Java Code Examples for org.ejml.ops.CommonOps#invert()

private void computeInverse() {

        if (DEBUG) {
            System.err.println("CachedMatrixInverse.computeInverse()");
        }
        
        if (EMJL) {
            // TODO Avoid multiple copies
            DenseMatrix64F source = new DenseMatrix64F(base.getParameterAsMatrix());
            DenseMatrix64F destination = new DenseMatrix64F(getColumnDimension(), getColumnDimension());
            CommonOps.invert(source, destination);
            inverse = new WrappedMatrix.WrappedDenseMatrix(destination);
        } else {
         inverse = new WrappedMatrix.ArrayOfArray(new Matrix(base.getParameterAsMatrix()).inverse().toComponents());
        }
    }
 

@Override
public void setDiffusionPrecision(int precisionIndex, final double[] matrix, double logDeterminant) {
    super.setDiffusionPrecision(precisionIndex, matrix, logDeterminant);

    assert (inverseDiffusions != null);

    final int offset = dimProcess * dimProcess * precisionIndex;
    DenseMatrix64F precision = wrap(diffusions, offset, dimProcess, dimProcess);
    DenseMatrix64F variance = new DenseMatrix64F(dimProcess, dimProcess);
    CommonOps.invert(precision, variance);
    unwrap(variance, inverseDiffusions, offset);

    if (DEBUG) {
        System.err.println("At precision index: " + precisionIndex);
        System.err.println("precision: " + precision);
        System.err.println("variance : " + variance);
    }
}
 

private void schurComplementInverse(final DenseMatrix64F A, final DenseMatrix64F D,
                                    final DenseMatrix64F C, final DenseMatrix64F B,
                                    final double[] destination, final int offset) {
    DenseMatrix64F invA = new DenseMatrix64F(dimProcess, dimProcess);
    CommonOps.invert(A, invA);
    DenseMatrix64F invMatD = getSchurInverseComplement(invA, D, C, B);

    DenseMatrix64F invAB = new DenseMatrix64F(dimProcess, dimProcess);
    CommonOps.mult(invA, B, invAB);
    DenseMatrix64F invMatB = new DenseMatrix64F(dimProcess, dimProcess);
    CommonOps.mult(-1.0, invAB, invMatD, invMatB);

    DenseMatrix64F CinvA = new DenseMatrix64F(dimProcess, dimProcess);
    CommonOps.mult(C, invA, CinvA);
    DenseMatrix64F invMatC = new DenseMatrix64F(dimProcess, dimProcess);
    CommonOps.mult(-1.0, invMatD, CinvA, invMatC);

    DenseMatrix64F invMatA = new DenseMatrix64F(dimProcess, dimProcess);
    CommonOps.mult(-1.0, invMatB, CinvA, invMatA);
    CommonOps.addEquals(invMatA, invA);

    blockUnwrap(invMatA, invMatD, invMatC, invMatB, destination, offset);
}
 

public static void symmPosDefInvert(DenseMatrix64F P, DenseMatrix64F P_inv) {
    LinearSolver<DenseMatrix64F> solver = LinearSolverFactory.symmPosDef(P.getNumCols());
    DenseMatrix64F Pbis = new DenseMatrix64F(P);
    if (!solver.setA(Pbis)) {
        CommonOps.invert(P, P_inv);
    } else {
        solver.invert(P_inv);
    }
}
 

private void computeTransformedMatrix() {
    DenseMatrix64F baseMatrix = wrapSpherical(offDiagonalParameter.getParameterValues(), 0, dim);
    DenseMatrix64F diagonalMatrix = wrapDiagonal(diagonalParameter.getParameterValues(), 0, dim);

    CommonOps.mult(baseMatrix, diagonalMatrix, temp);
    CommonOps.invert(baseMatrix);
    CommonOps.mult(temp, baseMatrix, transformedMatrix);

    compositionKnown = true;
}
 

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;
    }
}
 

public ConditionalPrecisionAndTransform2(final DenseMatrix64F precision,
                                         final int[] missingIndices,
                                         final int[] notMissingIndices) {

    assert (missingIndices.length + notMissingIndices.length == precision.getNumRows());
    assert (missingIndices.length + notMissingIndices.length == precision.getNumCols());

    this.missingIndices = missingIndices;
    this.notMissingIndices = notMissingIndices;

    DenseMatrix64F P11 = new DenseMatrix64F(missingIndices.length, missingIndices.length);
    MissingOps.gatherRowsAndColumns(precision, P11, missingIndices, missingIndices);

    P11Inv = new DenseMatrix64F(missingIndices.length, missingIndices.length);
    CommonOps.invert(P11, P11Inv);

    DenseMatrix64F P12 = new DenseMatrix64F(missingIndices.length, notMissingIndices.length);
    MissingOps.gatherRowsAndColumns(precision, P12, missingIndices, notMissingIndices);

    DenseMatrix64F P11InvP12 = new DenseMatrix64F(missingIndices.length, notMissingIndices.length);
    CommonOps.mult(P11Inv, P12, P11InvP12);

    this.affineTransform = P11InvP12;

    this.numMissing = missingIndices.length;
    this.numNotMissing = notMissingIndices.length;
}
 

