Search in sources :

Example 6 with CompoundDistribution

use of beast.core.util.CompoundDistribution 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 CompoundDistribution

use of beast.core.util.CompoundDistribution in project beast2 by CompEvol.

the class AlignmentListInputEditor method initTableData.

void initTableData() {
    this.likelihoods = new GenericTreeLikelihood[partitionCount];
    if (tableData == null) {
        tableData = new Object[partitionCount][NR_OF_COLUMNS];
    }
    CompoundDistribution likelihoods = (CompoundDistribution) doc.pluginmap.get("likelihood");
    for (int i = 0; i < partitionCount; i++) {
        Alignment data = alignments.get(i);
        // partition name
        tableData[i][NAME_COLUMN] = data;
        // alignment name
        if (data instanceof FilteredAlignment) {
            tableData[i][FILE_COLUMN] = ((FilteredAlignment) data).alignmentInput.get();
        } else {
            tableData[i][FILE_COLUMN] = data;
        }
        // # taxa
        tableData[i][TAXA_COLUMN] = data.getTaxonCount();
        // # sites
        tableData[i][SITES_COLUMN] = data.getSiteCount();
        // Data type
        tableData[i][TYPE_COLUMN] = data.getDataType();
        // site model
        GenericTreeLikelihood likelihood = (GenericTreeLikelihood) likelihoods.pDistributions.get().get(i);
        assert (likelihood != null);
        this.likelihoods[i] = likelihood;
        tableData[i][SITEMODEL_COLUMN] = getPartition(likelihood.siteModelInput);
        // clock model
        tableData[i][CLOCKMODEL_COLUMN] = getPartition(likelihood.branchRateModelInput);
        // tree
        tableData[i][TREE_COLUMN] = getPartition(likelihood.treeInput);
        // useAmbiguities
        tableData[i][USE_AMBIGUITIES_COLUMN] = null;
        try {
            if (hasUseAmbiguitiesInput(i)) {
                tableData[i][USE_AMBIGUITIES_COLUMN] = likelihood.getInputValue("useAmbiguities");
            }
        } catch (Exception e) {
        // ignore
        }
    }
}
Also used : CompoundDistribution(beast.core.util.CompoundDistribution) FilteredAlignment(beast.evolution.alignment.FilteredAlignment) Alignment(beast.evolution.alignment.Alignment) GenericTreeLikelihood(beast.evolution.likelihood.GenericTreeLikelihood) FilteredAlignment(beast.evolution.alignment.FilteredAlignment)

Example 8 with CompoundDistribution

use of beast.core.util.CompoundDistribution in project beast2 by CompEvol.

the class SiteModelInputEditor method customConnector.

public static boolean customConnector(BeautiDoc doc) {
    try {
        DeltaExchangeOperator operator = (DeltaExchangeOperator) doc.pluginmap.get("FixMeanMutationRatesOperator");
        if (operator == null) {
            return false;
        }
        List<RealParameter> parameters = operator.parameterInput.get();
        parameters.clear();
        // String weights = "";
        CompoundDistribution likelihood = (CompoundDistribution) doc.pluginmap.get("likelihood");
        boolean hasOneEstimatedRate = false;
        List<String> rateIDs = new ArrayList<>();
        List<Integer> weights = new ArrayList<>();
        for (Distribution d : likelihood.pDistributions.get()) {
            GenericTreeLikelihood treelikelihood = (GenericTreeLikelihood) d;
            Alignment data = treelikelihood.dataInput.get();
            int weight = data.getSiteCount();
            if (data.isAscertained) {
                weight -= data.getExcludedPatternCount();
            }
            if (treelikelihood.siteModelInput.get() instanceof SiteModel) {
                SiteModel siteModel = (SiteModel) treelikelihood.siteModelInput.get();
                RealParameter mutationRate = siteModel.muParameterInput.get();
                // clockRate.m_bIsEstimated.setValue(true, clockRate);
                if (mutationRate.isEstimatedInput.get()) {
                    hasOneEstimatedRate = true;
                    if (rateIDs.indexOf(mutationRate.getID()) == -1) {
                        parameters.add(mutationRate);
                        weights.add(weight);
                        rateIDs.add(mutationRate.getID());
                    } else {
                        int k = rateIDs.indexOf(mutationRate.getID());
                        weights.set(k, weights.get(k) + weight);
                    }
                }
            }
        }
        IntegerParameter weightParameter;
        if (weights.size() == 0) {
            weightParameter = new IntegerParameter();
        } else {
            String weightString = "";
            for (int k : weights) {
                weightString += k + " ";
            }
            weightParameter = new IntegerParameter(weightString);
            weightParameter.setID("weightparameter");
        }
        weightParameter.isEstimatedInput.setValue(false, weightParameter);
        operator.parameterWeightsInput.setValue(weightParameter, operator);
        return hasOneEstimatedRate;
    } catch (Exception e) {
    }
    return false;
}
Also used : IntegerParameter(beast.core.parameter.IntegerParameter) GenericTreeLikelihood(beast.evolution.likelihood.GenericTreeLikelihood) ArrayList(java.util.ArrayList) RealParameter(beast.core.parameter.RealParameter) SiteModel(beast.evolution.sitemodel.SiteModel) InvocationTargetException(java.lang.reflect.InvocationTargetException) CompoundDistribution(beast.core.util.CompoundDistribution) Alignment(beast.evolution.alignment.Alignment) CompoundDistribution(beast.core.util.CompoundDistribution) DeltaExchangeOperator(beast.evolution.operators.DeltaExchangeOperator)

