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

Example 81 with TreeModel

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

the class BeastCheckpointer method writeStateToFile.

private boolean writeStateToFile(File file, long state, double lnL, MarkovChain markovChain) {
    OperatorSchedule operatorSchedule = markovChain.getSchedule();
    OutputStream fileOut = null;
    try {
        fileOut = new FileOutputStream(file);
        PrintStream out = new PrintStream(fileOut);
        ArrayList<TreeParameterModel> traitModels = new ArrayList<TreeParameterModel>();
        int[] rngState = MathUtils.getRandomState();
        out.print("rng");
        for (int i = 0; i < rngState.length; i++) {
            out.print("\t");
            out.print(rngState[i]);
        }
        out.println();
        out.print("state\t");
        out.println(state);
        out.print("lnL\t");
        out.println(lnL);
        for (Parameter parameter : Parameter.CONNECTED_PARAMETER_SET) {
            out.print("parameter");
            out.print("\t");
            out.print(parameter.getParameterName());
            out.print("\t");
            out.print(parameter.getDimension());
            for (int dim = 0; dim < parameter.getDimension(); dim++) {
                out.print("\t");
                out.print(parameter.getParameterValue(dim));
            }
            out.println();
        }
        for (int i = 0; i < operatorSchedule.getOperatorCount(); i++) {
            MCMCOperator operator = operatorSchedule.getOperator(i);
            out.print("operator");
            out.print("\t");
            out.print(operator.getOperatorName());
            out.print("\t");
            out.print(operator.getAcceptCount());
            out.print("\t");
            out.print(operator.getRejectCount());
            if (operator instanceof CoercableMCMCOperator) {
                out.print("\t");
                out.print(((CoercableMCMCOperator) operator).getCoercableParameter());
            }
            out.println();
        }
        //check up front if there are any TreeParameterModel objects
        for (Model model : Model.CONNECTED_MODEL_SET) {
            if (model instanceof TreeParameterModel) {
                //System.out.println("\nDetected TreeParameterModel: " + ((TreeParameterModel) model).toString());
                traitModels.add((TreeParameterModel) model);
            }
        }
        for (Model model : Model.CONNECTED_MODEL_SET) {
            if (model instanceof TreeModel) {
                out.print("tree");
                out.print("\t");
                out.println(model.getModelName());
                //replace Newick format by printing general graph structure
                //out.println(((TreeModel) model).getNewick());
                out.println("#node height taxon");
                int nodeCount = ((TreeModel) model).getNodeCount();
                out.println(nodeCount);
                for (int i = 0; i < nodeCount; i++) {
                    out.print(((TreeModel) model).getNode(i).getNumber());
                    out.print("\t");
                    out.print(((TreeModel) model).getNodeHeight(((TreeModel) model).getNode(i)));
                    if (((TreeModel) model).isExternal(((TreeModel) model).getNode(i))) {
                        out.print("\t");
                        out.print(((TreeModel) model).getNodeTaxon(((TreeModel) model).getNode(i)).getId());
                    }
                    out.println();
                }
                out.println("#edges");
                out.println("#child-node parent-node L/R-child traits");
                out.println(nodeCount);
                for (int i = 0; i < nodeCount; i++) {
                    NodeRef parent = ((TreeModel) model).getParent(((TreeModel) model).getNode(i));
                    if (parent != null) {
                        out.print(((TreeModel) model).getNode(i).getNumber());
                        out.print("\t");
                        out.print(((TreeModel) model).getParent(((TreeModel) model).getNode(i)).getNumber());
                        out.print("\t");
                        if ((((TreeModel) model).getChild(parent, 0) == ((TreeModel) model).getNode(i))) {
                            //left child
                            out.print(0);
                        } else if ((((TreeModel) model).getChild(parent, 1) == ((TreeModel) model).getNode(i))) {
                            //right child
                            out.print(1);
                        } else {
                            throw new RuntimeException("Operation currently only supported for nodes with 2 children.");
                        }
                        for (TreeParameterModel tpm : traitModels) {
                            out.print("\t");
                            out.print(tpm.getNodeValue((TreeModel) model, ((TreeModel) model).getNode(i)));
                        }
                        out.println();
                    }
                }
            }
        }
        out.close();
        fileOut.close();
    } catch (IOException ioe) {
        System.err.println("Unable to write file: " + ioe.getMessage());
        return false;
    }
    if (DEBUG) {
        for (Likelihood likelihood : Likelihood.CONNECTED_LIKELIHOOD_SET) {
            System.err.println(likelihood.getId() + ": " + likelihood.getLogLikelihood());
        }
    }
    return true;
}

