@hidden
Perform LU decomposition

$$ A = P L U $$
351348
function lu_custom(A) {
    // Outer product LU decomposition with partial pivoting
    // Ref: Algo 3.4.1 Golub and Van Loan
    if (A.shape.length != 2) {
        throw new Error('Input is not a Matrix (2D)');
    }
    if (A.shape[0] !== A.shape[1]) {
        throw new Error('Input is not a Square Matrix');
    }
    var n = A.shape[0];
    var perm = new common_1.NDArray({ shape: [n] });
    for (var i = 0; i < n; i++) {
        perm.set(i, i);
    }
    function recordPermutation(ri, rj) {
        var tmp = perm.get(ri);
        perm.set(ri, perm.get(rj));
        perm.set(rj, tmp);
    }
    for (var k = 0; k < n - 1; k++) {
        // Find the maximum absolute entry in k'th column
        var ipivot = 0;
        var pivot = -Infinity;
        for (var i = k; i < n; i++) {
            var val = Math.abs(A.get(i, k)); // TODO:Complex
            if (val > pivot) {
                pivot = val;
                ipivot = i;
            }
        }
        // Swap rows k and ipivot
        A.swaprows(k, ipivot);
        recordPermutation(k, ipivot);
        if (common_1.iszero(pivot)) {
            throw new Error('Can\'t perform LU decomp. 0 on diagonal');
        }
        for (var i = k + 1; i < n; i++) {
            A.set(i, k, A.get(i, k) / pivot); // TODO:Complex
        }
        for (var i = k + 1; i < n; i++) {
            for (var j = n - 1; j > k; j--) {
                // TODO:Complex
                A.set(i, j, A.get(i, j) - A.get(i, k) * A.get(k, j));
            }
        }
    }
    return perm;
}
781430
public static Scsn cs_lu(Scs A, Scss S, float tol) {
        Scs L, U;
        Scsn N;
        float pivot, Lx[], Ux[], x[], a, t;
        int Lp[], Li[], Up[], Ui[], pinv[], xi[], q[], n, ipiv, k, top, p, i, col, lnz, unz;
        if (!Scs_util.CS_CSC(A) || S == null)
            return (null); /* check inputs */
        n = A.n;
        q = S.q;
        lnz = S.lnz;
        unz = S.unz;
        x = new float[n]; /* get float workspace */
        xi = new int[2 * n]; /* get int workspace */
        N = new Scsn(); /* allocate result */
        N.L = L = Scs_util.cs_spalloc(n, n, lnz, true, false); /* allocate result L */
        N.U = U = Scs_util.cs_spalloc(n, n, unz, true, false); /* allocate result U */
        N.pinv = pinv = new int[n]; /* allocate result pinv */
        Lp = L.p;
        Up = U.p;
        for (i = 0; i < n; i++)
            x[i] = 0; /* clear workspace */
        for (i = 0; i < n; i++)
            pinv[i] = -1; /* no rows pivotal yet */
        for (k = 0; k <= n; k++)
            Lp[k] = 0; /* no cols of L yet */
        lnz = unz = 0;
        for (k = 0; k < n; k++) /* compute L(:,k) and U(:,k) */
        {
            /* --- Triangular solve --------------------------------------------- */
            Lp[k] = lnz; /* L(:,k) starts here */
            Up[k] = unz; /* U(:,k) starts here */
            if (lnz + n > L.nzmax) {
                Scs_util.cs_sprealloc(L, 2 * L.nzmax + n);
            }
            if (unz + n > U.nzmax) {
                Scs_util.cs_sprealloc(U, 2 * U.nzmax + n);
            }
            Li = L.i;
            Lx = L.x;
            Ui = U.i;
            Ux = U.x;
            col = q != null ? (q[k]) : k;
            top = Scs_spsolve.cs_spsolve(L, A, col, xi, x, pinv, true); /* x = L\A(:,col) */
            /* --- Find pivot --------------------------------------------------- */
            ipiv = -1;
            a = -1;
            for (p = top; p < n; p++) {
                i = xi[p]; /* x(i) is nonzero */
                if (pinv[i] < 0) /* row i is not yet pivotal */
                {
                    if ((t = Math.abs(x[i])) > a) {
                        a = t; /* largest pivot candidate so far */
                        ipiv = i;
                    }
                } else /* x(i) is the entry U(pinv[i],k) */
                {
                    Ui[unz] = pinv[i];
                    Ux[unz++] = x[i];
                }
            }
            if (ipiv == -1 || a <= 0)
                return (null);
            if (pinv[col] < 0 && Math.abs(x[col]) >= a * tol)
                ipiv = col;
            /* --- Sivide by pivot ---------------------------------------------- */
            pivot = x[ipiv]; /* the chosen pivot */
            Ui[unz] = k; /* last entry in U(:,k) is U(k,k) */
            Ux[unz++] = pivot;
            pinv[ipiv] = k; /* ipiv is the kth pivot row */
            Li[lnz] = ipiv; /* first entry in L(:,k) is L(k,k) = 1 */
            Lx[lnz++] = 1;
            for (p = top; p < n; p++) /* L(k+1:n,k) = x / pivot */
            {
                i = xi[p];
                if (pinv[i] < 0) /* x(i) is an entry in L(:,k) */
                {
                    Li[lnz] = i; /* save unpermuted row in L */
                    Lx[lnz++] = x[i] / pivot; /* scale pivot column */
                }
                x[i] = 0; /* x [0.f.n-1] = 0 for next k */
            }
        }
        /* --- Finalize L and U ------------------------------------------------- */
        Lp[n] = lnz;
        Up[n] = unz;
        Li = L.i; /* fix row indices of L for final pinv */
        for (p = 0; p < lnz; p++)
            Li[p] = pinv[Li[p]];
        Scs_util.cs_sprealloc(L, 0); /* remove extra space from L and U */
        Scs_util.cs_sprealloc(U, 0);
        return N;
    }
