Examples with SpeciationLikelihood dr.evomodel.speciation.SpeciationLikelihood used on opensource projects

Example 6 with SpeciationLikelihood

use of dr.evomodel.speciation.SpeciationLikelihood in project beast-mcmc by beast-dev.

the class SpeciationLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    XMLObject cxo = xo.getChild(MODEL);
    final SpeciationModel specModel = (SpeciationModel) cxo.getChild(SpeciationModel.class);
    cxo = xo.getChild(TREE);
    final Tree tree = (Tree) cxo.getChild(Tree.class);
    Set<Taxon> excludeTaxa = null;
    if (xo.hasChildNamed(INCLUDE)) {
        excludeTaxa = new HashSet<Taxon>();
        for (int i = 0; i < tree.getTaxonCount(); i++) {
            excludeTaxa.add(tree.getTaxon(i));
        }
        cxo = xo.getChild(INCLUDE);
        for (int i = 0; i < cxo.getChildCount(); i++) {
            TaxonList taxonList = (TaxonList) cxo.getChild(i);
            for (int j = 0; j < taxonList.getTaxonCount(); j++) {
                excludeTaxa.remove(taxonList.getTaxon(j));
            }
        }
    }
    if (xo.hasChildNamed(EXCLUDE)) {
        excludeTaxa = new HashSet<Taxon>();
        cxo = xo.getChild(EXCLUDE);
        for (int i = 0; i < cxo.getChildCount(); i++) {
            TaxonList taxonList = (TaxonList) cxo.getChild(i);
            for (int j = 0; j < taxonList.getTaxonCount(); j++) {
                excludeTaxa.add(taxonList.getTaxon(j));
            }
        }
    }
    if (excludeTaxa != null) {
        Logger.getLogger("dr.evomodel").info("Speciation model excluding " + excludeTaxa.size() + " taxa from prior - " + (tree.getTaxonCount() - excludeTaxa.size()) + " taxa remaining.");
    }
    final XMLObject cal = xo.getChild(CALIBRATION);
    if (cal != null) {
        if (excludeTaxa != null) {
            throw new XMLParseException("Sorry, not implemented: internal calibration prior + excluded taxa");
        }
        List<Distribution> dists = new ArrayList<Distribution>();
        List<Taxa> taxa = new ArrayList<Taxa>();
        List<Boolean> forParent = new ArrayList<Boolean>();
        // (Statistic) cal.getChild(Statistic.class);
        Statistic userPDF = null;
        for (int k = 0; k < cal.getChildCount(); ++k) {
            final Object ck = cal.getChild(k);
            if (DistributionLikelihood.class.isInstance(ck)) {
                dists.add(((DistributionLikelihood) ck).getDistribution());
            } else if (Distribution.class.isInstance(ck)) {
                dists.add((Distribution) ck);
            } else if (Taxa.class.isInstance(ck)) {
                final Taxa tx = (Taxa) ck;
                taxa.add(tx);
                forParent.add(tx.getTaxonCount() == 1);
            } else if (Statistic.class.isInstance(ck)) {
                if (userPDF != null) {
                    throw new XMLParseException("more than one userPDF correction???");
                }
                userPDF = (Statistic) cal.getChild(Statistic.class);
            } else {
                XMLObject cko = (XMLObject) ck;
                assert cko.getChildCount() == 2;
                for (int i = 0; i < 2; ++i) {
                    final Object chi = cko.getChild(i);
                    if (DistributionLikelihood.class.isInstance(chi)) {
                        dists.add(((DistributionLikelihood) chi).getDistribution());
                    } else if (Distribution.class.isInstance(chi)) {
                        dists.add((Distribution) chi);
                    } else if (Taxa.class.isInstance(chi)) {
                        taxa.add((Taxa) chi);
                        boolean fp = ((Taxa) chi).getTaxonCount() == 1;
                        if (cko.hasAttribute(PARENT)) {
                            boolean ufp = cko.getBooleanAttribute(PARENT);
                            if (fp && !ufp) {
                                throw new XMLParseException("forParent==false for a single taxon?? (must be true)");
                            }
                            fp = ufp;
                        }
                        forParent.add(fp);
                    } else {
                        assert false;
                    }
                }
            }
        }
        if (dists.size() != taxa.size()) {
            throw new XMLParseException("Mismatch in number of distributions and taxa specs");
        }
        try {
            final String correction = cal.getAttribute(CORRECTION, EXACT);
            final CalibrationPoints.CorrectionType type = correction.equals(EXACT) ? CalibrationPoints.CorrectionType.EXACT : (correction.equals(APPROX) ? CalibrationPoints.CorrectionType.APPROXIMATED : (correction.equals(NONE) ? CalibrationPoints.CorrectionType.NONE : (correction.equals(PEXACT) ? CalibrationPoints.CorrectionType.PEXACT : null)));
            if (cal.hasAttribute(CORRECTION) && type == null) {
                throw new XMLParseException("correction type == " + correction + "???");
            }
            final CalibrationPoints calib = new CalibrationPoints(tree, specModel.isYule(), dists, taxa, forParent, userPDF, type);
            final SpeciationLikelihood speciationLikelihood = new SpeciationLikelihood(tree, specModel, null, calib);
            return speciationLikelihood;
        } catch (IllegalArgumentException e) {
            throw new XMLParseException(e.getMessage());
        }
    }
    return new SpeciationLikelihood(tree, specModel, excludeTaxa, null);
}

