org.ejml.data.DenseMatrix64F Java Examples

public PermutationIndices(DenseMatrix64F matrix) {

        this.matrix = matrix;
        dim = matrix.getNumCols();
        assert (dim == matrix.getNumRows());

        for (int i = 0; i < dim; ++i) {
            double diagonal = matrix.get(i, i);
            if (Double.isInfinite(diagonal)) {
                ++infiniteCount;
            } else if (diagonal == 0.0) {
                ++zeroCount;
            } else {
                ++nonZeroFiniteCount;
            }
        }
    }
 

public WrappedNormalSufficientStatistics(double[] buffer,
                                         int index,
                                         int dim,
                                         DenseMatrix64F Pd,
                                         PrecisionType precisionType) {

    int partialOffset = (dim + precisionType.getMatrixLength(dim)) * index;
    this.mean = new WrappedVector.Raw(buffer, partialOffset, dim);
    if (precisionType == PrecisionType.SCALAR) {
        this.precision = new WrappedMatrix.Raw(Pd.getData(), 0, dim, dim);
        this.precisionScalar = buffer[partialOffset + dim];
        this.variance = null;
    } else {
        this.precisionScalar = 1.0;
        this.precision = new WrappedMatrix.Raw(buffer, partialOffset + dim, dim, dim);
        this.variance = new WrappedMatrix.Raw(buffer, partialOffset + dim + dim * dim, dim, dim);
    }
}
 

public static double invertAndGetDeterminant(DenseMatrix64F mat, DenseMatrix64F result, boolean log) {

        final int numCol = mat.getNumCols();
        final int numRow = mat.getNumRows();
        if (numCol != numRow) {
            throw new IllegalArgumentException("Must be a square matrix.");
        }

        if (numCol <= 5) {

            if (numCol >= 2) {
                UnrolledInverseFromMinor.inv(mat, result);
            } else {
                result.set(0, 1.0D / mat.get(0));
            }

            double det = numCol >= 2 ?
                    UnrolledDeterminantFromMinor.det(mat) :
                    mat.get(0);
            return log ? Math.log(det) : det;

        } else {

            LUDecompositionAlt_D64 alg = new LUDecompositionAlt_D64();
            LinearSolverLu_D64 solver = new LinearSolverLu_D64(alg);
            if (solver.modifiesA()) {
                mat = mat.copy();
            }

            if (!solver.setA(mat)) {
                return Double.NaN;
            }

            solver.invert(result);

            return log ? computeLogDeterminant(alg) : alg.computeDeterminant().real;

        }
    }
 

private WrappedNormalSufficientStatistics getWeightedAverage(ReadableVector m1, ReadableMatrix p1,
                                                             ReadableVector m2, ReadableMatrix p2) {

    assert (m1.getDim() == m2.getDim());
    assert (p1.getDim() == p2.getDim());

    assert (m1.getDim() == p1.getMinorDim());
    assert (m1.getDim() == p1.getMajorDim());

    final WrappedVector m12 = new WrappedVector.Raw(new double[m1.getDim()], 0, dimFactors);
    final DenseMatrix64F p12 = new DenseMatrix64F(dimFactors, dimFactors);
    final DenseMatrix64F v12 = new DenseMatrix64F(dimFactors, dimFactors);

    final WrappedMatrix wP12 = new WrappedMatrix.WrappedDenseMatrix(p12);
    final WrappedMatrix wV12 = new WrappedMatrix.WrappedDenseMatrix(v12);

    MissingOps.add(p1, p2, wP12);
    safeInvert2(p12, v12, false);

    weightedAverage(m1, p1, m2, p2, m12, wV12, dimFactors);

    return new WrappedNormalSufficientStatistics(m12, wP12, wV12);
}
 

public static void gatherRowsAndColumns(final DenseMatrix64F source, final DenseMatrix64F destination,
                                        final int[] rowIndices, final int[] colIndices) {

    final int rowLength = rowIndices.length;
    final int colLength = colIndices.length;
    final double[] out = destination.getData();

    int index = 0;
    for (int i = 0; i < rowLength; ++i) {
        final int rowIndex = rowIndices[i];
        for (int j = 0; j < colLength; ++j) {
            out[index] = source.unsafe_get(rowIndex, colIndices[j]);
            ++index;
        }
    }
}
 

