package test.dr.evomodel.substmodel; import dr.evolution.datatype.OldHiddenNucleotides; import dr.evomodel.substmodel.*; import dr.inference.model.Parameter; import dr.oldevomodel.substmodel.CovarionHKY; import dr.oldevomodel.substmodel.FrequencyModel; import dr.oldevomodel.substmodel.SubstitutionModelUtils; import junit.framework.Test; import junit.framework.TestCase; import junit.framework.TestSuite; /** * CovarionHKY Tester. * * @author Alexei Drummond * @version 1.0 * @since <pre>08/20/2007</pre> */ public class CovarionHKYTest extends TestCase { public CovarionHKYTest(String name) { super(name); } public void setUp() throws Exception { super.setUp(); alpha = new Parameter.Default(0.0); switchingRate = new Parameter.Default(1.0); kappa = new Parameter.Default(2.0); dataType = new OldHiddenNucleotides(2); } public void tearDown() throws Exception { super.tearDown(); } public void testK2PTransitionProbabilities() { double[] frequencies = new double[] {0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125}; transitionProbabilitiesTester(frequencies, matLabPChange_K2P); } public void testHKYTransitionProbabilities() { double[] baseFrequencies = new double[]{0.15, 0.2, 0.3, 0.35}; double[] frequencies = new double[8]; for (int i = 0; i < frequencies.length; i++) { frequencies[i] = baseFrequencies[i % 4] / 2.0; } transitionProbabilitiesTester(frequencies, matLabPChange_HKY); } public void transitionProbabilitiesTester(double[] freqs, double[] expectedP) { Parameter frequencies = new Parameter.Default(freqs); FrequencyModel freqModel = new FrequencyModel(dataType, frequencies); model = new CovarionHKY(dataType, kappa, alpha, switchingRate, freqModel); alpha.setParameterValue(0, 0.0); switchingRate.setParameterValue(0, 1.0); kappa.setParameterValue(0, 2.0); model.setupMatrix(); System.out.println(SubstitutionModelUtils.toString(model.getQ(), dataType, 2)); double[] matrix = new double[64]; double[] pi = model.getFrequencyModel().getFrequencies(); int index = 0; for (double distance = 0.01; distance <= 1.005; distance += 0.01) { model.getTransitionProbabilities(distance, matrix); double pChange = 0.0; double pNoOrHiddenChange = 0.0; int k = 0; for (int i = 0; i < 8; i++) { for (int j = 0; j < 8; j++) { if ((i % 4) != (j % 4)) { pChange += matrix[k] * pi[i]; } else { pNoOrHiddenChange += matrix[k] * pi[i]; } k += 1; } } double totalP = pChange + pNoOrHiddenChange; assertEquals(1.0, totalP, 1e-10); System.out.print(distance + "\t" + "\t" + pChange + "\t"); if (index < 100) System.out.print(expectedP[index]); System.out.println(); //assertEquals(expectedP[index], pChange, 1e-14); index += 1; } } public void testGetRelativeDNARates() throws Exception { //TODO: Test goes here... } public void testSetGetKappa() throws Exception { //TODO: Test goes here... } public static Test suite() { return new TestSuite(CovarionHKYTest.class); } CovarionHKY model; OldHiddenNucleotides dataType; Parameter kappa; Parameter switchingRate; Parameter alpha; /* %The following Matlab code provided by David Bryant was used to calculate %the probabilities for a covarion K2P model with hidden rates of 0 and 1, %kappa of 2 and a switching rate of 1 format long; clear all; r=4; k=2; Q = [-2-k 1 k 1; 1 -2-k 1 k; k 1 -2-k 1; 1 k 1 -2-k]; piQ = diag([0.25 0.25 0.25 0.25]); %Next step necessary for unequal base cases. Commented out %to match Alexei's code. %Q = Q*piQ; G = [-1 1; 1 -1]; %Transition matrix for rate classes piG = diag([1/2 1/2]); D = diag([0 1]); I = eye(r,r); R = kron(D,Q) + kron(G,I); pi = kron(piG,piQ); %Matrix C is a 0 one matrix. Picks out where bases change C = kron(ones(2,2),ones(r,r)-eye(r)); %Normalise rate = sum(sum(C.*(pi*R))); R = R/rate; for i=1:1:100 t = i*0.01; P = expm(R*t); pchange(i) = sum(sum(C.*(pi*P))); end pchange' */ static final double[] matLabPChange_K2P = { 0.00986400785088, 0.01946196037798, 0.02880252532240, 0.03789407770533, 0.04674470983328, 0.05536224095885, 0.06375422660885, 0.07192796759140, 0.07989051869303, 0.08764869707661, 0.09520909039030, 0.10257806459771, 0.10976177153876, 0.11676615623066, 0.12359696391788, 0.13025974687993, 0.13675987100520, 0.14310252213901, 0.14929271221364, 0.15533528516793, 0.16123492266363, 0.16699614960557, 0.17262333947249, 0.17812071946489, 0.18349237547643, 0.18874225689480, 0.19387418123805, 0.19889183863200, 0.20379879613422, 0.20859850190993, 0.21329428926481, 0.21788938053981, 0.22238689087256, 0.22678983183007, 0.23110111491715, 0.23532355496481, 0.23945987340274, 0.24351270142001, 0.24748458301758, 0.25137797795665, 0.25519526460616, 0.25893874269311, 0.26261063595893, 0.26621309472519, 0.26974819837174, 0.27321795773026, 0.27662431739624, 0.27996915796208, 0.28325429817402, 0.28648149701564, 0.28965245572033, 0.29276881971518, 0.29583218049872, 0.29884407745474, 0.30180599960433, 0.30471938729832, 0.30758563385225, 0.31040608712556, 0.31318205104731, 0.31591478708987, 0.31860551569264, 0.32125541763742, 0.32386563537701, 0.32643727431884, 0.32897140406488, 0.33146905960964, 0.33393124249748, 0.33635892194069, 0.33875303589965, 0.34111449212649, 0.34344416917334, 0.34574291736640, 0.34801155974717, 0.35025089298173, 0.35246168823930, 0.35464469204110, 0.35680062708056, 0.35893019301571, 0.36103406723493, 0.36311290559680, 0.36516734314491, 0.36719799479871, 0.36920545602088, 0.37119030346234, 0.37315309558549, 0.37509437326636, 0.37701466037663, 0.37891446434595, 0.38079427670539, 0.38265457361258, 0.38449581635920, 0.38631845186142, 0.38812291313386, 0.38990961974765, 0.39167897827306, 0.39343138270730, 0.39516721488804, 0.39688684489292, 0.39859063142581, 0.40027892219004 }; static final double[] matLabPChange_HKY = { 0.00985954695768, 0.01944735733283, 0.02877658973544, 0.03785964647441, 0.04670821854965, 0.05533332794691, 0.06374536739747, 0.07195413775516, 0.07996888313383, 0.08779832394009, 0.09545068792763, 0.10293373939203, 0.11025480661789, 0.11742080768318, 0.12443827471949, 0.13131337672096, 0.13805194098895, 0.14465947329452, 0.15114117683550, 0.15750197006065, 0.16374650342878, 0.16987917516689, 0.17590414608721, 0.18182535351970, 0.18764652441305, 0.19337118765401, 0.19900268565191, 0.20454418523250, 0.20999868788237, 0.21536903938304, 0.22065793887109, 0.22586794735891, 0.23100149574820, 0.23606089236674, 0.24104833005692, 0.24596589284278, 0.25081556220095, 0.25559922295905, 0.26031866884393, 0.26497560770062, 0.26957166640174, 0.27410839546591, 0.27858727340240, 0.28300971079857, 0.28737705416541, 0.29169058955549, 0.29595154596716, 0.30016109854757, 0.30432037160662, 0.30843044145304, 0.31249233906332, 0.31650705259341, 0.32047552974251, 0.32439867997787, 0.32827737662883, 0.33211245885791, 0.33590473351622, 0.33965497689022, 0.34336393634610, 0.34703233187806, 0.35066085756610, 0.35425018294872, 0.35780095431559, 0.36131379592506, 0.36478931115072, 0.36822808356150, 0.37163067793909, 0.37499764123646, 0.37832950348097, 0.38162677862531, 0.38488996534947, 0.38811954781650, 0.39131599638496, 0.39447976828045, 0.39761130822876, 0.40071104905293, 0.40377941223623, 0.40681680845316, 0.40982363807038, 0.41280029161931, 0.41574715024200, 0.41866458611203, 0.42155296283168, 0.42441263580698, 0.42724395260180, 0.43004725327227, 0.43282287068266, 0.43557113080388, 0.43829235299547, 0.44098685027222, 0.44365492955622, 0.44629689191518, 0.44891303278793, 0.45150364219766, 0.45406900495387, 0.45660940084344, 0.45912510481158, 0.46161638713329, 0.46408351357576, 0.46652674555230 }; }