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

private static NormalSufficientStatistics computeJointLatent(NormalSufficientStatistics below,
                                                             NormalSufficientStatistics above,
                                                             int dim) {

    DenseMatrix64F mean = new DenseMatrix64F(dim, 1);
    DenseMatrix64F precision = new DenseMatrix64F(dim, dim);
    DenseMatrix64F variance = new DenseMatrix64F(dim, dim);

    CommonOps.add(below.getRawPrecision(), above.getRawPrecision(), precision);
    safeInvert2(precision, variance, false);

    safeWeightedAverage(
            new WrappedVector.Raw(below.getRawMean().getData(), 0, dim),
            below.getRawPrecision(),
            new WrappedVector.Raw(above.getRawMean().getData(), 0, dim),
            above.getRawPrecision(),
            new WrappedVector.Raw(mean.getData(), 0, dim),
            variance,
            dim);

    return new NormalSufficientStatistics(mean, precision, variance);
}
 

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

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

@Override
protected void simulateTraitForRoot(final int offsetSample, final int offsetPartial) {

    // TODO Bad programming -- should not need to know about internal layout
    // Layout, offset, dim
    // trait, 0, dT
    // precision, dT, dT * dT
    // variance, dT + dT * dT, dT * dT
    // scalar, dT + 2 * dT * dT, 1

    // Integrate out against prior
    final DenseMatrix64F rootPrec = wrap(partialNodeBuffer, offsetPartial + dimTrait, dimTrait, dimTrait);
    final DenseMatrix64F priorPrec = new DenseMatrix64F(dimTrait, dimTrait);
    MissingOps.safeMult(Pd, wrap(partialPriorBuffer, offsetPartial + dimTrait, dimTrait, dimTrait), priorPrec);

    final DenseMatrix64F totalPrec = new DenseMatrix64F(dimTrait, dimTrait);
    CommonOps.add(rootPrec, priorPrec, totalPrec);

    final DenseMatrix64F totalVar = new DenseMatrix64F(dimTrait, dimTrait);
    safeInvert2(totalPrec, totalVar, false);

    final double[] mean = new double[dimTrait];

    safeWeightedAverage(new WrappedVector.Raw(partialNodeBuffer, offsetPartial, dimTrait),
            rootPrec,
            new WrappedVector.Raw(partialPriorBuffer, offsetPartial, dimTrait),
            priorPrec,
            new WrappedVector.Raw(mean, 0, dimTrait),
            totalVar,
            dimTrait
            );

    final DenseMatrix64F cholesky = getCholeskyOfVariance(totalVar, dimTrait);

    MultivariateNormalDistribution.nextMultivariateNormalCholesky(
            new WrappedVector.Raw(mean, 0, dimTrait), // input mean
            new WrappedMatrix.Raw(cholesky.getData(), 0, dimTrait, dimTrait), 1.0, // input variance
            new WrappedVector.Raw(sample, offsetSample, dimTrait), // output sample
            tmpEpsilon);

    if (DEBUG) {
        System.err.println("Attempt to simulate root");
        System.err.println("Root  mean: " + new WrappedVector.Raw(partialNodeBuffer, offsetPartial, dimTrait));
        System.err.println("Root  prec: " + rootPrec);
        System.err.println("Prior mean: " + new WrappedVector.Raw(partialPriorBuffer, offsetPartial, dimTrait));
        System.err.println("Prior prec: " + priorPrec);
        System.err.println("Total prec: " + totalPrec);
        System.err.println("Total  var: " + totalVar);


        System.err.println("draw mean: " + new WrappedVector.Raw(mean, 0, dimTrait));
        System.err.println("sample: " + new WrappedVector.Raw(sample, offsetSample, dimTrait));
    }
}
 