public static double det(DenseMatrix64F mat) {
    int numCol = mat.getNumCols();
    int numRow = mat.getNumRows();
    if (numCol != numRow) {
        throw new IllegalArgumentException("Must be a square matrix.");
    } else if (numCol <= 6) {
        return numCol >= 2 ? UnrolledDeterminantFromMinor.det(mat) : mat.get(0);
    } else {
        LUDecompositionAlt_D64 alg = new LUDecompositionAlt_D64();
        if (alg.inputModified()) {
            mat = mat.copy();
        }

        return !alg.decompose(mat) ? 0.0D : alg.computeDeterminant().real;
    }
}
 

public static InversionResult safeSolve(DenseMatrix64F A,
                                        WrappedVector b,
                                        WrappedVector x,
                                        boolean getDeterminat) {
    final int dim = b.getDim();

    assert (A.getNumRows() == dim && A.getNumCols() == dim);

    final DenseMatrix64F B = wrap(b.getBuffer(), b.getOffset(), dim, 1);
    final DenseMatrix64F X = new DenseMatrix64F(dim, 1);

    InversionResult ir = safeSolve(A, B, X, getDeterminat);


    for (int row = 0; row < dim; ++row) {
        x.set(row, X.unsafe_get(row, 0));
    }

    return ir;
}
 

@Override
void computePartialPrecision(int ido, int jdo, int imo, int jmo,
                             DenseMatrix64F Pip, DenseMatrix64F Pjp, DenseMatrix64F Pk) {

    final double[] diagQdi = vectorDiagQdi;
    System.arraycopy(diagonal1mActualizations, ido, diagQdi, 0, dimTrait);
    oneMinus(diagQdi);
    final double[] diagQdj = vectorDiagQdj;
    System.arraycopy(diagonal1mActualizations, jdo, diagQdj, 0, dimTrait);
    oneMinus(diagQdj);

    final DenseMatrix64F QdiPipQdi = matrix0;
    final DenseMatrix64F QdjPjpQdj = matrix1;
    diagonalDoubleProduct(Pip, diagQdi, QdiPipQdi);
    diagonalDoubleProduct(Pjp, diagQdj, QdjPjpQdj);
    CommonOps.add(QdiPipQdi, QdjPjpQdj, Pk);

    if (DEBUG) {
        System.err.println("Qdi: " + Arrays.toString(diagQdi));
        System.err.println("\tQdiPipQdi: " + QdiPipQdi);
        System.err.println("\tQdj: " + Arrays.toString(diagQdj));
        System.err.println("\tQdjPjpQdj: " + QdjPjpQdj);
    }
}
 

public static void weightedSumActualized(final double[] ipartial,
                                         final int ibo,
                                         final DenseMatrix64F Pi,
                                         final double[] iactualization,
                                         final int ido,
                                         final double[] jpartial,
                                         final int jbo,
                                         final DenseMatrix64F Pj,
                                         final double[] jactualization,
                                         final int jdo,
                                         final int dimTrait,
                                         final double[] out) {
    for (int g = 0; g < dimTrait; ++g) {
        double sum = 0.0;
        for (int h = 0; h < dimTrait; ++h) {
            sum += iactualization[ido + g] * Pi.unsafe_get(g, h) * ipartial[ibo + h];
            sum += jactualization[jdo + g] * Pj.unsafe_get(g, h) * jpartial[jbo + h];
        }
        out[g] = sum;
    }
}
 

@Override
public void getBranchExpectation(double[] actualization, double[] parentValue, double[] displacement,
                                 double[] expectation) {

    assert (expectation != null);
    assert (expectation.length >= dimTrait);

    assert (actualization != null);
    assert (actualization.length >= dimTrait * dimTrait);

    assert (parentValue != null);
    assert (parentValue.length >= dimTrait);

    assert (displacement != null);
    assert (displacement.length >= dimTrait);

    DenseMatrix64F branchExpectationMatrix = new DenseMatrix64F(dimTrait, 1);
    CommonOps.mult(wrap(actualization, 0, dimTrait, dimTrait),
            wrap(parentValue, 0, dimTrait, 1),
            branchExpectationMatrix);
    CommonOps.addEquals(branchExpectationMatrix, wrap(displacement, 0, dimTrait, 1));

    unwrap(branchExpectationMatrix, expectation, 0);
}
 

