Examples with Tree dr.evolution.tree.Tree used on opensource projects

Example 46 with Tree

use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.

the class PLCoalescentSimulator method main.

public static void main(String[] arg) throws IOException {
    // READ DEMOGRAPHIC FUNCTION
    String filename = arg[0];
    BufferedReader reader = new BufferedReader(new FileReader(filename));
    double popSizeScale = 1.0;
    double generationTime = 1.0;
    if (arg.length > 2) {
        popSizeScale = Double.parseDouble(arg[2]);
    }
    if (arg.length > 3) {
        generationTime = Double.parseDouble(arg[3]);
    }
    PrintWriter populationFuncLogger = null;
    if (arg.length > 5) {
        String logFileName = arg[5];
        if (logFileName.equals("-")) {
            populationFuncLogger = new PrintWriter(System.out);
        } else {
            populationFuncLogger = new PrintWriter(new FileWriter(logFileName));
        }
    }
    List<Double> times = new ArrayList<Double>();
    String line = reader.readLine();
    String[] tokens = line.trim().split("[\t ]+");
    if (tokens.length < 2)
        throw new RuntimeException();
    ArrayList<ArrayList> popSizes = new ArrayList<ArrayList>();
    while (line != null) {
        double time = Double.parseDouble(tokens[0]) / generationTime;
        times.add(time);
        for (int i = 1; i < tokens.length; i++) {
            popSizes.add(new ArrayList<Double>());
            popSizes.get(i - 1).add(Double.parseDouble(tokens[i]));
        }
        line = reader.readLine();
        if (line != null) {
            tokens = line.trim().split("[\t ]+");
            if (tokens.length != popSizes.size() + 1)
                throw new RuntimeException();
        }
    }
    reader.close();
    // READ SAMPLE TIMES
    String samplesFilename = arg[1];
    reader = new BufferedReader(new FileReader(samplesFilename));
    line = reader.readLine();
    Taxa taxa = new Taxa();
    int id = 0;
    while (line != null) {
        if (!line.startsWith("#")) {
            tokens = line.split("[\t ]+");
            if (tokens.length == 4) {
                double t0 = Double.parseDouble(tokens[0]);
                double t1 = Double.parseDouble(tokens[1]);
                double dt = Double.parseDouble(tokens[2]);
                int k = Integer.parseInt(tokens[3]);
                for (double time = t0; time <= t1; time += dt) {
                    double sampleTime = time / generationTime;
                    for (int i = 0; i < k; i++) {
                        Taxon taxon = new Taxon("t" + id);
                        taxon.setAttribute(dr.evolution.util.Date.DATE, new Date(sampleTime, Units.Type.GENERATIONS, true));
                        taxa.addTaxon(taxon);
                        id += 1;
                    }
                }
            } else {
                // sample times are in the same units as simulation
                double sampleTime = Double.parseDouble(tokens[0]) / generationTime;
                int count = Integer.parseInt(tokens[1]);
                for (int i = 0; i < count; i++) {
                    Taxon taxon = new Taxon(id + "");
                    taxon.setAttribute(dr.evolution.util.Date.DATE, new Date(sampleTime, Units.Type.GENERATIONS, true));
                    taxa.addTaxon(taxon);
                    id += 1;
                }
            }
        }
        line = reader.readLine();
    }
    double minTheta = Double.MAX_VALUE;
    double maxTheta = 0.0;
    PrintWriter out;
    if (arg.length > 4) {
        out = new PrintWriter(new FileWriter(arg[4]));
    } else {
        out = new PrintWriter(System.out);
    }
    int pp = 0;
    for (List<Double> popSize : popSizes) {
        double[] thetas = new double[popSize.size()];
        double[] intervals = new double[times.size() - 1];
        if (populationFuncLogger != null) {
            populationFuncLogger.println("# " + pp);
            ++pp;
        }
        // must reverse the direction of the model
        for (int j = intervals.length; j > 0; j--) {
            intervals[intervals.length - j] = times.get(j) - times.get(j - 1);
            final double theta = popSize.get(j) * popSizeScale;
            thetas[intervals.length - j] = theta;
            if (theta < minTheta) {
                minTheta = theta;
            }
            if (theta > maxTheta) {
                maxTheta = theta;
            }
            final double t = times.get(intervals.length) - times.get(j);
            if (populationFuncLogger != null) {
                populationFuncLogger.println(t + "\t" + theta);
            }
        }
        if (debug != null) {
            debug.println("min theta = " + minTheta);
            debug.println("max theta = " + maxTheta);
        }
        PiecewiseLinearPopulation demo = new PiecewiseLinearPopulation(intervals, thetas, Units.Type.GENERATIONS);
        CoalescentSimulator simulator = new CoalescentSimulator();
        Tree tree = simulator.simulateTree(taxa, demo);
        out.println(TreeUtils.newick(tree));
        if (debug != null) {
            debug.println(TreeUtils.newick(tree));
        }
    }
    if (populationFuncLogger != null) {
        populationFuncLogger.flush();
        populationFuncLogger.close();
    }
    out.flush();
    out.close();
}

