diff --git a/src/test/substmodels/nucleotide/TIM1Test.java b/src/test/substmodels/nucleotide/TIM1Test.java new file mode 100644 index 0000000..755141d --- /dev/null +++ b/src/test/substmodels/nucleotide/TIM1Test.java @@ -0,0 +1,82 @@ +package test.substmodels.nucleotide; + +import beast.core.Description; +import beast.core.parameter.RealParameter; +import junit.framework.TestCase; +import substmodels.nucleotide.TIM1; + +@Description("Test TIM1 matrix exponentiation") +public class TIM1Test extends TestCase { + + /* + * import numpy as np + * from scipy.linalg import expm + * + * piQ = np.diag([.25, .25, .25, .25]) + * d = 0.3 + * # Q matrix with zeroed diagonal + * XQ = np.matrix([[0, 1, 2, 3], [1, 0, 3, 4], [2, 3, 0, 1], [3, 4, 1, 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[] {1.0, 2.0, 3.0, 4.0}; + } + + @Override + public double getDistance() { + return 0.3; + } + + @Override + public double[] getExpectedResult() { + return new double[]{ + 0.783537640428, 0.044577861463, 0.071332157416, 0.100552340692, + 0.044577861463, 0.7275631612 , 0.100552340692, 0.127306636645, + 0.071332157416, 0.100552340692, 0.783537640428, 0.044577861463, + 0.100552340692, 0.127306636645, 0.044577861463, 0.7275631612 + }; + } + }; + + + EqualBaseFrequencies[] all = {test0}; + + public void testTIM1() throws Exception { + for (EqualBaseFrequencies test : all) { + + TIM1 tim1 = new TIM1(); + RealParameter rates = new RealParameter(test.getRates()); + tim1.initByName("rates", rates); + tim1.printQ(System.out); // to obtain XQ for python script +// for (int i = 0; i < 6; ++i) +// System.out.println("Rate " + tim1.getSubstitution(i) + " : " + tim1.getRate(i)); + + // AC=GT + assertEquals(true, tim1.getRateAC()== tim1.getRateGT() && + tim1.getRateAC()!=tim1.getRateAT() && tim1.getRateAC()!=tim1.getRateAG()); + // AT=CG + assertEquals(true, tim1.getRateAT()== tim1.getRateCG() && + tim1.getRateAT()!=tim1.getRateAC() && tim1.getRateAT()!=tim1.getRateAG()); + // AG!=CT + assertEquals(true, tim1.getRateAC()!=tim1.getRateAG() && + tim1.getRateAC()!=tim1.getRateCT() && tim1.getRateAG()!=tim1.getRateCT()); + + double distance = test.getDistance(); + double[] mat = new double[4 * 4]; + tim1.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/TIM1ufTest.java b/src/test/substmodels/nucleotide/TIM1ufTest.java new file mode 100644 index 0000000..533f6f7 --- /dev/null +++ b/src/test/substmodels/nucleotide/TIM1ufTest.java @@ -0,0 +1,93 @@ +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.TIM1uf; + +@Description("Test TIM1uf matrix exponentiation") +public class TIM1ufTest extends TestCase { + + /* + * import numpy as np + * from scipy.linalg import expm + * + * piQ = np.diag([0.4, 0.3, 0.2, 0.1]) + * d = 0.3 + * # Q matrix with zeroed diagonal + * XQ = np.matrix([[0, 1, 2, 3], [1, 0, 3, 4], [2, 3, 0, 1], [3, 4, 1, 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[] {1.0, 2.0, 3.0, 4.0}; + } + + @Override + public double getDistance() { + return 0.3; + } + + @Override + public double[] getExpectedResult() { + return new double[]{ + 0.825871651293, 0.060874601338, 0.067043906735, 0.046209840634, + 0.081166135118, 0.765617592587, 0.094589810557, 0.058626461738, + 0.134087813471, 0.141884715835, 0.701565807625, 0.022461663069, + 0.184839362534, 0.175879385215, 0.044923326138, 0.594357926113 + }; + } + }; + + + UnequalBaseFrequencies[] all = {test0}; + + public void testTIM1uf() throws Exception { + for (UnequalBaseFrequencies test : all) { + + RealParameter f = new RealParameter(test.getPi()); + Frequencies freqs = new Frequencies(); + freqs.initByName("frequencies", f); // "estimate", true + + TIM1uf tim1uf = new TIM1uf(); + RealParameter rates = new RealParameter(test.getRates()); + tim1uf.initByName("rates", rates, "frequencies", freqs); + tim1uf.printQ(System.out); // to obtain XQ for python script +// for (int i = 0; i < 6; ++i) +// System.out.println("Rate " + tim1uf.getSubstitution(i) + " : " + tim1uf.getRate(i)); + + // AC=GT + assertEquals(true, tim1uf.getRateAC()== tim1uf.getRateGT() && + tim1uf.getRateAC()!=tim1uf.getRateAT() && tim1uf.getRateAC()!=tim1uf.getRateAG()); + // AT=CG + assertEquals(true, tim1uf.getRateAT()== tim1uf.getRateCG() && + tim1uf.getRateAT()!=tim1uf.getRateAC() && tim1uf.getRateAT()!=tim1uf.getRateAG()); + // AG!=CT + assertEquals(true, tim1uf.getRateAC()!=tim1uf.getRateAG() && + tim1uf.getRateAC()!=tim1uf.getRateCT() && tim1uf.getRateAG()!=tim1uf.getRateCT()); + + + double distance = test.getDistance(); + double[] mat = new double[4 * 4]; + tim1uf.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/TIM2Test.java b/src/test/substmodels/nucleotide/TIM2Test.java new file mode 100644 index 0000000..5561ccd --- /dev/null +++ b/src/test/substmodels/nucleotide/TIM2Test.java @@ -0,0 +1,82 @@ +package test.substmodels.nucleotide; + +import beast.core.Description; +import beast.core.parameter.RealParameter; +import junit.framework.TestCase; +import substmodels.nucleotide.TIM2; + +@Description("Test TIM2 matrix exponentiation") +public class TIM2Test extends TestCase { + + /* + * import numpy as np + * from scipy.linalg import expm + * + * piQ = np.diag([.25, .25, .25, .25]) + * d = 0.3 + * # Q matrix with zeroed diagonal + * XQ = np.matrix([[0, 1, 2, 1], [1, 0, 3, 4], [2, 3, 0, 3], [1, 4, 3, 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[] {1.0, 2.0, 3.0, 4.0}; + } + + @Override + public double getDistance() { + return 0.3; + } + + @Override + public double[] getExpectedResult() { + return new double[]{ + 0.847204213099, 0.040767943097, 0.071259900708, 0.040767943097, + 0.040767943097, 0.727703865901, 0.101751858319, 0.129776332684, + 0.071259900708, 0.101751858319, 0.725236382655, 0.101751858319, + 0.040767943097, 0.129776332684, 0.101751858319, 0.727703865901 + }; + } + }; + + + EqualBaseFrequencies[] all = {test0}; + + public void testTIM2() throws Exception { + for (EqualBaseFrequencies test : all) { + + TIM2 tim2 = new TIM2(); + RealParameter rates = new RealParameter(test.getRates()); + tim2.initByName("rates", rates); + tim2.printQ(System.out); // to obtain XQ for python script +// for (int i = 0; i < 6; ++i) +// System.out.println("Rate " + tim2.getSubstitution(i) + " : " + tim2.getRate(i)); + + // AC=AT + assertEquals(true, tim2.getRateAC()== tim2.getRateAT() && + tim2.getRateAC()!=tim2.getRateCG() && tim2.getRateAC()!=tim2.getRateAG()); + // CG=GT + assertEquals(true, tim2.getRateCG()== tim2.getRateGT() && + tim2.getRateCG()!=tim2.getRateCT() && tim2.getRateCG()!=tim2.getRateAG()); + // AG!=CT + assertEquals(true, tim2.getRateAC()!=tim2.getRateAG() && + tim2.getRateAC()!=tim2.getRateCT() && tim2.getRateAG()!=tim2.getRateCT()); + + double distance = test.getDistance(); + double[] mat = new double[4 * 4]; + tim2.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/TIM2ufTest.java b/src/test/substmodels/nucleotide/TIM2ufTest.java new file mode 100644 index 0000000..8c2ea71 --- /dev/null +++ b/src/test/substmodels/nucleotide/TIM2ufTest.java @@ -0,0 +1,93 @@ +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.TIM2uf; + +@Description("Test TIM2uf matrix exponentiation") +public class TIM2ufTest extends TestCase { + + /* + * import numpy as np + * from scipy.linalg import expm + * + * piQ = np.diag([0.4, 0.3, 0.2, 0.1]) + * d = 0.3 + * # Q matrix with zeroed diagonal + * XQ = np.matrix([[0, 1, 2, 1], [1, 0, 3, 4], [2, 3, 0, 3], [1, 4, 3, 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[] {1.0, 2.0, 3.0, 4.0}; + } + + @Override + public double getDistance() { + return 0.3; + } + + @Override + public double[] getExpectedResult() { + return new double[]{ + 0.848004762603, 0.061442861874, 0.070071421565, 0.020480953958, + 0.081923815832, 0.755054682896, 0.099180935215, 0.063840566057, + 0.140142843131, 0.148771402822, 0.661495286439, 0.049590467607, + 0.081923815832, 0.19152169817 , 0.099180935215, 0.627373550783 + }; + } + }; + + + UnequalBaseFrequencies[] all = {test0}; + + public void testTIM2uf() throws Exception { + for (UnequalBaseFrequencies test : all) { + + RealParameter f = new RealParameter(test.getPi()); + Frequencies freqs = new Frequencies(); + freqs.initByName("frequencies", f); // "estimate", true + + TIM2uf tim2uf = new TIM2uf(); + RealParameter rates = new RealParameter(test.getRates()); + tim2uf.initByName("rates", rates, "frequencies", freqs); + tim2uf.printQ(System.out); // to obtain XQ for python script +// for (int i = 0; i < 6; ++i) +// System.out.println("Rate " + tim2uf.getSubstitution(i) + " : " + tim2uf.getRate(i)); + + // AC=AT + assertEquals(true, tim2uf.getRateAC()== tim2uf.getRateAT() && + tim2uf.getRateAC()!=tim2uf.getRateCG() && tim2uf.getRateAC()!=tim2uf.getRateAG()); + // CG=GT + assertEquals(true, tim2uf.getRateCG()== tim2uf.getRateGT() && + tim2uf.getRateCG()!=tim2uf.getRateCT() && tim2uf.getRateCG()!=tim2uf.getRateAG()); + // AG!=CT + assertEquals(true, tim2uf.getRateAC()!=tim2uf.getRateAG() && + tim2uf.getRateAC()!=tim2uf.getRateCT() && tim2uf.getRateAG()!=tim2uf.getRateCT()); + + + double distance = test.getDistance(); + double[] mat = new double[4 * 4]; + tim2uf.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/TIM3Test.java b/src/test/substmodels/nucleotide/TIM3Test.java new file mode 100644 index 0000000..eff38b6 --- /dev/null +++ b/src/test/substmodels/nucleotide/TIM3Test.java @@ -0,0 +1,82 @@ +package test.substmodels.nucleotide; + +import beast.core.Description; +import beast.core.parameter.RealParameter; +import junit.framework.TestCase; +import substmodels.nucleotide.TIM3; + +@Description("Test TIM3 matrix exponentiation") +public class TIM3Test extends TestCase { + + /* + * import numpy as np + * from scipy.linalg import expm + * + * piQ = np.diag([.25, .25, .25, .25]) + * d = 0.3 + * # Q matrix with zeroed diagonal + * XQ = np.matrix([[0, 1, 2, 3], [1, 0, 1, 4], [2, 1, 0, 3], [3, 4, 3, 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[] {1.0, 2.0, 3.0, 4.0}; + } + + @Override + public double getDistance() { + return 0.3; + } + + @Override + public double[] getExpectedResult() { + return new double[]{ + 0.78353838664 , 0.043309922451, 0.073798790951, 0.099352899958, + 0.043309922451, 0.786005766388, 0.043309922451, 0.127374388711, + 0.073798790951, 0.043309922451, 0.78353838664 , 0.099352899958, + 0.099352899958, 0.127374388711, 0.099352899958, 0.673919811374 + }; + } + }; + + + EqualBaseFrequencies[] all = {test0}; + + public void testTIM3() throws Exception { + for (EqualBaseFrequencies test : all) { + + TIM3 tim3 = new TIM3(); + RealParameter rates = new RealParameter(test.getRates()); + tim3.initByName("rates", rates); + tim3.printQ(System.out); // to obtain XQ for python script +// for (int i = 0; i < 6; ++i) +// System.out.println("Rate " + tim3.getSubstitution(i) + " : " + tim3.getRate(i)); + + // AC=CG + assertEquals(true, tim3.getRateAC()== tim3.getRateCG() && + tim3.getRateAC()!=tim3.getRateAT() && tim3.getRateAC()!=tim3.getRateAG()); + // AT=GT + assertEquals(true, tim3.getRateAT()== tim3.getRateGT() && + tim3.getRateAT()!=tim3.getRateCT() && tim3.getRateCG()!=tim3.getRateAG()); + // AG!=CT + assertEquals(true, tim3.getRateAC()!=tim3.getRateAG() && + tim3.getRateAC()!=tim3.getRateCT() && tim3.getRateAG()!=tim3.getRateCT()); + + double distance = test.getDistance(); + double[] mat = new double[4 * 4]; + tim3.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/TIM3ufTest.java b/src/test/substmodels/nucleotide/TIM3ufTest.java new file mode 100644 index 0000000..fb1dcc0 --- /dev/null +++ b/src/test/substmodels/nucleotide/TIM3ufTest.java @@ -0,0 +1,92 @@ +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.TIM3uf; + +@Description("Test TIM3uf matrix exponentiation") +public class TIM3ufTest extends TestCase { + + /* + * import numpy as np + * from scipy.linalg import expm + * + * piQ = np.diag([0.4, 0.3, 0.2, 0.1]) + * d = 0.3 + * # Q matrix with zeroed diagonal + * XQ = np.matrix([[0, 1, 2, 3], [1, 0, 1, 4], [2, 1, 0, 3], [3, 4, 3, 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[] {1.0, 2.0, 3.0, 4.0}; + } + + @Override + public double getDistance() { + return 0.3; + } + + @Override + public double[] getExpectedResult() { + return new double[]{ + 0.808287239489, 0.065896148199, 0.076235614109, 0.049580998203, + 0.087861530931, 0.804818897665, 0.043930765466, 0.063388805938, + 0.152471228218, 0.065896148199, 0.73205162538 , 0.049580998203, + 0.198323992812, 0.190166417815, 0.099161996406, 0.512347592967 + }; + } + }; + + + UnequalBaseFrequencies[] all = {test0}; + + public void testTIM3uf() throws Exception { + for (UnequalBaseFrequencies test : all) { + + RealParameter f = new RealParameter(test.getPi()); + Frequencies freqs = new Frequencies(); + freqs.initByName("frequencies", f); // "estimate", true + + TIM3uf tim3uf = new TIM3uf(); + RealParameter rates = new RealParameter(test.getRates()); + tim3uf.initByName("rates", rates, "frequencies", freqs); + tim3uf.printQ(System.out); // to obtain XQ for python script +// for (int i = 0; i < 6; ++i) +// System.out.println("Rate " + tim3uf.getSubstitution(i) + " : " + tim3uf.getRate(i)); + + // AC=CG + assertEquals(true, tim3uf.getRateAC()== tim3uf.getRateCG() && + tim3uf.getRateAC()!=tim3uf.getRateAT() && tim3uf.getRateAC()!=tim3uf.getRateAG()); + // AT=GT + assertEquals(true, tim3uf.getRateAT()== tim3uf.getRateGT() && + tim3uf.getRateAT()!=tim3uf.getRateCT() && tim3uf.getRateCG()!=tim3uf.getRateAG()); + // AG!=CT + assertEquals(true, tim3uf.getRateAC()!=tim3uf.getRateAG() && + tim3uf.getRateAC()!=tim3uf.getRateCT() && tim3uf.getRateAG()!=tim3uf.getRateCT()); + + double distance = test.getDistance(); + double[] mat = new double[4 * 4]; + tim3uf.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