Examples with HKY dr.evomodel.substmodel.nucleotide.HKY used on opensource projects

Example 6 with HKY

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

the class MarkovJumpsSubstitutionModelTest method testMarginalRates.

public void testMarginalRates() {
    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];
    // A
    int from = 0;
    // C
    int to = 1;
    MarkovJumpsCore.fillRegistrationMatrix(r, from, to, states, 1.0);
    markovjumps.setRegistration(r);
    double marginalRate = markovjumps.getMarginalRate();
    System.out.println("Marginal rate = " + marginalRate);
    assertEquals(rMarkovMarginalRate, marginalRate, tolerance);
    MarkovJumpsCore.fillRegistrationMatrix(r, states);
    markovjumps.setRegistration(r);
    marginalRate = markovjumps.getMarginalRate();
    System.out.println("Marginal rate = " + marginalRate);
    assertEquals(1.0, marginalRate, tolerance);
}

Example 7 with HKY

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

the class CompleteHistorySimulatorTest method testHKYSimulation.

public void testHKYSimulation() {
    Parameter kappa = new Parameter.Default(1, 2.0);
    double[] pi = { 0.45, 0.05, 0.25, 0.25 };
    Parameter freqs = new Parameter.Default(pi);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    int stateCount = hky.getDataType().getStateCount();
    Parameter mu = new Parameter.Default(1, 0.5);
    Parameter alpha = new Parameter.Default(1, 0.5);
    GammaSiteRateModel siteModel = new GammaSiteRateModel("gammaModel", mu, alpha, 4, null);
    siteModel.setSubstitutionModel(hky);
    BranchRateModel branchRateModel = new DefaultBranchRateModel();
    double analyticResult = TreeUtils.getTreeLength(tree, tree.getRoot()) * mu.getParameterValue(0);
    int nSites = 200;
    double[] register1 = new double[stateCount * stateCount];
    double[] register2 = new double[stateCount * stateCount];
    // Count all jumps
    MarkovJumpsCore.fillRegistrationMatrix(register1, stateCount);
    // Move some jumps from 1 to 2
    register1[1 * stateCount + 2] = 0;
    register2[1 * stateCount + 2] = 1;
    register1[1 * stateCount + 3] = 0;
    register2[1 * stateCount + 3] = 1;
    register1[2 * stateCount + 3] = 0;
    register2[2 * stateCount + 3] = 1;
    runSimulation(tree, siteModel, branchRateModel, nSites, new double[][] { register1, register2 }, analyticResult);
}

Example 8 with HKY

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

the class HiddenMarkovRatesTest method testHiddenRates.

public void testHiddenRates() {
    final int index = 0;
    double[] freqs = getBinaryFreqs(index);
    FrequencyModel binaryFreqModel = new FrequencyModel(TwoStates.INSTANCE, freqs);
    int relativeTo = 0;
    Parameter ratesParameter = new Parameter.Default(0);
    GeneralSubstitutionModel binaryModel = new GeneralSubstitutionModel("binary", TwoStates.INSTANCE, binaryFreqModel, ratesParameter, relativeTo);
    UniformizedSubstitutionModel uSM = new UniformizedSubstitutionModel(binaryModel, MarkovJumpsType.REWARDS);
    uSM.setSaveCompleteHistory(true);
    double[] rewardRegister = new double[] { 0.0, 1.0 };
    uSM.setRegistration(rewardRegister);
    final double[] hkyFreqs = getHKYFreqs(index);
    FrequencyModel hkyFreqModel = new FrequencyModel(Nucleotides.INSTANCE, hkyFreqs);
    final double kappa = getKappa(index);
    final HKY hky = new HKY(kappa, hkyFreqModel);
    final double length = getLength(index);
    double[] resultCompleteHistory = new double[16];
    final int replicates = getNumberReplicates(index);
    double result = 0.0;
    for (int r = 0; r < replicates; ++r) {
        result += oneCompleteHistoryReplicate(resultCompleteHistory, hky, uSM, length);
    }
    result /= replicates;
    normalize(resultCompleteHistory, replicates);
    System.out.println("Averaged probabilities");
    System.out.println(result);
    System.out.println(new Vector(resultCompleteHistory));
    System.out.println();
    double[] intermediate = new double[16];
    hky.getTransitionProbabilities(result, intermediate);
    System.out.println("Intermediate using above average reward");
    System.out.println(result);
    System.out.println(new Vector(intermediate));
    System.out.println();
    double[] resultExpected = new double[16];
    UniformizedSubstitutionModel expectedUSM = new UniformizedSubstitutionModel(binaryModel, MarkovJumpsType.REWARDS, replicates);
    expectedUSM.setRegistration(rewardRegister);
    result = oneCompleteHistoryReplicate(resultExpected, hky, expectedUSM, length);
    System.out.println("Averaged reward");
    System.out.println(result);
    System.out.println(new Vector(resultExpected));
    System.out.println();
    double[] originalProbs = new double[16];
    hky.getTransitionProbabilities(length, originalProbs);
    System.out.println("Original probabilities");
    System.out.println(new Vector(originalProbs));
    System.out.println();
}

