Search in sources :

Example 31 with TaxonList

use of dr.evolution.util.TaxonList in project beast-mcmc by beast-dev.

the class RandomLocalClockTestProblem method testRandomLocalClock.

public void testRandomLocalClock() throws Exception {
    Parameter popSize = new Parameter.Default(ConstantPopulationModelParser.POPULATION_SIZE, 0.077, 0, Double.POSITIVE_INFINITY);
    ConstantPopulationModel constantModel = createRandomInitialTree(popSize);
    CoalescentLikelihood coalescent = new CoalescentLikelihood(treeModel, null, new ArrayList<TaxonList>(), constantModel);
    coalescent.setId("coalescent");
    // clock model
    Parameter ratesParameter = new Parameter.Default(RandomLocalClockModelParser.RATES, 10, 1);
    Parameter rateIndicatorParameter = new Parameter.Default(RandomLocalClockModelParser.RATE_INDICATORS, 10, 1);
    Parameter meanRateParameter = new Parameter.Default(RandomLocalClockModelParser.CLOCK_RATE, 1, 1.0);
    RandomLocalClockModel branchRateModel = new RandomLocalClockModel(treeModel, meanRateParameter, rateIndicatorParameter, ratesParameter, false, 0.5);
    SumStatistic rateChanges = new SumStatistic("rateChangeCount", true);
    rateChanges.addStatistic(rateIndicatorParameter);
    RateStatistic meanRate = new RateStatistic("meanRate", treeModel, branchRateModel, true, true, RateStatisticParser.MEAN);
    RateStatistic coefficientOfVariation = new RateStatistic(RateStatisticParser.COEFFICIENT_OF_VARIATION, treeModel, branchRateModel, true, true, RateStatisticParser.COEFFICIENT_OF_VARIATION);
    RateCovarianceStatistic covariance = new RateCovarianceStatistic("covariance", treeModel, branchRateModel);
    // Sub model
    Parameter freqs = new Parameter.Default(alignment.getStateFrequencies());
    Parameter kappa = new Parameter.Default(HKYParser.KAPPA, 1.0, 0, Double.POSITIVE_INFINITY);
    FrequencyModel f = new FrequencyModel(Nucleotides.INSTANCE, freqs);
    HKY hky = new HKY(kappa, f);
    //siteModel
    GammaSiteModel siteModel = new GammaSiteModel(hky);
    Parameter mu = new Parameter.Default(GammaSiteModelParser.MUTATION_RATE, 1.0, 0, Double.POSITIVE_INFINITY);
    siteModel.setMutationRateParameter(mu);
    //treeLikelihood
    SitePatterns patterns = new SitePatterns(alignment, null, 0, -1, 1, true);
    TreeLikelihood treeLikelihood = new TreeLikelihood(patterns, treeModel, siteModel, branchRateModel, null, false, false, true, false, false);
    treeLikelihood.setId(TreeLikelihoodParser.TREE_LIKELIHOOD);
    // Operators
    OperatorSchedule schedule = new SimpleOperatorSchedule();
    MCMCOperator operator = new ScaleOperator(kappa, 0.75);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    operator = new ScaleOperator(ratesParameter, 0.75);
    operator.setWeight(10.0);
    schedule.addOperator(operator);
    operator = new BitFlipOperator(rateIndicatorParameter, 15.0, true);
    schedule.addOperator(operator);
    operator = new ScaleOperator(popSize, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter rootHeight = treeModel.getRootHeightParameter();
    rootHeight.setId(TREE_HEIGHT);
    operator = new ScaleOperator(rootHeight, 0.75);
    operator.setWeight(3.0);
    schedule.addOperator(operator);
    Parameter internalHeights = treeModel.createNodeHeightsParameter(false, true, false);
    operator = new UniformOperator(internalHeights, 30.0);
    schedule.addOperator(operator);
    operator = new SubtreeSlideOperator(treeModel, 15.0, 0.0077, true, false, false, false, CoercionMode.COERCION_ON);
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.NARROW, treeModel, 15.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new ExchangeOperator(ExchangeOperator.WIDE, treeModel, 3.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    operator = new WilsonBalding(treeModel, 3.0);
    //        operator.doOperation();
    schedule.addOperator(operator);
    //CompoundLikelihood
    OneOnXPrior likelihood1 = new OneOnXPrior();
    likelihood1.addData(popSize);
    OneOnXPrior likelihood2 = new OneOnXPrior();
    likelihood2.addData(kappa);
    DistributionLikelihood likelihood3 = new DistributionLikelihood(new GammaDistribution(0.5, 2.0), 0.0);
    likelihood3.addData(ratesParameter);
    DistributionLikelihood likelihood4 = new DistributionLikelihood(new PoissonDistribution(1.0), 0.0);
    likelihood4.addData(rateChanges);
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(likelihood1);
    likelihoods.add(likelihood2);
    likelihoods.add(likelihood3);
    likelihoods.add(likelihood4);
    likelihoods.add(coalescent);
    Likelihood prior = new CompoundLikelihood(0, likelihoods);
    prior.setId(CompoundLikelihoodParser.PRIOR);
    likelihoods.clear();
    likelihoods.add(treeLikelihood);
    Likelihood likelihood = new CompoundLikelihood(-1, likelihoods);
    likelihoods.clear();
    likelihoods.add(prior);
    likelihoods.add(likelihood);
    Likelihood posterior = new CompoundLikelihood(0, likelihoods);
    posterior.setId(CompoundLikelihoodParser.POSTERIOR);
    // Log
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[2];
    loggers[0] = new MCLogger(formatter, 1000, false);
    loggers[0].add(posterior);
    loggers[0].add(prior);
    loggers[0].add(treeLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(kappa);
    //        loggers[0].add(meanRate);
    loggers[0].add(rateChanges);
    loggers[0].add(coefficientOfVariation);
    loggers[0].add(covariance);
    loggers[0].add(popSize);
    loggers[0].add(coalescent);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 10000, false);
    loggers[1].add(posterior);
    loggers[1].add(treeLikelihood);
    loggers[1].add(rootHeight);
    loggers[1].add(meanRate);
    loggers[1].add(rateChanges);
    // MCMC
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(1000000);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, posterior, schedule, loggers);
    mcmc.run();
    // time
    System.out.println(mcmc.getTimer().toString());
    // Tracer
    List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("RandomLocalClockTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    //        <expectation name="posterior" value="-1818.26"/>
    //        <expectation name="prior" value="-2.70143"/>
    //        <expectation name="likelihood" value="-1815.56"/>
    //        <expectation name="treeModel.rootHeight" value="6.363E-2"/>
    //        <expectation name="constant.popSize" value="9.67405E-2"/>
    //        <expectation name="hky.kappa" value="30.0394"/>
    //        <expectation name="coefficientOfVariation" value="7.02408E-2"/>
    //        covariance 0.47952
    //        <expectation name="rateChangeCount" value="0.40786"/>
    //        <expectation name="coalescent" value="7.29521"/>
    TraceCorrelation likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(CompoundLikelihoodParser.POSTERIOR));
    assertExpectation(CompoundLikelihoodParser.POSTERIOR, likelihoodStats, -1818.26);
    likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(CompoundLikelihoodParser.PRIOR));
    assertExpectation(CompoundLikelihoodParser.PRIOR, likelihoodStats, -2.70143);
    likelihoodStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TreeLikelihoodParser.TREE_LIKELIHOOD));
    assertExpectation(TreeLikelihoodParser.TREE_LIKELIHOOD, likelihoodStats, -1815.56);
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TREE_HEIGHT));
    assertExpectation(TREE_HEIGHT, treeHeightStats, 6.363E-2);
    TraceCorrelation kappaStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(HKYParser.KAPPA));
    assertExpectation(HKYParser.KAPPA, kappaStats, 30.0394);
    TraceCorrelation rateChangeStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("rateChangeCount"));
    assertExpectation("rateChangeCount", rateChangeStats, 0.40786);
    TraceCorrelation coefficientOfVariationStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(RateStatisticParser.COEFFICIENT_OF_VARIATION));
    assertExpectation(RateStatisticParser.COEFFICIENT_OF_VARIATION, coefficientOfVariationStats, 7.02408E-2);
    TraceCorrelation covarianceStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("covariance"));
    assertExpectation("covariance", covarianceStats, 0.47952);
    TraceCorrelation popStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(ConstantPopulationModelParser.POPULATION_SIZE));
    assertExpectation(ConstantPopulationModelParser.POPULATION_SIZE, popStats, 9.67405E-2);
    TraceCorrelation coalescentStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("coalescent"));
    assertExpectation("coalescent", coalescentStats, 7.29521);
}
Also used : FrequencyModel(dr.oldevomodel.substmodel.FrequencyModel) PoissonDistribution(dr.math.distributions.PoissonDistribution) TreeLikelihood(dr.oldevomodel.treelikelihood.TreeLikelihood) DistributionLikelihood(dr.inference.distribution.DistributionLikelihood) CoalescentLikelihood(dr.evomodel.coalescent.CoalescentLikelihood) TreeLikelihood(dr.oldevomodel.treelikelihood.TreeLikelihood) ExchangeOperator(dr.evomodel.operators.ExchangeOperator) ArrayList(java.util.ArrayList) MCMC(dr.inference.mcmc.MCMC) SubtreeSlideOperator(dr.evomodel.operators.SubtreeSlideOperator) RateCovarianceStatistic(dr.evomodel.tree.RateCovarianceStatistic) CoalescentLikelihood(dr.evomodel.coalescent.CoalescentLikelihood) MCMCOptions(dr.inference.mcmc.MCMCOptions) ArrayLogFormatter(dr.inference.loggers.ArrayLogFormatter) GammaDistribution(dr.math.distributions.GammaDistribution) WilsonBalding(dr.evomodel.operators.WilsonBalding) TraceCorrelation(dr.inference.trace.TraceCorrelation) SitePatterns(dr.evolution.alignment.SitePatterns) ConstantPopulationModel(dr.evomodel.coalescent.ConstantPopulationModel) TaxonList(dr.evolution.util.TaxonList) TabDelimitedFormatter(dr.inference.loggers.TabDelimitedFormatter) Trace(dr.inference.trace.Trace) GammaSiteModel(dr.oldevomodel.sitemodel.GammaSiteModel) RateStatistic(dr.evomodel.tree.RateStatistic) ArrayTraceList(dr.inference.trace.ArrayTraceList) RandomLocalClockModel(dr.evomodel.branchratemodel.RandomLocalClockModel) HKY(dr.oldevomodel.substmodel.HKY) DistributionLikelihood(dr.inference.distribution.DistributionLikelihood) MCLogger(dr.inference.loggers.MCLogger)