Example 47 with Tree

use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.

the class NewickImporter method importTrees.

/**
     * importTrees.
     */
public Tree[] importTrees(TaxonList taxonList) throws IOException, ImportException {
    boolean done = false;
    ArrayList<FlexibleTree> array = new ArrayList<FlexibleTree>();
    do {
        try {
            skipUntil("(");
            unreadCharacter('(');
            FlexibleNode root = readInternalNode(taxonList);
            FlexibleTree tree = new FlexibleTree(root, false, true);
            array.add(tree);
            if (taxonList == null) {
                taxonList = tree;
            }
            if (readCharacter() != ';') {
                throw new BadFormatException("Expecting ';' after tree");
            }
        } catch (EOFException e) {
            done = true;
        }
    } while (!done);
    Tree[] trees = new Tree[array.size()];
    array.toArray(trees);
    return trees;
}

Example 48 with Tree

use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.

the class PartitionedTreeModelParser method parseXMLObject.

/**
     * @return a tree object based on the XML element it was passed.
     */
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    Tree tree = (Tree) xo.getChild(Tree.class);
    AbstractOutbreak outbreak = (AbstractOutbreak) xo.getElementFirstChild(OUTBREAK);
    PartitionedTreeModel treeModel;
    if (xo.hasAttribute(STARTING_TT_FILE)) {
        treeModel = new PartitionedTreeModel(xo.getId(), tree, outbreak, xo.getStringAttribute(STARTING_TT_FILE));
    } else {
        treeModel = new PartitionedTreeModel(xo.getId(), tree, outbreak);
    }
    Logger.getLogger("dr.evomodel").info("Creating the partitioned tree model, '" + xo.getId() + "'");
    for (int i = 0; i < xo.getChildCount(); i++) {
        if (xo.getChild(i) instanceof XMLObject) {
            XMLObject cxo = (XMLObject) xo.getChild(i);
            if (cxo.getName().equals(ROOT_HEIGHT)) {
                ParameterParser.replaceParameter(cxo, treeModel.getRootHeightParameter());
            } else if (cxo.getName().equals(LEAF_HEIGHT)) {
                String taxonName;
                if (cxo.hasAttribute(TAXON)) {
                    taxonName = cxo.getStringAttribute(TAXON);
                } else {
                    throw new XMLParseException("taxa element missing from leafHeight element in treeModel element");
                }
                int index = treeModel.getTaxonIndex(taxonName);
                if (index == -1) {
                    throw new XMLParseException("taxon " + taxonName + " not found for leafHeight element in treeModel element");
                }
                NodeRef node = treeModel.getExternalNode(index);
                Parameter newParameter = treeModel.getLeafHeightParameter(node);
                ParameterParser.replaceParameter(cxo, newParameter);
                Taxon taxon = treeModel.getTaxon(index);
                setPrecisionBounds(newParameter, taxon);
            } else if (cxo.getName().equals(LEAF_HEIGHTS)) {
                // get a set of leaf height parameters out as a compound parameter...
                TaxonList taxa = (TaxonList) cxo.getChild(TaxonList.class);
                Parameter offsetParameter = (Parameter) cxo.getChild(Parameter.class);
                CompoundParameter leafHeights = new CompoundParameter("leafHeights");
                for (Taxon taxon : taxa) {
                    int index = treeModel.getTaxonIndex(taxon);
                    if (index == -1) {
                        throw new XMLParseException("taxon " + taxon.getId() + " not found for leafHeight element in treeModel element");
                    }
                    NodeRef node = treeModel.getExternalNode(index);
                    Parameter newParameter = treeModel.getLeafHeightParameter(node);
                    leafHeights.addParameter(newParameter);
                    setPrecisionBounds(newParameter, taxon);
                }
                ParameterParser.replaceParameter(cxo, leafHeights);
            } else if (cxo.getName().equals(NODE_HEIGHTS)) {
                boolean rootNode = cxo.getAttribute(ROOT_NODE, false);
                boolean internalNodes = cxo.getAttribute(INTERNAL_NODES, false);
                boolean leafNodes = cxo.getAttribute(LEAF_NODES, false);
                if (!rootNode && !internalNodes && !leafNodes) {
                    throw new XMLParseException("one or more of root, internal or leaf nodes must be selected for the nodeHeights element");
                }
                ParameterParser.replaceParameter(cxo, treeModel.createNodeHeightsParameter(rootNode, internalNodes, leafNodes));
            } else if (cxo.getName().equals(NODE_RATES)) {
                boolean rootNode = cxo.getAttribute(ROOT_NODE, false);
                boolean internalNodes = cxo.getAttribute(INTERNAL_NODES, false);
                boolean leafNodes = cxo.getAttribute(LEAF_NODES, false);
                double[] initialValues = null;
                if (cxo.hasAttribute(INITIAL_VALUE)) {
                    initialValues = cxo.getDoubleArrayAttribute(INITIAL_VALUE);
                }
                if (!rootNode && !internalNodes && !leafNodes) {
                    throw new XMLParseException("one or more of root, internal or leaf nodes must be selected for the nodeRates element");
                }
                ParameterParser.replaceParameter(cxo, treeModel.createNodeRatesParameter(initialValues, rootNode, internalNodes, leafNodes));
            } else if (cxo.getName().equals(NODE_TRAITS)) {
                boolean rootNode = cxo.getAttribute(ROOT_NODE, false);
                boolean internalNodes = cxo.getAttribute(INTERNAL_NODES, false);
                boolean leafNodes = cxo.getAttribute(LEAF_NODES, false);
                boolean fireTreeEvents = cxo.getAttribute(FIRE_TREE_EVENTS, false);
                String name = cxo.getAttribute(NAME, "trait");
                int dim = cxo.getAttribute(MULTIVARIATE_TRAIT, 1);
                double[] initialValues = null;
                if (cxo.hasAttribute(INITIAL_VALUE)) {
                    initialValues = cxo.getDoubleArrayAttribute(INITIAL_VALUE);
                }
                if (!rootNode && !internalNodes && !leafNodes) {
                    throw new XMLParseException("one or more of root, internal or leaf nodes must be selected for the nodeTraits element");
                }
                ParameterParser.replaceParameter(cxo, treeModel.createNodeTraitsParameter(name, dim, initialValues, rootNode, internalNodes, leafNodes, fireTreeEvents));
            } else if (cxo.getName().equals(LEAF_TRAIT)) {
                String name = cxo.getAttribute(NAME, "trait");
                String taxonName;
                if (cxo.hasAttribute(TAXON)) {
                    taxonName = cxo.getStringAttribute(TAXON);
                } else {
                    throw new XMLParseException("taxa element missing from leafTrait element in treeModel element");
                }
                int index = treeModel.getTaxonIndex(taxonName);
                if (index == -1) {
                    throw new XMLParseException("taxon '" + taxonName + "' not found for leafTrait element in treeModel element");
                }
                NodeRef node = treeModel.getExternalNode(index);
                Parameter parameter = treeModel.getNodeTraitParameter(node, name);
                if (parameter == null)
                    throw new XMLParseException("trait '" + name + "' not found for leafTrait (taxon, " + taxonName + ") element in treeModel element");
                ParameterParser.replaceParameter(cxo, parameter);
            } else {
                if (!cxo.getName().equals(OUTBREAK)) {
                    throw new XMLParseException("illegal child element in " + getParserName() + ": " + cxo.getName());
                }
            }
        } else if (xo.getChild(i) instanceof Tree) {
        // do nothing - already handled
        } else {
            throw new XMLParseException("illegal child element in  " + getParserName() + ": " + xo.getChildName(i) + " " + xo.getChild(i));
        }
    }
    // AR this is doubling up the number of bounds on each node.
    //        treeModel.setupHeightBounds();
    //System.err.println("done constructing treeModel");
    Logger.getLogger("dr.evomodel").info("  initial tree topology = " + TreeUtils.uniqueNewick(treeModel, treeModel.getRoot()));
    Logger.getLogger("dr.evomodel").info("  tree height = " + treeModel.getNodeHeight(treeModel.getRoot()));
    return treeModel;
}

