Examples with CompoundDistribution beast.core.util.CompoundDistribution used on opensource projects

Example 11 with CompoundDistribution

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

the class BeautiBase method assertPriorsEqual.

void assertPriorsEqual(String... ids) {
    System.err.println("assertPriorsEqual");
    CompoundDistribution prior = (CompoundDistribution) doc.pluginmap.get("prior");
    List<Distribution> priors = prior.pDistributions.get();
    for (String id : ids) {
        boolean found = false;
        for (BEASTObject node : priors) {
            if (node.getID().equals(id)) {
                found = true;
            }
        }
        assertThat(found).as("Could not find beastObject with ID " + id).isEqualTo(true);
    }
    List<String> extras = new ArrayList<>();
    for (BEASTObject node : priors) {
        boolean found = false;
        for (String id : ids) {
            if (node.getID().equals(id)) {
                found = true;
            }
        }
        if (!found) {
            extras.add(node.getID());
        }
    }
    if (extras.size() != 0) {
        System.err.println("Extra ids found: " + Arrays.toString(extras.toArray(new String[] {})));
    }
    assertThat(ids.length).as("list of beastObjects do not match").isEqualTo(priors.size());
    ;
}

Example 12 with CompoundDistribution

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

the class BeautiBase method assertParameterCountInPriorIs.

void assertParameterCountInPriorIs(int i) {
    // count nr of parameters in Prior objects in prior
    // including those for prior distributions (Normal, etc)
    // useful to make sure they do (or do not) get linked
    Set<Function> parameters = new LinkedHashSet<>();
    CompoundDistribution prior = (CompoundDistribution) doc.pluginmap.get("prior");
    for (Distribution p : prior.pDistributions.get()) {
        if (p instanceof Prior) {
            Prior p2 = (Prior) p;
            parameters.add(p2.m_x.get());
            for (BEASTInterface o : p2.distInput.get().listActiveBEASTObjects()) {
                if (o instanceof Parameter) {
                    parameters.add((Parameter<?>) o);
                }
            }
        }
    }
    System.err.println("Number of parameters in prior = " + parameters.size());
    if (i >= 0) {
        assertThat(parameters.size()).as("Expected " + i + " parameters in prior").isEqualTo(i);
    }
}

Example 13 with CompoundDistribution

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

the class CalibratedBirthDeathModel 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.initAndValidate();
                            cals.add(cal);
                            taxaSets.add(cal.taxa());
                            cal.taxa().initAndValidate();
                            calCount++;
                            calcCalibrations = false;
                        } else {
                            if (_MRCAPrior.isMonophyleticInput.get()) {
                                Log.warning.println("WARNING: MRCAPriors must have a distribution when monophyletic for Calibrated Yule prior");
                            }
                        }
                    }
                }
            }
        }
        xclades = new int[calCount][];
    }
    if (calCount == 0) {
        // 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;
        }
    }
    isYule = deathToBirthRatioInput.get() == null && sampleProbabilityInput.get() == null;
    userPDF = userMarInput.get();
    if (userPDF == null) {
        boolean needTables = false;
        if (type == Type.OVER_ALL_TOPOS) {
            if (calCount == 1 && isYule) {
            // 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 (isYule && calCount == 2 && orderedCalibrations[1].taxa().containsAll(orderedCalibrations[0].taxa())) {
                // closed form formulas
                } else {
                    needTables = true;
                    lastHeights = new double[calCount];
                }
            }
        } else if (type == Type.OVER_RANKED_COUNTS) {
            // setUpTables(tree.getLeafNodeCount() + 1);
            needTables = true;
        }
        if (needTables) {
            setUpTables(tree.getLeafNodeCount() + 1);
            linsIter = new CalibrationLineagesIterator(this.xclades, this.taxaPartialOrder, maximal, tree.getLeafNodeCount());
        }
    }
    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 Birth-Death Model does not handle dated tips correctly. " + "Consider using a coalescent prior instead.");
            break;
        }
    }
}

Example 14 with CompoundDistribution

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