Example 7 with SpeciationLikelihood

use of dr.evomodel.speciation.SpeciationLikelihood in project beast-mcmc by beast-dev.

the class TestCalibratedYuleModel method yuleTester.

private void yuleTester(TreeModel treeModel, OperatorSchedule schedule, Parameter brParameter, double S, int chainLength) throws IOException, TreeUtils.MissingTaxonException {
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(chainLength);
    TreeLengthStatistic tls = new TreeLengthStatistic(TL, treeModel);
    TreeHeightStatistic rootHeight = new TreeHeightStatistic(TREE_HEIGHT, treeModel);
    SpeciationModel speciationModel = new BirthDeathGernhard08Model("yule", brParameter, null, null, BirthDeathGernhard08Model.TreeType.UNSCALED, Units.Type.SUBSTITUTIONS, false);
    Likelihood speciationLikelihood = new SpeciationLikelihood(treeModel, speciationModel, "yule.like");
    Taxa halfTaxa = new Taxa();
    for (int i = 0; i < taxa.getTaxonCount() / 2; i++) {
        halfTaxa.addTaxon(new Taxon("T" + Integer.toString(i)));
    }
    TMRCAStatistic tmrca = new TMRCAStatistic("tmrca(halfTaxa)", treeModel, halfTaxa, false, false);
    DistributionLikelihood logNormalLikelihood = new DistributionLikelihood(new LogNormalDistribution(M, S), // meanInRealSpace="false"
    0);
    logNormalLikelihood.addData(tmrca);
    MonophylyStatistic monophylyStatistic = new MonophylyStatistic("monophyly(halfTaxa)", treeModel, halfTaxa, null);
    BooleanLikelihood booleanLikelihood = new BooleanLikelihood();
    booleanLikelihood.addData(monophylyStatistic);
    //CompoundLikelihood
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(speciationLikelihood);
    likelihoods.add(logNormalLikelihood);
    likelihoods.add(booleanLikelihood);
    Likelihood prior = new CompoundLikelihood(0, likelihoods);
    prior.setId(CompoundLikelihoodParser.PRIOR);
    ArrayLogFormatter logformatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[1];
    loggers[0] = new MCLogger(logformatter, (int) (options.getChainLength() / 10000), false);
    loggers[0].add(speciationLikelihood);
    loggers[0].add(rootHeight);
    loggers[0].add(tmrca);
    loggers[0].add(tls);
    loggers[0].add(brParameter);
    mcmc.setShowOperatorAnalysis(false);
    mcmc.init(options, prior, schedule, loggers);
    mcmc.run();
    List<Trace> traces = logformatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("yuleModelTest", traces, 1000);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    NumberFormatter formatter = new NumberFormatter(8);
    TraceCorrelation tlStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(TL));
    TraceCorrelation treeHeightStats = traceList.getCorrelationStatistics(traceList.getTraceIndex("tmrca(halfTaxa)"));
    //        out.write("tmrcaHeight = \t");
    out.write(formatter.format(treeHeightStats.getMean()));
    out.write("\t");
    double expectedNodeHeight = Math.pow(Math.E, (M + (Math.pow(S, 2) / 2)));
    //        out.write("expectation = \t");
    out.write(formatter.format(expectedNodeHeight));
    out.write("\t");
    double error = Math.abs((treeHeightStats.getMean() - expectedNodeHeight) / expectedNodeHeight);
    NumberFormat percentFormatter = NumberFormat.getNumberInstance();
    percentFormatter.setMinimumFractionDigits(5);
    percentFormatter.setMinimumFractionDigits(5);
    //        out.write("error = \t");
    out.write(percentFormatter.format(error));
    out.write("\t");
    //        out.write("tl.ess = \t");
    out.write(Double.toString(tlStats.getESS()));
    System.out.println("tmrcaHeight = " + formatter.format(treeHeightStats.getMean()) + ";  expectation = " + formatter.format(expectedNodeHeight) + ";  error = " + percentFormatter.format(error) + ";  tl.ess = " + tlStats.getESS());
}