private static void transformMatrixGeneral(DenseMatrix64F matrix, DenseMatrix64F rotation) {
    int dim = matrix.getNumRows();
    DenseMatrix64F tmp = new DenseMatrix64F(dim, dim);
    DenseMatrix64F rotationInverse = new DenseMatrix64F(dim, dim);
    CommonOps.invert(rotation, rotationInverse);
    CommonOps.mult(rotationInverse, matrix, tmp);
    CommonOps.multTransB(tmp, rotationInverse, matrix);
}
 

private DenseMatrix64F getSchurInverseComplement(final DenseMatrix64F invA, final DenseMatrix64F D,
                                    final DenseMatrix64F C, final DenseMatrix64F B) {
    DenseMatrix64F complement = new DenseMatrix64F(dimProcess, dimProcess);
    DenseMatrix64F tmp = new DenseMatrix64F(dimProcess, dimProcess);
    CommonOps.mult(invA, B, tmp);
    CommonOps.mult(-1.0, C, tmp, complement);
    CommonOps.addEquals(complement, D);
    CommonOps.invert(complement);
    return complement;
}
 

private double sampleNode(NodeRef node, WrappedNormalSufficientStatistics statistics) {

        final int thisNumber = node.getNumber();
        final int[] obsInds = maskDelegate.getObservedIndices(node);
        final int obsDim = obsInds.length;

        if (obsDim == dim) {
            return 0;
        }

        ReadableVector mean = statistics.getMean();
        ReadableMatrix thisP = statistics.getPrecision();
        double precisionScalar = statistics.getPrecisionScalar();
        for (int i = 0; i < mean.getDim(); ++i) {
            fcdMean[i] = mean.get(i);
        }

        for (int i = 0; i < mean.getDim(); ++i) {
            for (int j = 0; j < mean.getDim(); ++j) {
                fcdPrecision[i][j] = thisP.get(i, j) * precisionScalar;
            }
        }

        if (repeatedMeasuresModel != null) {
            //TODO: preallocate memory for all these matrices/vectors
            DenseMatrix64F Q = new DenseMatrix64F(fcdPrecision); //storing original fcd precision
            double[] tipPartial = repeatedMeasuresModel.getTipPartial(thisNumber, false);
            int offset = dim;
            DenseMatrix64F P = MissingOps.wrap(tipPartial, offset, dim, dim);

            DenseMatrix64F preOrderMean = new DenseMatrix64F(dim, 1);
            for (int i = 0; i < dim; i++) {
                preOrderMean.set(i, 0, fcdMean[i]);
            }
            DenseMatrix64F dataMean = new DenseMatrix64F(dim, 1);
            for (int i = 0; i < dim; i++) {
                dataMean.set(i, 0, tipPartial[i]);
            }

            D1Matrix64F bufferMean = new DenseMatrix64F(dim, 1);

            mult(Q, preOrderMean, bufferMean); //bufferMean = Q * preOrderMean
            multAdd(P, dataMean, bufferMean); //bufferMean = Q * preOderMean + P * dataMean


            CommonOps.addEquals(P, Q); //P = P + Q
            DenseMatrix64F V = new DenseMatrix64F(dim, dim);
            CommonOps.invert(P, V); //V = inv(P + Q)
            mult(V, bufferMean, dataMean); // dataMean = inv(P + Q) * (Q * preOderMean + P * dataMean)

            for (int i = 0; i < dim; i++) {
                fcdMean[i] = dataMean.get(i);
                for (int j = 0; j < dim; j++) {
                    fcdPrecision[i][j] = P.get(i, j);
                }

            }
        }

        MultivariateNormalDistribution fullDistribution = new MultivariateNormalDistribution(fcdMean, fcdPrecision); //TODO: should this not be declared until 'else' statement?
        MultivariateNormalDistribution drawDistribution;

        if (mask != null && obsDim > 0) {
            addMaskOnContiuousTraitsPrecisionSpace(thisNumber);
            drawDistribution = new MultivariateNormalDistribution(maskedMean, maskedPrecision);
        } else {
            drawDistribution = fullDistribution;
        }

        double[] oldValue = getNodeTrait(node);

        int attempt = 0;
        boolean validTip = false;

        while (!validTip & attempt < maxAttempts) {

            setNodeTrait(node, drawDistribution.nextMultivariateNormal());

            if (latentLiability.validTraitForTip(thisNumber)) {
                validTip = true;
            }
            attempt++;
        }

        if (attempt == maxAttempts) {
            return Double.NEGATIVE_INFINITY;
        }

        double[] newValue = getNodeTrait(node);

        if (latentLiability instanceof DummyLatentTruncationProvider) {
            return Double.POSITIVE_INFINITY;
        } else {
            return fullDistribution.logPdf(oldValue) - fullDistribution.logPdf(newValue);
        }
    }
 