the class BeautiDoc method addMRCAPrior.

public void addMRCAPrior(MRCAPrior mrcaPrior) {
    Tree tree = (Tree) pluginmap.get("Tree.t:" + alignments.get(0).getID());
    if (tree == null) {
        for (String key : pluginmap.keySet()) {
            if (key.startsWith("Tree.t:")) {
                tree = (Tree) pluginmap.get(key);
                break;
            }
        }
    }
    // make sure we have the appropriate tree:
    if (alignments.size() > 1) {
        String testTaxon = mrcaPrior.taxonsetInput.get().toString().split("\n")[1].trim();
        // = tree.getTaxaNames();
        String[] taxaNames;
        int index = -1;
        int j = 0;
        while (index < 0 && j++ < alignments.size()) {
            tree = (Tree) pluginmap.get("Tree.t:" + alignments.get(j - 1).getID());
            taxaNames = tree.getTaxaNames();
            for (int i = 0; i < taxaNames.length && index < 0; i++) if (testTaxon.equals(taxaNames[i]))
                index = i;
        }
    }
    CompoundDistribution prior = (CompoundDistribution) pluginmap.get("prior");
    mrcaPrior.treeInput.setValue(tree, mrcaPrior);
    ParametricDistribution distr = mrcaPrior.distInput.get();
    TaxonSet t = mrcaPrior.taxonsetInput.get();
    if (taxaset.keySet().contains(t.getID())) {
        Log.warning.println("taxonset " + t.getID() + " already exists: MRCAPrior " + mrcaPrior.getID() + " can not be added");
    } else {
        taxaset.put(t.getID(), t);
        // ensure TaxonSets are not duplicated
        List<Taxon> taxa = t.taxonsetInput.get();
        for (int i = 0; i < taxa.size(); i++) {
            if (taxaset.containsKey(taxa.get(i).getID())) {
                taxa.set(i, taxaset.get(taxa.get(i).getID()));
            } else {
                taxaset.put(taxa.get(i).getID(), taxa.get(i));
            }
        }
        if (distr instanceof Normal && (Double.isInfinite(((Normal) distr).sigmaInput.get().getValue()))) {
        // it is a 'fixed' calibration, no need to add a distribution
        } else {
            prior.pDistributions.setValue(mrcaPrior, prior);
            connect(mrcaPrior, "tracelog", "log");
        }
    }
    if (t.taxonsetInput.get().size() == 1 && distr != null) {
        // only add operators if it is NOT a 'fixed' calibration
        if (!(distr instanceof Normal && (Double.isInfinite(((Normal) distr).sigmaInput.get().getValue())))) {
            TipDatesRandomWalker operator = new TipDatesRandomWalker();
            t.initAndValidate();
            operator.initByName("taxonset", t, "weight", 1.0, "tree", tree, "windowSize", 1.0);
            operator.setID("TipDatesRandomWalker." + t.getID());
            MCMC mcmc = (MCMC) this.mcmc.get();
            mcmc.operatorsInput.setValue(operator, mcmc);
        }
        // set up date trait
        double date = distr.getMean();
        TraitSet dateTrait = null;
        for (TraitSet ts : tree.m_traitList.get()) {
            if (ts.isDateTrait()) {
                dateTrait = ts;
            }
        }
        if (dateTrait == null) {
            dateTrait = new TraitSet();
            dateTrait.setID("traitsetDate");
            registerPlugin(dateTrait);
            dateTrait.initByName("traitname", TraitSet.DATE_BACKWARD_TRAIT, "taxa", tree.getTaxonset(), "value", t.taxonsetInput.get().get(0).getID() + "=" + date);
            tree.m_traitList.setValue(dateTrait, tree);
            tree.initAndValidate();
        } else {
            dateTrait.traitsInput.setValue(dateTrait.traitsInput.get() + ",\n" + t.taxonsetInput.get().get(0).getID() + "=" + date, dateTrait);
        }
        mrcaPrior.onlyUseTipsInput.setValue(true, mrcaPrior);
    }
}