Example 8 with SpeciationLikelihood

use of dr.evomodel.speciation.SpeciationLikelihood in project beast-mcmc by beast-dev.

the class BirthDeathSSLikelihoodTest method testBirthDeathLikelihoodBEAST2.

public void testBirthDeathLikelihoodBEAST2() {
    System.out.println("RootHeight = " + tree2.getRootHeight());
    Variable<Double> origin = new Variable.D("origin", 6.0);
    final double birthRate = 2.0;
    final double deathRate = 1.0;
    // rate of sampling taxa through time
    final double psiRate = 0.5;
    // the proportion of taxa sampled, default to fix to 0
    final double sampleProbability = 0.0;
    final boolean hasFinalSample = false;
    Variable<Double> b = new Variable.D("b", birthRate);
    Variable<Double> d = new Variable.D("d", deathRate);
    Variable<Double> psi = new Variable.D("psi", psiRate);
    Variable<Double> p = new Variable.D("p", sampleProbability);
    // sampleBecomesNonInfectiousProb
    Variable<Double> r = new Variable.D("r", 0.0);
    SpeciationModel speciationModel = new BirthDeathSerialSamplingModel(b, d, psi, p, false, r, hasFinalSample, origin, Units.Type.YEARS);
    Likelihood likelihood = new SpeciationLikelihood(tree2, speciationModel, "bdss.like");
    assertEquals(-19.0198, likelihood.getLogLikelihood(), 1e-5);
}

Example 9 with SpeciationLikelihood

use of dr.evomodel.speciation.SpeciationLikelihood in project beast-mcmc by beast-dev.

the class YuleLikelihoodTest method yuleLikelihoodTester.

private void yuleLikelihoodTester(Tree tree, double birthRate, double logL) {
    Parameter b = new Parameter.Default("b", birthRate, 0.0, Double.MAX_VALUE);
    Parameter d = new Parameter.Default("d", 0.0, 0.0, Double.MAX_VALUE);
    SpeciationModel speciationModel = new BirthDeathGernhard08Model(b, d, null, BirthDeathGernhard08Model.TreeType.TIMESONLY, Units.Type.YEARS);
    Likelihood likelihood = new SpeciationLikelihood(tree, speciationModel, "yule.like");
    assertEquals(likelihood.getLogLikelihood(), logL);
}