Examples with EigenDecomposition dr.evomodel.substmodel.EigenDecomposition used on opensource projects

Example 1 with EigenDecomposition

use of dr.evomodel.substmodel.EigenDecomposition in project beast-mcmc by beast-dev.

the class SubstitutionModelDelegate method updateSubstitutionModels.

public void updateSubstitutionModels(Beagle beagle) {
    for (int i = 0; i < eigenCount; i++) {
        eigenBufferHelper.flipOffset(i);
        EigenDecomposition ed = substitutionModelList.get(i).getEigenDecomposition();
        beagle.setEigenDecomposition(eigenBufferHelper.getOffsetIndex(i), ed.getEigenVectors(), ed.getInverseEigenVectors(), ed.getEigenValues());
    }
}

Example 2 with EigenDecomposition

use of dr.evomodel.substmodel.EigenDecomposition in project beast-mcmc by beast-dev.

the class TN93Test method testTN93.

public void testTN93() {
    Parameter kappa1 = new Parameter.Default(5.0);
    Parameter kappa2 = new Parameter.Default(2.0);
    double[] pi = new double[] { 0.40, 0.20, 0.30, 0.10 };
    double time = 0.1;
    FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, pi);
    TN93 tn = new TN93(kappa1, kappa2, freqModel);
    EigenDecomposition decomp = tn.getEigenDecomposition();
    Vector eval = new Vector(decomp.getEigenValues());
    System.out.println("Eval = " + eval);
    double[] probs = new double[16];
    tn.getTransitionProbabilities(time, probs);
    System.out.println("new probs = " + new Vector(probs));
    // check against old implementation
    dr.oldevomodel.substmodel.FrequencyModel oldFreq = new dr.oldevomodel.substmodel.FrequencyModel(Nucleotides.INSTANCE, pi);
    dr.oldevomodel.substmodel.TN93 oldTN = new dr.oldevomodel.substmodel.TN93(kappa1, kappa2, oldFreq);
    double[] oldProbs = new double[16];
    oldTN.getTransitionProbabilities(time, oldProbs);
    System.out.println("old probs = " + new Vector(oldProbs));
    assertEquals(probs, oldProbs, 10E-6);
}

Example 3 with EigenDecomposition

use of dr.evomodel.substmodel.EigenDecomposition in project beast-mcmc by beast-dev.

the class MarkovModulatedGY94CodonModel method main.

public static void main(String[] args) {
    GY94CodonModel codonModel = new GY94CodonModel(Codons.UNIVERSAL, new Parameter.Default(1.0), new Parameter.Default(2.0), new FrequencyModel(Codons.UNIVERSAL, new Parameter.Default(61, 1.0 / 61.0)));
    EigenDecomposition ed1 = codonModel.getEigenDecomposition();
    //        double[][] q = codonModel.getQ();
    //        System.err.println("matrixQ = \n"+codonModel.printQ());// new Matrix(q));
    FrequencyModel freqModel = new FrequencyModel(HiddenCodons.UNIVERSAL_HIDDEN_2, new Parameter.Default(122, 1.0 / 122.0));
    System.err.println("freq = " + new Vector(freqModel.getFrequencies()));
    //        System.exit(-1);
    MarkovModulatedGY94CodonModel mmCodonModel = new MarkovModulatedGY94CodonModel(HiddenCodons.UNIVERSAL_HIDDEN_2, new Parameter.Default(2, 5.0), new Parameter.Default(2, 1.0), new Parameter.Default(2.0), freqModel);
    EigenDecomposition ed2 = mmCodonModel.getEigenDecomposition();
    // new Matrix(q));
    System.err.println("matrixQ = \n" + mmCodonModel.printQ());
}

Example 4 with EigenDecomposition

use of dr.evomodel.substmodel.EigenDecomposition in project beast-mcmc by beast-dev.

the class HKY method main.

