[med-svn] [Git][med-team/beast-mcmc][master] 7 commits: Preparing to use newer version of libejml-java. Performing the source...
Andrius Merkys
gitlab at salsa.debian.org
Tue May 21 10:29:20 BST 2019
Andrius Merkys pushed to branch master at Debian Med / beast-mcmc
Commits:
0f836acd by Andrius Merkys at 2019-05-21T08:30:15Z
Preparing to use newer version of libejml-java. Performing the source transition using convert_to_ejml31.py script from libejml-java.
- - - - -
6d7d3642 by Andrius Merkys at 2019-05-21T08:41:11Z
Migrating non-automigrated org.ejml.alg.dense.misc.UnrolledInverseFromMinor.
- - - - -
d65b3ffc by Andrius Merkys at 2019-05-21T09:03:00Z
Migrating non-automigrated org.ejml.alg.dense.decomposition.lu.LUDecompositionAlt_D64.
- - - - -
26c49758 by Andrius Merkys at 2019-05-21T09:18:29Z
Finishing migration.
- - - - -
a15aca9e by Andrius Merkys at 2019-05-21T09:25:27Z
Adding versioned b-dep on libejml-java.
- - - - -
8b6ddc1a by Andrius Merkys at 2019-05-21T09:27:11Z
Adding patch description for update-ejml.patch.
- - - - -
d5cab6e9 by Andrius Merkys at 2019-05-21T09:28:52Z
Preparing for the next upload.
- - - - -
5 changed files:
- debian/changelog
- debian/control
- debian/patches/fix_classpath_in_build_xml.patch
- debian/patches/series
- + debian/patches/update-ejml.patch
Changes:
=====================================
debian/changelog
=====================================
@@ -1,3 +1,10 @@
+beast-mcmc (1.10.4+dfsg-2) UNRELEASED; urgency=medium
+
+ * Team upload.
+ * Migrating to libejml-java 0.38.
+
+ -- Andrius Merkys <merkys at debian.org> Tue, 21 May 2019 05:27:43 -0400
+
beast-mcmc (1.10.4+dfsg-1) unstable; urgency=medium
[ Andreas Tille ]
=====================================
debian/control
=====================================
@@ -22,7 +22,7 @@ Build-Depends: debhelper (>= 11~),
junit4,
libmtj-java,
libitext1-java,
- libejml-java,
+ libejml-java (>= 0.38),
libjlapack-java
Standards-Version: 4.2.1
Vcs-Browser: https://salsa.debian.org/med-team/beast-mcmc
=====================================
debian/patches/fix_classpath_in_build_xml.patch
=====================================
@@ -13,7 +13,7 @@ Forwarded: no
<property name="dist" location="${build}/dist"/>
<property name="main_class_BEAST" value="dr.app.beast.BeastMain"/>
-@@ -48,6 +49,21 @@
+@@ -48,6 +49,31 @@
<path id="classpath">
<fileset dir="${lib}" includes="**/*.jar"/>
@@ -31,11 +31,21 @@ Forwarded: no
+ <fileset dir="${deblib}" includes="mtj.jar"/>
+ <fileset dir="${deblib}" includes="options.jar"/>
+ <fileset dir="${deblib}" includes="EJML.jar"/>
++ <fileset dir="${deblib}" includes="ejml-all.jar"/>
++ <fileset dir="${deblib}" includes="ejml-all.jar"/>
++ <fileset dir="${deblib}" includes="ejml-cdense.jar"/>
++ <fileset dir="${deblib}" includes="ejml-core.jar"/>
++ <fileset dir="${deblib}" includes="ejml-ddense.jar"/>
++ <fileset dir="${deblib}" includes="ejml-dsparse.jar"/>
++ <fileset dir="${deblib}" includes="ejml-experimental.jar"/>
++ <fileset dir="${deblib}" includes="ejml-fdense.jar"/>
++ <fileset dir="${deblib}" includes="ejml-simple.jar"/>
++ <fileset dir="${deblib}" includes="ejml-zdense.jar"/>
+ <fileset dir="${deblib}" includes="jlapack-lapack.jar"/>
</path>
<!-- start -->
-@@ -159,17 +175,6 @@
+@@ -161,17 +187,6 @@
<include name="org/virion/jam/**/*.png"/>
<include name="dr/**/*.properties"/>
</fileset>
@@ -53,7 +63,7 @@ Forwarded: no
</jar>
<!-- Put everything in ${build} into the beauti.jar file -->
-@@ -226,11 +231,6 @@
+@@ -228,11 +243,6 @@
<fileset dir="${src}">
<include name="dr/**/*.png"/>
</fileset>
@@ -67,7 +77,7 @@ Forwarded: no
<!-- Put everything in ${build} into the trace.jar file -->
--- a/.classpath
+++ b/.classpath
-@@ -1,20 +1,18 @@
+@@ -1,20 +1,27 @@
<?xml version="1.0" encoding="UTF-8"?>
<classpath>
<classpathentry kind="src" path="src"/>
@@ -98,6 +108,15 @@ Forwarded: no
+ <classpathentry kind="lib" path="/usr/share/java/mtj.jar"/>
+ <classpathentry kind="lib" path="/usr/share/java/commons-math.jar"/>
+ <classpathentry kind="lib" path="/usr/share/java/EJML.jar"/>
++ <classpathentry kind="lib" path="/usr/share/java/ejml-all.jar"/>
++ <classpathentry kind="lib" path="/usr/share/java/ejml-cdense.jar"/>
++ <classpathentry kind="lib" path="/usr/share/java/ejml-core.jar"/>
++ <classpathentry kind="lib" path="/usr/share/java/ejml-ddense.jar"/>
++ <classpathentry kind="lib" path="/usr/share/java/ejml-dsparse.jar"/>
++ <classpathentry kind="lib" path="/usr/share/java/ejml-experimental.jar"/>
++ <classpathentry kind="lib" path="/usr/share/java/ejml-fdense.jar"/>
++ <classpathentry kind="lib" path="/usr/share/java/ejml-simple.jar"/>
++ <classpathentry kind="lib" path="/usr/share/java/ejml-zdense.jar"/>
+ <classpathentry kind="lib" path="/usr/share/java/jlapack-lapack.jar"/>
<classpathentry kind="output" path="bin"/>
</classpath>
=====================================
debian/patches/series
=====================================
@@ -3,3 +3,4 @@ fix_classpath_in_build_xml.patch
create_soname.patch
fix_encoding.patch
ignore_mac_install.patch
+update-ejml.patch
=====================================
debian/patches/update-ejml.patch
=====================================
@@ -0,0 +1,1755 @@
+Author: Andrius Merkys <merkys at debian.org>
+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<String, Long> 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<String, Long> 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<PartiallyMissingInformation.HashedIntArray,
+@@ -221,7 +221,7 @@
+ double[][] diffusionPrecision = diffusionModel.getPrecisionmatrix();
+ diffusionVariance = getVectorizedVarianceFromPrecision(diffusionPrecision);
+ Vd = wrap(diffusionVariance, 0, dimTrait, dimTrait);
+- Pd = new DenseMatrix64F(diffusionPrecision);
++ Pd = new DMatrixRMaj(diffusionPrecision);
+ }
+ if (cholesky == null) {
+ cholesky = getCholeskyOfVariance(diffusionVariance, dimTrait);
+--- a/src/dr/inferencexml/operators/EllipticalSliceOperatorParser.java
++++ b/src/dr/inferencexml/operators/EllipticalSliceOperatorParser.java
+@@ -35,7 +35,7 @@
+ import dr.math.distributions.GaussianProcessRandomGenerator;
+ import dr.math.distributions.MultivariateNormalDistribution;
+ import dr.xml.*;
+-import org.ejml.ops.CommonOps;
++import org.ejml.dense.row.CommonOps_DDRM;
+
+ /**
+ */
+--- a/src/dr/math/matrixAlgebra/missingData/MissingOps.java
++++ b/src/dr/math/matrixAlgebra/missingData/MissingOps.java
+@@ -3,17 +3,17 @@
+ import dr.math.matrixAlgebra.Vector;
+ import dr.math.matrixAlgebra.WrappedVector;
+ import dr.util.Transform;
+-import org.ejml.alg.dense.decomposition.lu.LUDecompositionAlt_D64;
+-import org.ejml.alg.dense.decomposition.qr.QRColPivDecompositionHouseholderColumn_D64;
+-import org.ejml.alg.dense.linsol.lu.LinearSolverLu_D64;
+-import org.ejml.alg.dense.misc.UnrolledDeterminantFromMinor;
+-import org.ejml.alg.dense.misc.UnrolledInverseFromMinor;
+-import org.ejml.data.DenseMatrix64F;
+-import org.ejml.factory.LinearSolverFactory;
+-import org.ejml.interfaces.decomposition.QRPDecomposition;
+-import org.ejml.interfaces.decomposition.SingularValueDecomposition;
++import org.ejml.dense.row.decomposition.lu.LUDecompositionAlt_DDRM;
++import org.ejml.dense.row.decomposition.qr.QRColPivDecompositionHouseholderColumn_DDRM;
++import org.ejml.dense.row.linsol.lu.LinearSolverLu_DDRM;
++import org.ejml.dense.row.misc.UnrolledDeterminantFromMinor_DDRM;
++import org.ejml.dense.row.misc.UnrolledInverseFromMinor_DDRM;
++import org.ejml.data.DMatrixRMaj;
++import org.ejml.dense.row.factory.LinearSolverFactory_DDRM;
++import org.ejml.interfaces.decomposition.QRPDecomposition_F64;
++import org.ejml.interfaces.decomposition.SingularValueDecomposition_F64;
+ import org.ejml.interfaces.linsol.LinearSolver;
+-import org.ejml.ops.CommonOps;
++import org.ejml.dense.row.CommonOps_DDRM;
+
+ import java.util.Arrays;
+
+@@ -26,21 +26,21 @@
+ */
+ public class MissingOps {
+
+- public static DenseMatrix64F wrap(final double[] source, final int offset,
++ public static DMatrixRMaj wrap(final double[] source, final int offset,
+ final int numRows, final int numCols) {
+ double[] buffer = new double[numRows * numCols];
+ return wrap(source, offset, numRows, numCols, buffer);
+ }
+
+- public static DenseMatrix64F wrap(final double[] source, final int offset,
++ public static DMatrixRMaj wrap(final double[] source, final int offset,
+ final int numRows, final int numCols,
+ final double[] buffer) {
+ System.arraycopy(source, offset, buffer, 0, numRows * numCols);
+- return DenseMatrix64F.wrap(numRows, numCols, buffer);
++ return DMatrixRMaj.wrap(numRows, numCols, buffer);
+ }
+
+
+- public static void gatherRowsAndColumns(final DenseMatrix64F source, final DenseMatrix64F destination,
++ public static void gatherRowsAndColumns(final DMatrixRMaj source, final DMatrixRMaj destination,
+ final int[] rowIndices, final int[] colIndices) {
+
+ final int rowLength = rowIndices.length;
+@@ -57,7 +57,7 @@
+ }
+ }
+
+- public static void scatterRowsAndColumns(final DenseMatrix64F source, final DenseMatrix64F destination,
++ public static void scatterRowsAndColumns(final DMatrixRMaj source, final DMatrixRMaj destination,
+ final int[] rowIdices, final int[] colIndices, final boolean clear) {
+ if (clear) {
+ Arrays.fill(destination.getData(), 0.0);
+@@ -77,11 +77,11 @@
+ }
+ }
+
+- public static void unwrap(final DenseMatrix64F source, final double[] destination, final int offset) {
++ public static void unwrap(final DMatrixRMaj source, final double[] destination, final int offset) {
+ System.arraycopy(source.getData(), 0, destination, offset, source.getNumElements());
+ }
+
+- public static boolean anyDiagonalInfinities(DenseMatrix64F source) {
++ public static boolean anyDiagonalInfinities(DMatrixRMaj source) {
+ boolean anyInfinities = false;
+ for (int i = 0; i < source.getNumCols() && !anyInfinities; ++i) {
+ if (Double.isInfinite(source.unsafe_get(i, i))) {
+@@ -91,7 +91,7 @@
+ return anyInfinities;
+ }
+
+- public static boolean allFiniteDiagonals(DenseMatrix64F source) {
++ public static boolean allFiniteDiagonals(DMatrixRMaj source) {
+ boolean allFinite = true;
+
+ final int length = source.getNumCols();
+@@ -101,7 +101,7 @@
+ return allFinite;
+ }
+
+- public static int countFiniteDiagonals(DenseMatrix64F source) {
++ public static int countFiniteDiagonals(DMatrixRMaj source) {
+ final int length = source.getNumCols();
+
+ int count = 0;
+@@ -114,7 +114,7 @@
+ return count;
+ }
+
+- public static int countZeroDiagonals(DenseMatrix64F source) {
++ public static int countZeroDiagonals(DMatrixRMaj source) {
+ final int length = source.getNumCols();
+
+ int count = 0;
+@@ -127,7 +127,7 @@
+ return count;
+ }
+
+- public static void getFiniteDiagonalIndices(final DenseMatrix64F source, final int[] indices) {
++ public static void getFiniteDiagonalIndices(final DMatrixRMaj source, final int[] indices) {
+ final int length = source.getNumCols();
+
+ int index = 0;
+@@ -140,7 +140,7 @@
+ }
+ }
+
+- public static int countFiniteNonZeroDiagonals(DenseMatrix64F source) {
++ public static int countFiniteNonZeroDiagonals(DMatrixRMaj source) {
+ final int length = source.getNumCols();
+
+ int count = 0;
+@@ -153,7 +153,7 @@
+ return count;
+ }
+
+- public static void getFiniteNonZeroDiagonalIndices(final DenseMatrix64F source, final int[] indices) {
++ public static void getFiniteNonZeroDiagonalIndices(final DMatrixRMaj source, final int[] indices) {
+ final int length = source.getNumCols();
+
+ int index = 0;
+@@ -166,22 +166,22 @@
+ }
+ }
+
+- public static void addToDiagonal(DenseMatrix64F source, double increment) {
++ public static void addToDiagonal(DMatrixRMaj source, double increment) {
+ final int width = source.getNumRows();
+ for (int i = 0; i < width; ++i) {
+ source.unsafe_set(i,i, source.unsafe_get(i, i) + increment);
+ }
+ }
+
+- public static double det(DenseMatrix64F mat) {
++ public static double det(DMatrixRMaj mat) {
+ int numCol = mat.getNumCols();
+ int numRow = mat.getNumRows();
+ if(numCol != numRow) {
+ throw new IllegalArgumentException("Must be a square matrix.");
+ } else if(numCol <= 6) {
+- return numCol >= 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<DenseMatrix64F> solver = LinearSolverFactory.pseudoInverse(true);
++ LinearSolver<DMatrixRMaj, DMatrixRMaj> solver = LinearSolverFactory_DDRM.pseudoInverse(true);
+ solver.setA(source);
+
+- SingularValueDecomposition<DenseMatrix64F> svd = solver.getDecomposition();
++ SingularValueDecomposition_F64<DMatrixRMaj> 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<DenseMatrix64F> solver = LinearSolverFactory.pseudoInverse(true);
++ LinearSolver<DMatrixRMaj, DMatrixRMaj> solver = LinearSolverFactory_DDRM.pseudoInverse(true);
+ solver.setA(A);
+ solver.solve(B, X);
+
+@@ -306,7 +306,7 @@
+ double det = 1;
+
+ if (getDeterminant) {
+- SingularValueDecomposition<DenseMatrix64F> svd = solver.getDecomposition();
++ SingularValueDecomposition_F64<DMatrixRMaj> 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) {
View it on GitLab: https://salsa.debian.org/med-team/beast-mcmc/compare/0e8c48145b38c678c18a9d1cbd04eb7aa70df039...d5cab6e953b9b255fc9b227c1a34b0a1f026f1f4
--
View it on GitLab: https://salsa.debian.org/med-team/beast-mcmc/compare/0e8c48145b38c678c18a9d1cbd04eb7aa70df039...d5cab6e953b9b255fc9b227c1a34b0a1f026f1f4
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20190521/a5880b06/attachment-0001.html>
More information about the debian-med-commit
mailing list