Examples with CompoundLikelihood dr.inference.model.CompoundLikelihood used on opensource projects

Example 6 with CompoundLikelihood

use of dr.inference.model.CompoundLikelihood in project beast-mcmc by beast-dev.

the class LognormalPriorTest method testLognormalPrior.

public void testLognormalPrior() {
    //        ConstantPopulation constant = new ConstantPopulation(Units.Type.YEARS);
    //        constant.setN0(popSize); // popSize
    Parameter popSize = new Parameter.Default(6.0);
    popSize.setId(ConstantPopulationModelParser.POPULATION_SIZE);
    ConstantPopulationModel demo = new ConstantPopulationModel(popSize, Units.Type.YEARS);
    //Likelihood
    Likelihood dummyLikelihood = new DummyLikelihood(demo);
    // Operators
    OperatorSchedule schedule = new SimpleOperatorSchedule();
    MCMCOperator operator = new ScaleOperator(popSize, 0.75);
    operator.setWeight(1.0);
    schedule.addOperator(operator);
    // Log
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    MCLogger[] loggers = new MCLogger[2];
    loggers[0] = new MCLogger(formatter, 1000, false);
    //        loggers[0].add(treeLikelihood);
    loggers[0].add(popSize);
    loggers[1] = new MCLogger(new TabDelimitedFormatter(System.out), 100000, false);
    //        loggers[1].add(treeLikelihood);
    loggers[1].add(popSize);
    // MCMC
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(1000000);
    // meanInRealSpace="false"
    DistributionLikelihood logNormalLikelihood = new DistributionLikelihood(new LogNormalDistribution(1.0, 1.0), 0);
    logNormalLikelihood.addData(popSize);
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(logNormalLikelihood);
    Likelihood prior = new CompoundLikelihood(0, likelihoods);
    likelihoods.clear();
    likelihoods.add(dummyLikelihood);
    Likelihood likelihood = new CompoundLikelihood(-1, likelihoods);
    likelihoods.clear();
    likelihoods.add(prior);
    likelihoods.add(likelihood);
    Likelihood posterior = new CompoundLikelihood(0, likelihoods);
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, posterior, schedule, loggers);
    mcmc.run();
    // time
    System.out.println(mcmc.getTimer().toString());
    // Tracer
    List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("LognormalPriorTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    //      <expectation name="param" value="4.48168907"/>
    TraceCorrelation popSizeStats = traceList.getCorrelationStatistics(traceList.getTraceIndex(ConstantPopulationModelParser.POPULATION_SIZE));
    System.out.println("Expectation of Log-Normal(1,1) is e^(M+S^2/2) = e^(1.5) = " + Math.exp(1.5));
    assertExpectation(ConstantPopulationModelParser.POPULATION_SIZE, popSizeStats, Math.exp(1.5));
}

Example 7 with CompoundLikelihood

use of dr.inference.model.CompoundLikelihood in project beast-mcmc by beast-dev.

the class CompoundLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    // the default is -1 threads (automatic thread pool size) but an XML attribute can override it
    int threads = xo.getAttribute(THREADS, -1);
    // both the XML attribute and a system property can override it
    if (System.getProperty("thread.count") != null) {
        threads = Integer.parseInt(System.getProperty("thread.count"));
        if (threads < -1 || threads > 1000) {
            // put an upper limit here - may be unnecessary?
            threads = -1;
        }
    }
    //        }
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    for (int i = 0; i < xo.getChildCount(); i++) {
        final Object child = xo.getChild(i);
        if (child instanceof Likelihood) {
            if (likelihoods.contains(child)) {
                throw new XMLParseException("The likelihood element, '" + ((Likelihood) child).getId() + "', is already present in the likelihood or prior density.");
            }
            likelihoods.add((Likelihood) child);
        //            } else if (child instanceof BeagleBranchLikelihoods){
        //                
        //            	//TODO
        //            	likelihoods.addAll( ((BeagleBranchLikelihoods)child).getBranchLikelihoods());
        } else {
            throw new XMLParseException("An element (" + child + ") which is not a likelihood has been added to a " + COMPOUND_LIKELIHOOD + " element");
        }
    }
    CompoundLikelihood compoundLikelihood;
    if (xo.getName().equalsIgnoreCase(LIKELIHOOD)) {
        compoundLikelihood = new CompoundLikelihood(threads, likelihoods);
        switch(threads) {
            case -1:
                Logger.getLogger("dr.evomodel").info("\nLikelihood computation is using an auto sizing thread pool.");
                break;
            case 0:
                Logger.getLogger("dr.evomodel").info("\nLikelihood computation is using a single thread.");
                break;
            default:
                Logger.getLogger("dr.evomodel").info("\nLikelihood computation is using a pool of " + threads + " threads.");
                break;
        }
    } else {
        compoundLikelihood = new CompoundLikelihood(likelihoods);
    }
    return compoundLikelihood;
}