Example 32 with TaxonList

use of dr.evolution.util.TaxonList in project beast-mcmc by beast-dev.

the class IndependentCoalescentSampler method doOperation.

/**
	 * change the parameter and return the hastings ratio.
     */
public double doOperation() {
    CoalescentSimulator simulator = new CoalescentSimulator();
    List<TaxonList> taxonLists = new ArrayList<TaxonList>();
    double rootHeight = -1.0;
    double oldLikelihood = 0.0;
    double newLikelihood = 0.0;
    // should have one child that is node
    for (int i = 0; i < xo.getChildCount(); i++) {
        final Object child = xo.getChild(i);
        //careful: Trees are TaxonLists ... (AER); see OldCoalescentSimulatorParser
        if (child instanceof Tree) {
        //do nothing
        } else if (child instanceof TaxonList) {
            //taxonLists.add((TaxonList) child);
            taxonLists.add((Taxa) child);
            //taxa added
            break;
        }
    }
    try {
        Tree[] trees = new Tree[taxonLists.size()];
        // simulate each taxonList separately
        for (int i = 0; i < taxonLists.size(); i++) {
            trees[i] = simulator.simulateTree(taxonLists.get(i), demoModel);
        }
        oldLikelihood = coalescent.getLogLikelihood();
        SimpleTree simTree = simulator.simulateTree(trees, demoModel, rootHeight, trees.length != 1);
        //this would be the normal way to do it
        treeModel.beginTreeEdit();
        //now it's allowed to adjust the tree structure
        treeModel.adoptTreeStructure(simTree);
        //endTreeEdit() would then fire the events
        treeModel.endTreeEdit();
        newLikelihood = coalescent.getLogLikelihood();
    } catch (IllegalArgumentException iae) {
        try {
            throw new XMLParseException(iae.getMessage());
        } catch (XMLParseException e) {
            // TODO Auto-generated catch block
            e.printStackTrace();
        }
    }
    return oldLikelihood - newLikelihood;
}
Also used : TaxonList(dr.evolution.util.TaxonList) ArrayList(java.util.ArrayList) SimpleTree(dr.evolution.tree.SimpleTree) CoalescentSimulator(dr.evomodel.coalescent.CoalescentSimulator) Taxa(dr.evolution.util.Taxa) SimpleTree(dr.evolution.tree.SimpleTree) Tree(dr.evolution.tree.Tree) XMLObject(dr.xml.XMLObject) XMLParseException(dr.xml.XMLParseException)

