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

Example 6 with Likelihood

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

the class BirthDeathLikelihoodTest method birthDeathLikelihoodTester.

private void birthDeathLikelihoodTester(Tree tree, double birthRate, double deathRate, double logL) {
    Parameter b = new Parameter.Default("b", birthRate, 0.0, Double.MAX_VALUE);
    Parameter d = new Parameter.Default("d", deathRate, 0.0, Double.MAX_VALUE);
    SpeciationModel speciationModel = new BirthDeathGernhard08Model(b, d, null, BirthDeathGernhard08Model.TreeType.ORIENTED, Units.Type.YEARS);
    Likelihood likelihood = new SpeciationLikelihood(tree, speciationModel, "bd.like");
    assertEquals(logL, likelihood.getLogLikelihood(), 1e-14);
}

Example 7 with Likelihood

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

the class BirthDeathSSLikelihoodTest method likelihoodTester.

private void likelihoodTester(Tree tree, double birthRate, double deathRate, Variable<Double> origin, double logL) {
    Variable<Double> b = new Variable.D("b", birthRate);
    Variable<Double> d = new Variable.D("d", deathRate);
    Variable<Double> psi = new Variable.D("psi", this.psi);
    Variable<Double> p = new Variable.D("p", this.p);
    Variable<Double> r = new Variable.D("r", 0.5);
    Variable<Double> fTime = new Variable.D("time", 0.0);
    SpeciationModel speciationModel = new BirthDeathSerialSamplingModel(b, d, psi, p, false, r, true, origin, Units.Type.YEARS);
    Likelihood likelihood = new SpeciationLikelihood(tree, speciationModel, "bdss.like");
    assertEquals(logL, likelihood.getLogLikelihood());
}

Example 8 with Likelihood

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

the class TreeDataLikelihoodParser method createTreeDataLikelihood.

protected Likelihood createTreeDataLikelihood(List<PatternList> patternLists, List<BranchModel> branchModels, List<SiteRateModel> siteRateModels, TreeModel treeModel, BranchRateModel branchRateModel, TipStatesModel tipStatesModel, boolean useAmbiguities, PartialsRescalingScheme scalingScheme, boolean delayRescalingUntilUnderflow) throws XMLParseException {
    if (tipStatesModel != null) {
        throw new XMLParseException("Tip State Error models are not supported yet with TreeDataLikelihood");
    }
    List<Taxon> treeTaxa = treeModel.asList();
    List<Taxon> patternTaxa = patternLists.get(0).asList();
    if (!patternTaxa.containsAll(treeTaxa)) {
        throw new XMLParseException("TreeModel contains more taxa than the partition pattern list.");
    }
    if (!treeTaxa.containsAll(patternTaxa)) {
        throw new XMLParseException("TreeModel contains fewer taxa than the partition pattern list.");
    }
    //        DataLikelihoodDelegate dataLikelihoodDelegate = new BeagleDataLikelihoodDelegate(
    //                treeModel,
    //                patternLists.get(0),
    //                branchModel,
    //                siteRateModel,
    //                useAmbiguities,
    //                scalingScheme,
    //                delayRescalingUntilUnderflow);
    boolean useBeagle3 = Boolean.parseBoolean(System.getProperty("USE_BEAGLE3", "true"));
    boolean useJava = Boolean.parseBoolean(System.getProperty("java.only", "false"));
    if (useBeagle3 && MultiPartitionDataLikelihoodDelegate.IS_MULTI_PARTITION_COMPATIBLE() && !useJava) {
        DataLikelihoodDelegate dataLikelihoodDelegate = new MultiPartitionDataLikelihoodDelegate(treeModel, patternLists, branchModels, siteRateModels, useAmbiguities, scalingScheme, delayRescalingUntilUnderflow);
        return new TreeDataLikelihood(dataLikelihoodDelegate, treeModel, branchRateModel);
    } else {
        // The multipartition data likelihood isn't available so make a set of single partition data likelihoods
        List<Likelihood> treeDataLikelihoods = new ArrayList<Likelihood>();
        for (int i = 0; i < patternLists.size(); i++) {
            DataLikelihoodDelegate dataLikelihoodDelegate = new BeagleDataLikelihoodDelegate(treeModel, patternLists.get(i), branchModels.get(i), siteRateModels.get(i), useAmbiguities, scalingScheme, delayRescalingUntilUnderflow);
            treeDataLikelihoods.add(new TreeDataLikelihood(dataLikelihoodDelegate, treeModel, branchRateModel));
        }
        if (treeDataLikelihoods.size() == 1) {
            return treeDataLikelihoods.get(0);
        }
        return new CompoundLikelihood(treeDataLikelihoods);
    }
}