public static double weightedInnerProductOfDifferences(final double[] source1,
                                                       final int source1Offset,
                                                       final double[] source2,
                                                       final int source2Offset,
                                                       final DenseMatrix64F P,
                                                       final int dimTrait) {
    double SS = 0;
    for (int g = 0; g < dimTrait; ++g) {
        final double gDifference = source1[source1Offset + g] - source2[source2Offset + g];

        for (int h = 0; h < dimTrait; ++h) {
            final double hDifference = source1[source1Offset + h] - source2[source2Offset + h];

            SS += gDifference * P.unsafe_get(g, h) * hDifference;
        }
    }

    return SS;
}
 

private void computePartialsAndRemainders() {

        final DenseMatrix64F[] precisions = new DenseMatrix64F[taxonTaskPool.getNumThreads()];
        final DenseMatrix64F[] variances = new DenseMatrix64F[taxonTaskPool.getNumThreads()];

        for (int i = 0; i < taxonTaskPool.getNumThreads(); ++i) {
            precisions[i] = new DenseMatrix64F(numFactors, numFactors);
            variances[i] = new DenseMatrix64F(numFactors, numFactors);
        }

        if (USE_PRECISION_CACHE) {
            precisionMatrixMap.clear();
            if (DEBUG) {
                System.err.println("Hash CLEARED");
            }
        }

        if (TIMING) { // Do not use threads or lambda when timing
            for (int taxon = 0; taxon < numTaxa; ++taxon) {
                computePartialAndRemainderForOneTaxon(taxon, precisions[0], variances[0]);
            }
        } else {
            taxonTaskPool.fork((taxon, thread) ->
                    computePartialAndRemainderForOneTaxon(taxon, precisions[thread], variances[thread]));
        }
    }
 

@Override
void setMatrixRatio(DenseMatrix64F numeratorMatrix, DenseMatrix64F otherMatrix,
                    DenseMatrix64F destination) {
    int dim = destination.numRows;

    for (int i = 0; i < dim; i++) {
        for (int j = 0; j < dim; j++) {

            double n = Math.abs(numeratorMatrix.get(i, j));
            double d = Math.abs(otherMatrix.get(i, j));

            if (n == 0 && d == 0) {
                destination.set(i, j, 0);
            } else {
                destination.set(i, j, n / (n + d));
            }

        }
    }

}
 

@Override
void setMatrixRatio(DenseMatrix64F numeratorMatrix, DenseMatrix64F otherMatrix,
                    DenseMatrix64F destination) {

    for (int i = 0; i < destination.numRows; i++) {

        double val = numeratorMatrix.get(i, i) / (numeratorMatrix.get(i, i) + otherMatrix.get(i, i));
        destination.set(i, i, val);
        for (int j = i + 1; j < destination.numRows; j++) {

            double rg = numeratorMatrix.get(i, j);
            double vi = numeratorMatrix.get(i, i) + otherMatrix.get(i, i);
            double vj = numeratorMatrix.get(j, j) + otherMatrix.get(j, j);

            val = rg / sqrt(vi * vj);

            destination.set(i, j, val);
            destination.set(j, i, val);

        }
    }
}
 

public CompoundEigenMatrix(Parameter eigenValues, MatrixParameter eigenVectors) {
    super(eigenValues, eigenVectors);
    // Matrices
    temp = new DenseMatrix64F(dim, dim);
    transformedMatrix = new DenseMatrix64F(dim, dim);
    computeTransformedMatrix();
    savedTransformedMatrix = new double[dim * dim];
}
 

private void blockUnwrap(final DenseMatrix64F YY, final DenseMatrix64F XX, final DenseMatrix64F XY, final DenseMatrix64F YX, final double[] destination, final int offset) {
    for (int i = 0; i < dimProcess; i++) { // Rows
        for (int j = 0; j < dimProcess; j++) {
            destination[offset + i * dimTrait + j] = YY.get(i, j);
            destination[offset + (i + dimProcess) * dimTrait + j + dimProcess] = XX.get(i, j);
        }
        for (int j = 0; j < dimProcess; j++) {
            destination[offset + i * dimTrait + j + dimProcess] = YX.get(i, j);
            destination[offset + (i + dimProcess) * dimTrait + j] = XY.get(i, j);
        }
    }
}
 

