Create a matrix in wich each row is a tuple containing one of solutions or
    solution k-esima.
424597
def combine(items, k=None):
    
    length_items = len(items)
    lengths = [len(i) for i in items]
    length = reduce(lambda x, y: x * y, lengths)
    repeats = [reduce(lambda x, y: x * y, lengths[i:])
               for i in range(1, length_items)] + [1]
    if k is not None:
        k = k % length
        # Python division by default is integer division (~ floor(a/b))
        indices = [old_div((k % (lengths[i] * repeats[i])), repeats[i])
                   for i in range(length_items)]
        return [items[i][indices[i]] for i in range(length_items)]
    else:
        matrix = []
        for i, item in enumerate(items):
            row = []
            for subset in item:
                row.extend([subset] * repeats[i])
            times = old_div(length, len(row))
            matrix.append(row * times)
        # Transpose the matrix or return the columns instead rows
        return list(zip(*matrix))
542771
function(sbmA, sbmB) {
    // check if size compatible
    if (sbmA.y !== sbmB.x) {
      return
    }
    // create a new matrix
    return this.make((function (i, j) {
      return _.reduce(_.map(_.range(sbmA.y), (function(id) {
        return this.check(sbmA, i, id) && this.check(sbmB, id, j)
      }).bind(this)), function (sum, a) {
        return sum ^ a
      }) == 1
    }).bind(this), { x: sbmA.x, y: sbmB.y })
  }
230565
def fill
      @matrix[0][0] = 0
      # Set up the first column on each row.
      1.upto(@rows-1) {|i| @matrix[i][0] = @matrix[i-1][0] + @scoring.score_delete(@seq1[i])}
      # Set up the first row 
      1.upto(@cols-1) {|j| @matrix[0][j] = @matrix[0][j-1] + @scoring.score_insert(@seq2[j])}

      1.upto(@rows-1) do |i|
        prv_row = @matrix[i-1]
        cur_row = @matrix[i]

        1.upto(@cols-1) do |j|
          
          seq1_obj = @seq1[i-1]
          seq2_obj = @seq2[j-1]

          # Calculate the score.
          score_align = prv_row[j-1] + @scoring.score_align(seq1_obj, seq2_obj)
          score_delete = prv_row[j] + @scoring.score_delete(seq1_obj)
          score_insert = cur_row[j-1] + @scoring.score_insert(seq2_obj)
          max = max3(score_align, score_delete, score_insert)

          @matrix[i][j] = max
        end
      end
    end
565436
def _build_matrix(p, q, deriv):
    """"""
    A = [([1 for _ in range(-p, q+1)])]
    for i in range(1, p + q + 1):
        A.append([j**i for j in range(-p, q+1)])
    return np.array(A)
883931
def get_matrix(self):
        
        # Init matrix (remember, a matrix, not an array)
        a = np.zeros((4, 4), dtype=np.float32)
        w, x, y, z = self.w, self.x, self.y, self.z
        # First row
        a[0, 0] = - 2.0 * (y * y + z * z) + 1.0
        a[1, 0] = + 2.0 * (x * y + z * w)
        a[2, 0] = + 2.0 * (x * z - y * w)
        a[3, 0] = 0.0
        # Second row
        a[0, 1] = + 2.0 * (x * y - z * w)
        a[1, 1] = - 2.0 * (x * x + z * z) + 1.0
        a[2, 1] = + 2.0 * (z * y + x * w)
        a[3, 1] = 0.0
        # Third row
        a[0, 2] = + 2.0 * (x * z + y * w)
        a[1, 2] = + 2.0 * (y * z - x * w)
        a[2, 2] = - 2.0 * (x * x + y * y) + 1.0
        a[3, 2] = 0.0
        # Fourth row
        a[0, 3] = 0.0
        a[1, 3] = 0.0
        a[2, 3] = 0.0
        a[3, 3] = 1.0
        return a
625225
function matrix(rows, columns, values) {
    // console.log(values);
    const res = [];
    let row = [];
    let real;
    values.forEach((v, i) => {
        if (i % 2 == 0) {
            real = v;
        } else {
            row.push({
                real: real,
                imag: v
            });
            if ((i + 1) / 2 % columns == 0) {
                res.push(row);
                row = [];
            }
        }
    });
    return res;
}
569420
function() {
            var m = this.m,
                a = m.a,
                b = m.b,
                c = m.c,
                d = m.d,
                e = m.e,
                f = m.f,
                k, aa, bb, cc, dd, ee, ff;
            k = a * d - b * c;
            aa = d / k;
            bb = -b / k;
            cc = -c / k;
            dd = a / k;
            ee = (c * f - e * d) / k;
            ff = (b * e - a * f) / k;
            return new Matrix(aa, bb, cc, dd, ee, ff);
        }
