Examples with BranchModel dr.evomodel.branchmodel.BranchModel used on opensource projects

Example 6 with BranchModel

use of dr.evomodel.branchmodel.BranchModel in project beast-mcmc by beast-dev.

the class BeagleTreeLikelihood method main.

public static void main(String[] args) {
    try {
        MathUtils.setSeed(666);
        System.out.println("Test case 1: simulateOnePartition");
        int sequenceLength = 1000;
        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 branch model
        Parameter kappa1 = new Parameter.Default(1, 1);
        Parameter kappa2 = new Parameter.Default(1, 1);
        HKY hky1 = new HKY(kappa1, freqModel);
        HKY hky2 = new HKY(kappa2, freqModel);
        HomogeneousBranchModel homogenousBranchSubstitutionModel = new HomogeneousBranchModel(hky1);
        List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
        substitutionModels.add(hky1);
        substitutionModels.add(hky2);
        List<FrequencyModel> freqModels = new ArrayList<FrequencyModel>();
        freqModels.add(freqModel);
        Parameter epochTimes = new Parameter.Default(1, 20);
        // create branch rate model
        Parameter rate = new Parameter.Default(1, 0.001);
        BranchRateModel branchRateModel = new StrictClockBranchRates(rate);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        BranchModel homogeneousBranchModel = new HomogeneousBranchModel(hky1);
        BranchModel epochBranchModel = new EpochBranchModel(treeModel, substitutionModels, epochTimes);
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        homogenousBranchSubstitutionModel, //
        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(false, false);
        BeagleTreeLikelihood nbtl = new BeagleTreeLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
        System.out.println("nBTL(homogeneous) = " + nbtl.getLogLikelihood());
        nbtl = new BeagleTreeLikelihood(alignment, treeModel, epochBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
        System.out.println("nBTL(epoch) = " + nbtl.getLogLikelihood());
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 7 with BranchModel

use of dr.evomodel.branchmodel.BranchModel in project beast-mcmc by beast-dev.

the class PartitionData method createBranchModel.

public BranchModel createBranchModel() {
    BranchModel branchModel = null;
    if (this.substitutionModelIndex == 0) {
        // HKY
        Parameter kappa = new Parameter.Default(1, substitutionParameterValues[0]);
        FrequencyModel frequencyModel = this.createFrequencyModel();
        HKY hky = new HKY(kappa, frequencyModel);
        branchModel = new HomogeneousBranchModel(hky);
    } else if (this.substitutionModelIndex == 1) {
        // GTR
        Parameter ac = new Parameter.Default(1, substitutionParameterValues[1]);
        Parameter ag = new Parameter.Default(1, substitutionParameterValues[2]);
        Parameter at = new Parameter.Default(1, substitutionParameterValues[3]);
        Parameter cg = new Parameter.Default(1, substitutionParameterValues[4]);
        Parameter ct = new Parameter.Default(1, substitutionParameterValues[5]);
        Parameter gt = new Parameter.Default(1, substitutionParameterValues[6]);
        FrequencyModel frequencyModel = this.createFrequencyModel();
        GTR gtr = new GTR(ac, ag, at, cg, ct, gt, frequencyModel);
        branchModel = new HomogeneousBranchModel(gtr);
    } else if (this.substitutionModelIndex == 2) {
        // TN93
        Parameter kappa1 = new Parameter.Default(1, substitutionParameterValues[7]);
        Parameter kappa2 = new Parameter.Default(1, substitutionParameterValues[8]);
        FrequencyModel frequencyModel = this.createFrequencyModel();
        TN93 tn93 = new TN93(kappa1, kappa2, frequencyModel);
        branchModel = new HomogeneousBranchModel(tn93);
    } else if (this.substitutionModelIndex == 3) {
        // Yang Codon Model
        FrequencyModel frequencyModel = this.createFrequencyModel();
        Parameter kappa = new Parameter.Default(1, substitutionParameterValues[9]);
        Parameter omega = new Parameter.Default(1, substitutionParameterValues[10]);
        GY94CodonModel yangCodonModel = new GY94CodonModel(Codons.UNIVERSAL, omega, kappa, frequencyModel);
        branchModel = new HomogeneousBranchModel(yangCodonModel);
    } else if (this.substitutionModelIndex == 4) {
        // MG94CodonModel
        FrequencyModel frequencyModel = this.createFrequencyModel();
        Parameter alpha = new Parameter.Default(1, substitutionParameterValues[11]);
        Parameter beta = new Parameter.Default(1, substitutionParameterValues[12]);
        Parameter kappa = new Parameter.Default(1, substitutionParameterValues[13]);
        MG94HKYCodonModel mg94 = new MG94HKYCodonModel(Codons.UNIVERSAL, alpha, beta, kappa, frequencyModel);
        branchModel = new HomogeneousBranchModel(mg94);
    } else if (this.substitutionModelIndex == 5) {
        // Blosum62
        FrequencyModel frequencyModel = this.createFrequencyModel();
        EmpiricalRateMatrix rateMatrix = Blosum62.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, frequencyModel);
        branchModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
    } else if (this.substitutionModelIndex == 6) {
        // CPREV
        FrequencyModel frequencyModel = this.createFrequencyModel();
        EmpiricalRateMatrix rateMatrix = CPREV.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, frequencyModel);
        branchModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
    } else if (this.substitutionModelIndex == 7) {
        // Dayhoff
        FrequencyModel frequencyModel = this.createFrequencyModel();
        EmpiricalRateMatrix rateMatrix = Dayhoff.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, frequencyModel);
        branchModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
    } else if (this.substitutionModelIndex == 8) {
        // JTT
        FrequencyModel frequencyModel = this.createFrequencyModel();
        EmpiricalRateMatrix rateMatrix = JTT.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, frequencyModel);
        branchModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
    } else if (this.substitutionModelIndex == 9) {
        // LG
        FrequencyModel frequencyModel = this.createFrequencyModel();
        EmpiricalRateMatrix rateMatrix = LG.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, frequencyModel);
        branchModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
    } else if (this.substitutionModelIndex == 10) {
        // MTREV
        FrequencyModel frequencyModel = this.createFrequencyModel();
        EmpiricalRateMatrix rateMatrix = MTREV.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, frequencyModel);
        branchModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
    } else if (this.substitutionModelIndex == 11) {
        // WAG
        FrequencyModel frequencyModel = this.createFrequencyModel();
        EmpiricalRateMatrix rateMatrix = WAG.INSTANCE;
        EmpiricalAminoAcidModel empiricalAminoAcidModel = new EmpiricalAminoAcidModel(rateMatrix, frequencyModel);
        branchModel = new HomogeneousBranchModel(empiricalAminoAcidModel);
    } else {
        System.out.println("Not yet implemented");
    }
    return branchModel;
}

