Search in sources :

Example 6 with Distribution

use of beast.core.Distribution in project beast2 by CompEvol.

the class CalibratedYuleModel method initAndValidate.

@Override
public void initAndValidate() {
    super.initAndValidate();
    type = correctionTypeInput.get();
    final TreeInterface tree = treeInput.get();
    // shallow copy. we shall change cals later
    final List<CalibrationPoint> cals = new ArrayList<>(calibrationsInput.get());
    int calCount = cals.size();
    final List<TaxonSet> taxaSets = new ArrayList<>(calCount);
    if (cals.size() > 0) {
        xclades = new int[calCount][];
        // convenience
        for (final CalibrationPoint cal : cals) {
            taxaSets.add(cal.taxa());
        }
    } else {
        // find calibration points from prior
        for (final Object beastObject : getOutputs()) {
            if (beastObject instanceof CompoundDistribution) {
                final CompoundDistribution prior = (CompoundDistribution) beastObject;
                for (final Distribution distr : prior.pDistributions.get()) {
                    if (distr instanceof MRCAPrior) {
                        final MRCAPrior _MRCAPrior = (MRCAPrior) distr;
                        // make sure MRCAPrior is monophyletic
                        if (_MRCAPrior.distInput.get() != null) {
                            // make sure MRCAPrior is monophyletic
                            if (!_MRCAPrior.isMonophyleticInput.get()) {
                                throw new IllegalArgumentException("MRCAPriors must be monophyletic for Calibrated Yule prior");
                            }
                            // create CalibrationPoint from MRCAPrior
                            final CalibrationPoint cal = new CalibrationPoint();
                            cal.distInput.setValue(_MRCAPrior.distInput.get(), cal);
                            cal.taxonsetInput.setValue(_MRCAPrior.taxonsetInput.get(), cal);
                            cal.forParentInput.setValue(_MRCAPrior.useOriginateInput.get(), cal);
                            cal.initAndValidate();
                            cals.add(cal);
                            taxaSets.add(cal.taxa());
                            cal.taxa().initAndValidate();
                            calCount++;
                            calcCalibrations = false;
                        } else {
                            if (_MRCAPrior.isMonophyleticInput.get()) {
                                Log.warning.println("WARNING: MRCAPriors (" + _MRCAPrior.getID() + ") must have a distribution when monophyletic. Ignored for Calibrated Yule prior");
                            } else {
                                Log.warning.println("WARNING: MRCAPriors (" + _MRCAPrior.getID() + ") found that is not monophyletic. Ignored for Calibrated Yule prior");
                            }
                        }
                    }
                }
            }
        }
        xclades = new int[calCount][];
    }
    if (calCount == 0) {
        Log.warning.println("WARNING: Calibrated Yule prior could not find any properly configured calibrations. Expect this to crash in a BEAST run.");
        // assume we are in beauti, back off for now
        return;
    }
    for (int k = 0; k < calCount; ++k) {
        final TaxonSet tk = taxaSets.get(k);
        for (int i = k + 1; i < calCount; ++i) {
            final TaxonSet ti = taxaSets.get(i);
            if (ti.containsAny(tk)) {
                if (!(ti.containsAll(tk) || tk.containsAll(ti))) {
                    throw new IllegalArgumentException("Overlapping taxaSets??");
                }
            }
        }
    }
    orderedCalibrations = new CalibrationPoint[calCount];
    {
        int loc = taxaSets.size() - 1;
        while (loc >= 0) {
            assert loc == taxaSets.size() - 1;
            // place maximal taxaSets at end one at a time
            int k = 0;
            for (; /**/
            k < taxaSets.size(); ++k) {
                if (isMaximal(taxaSets, k)) {
                    break;
                }
            }
            final List<String> tk = taxaSets.get(k).asStringList();
            final int tkcount = tk.size();
            this.xclades[loc] = new int[tkcount];
            for (int nt = 0; nt < tkcount; ++nt) {
                final int taxonIndex = getTaxonIndex(tree, tk.get(nt));
                this.xclades[loc][nt] = taxonIndex;
                if (taxonIndex < 0) {
                    throw new IllegalArgumentException("Taxon not found in tree: " + tk.get(nt));
                }
            }
            orderedCalibrations[loc] = cals.remove(k);
            taxaSets.remove(k);
            // cals and taxaSets should match
            --loc;
        }
    }
    // tio[i] will contain all taxaSets contained in the i'th clade, in the form of thier index into orderedCalibrations
    @SuppressWarnings("unchecked") final List<Integer>[] tio = new List[orderedCalibrations.length];
    for (int k = 0; k < orderedCalibrations.length; ++k) {
        tio[k] = new ArrayList<>();
    }
    for (int k = 0; k < orderedCalibrations.length; ++k) {
        final TaxonSet txk = orderedCalibrations[k].taxa();
        for (int i = k + 1; i < orderedCalibrations.length; ++i) {
            if (orderedCalibrations[i].taxa().containsAll(txk)) {
                tio[i].add(k);
                break;
            }
        }
    }
    this.taxaPartialOrder = new int[orderedCalibrations.length][];
    for (int k = 0; k < orderedCalibrations.length; ++k) {
        final List<Integer> tiok = tio[k];
        this.taxaPartialOrder[k] = new int[tiok.size()];
        for (int j = 0; j < tiok.size(); ++j) {
            this.taxaPartialOrder[k][j] = tiok.get(j);
        }
    }
    // true if clade is not contained in any other clade
    final boolean[] maximal = new boolean[calCount];
    for (int k = 0; k < calCount; ++k) {
        maximal[k] = true;
    }
    for (int k = 0; k < calCount; ++k) {
        for (final int i : this.taxaPartialOrder[k]) {
            maximal[i] = false;
        }
    }
    userPDF = userMarInput.get();
    if (userPDF == null) {
        if (type == Type.OVER_ALL_TOPOS) {
            if (calCount == 1) {
            // closed form formula
            } else {
                boolean anyParent = false;
                for (final CalibrationPoint c : orderedCalibrations) {
                    if (c.forParentInput.get()) {
                        anyParent = true;
                    }
                }
                if (anyParent) {
                    throw new IllegalArgumentException("Sorry, not implemented: calibration on parent for more than one clade.");
                }
                if (calCount == 2 && orderedCalibrations[1].taxa().containsAll(orderedCalibrations[0].taxa())) {
                // closed form formulas
                } else {
                    setUpTables(tree.getLeafNodeCount() + 1);
                    linsIter = new CalibrationLineagesIterator(this.xclades, this.taxaPartialOrder, maximal, tree.getLeafNodeCount());
                    lastHeights = new double[calCount];
                }
            }
        } else if (type == Type.OVER_RANKED_COUNTS) {
            setUpTables(tree.getLeafNodeCount() + 1);
        }
    }
    final List<Node> leafs = tree.getExternalNodes();
    final double height = leafs.get(0).getHeight();
    for (final Node leaf : leafs) {
        if (Math.abs(leaf.getHeight() - height) > 1e-8) {
            Log.warning.println("WARNING: Calibrated Yule Model cannot handle dated tips. Use for example a coalescent prior instead.");
            break;
        }
    }
}
Also used : Node(beast.evolution.tree.Node) ArrayList(java.util.ArrayList) TaxonSet(beast.evolution.alignment.TaxonSet) TreeInterface(beast.evolution.tree.TreeInterface) CompoundDistribution(beast.core.util.CompoundDistribution) ParametricDistribution(beast.math.distributions.ParametricDistribution) CompoundDistribution(beast.core.util.CompoundDistribution) Distribution(beast.core.Distribution) MRCAPrior(beast.math.distributions.MRCAPrior) ArrayList(java.util.ArrayList) List(java.util.List)