799143
public static boolean isRowsLinearIndependent( DMatrixRMaj A )
    {
        // LU decomposition
        LUDecomposition<DMatrixRMaj> lu = DecompositionFactory_DDRM.lu(A.numRows,A.numCols);
        if( lu.inputModified() )
            A = A.copy();

        if( !lu.decompose(A))
            throw new RuntimeException("Decompositon failed?");

        // if they are linearly independent it should not be singular
        return !lu.isSingular();
    }
892904
def split_linear_constraints(A, l, u):
    
    ieq = []
    igt = []
    ilt = []
    ibx = []
    for i in range(len(l)):
        if abs(u[i] - l[i]) <= EPS:
            ieq.append(i)
        elif (u[i] > 1e10) and (l[i] > -1e10):
            igt.append(i)
        elif (l[i] <= -1e10) and (u[i] < 1e10):
            ilt.append(i)
        elif (abs(u[i] - l[i]) > EPS) and (u[i] < 1e10) and (l[i] > -1e10):
            ibx.append(i)
        else:
            raise ValueError

    Ae = A[ieq, :]
    Ai = sparse([A[ilt, :], -A[igt, :], A[ibx, :], -A[ibx, :]])
    be = u[ieq, :]
    bi = matrix([u[ilt], -l[igt], u[ibx], -l[ibx]])

    return Ae, be, Ai, bi
1303325
public static void symmRankNMinus_U(int blockLength ,
                                        DSubmatrixD1 A , DSubmatrixD1 B )
    {

        int heightB = B.row1-B.row0;
        if( heightB > blockLength )
            throw new IllegalArgumentException("Height of B cannot be greater than the block length");

        int N = B.col1-B.col0;

        if( A.col1-A.col0 != N )
            throw new IllegalArgumentException("A does not have the expected number of columns based on B's width");
        if( A.row1-A.row0 != N )
            throw new IllegalArgumentException("A does not have the expected number of rows based on B's width");


        for( int i = B.col0; i < B.col1; i += blockLength ) {

            int indexB_i = B.row0*B.original.numCols + i*heightB;
            int widthB_i = Math.min(blockLength,B.col1-i);

            int rowA = i-B.col0+A.row0;
            int heightA = Math.min( blockLength , A.row1 - rowA);

            for( int j = i; j < B.col1; j += blockLength ) {

                int widthB_j = Math.min(blockLength,B.col1-j);

                int indexA = rowA * A.original.numCols + (j-B.col0+A.col0)*heightA;
                int indexB_j = B.row0*B.original.numCols + j*heightB;

                if( i == j ) {
                    // only the upper portion of this block needs to be modified since it is along a diagonal
                    multTransABlockMinus_U( B.original.data,A.original.data,
                            indexB_i,indexB_j,indexA,heightB,widthB_i,widthB_j);
                } else {
                    multTransABlockMinus( B.original.data,A.original.data,
                            indexB_i,indexB_j,indexA,heightB,widthB_i,widthB_j);
                }
            }
        }
    }
