Examples with TaxonSet beast.evolution.alignment.TaxonSet used on opensource projects

Example 26 with TaxonSet

use of beast.evolution.alignment.TaxonSet in project beast2 by CompEvol.

the class TreeParser method initAndValidate.

/**
 * Ensure the class behaves properly, even when inputs are not specified.
 */
@Override
public void initAndValidate() {
    boolean sortNodesAlphabetically = false;
    if (dataInput.get() != null) {
        labels = dataInput.get().getTaxaNames();
    } else if (m_taxonset.get() != null) {
        if (labels == null) {
            labels = m_taxonset.get().asStringList();
        } else {
            // else labels were set by TreeParser c'tor
            sortNodesAlphabetically = true;
        }
    } else {
        if (isLabelledNewickInput.get()) {
            if (m_initial.get() != null) {
                labels = m_initial.get().getTaxonset().asStringList();
            } else {
                labels = new ArrayList<>();
                createUnrecognizedTaxa = true;
                sortNodesAlphabetically = true;
            }
        } else {
            if (m_initial.get() != null) {
                // try to pick up taxa from initial tree
                final Tree tree = m_initial.get();
                if (tree.m_taxonset.get() != null) {
                    labels = tree.m_taxonset.get().asStringList();
                } else {
                // m_sLabels = null;
                }
            } else {
            // m_sLabels = null;
            }
        }
    // m_bIsLabelledNewick = false;
    }
    final String newick = newickInput.get();
    if (newick == null || newick.equals("")) {
        // can happen while initalising Beauti
        final Node dummy = new Node();
        setRoot(dummy);
    } else {
        setRoot(parseNewick(newickInput.get()));
    }
    super.initAndValidate();
    m_sTaxaNames = null;
    if (sortNodesAlphabetically) {
        // correct for node ordering: ensure order is alphabetical
        for (int i = 0; i < getNodeCount() && i < labels.size(); i++) {
            m_nodes[i].setID(labels.get(i));
        }
        Node[] nodes = new Node[labels.size()];
        System.arraycopy(m_nodes, 0, nodes, 0, labels.size());
        Arrays.sort(nodes, (o1, o2) -> o1.getID().compareTo(o2.getID()));
        for (int i = 0; i < labels.size(); i++) {
            m_nodes[i] = nodes[i];
            nodes[i].setNr(i);
        }
    }
    if (m_initial.get() != null)
        processTraits(m_initial.get().m_traitList.get());
    else
        processTraits(m_traitList.get());
    if (timeTraitSet != null) {
        adjustTreeNodeHeights(root);
    } else if (adjustTipHeightsInput.get()) {
        double treeLength = TreeUtils.getTreeLength(this, getRoot());
        double extraTreeLength = 0.0;
        double maxTipHeight = 0.0;
        // all nodes should be at zero height if no date-trait is available
        for (int i = 0; i < getLeafNodeCount(); i++) {
            double height = getNode(i).getHeight();
            if (maxTipHeight < height) {
                maxTipHeight = height;
            }
            extraTreeLength += height;
            getNode(i).setHeight(0);
        }
        double scaleFactor = (treeLength + extraTreeLength) / treeLength;
        final double SCALE_FACTOR_THRESHOLD = 0.001;
        // if the change in total tree length is more than 0.1% then give the user a warning!
        if (scaleFactor > 1.0 + SCALE_FACTOR_THRESHOLD) {
            DecimalFormat format = new DecimalFormat("#.##");
            Log.info.println("WARNING: Adjust tip heights attribute set to 'true' in " + getClass().getSimpleName());
            Log.info.println("         has resulted in significant (>" + format.format(SCALE_FACTOR_THRESHOLD * 100.0) + "%) change in tree length.");
            Log.info.println("         Use " + adjustTipHeightsInput.getName() + "='false' to override this default.");
            Log.info.printf("  original max tip age = %8.3f\n", maxTipHeight);
            Log.info.printf("       new max tip age = %8.3f\n", 0.0);
            Log.info.printf("  original tree length = %8.3f\n", treeLength);
            Log.info.printf("       new tree length = %8.3f\n", treeLength + extraTreeLength);
            Log.info.printf("       TL scale factor = %8.3f\n", scaleFactor);
        }
    }
    if (m_taxonset.get() == null && labels != null && isLabelledNewickInput.get()) {
        m_taxonset.setValue(new TaxonSet(Taxon.createTaxonList(labels)), this);
    }
    initStateNodes();
}

Example 27 with TaxonSet

use of beast.evolution.alignment.TaxonSet 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 28 with TaxonSet

use of beast.evolution.alignment.TaxonSet in project beast2 by CompEvol.

the class FilteredAlignmentTest method testWeightedSitesReordered.

public void testWeightedSitesReordered() throws Exception {
    // reorder taxa
    Alignment data = getAlignmentNoTInHuman();
    data.setID("data");
    List<Taxon> taxa = new ArrayList<>();
    taxa.add(new Taxon("1chimp"));
    taxa.add(new Taxon("0human"));
    TaxonSet set = new TaxonSet(taxa);
    data.taxonSetInput.setValue(set, data);
    data.siteWeightsInput.setValue("11232, 2, 3, 4 ,1123,2,3,4,112,2,3,4,11,2,3,	4 ", data);
    data.initAndValidate();
    String weights = Arrays.toString(data.getWeights());
    System.out.println(weights + "\n0human\t" + alignmentToString(data, data.getTaxonIndex("0human")) + "\n1chimp\t" + alignmentToString(data, data.getTaxonIndex("1chimp")));
    assertEquals("[11243, 1123, 112, 4, 2, 2, 6, 3, 3, 8, 4, 4]", weights);
}

Example 29 with TaxonSet

use of beast.evolution.alignment.TaxonSet in project beast2 by CompEvol.

the class UnorderedAlignmentsTest method testUnorderedAlignment.

@Test
public void testUnorderedAlignment() throws Exception {
    TaxonSet taxa = getTaxa();
    Tree tree = getTree(taxa);
    SiteModel siteModel = getSiteModel();
    double logP = 0.0;
    double shuffledLogP = 0.0;
    for (int i = 0; i < 3; ++i) {
        Alignment data = getAlignment(taxa, tree, siteModel);
        // First calculate in order
        TreeLikelihood likelihood = new TreeLikelihood();
        likelihood.initByName("data", data, "tree", tree, "siteModel", siteModel);
        logP += likelihood.calculateLogP();
        // Now calculate again, with shuffled taxon order
        Collections.shuffle(data.sequenceInput.get());
        likelihood = new TreeLikelihood();
        likelihood.initByName("data", data, "tree", tree, "siteModel", siteModel);
        shuffledLogP += likelihood.calculateLogP();
    }
    assertEquals(logP, shuffledLogP, 1E-9);
}

Example 30 with TaxonSet

use of beast.evolution.alignment.TaxonSet in project beast2 by CompEvol.

the class UnorderedAlignmentsTest method getDates.

public static TraitSet getDates(TaxonSet taxa) throws Exception {
    TraitSet timeTrait = new TraitSet();
    String trait = String.join(",", (Iterable<String>) IntStream.range(0, 16).mapToObj(i -> taxa.getTaxonId(i) + "=" + i / 3.0)::iterator);
    timeTrait.initByName("traitname", "date-forward", "taxa", taxa, "value", trait);
    return timeTrait;
}