Examples with ParametricDistribution beast.math.distributions.ParametricDistribution used on opensource projects

Example 1 with ParametricDistribution

use of beast.math.distributions.ParametricDistribution in project beast2 by CompEvol.

the class RandomTree method initStateNodes.

// taxonset intersection test
// private boolean intersects(final BitSet bitSet, final BitSet bitSet2) {
// for (int k = bitSet.nextSetBit(0); k >= 0; k = bitSet.nextSetBit(k + 1)) {
// if (bitSet2.get(k)) {
// return true;
// }
// }
// return false;
// }
// returns true if bitSet is a subset of bitSet2
// private boolean isSubset(final BitSet bitSet, final BitSet bitSet2) {
// boolean isSubset = true;
// for (int k = bitSet.nextSetBit(0); isSubset && k >= 0; k = bitSet.nextSetBit(k + 1)) {
// isSubset = bitSet2.get(k);
// }
// return isSubset;
// }
@SuppressWarnings("unchecked")
@Override
public void initStateNodes() {
    // find taxon sets we are dealing with
    taxonSets = new ArrayList<>();
    m_bounds = new ArrayList<>();
    distributions = new ArrayList<>();
    taxonSetIDs = new ArrayList<>();
    lastMonophyletic = 0;
    if (taxaInput.get() != null) {
        taxa.addAll(taxaInput.get().getTaxaNames());
    } else {
        taxa.addAll(m_taxonset.get().asStringList());
    }
    // pick up constraints from outputs, m_inititial input tree and output tree, if any
    List<MRCAPrior> calibrations = new ArrayList<>();
    calibrations.addAll(calibrationsInput.get());
    // pick up constraints in m_initial tree
    for (final Object beastObject : getOutputs()) {
        if (beastObject instanceof MRCAPrior && !calibrations.contains(beastObject)) {
            calibrations.add((MRCAPrior) beastObject);
        }
    }
    if (m_initial.get() != null) {
        for (final Object beastObject : m_initial.get().getOutputs()) {
            if (beastObject instanceof MRCAPrior && !calibrations.contains(beastObject)) {
                calibrations.add((MRCAPrior) beastObject);
            }
        }
    }
    for (final MRCAPrior prior : calibrations) {
        final TaxonSet taxonSet = prior.taxonsetInput.get();
        if (taxonSet != null && !prior.onlyUseTipsInput.get()) {
            final Set<String> usedTaxa = new LinkedHashSet<>();
            if (taxonSet.asStringList() == null) {
                taxonSet.initAndValidate();
            }
            for (final String taxonID : taxonSet.asStringList()) {
                if (!taxa.contains(taxonID)) {
                    throw new IllegalArgumentException("Taxon <" + taxonID + "> could not be found in list of taxa. Choose one of " + taxa);
                }
                usedTaxa.add(taxonID);
            }
            final ParametricDistribution distr = prior.distInput.get();
            final Bound bounds = new Bound();
            if (distr != null) {
                List<BEASTInterface> beastObjects = new ArrayList<>();
                distr.getPredecessors(beastObjects);
                for (int i = beastObjects.size() - 1; i >= 0; i--) {
                    beastObjects.get(i).initAndValidate();
                }
                try {
                    bounds.lower = distr.inverseCumulativeProbability(0.0) + distr.offsetInput.get();
                    bounds.upper = distr.inverseCumulativeProbability(1.0) + distr.offsetInput.get();
                } catch (MathException e) {
                    Log.warning.println("At RandomTree::initStateNodes, bound on MRCAPrior could not be set " + e.getMessage());
                }
            }
            if (prior.isMonophyleticInput.get()) {
                // add any monophyletic constraint
                taxonSets.add(lastMonophyletic, usedTaxa);
                distributions.add(lastMonophyletic, distr);
                m_bounds.add(lastMonophyletic, bounds);
                taxonSetIDs.add(prior.getID());
                lastMonophyletic++;
            } else {
                // only calibrations with finite bounds are added
                if (!Double.isInfinite(bounds.lower) || !Double.isInfinite(bounds.upper)) {
                    taxonSets.add(usedTaxa);
                    distributions.add(distr);
                    m_bounds.add(bounds);
                    taxonSetIDs.add(prior.getID());
                }
            }
        }
    }
    // assume all calibration constraints are MonoPhyletic
    // TODO: verify that this is a reasonable assumption
    lastMonophyletic = taxonSets.size();
    // sort constraints such that if taxon set i is subset of taxon set j, then i < j
    for (int i = 0; i < lastMonophyletic; i++) {
        for (int j = i + 1; j < lastMonophyletic; j++) {
            Set<String> intersection = new LinkedHashSet<>(taxonSets.get(i));
            intersection.retainAll(taxonSets.get(j));
            if (intersection.size() > 0) {
                final boolean isSubset = taxonSets.get(i).containsAll(taxonSets.get(j));
                final boolean isSubset2 = taxonSets.get(j).containsAll(taxonSets.get(i));
                // o taxonset1 does not intersect taxonset2
                if (!(isSubset || isSubset2)) {
                    throw new IllegalArgumentException("333: Don't know how to generate a Random Tree for taxon sets that intersect, " + "but are not inclusive. Taxonset " + taxonSetIDs.get(i) + " and " + taxonSetIDs.get(j));
                }
                // swap i & j if b1 subset of b2
                if (isSubset) {
                    swap(taxonSets, i, j);
                    swap(distributions, i, j);
                    swap(m_bounds, i, j);
                    swap(taxonSetIDs, i, j);
                }
            }
        }
    }
    // build tree of mono constraints such that j is parent of i if i is a subset of j but i+1,i+2,...,j-1 are not.
    // The last one, standing for the virtual "root" of all monophyletic clades is not associated with an actual clade
    final int[] parent = new int[lastMonophyletic];
    children = new List[lastMonophyletic + 1];
    for (int i = 0; i < lastMonophyletic + 1; i++) {
        children[i] = new ArrayList<>();
    }
    for (int i = 0; i < lastMonophyletic; i++) {
        int j = i + 1;
        while (j < lastMonophyletic && !taxonSets.get(j).containsAll(taxonSets.get(i))) {
            j++;
        }
        parent[i] = j;
        children[j].add(i);
    }
    // make sure upper bounds of a child does not exceed the upper bound of its parent
    for (int i = lastMonophyletic - 1; i >= 0; --i) {
        if (parent[i] < lastMonophyletic) {
            if (m_bounds.get(i).upper > m_bounds.get(parent[i]).upper) {
                m_bounds.get(i).upper = m_bounds.get(parent[i]).upper - 1e-100;
            }
        }
    }
    final PopulationFunction popFunction = populationFunctionInput.get();
    simulateTree(taxa, popFunction);
    if (rootHeightInput.get() != null) {
        scaleToFit(rootHeightInput.get() / root.getHeight(), root);
    }
    nodeCount = 2 * taxa.size() - 1;
    internalNodeCount = taxa.size() - 1;
    leafNodeCount = taxa.size();
    HashMap<String, Integer> taxonToNR = null;
    // preserve node numbers where possible
    if (m_initial.get() != null) {
        if (leafNodeCount == m_initial.get().getLeafNodeCount()) {
            // dont ask me how the initial tree is rubbish  (i.e. 0:0.0)
            taxonToNR = new HashMap<>();
            for (Node n : m_initial.get().getExternalNodes()) {
                taxonToNR.put(n.getID(), n.getNr());
            }
        }
    } else {
        taxonToNR = new HashMap<>();
        String[] taxa = getTaxaNames();
        for (int k = 0; k < taxa.length; ++k) {
            taxonToNR.put(taxa[k], k);
        }
    }
    // multiple simulation tries may produce an excess of nodes with invalid nr's. reset those.
    setNodesNrs(root, 0, new int[1], taxonToNR);
    initArrays();
    if (m_initial.get() != null) {
        m_initial.get().assignFromWithoutID(this);
    }
    for (int k = 0; k < lastMonophyletic; ++k) {
        final MRCAPrior p = calibrations.get(k);
        if (p.isMonophyleticInput.get()) {
            final TaxonSet taxonSet = p.taxonsetInput.get();
            if (taxonSet == null) {
                throw new IllegalArgumentException("Something is wrong with constraint " + p.getID() + " -- a taxonset must be specified if a monophyletic constraint is enforced.");
            }
            final Set<String> usedTaxa = new LinkedHashSet<>();
            if (taxonSet.asStringList() == null) {
                taxonSet.initAndValidate();
            }
            usedTaxa.addAll(taxonSet.asStringList());
            /* int c = */
            traverse(root, usedTaxa, taxonSet.getTaxonCount(), new int[1]);
        // boolean b = c == nrOfTaxa + 127;
        }
    }
}