Example 33 with TaxonList

use of dr.evolution.util.TaxonList in project beast-mcmc by beast-dev.

the class HiddenLinkageTreeLogger method processTree.

/**
     * Create a tree of the same topology as the input tree, but with 
     * the linkage groups replaced by their constituent reads
     */
protected static SimpleTree processTree(Tree tree, HiddenLinkageModel hlm) {
    TaxonList reads = hlm.getData().getReadsTaxa();
    TaxonList reference = hlm.getData().getReferenceTaxa();
    // allocate space
    int nodeCount = tree.getTaxonCount() + reads.getTaxonCount();
    nodeCount = 2 * nodeCount - 1;
    SimpleNode[] nodes = new SimpleNode[nodeCount];
    for (int i = 0; i < nodes.length; i++) {
        nodes[i] = new SimpleNode();
        nodes[i].setNumber(i);
    }
    SimpleNode root = null;
    // copy the tree structure
    for (int i = 0; i < tree.getNodeCount(); i++) {
        NodeRef n = tree.getNode(i);
        for (int cI = 0; cI < tree.getChildCount(n); cI++) {
            NodeRef c1 = tree.getChild(n, cI);
            nodes[n.getNumber()].addChild(nodes[c1.getNumber()]);
        }
        nodes[n.getNumber()].setHeight(tree.getNodeHeight(n));
        nodes[n.getNumber()].setRate(tree.getNodeRate(n));
        nodes[n.getNumber()].setTaxon(tree.getNodeTaxon(n));
    }
    root = nodes[tree.getRoot().getNumber()];
    // now replace linkage groups with their constituent reads
    // first free up anything in the range of read leaf nodes
    int nextFree = tree.getNodeCount();
    int readI = 0;
    for (int i = reference.getTaxonCount(); i < reference.getTaxonCount() + reads.getTaxonCount(); i++) {
        SimpleNode tmp = nodes[nextFree];
        nodes[nextFree] = nodes[i];
        nodes[nextFree].setNumber(nextFree);
        nodes[i] = tmp;
        nodes[i].setNumber(i);
        nodes[i].setTaxon(reads.getTaxon(readI));
        readI++;
        nextFree++;
    }
    // if the linkage group has many reads, build a ladder
    for (int i = 0; i < nodes.length; i++) {
        SimpleNode n = nodes[i];
        if (n.getTaxon() == null)
            continue;
        if (reads.getTaxonIndex(n.getTaxon()) >= 0 || reference.getTaxonIndex(n.getTaxon()) >= 0)
            // not a linkage group
            continue;
        int gid = hlm.getTaxonIndex(n.getTaxon()) - reference.getTaxonCount();
        if (gid < 0) {
            System.err.println("big trouble, little china");
        }
        Set<Taxon> group = hlm.getGroup(gid);
        if (group.size() == 0) {
            // remove the group completely
            SimpleNode parent = n.getParent();
            parent.removeChild(n);
            if (parent.getChildCount() == 1) {
                SimpleNode grandparent = parent.getParent();
                SimpleNode child = parent.getChild(0);
                parent.removeChild(child);
                if (grandparent == null) {
                    // parent is root!  other child should become root.
                    root = child;
                } else {
                    grandparent.removeChild(parent);
                    grandparent.addChild(child);
                }
            }
        } else if (group.size() == 1) {
            // swap the group with the constituent read
            Taxon[] tax = new Taxon[group.size()];
            tax = (Taxon[]) group.toArray(tax);
            int rI = getTaxonNode(tax[0], nodes);
            SimpleNode parent = n.getParent();
            parent.removeChild(n);
            parent.addChild(nodes[rI]);
        } else {
            // create a star tree with the reads
            Taxon[] tax = new Taxon[group.size()];
            tax = (Taxon[]) group.toArray(tax);
            SimpleNode parent = n.getParent();
            parent.removeChild(n);
            parent.addChild(nodes[nextFree]);
            int tI = 0;
            for (; tI < tax.length - 2; tI++) {
                int rI = getTaxonNode(tax[tI], nodes);
                nodes[nextFree].addChild(nodes[rI]);
                nodes[nextFree].addChild(nodes[nextFree + 1]);
                nextFree++;
            }
            int rI = getTaxonNode(tax[tI], nodes);
            nodes[nextFree].addChild(nodes[rI]);
            int rJ = getTaxonNode(tax[tI + 1], nodes);
            nodes[nextFree].addChild(nodes[rJ]);
            nextFree++;
        }
    }
    SimpleTree st = new SimpleTree(root);
    return st;
}
Also used : NodeRef(dr.evolution.tree.NodeRef) TaxonList(dr.evolution.util.TaxonList) Taxon(dr.evolution.util.Taxon) SimpleTree(dr.evolution.tree.SimpleTree) SimpleNode(dr.evolution.tree.SimpleNode)

