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

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

    DenseMatrix64F gradQInvDiag = ((OUDiffusionModelDelegate) diffusionProcessDelegate).getGradientVarianceWrtAttenuation(node, cdi, statistics, gradQInv);

    if (DEBUG) {
        System.err.println("gradQ = " + NormalSufficientStatistics.toVectorizedString(gradQInv));
    }

    DenseMatrix64F gradNDiag = ((OUDiffusionModelDelegate) diffusionProcessDelegate).getGradientDisplacementWrtAttenuation(node, cdi, statistics, gradN);

    CommonOps.addEquals(gradQInvDiag, gradNDiag);

    return gradQInvDiag.getData();
}
 

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

private void computeIOUActualizedDisplacement(final double[] optimalRates,
                                              final int offset,
                                              final int pio,
                                              double branchLength,
                                              DenseMatrix64F inverseSelectionStrength) {
    DenseMatrix64F displacementOU = wrap(displacements, pio, dimProcess, 1);
    DenseMatrix64F optVal = wrap(optimalRates, offset, dimProcess, 1);
    DenseMatrix64F displacement = new DenseMatrix64F(dimProcess, 1);

    CommonOps.mult(inverseSelectionStrength, displacementOU, displacement);
    CommonOps.scale(-1.0, displacement);

    CommonOps.addEquals(displacement, branchLength, optVal);

    unwrap(displacement, displacements, pio + dimProcess);
}
 

private void computeIOUActualization(final int scaledOffset,
                                     DenseMatrix64F inverseSelectionStrength) {
    // YY
    DenseMatrix64F actualizationOU = wrap(actualizations, scaledOffset, dimProcess, dimProcess);

    // XX
    DenseMatrix64F temp = CommonOps.identity(dimProcess);
    CommonOps.addEquals(temp, -1.0, actualizationOU);
    DenseMatrix64F actualizationIOU = new DenseMatrix64F(dimProcess, dimProcess);
    CommonOps.mult(inverseSelectionStrength, temp, actualizationIOU);

    // YX and XX
    DenseMatrix64F actualizationYX = new DenseMatrix64F(dimProcess, dimProcess); // zeros
    DenseMatrix64F actualizationXX = CommonOps.identity(dimProcess);

    blockUnwrap(actualizationOU, actualizationXX, actualizationIOU, actualizationYX, actualizations, scaledOffset);
}
 

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

private DenseMatrix64F getGradientVarianceWrtAttenuationDiagonal(ContinuousDiffusionIntegrator cdi,
                                                                 BranchSufficientStatistics statistics,
                                                                 int nodeIndex, DenseMatrix64F gradient) {
    // wrt to q_i actualization
    DenseMatrix64F gradActualization = getGradientVarianceWrtActualizationDiagonal(cdi, statistics, nodeIndex, gradient);

    // wrt to Gamma stationary variance
    DenseMatrix64F gradStationary = getGradientBranchVarianceWrtAttenuationDiagonal(cdi, nodeIndex, gradient);

    CommonOps.addEquals(gradActualization, gradStationary);
    return gradActualization;

}
 

@Override
void computeOUActualizedDisplacement(final double[] optimalRates,
                                     final int offset,
                                     final int actualizationOffset,
                                     final int pio) {
    DenseMatrix64F actualization = wrap(actualizations, actualizationOffset, dimProcess, dimProcess);
    DenseMatrix64F optVal = wrap(optimalRates, offset, dimProcess, 1);
    DenseMatrix64F temp = CommonOps.identity(dimProcess);
    DenseMatrix64F displacement = new DenseMatrix64F(dimProcess, 1);

    CommonOps.addEquals(temp, -1.0, actualization);
    CommonOps.mult(temp, optVal, displacement);

    unwrap(displacement, displacements, pio);
}
 

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

public static void forceSymmetric(DenseMatrix64F P) {
    DenseMatrix64F Ptrans = new DenseMatrix64F(P);
    CommonOps.transpose(P, Ptrans);
    CommonOps.addEquals(P, Ptrans);
    CommonOps.scale(0.5, P);
}
 

