Author: Andrius Merkys Description: Migrating to libejml-java v0.38. --- a/src/dr/evomodel/treedatalikelihood/continuous/cdi/ContinuousDiffusionIntegrator.java +++ b/src/dr/evomodel/treedatalikelihood/continuous/cdi/ContinuousDiffusionIntegrator.java @@ -413,7 +413,7 @@ final double vi = branchLengths[imo]; final double vj = branchLengths[jmo]; // -// final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); +// final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); // if (DEBUG) { System.err.println("updatePreOrderPartial for node " + iBuffer); @@ -427,33 +427,33 @@ final double pk = prePartials[kbo + dimTrait]; final double pj = partials[jbo + dimTrait]; -// final DenseMatrix64F Pk = wrap(prePartials, kbo + dimTrait, dimTrait, dimTrait); -// // final DenseMatrix64F Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait); +// final DMatrixRMaj Pk = wrap(prePartials, kbo + dimTrait, dimTrait, dimTrait); +// // final DMatrixRMaj Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait); // -// // final DenseMatrix64F Vk = wrap(prePartials, kbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); -// final DenseMatrix64F Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); +// // final DMatrixRMaj Vk = wrap(prePartials, kbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); +// final DMatrixRMaj Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); // // B. Inflate variance along sibling branch using matrix inversion final double pjp = Double.isInfinite(pj) ? 1.0 / vj : pj / (1.0 + pj * vj); -// // final DenseMatrix64F Vjp = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Vjp = matrix0; -// CommonOps.add(Vj, vj, Vd, Vjp); +// // final DMatrixRMaj Vjp = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Vjp = matrix0; +// CommonOps_DDRM.add(Vj, vj, Vd, Vjp); // -// // final DenseMatrix64F Pjp = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Pjp = matrix1; +// // final DMatrixRMaj Pjp = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Pjp = matrix1; // InversionResult cj = safeInvert(Vjp, Pjp, false); // -// // final DenseMatrix64F Pip = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Pip = matrix2; -// CommonOps.add(Pk, Pjp, Pip); +// // final DMatrixRMaj Pip = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Pip = matrix2; +// CommonOps_DDRM.add(Pk, Pjp, Pip); final double pip = pjp + pk; // -// // final DenseMatrix64F Vip = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Vip = matrix3; +// // final DMatrixRMaj Vip = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Vip = matrix3; // InversionResult cip = safeInvert(Pip, Vip, false); // // // C. Compute prePartial mean @@ -484,11 +484,11 @@ final double pi = Double.isInfinite(pip) ? 1.0 / vi : pip / (1.0 + pip * vi); -// final DenseMatrix64F Vi = Vip; -// CommonOps.add(vi, Vd, Vip, Vi); +// final DMatrixRMaj Vi = Vip; +// CommonOps_DDRM.add(vi, Vd, Vip, Vi); // -// // final DenseMatrix64F Pi = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Pi = matrix4; +// // final DMatrixRMaj Pi = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Pi = matrix4; // InversionResult ci = safeInvert(Vi, Pi, false); // // // X. Store precision results for node --- a/src/dr/evomodel/treedatalikelihood/continuous/cdi/MultivariateIntegrator.java +++ b/src/dr/evomodel/treedatalikelihood/continuous/cdi/MultivariateIntegrator.java @@ -2,8 +2,8 @@ import dr.math.matrixAlgebra.WrappedVector; import dr.math.matrixAlgebra.missingData.InversionResult; -import org.ejml.data.DenseMatrix64F; -import org.ejml.ops.CommonOps; +import org.ejml.data.DMatrixRMaj; +import org.ejml.dense.row.CommonOps_DDRM; import java.util.HashMap; import java.util.Map; @@ -56,13 +56,13 @@ private final Map times; - DenseMatrix64F matrix0; - DenseMatrix64F matrix1; - DenseMatrix64F matrix2; - DenseMatrix64F matrix3; - DenseMatrix64F matrix4; - DenseMatrix64F matrix5; - DenseMatrix64F matrix6; + DMatrixRMaj matrix0; + DMatrixRMaj matrix1; + DMatrixRMaj matrix2; + DMatrixRMaj matrix3; + DMatrixRMaj matrix4; + DMatrixRMaj matrix5; + DMatrixRMaj matrix6; double[] vector0; @@ -70,13 +70,13 @@ inverseDiffusions = new double[dimTrait * dimTrait * diffusionCount]; vector0 = new double[dimTrait]; - matrix0 = new DenseMatrix64F(dimTrait, dimTrait); - matrix1 = new DenseMatrix64F(dimTrait, dimTrait); - matrix2 = new DenseMatrix64F(dimTrait, dimTrait); - matrix3 = new DenseMatrix64F(dimTrait, dimTrait); - matrix4 = new DenseMatrix64F(dimTrait, dimTrait); - matrix5 = new DenseMatrix64F(dimTrait, dimTrait); - matrix6 = new DenseMatrix64F(dimTrait, dimTrait); + matrix0 = new DMatrixRMaj(dimTrait, dimTrait); + matrix1 = new DMatrixRMaj(dimTrait, dimTrait); + matrix2 = new DMatrixRMaj(dimTrait, dimTrait); + matrix3 = new DMatrixRMaj(dimTrait, dimTrait); + matrix4 = new DMatrixRMaj(dimTrait, dimTrait); + matrix5 = new DMatrixRMaj(dimTrait, dimTrait); + matrix6 = new DMatrixRMaj(dimTrait, dimTrait); } @Override @@ -86,9 +86,9 @@ assert (inverseDiffusions != null); final int offset = dimTrait * dimTrait * precisionIndex; - DenseMatrix64F precision = wrap(diffusions, offset, dimTrait, dimTrait); - DenseMatrix64F variance = new DenseMatrix64F(dimTrait, dimTrait); - CommonOps.invert(precision, variance); + DMatrixRMaj precision = wrap(diffusions, offset, dimTrait, dimTrait); + DMatrixRMaj variance = new DMatrixRMaj(dimTrait, dimTrait); + CommonOps_DDRM.invert(precision, variance); unwrap(variance, inverseDiffusions, offset); if (DEBUG) { @@ -124,7 +124,7 @@ final double vi = branchLengths[imo]; final double vj = branchLengths[jmo]; - final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); + final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); if (DEBUG) { System.err.println("updatePreOrderPartial for node " + iBuffer); @@ -137,27 +137,27 @@ for (int trait = 0; trait < numTraits; ++trait) { // A. Get current precision of k and j - final DenseMatrix64F Pk = wrap(prePartials, kbo + dimTrait, dimTrait, dimTrait); -// final DenseMatrix64F Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Pk = wrap(prePartials, kbo + dimTrait, dimTrait, dimTrait); +// final DMatrixRMaj Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait); -// final DenseMatrix64F Vk = wrap(prePartials, kbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); - final DenseMatrix64F Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); +// final DMatrixRMaj Vk = wrap(prePartials, kbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); // B. Inflate variance along sibling branch using matrix inversion -// final DenseMatrix64F Vjp = new DenseMatrix64F(dimTrait, dimTrait); - final DenseMatrix64F Vjp = matrix0; - CommonOps.add(Vj, vj, Vd, Vjp); +// final DMatrixRMaj Vjp = new DMatrixRMaj(dimTrait, dimTrait); + final DMatrixRMaj Vjp = matrix0; + CommonOps_DDRM.add(Vj, vj, Vd, Vjp); -// final DenseMatrix64F Pjp = new DenseMatrix64F(dimTrait, dimTrait); - final DenseMatrix64F Pjp = matrix1; +// final DMatrixRMaj Pjp = new DMatrixRMaj(dimTrait, dimTrait); + final DMatrixRMaj Pjp = matrix1; InversionResult cj = safeInvert(Vjp, Pjp, false); -// final DenseMatrix64F Pip = new DenseMatrix64F(dimTrait, dimTrait); - final DenseMatrix64F Pip = matrix2; - CommonOps.add(Pk, Pjp, Pip); +// final DMatrixRMaj Pip = new DMatrixRMaj(dimTrait, dimTrait); + final DMatrixRMaj Pip = matrix2; + CommonOps_DDRM.add(Pk, Pjp, Pip); -// final DenseMatrix64F Vip = new DenseMatrix64F(dimTrait, dimTrait); - final DenseMatrix64F Vip = matrix3; +// final DMatrixRMaj Vip = new DMatrixRMaj(dimTrait, dimTrait); + final DMatrixRMaj Vip = matrix3; InversionResult cip = safeInvert(Pip, Vip, false); // C. Compute prePartial mean @@ -180,11 +180,11 @@ } // C. Inflate variance along node branch - final DenseMatrix64F Vi = Vip; - CommonOps.add(vi, Vd, Vip, Vi); + final DMatrixRMaj Vi = Vip; + CommonOps_DDRM.add(vi, Vd, Vip, Vi); -// final DenseMatrix64F Pi = new DenseMatrix64F(dimTrait, dimTrait); - final DenseMatrix64F Pi = matrix4; +// final DMatrixRMaj Pi = new DMatrixRMaj(dimTrait, dimTrait); + final DMatrixRMaj Pi = matrix4; InversionResult ci = safeInvert(Vi, Pi, false); // X. Store precision results for node @@ -242,7 +242,7 @@ final double vi = branchLengths[imo]; final double vj = branchLengths[jmo]; - final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); + final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); if (DEBUG) { System.err.println("variance diffusion: " + Vd); @@ -269,11 +269,11 @@ final double lpi = partials[ibo + dimTrait + 2 * dimTrait * dimTrait]; final double lpj = partials[jbo + dimTrait + 2 * dimTrait * dimTrait]; - final DenseMatrix64F Pi = wrap(partials, ibo + dimTrait, dimTrait, dimTrait); - final DenseMatrix64F Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Pi = wrap(partials, ibo + dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait); - final DenseMatrix64F Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); - final DenseMatrix64F Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); if (TIMING) { endTime("peel1"); @@ -286,23 +286,23 @@ final double lpjp = Double.isInfinite(lpj) ? 1.0 / vj : lpj / (1.0 + lpj * vj); -// final DenseMatrix64F Vip = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Vjp = new DenseMatrix64F(dimTrait, dimTrait); - final DenseMatrix64F Vip = matrix0; - final DenseMatrix64F Vjp = matrix1; +// final DMatrixRMaj Vip = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Vjp = new DMatrixRMaj(dimTrait, dimTrait); + final DMatrixRMaj Vip = matrix0; + final DMatrixRMaj Vjp = matrix1; - CommonOps.add(Vi, vi, Vd, Vip); - CommonOps.add(Vj, vj, Vd, Vjp); + CommonOps_DDRM.add(Vi, vi, Vd, Vip); + CommonOps_DDRM.add(Vj, vj, Vd, Vjp); if (TIMING) { endTime("peel2"); startTime("peel2a"); } -// final DenseMatrix64F Pip = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Pjp = new DenseMatrix64F(dimTrait, dimTrait); - final DenseMatrix64F Pip = matrix2; - final DenseMatrix64F Pjp = matrix3; +// final DMatrixRMaj Pip = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Pjp = new DMatrixRMaj(dimTrait, dimTrait); + final DMatrixRMaj Pip = matrix2; + final DMatrixRMaj Pjp = matrix3; InversionResult ci = safeInvert(Vip, Pip, true); InversionResult cj = safeInvert(Vjp, Pjp, true); @@ -317,13 +317,13 @@ // A. Partial precision and variance (for later use) using one matrix inversion final double lpk = lpip + lpjp; -// final DenseMatrix64F Pk = new DenseMatrix64F(dimTrait, dimTrait); - final DenseMatrix64F Pk = matrix4; +// final DMatrixRMaj Pk = new DMatrixRMaj(dimTrait, dimTrait); + final DMatrixRMaj Pk = matrix4; - CommonOps.add(Pip, Pjp, Pk); + CommonOps_DDRM.add(Pip, Pjp, Pk); -// final DenseMatrix64F Vk = new DenseMatrix64F(dimTrait, dimTrait); - final DenseMatrix64F Vk = matrix5; +// final DMatrixRMaj Vk = new DMatrixRMaj(dimTrait, dimTrait); + final DMatrixRMaj Vk = matrix5; InversionResult ck = safeInvert(Pk, Vk, true); // B. Partial mean @@ -433,9 +433,9 @@ } } -// final DenseMatrix64F Vt = new DenseMatrix64F(dimTrait, dimTrait); - final DenseMatrix64F Vt = matrix6; - CommonOps.add(Vip, Vjp, Vt); +// final DMatrixRMaj Vt = new DMatrixRMaj(dimTrait, dimTrait); + final DMatrixRMaj Vt = matrix6; + CommonOps_DDRM.add(Vip, Vjp, Vt); if (DEBUG) { System.err.println("Vt: " + Vt); @@ -445,7 +445,7 @@ - ck.getEffectiveDimension(); remainder += -dimensionChange * LOG_SQRT_2_PI - 0.5 * -// (Math.log(CommonOps.det(Vip)) + Math.log(CommonOps.det(Vjp)) - Math.log(CommonOps.det(Vk))) +// (Math.log(CommonOps_DDRM.det(Vip)) + Math.log(CommonOps_DDRM.det(Vjp)) - Math.log(CommonOps_DDRM.det(Vk))) (Math.log(ci.getDeterminant()) + Math.log(cj.getDeterminant()) + Math.log(ck.getDeterminant())) - 0.5 * (SSi + SSj - SSk); @@ -469,7 +469,7 @@ assert (false); - final DenseMatrix64F Pt = new DenseMatrix64F(dimTrait, dimTrait); + final DMatrixRMaj Pt = new DMatrixRMaj(dimTrait, dimTrait); InversionResult ct = safeInvert(Vt, Pt, false); int opo = dimTrait * dimTrait * trait; @@ -551,29 +551,29 @@ int rootOffset = dimPartial * rootBufferIndex; int priorOffset = dimPartial * priorBufferIndex; - final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); + final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); // TODO For each trait in parallel for (int trait = 0; trait < numTraits; ++trait) { - final DenseMatrix64F Proot = wrap(partials, rootOffset + dimTrait, dimTrait, dimTrait); - final DenseMatrix64F Pprior = wrap(partials, priorOffset + dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Proot = wrap(partials, rootOffset + dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Pprior = wrap(partials, priorOffset + dimTrait, dimTrait, dimTrait); - final DenseMatrix64F Vroot = wrap(partials, rootOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); - final DenseMatrix64F Vprior = wrap(partials, priorOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Vroot = wrap(partials, rootOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Vprior = wrap(partials, priorOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); // TODO Block below is for the conjugate prior ONLY { - final DenseMatrix64F Vtmp = new DenseMatrix64F(dimTrait, dimTrait); - CommonOps.mult(Vd, Vprior, Vtmp); + final DMatrixRMaj Vtmp = new DMatrixRMaj(dimTrait, dimTrait); + CommonOps_DDRM.mult(Vd, Vprior, Vtmp); Vprior.set(Vtmp); } - final DenseMatrix64F Vtotal = new DenseMatrix64F(dimTrait, dimTrait); - CommonOps.add(Vroot, Vprior, Vtotal); + final DMatrixRMaj Vtotal = new DMatrixRMaj(dimTrait, dimTrait); + CommonOps_DDRM.add(Vroot, Vprior, Vtotal); - final DenseMatrix64F Ptotal = new DenseMatrix64F(dimTrait, dimTrait); - CommonOps.invert(Vtotal, Ptotal); // TODO Can return determinant at same time to avoid extra QR decomp + final DMatrixRMaj Ptotal = new DMatrixRMaj(dimTrait, dimTrait); + CommonOps_DDRM.invert(Vtotal, Ptotal); // TODO Can return determinant at same time to avoid extra QR decomp double SS = 0; for (int g = 0; g < dimTrait; ++g) { @@ -586,7 +586,7 @@ } } - final double logLike = -dimTrait * LOG_SQRT_2_PI - 0.5 * Math.log(CommonOps.det(Vtotal)) - 0.5 * SS; + final double logLike = -dimTrait * LOG_SQRT_2_PI - 0.5 * Math.log(CommonOps_DDRM.det(Vtotal)) - 0.5 * SS; final double remainder = remainders[rootBufferIndex * numTraits + trait]; logLikelihoods[trait] = logLike + remainder; --- a/src/dr/evomodel/treedatalikelihood/continuous/cdi/SafeMultivariateIntegrator.java +++ b/src/dr/evomodel/treedatalikelihood/continuous/cdi/SafeMultivariateIntegrator.java @@ -2,8 +2,8 @@ import dr.math.matrixAlgebra.WrappedVector; import dr.math.matrixAlgebra.missingData.InversionResult; -import org.ejml.data.DenseMatrix64F; -import org.ejml.ops.CommonOps; +import org.ejml.data.DMatrixRMaj; +import org.ejml.dense.row.CommonOps_DDRM; import static dr.math.matrixAlgebra.missingData.InversionResult.Code.NOT_OBSERVED; import static dr.math.matrixAlgebra.missingData.MissingOps.*; @@ -53,13 +53,13 @@ // private final Map times; // -// private DenseMatrix64F matrix0; -// private DenseMatrix64F matrix1; -// private DenseMatrix64F matrix2; -// private DenseMatrix64F matrix3; -// private DenseMatrix64F matrix4; -// private DenseMatrix64F matrix5; -// private DenseMatrix64F matrix6; +// private DMatrixRMaj matrix0; +// private DMatrixRMaj matrix1; +// private DMatrixRMaj matrix2; +// private DMatrixRMaj matrix3; +// private DMatrixRMaj matrix4; +// private DMatrixRMaj matrix5; +// private DMatrixRMaj matrix6; // // private double[] vector0; @@ -67,13 +67,13 @@ // inverseDiffusions = new double[dimTrait * dimTrait * diffusionCount]; // // vector0 = new double[dimTrait]; -// matrix0 = new DenseMatrix64F(dimTrait, dimTrait); -// matrix1 = new DenseMatrix64F(dimTrait, dimTrait); -// matrix2 = new DenseMatrix64F(dimTrait, dimTrait); -// matrix3 = new DenseMatrix64F(dimTrait, dimTrait); -// matrix4 = new DenseMatrix64F(dimTrait, dimTrait); -// matrix5 = new DenseMatrix64F(dimTrait, dimTrait); -// matrix6 = new DenseMatrix64F(dimTrait, dimTrait); +// matrix0 = new DMatrixRMaj(dimTrait, dimTrait); +// matrix1 = new DMatrixRMaj(dimTrait, dimTrait); +// matrix2 = new DMatrixRMaj(dimTrait, dimTrait); +// matrix3 = new DMatrixRMaj(dimTrait, dimTrait); +// matrix4 = new DMatrixRMaj(dimTrait, dimTrait); +// matrix5 = new DMatrixRMaj(dimTrait, dimTrait); +// matrix6 = new DMatrixRMaj(dimTrait, dimTrait); // } // @Override @@ -83,9 +83,9 @@ // assert (inverseDiffusions != null); // // final int offset = dimTrait * dimTrait * precisionIndex; -// DenseMatrix64F precision = wrap(diffusions, offset, dimTrait, dimTrait); -// DenseMatrix64F variance = new DenseMatrix64F(dimTrait, dimTrait); -// CommonOps.invert(precision, variance); +// DMatrixRMaj precision = wrap(diffusions, offset, dimTrait, dimTrait); +// DMatrixRMaj variance = new DMatrixRMaj(dimTrait, dimTrait); +// CommonOps_DDRM.invert(precision, variance); // unwrap(variance, inverseDiffusions, offset); // // if (DEBUG) { @@ -121,7 +121,7 @@ // final double vi = branchLengths[imo]; // final double vj = branchLengths[jmo]; // -// final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); +// final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); // // if (DEBUG) { // System.err.println("updatePreOrderPartial for node " + iBuffer); @@ -134,27 +134,27 @@ // for (int trait = 0; trait < numTraits; ++trait) { // // // A. Get current precision of k and j -// final DenseMatrix64F Pk = wrap(prePartials, kbo + dimTrait, dimTrait, dimTrait); -//// final DenseMatrix64F Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait); +// final DMatrixRMaj Pk = wrap(prePartials, kbo + dimTrait, dimTrait, dimTrait); +//// final DMatrixRMaj Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait); // -//// final DenseMatrix64F Vk = wrap(prePartials, kbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); -// final DenseMatrix64F Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); +//// final DMatrixRMaj Vk = wrap(prePartials, kbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); +// final DMatrixRMaj Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); // // // B. Inflate variance along sibling branch using matrix inversion -//// final DenseMatrix64F Vjp = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Vjp = matrix0; -// CommonOps.add(Vj, vj, Vd, Vjp); +//// final DMatrixRMaj Vjp = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Vjp = matrix0; +// CommonOps_DDRM.add(Vj, vj, Vd, Vjp); // -//// final DenseMatrix64F Pjp = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Pjp = matrix1; +//// final DMatrixRMaj Pjp = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Pjp = matrix1; // InversionResult cj = safeInvert(Vjp, Pjp, false); // -//// final DenseMatrix64F Pip = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Pip = matrix2; -// CommonOps.add(Pk, Pjp, Pip); +//// final DMatrixRMaj Pip = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Pip = matrix2; +// CommonOps_DDRM.add(Pk, Pjp, Pip); // -//// final DenseMatrix64F Vip = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Vip = matrix3; +//// final DMatrixRMaj Vip = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Vip = matrix3; // InversionResult cip = safeInvert(Pip, Vip, false); // // // C. Compute prePartial mean @@ -177,11 +177,11 @@ // } // // // C. Inflate variance along node branch -// final DenseMatrix64F Vi = Vip; -// CommonOps.add(vi, Vd, Vip, Vi); +// final DMatrixRMaj Vi = Vip; +// CommonOps_DDRM.add(vi, Vd, Vip, Vi); // -//// final DenseMatrix64F Pi = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Pi = matrix4; +//// final DMatrixRMaj Pi = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Pi = matrix4; // InversionResult ci = safeInvert(Vi, Pi, false); // // // X. Store precision results for node @@ -239,8 +239,8 @@ final double vi = branchLengths[imo]; final double vj = branchLengths[jmo]; - final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); - final DenseMatrix64F Pd = wrap(diffusions, precisionOffset, dimTrait, dimTrait); + final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); + final DMatrixRMaj Pd = wrap(diffusions, precisionOffset, dimTrait, dimTrait); if (DEBUG) { System.err.println("variance diffusion: " + Vd); @@ -267,8 +267,8 @@ final double lpi = partials[ibo + dimTrait + 2 * dimTrait * dimTrait]; final double lpj = partials[jbo + dimTrait + 2 * dimTrait * dimTrait]; - final DenseMatrix64F Pi = wrap(partials, ibo + dimTrait, dimTrait, dimTrait); - final DenseMatrix64F Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Pi = wrap(partials, ibo + dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait); if (TIMING) { endTime("peel1"); @@ -284,8 +284,8 @@ InversionResult ci; InversionResult cj; - final DenseMatrix64F Pip = matrix2; - final DenseMatrix64F Pjp = matrix3; + final DMatrixRMaj Pip = matrix2; + final DMatrixRMaj Pjp = matrix3; // boolean useVariance = anyDiagonalInfinities(Pi) || anyDiagonalInfinities(Pj); final boolean useVariancei = anyDiagonalInfinities(Pi); @@ -293,77 +293,77 @@ if (useVariancei) { - final DenseMatrix64F Vip = matrix0; - final DenseMatrix64F Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); - CommonOps.add(Vi, vi, Vd, Vip); + final DMatrixRMaj Vip = matrix0; + final DMatrixRMaj Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); + CommonOps_DDRM.add(Vi, vi, Vd, Vip); ci = safeInvert(Vip, Pip, true); } else { - final DenseMatrix64F PiPlusPd = matrix0; - CommonOps.add(Pi, 1.0 / vi, Pd, PiPlusPd); - final DenseMatrix64F PiPlusPdInv = new DenseMatrix64F(dimTrait, dimTrait); + final DMatrixRMaj PiPlusPd = matrix0; + CommonOps_DDRM.add(Pi, 1.0 / vi, Pd, PiPlusPd); + final DMatrixRMaj PiPlusPdInv = new DMatrixRMaj(dimTrait, dimTrait); safeInvert(PiPlusPd, PiPlusPdInv, false); - CommonOps.mult(PiPlusPdInv, Pi, Pip); - CommonOps.mult(Pi, Pip, PiPlusPdInv); - CommonOps.add(Pi, -1, PiPlusPdInv, Pip); + CommonOps_DDRM.mult(PiPlusPdInv, Pi, Pip); + CommonOps_DDRM.mult(Pi, Pip, PiPlusPdInv); + CommonOps_DDRM.add(Pi, -1, PiPlusPdInv, Pip); ci = safeDeterminant(Pip, false); } if (useVariancej) { - final DenseMatrix64F Vjp = matrix1; - final DenseMatrix64F Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); - CommonOps.add(Vj, vj, Vd, Vjp); + final DMatrixRMaj Vjp = matrix1; + final DMatrixRMaj Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); + CommonOps_DDRM.add(Vj, vj, Vd, Vjp); cj = safeInvert(Vjp, Pjp, true); } else { - final DenseMatrix64F PjPlusPd = matrix1; - CommonOps.add(Pj, 1.0 / vj, Pd, PjPlusPd); - final DenseMatrix64F PjPlusPdInv = new DenseMatrix64F(dimTrait, dimTrait); + final DMatrixRMaj PjPlusPd = matrix1; + CommonOps_DDRM.add(Pj, 1.0 / vj, Pd, PjPlusPd); + final DMatrixRMaj PjPlusPdInv = new DMatrixRMaj(dimTrait, dimTrait); safeInvert(PjPlusPd, PjPlusPdInv, false); - CommonOps.mult(PjPlusPdInv, Pj, Pjp); - CommonOps.mult(Pj, Pjp, PjPlusPdInv); - CommonOps.add(Pj, -1, PjPlusPdInv, Pjp); + CommonOps_DDRM.mult(PjPlusPdInv, Pj, Pjp); + CommonOps_DDRM.mult(Pj, Pjp, PjPlusPdInv); + CommonOps_DDRM.add(Pj, -1, PjPlusPdInv, Pjp); cj = safeDeterminant(Pjp, false); } // if (useVariance) { // -// final DenseMatrix64F Vip = matrix0; -// final DenseMatrix64F Vjp = matrix1; +// final DMatrixRMaj Vip = matrix0; +// final DMatrixRMaj Vjp = matrix1; // -// final DenseMatrix64F Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); -// final DenseMatrix64F Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); +// final DMatrixRMaj Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); +// final DMatrixRMaj Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); // -// CommonOps.add(Vi, vi, Vd, Vip); -// CommonOps.add(Vj, vj, Vd, Vjp); +// CommonOps_DDRM.add(Vi, vi, Vd, Vip); +// CommonOps_DDRM.add(Vj, vj, Vd, Vjp); // // ci = safeInvert(Vip, Pip, true); // cj = safeInvert(Vjp, Pjp, true); // } else { // -// final DenseMatrix64F PiPlusPd = matrix0; -// final DenseMatrix64F PjPlusPd = matrix1; +// final DMatrixRMaj PiPlusPd = matrix0; +// final DMatrixRMaj PjPlusPd = matrix1; // -// CommonOps.add(Pi, 1.0 / vi, Pd, PiPlusPd); -// CommonOps.add(Pj, 1.0 / vj, Pd, PjPlusPd); +// CommonOps_DDRM.add(Pi, 1.0 / vi, Pd, PiPlusPd); +// CommonOps_DDRM.add(Pj, 1.0 / vj, Pd, PjPlusPd); // -// final DenseMatrix64F PiPlusPdInv = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F PjPlusPdInv = new DenseMatrix64F(dimTrait, dimTrait); +// final DMatrixRMaj PiPlusPdInv = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj PjPlusPdInv = new DMatrixRMaj(dimTrait, dimTrait); // // safeInvert(PiPlusPd, PiPlusPdInv, false); // safeInvert(PjPlusPd, PjPlusPdInv, false); // -// CommonOps.mult(PiPlusPdInv, Pi, Pip); -// CommonOps.mult(PjPlusPdInv, Pj, Pjp); +// CommonOps_DDRM.mult(PiPlusPdInv, Pi, Pip); +// CommonOps_DDRM.mult(PjPlusPdInv, Pj, Pjp); // -// CommonOps.mult(Pi, Pip, PiPlusPdInv); -// CommonOps.mult(Pj, Pjp, PjPlusPdInv); +// CommonOps_DDRM.mult(Pi, Pip, PiPlusPdInv); +// CommonOps_DDRM.mult(Pj, Pjp, PjPlusPdInv); // -// CommonOps.add(Pi, -1, PiPlusPdInv, Pip); -// CommonOps.add(Pj, -1, PjPlusPdInv, Pjp); +// CommonOps_DDRM.add(Pi, -1, PiPlusPdInv, Pip); +// CommonOps_DDRM.add(Pj, -1, PjPlusPdInv, Pjp); // // ci = safeDeterminant(Pip, false); // cj = safeDeterminant(Pjp, false); @@ -385,12 +385,12 @@ // A. Partial precision and variance (for later use) using one matrix inversion final double lpk = lpip + lpjp; -// final DenseMatrix64F Pk = new DenseMatrix64F(dimTrait, dimTrait); - final DenseMatrix64F Pk = matrix4; - CommonOps.add(Pip, Pjp, Pk); +// final DMatrixRMaj Pk = new DMatrixRMaj(dimTrait, dimTrait); + final DMatrixRMaj Pk = matrix4; + CommonOps_DDRM.add(Pip, Pjp, Pk); -// final DenseMatrix64F Vk = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Vk = matrix5; +// final DMatrixRMaj Vk = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Vk = matrix5; // if (useVariance) { //// InversionResult ck = @@ -524,9 +524,9 @@ } } -// final DenseMatrix64F Vt = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Vt = matrix6; -// CommonOps.add(Vip, Vjp, Vt); +// final DMatrixRMaj Vt = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Vt = matrix6; +// CommonOps_DDRM.add(Vip, Vjp, Vt); // if (DEBUG) { // System.err.println("Vt: " + Vt); @@ -536,16 +536,16 @@ - ck.getEffectiveDimension(); // System.err.println(ci.getDeterminant()); -// System.err.println(CommonOps.det(Vip)); +// System.err.println(CommonOps_DDRM.det(Vip)); // // System.err.println(cj.getDeterminant()); -// System.err.println(CommonOps.det(Vjp)); +// System.err.println(CommonOps_DDRM.det(Vjp)); // // System.err.println(1.0 / ck.getDeterminant()); -// System.err.println(CommonOps.det(Vk)); +// System.err.println(CommonOps_DDRM.det(Vk)); remainder += -dimensionChange * LOG_SQRT_2_PI - 0.5 * -// (Math.log(CommonOps.det(Vip)) + Math.log(CommonOps.det(Vjp)) - Math.log(CommonOps.det(Vk))) +// (Math.log(CommonOps_DDRM.det(Vip)) + Math.log(CommonOps_DDRM.det(Vjp)) - Math.log(CommonOps_DDRM.det(Vk))) (Math.log(ci.getDeterminant()) + Math.log(cj.getDeterminant()) + Math.log(ck.getDeterminant())) - 0.5 * (SSi + SSj - SSk); @@ -632,42 +632,42 @@ int rootOffset = dimPartial * rootBufferIndex; int priorOffset = dimPartial * priorBufferIndex; - final DenseMatrix64F Pd = wrap(diffusions, precisionOffset, dimTrait, dimTrait); -// final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); + final DMatrixRMaj Pd = wrap(diffusions, precisionOffset, dimTrait, dimTrait); +// final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); // TODO For each trait in parallel for (int trait = 0; trait < numTraits; ++trait) { - final DenseMatrix64F Proot = wrap(partials, rootOffset + dimTrait, dimTrait, dimTrait); - final DenseMatrix64F Pprior = wrap(partials, priorOffset + dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Proot = wrap(partials, rootOffset + dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Pprior = wrap(partials, priorOffset + dimTrait, dimTrait, dimTrait); -// final DenseMatrix64F Vroot = wrap(partials, rootOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); -// final DenseMatrix64F Vprior = wrap(partials, priorOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); +// final DMatrixRMaj Vroot = wrap(partials, rootOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); +// final DMatrixRMaj Vprior = wrap(partials, priorOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); // TODO Block below is for the conjugate prior ONLY { -// final DenseMatrix64F Vtmp = new DenseMatrix64F(dimTrait, dimTrait); -// CommonOps.mult(Vd, Vprior, Vtmp); +// final DMatrixRMaj Vtmp = new DMatrixRMaj(dimTrait, dimTrait); +// CommonOps_DDRM.mult(Vd, Vprior, Vtmp); // Vprior.set(Vtmp); - final DenseMatrix64F Ptmp = new DenseMatrix64F(dimTrait, dimTrait); - CommonOps.mult(Pd, Pprior, Ptmp); + final DMatrixRMaj Ptmp = new DMatrixRMaj(dimTrait, dimTrait); + CommonOps_DDRM.mult(Pd, Pprior, Ptmp); Pprior.set(Ptmp); // TODO What does this do? } - final DenseMatrix64F Vtotal = new DenseMatrix64F(dimTrait, dimTrait); -// CommonOps.add(Vroot, Vprior, Vtotal); + final DMatrixRMaj Vtotal = new DMatrixRMaj(dimTrait, dimTrait); +// CommonOps_DDRM.add(Vroot, Vprior, Vtotal); - final DenseMatrix64F Ptotal = new DenseMatrix64F(dimTrait, dimTrait); - CommonOps.invert(Vtotal, Ptotal); // TODO Can return determinant at same time to avoid extra QR decomp + final DMatrixRMaj Ptotal = new DMatrixRMaj(dimTrait, dimTrait); + CommonOps_DDRM.invert(Vtotal, Ptotal); // TODO Can return determinant at same time to avoid extra QR decomp - final DenseMatrix64F tmp1 = new DenseMatrix64F(dimTrait, dimTrait); - final DenseMatrix64F tmp2 = new DenseMatrix64F(dimTrait, dimTrait); - CommonOps.add(Proot, Pprior, Ptotal); - CommonOps.invert(Ptotal, Vtotal); - CommonOps.mult(Vtotal, Proot, tmp1); - CommonOps.mult(Proot, tmp1, tmp2); - CommonOps.add(Proot, -1.0, tmp2, Ptotal); + final DMatrixRMaj tmp1 = new DMatrixRMaj(dimTrait, dimTrait); + final DMatrixRMaj tmp2 = new DMatrixRMaj(dimTrait, dimTrait); + CommonOps_DDRM.add(Proot, Pprior, Ptotal); + CommonOps_DDRM.invert(Ptotal, Vtotal); + CommonOps_DDRM.mult(Vtotal, Proot, tmp1); + CommonOps_DDRM.mult(Proot, tmp1, tmp2); + CommonOps_DDRM.add(Proot, -1.0, tmp2, Ptotal); double SS = 0; for (int g = 0; g < dimTrait; ++g) { @@ -681,8 +681,8 @@ } final double logLike = -dimTrait * LOG_SQRT_2_PI -// - 0.5 * Math.log(CommonOps.det(Vtotal)) - + 0.5 * Math.log(CommonOps.det(Ptotal)) +// - 0.5 * Math.log(CommonOps_DDRM.det(Vtotal)) + + 0.5 * Math.log(CommonOps_DDRM.det(Ptotal)) - 0.5 * SS; final double remainder = remainders[rootBufferIndex * numTraits + trait]; --- a/src/dr/evomodel/treedatalikelihood/continuous/cdi/SafeMultivariateWithDriftIntegrator.java +++ b/src/dr/evomodel/treedatalikelihood/continuous/cdi/SafeMultivariateWithDriftIntegrator.java @@ -2,8 +2,8 @@ import dr.math.matrixAlgebra.WrappedVector; import dr.math.matrixAlgebra.missingData.InversionResult; -import org.ejml.data.DenseMatrix64F; -import org.ejml.ops.CommonOps; +import org.ejml.data.DMatrixRMaj; +import org.ejml.dense.row.CommonOps_DDRM; import static dr.math.matrixAlgebra.missingData.InversionResult.Code.NOT_OBSERVED; import static dr.math.matrixAlgebra.missingData.MissingOps.*; @@ -65,9 +65,9 @@ assert (inverseDiffusions != null); final int offset = dimTrait * dimTrait * precisionIndex; - DenseMatrix64F precision = wrap(diffusions, offset, dimTrait, dimTrait); - DenseMatrix64F variance = new DenseMatrix64F(dimTrait, dimTrait); - CommonOps.invert(precision, variance); + DMatrixRMaj precision = wrap(diffusions, offset, dimTrait, dimTrait); + DMatrixRMaj variance = new DMatrixRMaj(dimTrait, dimTrait); + CommonOps_DDRM.invert(precision, variance); unwrap(variance, inverseDiffusions, offset); if (DEBUG) { @@ -199,14 +199,14 @@ // final double vi = variances[imo]; // final double vj = variances[jmo]; - final DenseMatrix64F Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); -// final DenseMatrix64F Pd = wrap(diffusions, precisionOffset, dimTrait, dimTrait); + final DMatrixRMaj Vd = wrap(inverseDiffusions, precisionOffset, dimTrait, dimTrait); +// final DMatrixRMaj Pd = wrap(diffusions, precisionOffset, dimTrait, dimTrait); - final DenseMatrix64F Vdi = wrap(variances, imo, dimTrait, dimTrait); - final DenseMatrix64F Vdj = wrap(variances, jmo, dimTrait, dimTrait); + final DMatrixRMaj Vdi = wrap(variances, imo, dimTrait, dimTrait); + final DMatrixRMaj Vdj = wrap(variances, jmo, dimTrait, dimTrait); - final DenseMatrix64F Pdi = wrap(precisions, imo, dimTrait, dimTrait); // TODO Only if needed - final DenseMatrix64F Pdj = wrap(precisions, jmo, dimTrait, dimTrait); // TODO Only if needed + final DMatrixRMaj Pdi = wrap(precisions, imo, dimTrait, dimTrait); // TODO Only if needed + final DMatrixRMaj Pdj = wrap(precisions, jmo, dimTrait, dimTrait); // TODO Only if needed // TODO End fix @@ -237,8 +237,8 @@ // final double lpi = partials[ibo + dimTrait + 2 * dimTrait * dimTrait]; // final double lpj = partials[jbo + dimTrait + 2 * dimTrait * dimTrait]; - final DenseMatrix64F Pi = wrap(partials, ibo + dimTrait, dimTrait, dimTrait); - final DenseMatrix64F Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Pi = wrap(partials, ibo + dimTrait, dimTrait, dimTrait); + final DMatrixRMaj Pj = wrap(partials, jbo + dimTrait, dimTrait, dimTrait); if (TIMING) { endTime("peel1"); @@ -254,8 +254,8 @@ InversionResult ci; InversionResult cj; - final DenseMatrix64F Pip = matrix2; - final DenseMatrix64F Pjp = matrix3; + final DMatrixRMaj Pip = matrix2; + final DMatrixRMaj Pjp = matrix3; // boolean useVariance = anyDiagonalInfinities(Pi) || anyDiagonalInfinities(Pj); final boolean useVariancei = anyDiagonalInfinities(Pi); @@ -263,43 +263,43 @@ if (useVariancei) { - final DenseMatrix64F Vip = matrix0; - final DenseMatrix64F Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); -// CommonOps.add(Vi, vi, Vd, Vip); // TODO Fix - CommonOps.add(Vi, Vdi, Vip); + final DMatrixRMaj Vip = matrix0; + final DMatrixRMaj Vi = wrap(partials, ibo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); +// CommonOps_DDRM.add(Vi, vi, Vd, Vip); // TODO Fix + CommonOps_DDRM.add(Vi, Vdi, Vip); ci = safeInvert(Vip, Pip, true); } else { - final DenseMatrix64F PiPlusPd = matrix0; -// CommonOps.add(Pi, 1.0 / vi, Pd, PiPlusPd); // TODO Fix - CommonOps.add(Pi, Pdi, PiPlusPd); - final DenseMatrix64F PiPlusPdInv = new DenseMatrix64F(dimTrait, dimTrait); + final DMatrixRMaj PiPlusPd = matrix0; +// CommonOps_DDRM.add(Pi, 1.0 / vi, Pd, PiPlusPd); // TODO Fix + CommonOps_DDRM.add(Pi, Pdi, PiPlusPd); + final DMatrixRMaj PiPlusPdInv = new DMatrixRMaj(dimTrait, dimTrait); safeInvert(PiPlusPd, PiPlusPdInv, false); - CommonOps.mult(PiPlusPdInv, Pi, Pip); - CommonOps.mult(Pi, Pip, PiPlusPdInv); - CommonOps.add(Pi, -1, PiPlusPdInv, Pip); + CommonOps_DDRM.mult(PiPlusPdInv, Pi, Pip); + CommonOps_DDRM.mult(Pi, Pip, PiPlusPdInv); + CommonOps_DDRM.add(Pi, -1, PiPlusPdInv, Pip); ci = safeDeterminant(Pip, false); } if (useVariancej) { - final DenseMatrix64F Vjp = matrix1; - final DenseMatrix64F Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); -// CommonOps.add(Vj, vj, Vd, Vjp); // TODO Fix - CommonOps.add(Vj, Vdj, Vjp); + final DMatrixRMaj Vjp = matrix1; + final DMatrixRMaj Vj = wrap(partials, jbo + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); +// CommonOps_DDRM.add(Vj, vj, Vd, Vjp); // TODO Fix + CommonOps_DDRM.add(Vj, Vdj, Vjp); cj = safeInvert(Vjp, Pjp, true); } else { - final DenseMatrix64F PjPlusPd = matrix1; -// CommonOps.add(Pj, 1.0 / vj, Pd, PjPlusPd); // TODO Fix - CommonOps.add(Pj, Pdj, PjPlusPd); - final DenseMatrix64F PjPlusPdInv = new DenseMatrix64F(dimTrait, dimTrait); + final DMatrixRMaj PjPlusPd = matrix1; +// CommonOps_DDRM.add(Pj, 1.0 / vj, Pd, PjPlusPd); // TODO Fix + CommonOps_DDRM.add(Pj, Pdj, PjPlusPd); + final DMatrixRMaj PjPlusPdInv = new DMatrixRMaj(dimTrait, dimTrait); safeInvert(PjPlusPd, PjPlusPdInv, false); - CommonOps.mult(PjPlusPdInv, Pj, Pjp); - CommonOps.mult(Pj, Pjp, PjPlusPdInv); - CommonOps.add(Pj, -1, PjPlusPdInv, Pjp); + CommonOps_DDRM.mult(PjPlusPdInv, Pj, Pjp); + CommonOps_DDRM.mult(Pj, Pjp, PjPlusPdInv); + CommonOps_DDRM.add(Pj, -1, PjPlusPdInv, Pjp); cj = safeDeterminant(Pjp, false); } @@ -313,12 +313,12 @@ // A. Partial precision and variance (for later use) using one matrix inversion // final double lpk = lpip + lpjp; -// final DenseMatrix64F Pk = new DenseMatrix64F(dimTrait, dimTrait); - final DenseMatrix64F Pk = matrix4; - CommonOps.add(Pip, Pjp, Pk); +// final DMatrixRMaj Pk = new DMatrixRMaj(dimTrait, dimTrait); + final DMatrixRMaj Pk = matrix4; + CommonOps_DDRM.add(Pip, Pjp, Pk); -// final DenseMatrix64F Vk = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Vk = matrix5; +// final DMatrixRMaj Vk = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Vk = matrix5; // if (useVariance) { //// InversionResult ck = @@ -475,9 +475,9 @@ } } -// final DenseMatrix64F Vt = new DenseMatrix64F(dimTrait, dimTrait); -// final DenseMatrix64F Vt = matrix6; -// CommonOps.add(Vip, Vjp, Vt); +// final DMatrixRMaj Vt = new DMatrixRMaj(dimTrait, dimTrait); +// final DMatrixRMaj Vt = matrix6; +// CommonOps_DDRM.add(Vip, Vjp, Vt); // if (DEBUG) { // System.err.println("Vt: " + Vt); @@ -487,16 +487,16 @@ - ck.getEffectiveDimension(); // System.err.println(ci.getDeterminant()); -// System.err.println(CommonOps.det(Vip)); +// System.err.println(CommonOps_DDRM.det(Vip)); // // System.err.println(cj.getDeterminant()); -// System.err.println(CommonOps.det(Vjp)); +// System.err.println(CommonOps_DDRM.det(Vjp)); // // System.err.println(1.0 / ck.getDeterminant()); -// System.err.println(CommonOps.det(Vk)); +// System.err.println(CommonOps_DDRM.det(Vk)); remainder += -dimensionChange * LOG_SQRT_2_PI - 0.5 * -// (Math.log(CommonOps.det(Vip)) + Math.log(CommonOps.det(Vjp)) - Math.log(CommonOps.det(Vk))) +// (Math.log(CommonOps_DDRM.det(Vip)) + Math.log(CommonOps_DDRM.det(Vjp)) - Math.log(CommonOps_DDRM.det(Vk))) (Math.log(ci.getDeterminant()) + Math.log(cj.getDeterminant()) + Math.log(ck.getDeterminant())) - 0.5 * (SSi + SSj - SSk); --- a/src/dr/evomodel/treedatalikelihood/continuous/IntegratedFactorAnalysisLikelihood.java +++ b/src/dr/evomodel/treedatalikelihood/continuous/IntegratedFactorAnalysisLikelihood.java @@ -38,7 +38,7 @@ import dr.math.matrixAlgebra.WrappedVector; import dr.math.matrixAlgebra.missingData.InversionResult; import dr.xml.*; -import org.ejml.data.DenseMatrix64F; +import org.ejml.data.DMatrixRMaj; import java.util.ArrayList; import java.util.Arrays; @@ -229,7 +229,7 @@ private final double nuggetPrecision; - private void computePrecisionForTaxon(final DenseMatrix64F precision, final int taxon, + private void computePrecisionForTaxon(final DMatrixRMaj precision, final int taxon, final int numFactors) { final double[] observed = observedIndicators[taxon]; @@ -249,7 +249,7 @@ } } - private InversionResult fillInMeanForTaxon(final WrappedVector output, final DenseMatrix64F precision, final int taxon) { + private InversionResult fillInMeanForTaxon(final WrappedVector output, final DMatrixRMaj precision, final int taxon) { final double[] observed = observedIndicators[taxon]; final Parameter Y = traitParameter.getParameter(taxon); @@ -269,8 +269,8 @@ tmp[row] = sum; } - DenseMatrix64F B = DenseMatrix64F.wrap(numFactors, 1, tmp); - DenseMatrix64F X = DenseMatrix64F.wrap(numFactors, 1, tmp2); + DMatrixRMaj B = DMatrixRMaj.wrap(numFactors, 1, tmp); + DMatrixRMaj X = DMatrixRMaj.wrap(numFactors, 1, tmp2); InversionResult ci = safeSolve(precision, B, X, true); @@ -294,7 +294,7 @@ return sum; } - private double computeFactorInnerProduct(final WrappedVector mean, final DenseMatrix64F precision) { + private double computeFactorInnerProduct(final WrappedVector mean, final DMatrixRMaj precision) { // Compute \mu_i^t P_i \mu^t double sum = 0; for (int row = 0; row < numFactors; ++row) { @@ -334,7 +334,7 @@ return logDet; } - private void makeCompletedUnobserved(final DenseMatrix64F matrix, double diagonal) { + private void makeCompletedUnobserved(final DMatrixRMaj matrix, double diagonal) { for (int row = 0; row < numFactors; ++row) { for (int col = 0; col < numFactors; ++col) { double x = (row == col) ? diagonal : 0.0; @@ -345,8 +345,8 @@ private void computePartialsAndRemainders() { - final DenseMatrix64F precision = new DenseMatrix64F(numFactors, numFactors); - final DenseMatrix64F variance = new DenseMatrix64F(numFactors, numFactors); + final DMatrixRMaj precision = new DMatrixRMaj(numFactors, numFactors); + final DMatrixRMaj variance = new DMatrixRMaj(numFactors, numFactors); int partialsOffset = 0; for (int taxon = 0; taxon < numTaxa; ++taxon) { --- a/src/dr/evomodel/treedatalikelihood/continuous/RepeatedMeasuresTraitLikelihood.java +++ b/src/dr/evomodel/treedatalikelihood/continuous/RepeatedMeasuresTraitLikelihood.java @@ -39,7 +39,7 @@ import dr.math.matrixAlgebra.WrappedVector; import dr.math.matrixAlgebra.missingData.InversionResult; import dr.xml.*; -import org.ejml.data.DenseMatrix64F; +import org.ejml.data.DMatrixRMaj; import java.util.ArrayList; import java.util.Arrays; @@ -230,7 +230,7 @@ private final double nuggetPrecision; -// private void computePrecisionForTaxon(final DenseMatrix64F precision, final int taxon, +// private void computePrecisionForTaxon(final DMatrixRMaj precision, final int taxon, // final int numFactors) { // // final double[] observed = observedIndicators[taxon]; @@ -250,7 +250,7 @@ // } // } -// private InversionResult fillInMeanForTaxon(final WrappedVector output, final DenseMatrix64F precision, final int taxon) { +// private InversionResult fillInMeanForTaxon(final WrappedVector output, final DMatrixRMaj precision, final int taxon) { // // final double[] observed = observedIndicators[taxon]; // final Parameter Y = traitParameter.getParameter(taxon); @@ -270,8 +270,8 @@ // tmp[row] = sum; // } // -// DenseMatrix64F B = DenseMatrix64F.wrap(numFactors, 1, tmp); -// DenseMatrix64F X = DenseMatrix64F.wrap(numFactors, 1, tmp2); +// DMatrixRMaj B = DMatrixRMaj.wrap(numFactors, 1, tmp); +// DMatrixRMaj X = DMatrixRMaj.wrap(numFactors, 1, tmp2); // // InversionResult ci = safeSolve(precision, B, X, true); // @@ -295,7 +295,7 @@ // return sum; // } -// private double computeFactorInnerProduct(final WrappedVector mean, final DenseMatrix64F precision) { +// private double computeFactorInnerProduct(final WrappedVector mean, final DMatrixRMaj precision) { // // Compute \mu_i^t P_i \mu^t // double sum = 0; // for (int row = 0; row < numFactors; ++row) { @@ -320,7 +320,7 @@ // return det; // } -// private void makeCompletedUnobserved(final DenseMatrix64F matrix, double diagonal) { +// private void makeCompletedUnobserved(final DMatrixRMaj matrix, double diagonal) { // for (int row = 0; row < numFactors; ++row) { // for (int col = 0; col < numFactors; ++col) { // double x = (row == col) ? diagonal : 0.0; @@ -331,8 +331,8 @@ private void computePartialsAndRemainders() { -// final DenseMatrix64F precision = new DenseMatrix64F(numFactors, numFactors); -// final DenseMatrix64F variance = new DenseMatrix64F(numFactors, numFactors); +// final DMatrixRMaj precision = new DMatrixRMaj(numFactors, numFactors); +// final DMatrixRMaj variance = new DMatrixRMaj(numFactors, numFactors); // // int partialsOffset = 0; // for (int taxon = 0; taxon < numTaxa; ++taxon) { --- a/src/dr/evomodel/treedatalikelihood/preorder/AbstractValuesViaFullConditionalDelegate.java +++ b/src/dr/evomodel/treedatalikelihood/preorder/AbstractValuesViaFullConditionalDelegate.java @@ -5,7 +5,7 @@ import dr.evomodel.continuous.MultivariateDiffusionModel; import dr.evomodel.treedatalikelihood.continuous.*; import dr.math.matrixAlgebra.WrappedVector; -import org.ejml.data.DenseMatrix64F; +import org.ejml.data.DMatrixRMaj; import static dr.math.matrixAlgebra.missingData.MissingOps.wrap; @@ -63,7 +63,7 @@ System.err.println("Missing tip = " + node.getNumber() + " (" + nodeBuffer + "), trait = " + trait); final WrappedVector preMean = new WrappedVector.Raw(conditionalNodeBuffer, partialOffset, dimTrait); - final DenseMatrix64F preVar = wrap(conditionalNodeBuffer, partialOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); + final DMatrixRMaj preVar = wrap(conditionalNodeBuffer, partialOffset + dimTrait + dimTrait * dimTrait, dimTrait, dimTrait); final WrappedVector postObs = new WrappedVector.Raw(partialNodeBuffer, partialOffset, dimTrait); --- a/src/dr/evomodel/treedatalikelihood/preorder/ConditionalVarianceAndTransform2.java +++ b/src/dr/evomodel/treedatalikelihood/preorder/ConditionalVarianceAndTransform2.java @@ -1,8 +1,8 @@ package dr.evomodel.treedatalikelihood.preorder; import dr.math.matrixAlgebra.WrappedVector; -import org.ejml.data.DenseMatrix64F; -import org.ejml.ops.CommonOps; +import org.ejml.data.DMatrixRMaj; +import org.ejml.dense.row.CommonOps_DDRM; import static dr.math.matrixAlgebra.missingData.MissingOps.gatherRowsAndColumns; @@ -23,8 +23,8 @@ * \bar{\Sigma} = \Sigma_{11} - \Sigma_{12}\Sigma_{22}^1\Sigma{21} */ - final private DenseMatrix64F sBar; - final private DenseMatrix64F affineTransform; + final private DMatrixRMaj sBar; + final private DMatrixRMaj affineTransform; private final int[] missingIndices; private final int[] notMissingIndices; @@ -36,9 +36,9 @@ private static final boolean DEBUG = false; private double[][] cholesky = null; - private DenseMatrix64F sBarInv = null; + private DMatrixRMaj sBarInv = null; - ConditionalVarianceAndTransform2(final DenseMatrix64F variance, + ConditionalVarianceAndTransform2(final DMatrixRMaj variance, final int[] missingIndices, final int[] notMissingIndices) { assert (missingIndices.length + notMissingIndices.length == variance.getNumRows()); @@ -51,44 +51,44 @@ System.err.println("variance:\n" + variance); } - DenseMatrix64F S22 = new DenseMatrix64F(notMissingIndices.length, notMissingIndices.length); + DMatrixRMaj S22 = new DMatrixRMaj(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); + DMatrixRMaj S22Inv = new DMatrixRMaj(notMissingIndices.length, notMissingIndices.length); + CommonOps_DDRM.invert(S22, S22Inv); if (DEBUG) { System.err.println("S22Inv:\n" + S22Inv); } - DenseMatrix64F S12 = new DenseMatrix64F(missingIndices.length, notMissingIndices.length); + DMatrixRMaj S12 = new DMatrixRMaj(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); + DMatrixRMaj S12S22Inv = new DMatrixRMaj(missingIndices.length, notMissingIndices.length); + CommonOps_DDRM.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); + DMatrixRMaj S12S22InvS21 = new DMatrixRMaj(missingIndices.length, missingIndices.length); + CommonOps_DDRM.multTransB(S12S22Inv, S12, S12S22InvS21); if (DEBUG) { System.err.println("S12S22InvS21:\n" + S12S22InvS21); } - sBar = new DenseMatrix64F(missingIndices.length, missingIndices.length); + sBar = new DMatrixRMaj(missingIndices.length, missingIndices.length); gatherRowsAndColumns(variance, sBar, missingIndices, missingIndices); - CommonOps.subtract(sBar, S12S22InvS21, sBar); + CommonOps_DDRM.subtract(sBar, S12S22InvS21, sBar); if (DEBUG) { @@ -141,18 +141,18 @@ return cholesky; } -// public final DenseMatrix64F getAffineTransform() { +// public final DMatrixRMaj getAffineTransform() { // return affineTransform; // } - final DenseMatrix64F getConditionalVariance() { + final DMatrixRMaj getConditionalVariance() { return sBar; } - final DenseMatrix64F getConditionalPrecision() { + final DMatrixRMaj getConditionalPrecision() { if (sBarInv == null) { - sBarInv = new DenseMatrix64F(numMissing, numMissing); - CommonOps.invert(sBar, sBarInv); + sBarInv = new DMatrixRMaj(numMissing, numMissing); + CommonOps_DDRM.invert(sBar, sBarInv); } return sBarInv; } --- a/src/dr/evomodel/treedatalikelihood/preorder/MultivariateConditionalOnTipsRealizedDelegate.java +++ b/src/dr/evomodel/treedatalikelihood/preorder/MultivariateConditionalOnTipsRealizedDelegate.java @@ -6,8 +6,8 @@ import dr.math.distributions.MultivariateNormalDistribution; import dr.math.matrixAlgebra.Matrix; import dr.math.matrixAlgebra.WrappedVector; -import org.ejml.data.DenseMatrix64F; -import org.ejml.ops.CommonOps; +import org.ejml.data.DMatrixRMaj; +import org.ejml.dense.row.CommonOps_DDRM; import static dr.math.matrixAlgebra.missingData.MissingOps.*; @@ -41,14 +41,14 @@ // 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); - CommonOps.mult(Pd, wrap(partialPriorBuffer, offsetPartial + dimTrait, dimTrait, dimTrait), priorPrec); + final DMatrixRMaj rootPrec = wrap(partialNodeBuffer, offsetPartial + dimTrait, dimTrait, dimTrait); + final DMatrixRMaj priorPrec = new DMatrixRMaj(dimTrait, dimTrait); + CommonOps_DDRM.mult(Pd, wrap(partialPriorBuffer, offsetPartial + dimTrait, dimTrait, dimTrait), priorPrec); - final DenseMatrix64F totalPrec = new DenseMatrix64F(dimTrait, dimTrait); - CommonOps.add(rootPrec, priorPrec, totalPrec); + final DMatrixRMaj totalPrec = new DMatrixRMaj(dimTrait, dimTrait); + CommonOps_DDRM.add(rootPrec, priorPrec, totalPrec); - final DenseMatrix64F totalVar = new DenseMatrix64F(dimTrait, dimTrait); + final DMatrixRMaj totalVar = new DMatrixRMaj(dimTrait, dimTrait); safeInvert(totalPrec, totalVar, false); final double[] tmp = new double[dimTrait]; @@ -93,7 +93,7 @@ } } -// boolean extremeValue(final DenseMatrix64F x) { +// boolean extremeValue(final DMatrixRMaj x) { // return extremeValue(new WrappedVector.Raw(x.getData(), 0, x.getNumElements())); // } // @@ -130,7 +130,7 @@ final int offsetPartial, final double branchPrecision) { - final DenseMatrix64F P0 = wrap(partialNodeBuffer, offsetPartial + dimTrait, dimTrait, dimTrait); + final DMatrixRMaj P0 = wrap(partialNodeBuffer, offsetPartial + dimTrait, dimTrait, dimTrait); final int missingCount = countFiniteDiagonals(P0); if (missingCount == 0) { // Completely observed @@ -167,26 +167,26 @@ final int[] observed = indices.getComplement(); final int[] missing = indices.getArray(); - final DenseMatrix64F V1 = new DenseMatrix64F(dimTrait, dimTrait); - CommonOps.scale(1.0 / branchPrecision, Vd, V1); + final DMatrixRMaj V1 = new DMatrixRMaj(dimTrait, dimTrait); + CommonOps_DDRM.scale(1.0 / branchPrecision, Vd, V1); ConditionalVarianceAndTransform2 transform = new ConditionalVarianceAndTransform2( V1, missing, observed ); // TODO Cache (via delegated function) - final DenseMatrix64F cP0 = new DenseMatrix64F(missing.length, missing.length); + final DMatrixRMaj cP0 = new DMatrixRMaj(missing.length, missing.length); gatherRowsAndColumns(P0, cP0, missing, missing); final WrappedVector cM2 = transform.getConditionalMean( partialNodeBuffer, offsetPartial, // Tip value sample, offsetParent); // Parent value - final DenseMatrix64F cP1 = transform.getConditionalPrecision(); + final DMatrixRMaj cP1 = transform.getConditionalPrecision(); - final DenseMatrix64F cP2 = new DenseMatrix64F(missing.length, missing.length); - final DenseMatrix64F cV2 = new DenseMatrix64F(missing.length, missing.length); - CommonOps.add(cP0, cP1, cP2); + final DMatrixRMaj cP2 = new DMatrixRMaj(missing.length, missing.length); + final DMatrixRMaj cV2 = new DMatrixRMaj(missing.length, missing.length); + CommonOps_DDRM.add(cP0, cP1, cP2); safeInvert(cP2, cV2, false); double[][] cC2 = getCholeskyOfVariance(cV2.getData(), missing.length); @@ -205,8 +205,8 @@ final WrappedVector M0 = new WrappedVector.Raw(partialNodeBuffer, offsetPartial, dimTrait); final WrappedVector M1 = new WrappedVector.Raw(sample, offsetParent, dimTrait); - final DenseMatrix64F P1 = new DenseMatrix64F(dimTrait, dimTrait); - CommonOps.scale(branchPrecision, Pd, P1); + final DMatrixRMaj P1 = new DMatrixRMaj(dimTrait, dimTrait); + CommonOps_DDRM.scale(branchPrecision, Pd, P1); final WrappedVector newSample = new WrappedVector.Raw(sample, offsetSample, dimTrait); @@ -246,25 +246,25 @@ if (!Double.isInfinite(branchPrecision)) { final WrappedVector M0 = new WrappedVector.Raw(partialNodeBuffer, offsetPartial, dimTrait); - final DenseMatrix64F P0 = wrap(partialNodeBuffer, offsetPartial + dimTrait, dimTrait, dimTrait); + final DMatrixRMaj P0 = wrap(partialNodeBuffer, offsetPartial + dimTrait, dimTrait, dimTrait); final WrappedVector M1; - final DenseMatrix64F P1; + final DMatrixRMaj P1; if (hasNoDrift) { M1 = new WrappedVector.Raw(sample, offsetParent, dimTrait); - P1 = new DenseMatrix64F(dimTrait, dimTrait); - CommonOps.scale(branchPrecision, Pd, P1); + P1 = new DMatrixRMaj(dimTrait, dimTrait); + CommonOps_DDRM.scale(branchPrecision, Pd, P1); } else { M1 = getMeanWithDrift(sample, offsetParent, displacementBuffer, dimTrait); - P1 = DenseMatrix64F.wrap(dimTrait, dimTrait, precisionBuffer); + P1 = DMatrixRMaj.wrap(dimTrait, dimTrait, precisionBuffer); } // boolean DEBUG_PRECISION = false; // // if (DEBUG_PRECISION) { -// DenseMatrix64F tP1 = new DenseMatrix64F(dimTrait, dimTrait); -// CommonOps.scale(branchPrecision, Pd, tP1); +// DMatrixRMaj tP1 = new DMatrixRMaj(dimTrait, dimTrait); +// CommonOps_DDRM.scale(branchPrecision, Pd, tP1); // // for (int i = 0; i < dimTrait; ++i) { // for (int j = 0; j < dimTrait; ++j) { @@ -277,10 +277,10 @@ // } final WrappedVector M2 = new WrappedVector.Raw(tmpMean, 0, dimTrait); - final DenseMatrix64F P2 = new DenseMatrix64F(dimTrait, dimTrait); - final DenseMatrix64F V2 = new DenseMatrix64F(dimTrait, dimTrait); + final DMatrixRMaj P2 = new DMatrixRMaj(dimTrait, dimTrait); + final DMatrixRMaj V2 = new DMatrixRMaj(dimTrait, dimTrait); - CommonOps.add(P0, P1, P2); + CommonOps_DDRM.add(P0, P1, P2); safeInvert(P2, V2, false); weightedAverage(M0, P0, M1, P1, M2, V2, dimTrait); --- a/src/dr/evomodel/treedatalikelihood/preorder/ProcessSimulationDelegate.java +++ b/src/dr/evomodel/treedatalikelihood/preorder/ProcessSimulationDelegate.java @@ -35,7 +35,7 @@ import dr.inference.model.Model; import dr.inference.model.ModelListener; import dr.math.matrixAlgebra.*; -import org.ejml.data.DenseMatrix64F; +import org.ejml.data.DMatrixRMaj; import java.util.List; import java.util.Map; @@ -161,8 +161,8 @@ final ContinuousDataLikelihoodDelegate likelihoodDelegate; double[] diffusionVariance; - DenseMatrix64F Vd; - DenseMatrix64F Pd; + DMatrixRMaj Vd; + DMatrixRMaj Pd; double[][] cholesky; Map= 2? UnrolledDeterminantFromMinor.det(mat):mat.get(0); + return numCol >= 2? UnrolledDeterminantFromMinor_DDRM.det(mat):mat.get(0); } else { - LUDecompositionAlt_D64 alg = new LUDecompositionAlt_D64(); + LUDecompositionAlt_DDRM alg = new LUDecompositionAlt_DDRM(); if(alg.inputModified()) { mat = mat.copy(); } @@ -190,7 +190,7 @@ } } - public static double invertAndGetDeterminant(DenseMatrix64F mat, DenseMatrix64F result) { + public static double invertAndGetDeterminant(DMatrixRMaj mat, DMatrixRMaj result) { final int numCol = mat.getNumCols(); final int numRow = mat.getNumRows(); @@ -201,18 +201,18 @@ if (numCol <= 5) { if (numCol >= 2) { - UnrolledInverseFromMinor.inv(mat, result); + UnrolledInverseFromMinor_DDRM.inv(mat, result); } else { result.set(0, 1.0D / mat.get(0)); } return numCol >= 2 ? - UnrolledDeterminantFromMinor.det(mat) : + UnrolledDeterminantFromMinor_DDRM.det(mat) : mat.get(0); } else { - LUDecompositionAlt_D64 alg = new LUDecompositionAlt_D64(); - LinearSolverLu_D64 solver = new LinearSolverLu_D64(alg); + LUDecompositionAlt_DDRM alg = new LUDecompositionAlt_DDRM(); + LinearSolverLu_DDRM solver = new LinearSolverLu_DDRM(alg); if (solver.modifiesA()) { mat = mat.copy(); } @@ -228,7 +228,7 @@ } } - public static InversionResult safeDeterminant(DenseMatrix64F source, boolean invert) { + public static InversionResult safeDeterminant(DMatrixRMaj source, boolean invert) { final int finiteCount = countFiniteNonZeroDiagonals(source); InversionResult result; @@ -236,10 +236,10 @@ if (finiteCount == 0) { result = new InversionResult(NOT_OBSERVED, 0, 0); } else { - LinearSolver solver = LinearSolverFactory.pseudoInverse(true); + LinearSolver solver = LinearSolverFactory_DDRM.pseudoInverse(true); solver.setA(source); - SingularValueDecomposition svd = solver.getDecomposition(); + SingularValueDecomposition_F64 svd = solver.getDecomposition(); double[] values = svd.getSingularValues(); @@ -267,7 +267,7 @@ return result; } - public static InversionResult safeSolve(DenseMatrix64F A, + public static InversionResult safeSolve(DMatrixRMaj A, WrappedVector b, WrappedVector x, boolean getDeterminat) { @@ -275,8 +275,8 @@ assert(A.getNumRows() == dim && A.getNumCols() == dim); - final DenseMatrix64F B = wrap(b.getBuffer(), b.getOffset(), dim, 1); - final DenseMatrix64F X = new DenseMatrix64F(dim, 1); + final DMatrixRMaj B = wrap(b.getBuffer(), b.getOffset(), dim, 1); + final DMatrixRMaj X = new DMatrixRMaj(dim, 1); InversionResult ir = safeSolve(A, B, X, getDeterminat); @@ -288,7 +288,7 @@ return ir; } - public static InversionResult safeSolve(DenseMatrix64F A, DenseMatrix64F B, DenseMatrix64F X, boolean getDeterminant) { + public static InversionResult safeSolve(DMatrixRMaj A, DMatrixRMaj B, DMatrixRMaj X, boolean getDeterminant) { final int finiteCount = countFiniteNonZeroDiagonals(A); @@ -298,7 +298,7 @@ result = new InversionResult(NOT_OBSERVED, 0, 0); } else { - LinearSolver solver = LinearSolverFactory.pseudoInverse(true); + LinearSolver solver = LinearSolverFactory_DDRM.pseudoInverse(true); solver.setA(A); solver.solve(B, X); @@ -306,7 +306,7 @@ double det = 1; if (getDeterminant) { - SingularValueDecomposition svd = solver.getDecomposition(); + SingularValueDecomposition_F64 svd = solver.getDecomposition(); double[] values = svd.getSingularValues(); for (int i = 0; i < values.length; ++i) { @@ -325,7 +325,7 @@ } - public static InversionResult safeInvert(DenseMatrix64F source, DenseMatrix64F destination, boolean getDeterminant) { + public static InversionResult safeInvert(DMatrixRMaj source, DMatrixRMaj destination, boolean getDeterminant) { final int dim = source.getNumCols(); final int finiteCount = countFiniteNonZeroDiagonals(source); @@ -335,7 +335,7 @@ if (getDeterminant) { det = invertAndGetDeterminant(source, destination); } else { - CommonOps.invert(source, destination); + CommonOps_DDRM.invert(source, destination); } return new InversionResult(FULLY_OBSERVED, dim, det); } else { @@ -346,14 +346,14 @@ final int[] finiteIndices = new int[finiteCount]; getFiniteNonZeroDiagonalIndices(source, finiteIndices); - final DenseMatrix64F subSource = new DenseMatrix64F(finiteCount, finiteCount); + final DMatrixRMaj subSource = new DMatrixRMaj(finiteCount, finiteCount); gatherRowsAndColumns(source, subSource, finiteIndices, finiteIndices); - final DenseMatrix64F inverseSubSource = new DenseMatrix64F(finiteCount, finiteCount); + final DMatrixRMaj inverseSubSource = new DMatrixRMaj(finiteCount, finiteCount); if (getDeterminant) { det = invertAndGetDeterminant(subSource, inverseSubSource); } else { - CommonOps.invert(subSource, inverseSubSource); + CommonOps_DDRM.invert(subSource, inverseSubSource); } scatterRowsAndColumns(inverseSubSource, destination, finiteIndices, finiteIndices, true); @@ -364,11 +364,11 @@ } public static void weightedAverage(final WrappedVector mi, - final DenseMatrix64F Pi, + final DMatrixRMaj Pi, final WrappedVector mj, - final DenseMatrix64F Pj, + final DMatrixRMaj Pj, final WrappedVector mk, - final DenseMatrix64F Vk, + final DMatrixRMaj Vk, final int dimTrait) { final double[] tmp = new double[dimTrait]; for (int g = 0; g < dimTrait; ++g) { @@ -390,13 +390,13 @@ public static void weightedAverage(final double[] ipartial, final int ibo, - final DenseMatrix64F Pi, + final DMatrixRMaj Pi, final double[] jpartial, final int jbo, - final DenseMatrix64F Pj, + final DMatrixRMaj Pj, final double[] kpartial, final int kbo, - final DenseMatrix64F Vk, + final DMatrixRMaj Vk, final int dimTrait) { final double[] tmp = new double[dimTrait]; for (int g = 0; g < dimTrait; ++g) {