Examples with CompleteHistorySimulator dr.app.beagle.tools.CompleteHistorySimulator used on opensource projects

Example 1 with CompleteHistorySimulator

use of dr.app.beagle.tools.CompleteHistorySimulator in project beast-mcmc by beast-dev.

the class CompleteHistorySimulatorParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    int nReplications = xo.getIntegerAttribute(REPLICATIONS);
    Tree tree = (Tree) xo.getChild(Tree.class);
    GammaSiteRateModel siteModel = (GammaSiteRateModel) xo.getChild(GammaSiteRateModel.class);
    BranchRateModel rateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    if (rateModel == null)
        rateModel = new DefaultBranchRateModel();
    DataType dataType = siteModel.getSubstitutionModel().getDataType();
    String jumpTag = xo.getAttribute(JUMP_TAG_NAME, JUMP_TAG);
    boolean sumAcrossSites = xo.getAttribute(SUM_SITES, false);
    Parameter branchSpecificParameter = null;
    Parameter variableValueParameter = null;
    if (xo.hasChildNamed(BRANCH_SPECIFIC_SPECIFICATION)) {
        XMLObject cxo = xo.getChild(BRANCH_SPECIFIC_SPECIFICATION);
        branchSpecificParameter = (Parameter) cxo.getChild(BRANCH_VARIABLE_PARAMETER).getChild(Parameter.class);
        variableValueParameter = (Parameter) cxo.getChild(VARIABLE_VALUE_PARAMETER).getChild(Parameter.class);
    }
    CompleteHistorySimulator history = new CompleteHistorySimulator(tree, siteModel, rateModel, nReplications, sumAcrossSites, branchSpecificParameter, variableValueParameter);
    XMLObject cxo = xo.getChild(COUNTS);
    if (cxo != null) {
        MarkovJumpsTreeLikelihoodParser.parseAllChildren(cxo, history, dataType.getStateCount(), jumpTag, MarkovJumpsType.COUNTS, false);
    }
    cxo = xo.getChild(REWARDS);
    if (cxo != null) {
        MarkovJumpsTreeLikelihoodParser.parseAllChildren(cxo, history, dataType.getStateCount(), jumpTag, MarkovJumpsType.REWARDS, false);
    }
    if (dataType instanceof Codons) {
        Codons codons = (Codons) dataType;
        if (xo.getAttribute(SYN_JUMPS, false)) {
            // use base 61
            double[] synRegMatrix = CodonLabeling.getRegisterMatrix(CodonLabeling.SYN, codons, false);
            Parameter registerParameter = new Parameter.Default(synRegMatrix);
            registerParameter.setId("S");
            history.addRegister(registerParameter, MarkovJumpsType.COUNTS, false);
        }
        if (xo.getAttribute(NON_SYN_JUMPS, false)) {
            // use base 61
            double[] nonSynRegMatrix = CodonLabeling.getRegisterMatrix(CodonLabeling.NON_SYN, codons, false);
            Parameter registerParameter = new Parameter.Default(nonSynRegMatrix);
            registerParameter.setId("N");
            history.addRegister(registerParameter, MarkovJumpsType.COUNTS, false);
        }
    }
    if (xo.getAttribute(ANNOTATE_WITH_ALIGNMENT, false)) {
        history.addAlignmentTrait();
    }
    boolean alignmentOnly = xo.getAttribute(ALIGNMENT_ONLY, false);
    if (dataType instanceof Codons && !alignmentOnly) {
        System.out.println("Codon models give exception when count statistics are done on them. " + "You can supress this by setting alignmentOnly to true.");
    }
    if (alignmentOnly) {
        history.setAlignmentOnly();
    }
    history.simulate();
    return history;
}

Example 2 with CompleteHistorySimulator

use of dr.app.beagle.tools.CompleteHistorySimulator in project beast-mcmc by beast-dev.

the class CompleteHistorySimulatorTest method runSimulation.

protected void runSimulation(int N, Tree tree, GammaSiteRateModel siteModel, BranchRateModel branchRateModel, int nSites, double[][] registers, double analyticResult, Parameter variableParam, Parameter valuesParam) {
    CompleteHistorySimulator simulator = new CompleteHistorySimulator(tree, siteModel, branchRateModel, nSites, false, variableParam, valuesParam);
    for (int r = 0; r < registers.length; r++) {
        Parameter registerParameter = new Parameter.Default(registers[r]);
        registerParameter.setId("set" + (r + 1));
        simulator.addRegister(registerParameter, MarkovJumpsType.COUNTS, false);
    }
    int nJumps = simulator.getNumberOfJumpProcess();
    double[] results = new double[nJumps];
    for (int rep = 0; rep < N; rep++) {
        // Generate new history
        simulator.simulate();
        for (int jump = 0; jump < nJumps; jump++) {
            results[jump] += accumulateJumpsOverTree(tree, simulator, jump);
        }
    }
    for (int jump = 0; jump < nJumps; jump++) {
        results[jump] /= (double) N * (double) nSites;
    }
    if (nJumps > 0) {
        System.out.println("simulations = " + new Vector(results));
    } else {
        System.out.println("No jump processes!");
    }
    double simulationResult = accumulate(results);
    assertEquals(analyticResult, simulationResult, 1E-2);
}