Example 49 with Tree

use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.

the class ModelAveragingSpeciationLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    List<Tree> trees = new ArrayList<Tree>();
    List<MaskableSpeciationModel> models = new ArrayList<MaskableSpeciationModel>();
    Variable<Integer> index;
    System.out.println("id = " + xo.getId());
    XMLObject cxo = xo.getChild(MODEL);
    for (int m = 0; m < cxo.getChildCount(); m++) {
        final MaskableSpeciationModel specModel = (MaskableSpeciationModel) cxo.getChild(m);
        models.add(specModel);
    }
    cxo = xo.getChild(TREE);
    for (int t = 0; t < cxo.getChildCount(); t++) {
        final Tree tree = (Tree) cxo.getChild(t);
        trees.add(tree);
    }
    //        cxo = xo.getChild(INDEX);
    // integer index parameter size = real size - 1
    index = (Variable<Integer>) xo.getElementFirstChild(INDEX);
    Parameter maxIndex = (Parameter) xo.getElementFirstChild(MAX_INDEX);
    //        System.out.println(index.getClass());
    //        for (int i=0; i<index.getSize(); i++) {
    //            System.out.println(index.getValue(i).getClass());
    //            index.setValue(i, 0);
    //        }
    int indexLength = models.size();
    if (indexLength < 1 || trees.size() < 1) {
        throw new XMLParseException("It requires at least one tree or one speciation model.");
    } else if (indexLength != trees.size()) {
        throw new XMLParseException("The number of trees and the number of speciation models should be equal.");
    } else if (indexLength != index.getSize() + 1) {
        // integer index parameter size = real size - 1
        throw new XMLParseException("Index parameter must be same size as the number of trees.");
    }
    Logger.getLogger("dr.evomodel").info("Speciation model excluding " + " taxa remaining.");
    return new ModelAveragingSpeciationLikelihood(trees, models, index, maxIndex, xo.getId());
}

