Examples with BranchRateModel dr.evomodel.branchratemodel.BranchRateModel used on opensource projects

Example 21 with BranchRateModel

use of dr.evomodel.branchratemodel.BranchRateModel in project beast-mcmc by beast-dev.

the class MsatSamplingTreeLikelihoodTest method setUp.

public void setUp() throws Exception {
    super.setUp();
    //taxa
    ArrayList<Taxon> taxonList3 = new ArrayList<Taxon>();
    Collections.addAll(taxonList3, new Taxon("Taxon1"), new Taxon("Taxon2"), new Taxon("Taxon3"), new Taxon("Taxon4"), new Taxon("Taxon5"), new Taxon("Taxon6"), new Taxon("Taxon7"));
    Taxa taxa3 = new Taxa(taxonList3);
    //msat datatype
    Microsatellite msat = new Microsatellite(1, 6);
    Patterns msatPatterns = new Patterns(msat, taxa3);
    //pattern in the correct code form.
    msatPatterns.addPattern(new int[] { 0, 1, 3, 2, 4, 5, 1 });
    //create tree
    NewickImporter importer = new NewickImporter("(((Taxon1:0.3,Taxon2:0.3):0.6,Taxon3:0.9):0.9,((Taxon4:0.5,Taxon5:0.5):0.3,(Taxon6:0.7,Taxon7:0.7):0.1):1.0);");
    Tree tree = importer.importTree(null);
    //treeModel
    TreeModel treeModel = new TreeModel(tree);
    //msatsubstModel
    AsymmetricQuadraticModel eu1 = new AsymmetricQuadraticModel(msat, null);
    //create msatSamplerTreeModel
    Parameter internalVal = new Parameter.Default(new double[] { 2, 3, 4, 2, 1, 5 });
    int[] externalValues = msatPatterns.getPattern(0);
    HashMap<String, Integer> taxaMap = new HashMap<String, Integer>(externalValues.length);
    boolean internalValuesProvided = true;
    for (int i = 0; i < externalValues.length; i++) {
        taxaMap.put(msatPatterns.getTaxonId(i), i);
    }
    MicrosatelliteSamplerTreeModel msatTreeModel = new MicrosatelliteSamplerTreeModel("JUnitTestEx", treeModel, internalVal, msatPatterns, externalValues, taxaMap, internalValuesProvided);
    //create msatSamplerTreeLikelihood
    BranchRateModel branchRateModel = new StrictClockBranchRates(new Parameter.Default(1.0));
    eu1Likelihood = new MicrosatelliteSamplerTreeLikelihood(msatTreeModel, eu1, branchRateModel);
    //eu2
    TwoPhaseModel eu2 = new TwoPhaseModel(msat, null, eu1, new Parameter.Default(0.0), new Parameter.Default(0.4), null, false);
    eu2Likelihood = new MicrosatelliteSamplerTreeLikelihood(msatTreeModel, eu2, branchRateModel);
    //ec1
    LinearBiasModel ec1 = new LinearBiasModel(msat, null, eu1, new Parameter.Default(0.48), new Parameter.Default(0.0), false, false, false);
    ec1Likelihood = new MicrosatelliteSamplerTreeLikelihood(msatTreeModel, ec1, branchRateModel);
    //ec2
    TwoPhaseModel ec2 = new TwoPhaseModel(msat, null, ec1, new Parameter.Default(0.0), new Parameter.Default(0.4), null, false);
    ec2Likelihood = new MicrosatelliteSamplerTreeLikelihood(msatTreeModel, ec2, branchRateModel);
    //el1
    LinearBiasModel el1 = new LinearBiasModel(msat, null, eu1, new Parameter.Default(0.2), new Parameter.Default(-0.018), true, false, false);
    el1Likelihood = new MicrosatelliteSamplerTreeLikelihood(msatTreeModel, el1, branchRateModel);
    AsymmetricQuadraticModel pu1 = new AsymmetricQuadraticModel(msat, null, new Parameter.Default(1.0), new Parameter.Default(0.015), new Parameter.Default(0.0), new Parameter.Default(1.0), new Parameter.Default(0.015), new Parameter.Default(0.0), false);
    pu1Likelihood = new MicrosatelliteSamplerTreeLikelihood(msatTreeModel, pu1, branchRateModel);
    //ec2
    TwoPhaseModel pu2 = new TwoPhaseModel(msat, null, pu1, new Parameter.Default(0.0), new Parameter.Default(0.4), null, false);
    pu2Likelihood = new MicrosatelliteSamplerTreeLikelihood(msatTreeModel, pu2, branchRateModel);
    //ec1
    LinearBiasModel pc1 = new LinearBiasModel(msat, null, pu1, new Parameter.Default(0.48), new Parameter.Default(0.0), false, false, false);
    pc1Likelihood = new MicrosatelliteSamplerTreeLikelihood(msatTreeModel, pc1, branchRateModel);
}

Example 22 with BranchRateModel

use of dr.evomodel.branchratemodel.BranchRateModel 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 23 with BranchRateModel

