Examples with RealParameter beast.core.parameter.RealParameter used on opensource projects

Example 6 with RealParameter

use of beast.core.parameter.RealParameter in project beast2 by CompEvol.

the class WilsonBaldingTest method topologyDistribution.

/**
 * Test topology distribution.
 * @throws Exception
 */
@Test
public void topologyDistribution() throws Exception {
    // Fix seed: will hopefully ensure success of test unless something
    // goes terribly wrong.
    Randomizer.setSeed(42);
    // Assemble model:
    ConstantPopulation constantPop = new ConstantPopulation();
    constantPop.initByName("popSize", new RealParameter("10000.0"));
    List<Object> alignmentInitArgs = new ArrayList<Object>();
    for (int i = 0; i < 4; i++) {
        Sequence thisSeq = new Sequence();
        thisSeq.initByName("taxon", String.valueOf(i), "value", "?");
        alignmentInitArgs.add("sequence");
        alignmentInitArgs.add(thisSeq);
    }
    Alignment alignment = new Alignment();
    alignment.initByName(alignmentInitArgs.toArray());
    Tree tree = new RandomTree();
    tree.initByName("taxa", alignment, "populationModel", constantPop);
    TreeIntervals treeIntervals = new TreeIntervals();
    treeIntervals.initByName("tree", tree);
    Coalescent coalescentDistrib = new Coalescent();
    coalescentDistrib.initByName("treeIntervals", treeIntervals, "populationModel", constantPop);
    // Set up state:
    State state = new State();
    state.initByName("stateNode", tree);
    // Set up operator:
    WilsonBalding wilsonBalding = new WilsonBalding();
    wilsonBalding.initByName("weight", "1", "tree", tree);
    // Set up logger:
    TreeReport treeReport = new TreeReport();
    treeReport.initByName("logEvery", "100", "burnin", "200000", "credibleSetPercentage", "95.0", "log", tree, "silent", true);
    // Set up MCMC:
    MCMC mcmc = new MCMC();
    mcmc.initByName("chainLength", "2000000", "state", state, "distribution", coalescentDistrib, "operator", wilsonBalding, "logger", treeReport);
    // Run MCMC:
    mcmc.run();
    // Obtain analysis results:
    TreeTraceAnalysis analysis = treeReport.getAnalysis();
    Map<String, Integer> topologyCounts = analysis.getTopologyCounts();
    int totalTreesUsed = analysis.getNTrees();
    // Test topology distribution against ideal:
    double tol = 0.005;
    for (int i = 0; i < topologies.length; i++) {
        double thisProb = topologyCounts.get(topologies[i]) / (double) totalTreesUsed;
        boolean withinTol = (thisProb > probs[i] - tol && thisProb < probs[i] + tol);
        Assert.assertTrue(withinTol);
        System.err.format("Topology %s rel. freq. %.3f", topologies[i], thisProb);
        if (withinTol)
            System.err.println(" (Within tolerance " + tol + " of " + String.valueOf(probs[i]) + ")");
        else
            System.err.println(" (FAILURE: outside tolerance " + tol + " of " + String.valueOf(probs[i]) + ")");
    }
}

Example 7 with RealParameter

use of beast.core.parameter.RealParameter in project beast2 by CompEvol.

the class YuleModelTest method testYule.

@Test
public void testYule() throws Exception {
    YuleModel bdssm = new YuleModel();
    Tree tree1 = new TreeParser("((A:1.0,B:1.0):1.0,C:2.0);", false);
    bdssm.setInputValue("tree", tree1);
    bdssm.setInputValue("birthDiffRate", new RealParameter("10."));
    bdssm.setInputValue("originHeight", new RealParameter("10."));
    bdssm.initAndValidate();
    double logP1 = bdssm.calculateTreeLogLikelihood(tree1);
    Tree tree = new TreeParser("((A:1.0,B:1.0):2.0,C:3.0);", false);
    bdssm.setInputValue("tree", tree);
    bdssm.setInputValue("birthDiffRate", new RealParameter("10."));
    bdssm.setInputValue("originHeight", new RealParameter("10."));
    bdssm.initAndValidate();
    double logP2 = bdssm.calculateTreeLogLikelihood(tree);
    assertEquals(logP1 - logP2, 10.0);
}