559212
function solveIntersection (a, b, c, d) {
  var ba = rvSub(b, a)
  var dc = rvSub(d, c)

  var baXdc = ratPerp(ba, dc)

  if (ratSign(baXdc) === 0) {
    return null
  }

  var ac = rvSub(a, c)
  var dcXac = ratPerp(dc, ac)

  var t = ratDiv(dcXac, baXdc)
  var s = rvMuls(ba, t)
  var r = rvAdd(a, s)

  return r
}
1108314
function pnpoly (x,y,coords) {
    var vert = [ [0,0] ]

    for (var i = 0; i < coords.length; i++) {
      for (var j = 0; j < coords[i].length; j++) {
        vert.push(coords[i][j])
      }
      vert.push([0,0])
    }

    var inside = false
    for (var i = 0, j = vert.length - 1; i < vert.length; j = i++) {
      if (((vert[i][0] > y) != (vert[j][0] > y)) && (x < (vert[j][1] - vert[i][1]) * (y - vert[i][0]) / (vert[j][0] - vert[i][0]) + vert[i][1])) inside = !inside
    }

    return inside
  }
611930
def isoncurve(self, p):
        
        return p.iszero() or p.y ** 2 == p.x ** 3 + self.a * p.x + self.b
494766
def _u0Eq(logu,delta,pot,E,Lz22):
    """"""
    u= numpy.exp(logu)
    sinh2u= numpy.sinh(u)**2.
    cosh2u= numpy.cosh(u)**2.
    dU= cosh2u*actionAngleStaeckel.potentialStaeckel(u,numpy.pi/2.,pot,delta)
    return -(E*sinh2u-dU-Lz22/delta**2./sinh2u)
1195902
def lab2xyz(l__, a__, b__):
    


    def f_inv(arr):
        """Inverse of the f function.
        """
        return np.where(arr > 6.0/29.0,
                        arr ** 3,
                        3 * (6.0 / 29.0) * (6.0 / 29.0) * (arr - 4.0 / 29.0))

    new_l = (l__ + 16.0) / 116.0
    x__ = 95.047 * f_inv(new_l + a__ / 500.0)
    y__ = 100.0 * f_inv(new_l)
    z__ = 108.883 * f_inv(new_l - b__ / 200.0)

    return x__, y__, z__
169894
def relaxation_operators(p):
    
    k0 = np.array([[1.0, 0.0], [0.0, np.sqrt(1 - p)]])
    k1 = np.array([[0.0, np.sqrt(p)], [0.0, 0.0]])
    return k0, k1
622077
function vec3lerp(out, a, b, t, oa, ob) {
    var ax = a[oa];
    var ay = a[oa + 1];
    var az = a[oa + 2];
    out[0] = ax + t * (b[ob] - ax);
    out[1] = ay + t * (b[ob + 1] - ay);
    out[2] = az + t * (b[ob + 2] - az);

    return out;
}
1212879
def linsolve(self, A, b):
        
        if self.sparselib == 'umfpack':
            return umfpack.linsolve(A, b)

        elif self.sparselib == 'klu':
            return klu.linsolve(A, b)