Example 82 with TreeModel

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

the class BitFlipOperatorParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    double weight = xo.getDoubleAttribute(MCMCOperator.WEIGHT);
    Parameter parameter = (Parameter) xo.getChild(Parameter.class);
    boolean usesPriorOnSum = xo.getAttribute(USES_SUM_PRIOR, true);
    // boolean forDrift = xo.getAttribute(FOR_DRIFT,false);
    TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
    return new BitFlipOperator(parameter, weight, usesPriorOnSum, treeModel);
}

Example 83 with TreeModel

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

the class AdvancedTreeLikelihood method handleModelChangedEvent.

// **************************************************************
// ModelListener IMPLEMENTATION
// **************************************************************
/**
     * Handles model changed events from the submodels.
     */
protected void handleModelChangedEvent(Model model, Object object, int index) {
    if (model == treeModel) {
        if (object instanceof TreeModel.TreeChangedEvent) {
            if (((TreeModel.TreeChangedEvent) object).isNodeChanged()) {
                updateNodeAndChildren(((TreeModel.TreeChangedEvent) object).getNode());
            } else {
                updateAllNodes();
                commonAncestorsKnown = false;
            }
        }
    } else if (model == branchRateModel) {
        updateAllNodes();
    } else if (model == frequencyModel) {
        updateAllNodes();
    } else if (model instanceof SiteModel) {
        if (model == siteModel) {
            updateAllNodes();
        } else if (model == tipsSiteModel) {
            updateAllNodes();
        } else {
            // find the siteModel in the additional siteModel list
            NodeRef node = null;
            for (int i = 0, n = cladeSiteModels.size(); i < n; i++) {
                Clade clade = cladeSiteModels.get(i);
                if (!commonAncestorsKnown) {
                    clade.findMRCA();
                }
                if (clade.getSiteModel() == model) {
                    node = treeModel.getNode(clade.getNode());
                }
            }
            commonAncestorsKnown = true;
            updateNodeAndDescendents(node);
        }
    } else {
        throw new RuntimeException("Unknown componentChangedEvent");
    }
    super.handleModelChangedEvent(model, object, index);
}

Example 84 with TreeModel

use of dr.evomodel.tree.TreeModel 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);
}

Example 85 with TreeModel

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

the class TreeLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    boolean useAmbiguities = xo.getAttribute(USE_AMBIGUITIES, false);
    boolean allowMissingTaxa = xo.getAttribute(ALLOW_MISSING_TAXA, false);
    boolean storePartials = xo.getAttribute(STORE_PARTIALS, true);
    boolean forceJavaCore = xo.getAttribute(FORCE_JAVA_CORE, false);
    if (Boolean.valueOf(System.getProperty("java.only"))) {
        forceJavaCore = true;
    }
    PatternList patternList = (PatternList) xo.getChild(PatternList.class);
    TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
    SiteModel siteModel = (SiteModel) xo.getChild(SiteModel.class);
    BranchRateModel branchRateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    TipStatesModel tipStatesModel = (TipStatesModel) xo.getChild(TipStatesModel.class);
    if (tipStatesModel != null && tipStatesModel.getPatternList() != null) {
        throw new XMLParseException("The same sequence error model cannot be used for multiple partitions");
    }
    if (tipStatesModel != null && tipStatesModel.getModelType() == TipStatesModel.Type.STATES) {
        throw new XMLParseException("The state emitting TipStateModel requires BEAGLE");
    }
    boolean forceRescaling = xo.getAttribute(FORCE_RESCALING, false);
    return new TreeLikelihood(patternList, treeModel, siteModel, branchRateModel, tipStatesModel, useAmbiguities, allowMissingTaxa, storePartials, forceJavaCore, forceRescaling);
}