Example 2 with ParametricDistribution

use of beast.math.distributions.ParametricDistribution in project beast2 by CompEvol.

the class ParameterInputEditor method addComboBox.

@Override
protected void addComboBox(JComponent box, Input<?> input, BEASTInterface beastObject) {
    Box paramBox = Box.createHorizontalBox();
    Parameter.Base<?> parameter = null;
    if (itemNr >= 0) {
        parameter = (Parameter.Base<?>) ((List<?>) input.get()).get(itemNr);
    } else {
        parameter = (Parameter.Base<?>) input.get();
    }
    if (parameter == null) {
        super.addComboBox(box, input, beastObject);
    } else {
        setUpEntry();
        paramBox.add(m_entry);
        if (doc.allowLinking) {
            boolean isLinked = doc.isLinked(m_input);
            if (isLinked || doc.suggestedLinks((BEASTInterface) m_input.get()).size() > 0) {
                JButton linkbutton = new JButton(Utils.getIcon(BeautiPanel.ICONPATH + (isLinked ? "link.png" : "unlink.png")));
                linkbutton.setBorder(BorderFactory.createEmptyBorder());
                linkbutton.setToolTipText("link/unlink this parameter with another compatible parameter");
                linkbutton.addActionListener(e -> {
                    if (doc.isLinked(m_input)) {
                        // unlink
                        try {
                            BEASTInterface candidate = doc.getUnlinkCandidate(m_input, m_beastObject);
                            m_input.setValue(candidate, m_beastObject);
                            doc.deLink(m_input);
                        } catch (RuntimeException e2) {
                            e2.printStackTrace();
                            JOptionPane.showMessageDialog(this, "Could not unlink: " + e2.getMessage());
                        }
                    } else {
                        // create a link
                        List<BEASTInterface> candidates = doc.suggestedLinks((BEASTInterface) m_input.get());
                        JComboBox<BEASTInterface> jcb = new JComboBox<>(candidates.toArray(new BEASTInterface[] {}));
                        JOptionPane.showMessageDialog(null, jcb, "select parameter to link with", JOptionPane.QUESTION_MESSAGE);
                        BEASTInterface candidate = (BEASTInterface) jcb.getSelectedItem();
                        if (candidate != null) {
                            try {
                                m_input.setValue(candidate, m_beastObject);
                                doc.addLink(m_input);
                            } catch (Exception e2) {
                                e2.printStackTrace();
                            }
                        }
                    }
                    refreshPanel();
                });
                paramBox.add(linkbutton);
            }
        }
        paramBox.add(Box.createHorizontalGlue());
        m_isEstimatedBox = new JCheckBox(doc.beautiConfig.getInputLabel(parameter, parameter.isEstimatedInput.getName()));
        m_isEstimatedBox.setName(input.getName() + ".isEstimated");
        if (input.get() != null) {
            m_isEstimatedBox.setSelected(parameter.isEstimatedInput.get());
        }
        m_isEstimatedBox.setToolTipText(parameter.isEstimatedInput.getHTMLTipText());
        boolean isClockRate = false;
        for (Object output : parameter.getOutputs()) {
            if (output instanceof BranchRateModel.Base) {
                isClockRate |= ((BranchRateModel.Base) output).meanRateInput.get() == parameter;
            }
        }
        m_isEstimatedBox.setEnabled(!isClockRate || !getDoc().autoSetClockRate);
        m_isEstimatedBox.addActionListener(e -> {
            try {
                Parameter.Base<?> parameter2 = (Parameter.Base<?>) m_input.get();
                parameter2.isEstimatedInput.setValue(m_isEstimatedBox.isSelected(), parameter2);
                if (isParametricDistributionParameter) {
                    String id = parameter2.getID();
                    if (id.startsWith("RealParameter")) {
                        ParametricDistribution parent = null;
                        for (Object beastObject2 : parameter2.getOutputs()) {
                            if (beastObject2 instanceof ParametricDistribution) {
                                parent = (ParametricDistribution) beastObject2;
                                break;
                            }
                        }
                        Distribution grandparent = null;
                        for (Object beastObject2 : parent.getOutputs()) {
                            if (beastObject2 instanceof Distribution) {
                                grandparent = (Distribution) beastObject2;
                                break;
                            }
                        }
                        id = "parameter.hyper" + parent.getClass().getSimpleName() + "-" + m_input.getName() + "-" + grandparent.getID();
                        doc.pluginmap.remove(parameter2.getID());
                        parameter2.setID(id);
                        doc.addPlugin(parameter2);
                    }
                    PartitionContext context = new PartitionContext(id.substring("parameter.".length()));
                    Log.warning.println(context + " " + id);
                    doc.beautiConfig.hyperPriorTemplate.createSubNet(context, true);
                }
                refreshPanel();
            } catch (Exception ex) {
                Log.err.println("ParameterInputEditor " + ex.getMessage());
            }
        });
        paramBox.add(m_isEstimatedBox);
        // only show the estimate flag if there is an operator that works on this parameter
        m_isEstimatedBox.setVisible(doc.isExpertMode());
        m_isEstimatedBox.setToolTipText("Estimate value of this parameter in the MCMC chain");
        // m_bAddButtons = false;
        if (itemNr < 0) {
            for (Object beastObject2 : ((BEASTInterface) m_input.get()).getOutputs()) {
                if (beastObject2 instanceof ParametricDistribution) {
                    m_isEstimatedBox.setVisible(true);
                    isParametricDistributionParameter = true;
                    break;
                }
            }
            for (Object beastObject2 : ((BEASTInterface) m_input.get()).getOutputs()) {
                if (beastObject2 instanceof Operator) {
                    m_isEstimatedBox.setVisible(true);
                    // m_editPluginButton.setVisible(true);
                    break;
                }
            }
        } else {
            for (Object beastObject2 : ((BEASTInterface) ((List<?>) m_input.get()).get(itemNr)).getOutputs()) {
                if (beastObject2 instanceof Operator) {
                    m_isEstimatedBox.setVisible(true);
                    // m_editPluginButton.setVisible(true);
                    break;
                }
            }
        }
        box.add(paramBox);
    }
}