void getSufficientStatistics(BranchSufficientStatistics statistics) {
    final NormalSufficientStatistics below = statistics.getBelow();
    final NormalSufficientStatistics above = statistics.getAbove();

    // Sampling Parameters: Statistic data model
    DenseMatrix64F samplingVariance = dataModel.getExtensionVariance();

    // One more pre-order step
    // TODO This is just one more pre-order step. Should maybe be moved elsewhere ?
    // TODO Below only works for one data point (no repetition)
    // above data
    DenseMatrix64F aboveDataVariance = new DenseMatrix64F(dim, dim);
    DenseMatrix64F aboveDataPrecision = new DenseMatrix64F(dim, dim);
    CommonOps.add(above.getRawVariance(), samplingVariance, aboveDataVariance);
    safeInvert2(aboveDataVariance, aboveDataPrecision, false);
    final NormalSufficientStatistics aboveData = new NormalSufficientStatistics(
            above.getRawMeanCopy(),
            aboveDataPrecision,
            aboveDataVariance);
    // Below data
    int[] missings = statistics.getMissing();
    DenseMatrix64F precisionData = new DenseMatrix64F(dim, dim);
    for (int i = 0; i < dim; i++) {
        precisionData.unsafe_set(i, i, Double.POSITIVE_INFINITY);
    }
    for (int m : missings) {
        precisionData.unsafe_set(m, m, 0.0);
    }
    final NormalSufficientStatistics belowData = new NormalSufficientStatistics(
            below.getRawMeanCopy(),
            precisionData);

    // Joint data
    NormalSufficientStatistics jointStatisticsData =
            BranchRateGradient.ContinuousTraitGradientForBranch.Default.computeJointStatistics(
                    belowData, aboveData, dim
            );

    // Gradient
    matrixQ = aboveData.getRawPrecision();
    matrixW = aboveData.getRawVariance();
    matrixV = jointStatisticsData.getRawVariance();

    // Delta
    for (int row = 0; row < dim; ++row) {
        matrixDelta.unsafe_set(row, 0,
                jointStatisticsData.getMean(row) - aboveData.getMean(row)
        );
    }
}
 

private void inflateBranch(DenseMatrix64F Vj, DenseMatrix64F Vdj, DenseMatrix64F Vjp) {
    CommonOps.add(Vj, Vdj, Vjp);
}
 

private InversionResult increaseVariances(int ibo,
                                          int iBuffer,
                                          final DenseMatrix64F Vdi,
                                          final DenseMatrix64F Pdi,
                                          final DenseMatrix64F Pip,
                                          final boolean getDeterminant) {

    if (TIMING) {
        startTime("peel1");
    }

    // A. Get current precision of i and j
    final DenseMatrix64F Pi = wrap(partials, ibo + dimTrait, dimTrait, dimTrait);

    if (TIMING) {
        endTime("peel1");
        startTime("peel2");
    }

    // B. Integrate along branch using two matrix inversions

    final boolean useVariancei = anyDiagonalInfinities(Pi);
    InversionResult ci = null;

    if (useVariancei) {

        final DenseMatrix64F Vip = matrix0;
        final DenseMatrix64F Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait);
        CommonOps.add(Vi, Vdi, Vip);
        if (allZeroOrInfinite(Vip)) {
            throw new RuntimeException("Zero-length branch on data is not allowed.");
        }
        ci = safeInvert2(Vip, Pip, getDeterminant);

    } else {

        final DenseMatrix64F tmp1 = matrix0;
        CommonOps.add(Pi, Pdi, tmp1);
        final DenseMatrix64F tmp2 = matrix1;
        safeInvert2(tmp1, tmp2, false);
        CommonOps.mult(tmp2, Pi, tmp1);
        idMinusA(tmp1);
        if (getDeterminant) ci = safeDeterminant(tmp1, true);
        CommonOps.mult(Pi, tmp1, Pip);
        if (getDeterminant && getEffectiveDimension(iBuffer) > 0) {
            InversionResult cP = safeDeterminant(Pi, true);
            ci = mult(ci, cP);
        }
    }

    if (TIMING) {
        endTime("peel2");
    }

    return ci;
}
 

void computePartialPrecision(int ido, int jdo, int imo, int jmo,
                             DenseMatrix64F Pip, DenseMatrix64F Pjp, DenseMatrix64F Pk) {
    CommonOps.add(Pip, Pjp, Pk);
}