Examples with EpochBranchModel dr.evomodel.branchmodel.EpochBranchModel used on opensource projects

Example 1 with EpochBranchModel

use of dr.evomodel.branchmodel.EpochBranchModel in project beast-mcmc by beast-dev.

the class BeagleTreeLikelihood method main.

public static void main(String[] args) {
    try {
        MathUtils.setSeed(666);
        System.out.println("Test case 1: simulateOnePartition");
        int sequenceLength = 1000;
        ArrayList<Partition> partitionsList = new ArrayList<Partition>();
        // create tree
        NewickImporter importer = new NewickImporter("(SimSeq1:73.7468,(SimSeq2:25.256989999999995,SimSeq3:45.256989999999995):18.48981);");
        Tree tree = importer.importTree(null);
        TreeModel treeModel = new TreeModel(tree);
        // create Frequency Model
        Parameter freqs = new Parameter.Default(new double[] { 0.25, 0.25, 0.25, 0.25 });
        FrequencyModel freqModel = new FrequencyModel(Nucleotides.INSTANCE, freqs);
        // create branch model
        Parameter kappa1 = new Parameter.Default(1, 1);
        Parameter kappa2 = new Parameter.Default(1, 1);
        HKY hky1 = new HKY(kappa1, freqModel);
        HKY hky2 = new HKY(kappa2, freqModel);
        HomogeneousBranchModel homogenousBranchSubstitutionModel = new HomogeneousBranchModel(hky1);
        List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
        substitutionModels.add(hky1);
        substitutionModels.add(hky2);
        List<FrequencyModel> freqModels = new ArrayList<FrequencyModel>();
        freqModels.add(freqModel);
        Parameter epochTimes = new Parameter.Default(1, 20);
        // create branch rate model
        Parameter rate = new Parameter.Default(1, 0.001);
        BranchRateModel branchRateModel = new StrictClockBranchRates(rate);
        // create site model
        GammaSiteRateModel siteRateModel = new GammaSiteRateModel("siteModel");
        BranchModel homogeneousBranchModel = new HomogeneousBranchModel(hky1);
        BranchModel epochBranchModel = new EpochBranchModel(treeModel, substitutionModels, epochTimes);
        // create partition
        Partition partition1 = new //
        Partition(//
        treeModel, //
        homogenousBranchSubstitutionModel, //
        siteRateModel, //
        branchRateModel, //
        freqModel, // from
        0, // to
        sequenceLength - 1, // every
        1);
        partitionsList.add(partition1);
        // feed to sequence simulator and generate data
        BeagleSequenceSimulator simulator = new BeagleSequenceSimulator(partitionsList);
        Alignment alignment = simulator.simulate(false, false);
        BeagleTreeLikelihood nbtl = new BeagleTreeLikelihood(alignment, treeModel, homogeneousBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
        System.out.println("nBTL(homogeneous) = " + nbtl.getLogLikelihood());
        nbtl = new BeagleTreeLikelihood(alignment, treeModel, epochBranchModel, siteRateModel, branchRateModel, null, false, PartialsRescalingScheme.DEFAULT, false);
        System.out.println("nBTL(epoch) = " + nbtl.getLogLikelihood());
    } catch (Exception e) {
        e.printStackTrace();
        System.exit(-1);
    }
// END: try-catch block
}

Example 2 with EpochBranchModel

use of dr.evomodel.branchmodel.EpochBranchModel in project beast-mcmc by beast-dev.

the class MarkovJumpsBeagleTreeLikelihood method hookCalculation.

protected void hookCalculation(Tree tree, NodeRef parentNode, NodeRef childNode, int[] parentStates, int[] childStates, double[] inProbabilities, int[] rateCategory) {
    final int childNum = childNode.getNumber();
    double[] probabilities = inProbabilities;
    if (probabilities == null) {
        // Leaf will call this hook with a null
        getMatrix(childNum, tmpProbabilities);
        probabilities = tmpProbabilities;
    }
    final double branchRate = branchRateModel.getBranchRate(tree, childNode);
    final double parentTime = tree.getNodeHeight(parentNode);
    final double childTime = tree.getNodeHeight(childNode);
    final double substTime = parentTime - childTime;
    for (int r = 0; r < markovjumps.size(); r++) {
        MarkovJumpsSubstitutionModel thisMarkovJumps = markovjumps.get(r);
        final int modelNumberFromrRegistry = branchModelNumber.get(r);
        // int dummy = 0;
        // final int modelNumberFromTree = branchSubstitutionModel.getBranchIndex(tree, childNode, dummy);
        // @todo AR - not sure about this - if this is an epoch this is just going to get the most
        // @todo tipward model for the branch. I think this was what was happening before (in comment,
        // @todo above).
        BranchModel.Mapping mapping = branchModel.getBranchModelMapping(childNode);
        if (modelNumberFromrRegistry == mapping.getOrder()[0]) {
            if (useUniformization) {
                computeSampledMarkovJumpsForBranch(((UniformizedSubstitutionModel) thisMarkovJumps), substTime, branchRate, childNum, parentStates, childStates, parentTime, childTime, probabilities, scaleByTime[r], expectedJumps.get(r), rateCategory, (branchModel instanceof EpochBranchModel) || r == historyRegisterNumber);
            } else {
                computeIntegratedMarkovJumpsForBranch(thisMarkovJumps, substTime, branchRate, childNum, parentStates, childStates, probabilities, condJumps, scaleByTime[r], expectedJumps.get(r), rateCategory);
            }
        } else {
            // Fill with zeros
            double[] result = expectedJumps.get(r)[childNum];
            Arrays.fill(result, 0.0);
        }
    }
}

Example 3 with EpochBranchModel

use of dr.evomodel.branchmodel.EpochBranchModel in project beast-mcmc by beast-dev.

the class MarkovJumpsBeagleTreeLikelihood method addRegister.

public void addRegister(Parameter addRegisterParameter, MarkovJumpsType type, boolean scaleByTime) {
    if ((type == MarkovJumpsType.COUNTS && addRegisterParameter.getDimension() != stateCount * stateCount) || (type == MarkovJumpsType.REWARDS && addRegisterParameter.getDimension() != stateCount)) {
        throw new RuntimeException("Register parameter of wrong dimension");
    }
    addVariable(addRegisterParameter);
    final String tag = addRegisterParameter.getId();
    for (int i = 0; i < substitutionModelDelegate.getSubstitutionModelCount(); ++i) {
        boolean isEpochModel = branchModel instanceof EpochBranchModel;
        registerParameter.add(addRegisterParameter);
        MarkovJumpsSubstitutionModel mjModel;
        SubstitutionModel substitutionModel = substitutionModelDelegate.getSubstitutionModel(i);
        if (useUniformization) {
            mjModel = new UniformizedSubstitutionModel(substitutionModel, type, nSimulants);
        } else {
            if (type == MarkovJumpsType.HISTORY) {
                throw new RuntimeException("Can only report complete history using uniformization");
            }
            mjModel = new MarkovJumpsSubstitutionModel(substitutionModel, type);
        }
        markovjumps.add(mjModel);
        branchModelNumber.add(i);
        addModel(mjModel);
        setupRegistration(numRegisters);
        String traitName;
        if (substitutionModelDelegate.getSubstitutionModelCount() == 1) {
            traitName = tag;
        } else {
            traitName = tag + i;
        }
        jumpTag.add(traitName);
        expectedJumps.add(new double[treeModel.getNodeCount()][patternCount]);
        // storedExpectedJumps.add(new double[treeModel.getNodeCount()][patternCount]);
        boolean[] oldScaleByTime = this.scaleByTime;
        int oldScaleByTimeLength = (oldScaleByTime == null ? 0 : oldScaleByTime.length);
        this.scaleByTime = new boolean[oldScaleByTimeLength + 1];
        if (oldScaleByTimeLength > 0) {
            System.arraycopy(oldScaleByTime, 0, this.scaleByTime, 0, oldScaleByTimeLength);
        }
        this.scaleByTime[oldScaleByTimeLength] = scaleByTime;
        if (type != MarkovJumpsType.HISTORY) {
            TreeTrait.DA da = new TreeTrait.DA() {

                final int registerNumber = numRegisters;

                public String getTraitName() {
                    return tag;
                }

                public Intent getIntent() {
                    return Intent.BRANCH;
                }

                public double[] getTrait(Tree tree, NodeRef node) {
                    return getMarkovJumpsForNodeAndRegister(tree, node, registerNumber);
                }
            };
            treeTraits.addTrait(traitName + "_base", da);
            treeTraits.addTrait(addRegisterParameter.getId(), new TreeTrait.SumAcrossArrayD(new TreeTrait.SumOverTreeDA(da)));
        } else {
            if (i == 0 || !isEpochModel) {
                if (histories == null) {
                    histories = new String[treeModel.getNodeCount()][patternCount];
                } else {
                    throw new RuntimeException("Only one complete history per markovJumpTreeLikelihood is allowed");
                }
                if (nSimulants > 1) {
                    throw new RuntimeException("Only one simulant allowed when saving complete history");
                }
                // Add total number of changes over tree trait
                TreeTrait da = new TreeTrait.DA() {

                    final int registerNumber = numRegisters;

                    public String getTraitName() {
                        return tag;
                    }

                    public Intent getIntent() {
                        return Intent.BRANCH;
                    }

                    public double[] getTrait(Tree tree, NodeRef node) {
                        return getMarkovJumpsForNodeAndRegister(tree, node, registerNumber);
                    }
                };
                treeTraits.addTrait(addRegisterParameter.getId(), new TreeTrait.SumOverTreeDA(da));
                // Record the complete history for this register
                historyRegisterNumber = numRegisters;
                ((UniformizedSubstitutionModel) mjModel).setSaveCompleteHistory(true);
                if (useCompactHistory && logHistory) {
                    treeTraits.addTrait(ALL_HISTORY, new TreeTrait.SA() {

                        public String getTraitName() {
                            return ALL_HISTORY;
                        }

                        public Intent getIntent() {
                            return Intent.BRANCH;
                        }

                        public boolean getFormatAsArray() {
                            return true;
                        }

                        public String[] getTrait(Tree tree, NodeRef node) {
                            List<String> events = new ArrayList<String>();
                            for (int i = 0; i < patternCount; i++) {
                                String eventString = getHistoryForNode(tree, node, i);
                                if (eventString != null && eventString.compareTo("{}") != 0) {
                                    eventString = eventString.substring(1, eventString.length() - 1);
                                    if (eventString.contains("},{")) {
                                        // There are multiple events
                                        String[] elements = eventString.split("(?<=\\}),(?=\\{)");
                                        for (String e : elements) {
                                            events.add(e);
                                        }
                                    } else {
                                        events.add(eventString);
                                    }
                                }
                            }
                            String[] array = new String[events.size()];
                            events.toArray(array);
                            return array;
                        }

                        public boolean getLoggable() {
                            return true;
                        }
                    });
                }
                for (int site = 0; site < patternCount; ++site) {
                    final String anonName = (patternCount == 1) ? HISTORY : HISTORY + "_" + (site + 1);
                    final int anonSite = site;
                    treeTraits.addTrait(anonName, new TreeTrait.S() {

                        public String getTraitName() {
                            return anonName;
                        }

                        public Intent getIntent() {
                            return Intent.BRANCH;
                        }

                        public String getTrait(Tree tree, NodeRef node) {
                            String history = getHistoryForNode(tree, node, anonSite);
                            // Return null if empty
                            return (history.compareTo("{}") != 0) ? history : null;
                        }

                        public boolean getLoggable() {
                            return logHistory && !useCompactHistory;
                        }
                    });
                }
            }
            if (isEpochModel) {
                for (int j = 0; j < markovjumps.size(); ++j) {
                    ((UniformizedSubstitutionModel) markovjumps.get(j)).setSaveCompleteHistory(true);
                }
            }
        }
        numRegisters++;
    }
// End of loop over branch models
}

Example 4 with EpochBranchModel

use of dr.evomodel.branchmodel.EpochBranchModel 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);
    MutableTreeModel treeModel = (MutableTreeModel) xo.getChild(MutableTreeModel.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);
    if (branchModel instanceof EpochBranchModel && rootFreqModel != null) {
        EpochBranchModel epochBranchModel = (EpochBranchModel) branchModel;
        epochBranchModel.setRootFrequencyModel(rootFreqModel);
    }
    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");
    }
    int beagleThreadCount = -1;
    if (System.getProperty(BEAGLE_THREAD_COUNT) != null) {
        beagleThreadCount = Integer.parseInt(System.getProperty(BEAGLE_THREAD_COUNT));
    }
    if (beagleThreadCount == -1) {
        // the default is -1 threads (automatic thread pool size) but an XML attribute can override it
        int threadCount = xo.getAttribute(THREADS, -1);
        if (System.getProperty(THREAD_COUNT) != null) {
            threadCount = Integer.parseInt(System.getProperty(THREAD_COUNT));
        }
        // Todo: allow for different number of threads per beagle instance according to pattern counts
        if (threadCount >= 0) {
            System.setProperty(BEAGLE_THREAD_COUNT, Integer.toString(threadCount / instanceCount));
        }
    }
    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);
}