Example 3 with ParametricDistribution

use of beast.math.distributions.ParametricDistribution in project beast2 by CompEvol.

the class SampleOffValues method proposal.

@Override
public double proposal() {
    final BooleanParameter indicators = indicatorsInput.get(this);
    final RealParameter data = valuesInput.get(this);
    final ParametricDistribution distribution = distInput.get();
    final int idim = indicators.getDimension();
    final int offset = (data.getDimension() - 1) == idim ? 1 : 0;
    assert offset == 1 || data.getDimension() == idim : "" + idim + " (?+1) != " + data.getDimension();
    double hr = Double.NEGATIVE_INFINITY;
    if (scaleAll.get()) {
        for (int i = offset; i < idim; ++i) {
            if (!indicators.getValue(i - offset)) {
                try {
                    final double val = distribution.inverseCumulativeProbability(Randomizer.nextDouble());
                    hr += distribution.logDensity(data.getValue(i));
                    data.setValue(i, val);
                } catch (Exception e) {
                    // some distributions fail on extreme values - currently gamma
                    return Double.NEGATIVE_INFINITY;
                }
            }
        }
    } else {
        // available locations for direct sampling
        int[] loc = new int[idim];
        int locIndex = 0;
        for (int i = 0; i < idim; ++i) {
            if (!indicators.getValue(i)) {
                loc[locIndex] = i + offset;
                ++locIndex;
            }
        }
        if (locIndex > 0) {
            final int index = loc[Randomizer.nextInt(locIndex)];
            try {
                final double val = distribution.inverseCumulativeProbability(Randomizer.nextDouble());
                hr = distribution.logDensity(data.getValue(index));
                data.setValue(index, val);
            } catch (Exception e) {
                // some distributions fail on extreme values - currently gamma
                return Double.NEGATIVE_INFINITY;
            // throw new OperatorFailedException(e.getMessage());
            }
        } else {
        // no non-active indicators
        // return Double.NEGATIVE_INFINITY;
        }
    }
    return hr;
}