private void computeVariance(DenseMatrix64F matrix, double[] data, int numRows, int numCols) {

        double[] buffer = new double[numRows];
        DenseMatrix64F sumVec = new DenseMatrix64F(numCols, 1);
        DenseMatrix64F matrixBuffer = new DenseMatrix64F(numCols, numCols);

        Arrays.fill(matrix.getData(), 0);

        for (int i = 0; i < numRows; i++) {
            int offset = numCols * i;

            DenseMatrix64F wrapper = MissingOps.wrap(data, offset, numCols, 1, buffer);

            CommonOps.multTransB(wrapper, wrapper, matrixBuffer);
            CommonOps.addEquals(matrix, matrixBuffer);
            CommonOps.addEquals(sumVec, wrapper);

        }

        CommonOps.multTransB(sumVec, sumVec, matrixBuffer);
        CommonOps.addEquals(matrix, -1.0 / numRows, matrixBuffer);
        CommonOps.scale(1.0 / numRows, matrix);


    }
 

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

private DenseMatrix64F getGradientVarianceWrtActualizationDiagonal(ContinuousDiffusionIntegrator cdi, BranchSufficientStatistics statistics,
                                                                       int nodeIndex, DenseMatrix64F gradient) {
        // 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);
        // Q_i^-
        DenseMatrix64F Wi = statistics.getAbove().getRawVarianceCopy();
        // Gamma
//        DenseMatrix64F Gamma = wrap(
//                ((SafeMultivariateDiagonalActualizedWithDriftIntegrator) cdi).getStationaryVariance(getEigenBufferOffsetIndex(0)),
//                0, dim, dim);
        // Branch variance
        double[] branchVariance = new double[dim * dim];
        cdi.getBranchVariance(getMatrixBufferOffsetIndex(nodeIndex), getEigenBufferOffsetIndex(0) , branchVariance);
        DenseMatrix64F Sigma_i = wrap(branchVariance, 0, dim, dim);

        // factor
        DenseMatrix64F res = new DenseMatrix64F(dim, dim);
        CommonOps.addEquals(Wi, -1, Sigma_i);
//        MissingOps.diagDiv(qi, Wi);
//        CommonOps.multTransA(qiInv, tmp, res);
        CommonOps.multTransB(Wi, gradient, res);

        // temp + temp^T
//        MissingOps.addTransEquals(res); No need for diagonal case
        CommonOps.scale(2.0, res);

        // Get diagonal
        DenseMatrix64F resDiag = new DenseMatrix64F(dim, 1);
        CommonOps.extractDiag(res, resDiag);

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

        return resDiag;

    }
 

static void getGradientQInvForBranch(DenseMatrix64F Qi, DenseMatrix64F Wi,
                                     DenseMatrix64F Vi, DenseMatrix64F delta,
                                     DenseMatrix64F grad) {

    CommonOps.scale(0.5, Wi, grad);

    CommonOps.multAddTransB(-0.5, delta, delta, grad);

    CommonOps.addEquals(grad, -0.5, Vi);

    if (DEBUG) {
        System.err.println("\tgradientQi = " + NormalSufficientStatistics.toVectorizedString(grad));
    }

    MultivariateChainRule ruleI = new MultivariateChainRule.InverseGeneral(Qi);
    ruleI.chainGradient(grad);

    if (DEBUG) {
        System.err.println("\tgradientQiInv = " + NormalSufficientStatistics.toVectorizedString(grad));
    }

}
 