238344
private DoubleMatrix solveEquationSystem(){
		DoubleMatrix R=new DoubleMatrix(partition.getLength());
		DoubleMatrix M=new DoubleMatrix(partition.getLength(),partition.getLength());
		DoubleMatrix partitionAsVector=new DoubleMatrix(partition.getPoints());
		DoubleMatrix shiftedPartition=new DoubleMatrix(partition.getLength());
		for(int j=1; j<shiftedPartition.length;j++) {
			shiftedPartition.put(j, partition.getPoint(j-1));
		}
		DoubleMatrix partitionIncrements= partitionAsVector.sub(shiftedPartition).put(0,1);
		DoubleMatrix kernelValues=new DoubleMatrix(partition.getLength()-1);
		DoubleMatrix M1_1= new DoubleMatrix(1);
		DoubleMatrix MFirstCol= new DoubleMatrix(partition.getLength()-1);
		DoubleMatrix MSubDiagonal= new DoubleMatrix(partition.getLength()-1);
		DoubleMatrix MSubMatrix= new DoubleMatrix(partition.getLength()-1,partition.getLength()-1);
		DoubleMatrix MSubMatrixSum= new DoubleMatrix(partition.getLength()-1);

		for(int i=0;i<independentValues.length;i++){

			DoubleMatrix oneZeroVector= new DoubleMatrix(partition.getLength());
			DoubleMatrix kernelSum= new DoubleMatrix(partition.getLength());
			DoubleMatrix shiftedKernelVector= new DoubleMatrix(partition.getLength());


			for(int r=0;r<partition.getLength()-1;r++){
				oneZeroVector.put(r, 1);
				kernelValues.put( r,kernel.density((partition.getIntervalReferencePoint(r)-independentValues[i])/bandwidth));
				shiftedKernelVector.put(r+1,kernelValues.get( r) );
				kernelSum=kernelSum.add(oneZeroVector.mmul(kernelValues.get(r)));
			}

			R=R.add(shiftedPartition.neg().add(independentValues[i]).mul(shiftedKernelVector)
					.add(partitionIncrements.mul(kernelSum)).mul(dependentValues[i]));

			M1_1=M1_1.add( kernelSum.get(0));

			MFirstCol=MFirstCol.add(
					partitionAsVector.getRange(0,partitionAsVector.length-1).neg().add(independentValues[i])
					.mul(kernelValues).add(
							partitionIncrements.getRange(1,partitionAsVector.length)
							.mul(kernelSum.getRange(1, kernelSum.length))));

			MSubDiagonal=MSubDiagonal.add(
					partitionAsVector.getRange(0,partitionAsVector.length-1).neg().add(independentValues[i])
					.mul(partitionAsVector.getRange(0,partitionAsVector.length-1).neg().add(independentValues[i]))
					.mul(kernelValues).add(
							partitionIncrements.getRange(1,partitionAsVector.length)
							.mul(partitionIncrements.getRange(1,partitionAsVector.length)
									.mul(kernelSum.getRange(1, kernelSum.length)))));

			MSubMatrixSum=MSubMatrixSum.add(
					partitionAsVector.getRange(0, partitionAsVector.length-1).neg().add(independentValues[i])
					.mul(kernelValues).add(
							partitionIncrements.getRange(1, partitionIncrements.length)
							.mul(kernelSum.getRange(1, kernelSum.length))));
		}

		DoubleMatrix partitionIncrementMatrix= new DoubleMatrix(partition.getLength()-1,partition.getLength()-1);
		DoubleMatrix matrixDefine= DoubleMatrix.ones(partition.getLength()-1);
		for(int m=0;m<matrixDefine.length-1;m++) {
			matrixDefine.put(m, 0);
			partitionIncrementMatrix.putColumn(m, matrixDefine.mul(partitionIncrements.get(m+1)));
		}

		MSubMatrix=partitionIncrementMatrix.mulColumnVector(MSubMatrixSum);
		MSubMatrix=MSubMatrix.add(MSubMatrix.transpose()).add(DoubleMatrix.diag(MSubDiagonal));

		int[] rowColIndex =new int[partition.getLength()-1];
		for(int n=0;n<rowColIndex.length;n++) {
			rowColIndex[n]=n+1;
		}

		M.put(0,0,M1_1.get(0));
		M.put(rowColIndex, 0, MFirstCol);
		M.put(0, rowColIndex, MFirstCol.transpose());
		M.put(rowColIndex, rowColIndex, MSubMatrix);

		return Solve.solve(M, R);
	}
874130
def kmask(matrix, k=7, dists=None, logic='or'):
    

    dists = (kdists(matrix, k=k) if dists is None else dists)
    mask = (matrix <= dists)
    if logic == 'or' or logic == '|':
        return mask | mask.T
    elif logic == 'and' or logic == '&':
        return mask & mask.T
    return mask