Example 8 with CompoundLikelihood

use of dr.inference.model.CompoundLikelihood in project beast-mcmc by beast-dev.

the class ARGAddRemoveOperatorTest method flatPriorTester.

private void flatPriorTester(ARGModel arg, int chainLength, int sampleTreeEvery, double nodeCountSetting, double rootHeightAlpha, double rootHeightBeta, int maxCount) throws IOException, Importer.ImportException {
    MCMC mcmc = new MCMC("mcmc1");
    MCMCOptions options = new MCMCOptions(chainLength);
    //        double nodeCountSetting = 2.0;
    //        double rootHeightAlpha = 100;
    //        double rootHeightBeta = 0.5;
    OperatorSchedule schedule = getSchedule(arg);
    ARGUniformPrior uniformPrior = new ARGUniformPrior(arg, maxCount, arg.getExternalNodeCount());
    PoissonDistribution poisson = new PoissonDistribution(nodeCountSetting);
    DistributionLikelihood nodeCountPrior = new DistributionLikelihood(poisson, 0.0);
    ARGReassortmentNodeCountStatistic nodeCountStatistic = new ARGReassortmentNodeCountStatistic("nodeCount", arg);
    nodeCountPrior.addData(nodeCountStatistic);
    DistributionLikelihood rootPrior = new DistributionLikelihood(new GammaDistribution(rootHeightAlpha, rootHeightBeta), 0.0);
    CompoundParameter rootHeight = (CompoundParameter) arg.createNodeHeightsParameter(true, false, false);
    rootPrior.addData(rootHeight);
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    likelihoods.add(uniformPrior);
    likelihoods.add(rootPrior);
    likelihoods.add(nodeCountPrior);
    CompoundLikelihood compoundLikelihood = new CompoundLikelihood(1, likelihoods);
    compoundLikelihood.setId("likelihood1");
    MCLogger[] loggers = new MCLogger[3];
    loggers[0] = new MCLogger(new TabDelimitedFormatter(System.out), 10000, false);
    loggers[0].add(compoundLikelihood);
    loggers[0].add(arg);
    File file = new File("test.args");
    file.deleteOnExit();
    FileOutputStream out = new FileOutputStream(file);
    loggers[1] = new ARGLogger(arg, new TabDelimitedFormatter(out), sampleTreeEvery, "test");
    ArrayLogFormatter formatter = new ArrayLogFormatter(false);
    loggers[2] = new MCLogger(formatter, sampleTreeEvery, false);
    loggers[2].add(arg);
    arg.getRootHeightParameter().setId("root");
    loggers[2].add(arg.getRootHeightParameter());
    mcmc.setShowOperatorAnalysis(true);
    mcmc.init(options, compoundLikelihood, schedule, loggers);
    mcmc.run();
    out.flush();
    out.close();
    List<Trace> traces = formatter.getTraces();
    //        Set<String> uniqueTrees = new HashSet<String>();
    //
    //        NexusImporter importer = new NexusImporter(new FileReader(file));
    //        while (importer.hasTree()) {
    //            Tree t = importer.importNextTree();
    //            uniqueTrees.add(Tree.Utils.uniqueNewick(t, t.getRoot()));
    //        }
    //
    //        TestCase.assertEquals(numTopologies, uniqueTrees.size());            List<Trace> traces = formatter.getTraces();
    ArrayTraceList traceList = new ArrayTraceList("ARGTest", traces, 0);
    for (int i = 1; i < traces.size(); i++) {
        traceList.analyseTrace(i);
    }
    TraceCorrelation nodeCountStats = traceList.getCorrelationStatistics(1);
    TraceCorrelation rootHeightStats = traceList.getCorrelationStatistics(4);
    assertExpectation("nodeCount", nodeCountStats, poisson.truncatedMean(maxCount));
    assertExpectation(TreeModelParser.ROOT_HEIGHT, rootHeightStats, rootHeightAlpha * rootHeightBeta);
}

Example 9 with CompoundLikelihood

use of dr.inference.model.CompoundLikelihood 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 10 with CompoundLikelihood

use of dr.inference.model.CompoundLikelihood in project beast-mcmc by beast-dev.

the class MarkovChain method runChain.

/**
     * Run the chain for a given number of states.
     *
     * @param length number of states to run the chain.
     */