@Override
public double[] chainRuleRoot(ContinuousDiffusionIntegrator cdi,
                              DiffusionProcessDelegate diffusionProcessDelegate,
                              ContinuousDataLikelihoodDelegate likelihoodDelegate,
                              BranchSufficientStatistics statistics, NodeRef node,
                              final DenseMatrix64F gradQInv, final DenseMatrix64F gradN) {
    return WRT_ROOT_MEAN.chainRuleRoot(cdi, diffusionProcessDelegate, likelihoodDelegate, statistics, node, gradQInv, gradN);
}
 

private DenseMatrix64F getGradientDisplacementWrtAttenuationDiagonal(ContinuousDiffusionIntegrator cdi, BranchSufficientStatistics statistics,
                                                                         NodeRef node, DenseMatrix64F gradient) {
        int nodeIndex = node.getNumber();
        // q_i
//        double[] qi = new double[dim];
//        cdi.getBranchActualization(getMatrixBufferOffsetIndex(nodeIndex), qi);
//        DenseMatrix64F qi = wrapDiagonal(actualization, 0, dim);
        // q_i^-1
//        DenseMatrix64F qiInv = wrapDiagonalInverse(actualization, 0, dim);
        // ni
        DenseMatrix64F ni = statistics.getAbove().getRawMean();
        // beta_i
        DenseMatrix64F betai = wrap(getDriftRate(node), 0, dim, 1);

        // factor
//        DenseMatrix64F tmp = new DenseMatrix64F(dim, 1);
        DenseMatrix64F resFull = new DenseMatrix64F(dim, dim);
        DenseMatrix64F resDiag = new DenseMatrix64F(dim, 1);
        CommonOps.add(ni, -1, betai, resDiag);
//        MissingOps.diagDiv(qi, resDiag);
        CommonOps.multTransB(gradient, resDiag, resFull);

        // Extract diag
        CommonOps.extractDiag(resFull, resDiag);

        // Wrt attenuation
        double ti = cdi.getBranchLength(getMatrixBufferOffsetIndex(nodeIndex));
        chainRuleActualizationWrtAttenuationDiagonal(ti, resDiag);

        return resDiag;

    }
 

@Override
public double[] chainRuleRoot(ContinuousDiffusionIntegrator cdi,
                              DiffusionProcessDelegate diffusionProcessDelegate,
                              ContinuousDataLikelihoodDelegate likelihoodDelegate,
                              BranchSufficientStatistics statistics, NodeRef node,
                              final DenseMatrix64F gradQInv, final DenseMatrix64F gradN) {
    return chainRule(cdi, diffusionProcessDelegate, likelihoodDelegate, statistics, node, gradQInv, gradN);
}
 

MatrixSufficientStatistics(double[] displacement,
                           double[] precision,
                           double[] actualization,
                           int index,
                           int dim,
                           DenseMatrix64F Pd,
                           PrecisionType precisionType) {

    super(displacement, precision, index, dim, Pd, precisionType);

    int actualizationOffset = (dim * dim) * index;
    this.actualization = MissingOps.wrap(actualization, actualizationOffset, dim, dim);

}
 

public static void transformMatrix(DenseMatrix64F matrix, DenseMatrix64F rotation, Boolean isSymmetric) {
    if (isSymmetric) {
        transformMatrixSymmetric(matrix, rotation);
    } else {
        transformMatrixGeneral(matrix, rotation);
    }
}
 

DenseMatrix64F getVarianceBranch(double branchPrecision){
    if (!hasDrift) {
        final DenseMatrix64F V1 = new DenseMatrix64F(dimTrait, dimTrait);
        CommonOps.scale(1.0 / branchPrecision, Vd, V1);
        return V1;
    } else {
        DenseMatrix64F P = getPrecisionBranch(branchPrecision);
        DenseMatrix64F V = new DenseMatrix64F(dimTrait, dimTrait);
        CommonOps.invert(P, V);
        return V;
    }
}
 

@Override
public EigenDecomposition decomposeStrenghtOfSelection(MatrixParameterInterface AParam, int dim, boolean isSymmetric) {
    DenseMatrix64F A = MissingOps.wrap(AParam);
    org.ejml.interfaces.decomposition.EigenDecomposition eigA = DecompositionFactory.eig(dim, true, isSymmetric);
    if (!eigA.decompose(A)) throw new RuntimeException("Eigen decomposition failed.");
    return new EigenDecomposition(eigenVectorsMatrix(eigA),
            null,
            eigenValuesMatrix(eigA, dim));
}
 