Example 34 with TaxonList

use of dr.evolution.util.TaxonList in project beast-mcmc by beast-dev.

the class OldCoalescentSimulatorParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    CoalescentSimulator simulator = new CoalescentSimulator();
    DemographicModel demoModel = (DemographicModel) xo.getChild(DemographicModel.class);
    List<TaxonList> taxonLists = new ArrayList<TaxonList>();
    List<Tree> subtrees = new ArrayList<Tree>();
    double rootHeight = xo.getAttribute(ROOT_HEIGHT, -1.0);
    if (xo.hasAttribute(RESCALE_HEIGHT)) {
        rootHeight = xo.getDoubleAttribute(RESCALE_HEIGHT);
    }
    // should have one child that is node
    for (int i = 0; i < xo.getChildCount(); i++) {
        final Object child = xo.getChild(i);
        // AER - swapped the order of these round because Trees are TaxonLists...
        if (child instanceof Tree) {
            subtrees.add((Tree) child);
        } else if (child instanceof TaxonList) {
            taxonLists.add((TaxonList) child);
        } else if (xo.getChildName(i).equals(CONSTRAINED_TAXA)) {
            XMLObject constrainedTaxa = (XMLObject) child;
            // all taxa
            final TaxonList taxa = (TaxonList) constrainedTaxa.getChild(TaxonList.class);
            List<CoalescentSimulator.TaxaConstraint> constraints = new ArrayList<CoalescentSimulator.TaxaConstraint>();
            final String setsNotCompatibleMessage = "taxa sets not compatible";
            for (int nc = 0; nc < constrainedTaxa.getChildCount(); ++nc) {
                final Object object = constrainedTaxa.getChild(nc);
                if (object instanceof XMLObject) {
                    final XMLObject constraint = (XMLObject) object;
                    if (constraint.getName().equals(TMRCA_CONSTRAINT)) {
                        TaxonList taxaSubSet = (TaxonList) constraint.getChild(TaxonList.class);
                        ParametricDistributionModel dist = (ParametricDistributionModel) constraint.getChild(ParametricDistributionModel.class);
                        boolean isMono = constraint.getAttribute(IS_MONOPHYLETIC, true);
                        final CoalescentSimulator.TaxaConstraint taxaConstraint = new CoalescentSimulator.TaxaConstraint(taxaSubSet, dist, isMono);
                        int insertPoint;
                        for (insertPoint = 0; insertPoint < constraints.size(); ++insertPoint) {
                            // if new <= constraints[insertPoint] insert before insertPoint
                            final CoalescentSimulator.TaxaConstraint iConstraint = constraints.get(insertPoint);
                            if (iConstraint.isMonophyletic) {
                                if (!taxaConstraint.isMonophyletic) {
                                    continue;
                                }
                                final TaxonList taxonsip = iConstraint.taxons;
                                final int nIn = simulator.sizeOfIntersection(taxonsip, taxaSubSet);
                                if (nIn == taxaSubSet.getTaxonCount()) {
                                    break;
                                }
                                if (nIn > 0 && nIn != taxonsip.getTaxonCount()) {
                                    throw new XMLParseException(setsNotCompatibleMessage);
                                }
                            } else {
                                // reached non mono area
                                if (!taxaConstraint.isMonophyletic) {
                                    if (iConstraint.upper >= taxaConstraint.upper) {
                                        break;
                                    }
                                } else {
                                    break;
                                }
                            }
                        }
                        constraints.add(insertPoint, taxaConstraint);
                    }
                }
            }
            final int nConstraints = constraints.size();
            if (nConstraints == 0) {
                if (taxa != null) {
                    taxonLists.add(taxa);
                }
            } else {
                for (int nc = 0; nc < nConstraints; ++nc) {
                    CoalescentSimulator.TaxaConstraint cnc = constraints.get(nc);
                    if (!cnc.isMonophyletic) {
                        for (int nc1 = nc - 1; nc1 >= 0; --nc1) {
                            CoalescentSimulator.TaxaConstraint cnc1 = constraints.get(nc1);
                            int x = simulator.sizeOfIntersection(cnc.taxons, cnc1.taxons);
                            if (x > 0) {
                                Taxa combinedTaxa = new Taxa(cnc.taxons);
                                combinedTaxa.addTaxa(cnc1.taxons);
                                cnc = new CoalescentSimulator.TaxaConstraint(combinedTaxa, cnc.lower, cnc.upper, cnc.isMonophyletic);
                                constraints.set(nc, cnc);
                            }
                        }
                    }
                }
                // determine upper bound for each set.
                double[] upper = new double[nConstraints];
                for (int nc = nConstraints - 1; nc >= 0; --nc) {
                    final CoalescentSimulator.TaxaConstraint cnc = constraints.get(nc);
                    if (cnc.realLimits()) {
                        upper[nc] = cnc.upper;
                    } else {
                        upper[nc] = Double.POSITIVE_INFINITY;
                    }
                }
                for (int nc = nConstraints - 1; nc >= 0; --nc) {
                    final CoalescentSimulator.TaxaConstraint cnc = constraints.get(nc);
                    if (upper[nc] < Double.POSITIVE_INFINITY) {
                        for (int nc1 = nc - 1; nc1 >= 0; --nc1) {
                            final CoalescentSimulator.TaxaConstraint cnc1 = constraints.get(nc1);
                            if (simulator.contained(cnc1.taxons, cnc.taxons)) {
                                upper[nc1] = Math.min(upper[nc1], upper[nc]);
                                if (cnc1.realLimits() && cnc1.lower > upper[nc1]) {
                                    throw new XMLParseException(setsNotCompatibleMessage);
                                }
                                break;
                            }
                        }
                    }
                }
                // collect subtrees here
                List<Tree> st = new ArrayList<Tree>();
                for (int nc = 0; nc < constraints.size(); ++nc) {
                    final CoalescentSimulator.TaxaConstraint nxt = constraints.get(nc);
                    // collect all previously built subtrees which are a subset of taxa set to be added
                    List<Tree> subs = new ArrayList<Tree>();
                    Taxa newTaxons = new Taxa(nxt.taxons);
                    for (int k = 0; k < st.size(); ++k) {
                        final Tree stk = st.get(k);
                        int x = simulator.sizeOfIntersection(stk, nxt.taxons);
                        if (x == st.get(k).getTaxonCount()) {
                            final Tree tree = st.remove(k);
                            --k;
                            subs.add(tree);
                            newTaxons.removeTaxa(tree);
                        }
                    }
                    SimpleTree tree = simulator.simulateTree(newTaxons, demoModel);
                    final double lower = nxt.realLimits() ? nxt.lower : 0;
                    if (upper[nc] < Double.MAX_VALUE) {
                        simulator.attemptToScaleTree(tree, (lower + upper[nc]) / 2);
                    }
                    if (subs.size() > 0) {
                        if (tree.getTaxonCount() > 0)
                            subs.add(tree);
                        double h = -1;
                        if (upper[nc] < Double.MAX_VALUE) {
                            for (Tree t : subs) {
                                // protect against 1 taxa tree with height 0
                                final double rh = t.getNodeHeight(t.getRoot());
                                h = Math.max(h, rh > 0 ? rh : (lower + upper[nc]) / 2);
                            }
                            h = (h + upper[nc]) / 2;
                        }
                        tree = simulator.simulateTree(subs.toArray(new Tree[subs.size()]), demoModel, h, true);
                    }
                    st.add(tree);
                }
                // add a taxon list for remaining taxa
                if (taxa != null) {
                    final Taxa list = new Taxa();
                    for (int j = 0; j < taxa.getTaxonCount(); ++j) {
                        Taxon taxonj = taxa.getTaxon(j);
                        for (Tree aSt : st) {
                            if (aSt.getTaxonIndex(taxonj) >= 0) {
                                taxonj = null;
                                break;
                            }
                        }
                        if (taxonj != null) {
                            list.addTaxon(taxonj);
                        }
                    }
                    if (list.getTaxonCount() > 0) {
                        taxonLists.add(list);
                    }
                }
                if (st.size() > 1) {
                    final Tree t = simulator.simulateTree(st.toArray(new Tree[st.size()]), demoModel, -1, false);
                    subtrees.add(t);
                } else {
                    subtrees.add(st.get(0));
                }
            }
        }
    }
    if (taxonLists.size() == 0) {
        if (subtrees.size() == 1) {
            return subtrees.get(0);
        }
        throw new XMLParseException("Expected at least one taxonList or two subtrees in " + getParserName() + " element.");
    }
    try {
        Tree[] trees = new Tree[taxonLists.size() + subtrees.size()];
        // simulate each taxonList separately
        for (int i = 0; i < taxonLists.size(); i++) {
            trees[i] = simulator.simulateTree(taxonLists.get(i), demoModel);
        }
        // add the preset trees
        for (int i = 0; i < subtrees.size(); i++) {
            trees[i + taxonLists.size()] = subtrees.get(i);
        }
        return simulator.simulateTree(trees, demoModel, rootHeight, trees.length != 1);
    } catch (IllegalArgumentException iae) {
        throw new XMLParseException(iae.getMessage());
    }
}
Also used : TaxonList(dr.evolution.util.TaxonList) Taxon(dr.evolution.util.Taxon) ArrayList(java.util.ArrayList) DemographicModel(dr.evomodel.coalescent.DemographicModel) SimpleTree(dr.evolution.tree.SimpleTree) CoalescentSimulator(dr.evomodel.coalescent.CoalescentSimulator) Taxa(dr.evolution.util.Taxa) ParametricDistributionModel(dr.inference.distribution.ParametricDistributionModel) Tree(dr.evolution.tree.Tree) SimpleTree(dr.evolution.tree.SimpleTree)

