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;
}
}
}
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;
}
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;
}
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();
}
}
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;
}
}
}
}
}
}
Aggregations