/* * ColourChangeMatrix.java * * Copyright (c) 2002-2015 Alexei Drummond, Andrew Rambaut and Marc Suchard * * This file is part of BEAST. * See the NOTICE file distributed with this work for additional * information regarding copyright ownership and licensing. * * BEAST is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * BEAST is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with BEAST; if not, write to the * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, * Boston, MA 02110-1301 USA */ package dr.evolution.colouring; /** * MigrationMatrix.java * * Represents a migration matrix, both backward and forward in time. * * @author Gerton Lunter * * @version 1 */ public class ColourChangeMatrix { public ColourChangeMatrix() { } public ColourChangeMatrix( double[] changeRates, int numColours ) { // For now, assume two colours if (numColours != 2) { throw new IllegalArgumentException("Only 2 colours supported"); } this.numColours = numColours; this.bwMatrix = new double[numColours][numColours]; this.equilibrium = new double[numColours]; // Populate matrix bwMatrix[0][1] = changeRates[0]; bwMatrix[1][0] = changeRates[1]; calculateExitRates(); // Calculate equilibrium distribution calculateEquilibrium(); } /* * populate diagonal of bw migration matrix */ private void calculateExitRates() { for (int i=0; i<numColours; i++) { double exitRate = 0.0; for (int j=0; j<numColours; j++) { if (i!=j) { exitRate += bwMatrix[i][j]; } } bwMatrix[i][i] = -exitRate; } } private void calculateEquilibrium() { if (numColours == 2) { double r01 = bwMatrix[0][1]; double r10 = bwMatrix[1][0]; equilibrium[0] = r10/(r01+r10); equilibrium[1] = r01/(r01+r10); } else { throw new Error("Only 2 colours supported"); } } /* * Returns backward migration rate from i (child) to j (parent) * */ public double getBackwardRate(int i, int j) { return bwMatrix[i][j]; } /* * Returns forward migration rate from (parent) i to (child) j */ public double getForwardRate(int i, int j) { return equilibrium[j]*bwMatrix[j][i]/equilibrium[i]; } /** * The probability of ending in state y after time t, conditional on starting in state x, * according to the forward evolution matrix. * * @param x starting state * @param y end state * @param t elapsed time * @return probability of starting from state x and ending in state y after time t. */ public double forwardTimeEvolution(int x, int y, double t) { if (t<0) { throw new IllegalArgumentException("Cannot go backwards in time: t="+t); } double m[] = {getForwardRate(0,1), getForwardRate(1,0)}; double mt = m[0]+m[1]; if (y==x) { return (m[x]*Math.exp(-mt*t) + m[1-x])/mt; } else { return m[x]*(1.0-Math.exp(-mt*t))/mt; } } /** * The probability of ending in state y after time t, conditional on starting in state x, * according to the backward evolution matrix. * * @param x starting state (child) * @param y end state (parent) * @param t elapsed time (backwards in time) * @return probability of starting from state x and ending in state y after time t. */ public double backwardTimeEvolution(int x, int y, double t) { if (t<0) { throw new IllegalArgumentException("Cannot go backwards in time: t="+t); } double m[] = {getBackwardRate(0,1), getBackwardRate(1,0)}; double mt = m[0]+m[1]; if (y==x) { // return (m[1-y]*Math.exp(-mt*t) + m[y])/mt; return (m[x]*Math.exp(-mt*t) + m[1-x])/mt; } else { // return m[y]*(1.0-Math.exp(-mt*t))/mt; return m[x]*(1.0-Math.exp(-mt*t))/mt; } } /** * * @return equilibrium distribution */ public double[] getEquilibrium() { return equilibrium.clone(); } public double getEquilibrium(int i) { return equilibrium[i]; } public static void main(String[] args) { double[] pars = {1.0, 0.5}; ColourChangeMatrix mm = new ColourChangeMatrix( pars, 2 ); System.out.println( "BW 0->0, matrix:"+mm.getBackwardRate(0,0)); System.out.println( "BW 0->1, matrix:"+mm.getBackwardRate(0,1)); System.out.println( "BW 1->0, matrix:"+mm.getBackwardRate(1,0)); System.out.println( "BW 1->1, matrix:"+mm.getBackwardRate(1,1)); System.out.println( "FW 0->0, matrix:"+mm.getForwardRate(0,0)); System.out.println( "FW 0->1, matrix:"+mm.getForwardRate(0,1)); System.out.println( "FW 1->0, matrix:"+mm.getForwardRate(1,0)); System.out.println( "FW 1->1, matrix:"+mm.getForwardRate(1,1)); System.out.println( "equilibrium="+mm.equilibrium[0]+","+mm.equilibrium[1] ); System.out.println( "BW 0->0, t=1:"+mm.backwardTimeEvolution(0,0,1.0)); System.out.println( "BW 0->1, t=1:"+mm.backwardTimeEvolution(0,1,1.0)); System.out.println( "BW 1->0, t=1:"+mm.backwardTimeEvolution(1,0,1.0)); System.out.println( "BW 1->1, t=1:"+mm.backwardTimeEvolution(1,1,1.0)); System.out.println( "FW 0->0, t=1:"+mm.forwardTimeEvolution(0,0,1.0)); System.out.println( "FW 0->1, t=1:"+mm.forwardTimeEvolution(0,1,1.0)); System.out.println( "FW 1->0, t=1:"+mm.forwardTimeEvolution(1,0,1.0)); System.out.println( "FW 1->1, t=1:"+mm.forwardTimeEvolution(1,1,1.0)); System.out.println( "BW 0->0, t=infty:"+mm.backwardTimeEvolution(0,0,1000.0)); System.out.println( "BW 0->1, t=infty:"+mm.backwardTimeEvolution(0,1,1000.0)); System.out.println( "BW 1->0, t=infty:"+mm.backwardTimeEvolution(1,0,1000.0)); System.out.println( "BW 1->1, t=infty:"+mm.backwardTimeEvolution(1,1,1000.0)); System.out.println( "FW 0->0, t=infty:"+mm.forwardTimeEvolution(0,0,1000.0)); System.out.println( "FW 0->1, t=infty:"+mm.forwardTimeEvolution(0,1,1000.0)); System.out.println( "FW 1->0, t=infty:"+mm.forwardTimeEvolution(1,0,1000.0)); System.out.println( "FW 1->1, t=infty:"+mm.forwardTimeEvolution(1,1,1000.0)); } // private stuff private int numColours; private double[][] bwMatrix; private double[] equilibrium; }