324694
public Matrix solve(Matrix rhs) {
        if (M != N || rhs.M != N || rhs.N != 1) {
            throw new RuntimeException("Illegal matrix dimensions.");
        }

        // create copies of the data
        Matrix A = new Matrix(this);
        Matrix b = new Matrix(rhs);

        // Gaussian elimination with partial pivoting
        for (int i = 0; i < N; i++) {

            // find pivot row and swap
            int max = i;
            for (int j = i + 1; j < N; j++) {
                if (Math.abs(A.data[j][i]) > Math.abs(A.data[max][i])) {
                    max = j;
                }
            }
            A.swap(i, max);
            b.swap(i, max);

            // singular
            if (A.data[i][i] == 0.0) {
                throw new RuntimeException("Matrix is singular.");
            }

            // pivot within b
            for (int j = i + 1; j < N; j++) {
                b.data[j][0] -= b.data[i][0] * A.data[j][i] / A.data[i][i];
            }

            // pivot within A
            for (int j = i + 1; j < N; j++) {
                double m = A.data[j][i] / A.data[i][i];
                for (int k = i + 1; k < N; k++) {
                    A.data[j][k] -= A.data[i][k] * m;
                }
                A.data[j][i] = 0.0;
            }
        }

        // back substitution
        Matrix x = new Matrix(N, 1);
        for (int j = N - 1; j >= 0; j--) {
            double t = 0.0;
            for (int k = j + 1; k < N; k++) {
                t += A.data[j][k] * x.data[k][0];
            }
            x.data[j][0] = (b.data[j][0] - t) / A.data[j][j];
        }
        return x;

    }
1051600
public DenseMatrix solve(DenseMatrix B) throws MatrixNotSPDException {
        if (notspd)
            throw new MatrixNotSPDException();
        if (B.numRows() != n)
            throw new IllegalArgumentException("B.numRows() != n");

        intW info = new intW(0);
        if (upper)
            LAPACK.getInstance().dpbtrs(UpLo.Upper.netlib(), n, kd,
                    B.numColumns(), Cu.getData(), Matrices.ld(kd + 1),
                    B.getData(), Matrices.ld(n), info);
        else
            LAPACK.getInstance().dpbtrs(UpLo.Lower.netlib(), n, kd,
                    B.numColumns(), Cl.getData(), Matrices.ld(kd + 1),
                    B.getData(), Matrices.ld(n), info);

        if (info.val < 0)
            throw new IllegalArgumentException();

        return B;
    }
37046
function() {
                var rows = this.elements.length, i, cols = this.elements[0].length, j;
            	var elements = [], i = cols;
            	while (i--) {
            	    j = rows;
            	    elements[i] = [];
            	    while (j--) {
            		elements[i][j] = this.elements[j][i];
            	    }
            	}
            	return Matrix$b.create(elements);
                }
927880
public void solve(DenseMatrix B) {
        if (B.nrows() != L.nrows()) {
            throw new IllegalArgumentException(String.format("Row dimensions do not agree: A is %d x %d, but B is %d x %d", L.nrows(), L.ncols(), B.nrows(), B.ncols()));
        }

        int n = B.nrows();
        int nrhs = B.ncols();

        // Solve L*Y = B;
        for (int k = 0; k < n; k++) {
            for (int j = 0; j < nrhs; j++) {
                for (int i = 0; i < k; i++) {
                    B.sub(k, j, B.get(i, j) * L.get(k, i));
                }
                B.div(k, j, L.get(k, k));
            }
        }

        // Solve L'*X = Y;
        for (int k = n - 1; k >= 0; k--) {
            for (int j = 0; j < nrhs; j++) {
                for (int i = k + 1; i < n; i++) {
                    B.sub(k, j, B.get(i, j) * L.get(i, k));
                }
                B.div(k, j, L.get(k, k));
            }
        }
    }
600535
def Ry_matrix(theta):
    """"""
    return np.array([
        [np.cos(theta), 0, np.sin(theta)],
        [0, 1, 0],
        [-np.sin(theta), 0, np.cos(theta)]
    ])
646285
def matrix_multiply(m0, m1):
    new_matrix = numpy.zeros((m0.shape[0],m1.shape[1]))
    ""
    for i in xrange(m0.shape[0]):
        for j in xrange(m1.shape[1]):
            r=0
            for k in xrange(m1.shape[0]):
                r += m0[i,k]*m1[k,j]
            new_matrix[i,j]=r
    return new_matrix
486507
@Override
    public void solve(DMatrixRMaj B, DMatrixRMaj X) {
        if( B.numRows != numRows )
            throw new IllegalArgumentException("Unexpected dimensions for X: X rows = "+X.numRows+" expected = "+numCols);
        X.reshape(numCols,B.numCols);

        int BnumCols = B.numCols;

        Y.reshape(numRows,1, false);
        Z.reshape(numRows,1, false);

        // solve each column one by one
        for( int colB = 0; colB < BnumCols; colB++ ) {

            // make a copy of this column in the vector
            for( int i = 0; i < numRows; i++ ) {
                Y.data[i] = B.get(i,colB);
            }

            // Solve Qa=b
            // a = Q'b
            CommonOps_DDRM.multTransA(Q,Y,Z);

            // solve for Rx = b using the standard upper triangular solver
            TriangularSolver_DDRM.solveU(R.data,Z.data,numCols);

            // save the results
            for( int i = 0; i < numCols; i++ ) {
                X.set(i,colB,Z.data[i]);
            }
        }
    }
728440
def Tmatrix(X):
    
    T = [[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]]
    for row in X:
        for k in range(3):
            for l in range(3):
                T[k][l] += row[k] * row[l]
    return T