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

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;

    }
 

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

@Override
void computeOUVarianceBranch(final int sourceOffset,
                             final int destinationOffset,
                             final int destinationOffsetDiagonal,
                             final double edgeLength) {
    DenseMatrix64F actualization = wrap(actualizations, destinationOffset, dimProcess, dimProcess);
    DenseMatrix64F variance = wrap(stationaryVariances, sourceOffset, dimProcess, dimProcess);
    DenseMatrix64F temp = new DenseMatrix64F(dimProcess, dimProcess);

    CommonOps.multTransB(variance, actualization, temp);
    CommonOps.multAdd(-1.0, actualization, temp, variance);

    unwrap(variance, variances, destinationOffset);
}
 

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


    }
 

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

    }
 

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

public static void transformMatrixBack(DenseMatrix64F matrix, DenseMatrix64F rotation) {
    int dim = matrix.getNumRows();
    DenseMatrix64F tmp = new DenseMatrix64F(dim, dim);
    CommonOps.multTransB(matrix, rotation, tmp);
    CommonOps.mult(rotation, tmp, matrix);
}
 

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 scaleVariance(DenseMatrix64F Q, DenseMatrix64F P,
                           DenseMatrix64F QtP, DenseMatrix64F QtPQ) {
    CommonOps.mult(Q, P, QtP);
    CommonOps.multTransB(QtP, Q, QtPQ);
}