Example 9 with Likelihood

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

the class BalancedBeagleTreeLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    boolean useAmbiguities = xo.getAttribute(BeagleTreeLikelihoodParser.USE_AMBIGUITIES, false);
    /*int instanceCount = xo.getAttribute(INSTANCE_COUNT, 1);
        if (instanceCount < 1) {
            instanceCount = 1;
        }

        String ic = System.getProperty(BEAGLE_INSTANCE_COUNT);
        if (ic != null && ic.length() > 0) {
            instanceCount = Integer.parseInt(ic);
        }*/
    PatternList patternList = (PatternList) xo.getChild(PatternList.class);
    TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
    GammaSiteRateModel siteRateModel = (GammaSiteRateModel) xo.getChild(GammaSiteRateModel.class);
    FrequencyModel rootFreqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
    BranchModel branchModel = (BranchModel) xo.getChild(BranchModel.class);
    if (branchModel == null) {
        SubstitutionModel substitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
        if (substitutionModel == null) {
            substitutionModel = siteRateModel.getSubstitutionModel();
        }
        if (substitutionModel == null) {
            throw new XMLParseException("No substitution model available for TreeLikelihood: " + xo.getId());
        }
        branchModel = new HomogeneousBranchModel(substitutionModel, rootFreqModel);
    }
    BranchRateModel branchRateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    TipStatesModel tipStatesModel = (TipStatesModel) xo.getChild(TipStatesModel.class);
    //        if (xo.getChild(TipStatesModel.class) != null) {
    //            throw new XMLParseException("Sequence Error Models are not supported under BEAGLE yet. Please use Native BEAST Likelihood.");
    //        }
    PartialsRescalingScheme scalingScheme = PartialsRescalingScheme.DEFAULT;
    if (xo.hasAttribute(BeagleTreeLikelihoodParser.SCALING_SCHEME)) {
        //            scalingScheme = PartialsRescalingScheme.parseFromString(xo.getStringAttribute(BeagleTreeLikelihoodParser.SCALING_SCHEME));
        if (scalingScheme == null)
            throw new XMLParseException("Unknown scaling scheme '" + xo.getStringAttribute(BeagleTreeLikelihoodParser.SCALING_SCHEME) + "' in " + "OldBeagleTreeLikelihood object '" + xo.getId());
    }
    boolean delayScaling = true;
    Map<Set<String>, Parameter> partialsRestrictions = null;
    if (xo.hasChildNamed(PARTIALS_RESTRICTION)) {
        XMLObject cxo = xo.getChild(PARTIALS_RESTRICTION);
        TaxonList taxonList = (TaxonList) cxo.getChild(TaxonList.class);
        //            Parameter parameter = (Parameter) cxo.getChild(Parameter.class);
        try {
            TreeUtils.getLeavesForTaxa(treeModel, taxonList);
        } catch (TreeUtils.MissingTaxonException e) {
            throw new XMLParseException("Unable to parse taxon list: " + e.getMessage());
        }
        throw new XMLParseException("Restricting internal nodes is not yet implemented.  Contact Marc");
    }
    /*if (instanceCount == 1 || patternList.getPatternCount() < instanceCount) {
            return createTreeLikelihood(
                    patternList,
                    treeModel,
                    branchModel,
                    siteRateModel,
                    branchRateModel,
                    tipStatesModel,
                    useAmbiguities,
                    scalingScheme,
                    partialsRestrictions,
                    xo
            );
        }*/
    //first run a test for instanceCount == 1
    System.err.println("\nTesting instanceCount == 1");
    Likelihood baseLikelihood = createTreeLikelihood(patternList, treeModel, branchModel, siteRateModel, branchRateModel, tipStatesModel, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, xo);
    double start = System.nanoTime();
    for (int i = 0; i < TEST_RUNS; i++) {
        baseLikelihood.makeDirty();
        baseLikelihood.getLogLikelihood();
    }
    double end = System.nanoTime();
    double baseResult = end - start;
    System.err.println("Evaluation took: " + baseResult);
    if (!(patternList instanceof SitePatterns)) {
        throw new XMLParseException("BEAGLE_INSTANCES option cannot be used with BEAUti-selected codon partitioning.");
    }
    if (tipStatesModel != null) {
        throw new XMLParseException("BEAGLE_INSTANCES option cannot be used with a TipStateModel (i.e., a sequence error model).");
    }
    //List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    List<Likelihood> likelihoods = null;
    CompoundLikelihood compound = null;
    int instanceCount = 2;
    boolean optimal = false;
    while (optimal == false) {
        System.err.println("\nCreating instanceCount == " + instanceCount);
        likelihoods = new ArrayList<Likelihood>();
        for (int i = 0; i < instanceCount; i++) {
            Patterns subPatterns = new Patterns((SitePatterns) patternList, 0, 0, 1, i, instanceCount);
            AbstractTreeLikelihood treeLikelihood = createTreeLikelihood(subPatterns, treeModel, branchModel, siteRateModel, branchRateModel, null, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, xo);
            treeLikelihood.setId(xo.getId() + "_" + instanceCount);
            likelihoods.add(treeLikelihood);
        }
        //construct compoundLikelihood
        compound = new CompoundLikelihood(instanceCount, likelihoods);
        //test timings 
        System.err.println("\nTesting instanceCount == " + instanceCount);
        start = System.nanoTime();
        for (int i = 0; i < TEST_RUNS; i++) {
            compound.makeDirty();
            compound.getLogLikelihood();
        }
        end = System.nanoTime();
        double newResult = end - start;
        System.err.println("Evaluation took: " + newResult);
        if (baseResult / newResult > TEST_CUTOFF) {
            instanceCount++;
            baseResult = newResult;
        } else {
            optimal = true;
            instanceCount--;
            System.err.println("\nCreating final BeagleTreeLikelihood with instanceCount: " + instanceCount);
            likelihoods = new ArrayList<Likelihood>();
            for (int i = 0; i < instanceCount; i++) {
                Patterns subPatterns = new Patterns((SitePatterns) patternList, 0, 0, 1, i, instanceCount);
                AbstractTreeLikelihood treeLikelihood = createTreeLikelihood(subPatterns, treeModel, branchModel, siteRateModel, branchRateModel, null, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, xo);
                treeLikelihood.setId(xo.getId() + "_" + instanceCount);
                likelihoods.add(treeLikelihood);
            }
            //construct compoundLikelihood
            compound = new CompoundLikelihood(instanceCount, likelihoods);
        }
    }
    return compound;