Example 8 with RealParameter

use of beast.core.parameter.RealParameter in project beast2 by CompEvol.

the class LoggerTest method isLoggingToStdout.

@Test
public void isLoggingToStdout() throws Exception {
    logger = new Logger();
    logger.initByName("fileName", null, "log", new RealParameter());
    assertTrue("fileName == null", logger.isLoggingToStdout());
    logger = new Logger();
    logger.initByName("fileName", "", "log", new RealParameter());
    assertTrue("fileName.length() == 0", logger.isLoggingToStdout());
    logger = new Logger();
    logger.initByName("fileName", "beast.log", "log", new RealParameter());
    assertFalse("fileName = \"beast.log\"", logger.isLoggingToStdout());
}

Example 9 with RealParameter

use of beast.core.parameter.RealParameter in project beast2 by CompEvol.

the class NexusParser method getMRCAPrior.

// TODO mv to somewhere easy to access ?
/**
 * get a MRCAPrior object for given taxon set,
 * from a string array which determines the distribution
 * @param taxonset
 * @param strs3 [0] is distribution name,
 *              [1]-[3] for values to determine the distribution
 * @return a MRCAPrior object
 * @throws RuntimeException
 */
public MRCAPrior getMRCAPrior(TaxonSet taxonset, String[] strs3) throws RuntimeException {
    RealParameter[] param = new RealParameter[strs3.length];
    for (int i = 1; i < strs3.length; i++) {
        try {
            param[i] = new RealParameter(strs3[i]);
            param[i].setID("param." + i);
        } catch (Exception e) {
        // ignore parsing errors
        }
    }
    ParametricDistribution distr = null;
    switch(strs3[0]) {
        case "normal":
            distr = new Normal();
            distr.initByName("mean", param[1], "sigma", param[2]);
            distr.setID("Normal.0");
            break;
        case "uniform":
            distr = new Uniform();
            distr.initByName("lower", strs3[1], "upper", strs3[2]);
            distr.setID("Uniform.0");
            break;
        case "fixed":
            // uniform with lower == upper
            distr = new Normal();
            distr.initByName("mean", param[1], "sigma", "+Infinity");
            distr.setID("Normal.0");
            break;
        case "offsetlognormal":
            distr = new LogNormalDistributionModel();
            distr.initByName("offset", strs3[1], "M", param[2], "S", param[3], "meanInRealSpace", true);
            distr.setID("LogNormalDistributionModel.0");
            break;
        case "lognormal":
            distr = new LogNormalDistributionModel();
            distr.initByName("M", param[1], "S", param[2], "meanInRealSpace", true);
            distr.setID("LogNormalDistributionModel.0");
            break;
        case "offsetexponential":
            distr = new Exponential();
            distr.initByName("offset", strs3[1], "mean", param[2]);
            distr.setID("Exponential.0");
            break;
        case "gamma":
            distr = new Gamma();
            distr.initByName("alpha", param[1], "beta", param[2]);
            distr.setID("Gamma.0");
            break;
        case "offsetgamma":
            distr = new Gamma();
            distr.initByName("offset", strs3[1], "alpha", param[2], "beta", param[3]);
            distr.setID("Gamma.0");
            break;
        default:
            throw new RuntimeException("Unknwon distribution " + strs3[0]);
    }
    MRCAPrior prior = new MRCAPrior();
    prior.isMonophyleticInput.setValue(true, prior);
    prior.distInput.setValue(distr, prior);
    prior.taxonsetInput.setValue(taxonset, prior);
    prior.setID(taxonset.getID() + ".prior");
    return prior;
}

Example 10 with RealParameter

use of beast.core.parameter.RealParameter in project beast2 by CompEvol.

the class ScaleOperator method proposal.

/**
 * override this for proposals,
 *
 * @return log of Hastings Ratio, or Double.NEGATIVE_INFINITY if proposal should not be accepted *
 */