public static void main(String[] args) {
    //        double kappa = 2.0;
    //        double[] pi = new double[]{0.15,0.30,0.20,0.35};
    //        double time = 0.1;
    double kappa = 1.0;
    double[] pi = new double[] { 0.25, 0.25, 0.25, 0.25 };
    double time = 0.1;
    FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, pi);
    HKY hky = new HKY(kappa, freqModel);
    EigenDecomposition decomp = hky.getEigenDecomposition();
    //        Matrix evec = new Matrix(decomp.getEigenVectors());
    //        Matrix ivec = new Matrix(decomp.getInverseEigenVectors());
    //        System.out.println("Evec =\n"+evec);
    //        System.out.println("Ivec =\n"+ivec);
    Vector eval = new Vector(decomp.getEigenValues());
    System.out.println("Eval = " + eval);
    double[] probs = new double[16];
    hky.getTransitionProbabilities(time, probs);
    System.out.println("new probs = " + new Vector(probs));
    // check against old implementation
    dr.oldevomodel.substmodel.FrequencyModel oldFreq = new dr.oldevomodel.substmodel.FrequencyModel(Nucleotides.INSTANCE, pi);
    dr.oldevomodel.substmodel.HKY oldHKY = new dr.oldevomodel.substmodel.HKY(kappa, oldFreq);
    oldHKY.setKappa(kappa);
    oldHKY.getTransitionProbabilities(time, probs);
    System.out.println("old probs = " + new Vector(probs));
}

Example 5 with EigenDecomposition

use of dr.evomodel.substmodel.EigenDecomposition in project beast-mcmc by beast-dev.

the class HKY method getEigenDecomposition.

public synchronized EigenDecomposition getEigenDecomposition() {
    if (eigenDecomposition == null) {
        double[] evec = new double[stateCount * stateCount];
        double[] ivec = new double[stateCount * stateCount];
        double[] eval = new double[stateCount];
        eigenDecomposition = new EigenDecomposition(evec, ivec, eval);
        // left eigenvectors 3 = (0,1,0,-1); 4 = (1,0,-1,0)
        ivec[2 * stateCount + 1] = 1;
        ivec[2 * stateCount + 3] = -1;
        ivec[3 * stateCount + 0] = 1;
        ivec[3 * stateCount + 2] = -1;
        // right eigenvector 1 = (1,1,1,1)'
        evec[0 * stateCount + 0] = 1;
        evec[1 * stateCount + 0] = 1;
        evec[2 * stateCount + 0] = 1;
        evec[3 * stateCount + 0] = 1;
    }
    if (updateMatrix) {
        double[] evec = eigenDecomposition.getEigenVectors();
        double[] ivec = eigenDecomposition.getInverseEigenVectors();
        double[] pi = freqModel.getFrequencies();
        double piR = pi[0] + pi[2];
        double piY = pi[1] + pi[3];
        // left eigenvector #1
        // or, evec[0] = pi;
        ivec[0 * stateCount + 0] = pi[0];
        ivec[0 * stateCount + 1] = pi[1];
        ivec[0 * stateCount + 2] = pi[2];
        ivec[0 * stateCount + 3] = pi[3];
        // left eigenvector #2
        ivec[1 * stateCount + 0] = pi[0] * piY;
        ivec[1 * stateCount + 1] = -pi[1] * piR;
        ivec[1 * stateCount + 2] = pi[2] * piY;
        ivec[1 * stateCount + 3] = -pi[3] * piR;
        // right eigenvector #2
        evec[0 * stateCount + 1] = 1.0 / piR;
        evec[1 * stateCount + 1] = -1.0 / piY;
        evec[2 * stateCount + 1] = 1.0 / piR;
        evec[3 * stateCount + 1] = -1.0 / piY;
        // right eigenvector #3
        evec[1 * stateCount + 2] = pi[3] / piY;
        evec[3 * stateCount + 2] = -pi[1] / piY;
        // right eigenvector #4
        evec[0 * stateCount + 3] = pi[2] / piR;
        evec[2 * stateCount + 3] = -pi[0] / piR;
        // eigenvectors
        double[] eval = eigenDecomposition.getEigenValues();
        final double kappa = getKappa();
        final double beta = -1.0 / (2.0 * (piR * piY + kappa * (pi[0] * pi[2] + pi[1] * pi[3])));
        final double A_R = 1.0 + piR * (kappa - 1);
        final double A_Y = 1.0 + piY * (kappa - 1);
        eval[1] = beta;
        eval[2] = beta * A_Y;
        eval[3] = beta * A_R;
        updateMatrix = false;
    }
    return eigenDecomposition;
}