public ConditionalVarianceAndTransform2(final DenseMatrix64F variance,
                                        final int[] missingIndices, final int[] notMissingIndices) {

    assert (missingIndices.length + notMissingIndices.length == variance.getNumRows());
    assert (missingIndices.length + notMissingIndices.length == variance.getNumCols());

    this.missingIndices = missingIndices;
    this.notMissingIndices = notMissingIndices;

    if (DEBUG) {
        System.err.println("variance:\n" + variance);
    }

    DenseMatrix64F S22 = new DenseMatrix64F(notMissingIndices.length, notMissingIndices.length);
    gatherRowsAndColumns(variance, S22, notMissingIndices, notMissingIndices);

    if (DEBUG) {
        System.err.println("S22:\n" + S22);
    }

    DenseMatrix64F S22Inv = new DenseMatrix64F(notMissingIndices.length, notMissingIndices.length);
    CommonOps.invert(S22, S22Inv);

    if (DEBUG) {
        System.err.println("S22Inv:\n" + S22Inv);
    }

    DenseMatrix64F S12 = new DenseMatrix64F(missingIndices.length, notMissingIndices.length);
    gatherRowsAndColumns(variance, S12, missingIndices, notMissingIndices);

    if (DEBUG) {
        System.err.println("S12:\n" + S12);
    }

    DenseMatrix64F S12S22Inv = new DenseMatrix64F(missingIndices.length, notMissingIndices.length);
    CommonOps.mult(S12, S22Inv, S12S22Inv);

    if (DEBUG) {
        System.err.println("S12S22Inv:\n" + S12S22Inv);
    }

    DenseMatrix64F S12S22InvS21 = new DenseMatrix64F(missingIndices.length, missingIndices.length);
    CommonOps.multTransB(S12S22Inv, S12, S12S22InvS21);

    if (DEBUG) {
        System.err.println("S12S22InvS21:\n" + S12S22InvS21);
    }

    sBar = new DenseMatrix64F(missingIndices.length, missingIndices.length);
    gatherRowsAndColumns(variance, sBar, missingIndices, missingIndices);
    CommonOps.subtract(sBar, S12S22InvS21, sBar);


    if (DEBUG) {
        System.err.println("sBar:\n" + sBar);
    }


    this.affineTransform = S12S22Inv;
    this.tempStorage = new double[missingIndices.length];

    this.numMissing = missingIndices.length;
    this.numNotMissing = notMissingIndices.length;

}
 