public long runChain(long length, boolean disableCoerce) {
    likelihood.makeDirty();
    currentScore = evaluate(likelihood);
    long currentState = currentLength;
    final Model currentModel = likelihood.getModel();
    if (currentState == 0) {
        initialScore = currentScore;
        bestScore = currentScore;
        fireBestModel(currentState, currentModel);
    }
    if (currentScore == Double.NEGATIVE_INFINITY) {
        // identify which component of the score is zero...
        String message = "The initial likelihood is zero";
        if (likelihood instanceof CompoundLikelihood) {
            message += ": " + ((CompoundLikelihood) likelihood).getDiagnosis();
        } else if (likelihood instanceof PathLikelihood) {
            message += ": " + ((CompoundLikelihood) ((PathLikelihood) likelihood).getSourceLikelihood()).getDiagnosis();
            message += ": " + ((CompoundLikelihood) ((PathLikelihood) likelihood).getDestinationLikelihood()).getDiagnosis();
        } else {
            message += ".";
        }
        throw new IllegalArgumentException(message);
    } else if (currentScore == Double.POSITIVE_INFINITY || Double.isNaN(currentScore)) {
        String message = "A likelihood returned with a numerical error";
        if (likelihood instanceof CompoundLikelihood) {
            message += ": " + ((CompoundLikelihood) likelihood).getDiagnosis();
        } else {
            message += ".";
        }
        throw new IllegalArgumentException(message);
    }
    pleaseStop = false;
    isStopped = false;
    //int otfcounter = onTheFlyOperatorWeights > 0 ? onTheFlyOperatorWeights : 0;
    double[] logr = { 0.0 };
    boolean usingFullEvaluation = true;
    // set ops count in mcmc element instead
    if (// Temporary solution until full code review
    fullEvaluationCount == 0)
        usingFullEvaluation = false;
    boolean fullEvaluationError = false;
    while (!pleaseStop && (currentState < (currentLength + length))) {
        String diagnosticStart = "";
        // periodically log states
        fireCurrentModel(currentState, currentModel);
        if (pleaseStop) {
            isStopped = true;
            break;
        }
        // Get the operator
        final int op = schedule.getNextOperatorIndex();
        final MCMCOperator mcmcOperator = schedule.getOperator(op);
        double oldScore = currentScore;
        if (usingFullEvaluation) {
            diagnosticStart = likelihood instanceof CompoundLikelihood ? ((CompoundLikelihood) likelihood).getDiagnosis() : "";
        }
        // The current model is stored here in case the proposal fails
        if (currentModel != null) {
            currentModel.storeModelState();
        }
        // assert Profiler.stopProfile("Store");
        boolean operatorSucceeded = true;
        double hastingsRatio = 1.0;
        boolean accept = false;
        logr[0] = -Double.MAX_VALUE;
        long elaspedTime = 0;
        if (PROFILE) {
            elaspedTime = System.currentTimeMillis();
        }
        if (DEBUG) {
            System.out.println("\n>> Iteration: " + currentState);
            System.out.println("\n&& Operator: " + mcmcOperator.getOperatorName());
        }
        if (mcmcOperator instanceof GeneralOperator) {
            hastingsRatio = ((GeneralOperator) mcmcOperator).operate(likelihood);
        } else {
            hastingsRatio = mcmcOperator.operate();
        }
        // assert Profiler.stopProfile("Operate");
        if (hastingsRatio == Double.NEGATIVE_INFINITY) {
            // Should the evaluation be short-cutted?
            // Previously this was set to false if OperatorFailedException was thrown.
            // Now a -Inf HR is returned.
            operatorSucceeded = false;
        }
        if (PROFILE) {
            long duration = System.currentTimeMillis() - elaspedTime;
            if (DEBUG) {
                System.out.println("Time: " + duration);
            }
            mcmcOperator.addEvaluationTime(duration);
        }
        double score = Double.NaN;
        double deviation = Double.NaN;
        //    System.err.print("" + currentState + ": ");
        if (operatorSucceeded) {
            if (DEBUG) {
                System.out.println("** Evaluate");
            }
            long elapsedTime = 0;
            if (PROFILE) {
                elapsedTime = System.currentTimeMillis();
            }
            // The new model is evaluated
            score = evaluate(likelihood);
            if (PROFILE) {
                long duration = System.currentTimeMillis() - elapsedTime;
                if (DEBUG) {
                    System.out.println("Time: " + duration);
                }
                mcmcOperator.addEvaluationTime(duration);
            }
            String diagnosticOperator = "";
            if (usingFullEvaluation) {
                diagnosticOperator = likelihood instanceof CompoundLikelihood ? ((CompoundLikelihood) likelihood).getDiagnosis() : "";
            }
            if (score == Double.NEGATIVE_INFINITY && mcmcOperator instanceof GibbsOperator) {
                if (!(mcmcOperator instanceof GibbsIndependentNormalDistributionOperator) && !(mcmcOperator instanceof GibbsIndependentGammaOperator) && !(mcmcOperator instanceof GibbsIndependentCoalescentOperator) && !(mcmcOperator instanceof GibbsIndependentJointNormalGammaOperator)) {
                    Logger.getLogger("error").severe("State " + currentState + ": A Gibbs operator, " + mcmcOperator.getOperatorName() + ", returned a state with zero likelihood.");
                }
            }
            if (score == Double.POSITIVE_INFINITY || Double.isNaN(score)) {
                if (likelihood instanceof CompoundLikelihood) {
                    Logger.getLogger("error").severe("State " + currentState + ": A likelihood returned with a numerical error:\n" + ((CompoundLikelihood) likelihood).getDiagnosis());
                } else {
                    Logger.getLogger("error").severe("State " + currentState + ": A likelihood returned with a numerical error.");
                }
                // If the user has chosen to ignore this error then we transform it
                // to a negative infinity so the state is rejected.
                score = Double.NEGATIVE_INFINITY;
            }
            if (usingFullEvaluation) {
                // This is a test that the state was correctly evaluated. The
                // likelihood of all components of the model are flagged as
                // needing recalculation, then the full likelihood is calculated
                // again and compared to the first result. This checks that the
                // BEAST is aware of all changes that the operator induced.
                likelihood.makeDirty();
                final double testScore = evaluate(likelihood);
                final String d2 = likelihood instanceof CompoundLikelihood ? ((CompoundLikelihood) likelihood).getDiagnosis() : "";
                if (Math.abs(testScore - score) > evaluationTestThreshold) {
                    Logger.getLogger("error").severe("State " + currentState + ": State was not correctly calculated after an operator move.\n" + "Likelihood evaluation: " + score + "\nFull Likelihood evaluation: " + testScore + "\n" + "Operator: " + mcmcOperator + " " + mcmcOperator.getOperatorName() + (diagnosticOperator.length() > 0 ? "\n\nDetails\nBefore: " + diagnosticOperator + "\nAfter: " + d2 : "") + "\n\n");
                    fullEvaluationError = true;
                }
            }
            if (score > bestScore) {
                bestScore = score;
                fireBestModel(currentState, currentModel);
            }
            accept = mcmcOperator instanceof GibbsOperator || acceptor.accept(oldScore, score, hastingsRatio, logr);
            deviation = score - oldScore;
        }
        // The new model is accepted or rejected
        if (accept) {
            if (DEBUG) {
                System.out.println("** Move accepted: new score = " + score + ", old score = " + oldScore);
            }
            mcmcOperator.accept(deviation);
            currentModel.acceptModelState();
            currentScore = score;
        } else {
            if (DEBUG) {
                System.out.println("** Move rejected: new score = " + score + ", old score = " + oldScore + " (logr = " + logr[0] + ")");
            }
            mcmcOperator.reject();
            // assert Profiler.startProfile("Restore");
            currentModel.restoreModelState();
            if (usingFullEvaluation) {
                // This is a test that the state is correctly restored. The
                // restored state is fully evaluated and the likelihood compared with
                // that before the operation was made.
                likelihood.makeDirty();
                final double testScore = evaluate(likelihood);
                final String d2 = likelihood instanceof CompoundLikelihood ? ((CompoundLikelihood) likelihood).getDiagnosis() : "";
                if (Math.abs(testScore - oldScore) > evaluationTestThreshold) {
                    final Logger logger = Logger.getLogger("error");
                    logger.severe("State " + currentState + ": State was not correctly restored after reject step.\n" + "Likelihood before: " + oldScore + " Likelihood after: " + testScore + "\n" + "Operator: " + mcmcOperator + " " + mcmcOperator.getOperatorName() + (diagnosticStart.length() > 0 ? "\n\nDetails\nBefore: " + diagnosticStart + "\nAfter: " + d2 : "") + "\n\n");
                    fullEvaluationError = true;
                }
            }
        }
        if (!disableCoerce && mcmcOperator instanceof CoercableMCMCOperator) {
            coerceAcceptanceProbability((CoercableMCMCOperator) mcmcOperator, logr[0]);
        }
        if (usingFullEvaluation) {
            if (schedule.getMinimumAcceptAndRejectCount() >= minOperatorCountForFullEvaluation && currentState >= fullEvaluationCount) {
                // full evaluation is only switched off when each operator has done a
                // minimum number of operations (currently 1) and fullEvalationCount
                // operations in total.
                usingFullEvaluation = false;
                if (fullEvaluationError) {
                    // If there has been an error then stop with an error
                    throw new RuntimeException("One or more evaluation errors occurred during the test phase of this\n" + "run. These errors imply critical errors which may produce incorrect\n" + "results.");
                }
            }
        }
        fireEndCurrentIteration(currentState);
        currentState += 1;
    }
    currentLength = currentState;
    return currentLength;
}