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

Example 1 with ComplexSubstitutionModel

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

the class TwoStateMarkovRewardsTest method testTwoStateRewards.

public void testTwoStateRewards() {
    DataType dataType = TwoStates.INSTANCE;
    FrequencyModel freqModel = new FrequencyModel(TwoStates.INSTANCE, new double[] { 0.5, 0.5 });
    Parameter rates = new Parameter.Default(new double[] { 4.0, 6.0 });
    ComplexSubstitutionModel twoStateModel = new ComplexSubstitutionModel("two", dataType, freqModel, rates) {
    };
    twoStateModel.setNormalization(false);
    MarkovJumpsSubstitutionModel markovRewards = new MarkovJumpsSubstitutionModel(twoStateModel, MarkovJumpsType.REWARDS);
    double[] r = new double[2];
    double[] q = new double[4];
    double[] c = new double[4];
    int mark = 0;
    double weight = 1.0;
    r[mark] = weight;
    markovRewards.setRegistration(r);
    twoStateModel.getInfinitesimalMatrix(q);
    System.out.println("Q = " + new Vector(q));
    System.out.println("Reward for state 0");
    double time = 1.0;
    markovRewards.computeCondStatMarkovJumps(time, c);
    System.out.println("Reward conditional on X(0) = i, X(t) = j: " + new Vector(c));
    double endTime = 10.0;
    int steps = 10;
    for (time = 0.0; time < endTime; time += (endTime / steps)) {
        markovRewards.computeCondStatMarkovJumps(time, c);
        // start = 0, end = 0
        System.out.println(time + "," + c[0]);
    }
}

Example 2 with ComplexSubstitutionModel

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

the class EigendecompositionTest method doTest.

private void doTest(Parameter.Default rateMatrix) {
    normalise(rateMatrix);
    ComplexSubstitutionModel substModel = new ComplexSubstitutionModel("test", dataType, freqModel, rateMatrix);
    boolean eigenDecompOK = testEigenDecomposition(substModel);
    double[] probs = new double[dim * dim];
    substModel.getTransitionProbabilities(1, probs);
    boolean connCheckResult = true;
    String eSuccess = eigenDecompOK ? "Succeeded" : "Failed (all zero eigenvectors returned)";
    String cSuccess = connCheckResult ? "Passed" : "Failed";
    System.out.println("Eigendecomposition: " + eSuccess);
    System.out.println("Connectivity check: " + cSuccess);
    System.out.println("Transition probability matrix:");
    for (int row = 0; row < dim; row++) {
        for (int col = 0; col < dim; col++) {
            System.out.print(probs[row * dim + col] + " ");
        }
        System.out.println();
    }
}

Example 3 with ComplexSubstitutionModel

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

the class EigendecompositionTest method doRandomTests.

private void doRandomTests() {
    ComplexSubstitutionModel substModel = new ComplexSubstitutionModel("test", dataType, freqModel, rateMatrix);
    ArrayList<Integer[]> indicators = makeMatrices(getAllPossibilities(dim), columnNumberLookup);
    for (Integer[] ind : indicators) {
        int successED = 0;
        int[] successTP = new int[dim];
        int[] zeroRows = new int[dim];
        for (int row = 0; row < dim; row++) {
            boolean zeroRow = true;
            for (int col = 0; col < dim; col++) {
                if (row != col && ind[columnNumberLookup[row][col]] != 0) {
                    zeroRow = false;
                }
            }
            zeroRows[row] = zeroRow ? 1 : 0;
        }
        for (int run = 0; run < 100; run++) {
            for (int i = 0; i < dim * (dim - 1); i++) {
                rateMatrix.setParameterValue(i, MathUtils.nextDouble() * ind[i]);
            }
            normalise(rateMatrix);
            boolean EigenDecompOK = testEigenDecomposition(substModel);
            if (!EigenDecompOK) {
                for (int row = 0; row < dim; row++) {
                    for (int col = 0; col < dim; col++) {
                        if (row == col) {
                            System.out.print("- ");
                        } else {
                            System.out.print(rateMatrix.getParameterValue(columnNumberLookup[row][col]) + " ");
                        }
                    }
                    System.out.println();
                }
                System.out.println();
            }
            Boolean[] transProbsOK = new Boolean[dim];
            for (int i = 0; i < dim; i++) {
                transProbsOK[i] = zeroRows[i] != 1 || testProbMatrix(substModel, i, dim, tolerance);
            }
            if (EigenDecompOK) {
                successED++;
            }
            for (int i = 0; i < dim; i++) if (EigenDecompOK && transProbsOK[i]) {
                successTP[i]++;
            }
        }
        if (successED != 100) {
            System.out.print("Indicators: ");
            for (int i = 0; i < dim * (dim - 1); i++) {
                System.out.print(ind[i] + " ");
            }
            System.out.println();
            System.out.println("Eigendecomposition succeeded " + successED + " times out of 100");
            for (int i = 0; i < dim; i++) {
                if (zeroRows[i] == 1) {
                    System.out.println("Row " + (i + 1) + ": Connectivity check would work " + successTP[i] + " times out of " + successED);
                }
            }
            System.out.println();
        }
    }
}