Example 50 with Tree

use of dr.evolution.tree.Tree in project beast-mcmc by beast-dev.

the class SpeciationLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    XMLObject cxo = xo.getChild(MODEL);
    final SpeciationModel specModel = (SpeciationModel) cxo.getChild(SpeciationModel.class);
    cxo = xo.getChild(TREE);
    final Tree tree = (Tree) cxo.getChild(Tree.class);
    Set<Taxon> excludeTaxa = null;
    if (xo.hasChildNamed(INCLUDE)) {
        excludeTaxa = new HashSet<Taxon>();
        for (int i = 0; i < tree.getTaxonCount(); i++) {
            excludeTaxa.add(tree.getTaxon(i));
        }
        cxo = xo.getChild(INCLUDE);
        for (int i = 0; i < cxo.getChildCount(); i++) {
            TaxonList taxonList = (TaxonList) cxo.getChild(i);
            for (int j = 0; j < taxonList.getTaxonCount(); j++) {
                excludeTaxa.remove(taxonList.getTaxon(j));
            }
        }
    }
    if (xo.hasChildNamed(EXCLUDE)) {
        excludeTaxa = new HashSet<Taxon>();
        cxo = xo.getChild(EXCLUDE);
        for (int i = 0; i < cxo.getChildCount(); i++) {
            TaxonList taxonList = (TaxonList) cxo.getChild(i);
            for (int j = 0; j < taxonList.getTaxonCount(); j++) {
                excludeTaxa.add(taxonList.getTaxon(j));
            }
        }
    }
    if (excludeTaxa != null) {
        Logger.getLogger("dr.evomodel").info("Speciation model excluding " + excludeTaxa.size() + " taxa from prior - " + (tree.getTaxonCount() - excludeTaxa.size()) + " taxa remaining.");
    }
    final XMLObject cal = xo.getChild(CALIBRATION);
    if (cal != null) {
        if (excludeTaxa != null) {
            throw new XMLParseException("Sorry, not implemented: internal calibration prior + excluded taxa");
        }
        List<Distribution> dists = new ArrayList<Distribution>();
        List<Taxa> taxa = new ArrayList<Taxa>();
        List<Boolean> forParent = new ArrayList<Boolean>();
        // (Statistic) cal.getChild(Statistic.class);
        Statistic userPDF = null;
        for (int k = 0; k < cal.getChildCount(); ++k) {
            final Object ck = cal.getChild(k);
            if (DistributionLikelihood.class.isInstance(ck)) {
                dists.add(((DistributionLikelihood) ck).getDistribution());
            } else if (Distribution.class.isInstance(ck)) {
                dists.add((Distribution) ck);
            } else if (Taxa.class.isInstance(ck)) {
                final Taxa tx = (Taxa) ck;
                taxa.add(tx);
                forParent.add(tx.getTaxonCount() == 1);
            } else if (Statistic.class.isInstance(ck)) {
                if (userPDF != null) {
                    throw new XMLParseException("more than one userPDF correction???");
                }
                userPDF = (Statistic) cal.getChild(Statistic.class);
            } else {
                XMLObject cko = (XMLObject) ck;
                assert cko.getChildCount() == 2;
                for (int i = 0; i < 2; ++i) {
                    final Object chi = cko.getChild(i);
                    if (DistributionLikelihood.class.isInstance(chi)) {
                        dists.add(((DistributionLikelihood) chi).getDistribution());
                    } else if (Distribution.class.isInstance(chi)) {
                        dists.add((Distribution) chi);
                    } else if (Taxa.class.isInstance(chi)) {
                        taxa.add((Taxa) chi);
                        boolean fp = ((Taxa) chi).getTaxonCount() == 1;
                        if (cko.hasAttribute(PARENT)) {
                            boolean ufp = cko.getBooleanAttribute(PARENT);
                            if (fp && !ufp) {
                                throw new XMLParseException("forParent==false for a single taxon?? (must be true)");
                            }
                            fp = ufp;
                        }
                        forParent.add(fp);
                    } else {
                        assert false;
                    }
                }
            }
        }
        if (dists.size() != taxa.size()) {
            throw new XMLParseException("Mismatch in number of distributions and taxa specs");
        }
        try {
            final String correction = cal.getAttribute(CORRECTION, EXACT);
            final CalibrationPoints.CorrectionType type = correction.equals(EXACT) ? CalibrationPoints.CorrectionType.EXACT : (correction.equals(APPROX) ? CalibrationPoints.CorrectionType.APPROXIMATED : (correction.equals(NONE) ? CalibrationPoints.CorrectionType.NONE : (correction.equals(PEXACT) ? CalibrationPoints.CorrectionType.PEXACT : null)));
            if (cal.hasAttribute(CORRECTION) && type == null) {
                throw new XMLParseException("correction type == " + correction + "???");
            }
            final CalibrationPoints calib = new CalibrationPoints(tree, specModel.isYule(), dists, taxa, forParent, userPDF, type);
            final SpeciationLikelihood speciationLikelihood = new SpeciationLikelihood(tree, specModel, null, calib);
            return speciationLikelihood;
        } catch (IllegalArgumentException e) {
            throw new XMLParseException(e.getMessage());
        }
    }
    return new SpeciationLikelihood(tree, specModel, excludeTaxa, null);
}