Example 15 with CompoundDistribution

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

the class BeautiDoc method determinePartitions.

public void determinePartitions() {
    CompoundDistribution likelihood = (CompoundDistribution) pluginmap.get("likelihood");
    if (likelihood == null) {
        return;
    }
    partitionNames.clear();
    possibleContexts.clear();
    for (Distribution distr : likelihood.pDistributions.get()) {
        if (distr instanceof GenericTreeLikelihood) {
            GenericTreeLikelihood treeLikelihood = (GenericTreeLikelihood) distr;
            alignments.add(treeLikelihood.dataInput.get());
            PartitionContext context = new PartitionContext(treeLikelihood);
            partitionNames.add(context);
            boolean found = false;
            for (PartitionContext context2 : possibleContexts) {
                if (context.equals(context2)) {
                    found = true;
                }
            }
            if (!found) {
                possibleContexts.add(context);
            }
        }
    }
    alignments.clear();
    for (int i = 0; i < 3; i++) {
        pPartitionByAlignments[i].clear();
        pPartition[i].clear();
        currentPartitions[i].clear();
    }
    List<GenericTreeLikelihood> treeLikelihoods = new ArrayList<>();
    for (Distribution distr : likelihood.pDistributions.get()) {
        if (distr instanceof GenericTreeLikelihood) {
            GenericTreeLikelihood treeLikelihood = (GenericTreeLikelihood) distr;
            alignments.add(treeLikelihood.dataInput.get());
            treeLikelihoods.add(treeLikelihood);
        }
    }
    for (Distribution distr : likelihood.pDistributions.get()) {
        if (distr instanceof GenericTreeLikelihood) {
            GenericTreeLikelihood treeLikelihood = (GenericTreeLikelihood) distr;
            try {
                // sync SiteModel, ClockModel and Tree to any changes that
                // may have occurred
                // this should only affect the clock model in practice
                int partition = getPartitionNr((BEASTInterface) treeLikelihood.siteModelInput.get());
                GenericTreeLikelihood treeLikelihood2 = treeLikelihoods.get(partition);
                treeLikelihood.siteModelInput.setValue(treeLikelihood2.siteModelInput.get(), treeLikelihood);
                currentPartitions[0].add(partition);
                BranchRateModel rateModel = treeLikelihood.branchRateModelInput.get();
                if (rateModel != null) {
                    partition = getPartitionNr((BEASTInterface) rateModel);
                    treeLikelihood2 = treeLikelihoods.get(partition);
                    treeLikelihood.branchRateModelInput.setValue(treeLikelihood2.branchRateModelInput.get(), treeLikelihood);
                    currentPartitions[1].add(partition);
                } else {
                    currentPartitions[1].add(0);
                }
                partition = getPartitionNr((BEASTInterface) treeLikelihood.treeInput.get());
                treeLikelihood2 = treeLikelihoods.get(partition);
                treeLikelihood.treeInput.setValue(treeLikelihood2.treeInput.get(), treeLikelihood);
                currentPartitions[2].add(partition);
            } catch (Exception e) {
                e.printStackTrace();
            }
            pPartitionByAlignments[0].add(treeLikelihood);
            pPartitionByAlignments[1].add(treeLikelihood);
            pPartitionByAlignments[2].add(treeLikelihood);
        }
    }
    int partitionCount = partitionNames.size();
    for (int i = 0; i < 3; i++) {
        boolean[] usedPartition = new boolean[partitionCount];
        for (int j = 0; j < partitionCount; j++) {
            // getPartitionNr(m_pPartitionByAlignments[i].get(j));
            int partitionIndex = currentPartitions[i].get(j);
            usedPartition[partitionIndex] = true;
        }
        for (int j = 0; j < partitionCount; j++) {
            if (usedPartition[j]) {
                pPartition[i].add(pPartitionByAlignments[i].get(j));
            }
        }
    }
    Log.warning.println("PARTITIONS0:\n");
    Log.warning.println(Arrays.toString(currentPartitions));
}