Example 4 with ParametricDistribution

use of beast.math.distributions.ParametricDistribution in project beast2 by CompEvol.

the class CalibratedYuleModel method compatibleInitialTree.

public Tree compatibleInitialTree() throws MathException {
    final int calCount = orderedCalibrations.length;
    final double[] lowBound = new double[calCount];
    final double[] cladeHeight = new double[calCount];
    // get lower  bound: max(lower bound of dist , bounds of nested clades)
    for (int k = 0; k < calCount; ++k) {
        final CalibrationPoint cal = orderedCalibrations[k];
        final ParametricDistribution dist = cal.dist();
        // final double offset = dist.getOffset();
        lowBound[k] = dist.inverseCumulativeProbability(0);
        // those are node heights
        if (lowBound[k] < 0) {
            lowBound[k] = 0;
        }
        for (final int i : taxaPartialOrder[k]) {
            lowBound[k] = Math.max(lowBound[k], lowBound[i]);
        }
        cladeHeight[k] = dist.inverseCumulativeProbability(1);
    // if (! Double.isInfinite(cladeHeight[k])) {
    // cladeHeight[k] += offset;
    // }
    }
    for (int k = calCount - 1; k >= 0; --k) {
        // cladeHeight[k] should be the upper bound of k
        double upper = cladeHeight[k];
        if (Double.isInfinite(upper)) {
            upper = lowBound[k] + 1;
        }
        cladeHeight[k] = (upper + lowBound[k]) / 2.0;
        for (final int i : taxaPartialOrder[k]) {
            cladeHeight[i] = Math.min(cladeHeight[i], cladeHeight[k]);
        }
    }
    final TreeInterface tree = treeInput.get();
    final int nodeCount = tree.getLeafNodeCount();
    final boolean[] used = new boolean[nodeCount];
    int curLeaf = -1;
    int curInternal = nodeCount - 1;
    final Node[] subTree = new Node[calCount];
    for (int k = 0; k < calCount; ++k) {
        final List<Integer> freeTaxa = new ArrayList<>();
        for (final int ti : xclades[k]) {
            freeTaxa.add(ti);
        }
        for (final int i : taxaPartialOrder[k]) {
            for (final int u : xclades[i]) {
                freeTaxa.remove(new Integer(u));
            }
        }
        final List<Node> sbs = new ArrayList<>();
        for (final int i : freeTaxa) {
            final Node n = new Node(tree.getNode(i).getID());
            n.setNr(++curLeaf);
            n.setHeight(0.0);
            sbs.add(n);
            used[i] = true;
        }
        for (final int i : taxaPartialOrder[k]) {
            sbs.add(subTree[i]);
            subTree[i] = null;
        }
        final double base = sbs.get(sbs.size() - 1).getHeight();
        final double step = (cladeHeight[k] - base) / (sbs.size() - 1);
        Node tr = sbs.get(0);
        for (int i = 1; i < sbs.size(); ++i) {
            tr = Node.connect(tr, sbs.get(i), base + i * step);
            tr.setNr(++curInternal);
        }
        subTree[k] = tr;
    }
    Node finalTree = subTree[calCount - 1];
    double h = cladeHeight[calCount - 1];
    for (int k = 0; k < calCount - 1; ++k) {
        final Node s = subTree[k];
        if (s != null) {
            h = Math.max(h, cladeHeight[k]) + 1;
            finalTree = Node.connect(finalTree, s, h);
            finalTree.setNr(++curInternal);
        }
    }
    for (int k = 0; k < used.length; ++k) {
        if (!used[k]) {
            final String tx = tree.getNode(k).getID();
            final Node n = new Node(tx);
            n.setHeight(0.0);
            n.setNr(++curLeaf);
            finalTree = Node.connect(finalTree, n, h + 1);
            finalTree.setNr(++curInternal);
            h += 1;
        }
    }
    final Tree t = new Tree();
    t.setRoot(finalTree);
    t.initAndValidate();
    return t;
}

Example 5 with ParametricDistribution

use of beast.math.distributions.ParametricDistribution 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);
    }
}