Example 7 with Distribution

use of beast.core.Distribution in project beast2 by CompEvol.

the class CompoundDistribution method getConditions.

@Override
public List<String> getConditions() {
    List<String> conditions = new ArrayList<>();
    for (Distribution distribution : pDistributions.get()) {
        conditions.addAll(distribution.getConditions());
    }
    conditions.removeAll(getArguments());
    return conditions;
}
Also used : Distribution(beast.core.Distribution) ArrayList(java.util.ArrayList)

Example 8 with Distribution

use of beast.core.Distribution in project beast2 by CompEvol.

the class CompoundDistribution method calculateLogP.

/**
 * Distribution implementation follows *
 */
@Override
public double calculateLogP() {
    logP = 0;
    if (ignore) {
        return logP;
    }
    int workAvailable = 0;
    if (useThreads) {
        for (Distribution dists : pDistributions.get()) {
            if (dists.isDirtyCalculation()) {
                workAvailable++;
            }
        }
    }
    if (useThreads && workAvailable > 1) {
        logP = calculateLogPUsingThreads();
    } else {
        for (Distribution dists : pDistributions.get()) {
            if (dists.isDirtyCalculation()) {
                logP += dists.calculateLogP();
            } else {
                logP += dists.getCurrentLogP();
            }
            if (Double.isInfinite(logP) || Double.isNaN(logP)) {
                return logP;
            }
        }
    }
    return logP;
}
Also used : Distribution(beast.core.Distribution)