Example 9 with CompoundDistribution

use of beast.core.util.CompoundDistribution 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)

Example 10 with CompoundDistribution

use of beast.core.util.CompoundDistribution in project beast2 by CompEvol.

the class BeautiPanel method cloneFrom.

/**
 * Clones partition identified by sourceID to targetID and type (Site/Clock/Tree model)
 * as stored in config.
 * @param sourceID
 * @param targetID
 */
public void cloneFrom(String sourceID, String targetID) {
    if (sourceID.equals(targetID)) {
        return;
    }
    String type = config.hasPartitionsInput.get().toString();
    java.util.List<BEASTInterface> list = doc.getPartitions(type);
    int source = -1, target = -1;
    for (int i = 0; i < list.size(); i++) {
        BEASTInterface partition = list.get(i);
        if (type.equals("SiteModel")) {
            partition = (BEASTInterface) ((GenericTreeLikelihood) partition).siteModelInput.get();
        } else if (type.equals("ClockModel")) {
            partition = ((GenericTreeLikelihood) partition).branchRateModelInput.get();
        } else if (type.equals("Tree")) {
            partition = (BEASTInterface) ((GenericTreeLikelihood) partition).treeInput.get();
        }
        String partitionID = partition.getID();
        partitionID = partitionID.substring(partitionID.lastIndexOf('.') + 1);
        if (partitionID.length() > 1 && partitionID.charAt(1) == ':') {
            partitionID = partitionID.substring(2);
        }
        if (partitionID.equals(sourceID)) {
            source = i;
        }
        if (partitionID.equals(targetID)) {
            target = i;
        }
    }
    if (target == -1) {
        throw new RuntimeException("Programmer error: sourceID and targetID should be in list");
    }
    CompoundDistribution likelihoods = (CompoundDistribution) doc.pluginmap.get("likelihood");
    GenericTreeLikelihood likelihoodSource = (GenericTreeLikelihood) likelihoods.pDistributions.get().get(source);
    GenericTreeLikelihood likelihood = (GenericTreeLikelihood) likelihoods.pDistributions.get().get(target);
    PartitionContext oldContext = doc.getContextFor(likelihoodSource);
    PartitionContext newContext = doc.getContextFor(likelihood);
    // this ensures the config.sync does not set any input value
    config._input.setValue(null, config);
    if (type.equals("SiteModel")) {
        SiteModelInterface siteModelSource = likelihoodSource.siteModelInput.get();
        SiteModelInterface siteModel = null;
        try {
            siteModel = (SiteModel.Base) BeautiDoc.deepCopyPlugin((BEASTInterface) siteModelSource, likelihood, (MCMC) doc.mcmc.get(), oldContext, newContext, doc, null);
        } catch (RuntimeException e) {
            JOptionPane.showMessageDialog(this, "Could not clone " + sourceID + " to " + targetID + " " + e.getMessage());
            return;
        }
        likelihood.siteModelInput.setValue(siteModel, likelihood);
        return;
    } else if (type.equals("ClockModel")) {
        BranchRateModel clockModelSource = likelihoodSource.branchRateModelInput.get();
        BranchRateModel clockModel = null;
        try {
            clockModel = (BranchRateModel) BeautiDoc.deepCopyPlugin((BEASTInterface) clockModelSource, likelihood, (MCMC) doc.mcmc.get(), oldContext, newContext, doc, null);
        } catch (Exception e) {
            JOptionPane.showMessageDialog(this, "Could not clone " + sourceID + " to " + targetID + " " + e.getMessage());
            return;
        }
        // make sure that *if* the clock model has a tree as input, it is
        // the same as for the likelihood
        TreeInterface tree = null;
        try {
            for (Input<?> input : ((BEASTInterface) clockModel).listInputs()) {
                if (input.getName().equals("tree")) {
                    tree = (TreeInterface) input.get();
                }
            }
        } catch (IllegalArgumentException e) {
            // TODO Auto-generated catch block
            e.printStackTrace();
        }
        if (tree != null && tree != likelihood.treeInput.get()) {
            // likelihood.treeInput.setValue(tree, likelihood);
            JOptionPane.showMessageDialog(null, "Cannot clone clock model with different trees");
            return;
        }
        likelihood.branchRateModelInput.setValue(clockModel, likelihood);
        return;
    } else if (type.equals("Tree")) {
        TreeInterface tree = null;
        TreeInterface treeSource = likelihoodSource.treeInput.get();
        try {
            tree = (TreeInterface) BeautiDoc.deepCopyPlugin((BEASTInterface) treeSource, likelihood, (MCMC) doc.mcmc.get(), oldContext, newContext, doc, null);
        } catch (Exception e) {
            JOptionPane.showMessageDialog(this, "Could not clone " + sourceID + " to " + targetID + " " + e.getMessage());
            return;
        }
        // sanity check: make sure taxon sets are compatible
        Taxon.assertSameTaxa(tree.getID(), tree.getTaxonset().getTaxaNames(), likelihood.dataInput.get().getID(), likelihood.dataInput.get().getTaxaNames());
        likelihood.treeInput.setValue(tree, likelihood);
        return;
    } else {
        throw new RuntimeException("Programmer error calling cloneFrom: Should only clone Site/Clock/Tree model");
    }
}
Also used : GenericTreeLikelihood(beast.evolution.likelihood.GenericTreeLikelihood) MCMC(beast.core.MCMC) SiteModel(beast.evolution.sitemodel.SiteModel) InvocationTargetException(java.lang.reflect.InvocationTargetException) TreeInterface(beast.evolution.tree.TreeInterface) CompoundDistribution(beast.core.util.CompoundDistribution) Input(beast.core.Input) BranchRateModel(beast.evolution.branchratemodel.BranchRateModel) BEASTInterface(beast.core.BEASTInterface) SiteModelInterface(beast.evolution.sitemodel.SiteModelInterface)

Aggregations

CompoundDistribution (beast.core.util.CompoundDistribution)19 Distribution (beast.core.Distribution)12 GenericTreeLikelihood (beast.evolution.likelihood.GenericTreeLikelihood)9 BEASTInterface (beast.core.BEASTInterface)8 ParametricDistribution (beast.math.distributions.ParametricDistribution)8 ArrayList (java.util.ArrayList)7 Tree (beast.evolution.tree.Tree)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 RealParameter (beast.core.parameter.RealParameter)4 TaxonSet (beast.evolution.alignment.TaxonSet)4 BranchRateModel (beast.evolution.branchratemodel.BranchRateModel)4 BEASTObject (beast.core.BEASTObject)3 MCMC (beast.core.MCMC)3 Alignment (beast.evolution.alignment.Alignment)3 SiteModel (beast.evolution.sitemodel.SiteModel)3