/*for (int i = 0; i < instanceCount; i++) {

            Patterns subPatterns = new Patterns((SitePatterns)patternList, 0, 0, 1, i, instanceCount);

            AbstractTreeLikelihood treeLikelihood = createTreeLikelihood(
                    subPatterns,
                    treeModel,
                    branchModel,
                    siteRateModel,
                    branchRateModel,
                    null,
                    useAmbiguities,
                    scalingScheme,
                    partialsRestrictions,
                    xo);
            treeLikelihood.setId(xo.getId() + "_" + instanceCount);
            likelihoods.add(treeLikelihood);
        }

        return new CompoundLikelihood(likelihoods);*/
}

Example 10 with Likelihood

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

the class BeagleTreeLikelihoodParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    boolean useAmbiguities = xo.getAttribute(USE_AMBIGUITIES, false);
    int instanceCount = xo.getAttribute(INSTANCE_COUNT, 1);
    if (instanceCount < 1) {
        instanceCount = 1;
    }
    String ic = System.getProperty(BEAGLE_INSTANCE_COUNT);
    if (ic != null && ic.length() > 0) {
        instanceCount = Integer.parseInt(ic);
    }
    PatternList patternList = (PatternList) xo.getChild(PatternList.class);
    TreeModel treeModel = (TreeModel) xo.getChild(TreeModel.class);
    GammaSiteRateModel siteRateModel = (GammaSiteRateModel) xo.getChild(GammaSiteRateModel.class);
    FrequencyModel rootFreqModel = (FrequencyModel) xo.getChild(FrequencyModel.class);
    BranchModel branchModel = (BranchModel) xo.getChild(BranchModel.class);
    if (branchModel == null) {
        SubstitutionModel substitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
        if (substitutionModel == null) {
            substitutionModel = siteRateModel.getSubstitutionModel();
        }
        if (substitutionModel == null) {
            throw new XMLParseException("No substitution model available for TreeLikelihood: " + xo.getId());
        }
        branchModel = new HomogeneousBranchModel(substitutionModel, rootFreqModel);
    }
    BranchRateModel branchRateModel = (BranchRateModel) xo.getChild(BranchRateModel.class);
    TipStatesModel tipStatesModel = (TipStatesModel) xo.getChild(TipStatesModel.class);
    //        if (xo.getChild(TipStatesModel.class) != null) {
    //            throw new XMLParseException("Sequence Error Models are not supported under BEAGLE yet. Please use Native BEAST Likelihood.");
    //        }
    PartialsRescalingScheme scalingScheme = PartialsRescalingScheme.DEFAULT;
    boolean delayScaling = true;
    if (xo.hasAttribute(SCALING_SCHEME)) {
        scalingScheme = PartialsRescalingScheme.parseFromString(xo.getStringAttribute(SCALING_SCHEME));
        if (scalingScheme == null)
            throw new XMLParseException("Unknown scaling scheme '" + xo.getStringAttribute(SCALING_SCHEME) + "' in " + "OldBeagleTreeLikelihood object '" + xo.getId());
    }
    if (xo.hasAttribute(DELAY_SCALING)) {
        delayScaling = xo.getBooleanAttribute(DELAY_SCALING);
    }
    Map<Set<String>, Parameter> partialsRestrictions = null;
    if (xo.hasChildNamed(PARTIALS_RESTRICTION)) {
        XMLObject cxo = xo.getChild(PARTIALS_RESTRICTION);
        TaxonList taxonList = (TaxonList) cxo.getChild(TaxonList.class);
        //            Parameter parameter = (Parameter) cxo.getChild(Parameter.class);
        try {
            TreeUtils.getLeavesForTaxa(treeModel, taxonList);
        } catch (TreeUtils.MissingTaxonException e) {
            throw new XMLParseException("Unable to parse taxon list: " + e.getMessage());
        }
        throw new XMLParseException("Restricting internal nodes is not yet implemented.  Contact Marc");
    }
    if (instanceCount == 1 || patternList.getPatternCount() < instanceCount) {
        return createTreeLikelihood(patternList, treeModel, branchModel, siteRateModel, branchRateModel, tipStatesModel, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, xo);
    }
    if (tipStatesModel != null) {
        throw new XMLParseException("BEAGLE_INSTANCES option cannot be used with a TipStateModel (i.e., a sequence error model).");
    }
    List<Likelihood> likelihoods = new ArrayList<Likelihood>();
    for (int i = 0; i < instanceCount; i++) {
        Patterns subPatterns = new Patterns(patternList, i, instanceCount);
        AbstractTreeLikelihood treeLikelihood = createTreeLikelihood(subPatterns, treeModel, branchModel, siteRateModel, branchRateModel, null, useAmbiguities, scalingScheme, delayScaling, partialsRestrictions, xo);
        treeLikelihood.setId(xo.getId() + "_" + instanceCount);
        likelihoods.add(treeLikelihood);
    }
    return new CompoundLikelihood(likelihoods);
}