Examples with TreeModel dr.evomodel.tree.TreeModel used on opensource projects

Example 86 with TreeModel

use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.

the class BeagleBranchLikelihood method main.

// END: LikelihoodColumn class
// ////////////
// ---TEST---//
// ////////////
public static void main(String[] args) {
    try {
        MathUtils.setSeed(666);
        int sequenceLength = 1000;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("((SimSeq1:22.0,SimSeq2:22.0):12.0,(SimSeq3:23.1,SimSeq4:23.1):10.899999999999999);");
        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);
        HKY hky1 = new HKY(kappa1, freqModel);
        BranchModel homogeneousBranchModel = new HomogeneousBranchModel(hky1);
        List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
        substitutionModels.add(hky1);
        List<FrequencyModel> freqModels = new ArrayList<FrequencyModel>();
        freqModels.add(freqModel);
        // create branch rate model
        Parameter rate = new Parameter.Default(1, 1.000);
        BranchRateModel branchRateModel = new StrictClockBranchRates(rate);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        homogeneousBranchModel, //
        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);
        System.out.println(alignment);
        BeagleTreeLikelihood btl = new BeagleTreeLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, true);
        System.out.println("BTL(homogeneous) = " + btl.getLogLikelihood());
        BeagleBranchLikelihood bbl = new BeagleBranchLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, freqModel, branchRateModel);
        int branchIndex = 4;
        System.out.println(bbl.getBranchLogLikelihood(branchIndex));
        bbl.finalizeBeagle();
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    } catch (Throwable e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 87 with TreeModel

use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.

the class LineageSpecificBranchModel method main.

// END: acceptState
public static void main(String[] args) {
    try {
        // the seed of the BEAST
        MathUtils.setSeed(666);
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        TreeModel tree = new TreeModel(importer.importTree(null));
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        // create branch rate model
        BranchRateModel branchRateModel = new DefaultBranchRateModel();
        int sequenceLength = 10;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { //
        0.0163936, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344, //
        0.01639344 });
        FrequencyModel freqModel = new FrequencyModel(Codons.UNIVERSAL, freqs);
        // create substitution model
        Parameter alpha = new Parameter.Default(1, 10);
        Parameter beta = new Parameter.Default(1, 5);
        MG94CodonModel mg94 = new MG94CodonModel(Codons.UNIVERSAL, alpha, beta, freqModel);
        HomogeneousBranchModel substitutionModel = new HomogeneousBranchModel(mg94);
        // create partition
        Partition partition1 = new //
        Partition(//
        tree, //
        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(false, false);
        ConvertAlignment convert = new ConvertAlignment(Nucleotides.INSTANCE, GeneticCode.UNIVERSAL, alignment);
        List<SubstitutionModel> substModels = new ArrayList<SubstitutionModel>();
        for (int i = 0; i < 2; i++) {
            //				alpha = new Parameter.Default(1, 10 );
            //				beta = new Parameter.Default(1, 5 );
            //				mg94 = new MG94CodonModel(Codons.UNIVERSAL, alpha, beta,
            //						freqModel);
            substModels.add(mg94);
        }
        Parameter uCategories = new Parameter.Default(2, 0);
        //            CountableBranchCategoryProvider provider = new CountableBranchCategoryProvider.IndependentBranchCategoryModel(tree, uCategories);
        LineageSpecificBranchModel branchSpecific = new //provider, 
        LineageSpecificBranchModel(//provider, 
        tree, //provider, 
        freqModel, //provider, 
        substModels, uCategories);
        BeagleTreeLikelihood like = new //
        BeagleTreeLikelihood(//
        convert, //
        tree, //
        branchSpecific, //
        siteRateModel, //
        branchRateModel, //
        null, //
        false, PartialsRescalingScheme.DEFAULT, true);
        BeagleTreeLikelihood gold = new //
        BeagleTreeLikelihood(//
        convert, //
        tree, //
        substitutionModel, //
        siteRateModel, //
        branchRateModel, //
        null, //
        false, PartialsRescalingScheme.DEFAULT, true);
        System.out.println("likelihood (gold) = " + gold.getLogLikelihood());
        System.out.println("likelihood = " + like.getLogLikelihood());
    } catch (Exception e) {
        e.printStackTrace();
    }
}

Example 88 with TreeModel

use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.

the class ContinuousTraitBranchRateModel method getBranchRate.

public double getBranchRate(final Tree tree, final NodeRef node) {
    NodeRef parent = tree.getParent(node);
    if (parent == null) {
        throw new IllegalArgumentException("Root does not have a valid rate");
    }
    double rate = 1.0;
    TreeModel treeModel = (TreeModel) tree;
    if (rateParameter != null) {
        double scale = 1.0;
        double ratio = 1.0;
        if (rateParameter != null) {
            scale = rateParameter.getParameterValue(0);
        }
        if (ratioParameter != null) {
            ratio = ratioParameter.getParameterValue(0);
        }
        // get the log rate for the node and its parent
        double rate1 = ratio * treeModel.getMultivariateNodeTrait(node, trait)[0];
        double rate2 = ratio * treeModel.getMultivariateNodeTrait(parent, trait)[0];
        if (rate1 == rate2) {
            return scale * Math.exp(rate1);
        }
        rate = scale * (Math.exp(rate2) - Math.exp(rate1)) / (rate2 - rate1);
    } else {
        double rate1 = treeModel.getMultivariateNodeTrait(node, trait)[dimension];
        double rate2 = treeModel.getMultivariateNodeTrait(parent, trait)[dimension];
        if (rate1 == rate2) {
            return Math.exp(rate1);
        }
        // TODO Should this not be averaged on the log-scale?
        rate = (Math.exp(rate2) - Math.exp(rate1)) / (rate2 - rate1);
    }
    return rate;
}

Example 89 with TreeModel

use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.

the class MicrosatelliteSingleAncestralStateGibbsOperator method doOperation.

public double doOperation() {
    TreeModel tree = msatSamplerTreeModel.getTreeModel();
    int index = MathUtils.nextInt(parameter.getDimension());
    //double logq=0.0;
    //getLikelihoodGiven the children
    NodeRef node = tree.getNode(msatSamplerTreeModel.getNodeNumberFromParameterIndex(index));
    NodeRef lc = tree.getChild(node, 0);
    NodeRef rc = tree.getChild(node, 1);
    int lcState = msatSamplerTreeModel.getNodeValue(lc);
    int rcState = msatSamplerTreeModel.getNodeValue(rc);
    double branchLeftLength = tree.getBranchLength(lc) * branchRateModel.getBranchRate(tree, lc);
    double branchRightLength = tree.getBranchLength(rc) * branchRateModel.getBranchRate(tree, rc);
    double[] probLbranch = msatModel.getColTransitionProbabilities(branchLeftLength, lcState);
    double[] probRbranch = msatModel.getColTransitionProbabilities(branchRightLength, rcState);
    double[] lik = new double[msatModel.getDataType().getStateCount()];
    if (tree.isRoot(node)) {
        //likelihood of root state also depends on the stationary distribution
        double[] statDist = msatModel.getStationaryDistribution();
        for (int j = 0; j < lik.length; j++) {
            lik[j] = probLbranch[j] * probRbranch[j] * statDist[j];
        }
    } else {
        NodeRef parent = tree.getParent(node);
        int pState = msatSamplerTreeModel.getNodeValue(parent);
        double branchParentLength = tree.getBranchLength(node) * branchRateModel.getBranchRate(tree, node);
        double[] probPbranch = msatModel.getRowTransitionProbabilities(branchParentLength, pState);
        for (int j = 0; j < lik.length; j++) {
            lik[j] = probLbranch[j] * probRbranch[j] * probPbranch[j];
        }
    }
    int sampledState = MathUtils.randomChoicePDF(lik);
    //logq = logq + Math.log(lik[currState]) - Math.log(lik[sampledState]);
    parameter.setParameterValue(index, sampledState);
    return 0.0;
}

Example 90 with TreeModel

use of dr.evomodel.tree.TreeModel in project beast-mcmc by beast-dev.

the class MicrosatelliteFullAncestryImportanceSamplingOperator method doOperation.

public double doOperation() {
    TreeModel tree = msatSamplerTreeModel.getTreeModel();
    //get postOrder
    int[] postOrder = new int[tree.getNodeCount()];
    TreeUtils.postOrderTraversalList(tree, postOrder);
    int extNodeCount = tree.getExternalNodeCount();
    double logq = 0.0;
    for (int i = 0; i < postOrder.length; i++) {
        //if it's an internal node
        if (postOrder[i] >= extNodeCount) {
            //getLikelihoodGiven the children
            NodeRef node = tree.getNode(postOrder[i]);
            NodeRef lc = tree.getChild(node, 0);
            NodeRef rc = tree.getChild(node, 1);
            int lcState = msatSamplerTreeModel.getNodeValue(lc);
            int rcState = msatSamplerTreeModel.getNodeValue(rc);
            double branchLeftLength = tree.getBranchLength(lc) * branchRateModel.getBranchRate(tree, lc);
            double branchRightLength = tree.getBranchLength(rc) * branchRateModel.getBranchRate(tree, rc);
            double[] probLbranch = msatModel.getColTransitionProbabilities(branchLeftLength, lcState);
            double[] probRbranch = msatModel.getColTransitionProbabilities(branchRightLength, rcState);
            double[] lik = new double[msatModel.getDataType().getStateCount()];
            int currState = (int) parameter.getParameterValue(msatSamplerTreeModel.getParameterIndexFromNodeNumber(postOrder[i]));
            //if node = root node
            if (i == postOrder.length - 1) {
                //likelihood of root state also depends on the stationary distribution
                double[] statDist = msatModel.getStationaryDistribution();
                for (int j = 0; j < lik.length; j++) {
                    lik[j] = probLbranch[j] * probRbranch[j] * statDist[j];
                }
            } else {
                for (int j = 0; j < lik.length; j++) {
                    lik[j] = probLbranch[j] * probRbranch[j];
                }
            }
            int sampledState = MathUtils.randomChoicePDF(lik);
            logq = logq + Math.log(lik[currState]) - Math.log(lik[sampledState]);
            parameter.setParameterValue(msatSamplerTreeModel.getParameterIndexFromNodeNumber(postOrder[i]), sampledState);
        }
    }
    return logq;
}