@Override
public double proposal() {
    try {
        double hastingsRatio;
        final double scale = getScaler();
        if (m_bIsTreeScaler) {
            final Tree tree = treeInput.get(this);
            if (rootOnlyInput.get()) {
                final Node root = tree.getRoot();
                final double newHeight = root.getHeight() * scale;
                if (newHeight < Math.max(root.getLeft().getHeight(), root.getRight().getHeight())) {
                    return Double.NEGATIVE_INFINITY;
                }
                root.setHeight(newHeight);
                return -Math.log(scale);
            } else {
                // scale the beast.tree
                final int internalNodes = tree.scale(scale);
                return Math.log(scale) * (internalNodes - 2);
            }
        }
        // not a tree scaler, so scale a parameter
        final boolean scaleAll = scaleAllInput.get();
        final int specifiedDoF = degreesOfFreedomInput.get();
        final boolean scaleAllIndependently = scaleAllIndependentlyInput.get();
        final RealParameter param = parameterInput.get(this);
        assert param.getLower() != null && param.getUpper() != null;
        final int dim = param.getDimension();
        if (scaleAllIndependently) {
            // update all dimensions independently.
            hastingsRatio = 0;
            final BooleanParameter indicators = indicatorInput.get();
            if (indicators != null) {
                final int dimCount = indicators.getDimension();
                final Boolean[] indicator = indicators.getValues();
                final boolean impliedOne = dimCount == (dim - 1);
                for (int i = 0; i < dim; i++) {
                    if ((impliedOne && (i == 0 || indicator[i - 1])) || (!impliedOne && indicator[i])) {
                        final double scaleOne = getScaler();
                        final double newValue = scaleOne * param.getValue(i);
                        hastingsRatio -= Math.log(scaleOne);
                        if (outsideBounds(newValue, param)) {
                            return Double.NEGATIVE_INFINITY;
                        }
                        param.setValue(i, newValue);
                    }
                }
            } else {
                for (int i = 0; i < dim; i++) {
                    final double scaleOne = getScaler();
                    final double newValue = scaleOne * param.getValue(i);
                    hastingsRatio -= Math.log(scaleOne);
                    if (outsideBounds(newValue, param)) {
                        return Double.NEGATIVE_INFINITY;
                    }
                    param.setValue(i, newValue);
                }
            }
        } else if (scaleAll) {
            // update all dimensions
            // hasting ratio is dim-2 times of 1dim case. would be nice to have a reference here
            // for the proof. It is supposed to be somewhere in an Alexei/Nicholes article.
            // all Values assumed independent!
            final int computedDoF = param.scale(scale);
            final int usedDoF = (specifiedDoF > 0) ? specifiedDoF : computedDoF;
            hastingsRatio = (usedDoF - 2) * Math.log(scale);
        } else {
            hastingsRatio = -Math.log(scale);
            // which position to scale
            final int index;
            final BooleanParameter indicators = indicatorInput.get();
            if (indicators != null) {
                final int dimCount = indicators.getDimension();
                final Boolean[] indicator = indicators.getValues();
                final boolean impliedOne = dimCount == (dim - 1);
                // available bit locations. there can be hundreds of them. scan list only once.
                final int[] loc = new int[dimCount + 1];
                int locIndex = 0;
                if (impliedOne) {
                    loc[locIndex] = 0;
                    ++locIndex;
                }
                for (int i = 0; i < dimCount; i++) {
                    if (indicator[i]) {
                        loc[locIndex] = i + (impliedOne ? 1 : 0);
                        ++locIndex;
                    }
                }
                if (locIndex > 0) {
                    final int rand = Randomizer.nextInt(locIndex);
                    index = loc[rand];
                } else {
                    // no active indicators
                    return Double.NEGATIVE_INFINITY;
                }
            } else {
                // any is good
                index = Randomizer.nextInt(dim);
            }
            final double oldValue = param.getValue(index);
            if (oldValue == 0) {
                // Error: parameter has value 0 and cannot be scaled
                return Double.NEGATIVE_INFINITY;
            }
            final double newValue = scale * oldValue;
            if (outsideBounds(newValue, param)) {
                // reject out of bounds scales
                return Double.NEGATIVE_INFINITY;
            }
            param.setValue(index, newValue);
        // provides a hook for subclasses
        // cleanupOperation(newValue, oldValue);
        }
        return hastingsRatio;
    } catch (Exception e) {
        // whatever went wrong, we want to abort this operation...
        return Double.NEGATIVE_INFINITY;
    }
}