Example 8 with BranchModel

use of dr.evomodel.branchmodel.BranchModel in project beast-mcmc by beast-dev.

the class BranchSpecificTraitParser method parseXMLObject.

@Override
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    BranchModel branchModel = (BranchModel) xo.getChild(BranchModel.class);
    //		CompoundParameter parameter = (CompoundParameter) xo.getChild(CompoundParameter.class);
    TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
    return new BranchSpecificTrait(treeModel, branchModel, xo.getId());
}

Example 9 with BranchModel

use of dr.evomodel.branchmodel.BranchModel 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 10 with BranchModel

use of dr.evomodel.branchmodel.BranchModel in project beast-mcmc by beast-dev.

the class AncestralStateBeagleTreeLikelihoodTest method testJointLikelihood.

public void testJointLikelihood() {
    TreeModel treeModel = new TreeModel("treeModel", tree);
    Sequence[] sequence = new Sequence[3];
    sequence[0] = new Sequence(new Taxon("0"), "A");
    sequence[1] = new Sequence(new Taxon("1"), "C");
    sequence[2] = new Sequence(new Taxon("2"), "C");
    Taxa taxa = new Taxa();
    for (Sequence s : sequence) {
        taxa.addTaxon(s.getTaxon());
    }
    SimpleAlignment alignment = new SimpleAlignment();
    for (Sequence s : sequence) {
        alignment.addSequence(s);
    }
    Parameter mu = new Parameter.Default(1, 1.0);
    Parameter kappa = new Parameter.Default(1, 1.0);
    double[] pi = { 0.25, 0.25, 0.25, 0.25 };
    Parameter freqs = new Parameter.Default(pi);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    GammaSiteRateModel siteRateModel = new GammaSiteRateModel("gammaModel", mu, null, -1, null);
    siteRateModel.setSubstitutionModel(hky);
    BranchModel branchModel = new HomogeneousBranchModel(siteRateModel.getSubstitutionModel());
    BranchRateModel branchRateModel = null;
    AncestralStateBeagleTreeLikelihood treeLikelihood = new AncestralStateBeagleTreeLikelihood(alignment, treeModel, branchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, true, null, hky.getDataType(), "stateTag", // useMap = true
    true, false);
    double logLike = treeLikelihood.getLogLikelihood();
    StringBuffer buffer = new StringBuffer();
    //        Tree.Utils.newick(treeModel, treeModel.getRoot(), false, Tree.BranchLengthType.LENGTHS_AS_TIME,
    //                null, null, new NodeAttributeProvider[]{treeLikelihood}, null, null, buffer);
    TreeUtils.newick(treeModel, treeModel.getRoot(), false, TreeUtils.BranchLengthType.LENGTHS_AS_TIME, null, null, new TreeTraitProvider[] { treeLikelihood }, null, buffer);
    System.out.println(buffer);
    System.out.println("t_CA(2) = " + t(false, 2.0));
    System.out.println("t_CC(1) = " + t(true, 1.0));
    double trueValue = 0.25 * t(false, 2.0) * Math.pow(t(true, 1.0), 3.0);
    assertEquals(logLike, Math.log(trueValue), 1e-6);
}