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

Example 6 with EigenDecomposition

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

the class TN93 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 kappa1 = getKappa1();
        final double kappa2 = getKappa2();
        final double beta = -1.0 / (2.0 * (piR * piY + kappa1 * pi[0] * pi[2] + kappa2 * pi[1] * pi[3]));
        final double A_R = 1.0 + piR * (kappa1 - 1);
        final double A_Y = 1.0 + piY * (kappa2 - 1);
        eval[1] = beta;
        eval[2] = beta * A_Y;
        eval[3] = beta * A_R;
        updateMatrix = false;
    }
    return eigenDecomposition;
}

Example 7 with EigenDecomposition

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

the class SubstitutionModelDelegate method updateSubstitutionModels.

// END: getStateFrequencies
@Override
public void updateSubstitutionModels(Beagle beagle, boolean flipBuffers) {
    for (int i = 0; i < eigenCount; i++) {
        if (flipBuffers) {
            eigenBufferHelper.flipOffset(i);
        }
        EigenDecomposition ed = substitutionModelList.get(i).getEigenDecomposition();
        beagle.setEigenDecomposition(eigenBufferHelper.getOffsetIndex(i), ed.getEigenVectors(), ed.getInverseEigenVectors(), ed.getEigenValues());
    }
}

Example 8 with EigenDecomposition

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

the class HomogenousSubstitutionModelDelegate method updateSubstitutionModels.

@Override
public void updateSubstitutionModels(Beagle beagle, boolean flip) {
    if (flip) {
        eigenBufferHelper.flipOffset(0);
    }
    EigenDecomposition ed = substitutionModel.getEigenDecomposition();
    beagle.setEigenDecomposition(eigenBufferHelper.getOffsetIndex(0), ed.getEigenVectors(), ed.getInverseEigenVectors(), ed.getEigenValues());
}