Example 9 with HKY

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

the class BeagleSeqSimTest method simulateThreePartitions.

// END: simulateTwoPartitions
static void simulateThreePartitions(int i, int N) {
    try {
        MathUtils.setSeed(666);
        System.out.println("Test case 3: simulateThreePartitions");
        int sequenceLength = 100000;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        Tree tree = importer.importTree(null);
        TreeModel treeModel = new TreeModel(tree);
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
        FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqs);
        // create substitution model
        Parameter kappa = new Parameter.Default(1, 10);
        HKY hky = new HKY(kappa, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(hky);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        sequenceLength - 1, // every
        3);
        // create partition
        Partition Partition = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        1, // to
        sequenceLength - 1, // every
        3);
        // create partition
        Partition partition3 = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        2, // to
        sequenceLength - 1, // every
        3);
        partitionsList.add(partition1);
        partitionsList.add(Partition);
        partitionsList.add(partition3);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        if (i == (N - 1)) {
            System.out.println(simulator.simulate(simulateInPar, false).toString());
        } else {
            simulator.simulate(simulateInPar, false);
        }
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 10 with HKY

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

the class BeagleSeqSimTest method simulateTopology.

// END: annotateTree
static void simulateTopology() {
    try {
        System.out.println("Test case 1: simulateTopology");
        MathUtils.setSeed(666);
        int sequenceLength = 10;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        Tree tree = importer.importTree(null);
        // set demographic function
        ExponentialGrowth exponentialGrowth = new ExponentialGrowth(Units.Type.YEARS);
        exponentialGrowth.setN0(10);
        exponentialGrowth.setGrowthRate(0.5);
        Taxa taxa = new Taxa();
        for (Taxon taxon : tree.asList()) {
            double absoluteHeight = Utils.getAbsoluteTaxonHeight(taxon, tree);
            taxon.setAttribute(Utils.ABSOLUTE_HEIGHT, absoluteHeight);
            // taxon.setAttribute("date", new Date(absoluteHeight,
            // Units.Type.YEARS, true));
            taxa.addTaxon(taxon);
        }
        // END: taxon loop
        CoalescentSimulator topologySimulator = new CoalescentSimulator();
        TreeModel treeModel = new TreeModel(topologySimulator.simulateTree(taxa, exponentialGrowth));
        System.out.println(treeModel.toString());
        Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
        FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqs);
        // create substitution model
        Parameter kappa = new Parameter.Default(1, 10);
        HKY hky = new HKY(kappa, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(hky);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        sequenceLength - 1, // every
        1);
        partitionsList.add(partition1);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        System.out.println(simulator.simulate(simulateInPar, false).toString());
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}