private DenseMatrix64F createMatrixA(int numRows, int numCols) {
    double[] data;
    if (hasIntercept) {
        data = fill(numRows, 1.0);
    } else {
        data = arrayFrom();
    }
    for (double[] predictor : predictors) {
        data = combine(data, arrayFrom(predictor));
    }
    boolean isRowMajor = false;
    return new DenseMatrix64F(numRows, numCols, isRowMajor, data);
}
 

DenseMatrix64F getPrecisionBranch(double branchPrecision){
    if (!hasDrift) {
        DenseMatrix64F P1 = new DenseMatrix64F(dimTrait, dimTrait);
        CommonOps.scale(branchPrecision, Pd, P1);
        return P1;
    } else {
        return DenseMatrix64F.wrap(dimTrait, dimTrait, precisionBuffer);
    }
}
 

private void solveSystem(int numRows, int numCols) {
    LinearSolver<DenseMatrix64F> qrSolver = LinearSolverFactory.qr(numRows, numCols);
    QRDecomposition<DenseMatrix64F> decomposition = qrSolver.getDecomposition();
    qrSolver.setA(X);
    y.setData(response);
    qrSolver.solve(this.y, this.b);
    DenseMatrix64F R = decomposition.getR(null, true);
    LinearSolver<DenseMatrix64F> linearSolver = LinearSolverFactory.linear(numCols);
    linearSolver.setA(R);
    DenseMatrix64F Rinverse = new DenseMatrix64F(numCols, numCols);
    linearSolver.invert(Rinverse); // stores solver's solution inside of Rinverse.
    CommonOps.multOuter(Rinverse, this.XtXInv);
}
 

public static int countFiniteNonZeroDiagonals(DenseMatrix64F source) {
    final int length = source.getNumCols();

    int count = 0;
    for (int i = 0; i < length; ++i) {
        final double d = source.unsafe_get(i, i);
        if (!Double.isInfinite(d) && d != 0.0) {
            ++count;
        }
    }
    return count;
}
 

@Override
    void computePartialPrecision(int ido, int jdo, int imo, int jmo,
                                 DenseMatrix64F Pip, DenseMatrix64F Pjp, DenseMatrix64F Pk) {

        final DenseMatrix64F Qdi = wrap(actualizations, imo, dimTrait, dimTrait);
        final DenseMatrix64F Qdj = wrap(actualizations, jmo, dimTrait, dimTrait);

        final DenseMatrix64F QdiPip = matrixQdiPip;
        final DenseMatrix64F QdiPipQdi = matrix0;
        scalePrecision(Qdi, Pip, QdiPip, QdiPipQdi);

        final DenseMatrix64F QdjPjpQdj = matrix1;
        final DenseMatrix64F QdjPjp = matrixQdjPjp;
        scalePrecision(Qdj, Pjp, QdjPjp, QdjPjpQdj);

        CommonOps.add(QdiPipQdi, QdjPjpQdj, Pk);

//        forceSymmetric(Pk);

        if (DEBUG) {
            System.err.println("Qdi: " + Qdi);
            System.err.println("\tQdiPip: " + QdiPip);
            System.err.println("\tQdiPipQdi: " + QdiPipQdi);
            System.err.println("\tQdj: " + Qdj);
            System.err.println("\tQdjPjp: " + QdjPjp);
            System.err.println("\tQdjPjpQdj: " + QdjPjpQdj);
        }
    }
 

private void makeCompletedUnobserved(final DenseMatrix64F matrix, double diagonal) {
    for (int row = 0; row < numFactors; ++row) {
        for (int col = 0; col < numFactors; ++col) {
            double x = (row == col) ? diagonal : 0.0;
            matrix.unsafe_set(row, col, x);
        }
    }
}
 

public static void weightedAverage(final double[] ipartial,
                                   final int ibo,
                                   final DenseMatrix64F Pi,
                                   final double[] jpartial,
                                   final int jbo,
                                   final DenseMatrix64F Pj,
                                   final double[] kpartial,
                                   final int kbo,
                                   final DenseMatrix64F Vk,
                                   final int dimTrait) {
    final double[] tmp = new double[dimTrait];
    weightedAverage(ipartial, ibo, Pi, jpartial, jbo, Pj, kpartial, kbo, Vk, dimTrait, tmp);
}