Example 35 with TaxonList

use of dr.evolution.util.TaxonList in project beast-mcmc by beast-dev.

the class CoalescentSimulatorParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    CoalescentSimulator simulator = new CoalescentSimulator();
    DemographicModel demoModel = (DemographicModel) xo.getChild(DemographicModel.class);
    List<TaxonList> taxonLists = new ArrayList<TaxonList>();
    List<Tree> subtrees = new ArrayList<Tree>();
    double height = xo.getAttribute(HEIGHT, Double.NaN);
    // should have one child that is node
    for (int i = 0; i < xo.getChildCount(); i++) {
        final Object child = xo.getChild(i);
        // AER - swapped the order of these round because Trees are TaxonLists...
        if (child instanceof Tree) {
            subtrees.add((Tree) child);
        } else if (child instanceof TaxonList) {
            taxonLists.add((TaxonList) child);
        }
    }
    if (taxonLists.size() == 0) {
        if (subtrees.size() == 1) {
            return subtrees.get(0);
        }
        throw new XMLParseException("Expected at least one taxonList or two subtrees in " + getParserName() + " element.");
    }
    Taxa remainingTaxa = new Taxa();
    for (int i = 0; i < taxonLists.size(); i++) {
        remainingTaxa.addTaxa(taxonLists.get(i));
    }
    for (int i = 0; i < subtrees.size(); i++) {
        remainingTaxa.removeTaxa(subtrees.get(i));
    }
    try {
        Tree[] trees = new Tree[subtrees.size() + remainingTaxa.getTaxonCount()];
        // add the preset trees
        for (int i = 0; i < subtrees.size(); i++) {
            trees[i] = subtrees.get(i);
        }
        // add all the remaining taxa in as single tip trees...
        for (int i = 0; i < remainingTaxa.getTaxonCount(); i++) {
            Taxa tip = new Taxa();
            tip.addTaxon(remainingTaxa.getTaxon(i));
            trees[i + subtrees.size()] = simulator.simulateTree(tip, demoModel);
        }
        return simulator.simulateTree(trees, demoModel, height, trees.length != 1);
    } catch (IllegalArgumentException iae) {
        throw new XMLParseException(iae.getMessage());
    }
}
Also used : TaxonList(dr.evolution.util.TaxonList) ArrayList(java.util.ArrayList) DemographicModel(dr.evomodel.coalescent.DemographicModel) CoalescentSimulator(dr.evomodel.coalescent.CoalescentSimulator) Taxa(dr.evolution.util.Taxa) SimpleTree(dr.evolution.tree.SimpleTree) Tree(dr.evolution.tree.Tree)

