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

Example 6 with MarkovJumpsSubstitutionModel

use of dr.evomodel.substmodel.MarkovJumpsSubstitutionModel 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 7 with MarkovJumpsSubstitutionModel

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

the class MarkovJumpsSubstitutionModelTest method testMarkovJumpsCounts.

public void testMarkovJumpsCounts() {
    HKY substModel = new HKY(2.0, new FrequencyModel(Nucleotides.INSTANCE, // A,C,G,T
    new double[] { 0.3, 0.2, 0.25, 0.25 }));
    int states = substModel.getDataType().getStateCount();
    MarkovJumpsSubstitutionModel markovjumps = new MarkovJumpsSubstitutionModel(substModel, MarkovJumpsType.COUNTS);
    double[] r = new double[states * states];
    double[] q = new double[states * states];
    double[] j = new double[states * states];
    double[] c = new double[states * states];
    double[] p = new double[states * states];
    double time = 1.0;
    // A
    int from = 0;
    // C
    int to = 1;
    MarkovJumpsCore.fillRegistrationMatrix(r, from, to, states, 1.0);
    markovjumps.setRegistration(r);
    substModel.getInfinitesimalMatrix(q);
    substModel.getTransitionProbabilities(time, p);
    markovjumps.computeJointStatMarkovJumps(time, j);
    markovjumps.computeCondStatMarkovJumps(time, c);
    MarkovJumpsCore.makeComparableToRPackage(q);
    System.out.println("Q = " + new Vector(q));
    MarkovJumpsCore.makeComparableToRPackage(p);
    System.out.println("P = " + new Vector(p));
    System.out.println("Counts:");
    MarkovJumpsCore.makeComparableToRPackage(r);
    System.out.println("R = " + new Vector(r));
    MarkovJumpsCore.makeComparableToRPackage(j);
    System.out.println("J = " + new Vector(j));
    assertEquals(rMarkovJumpsJ, j, tolerance);
    MarkovJumpsCore.makeComparableToRPackage(c);
    System.err.println("C = " + new Vector(c));
    assertEquals(rMarkovJumpsC, c, tolerance);
}

Example 8 with MarkovJumpsSubstitutionModel

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

the class MarkovJumpsSubstitutionModelTest method testMarkovJumpsReward.

public void testMarkovJumpsReward() {
    HKY substModel = new HKY(2.0, new FrequencyModel(Nucleotides.INSTANCE, // A,C,G,T
    new double[] { 0.3, 0.2, 0.25, 0.25 }));
    int states = substModel.getDataType().getStateCount();
    double[] j = new double[states * states];
    double[] c = new double[states * states];
    double time = 1.0;
    MarkovJumpsSubstitutionModel markovjumps = new MarkovJumpsSubstitutionModel(substModel, MarkovJumpsType.REWARDS);
    double[] rewards = { 1.0, 1.0, 1.0, 1.0 };
    markovjumps.setRegistration(rewards);
    markovjumps.computeJointStatMarkovJumps(time, j);
    markovjumps.computeCondStatMarkovJumps(time, c);
    System.out.println("Rewards:");
    MarkovJumpsCore.makeComparableToRPackage(rewards);
    System.out.println("R = " + new Vector(rewards));
    MarkovJumpsCore.makeComparableToRPackage(j);
    System.out.println("J = " + new Vector(j));
    assertEquals(rMarkovRewardsJ, j, tolerance);
    MarkovJumpsCore.makeComparableToRPackage(c);
    System.out.println("C = " + new Vector(c));
    assertEquals(rMarkovRewardsC, c, tolerance);
}

Example 9 with MarkovJumpsSubstitutionModel

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

the class UniformizedStateHistoryTest method testStateHistorySimulationForRewards.

public void testStateHistorySimulationForRewards() {
    try {
        double startingTime = 1.0;
        double endingTime = 3.0;
        int startingState = 1;
        int endingState = 3;
        int N = 1000000;
        double[] tmp = new double[stateCount * stateCount];
        hky.getTransitionProbabilities(endingTime - startingTime, tmp);
        double transitionProbability = tmp[startingState * stateCount + endingState];
        double[][] registers = new double[3][stateCount];
        // Reward all states
        Arrays.fill(registers[0], 1.0);
        // Reward just one state!
        registers[1][0] = 1.0;
        // Reward just one state!
        registers[2][3] = 1.0;
        double[] expectations = new double[registers.length];
        for (int i = 0; i < N; i++) {
            StateHistory history = UniformizedStateHistory.simulateConditionalOnEndingState(startingTime, startingState, endingTime, endingState, transitionProbability, stateCount, process);
            for (int j = 0; j < registers.length; j++) {
                expectations[j] += history.getTotalReward(registers[j]);
            }
        }
        // Determine analytic solution
        MarkovJumpsSubstitutionModel markovjumps = new MarkovJumpsSubstitutionModel(hky, MarkovJumpsType.REWARDS);
        double[] mjExpectations = new double[stateCount * stateCount];
        for (int j = 0; j < registers.length; j++) {
            expectations[j] /= (double) N;
            System.out.println("Expected reward for register[" + j + "] = " + expectations[j]);
            markovjumps.setRegistration(registers[j]);
            markovjumps.computeCondStatMarkovJumps(endingTime - startingTime, mjExpectations);
            assertEquals(mjExpectations[startingState * stateCount + endingState], expectations[j], 1E-2);
        }
    } catch (SubordinatedProcess.Exception e) {
        throw new RuntimeException("Subordinated process exception");
    }
}