Examples with MarkovJumpsSubstitutionModel dr.evomodel.substmodel.MarkovJumpsSubstitutionModel used on opensource projects

Example 1 with MarkovJumpsSubstitutionModel

use of dr.evomodel.substmodel.MarkovJumpsSubstitutionModel 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) {
        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 (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;
                    }
                });
            }
        }
        numRegisters++;
    }
// End of loop over branch models
}

Example 2 with MarkovJumpsSubstitutionModel

use of dr.evomodel.substmodel.MarkovJumpsSubstitutionModel 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, 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 MarkovJumpsSubstitutionModel

use of dr.evomodel.substmodel.MarkovJumpsSubstitutionModel in project beast-mcmc by beast-dev.

the class UniformizedStateHistoryTest method testStateHistorySimulationForJumps.

public void testStateHistorySimulationForJumps() {
    try {
        double startingTime = 1.0;
        double endingTime = 3.0;
        int startingState = 1;
        int endingState = 3;
        int N = 1000000;
        double[] tmp = new double[stateCount * stateCount];
        hky.getTransitionProbabilities(endingTime - startingTime, tmp);
        double transitionProbability = tmp[startingState * stateCount + endingState];
        double[][] registers = new double[2][stateCount * stateCount];
        // Count all jumps
        MarkovJumpsCore.fillRegistrationMatrix(registers[0], stateCount);
        // Mark just one state!
        registers[1][2 * stateCount + 1] = 1.0;
        double[] expectations = new double[registers.length];
        for (int i = 0; i < N; i++) {
            StateHistory history = UniformizedStateHistory.simulateConditionalOnEndingState(startingTime, startingState, endingTime, endingState, transitionProbability, stateCount, process);
            for (int j = 0; j < registers.length; j++) {
                expectations[j] += history.getTotalRegisteredCounts(registers[j]);
            }
        }
        // Determine analytic solution
        MarkovJumpsSubstitutionModel markovjumps = new MarkovJumpsSubstitutionModel(hky);
        double[] mjExpectations = new double[stateCount * stateCount];
        for (int j = 0; j < registers.length; j++) {
            expectations[j] /= (double) N;
            System.out.println("Expected number for register = " + expectations[j]);
            markovjumps.setRegistration(registers[j]);
            markovjumps.computeCondStatMarkovJumps(endingTime - startingTime, mjExpectations);
            assertEquals(mjExpectations[startingState * stateCount + endingState], expectations[j], 1E-2);
        }
    } catch (SubordinatedProcess.Exception e) {
        throw new RuntimeException("Subordinated process exception");
    }
}

Example 4 with MarkovJumpsSubstitutionModel

use of dr.evomodel.substmodel.MarkovJumpsSubstitutionModel in project beast-mcmc by beast-dev.

the class StateHistoryTest method setUp.

public void setUp() {
    MathUtils.setSeed(666);
    freqModel = new FrequencyModel(Nucleotides.INSTANCE, new double[] { 0.45, 0.25, 0.05, 0.25 });
    baseModel = new HKY(2.0, freqModel);
    stateCount = baseModel.getDataType().getStateCount();
    lambda = new double[stateCount * stateCount];
    baseModel.getInfinitesimalMatrix(lambda);
    System.out.println("lambda = " + new Vector(lambda));
    markovjumps = new MarkovJumpsSubstitutionModel(baseModel);
}

Example 5 with MarkovJumpsSubstitutionModel

use of dr.evomodel.substmodel.MarkovJumpsSubstitutionModel in project beast-mcmc by beast-dev.

the class MarkovJumpsSubstitutionModelTest method testMarginalRates.

public void testMarginalRates() {
    HKY substModel = new HKY(2.0, new FrequencyModel(Nucleotides.INSTANCE, // A,C,G,T
    new double[] { 0.3, 0.2, 0.25, 0.25 }));
    int states = substModel.getDataType().getStateCount();
    MarkovJumpsSubstitutionModel markovjumps = new MarkovJumpsSubstitutionModel(substModel, MarkovJumpsType.COUNTS);
    double[] r = new double[states * states];
    // A
    int from = 0;
    // C
    int to = 1;
    MarkovJumpsCore.fillRegistrationMatrix(r, from, to, states, 1.0);
    markovjumps.setRegistration(r);
    double marginalRate = markovjumps.getMarginalRate();
    System.out.println("Marginal rate = " + marginalRate);
    assertEquals(rMarkovMarginalRate, marginalRate, tolerance);
    MarkovJumpsCore.fillRegistrationMatrix(r, states);
    markovjumps.setRegistration(r);
    marginalRate = markovjumps.getMarginalRate();
    System.out.println("Marginal rate = " + marginalRate);
    assertEquals(1.0, marginalRate, tolerance);
}