private void computeIOUVarianceBranch(final int sourceOffset,
                                      final int destinationOffset,
                                      double branchLength,
                                      DenseMatrix64F inverseSelectionStrength) {
    DenseMatrix64F actualization = wrap(actualizations, destinationOffset, dimProcess, dimProcess);
    DenseMatrix64F stationaryVariance = wrap(stationaryVariances, sourceOffset, dimProcess, dimProcess);

    DenseMatrix64F invAS = new DenseMatrix64F(dimProcess, dimProcess);
    CommonOps.mult(inverseSelectionStrength, stationaryVariance, invAS);

    //// Variance YY
    DenseMatrix64F varianceYY = wrap(variances, destinationOffset, dimProcess, dimProcess);

    //// Variance XX
    DenseMatrix64F varianceXX = new DenseMatrix64F(dimProcess, dimProcess);
    // Variance 1
    CommonOps.multTransB(invAS, inverseSelectionStrength, varianceXX);
    DenseMatrix64F temp = new DenseMatrix64F(dimProcess, dimProcess);
    CommonOps.multTransB(varianceXX, actualization, temp);
    CommonOps.multAdd(-1.0, actualization, temp, varianceXX);
    // Delta
    DenseMatrix64F delta = new DenseMatrix64F(dimProcess, dimProcess);
    addTrans(invAS, delta);
    // Variance 2
    CommonOps.addEquals(varianceXX, branchLength, delta);
    // Variance 3
    DenseMatrix64F temp2 = CommonOps.identity(dimProcess);
    CommonOps.addEquals(temp2, -1.0, actualization);
    DenseMatrix64F temp3 = new DenseMatrix64F(dimProcess, dimProcess);
    CommonOps.mult(temp2, inverseSelectionStrength, temp3);
    CommonOps.mult(temp3, delta, temp2);
    addTrans(temp2, temp);
    // All
    CommonOps.addEquals(varianceXX, -1.0, temp);

    //// Variance XY
    DenseMatrix64F varianceXY = new DenseMatrix64F(dimProcess, dimProcess);
    // Variance 1
    CommonOps.multTransB(stationaryVariance, temp3, varianceXY);
    // Variance 2
    CommonOps.mult(temp3, stationaryVariance, temp);
    CommonOps.multTransB(temp, actualization, temp2);
    // All
    CommonOps.addEquals(varianceXY, -1.0, temp2);

    //// Variance YX
    DenseMatrix64F varianceYX = new DenseMatrix64F(dimProcess, dimProcess);
    CommonOps.transpose(varianceXY, varianceYX);

    blockUnwrap(varianceYY, varianceXX, varianceXY, varianceYX, variances, destinationOffset);
    schurComplementInverse(varianceYY, varianceXX, varianceXY, varianceYX, precisions, destinationOffset);
}
 

private void addTrans(DenseMatrix64F A, DenseMatrix64F B) {
    CommonOps.transpose(A, B);
    CommonOps.addEquals(B, A);
}
 

private void updateVarianceComponents() {

        if (useEmpiricalVariance) {

            if (forceResample) {
                likelihoodDelegate.fireModelChanged(); //Forces new sample
            }
            double[] tipTraits = (double[]) treeTrait.getTrait(treeLikelihood.getTree(), null);
            double[] data = extensionDelegate.getExtendedValues(tipTraits);

            int nTaxa = tree.getExternalNodeCount();

            computeVariance(diffusionComponent, tipTraits, nTaxa, dimTrait);
            computeVariance(samplingComponent, data, nTaxa, dimTrait);

            CommonOps.addEquals(samplingComponent, -1, diffusionComponent);


        } else {
            double n = tree.getExternalNodeCount();

            double diffusionScale = (treeSums.diagonalSum / n - treeSums.totalSum / (n * n));
            double samplingScale = (n - 1) / n;

            for (int i = 0; i < dimTrait; i++) {

                diffusionComponent.set(i, i, diffusionScale * diffusionVariance.component(i, i));
                samplingComponent.set(i, i, samplingScale * samplingVariance.component(i, i));

                for (int j = i + 1; j < dimTrait; j++) {

                    double diffValue = diffusionScale * diffusionVariance.component(i, j);
                    double sampValue = samplingScale * samplingVariance.component(i, j);

                    diffusionComponent.set(i, j, diffValue);
                    samplingComponent.set(i, j, sampValue);

                    diffusionComponent.set(j, i, diffValue);
                    samplingComponent.set(j, i, sampValue);

                }
            }

        }
    }