@Override
public double[][] getJointVariance(final double priorSampleSize,
                                   final double[][] treeVariance, final double[][] treeSharedLengths,
                                   final double[][] traitVariance) {

    double[] eigVals = this.getEigenValuesStrengthOfSelection();
    DenseMatrix64F V = wrap(this.getEigenVectorsStrengthOfSelection(), 0, dim, dim);
    DenseMatrix64F Vinv = new DenseMatrix64F(dim, dim);
    CommonOps.invert(V, Vinv);

    DenseMatrix64F transTraitVariance = new DenseMatrix64F(traitVariance);

    DenseMatrix64F tmp = new DenseMatrix64F(dim, dim);
    CommonOps.mult(Vinv, transTraitVariance, tmp);
    CommonOps.multTransB(tmp, Vinv, transTraitVariance);

    // Computation of matrix
    int ntaxa = tree.getExternalNodeCount();
    double ti;
    double tj;
    double tij;
    double ep;
    double eq;
    double var;
    DenseMatrix64F varTemp = new DenseMatrix64F(dim, dim);
    double[][] jointVariance = new double[dim * ntaxa][dim * ntaxa];
    for (int i = 0; i < ntaxa; ++i) {
        for (int j = 0; j < ntaxa; ++j) {
            ti = treeSharedLengths[i][i];
            tj = treeSharedLengths[j][j];
            tij = treeSharedLengths[i][j];
            for (int p = 0; p < dim; ++p) {
                for (int q = 0; q < dim; ++q) {
                    ep = eigVals[p];
                    eq = eigVals[q];
                    var = tij / ep / eq;
                    var += (1 - Math.exp(ep * tij)) * Math.exp(-ep * ti) / ep / ep / eq;
                    var += (1 - Math.exp(eq * tij)) * Math.exp(-eq * tj) / ep / eq / eq;
                    var -= (1 - Math.exp((ep + eq) * tij)) * Math.exp(-ep * ti) * Math.exp(-eq * tj) / ep / eq / (ep + eq);
                    var += (1 - Math.exp(-ep * ti)) * (1 - Math.exp(-eq * tj)) / ep / eq / priorSampleSize;
                    var += 1 / priorSampleSize;
                    varTemp.set(p, q, var * transTraitVariance.get(p, q));
                }
            }
            CommonOps.mult(V, varTemp, tmp);
            CommonOps.multTransB(tmp, V, varTemp);
            for (int p = 0; p < dim; ++p) {
                for (int q = 0; q < dim; ++q) {
                    jointVariance[i * dim + p][j * dim + q] = varTemp.get(p, q);
                }
            }
        }
    }
    return jointVariance;
}
 

private double[][] getJointVarianceFull(final double priorSampleSize,
                                        final double[][] treeVariance, final double[][] treeSharedLengths,
                                        final double[][] traitVariance) {

    double[] eigVals = this.getEigenValuesStrengthOfSelection();
    DenseMatrix64F V = wrap(this.getEigenVectorsStrengthOfSelection(), 0, dim, dim);
    DenseMatrix64F Vinv = new DenseMatrix64F(dim, dim);
    CommonOps.invert(V, Vinv);

    DenseMatrix64F transTraitVariance = new DenseMatrix64F(traitVariance);

    DenseMatrix64F tmp = new DenseMatrix64F(dim, dim);
    CommonOps.mult(Vinv, transTraitVariance, tmp);
    CommonOps.multTransB(tmp, Vinv, transTraitVariance);

    // inverse of eigenvalues
    double[][] invEigVals = new double[dim][dim];
    for (int p = 0; p < dim; ++p) {
        for (int q = 0; q < dim; ++q) {
            invEigVals[p][q] = 1 / (eigVals[p] + eigVals[q]);
        }
    }

    // Computation of matrix
    int ntaxa = tree.getExternalNodeCount();
    double ti;
    double tj;
    double tij;
    double ep;
    double eq;
    DenseMatrix64F varTemp = new DenseMatrix64F(dim, dim);
    double[][] jointVariance = new double[dim * ntaxa][dim * ntaxa];
    for (int i = 0; i < ntaxa; ++i) {
        for (int j = 0; j < ntaxa; ++j) {
            ti = treeSharedLengths[i][i];
            tj = treeSharedLengths[j][j];
            tij = treeSharedLengths[i][j];
            for (int p = 0; p < dim; ++p) {
                for (int q = 0; q < dim; ++q) {
                    ep = eigVals[p];
                    eq = eigVals[q];
                    varTemp.set(p, q, Math.exp(-ep * ti) * Math.exp(-eq * tj) * (invEigVals[p][q] * (Math.exp((ep + eq) * tij) - 1) + 1 / priorSampleSize) * transTraitVariance.get(p, q));
                }
            }
            CommonOps.mult(V, varTemp, tmp);
            CommonOps.multTransB(tmp, V, varTemp);
            for (int p = 0; p < dim; ++p) {
                for (int q = 0; q < dim; ++q) {
                    jointVariance[i * dim + p][j * dim + q] = varTemp.get(p, q);
                }
            }
        }
    }
    return jointVariance;
}
 

@Override
public void getBranchVariance(int bufferIndex, int precisionIndex, double[] variance) {
    getBranchPrecision(bufferIndex, precisionIndex, variance);
    DenseMatrix64F Var = DenseMatrix64F.wrap(dimTrait, dimTrait, variance);
    CommonOps.invert(Var);
}
 

private void transformMatrixBaseGeneral(DenseMatrix64F matrix, DenseMatrix64F rotation) {
    DenseMatrix64F tmp = new DenseMatrix64F(dimProcess, dimProcess);
    CommonOps.mult(rotation, matrix, tmp);
    CommonOps.invert(rotation); // Warning: side effect on rotation matrix.
    CommonOps.mult(tmp, rotation, matrix);
}