860019
public void decompose( DMatrixRMaj A ) {

        this.QR = A.copy();

        int N = Math.min(A.numCols,A.numRows);

        gammas = new double[ A.numCols ];

        DMatrixRMaj A_small = new DMatrixRMaj(A.numRows,A.numCols);
        DMatrixRMaj A_mod = new DMatrixRMaj(A.numRows,A.numCols);
        DMatrixRMaj v = new DMatrixRMaj(A.numRows,1);
        DMatrixRMaj Q_k = new DMatrixRMaj(A.numRows,A.numRows);

        for( int i = 0; i < N; i++ ) {
            // reshape temporary variables
            A_small.reshape(QR.numRows-i,QR.numCols-i);
            v.reshape(A_small.numRows,1);
            Q_k.reshape(v.getNumElements(),v.getNumElements());

            // use extract matrix to get the column that is to be zeroed
            CommonOps_DDRM.extract(QR,i,QR.numRows,i,i+1,v,0,0);

            double max = CommonOps_DDRM.elementMaxAbs(v);

            if( max > 0 && v.getNumElements() > 1 ) {
                // normalize to reduce overflow issues
                CommonOps_DDRM.divide(v,max);

                // compute the magnitude of the vector
                double tau = NormOps_DDRM.normF(v);

                if( v.get(0) < 0 )
                    tau *= -1.0;

                double u_0 = v.get(0) + tau;
                double gamma = u_0/tau;

                CommonOps_DDRM.divide(v,u_0);
                v.set(0,1.0);

                // extract the submatrix of A which is being operated on
                CommonOps_DDRM.extract(QR,i,QR.numRows,i,QR.numCols,A_small,0,0);

                // A = (I - &gamma;*u*u<sup>T</sup>)A
                CommonOps_DDRM.setIdentity(Q_k);
                CommonOps_DDRM.multAddTransB(-gamma,v,v,Q_k);
                CommonOps_DDRM.mult(Q_k,A_small,A_mod);

                // save the results
                CommonOps_DDRM.insert(A_mod, QR, i,i);
                CommonOps_DDRM.insert(v, QR, i,i);
                QR.unsafe_set(i,i,-tau*max);

                // save gamma for recomputing Q later on
                gammas[i] = gamma;
            }
        }
    }
262917
def solve_factorized_aug(z, Fval, LU, G, A):
    M, N=G.shape
    P, N=A.shape

    """"""
    m=M    

    """Primal variable"""
    x=z[0:N]

    """Multiplier for equality constraints"""
    nu=z[N:N+P]

    """Multiplier for inequality constraints"""
    l=z[N+P:N+P+M]

    """Slacks"""
    s=z[N+P+M:]

    """Dual infeasibility"""
    rd = Fval[0:N]
    
    """Primal infeasibility"""
    rp1 = Fval[N:N+P]
    rp2 = Fval[N+P:N+P+M]

    """Centrality"""
    rc = Fval[N+P+M:]

    """Sigma matrix"""
    SIG = diags(l/s, 0)

    """LU is actually the augmented system W"""
    W = LU

    b1 = -rd - mydot(G.T, mydot(SIG, rp2)) + mydot(G.T, rc/s)
    b2 = -rp1
    b = np.hstack((b1, b2))

    """Prepare iterative solve via MINRES"""
    sign = np.zeros(N+P)
    sign[0:N//2] = 1.0
    sign[N//2:] = -1.0
    T = diags(sign, 0)        
   
    """Change rhs"""
    b_new = mydot(T, b)

    dW = np.abs(W.diagonal())
    dPc = np.ones(W.shape[0])
    ind = (dW > 0.0)
    dPc[ind] = 1.0/dW[ind]
    Pc = diags(dPc, 0)    
    dxnu, info = minres(W, b_new, tol=1e-10, M=Pc)
    
    # dxnu = solve(J, b)
    dx = dxnu[0:N]
    dnu = dxnu[N:]

    """Obtain search directions for l and s"""
    ds = -rp2 - mydot(G, dx)
    # ds = s*ds
    # SIG = np.diag(l/s)
    dl = -mydot(SIG, ds) - rc/s

    dz = np.hstack((dx, dnu, dl, ds))
    return dz
120968
def _solve(self, x0, A, l, u, xmin, xmax):
        
        s = pips(self._costfcn, x0, A, l, u, xmin, xmax,
                 self._consfcn, self._hessfcn, self.opt)
        return s
440898
function ecp_mul(R, n1, N2, a, p, l) {
	var add = mpn_one_p(N2[2], l) ? ecp_add_aff : ecp_add;
	var active = false, swaps = 0, S = _ecp_S, T;
	for (var i = l; i > 0;) {
		var w = n1[--i];
		for (var j = 31; j > 0; --j) {
			if (active) {
				ecp_dbl(S, R, a, p, l);
				T = S, S = R, R = T, ++swaps;
			}
			if ((w <<= 1) < 0) {
				if (active) {
					add(S, R, N2, a, p, l);
					T = S, S = R, R = T, ++swaps;
				}
				else {
					ecp_copy(R, N2, l);
					active = true;
				}
			}
		}
	}
	if (swaps & 1) {
		return ecp_copy(S, R, l);
	}
	return R;
}