    001    package calhoun.analysis.crf.features.supporting.phylogenetic;
    002
    003    import java.io.Serializable;
    004
    005    import org.apache.commons.logging.Log;
    006    import org.apache.commons.logging.LogFactory;
    007
    008    import calhoun.util.Assert;
    009    import calhoun.util.ColtUtil;
    010    import cern.colt.matrix.DoubleMatrix2D;
    011    import cern.colt.matrix.impl.DenseDoubleMatrix2D;
    012
    013    public abstract class NucleotideSubstitutionModel implements Serializable{
    014            private static final Log log = LogFactory.getLog(NucleotideSubstitutionModel.class);
    015
    016            /* References:
    017             *
    018             * http://umber.sbs.man.ac.uk/resources/phase/manual/node65.html
    019             * (a description of nucleotide substitution models implemented by the program PHASE by Gowri-Shankar Vivek
    020             *
    021             * Whelan, S., P. Lio, and N. Goldman
    022             * 2001. Molecular phylogenetics: state-of-the art methods for looking into the past.
    023             * TRENDS in Genetics, 17(5):262-272.
    024             */
    025
    026            DoubleMatrix2D R = new DenseDoubleMatrix2D(4,4); // nucleotide substitution rate matrix
    027
    028            public abstract String getEvolutionaryModelName();
    029
    030            public abstract double[] getParameters();
    031
    032            public abstract void setParameters(double[] parms); // responsible for setting the transition rate matrix.
    033
    034            public DoubleMatrix2D getRateMatrix() {
    035                    return R;
    036            }
    037
    038            public  DoubleMatrix2D transitionMatrix(double t) {
    039                    //return exp(t*R) is what we want to do, but isit really that simple?
    040
    041                    DoubleMatrix2D S = new DenseDoubleMatrix2D(4,4);
    042                    S.assign(R);
    043                    ColtUtil.scalarmultiply(S,t);
    044
    045                    DoubleMatrix2D T = new DenseDoubleMatrix2D(4,4);
    046                    //ColtUtil.exponentiate_symmetric_matrix(S,T);
    047                    ColtUtil.exponentiate_real_matrix(S,T,100);
    048
    049                    // Now let us verify that the row sums of T are all zero.
    050                    for (int i=0; i<4; i++) {
    051                            double sumi = 0;
    052                            for (int j=0; j<4; j++) {
    053                                    sumi += T.getQuick(i,j);
    054                            }
    055                            if ((sumi<0.999) || (sumi>1.001) ) {
    056                                    log.warn("Assertion failure imminent.  The matrix S is:");
    057                                    log.warn(ColtUtil.format(S));
    058                                    log.warn("The matrix T (which should be exponential of S) is:");
    059                                    log.warn(ColtUtil.format(T));
    060                                    Assert.a(false);
    061                            }
    062                    }
    063
    064                    return T;
    065            }
    066
    067            public abstract void summarize();
    068
    069    }