add lphy and tests #406

examples/language/binaryCovarion.lphy

@@ -0,0 +1,15 @@
+// TODO better parameter values
+data {
+  L = 100;
+  n = 10;
+  birthRate = 10;
+}
+model {
+  alpha ~ Uniform(lower=1.0E-4, upper=1.0);
+  switchRate ~ Gamma(shape=0.05, scale=10.0);
+  // vfreq is the frequencies of the visible states, hfreq is the frequencies of the hidden states.
+  Q = binaryCovarion(alpha=alpha, s=switchRate, vfreq=[0.2, 0.8], hfreq=[0.5, 0.5]);
+
+  ψ ~ Yule(lambda=birthRate, n=n);
+  D ~ PhyloCTMC(tree=ψ, L=L, Q=Q, dataType=binaryDataType());
+}

lphy-base/src/main/java/lphy/base/evolution/likelihood/AbstractPhyloCTMC.java

@@ -79,12 +79,6 @@ public AbstractPhyloCTMC(Value<TimeTree> tree, Value<Number> clockRate, Value<Do
     // return Q matrix
     protected abstract Double[][] getQ();
 
-    // setup() before sample()
-    protected void setup() {
-        // overwrite the default if more setup
-        computePAndRootFreqs();
-    }
-
     // shared code in setup()
     protected void computePAndRootFreqs() {
         idMap.clear();
@@ -150,7 +144,13 @@ protected void traverseTree(TimeTreeNode node, int nodeState, SimpleAlignment al
         }
     }
 
-    //+++ getter +++//
+    //+++ public and getter +++//
+
+    // setup() before sample()
+    public void setup() {
+        // overwrite the default if more setup
+        computePAndRootFreqs();
+    }
 
     public Value<Double[]> getBranchRates() {
         return branchRates;
@@ -169,6 +169,33 @@ public SequenceType getDataType() {
         return dataType.value();
     }
 
+    // for unit test
+    public void getTransitionProbabilities(double branchLength, double[][] transProbs) {
+
+        int i, j, k;
+        double temp;
+
+        final int numStates = transProbs.length; // getQ().length ?
+        // inverse Eigen vectors
+        // Eigen values
+        for (i = 0; i < numStates; i++) {
+            temp = FastMath.exp(branchLength * Eval[i]);
+            for (j = 0; j < numStates; j++) {
+                iexp[i][j] = Ievc[i][j] * temp;
+            }
+        }
+
+        for (i = 0; i < numStates; i++) {
+            for (j = 0; j < numStates; j++) {
+                temp = 0.0;
+                for (k = 0; k < numStates; k++) {
+                    temp += Evec[i][k] * iexp[k][j];
+                }
+                transProbs[i][j] = FastMath.abs(temp);
+            }
+        }
+    }
+
     //+++ private methods +++//
 
     private Value<Double[]> computeEquilibrium(double[][] transProb) {
@@ -219,31 +246,6 @@ private int drawState(double[] p) {
         throw new RuntimeException("p vector should add to 1.0 but adds to " + totalP +  " instead.");
     }
 
-    private void getTransitionProbabilities(double branchLength, double[][] transProbs) {
-
-        int i, j, k;
-        double temp;
-
-        final int numStates = transProbs.length; // getQ().length ?
-        // inverse Eigen vectors
-        // Eigen values
-        for (i = 0; i < numStates; i++) {
-            temp = FastMath.exp(branchLength * Eval[i]);
-            for (j = 0; j < numStates; j++) {
-                iexp[i][j] = Ievc[i][j] * temp;
-            }
-        }
-
-        for (i = 0; i < numStates; i++) {
-            for (j = 0; j < numStates; j++) {
-                temp = 0.0;
-                for (k = 0; k < numStates; k++) {
-                    temp += Evec[i][k] * iexp[k][j];
-                }
-                transProbs[i][j] = FastMath.abs(temp);
-            }
-        }
-    }
 
     private static double EPSILON = 2.220446049250313E-16;
 

lphy-base/src/main/java/lphy/base/evolution/likelihood/PhyloCTMCSiteModel.java

@@ -110,7 +110,7 @@ public void setParam(String paramName, Value value) {
         }
     }
 
-    protected void setup() {
+    public void setup() {
 
         siteCount = getSiteCount();
 

lphy-base/src/main/java/lphy/base/evolution/substitutionmodel/BinaryCovarion.java

@@ -52,29 +52,13 @@ public Value<Double[][]> apply() {
         double[] vfreq = ValueUtils.doubleArrayValue(getParams().get(vfreqParamName));
         double[] hfreq = ValueUtils.doubleArrayValue(getParams().get(hfreqParamName));
 
-        Double[][] unnormalizedQ = setupUnnormalizedQMatrix(alpha, switchRate, vfreq, hfreq);
+        Double[][] Q = setupUnnormalizedQMatrix(alpha, switchRate, vfreq, hfreq);
         Double[] freqs = getFrequencies(vfreq, hfreq);
 
-        Double[][] Q = binaryCovarion(unnormalizedQ, freqs);
-
-        return new DoubleArray2DValue(Q, this);
-    }
-
-    private Double[][] binaryCovarion(Double[][] Q, Double[] freqs) {
-
-        // set up diagonal
-        for (int i = 0; i < NumOfStates; i++) {
-            double sum = 0.0;
-            for (int j = 0; j < NumOfStates; j++) {
-                if (i != j)
-                    sum += Q[i][j];
-            }
-            Q[i][i] = -sum;
-        }
         // normalise rate matrix to one expected substitution per unit time
         normalize(freqs, Q);
 
-        return Q;
+        return new DoubleArray2DValue(Q, this);
     }
 
     private Double[][] setupUnnormalizedQMatrix(double a, double s, double[] vf, double[] hf) {
@@ -104,6 +88,16 @@ private Double[][] setupUnnormalizedQMatrix(double a, double s, double[] vf, dou
         unnormalizedQ[3][1] = s;
         unnormalizedQ[3][2] = p0;
 
+        // set up diagonal
+        for (int i = 0; i < NumOfStates; i++) {
+            double sum = 0.0;
+            for (int j = 0; j < NumOfStates; j++) {
+                if (i != j)
+                    sum += unnormalizedQ[i][j];
+            }
+            unnormalizedQ[i][i] = -sum;
+        }
+
         return unnormalizedQ;
     }
 
@@ -138,7 +132,7 @@ void normalize(Double[] freqs, Double[][] Q) {
         }
     }
 
-    private Double[] getFrequencies(double[] vf, double[] hf) {
+    public Double[] getFrequencies(double[] vf, double[] hf) {
         Double[] freqs = new Double[NumOfStates];
         freqs[0] = vf[0] * hf[0];
         freqs[1] = vf[1] * hf[0];

lphy-base/src/test/java/lphy/base/evolution/substitutionmodel/BinaryCovarionTest.java

@@ -0,0 +1,136 @@
+package lphy.base.evolution.substitutionmodel;
+
+import lphy.base.evolution.coalescent.Coalescent;
+import lphy.base.evolution.likelihood.PhyloCTMC;
+import lphy.base.evolution.tree.TimeTree;
+import lphy.core.model.RandomVariable;
+import lphy.core.model.Value;
+import lphy.core.model.ValueUtils;
+import lphy.core.simulator.RandomUtils;
+import org.apache.commons.lang3.ArrayUtils;
+import org.junit.jupiter.api.BeforeEach;
+import org.junit.jupiter.api.Test;
+
+import java.util.Arrays;
+
+import static org.junit.jupiter.api.Assertions.assertArrayEquals;
+import static org.junit.jupiter.api.Assertions.assertEquals;
+
+/**
+ * mode = BEAST -- using classic BEAST implementation, reversible iff hfrequencies = (0.5, 0.5)
+ * FLAGS: reversible = false, TSParameterisation = false
+ *
+ * [ -(a*p1)-s ,   a*p1    ,    s   ,   0   ]
+ * [   a*p0    , -(a*p0)-s ,    0   ,   s   ]
+ * [    s      ,     0     ,  -p1-s ,  p1   ]
+ * [    0      ,     s     ,    p0  , -p0-s ]
+ *
+ * equilibrium frequencies
+ * [ p0 * f0, p1, * f0, p0 * f1, p1, * f1 ]
+ */
+class BinaryCovarionTest {
+
+    RandomVariable<TimeTree> tree;
+
+    Value<Number> alpha = new Value<>("alpha", 0.01);
+    Value<Number> s = new Value<>("switchRate", 0.1);
+
+    @BeforeEach
+    void setUp() {
+        // tree
+        Coalescent coalescent = new Coalescent(new Value<>("theta", 10.0),
+                new Value<>("n", 10), null);
+        tree = coalescent.sample();
+        // only tip labels are used in the follow code
+    }
+
+    @Test
+    void testQAndFreqs() {
+        Value<Number[]> vfreq = new Value<>("fVisibleStates", new Number[]{0.2, 0.8});
+        Value<Number[]> hfreq = new Value<>("fHiddenStates", new Number[]{0.5, 0.5});
+
+        BinaryCovarion binaryCovarion = new BinaryCovarion(alpha, s, vfreq, hfreq,
+                new Value<>("meanRate", 1.0));
+        Value<Double[][]> Q = binaryCovarion.apply();
+        System.out.println("Q matrix: " + Arrays.deepToString(Q.value()));
+
+        // 1st row of Q matrix
+        assertArrayEquals(new double[]{-0.6683, 0.0495, 0.6188, 0.0},
+                ArrayUtils.toPrimitive(Q.value()[0]), 1e-3);
+
+        // test freqs
+        Double[] freqs = binaryCovarion.getFrequencies(ValueUtils.doubleArrayValue(vfreq),
+                ValueUtils.doubleArrayValue(hfreq));
+        assertArrayEquals(new double[]{0.1, 0.4, 0.1, 0.4}, ArrayUtils.toPrimitive(freqs));
+    }
+
+
+    @Test
+    void testTransProbWithEqualHFreqs() {
+        double d = 0.05 + RandomUtils.getRandom().nextDouble() * 0.9;
+        Value<Number[]> vfreq = new Value<>("fVisibleStates", new Number[]{d, 1.0 - d});
+        Value<Number[]> hfreq = new Value<>("fHiddenStates", new Number[]{0.5, 0.5});
+
+        BinaryCovarion binaryCovarion = new BinaryCovarion(alpha, s, vfreq, hfreq,
+                new Value<>("meanRate", 1.0));
+        Value<Double[][]> Q = binaryCovarion.apply();
+        System.out.println("Q matrix: " + Arrays.deepToString(Q.value()));
+
+        PhyloCTMC phyloCTMC = new PhyloCTMC(tree, null, null, Q,
+                null, null, new Value<>("L", 100), null, null);
+        // not simulate sequences, so ignore sample(), then tree stats and L should not affect the result
+        phyloCTMC.setup();
+
+        double[][] p = new double[4][4];
+        // branchLength = 1000
+        phyloCTMC.getTransitionProbabilities(100, p);
+        System.out.println("Actual: " + Arrays.deepToString(p));
+
+        // freqs
+        Double[] freqs = binaryCovarion.getFrequencies(ValueUtils.doubleArrayValue(vfreq),
+                ValueUtils.doubleArrayValue(hfreq));
+        System.out.println("Expected: " + Arrays.toString(freqs));
+
+        // 1st row of p matrix
+        for (int j = 0; j < 4; j++) {
+            assertEquals(freqs[j], p[0][j], 1e-3);
+        }
+
+    }
+
+//TODO not working
+//    @Test
+//    void testTransProbWithUnEqualHFreqs() {
+//        double d = 0.05 + RandomUtils.getRandom().nextDouble() * 0.9;
+//        Value<Number[]> vfreq = new Value<>("fVisibleStates", new Number[]{d, 1.0 - d});
+//        d = 0.05 + RandomUtils.getRandom().nextDouble() * 0.9;
+//        Value<Number[]> hfreq = new Value<>("fHiddenStates", new Number[]{d, 1.0 - d});
+//
+//        BinaryCovarion binaryCovarion = new BinaryCovarion(alpha, s, vfreq, hfreq,
+//                new Value<>("meanRate", 1.0));
+//        Value<Double[][]> Q = binaryCovarion.apply();
+//        System.out.println("Q matrix: " + Arrays.deepToString(Q.value()));
+//
+//        PhyloCTMC phyloCTMC = new PhyloCTMC(tree, null, null, Q,
+//                null, null, new Value<>("L", 100), null, null);
+//        // not simulate sequences, so ignore sample(), then tree stats and L should not affect the result
+//        phyloCTMC.setup();
+//
+//        double[][] p = new double[4][4];
+//        // branchLength = 1000
+//        phyloCTMC.getTransitionProbabilities(100, p);
+//        System.out.println("Actual: " + Arrays.deepToString(p));
+//
+//        // freqs
+//        Double[] freqs = binaryCovarion.getFrequencies(ValueUtils.doubleArrayValue(vfreq),
+//                ValueUtils.doubleArrayValue(hfreq));
+//        System.out.println("Expected: " + Arrays.toString(freqs));
+//
+//        // 1st row of p matrix
+//        for (int j = 0; j < 4; j++) {
+//            assertEquals(freqs[j], p[0][j], 1e-3);
+//        }
+//
+//    }
+
+}