Example 9 with Distribution

use of beast.core.Distribution in project beast2 by CompEvol.

the class CompoundDistribution method calculateLogPUsingThreads.

private double calculateLogPUsingThreads() {
    try {
        int dirtyDistrs = 0;
        for (Distribution dists : pDistributions.get()) {
            if (dists.isDirtyCalculation()) {
                dirtyDistrs++;
            }
        }
        countDown = new CountDownLatch(dirtyDistrs);
        // kick off the threads
        for (Distribution dists : pDistributions.get()) {
            if (dists.isDirtyCalculation()) {
                CoreRunnable coreRunnable = new CoreRunnable(dists);
                exec.execute(coreRunnable);
            }
        }
        countDown.await();
        logP = 0;
        for (Distribution distr : pDistributions.get()) {
            logP += distr.getCurrentLogP();
        }
        return logP;
    } catch (RejectedExecutionException | InterruptedException e) {
        useThreads = false;
        Log.err.println("Stop using threads: " + e.getMessage());
        return calculateLogP();
    }
}
Also used : Distribution(beast.core.Distribution) CountDownLatch(java.util.concurrent.CountDownLatch) RejectedExecutionException(java.util.concurrent.RejectedExecutionException)

Example 10 with Distribution

use of beast.core.Distribution in project beast2 by CompEvol.

the class MRCAPriorInputEditor method customConnector.

public static void customConnector(BeautiDoc doc) {
    Object o0 = doc.pluginmap.get("prior");
    if (o0 != null && o0 instanceof CompoundDistribution) {
        CompoundDistribution p = (CompoundDistribution) o0;
        for (Distribution p0 : p.pDistributions.get()) {
            if (p0 instanceof MRCAPrior) {
                MRCAPrior prior = (MRCAPrior) p0;
                if (prior.treeInput.get() != null) {
                    boolean isInState = false;
                    for (BEASTInterface o : prior.treeInput.get().getOutputs()) {
                        if (o instanceof State) {
                            isInState = true;
                            break;
                        }
                    }
                    if (!isInState) {
                        doc.disconnect(prior, "prior", "distribution");
                        doc.disconnect(prior, "tracelog", "log");
                        if (prior.onlyUseTipsInput.get()) {
                            disableTipSampling(prior, doc);
                        }
                        doc.unregisterPlugin(prior);
                        return;
                    }
                }
            }
        }
    }
}
Also used : CompoundDistribution(beast.core.util.CompoundDistribution) State(beast.core.State) CompoundDistribution(beast.core.util.CompoundDistribution) Distribution(beast.core.Distribution) MRCAPrior(beast.math.distributions.MRCAPrior) BEASTInterface(beast.core.BEASTInterface)

Aggregations

Distribution (beast.core.Distribution)17 CompoundDistribution (beast.core.util.CompoundDistribution)12 BEASTInterface (beast.core.BEASTInterface)9 ParametricDistribution (beast.math.distributions.ParametricDistribution)8 ArrayList (java.util.ArrayList)7 GenericTreeLikelihood (beast.evolution.likelihood.GenericTreeLikelihood)5 MRCAPrior (beast.math.distributions.MRCAPrior)5 XMLParserException (beast.util.XMLParserException)5 IOException (java.io.IOException)5 ParserConfigurationException (javax.xml.parsers.ParserConfigurationException)5 TransformerException (javax.xml.transform.TransformerException)5 SAXException (org.xml.sax.SAXException)5 BranchRateModel (beast.evolution.branchratemodel.BranchRateModel)4 Tree (beast.evolution.tree.Tree)4 List (java.util.List)4 BEASTObject (beast.core.BEASTObject)3 Parameter (beast.core.parameter.Parameter)3 TaxonSet (beast.evolution.alignment.TaxonSet)3 StateNode (beast.core.StateNode)2 RealParameter (beast.core.parameter.RealParameter)2