Java Code Examples for org.ejml.ops.CommonOps#multTransB()
The following examples show how to use
org.ejml.ops.CommonOps#multTransB() .
You can vote up the ones you like or vote down the ones you don't like,
and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
| Source File: OUDiffusionModelDelegate.java From beast-mcmc with GNU Lesser General Public License v2.1 | 5 votes |
|
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;
}
Example 2
| Source File: SafeMultivariateActualizedWithDriftIntegrator.java From beast-mcmc with GNU Lesser General Public License v2.1 | 5 votes |
|
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);
}
Example 3
| Source File: SafeMultivariateActualizedWithDriftIntegrator.java From beast-mcmc with GNU Lesser General Public License v2.1 | 5 votes |
|
@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);
}
Example 4
| Source File: VarianceProportionStatistic.java From beast-mcmc with GNU Lesser General Public License v2.1 | 4 votes |
|
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);
}
Example 5
| Source File: ConditionalVarianceAndTransform2.java From beast-mcmc with GNU Lesser General Public License v2.1 | 4 votes |
|
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;
}
Example 6
| Source File: IntegratedOUDiffusionModelDelegate.java From beast-mcmc with GNU Lesser General Public License v2.1 | 4 votes |
|
@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;
}
Example 7
| Source File: OUDiffusionModelDelegate.java From beast-mcmc with GNU Lesser General Public License v2.1 | 4 votes |
|
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;
}
Example 8
| Source File: OUDiffusionModelDelegate.java From beast-mcmc with GNU Lesser General Public License v2.1 | 4 votes |
|
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;
}
Example 9
| Source File: SafeMultivariateActualizedWithDriftIntegrator.java From beast-mcmc with GNU Lesser General Public License v2.1 | 4 votes |
|
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);
}
Example 10
| Source File: SafeMultivariateActualizedWithDriftIntegrator.java From beast-mcmc with GNU Lesser General Public License v2.1 | 4 votes |
|
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);
}
Example 11
| Source File: SafeMultivariateActualizedWithDriftIntegrator.java From beast-mcmc with GNU Lesser General Public License v2.1 | 4 votes |
|
private void scaleVariance(DenseMatrix64F Q, DenseMatrix64F P,
DenseMatrix64F QtP, DenseMatrix64F QtPQ) {
CommonOps.mult(Q, P, QtP);
CommonOps.multTransB(QtP, Q, QtPQ);
}