Example 5 with EpochBranchModel

use of dr.evomodel.branchmodel.EpochBranchModel in project beast-mcmc by beast-dev.

the class EpochBranchModelParser method parseXMLObject.

public Object parseXMLObject(XMLObject xo) throws XMLParseException {
    Logger.getLogger("dr.evomodel").info("\nUsing multi-epoch branch model.");
    MutableTreeModel treeModel = (MutableTreeModel) xo.getChild(MutableTreeModel.class);
    SubstitutionModel ancestralSubstitutionModel = (SubstitutionModel) xo.getChild(SubstitutionModel.class);
    List<Epoch> epochs = new ArrayList<Epoch>();
    for (int i = 0; i < xo.getChildCount(); i++) {
        if (xo.getChild(i) instanceof XMLObject) {
            XMLObject xoc = (XMLObject) xo.getChild(i);
            if (xoc.getName().equals(EPOCH)) {
                Parameter tt = null;
                if (xoc.hasAttribute(TRANSITION_TIME)) {
                    double t = xoc.getAttribute(TRANSITION_TIME, 0.0);
                    tt = new Parameter.Default(1, t);
                }
                SubstitutionModel s = (SubstitutionModel) xoc.getChild(SubstitutionModel.class);
                if (xoc.hasChildNamed(TRANSITION_TIME)) {
                    if (tt != null) {
                        throw new XMLParseException("An epoch cannot have a transitionTime attribute and a parameter");
                    }
                    tt = (Parameter) xoc.getElementFirstChild(TRANSITION_TIME);
                }
                epochs.add(new Epoch(s, tt));
            }
        }
    }
    Collections.sort(epochs);
    List<SubstitutionModel> substitutionModels = new ArrayList<SubstitutionModel>();
    CompoundParameter transitionTimes = new CompoundParameter("epochTimes");
    for (Epoch epoch : epochs) {
        substitutionModels.add(epoch.substitutionModel);
        transitionTimes.addParameter(epoch.timeParameter);
    }
    substitutionModels.add(ancestralSubstitutionModel);
    return new EpochBranchModel(treeModel, substitutionModels, transitionTimes);
}