Aggregations

TaxonList (dr.evolution.util.TaxonList)59 Tree (dr.evolution.tree.Tree)18 Taxon (dr.evolution.util.Taxon)18 TreeUtils (dr.evolution.tree.TreeUtils)15 Taxa (dr.evolution.util.Taxa)13 Parameter (dr.inference.model.Parameter)13 ArrayList (java.util.ArrayList)13 TreeModel (dr.evomodel.tree.TreeModel)11 Alignment (dr.evolution.alignment.Alignment)6 SitePatterns (dr.evolution.alignment.SitePatterns)5 SimpleTree (dr.evolution.tree.SimpleTree)5 SimpleAlignment (dr.evolution.alignment.SimpleAlignment)4 Importer (dr.evolution.io.Importer)4 ImportException (dr.evolution.io.Importer.ImportException)4 NexusImporter (dr.evolution.io.NexusImporter)4 CoalescentLikelihood (dr.evomodel.coalescent.CoalescentLikelihood)4 CoalescentSimulator (dr.evomodel.coalescent.CoalescentSimulator)4 NodeRef (dr.evolution.tree.NodeRef)3 SubstitutionModel (dr.evomodel.substmodel.SubstitutionModel)3 CompoundLikelihood (dr.inference.model.CompoundLikelihood)3