use of dr.evomodel.branchratemodel.BranchRateModel in project beast-mcmc by beast-dev.

the class MicrosatelliteSamplerTreeLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    MicrosatelliteSamplerTreeModel mstModel = (MicrosatelliteSamplerTreeModel) xo.getChild(MicrosatelliteSamplerTreeModel.class);
    MicrosatelliteModel microsatelliteModel = (MicrosatelliteModel) xo.getChild(MicrosatelliteModel.class);
    BranchRateModel branchRateModel;
    Object cxo = xo.getChild(BranchRateModel.class);
    if (xo.getChild(BranchRateModel.class) != null) {
        branchRateModel = (BranchRateModel) cxo;
        System.out.println("BranchRateModel provided to MicrosatelliteSamplerTreeLikelihood");
    } else if (xo.hasChildNamed(MUTATION_RATE)) {
        Parameter muRate = (Parameter) xo.getElementFirstChild(MUTATION_RATE);
        branchRateModel = new StrictClockBranchRates(muRate);
        System.out.println("mutation rate provided to MicrosatelliteSamplerTreeLikelihood");
    } else {
        branchRateModel = new StrictClockBranchRates(new Parameter.Default(1.0));
    }
    return new MicrosatelliteSamplerTreeLikelihood(mstModel, microsatelliteModel, branchRateModel);
}

Example 24 with BranchRateModel

use of dr.evomodel.branchratemodel.BranchRateModel in project beast-mcmc by beast-dev.

the class PartitionParser method parseXMLObject.

@Override
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    int from = 0;
    int to = -1;
    int every = xo.getAttribute(EVERY, 1);
    if (xo.hasAttribute(FROM)) {
        from = xo.getIntegerAttribute(FROM) - 1;
        if (from < 0) {
            throw new XMLParseException("Illegal 'from' attribute in patterns element");
        }
    }
    if (xo.hasAttribute(TO)) {
        to = xo.getIntegerAttribute(TO) - 1;
        if (to < 0 || to < from) {
            throw new XMLParseException("Illegal 'to' attribute in patterns element");
        }
    }
    if (every <= 0) {
        throw new XMLParseException("Illegal 'every' attribute in patterns element");
    }
    // END: every check
    TreeModel tree = (TreeModel) xo.getChild(TreeModel.class);
    GammaSiteRateModel siteModel = (GammaSiteRateModel) xo.getChild(GammaSiteRateModel.class);
    FrequencyModel freqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
    Sequence rootSequence = (Sequence) xo.getChild(Sequence.class);
    BranchRateModel rateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    if (rateModel == null) {
        rateModel = new DefaultBranchRateModel();
    }
    BranchModel branchModel = (BranchModel) xo.getChild(BranchModel.class);
    if (branchModel == null) {
        SubstitutionModel substitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
        branchModel = new HomogeneousBranchModel(substitutionModel);
    }
    Partition partition = new Partition(tree, branchModel, siteModel, rateModel, freqModel, from, to, every);
    if (rootSequence != null) {
        partition.setRootSequence(rootSequence);
    }
    return partition;
}

Example 25 with BranchRateModel

use of dr.evomodel.branchratemodel.BranchRateModel in project beast-mcmc by beast-dev.

the class BeagleSeqSimTest method simulateCodon.

// END: simulateAminoAcid
static void simulateCodon() {
    try {
        boolean calculateLikelihood = true;
        System.out.println("Test case 6: simulate codons");
        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);
        TreeModel treeModel = new TreeModel(tree);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        // create Frequency Model
        Parameter freqs = new Parameter.Default(Utils.UNIFORM_CODON_FREQUENCIES);
        FrequencyModel freqModel = new FrequencyModel(Codons.UNIVERSAL, freqs);
        // create substitution model
        Parameter alpha = new Parameter.Default(1, 10);
        Parameter beta = new Parameter.Default(1, 5);
        //			Parameter kappa = new Parameter.Default(1, 1);
        MG94CodonModel mg94 = new MG94CodonModel(Codons.UNIVERSAL, alpha, beta, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(mg94);
        // 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);
        Alignment alignment = simulator.simulate(simulateInPar, false);
        System.out.println(alignment.toString());
        if (calculateLikelihood) {
            // NewBeagleSequenceLikelihood nbtl = new
            // NewBeagleSequenceLikelihood(alignment, treeModel,
            // substitutionModel, (SiteModel) siteRateModel,
            // branchRateModel, null, false,
            // PartialsRescalingScheme.DEFAULT);
            ConvertAlignment convert = new ConvertAlignment(Nucleotides.INSTANCE, GeneticCode.UNIVERSAL, alignment);
            BeagleTreeLikelihood nbtl = new //
            BeagleTreeLikelihood(//
            convert, //
            treeModel, //
            substitutionModel, //
            siteRateModel, //
            branchRateModel, //
            null, //
            false, PartialsRescalingScheme.DEFAULT, true);
            System.out.println("likelihood = " + nbtl.getLogLikelihood());
        }
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch
}