diff --git a/src/test/substmodels/nucleotide/SYMTest.java b/src/test/substmodels/nucleotide/SYMTest.java new file mode 100644 index 0000000..0392b47 --- /dev/null +++ b/src/test/substmodels/nucleotide/SYMTest.java @@ -0,0 +1,89 @@ +package test.substmodels.nucleotide; + + +import beast.core.Description; +import beast.core.parameter.RealParameter; +import junit.framework.TestCase; +import substmodels.nucleotide.SYM; + +/** + * Test SYM matrix exponentiation + * + */ +@Description("Test SYM matrix exponentiation") +public class SYMTest extends TestCase { + + /* + * import numpy as np + * from scipy.linalg import expm + * + * piQ = np.diag([.25, .25, .25, .25]) + * d = 0.1 + * # Q matrix with zeroed diagonal + * XQ = np.matrix([[0, 0.2, 10, .3], [0.2, 0, 0.4, 5], [10, 0.4, 0, 0.5], [0.3, 5, 0.5, 0]]) + * + * xx = XQ * piQ + * + * # fill diagonal and normalize by total substitution rate + * q0 = (xx + np.diag(np.squeeze(np.asarray(-np.sum(xx, axis=1))))) / np.sum(piQ * np.sum(xx, axis=1)) + * expm(q0 * d) + */ + protected EqualBaseFrequencies test0 = new EqualBaseFrequencies() { + @Override + public Double [] getRates() { + return new Double[] {0.2, 10.0, 0.3, 0.4, 5.0, 0.5}; + } + + @Override + public double getDistance() { + return 0.1; + } + + @Override + public double[] getExpectedResult() { + return new double[]{ + 0.886360401447, 0.002594129576, 0.107315348219, 0.003730120758, + 0.002594129576, 0.93573730447 , 0.004733198723, 0.056935367232, + 0.107315348219, 0.004733198723, 0.882087475436, 0.005863977622, + 0.003730120758, 0.056935367232, 0.005863977622, 0.933470534387 + }; + } + }; + + + EqualBaseFrequencies[] all = {test0}; + + public void testSYM() throws Exception { + for (EqualBaseFrequencies test : all) { + + SYM sym = new SYM(); + RealParameter symRates = new RealParameter(test.getRates()); + sym.initByName("rates", symRates); + sym.printQ(System.out); // to obtain XQ for python script +// for (int i = 0; i < 6; ++i) +// System.out.println("Rate " + sym.getSubstitution(i) + " : " + sym.getRate(i)); + + assertEquals(false, sym.getRateAC()== sym.getRateAG()); + assertEquals(false, sym.getRateAC()== sym.getRateAT()); + assertEquals(false, sym.getRateAC()== sym.getRateCG()); + assertEquals(false, sym.getRateAC()== sym.getRateCT()); + assertEquals(false, sym.getRateAC()== sym.getRateGT()); + assertEquals(false, sym.getRateAG()== sym.getRateAT()); + assertEquals(false, sym.getRateAG()== sym.getRateCG()); + assertEquals(false, sym.getRateAG()== sym.getRateCT()); + assertEquals(false, sym.getRateAG()== sym.getRateGT()); + assertEquals(false, sym.getRateCG()== sym.getRateGT()); + assertEquals(false, sym.getRateCT()== sym.getRateGT()); + + double distance = test.getDistance(); + double[] mat = new double[4 * 4]; + sym.getTransitionProbabilities(null, distance, 0, 1, mat); + + final double[] result = test.getExpectedResult(); + for (int k = 0; k < mat.length; ++k) { + assertEquals(mat[k], result[k], 1e-10); + System.out.println(k + " : " + (mat[k] - result[k])); + } + } + } +} \ No newline at end of file diff --git a/src/test/substmodels/nucleotide/TVMTest.java b/src/test/substmodels/nucleotide/TVMTest.java new file mode 100644 index 0000000..1656052 --- /dev/null +++ b/src/test/substmodels/nucleotide/TVMTest.java @@ -0,0 +1,100 @@ +package test.substmodels.nucleotide; + + +import beast.core.Description; +import beast.core.parameter.RealParameter; +import beast.evolution.substitutionmodel.Frequencies; +import junit.framework.TestCase; +import substmodels.nucleotide.TVM; + +/** + * Test TVM matrix exponentiation + * + */ +@Description("Test TVM matrix exponentiation") +public class TVMTest extends TestCase { + + /* + * import numpy as np + * from scipy.linalg import expm + * + * piQ = np.diag([0.4, 0.3, 0.2, 0.1]) + * d = 0.1 + * # Q matrix with zeroed diagonal + * XQ = np.matrix([[0, 0.2, 10, .3], [0.2, 0, 0.4, 10], [10, 0.4, 0, 5], [0.3, 10, 5, 0]]) + * + * xx = XQ * piQ + * + * # fill diagonal and normalize by total substitution rate + * q0 = (xx + np.diag(np.squeeze(np.asarray(-np.sum(xx, axis=1))))) / np.sum(piQ * np.sum(xx, axis=1)) + * expm(q0 * d) + */ + protected UnequalBaseFrequencies test0 = new UnequalBaseFrequencies() { + @Override + public Double[] getPi() { + return new Double[]{0.4, 0.3, 0.2, 0.1}; + } + + @Override + public Double [] getRates() { + return new Double[] {0.2, 10.0, 0.3, 0.4, 5.0}; + } + + @Override + public double getDistance() { + return 0.1; + } + + @Override + public double[] getExpectedResult() { + return new double[]{ + 0.926032879344, 0.002501473304, 0.069681075145, 0.001784572207, + 0.003335297739, 0.957195914822, 0.00364943052 , 0.035819356919, + 0.13936215029 , 0.005474145781, 0.838262134317, 0.016901569612, + 0.00713828883 , 0.107458070758, 0.033803139225, 0.851600501187 + }; + } + }; + + + UnequalBaseFrequencies[] all = {test0}; + + public void testTVM() throws Exception { + for (UnequalBaseFrequencies test : all) { + + RealParameter f = new RealParameter(test.getPi()); + Frequencies freqs = new Frequencies(); + freqs.initByName("frequencies", f); // "estimate", true + + TVM tvm = new TVM(); + RealParameter tvmRates = new RealParameter(test.getRates()); + tvm.initByName("rates", tvmRates, "frequencies", freqs); + tvm.printQ(System.out); // to obtain XQ for python script +// for (int i = 0; i < 6; ++i) +// System.out.println("Rate " + tvm.getSubstitution(i) + " : " + tvm.getRate(i)); + + assertEquals(false, tvm.getRateAC()== tvm.getRateAG()); + assertEquals(false, tvm.getRateAC()== tvm.getRateAT()); + assertEquals(false, tvm.getRateAC()== tvm.getRateCG()); + assertEquals(false, tvm.getRateAC()== tvm.getRateCT()); + assertEquals(false, tvm.getRateAC()== tvm.getRateGT()); + assertEquals(false, tvm.getRateAG()== tvm.getRateAT()); + assertEquals(false, tvm.getRateAG()== tvm.getRateCG()); + assertEquals(false, tvm.getRateAG()== tvm.getRateGT()); + assertEquals(false, tvm.getRateCG()== tvm.getRateGT()); + assertEquals(false, tvm.getRateCT()== tvm.getRateGT()); + // AG=CT + assertEquals(true, tvm.getRateAG()== tvm.getRateCT()); + + double distance = test.getDistance(); + double[] mat = new double[4 * 4]; + tvm.getTransitionProbabilities(null, distance, 0, 1, mat); + + final double[] result = test.getExpectedResult(); + for (int k = 0; k < mat.length; ++k) { + assertEquals(mat[k], result[k], 1e-10); + System.out.println(k + " : " + (mat[k] - result[k])); + } + } + } +} \ No newline at end of file diff --git a/src/test/substmodels/nucleotide/TVMefTest.java b/src/test/substmodels/nucleotide/TVMefTest.java new file mode 100644 index 0000000..2c63814 --- /dev/null +++ b/src/test/substmodels/nucleotide/TVMefTest.java @@ -0,0 +1,90 @@ +package test.substmodels.nucleotide; + + +import beast.core.Description; +import beast.core.parameter.RealParameter; +import junit.framework.TestCase; +import substmodels.nucleotide.TVMef; + +/** + * Test TVMef matrix exponentiation + * + */ +@Description("Test TVMef matrix exponentiation") +public class TVMefTest extends TestCase { + + /* + * import numpy as np + * from scipy.linalg import expm + * + * piQ = np.diag([.25, .25, .25, .25]) + * d = 0.1 + * # Q matrix with zeroed diagonal + * XQ = np.matrix([[0, 0.2, 10, .3], [0.2, 0, 0.4, 10], [10, 0.4, 0, 5], [0.3, 10, 5, 0]]) + * + * xx = XQ * piQ + * + * # fill diagonal and normalize by total substitution rate + * q0 = (xx + np.diag(np.squeeze(np.asarray(-np.sum(xx, axis=1))))) / np.sum(piQ * np.sum(xx, axis=1)) + * expm(q0 * d) + */ + protected EqualBaseFrequencies test0 = new EqualBaseFrequencies() { + @Override + public Double [] getRates() { + return new Double[] {0.2, 10.0, 0.3, 0.4, 5.0}; + } + + @Override + public double getDistance() { + return 0.1; + } + + @Override + public double[] getExpectedResult() { + return new double[]{ + 0.924841386585, 0.001651397824, 0.070005984376, 0.003501231215, + 0.001651397824, 0.924131486211, 0.004198599114, 0.07001851685, + 0.070005984376, 0.004198599114, 0.891234320911, 0.034561095599, + 0.003501231215, 0.07001851685 , 0.034561095599, 0.891919156335 + }; + } + }; + + + EqualBaseFrequencies[] all = {test0}; + + public void testTVMef() throws Exception { + for (EqualBaseFrequencies test : all) { + + TVMef tvmef = new TVMef(); + RealParameter tvmefRates = new RealParameter(test.getRates()); + tvmef.initByName("rates", tvmefRates); + tvmef.printQ(System.out); // to obtain XQ for python script +// for (int i = 0; i < 6; ++i) +// System.out.println("Rate " + tvmef.getSubstitution(i) + " : " + tvmef.getRate(i)); + + assertEquals(false, tvmef.getRateAC()== tvmef.getRateAG()); + assertEquals(false, tvmef.getRateAC()== tvmef.getRateAT()); + assertEquals(false, tvmef.getRateAC()== tvmef.getRateCG()); + assertEquals(false, tvmef.getRateAC()== tvmef.getRateCT()); + assertEquals(false, tvmef.getRateAC()== tvmef.getRateGT()); + assertEquals(false, tvmef.getRateAG()== tvmef.getRateAT()); + assertEquals(false, tvmef.getRateAG()== tvmef.getRateCG()); + assertEquals(false, tvmef.getRateAG()== tvmef.getRateGT()); + assertEquals(false, tvmef.getRateCG()== tvmef.getRateGT()); + assertEquals(false, tvmef.getRateCT()== tvmef.getRateGT()); + // AG=CT + assertEquals(true, tvmef.getRateAG()== tvmef.getRateCT()); + + double distance = test.getDistance(); + double[] mat = new double[4 * 4]; + tvmef.getTransitionProbabilities(null, distance, 0, 1, mat); + + final double[] result = test.getExpectedResult(); + for (int k = 0; k < mat.length; ++k) { + assertEquals(mat[k], result[k], 1e-10); + System.out.println(k + " : " + (mat[k] - result[k])); + } + } + } +} \ No newline at end of file