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

Example 1 with UniformizedSubstitutionModel

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

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

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

the class StateHistoryTest method testGetLogLikelihood.

public void testGetLogLikelihood() {
    System.out.println("Start of getLogLikelihood test");
    int startingState = 1;
    int endingState = 1;
    double duration = 0.5;
    int iterations = 5000000;
    //        int iterations = 1;
    double[] probs = new double[16];
    baseModel.getTransitionProbabilities(duration, probs);
    double trueProb = probs[startingState * 4 + endingState];
    System.out.println("Tru prob = " + trueProb);
    UniformizedSubstitutionModel uSM = new UniformizedSubstitutionModel(baseModel);
    uSM.setSaveCompleteHistory(true);
    double logProb = Double.NEGATIVE_INFINITY;
    double prob = 0.0;
    double condProb = 0.0;
    for (int i = 0; i < iterations; ++i) {
        uSM.computeCondStatMarkovJumps(startingState, endingState, duration, trueProb);
        StateHistory history = uSM.getStateHistory();
        //            System.out.println(history.getEndingTime() - history.getStartingTime());
        assertEquals(history.getEndingTime() - history.getStartingTime(), duration, 10E-3);
        assertEquals(history.getStartingState(), startingState);
        assertEquals(history.getEndingState(), endingState);
        double logLikelihood = history.getLogLikelihood(lambda, 4);
        prob += Math.exp(logLikelihood);
        logProb = LogTricks.logSum(logProb, logLikelihood);
        condProb += Math.exp(-logLikelihood);
    //            System.err.println(logLikelihood);
    }
    logProb = Math.exp(logProb - Math.log(iterations));
    prob /= iterations;
    condProb /= iterations;
    System.out.println("Sim prob = " + prob);
    System.out.println("Inv prob = " + (1.0 / condProb));
    //        System.out.println("log prob = " + logProb);
    //   System.exit(-1);
    System.out.println();
    System.out.println();
    // Try using unconditioned simulation
    double marginalProb = 0.0;
    double mcProb = 0.0;
    double invMcProb = 0.0;
    int totalTries = 0;
    int i = 0;
    while (i < iterations) {
        startingState = MathUtils.randomChoicePDF(freqModel.getFrequencies());
        StateHistory history = StateHistory.simulateUnconditionalOnEndingState(0, startingState, duration, lambda, 4);
        if (//                    history.getNumberOfJumps() == randomChoice1
        true) {
            marginalProb += 1.0;
            i++;
            double logLike = history.getLogLikelihood(lambda, 4);
            mcProb += Math.exp(logLike);
            invMcProb += Math.exp(-logLike);
        //                if (i % 100000 == 0) System.out.println(i);
        }
        totalTries++;
    }
    marginalProb /= totalTries;
    mcProb /= iterations;
    invMcProb /= iterations;
    System.out.println("Sim uncd = " + marginalProb);
    System.out.println("mc  prob = " + mcProb);
    System.out.println("m2  prob = " + (1.0 / invMcProb));
    assertEquals(prob, trueProb);
}

Example 4 with UniformizedSubstitutionModel

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

the class HiddenMarkovRatesTest method testHiddenRates.

public void testHiddenRates() {
    final int index = 0;
    double[] freqs = getBinaryFreqs(index);
    FrequencyModel binaryFreqModel = new FrequencyModel(TwoStates.INSTANCE, freqs);
    int relativeTo = 0;
    Parameter ratesParameter = new Parameter.Default(0);
    GeneralSubstitutionModel binaryModel = new GeneralSubstitutionModel("binary", TwoStates.INSTANCE, binaryFreqModel, ratesParameter, relativeTo);
    UniformizedSubstitutionModel uSM = new UniformizedSubstitutionModel(binaryModel, MarkovJumpsType.REWARDS);
    uSM.setSaveCompleteHistory(true);
    double[] rewardRegister = new double[] { 0.0, 1.0 };
    uSM.setRegistration(rewardRegister);
    final double[] hkyFreqs = getHKYFreqs(index);
    FrequencyModel hkyFreqModel = new FrequencyModel(Nucleotides.INSTANCE, hkyFreqs);
    final double kappa = getKappa(index);
    final HKY hky = new HKY(kappa, hkyFreqModel);
    final double length = getLength(index);
    double[] resultCompleteHistory = new double[16];
    final int replicates = getNumberReplicates(index);
    double result = 0.0;
    for (int r = 0; r < replicates; ++r) {
        result += oneCompleteHistoryReplicate(resultCompleteHistory, hky, uSM, length);
    }
    result /= replicates;
    normalize(resultCompleteHistory, replicates);
    System.out.println("Averaged probabilities");
    System.out.println(result);
    System.out.println(new Vector(resultCompleteHistory));
    System.out.println();
    double[] intermediate = new double[16];
    hky.getTransitionProbabilities(result, intermediate);
    System.out.println("Intermediate using above average reward");
    System.out.println(result);
    System.out.println(new Vector(intermediate));
    System.out.println();
    double[] resultExpected = new double[16];
    UniformizedSubstitutionModel expectedUSM = new UniformizedSubstitutionModel(binaryModel, MarkovJumpsType.REWARDS, replicates);
    expectedUSM.setRegistration(rewardRegister);
    result = oneCompleteHistoryReplicate(resultExpected, hky, expectedUSM, length);
    System.out.println("Averaged reward");
    System.out.println(result);
    System.out.println(new Vector(resultExpected));
    System.out.println();
    double[] originalProbs = new double[16];
    hky.getTransitionProbabilities(length, originalProbs);
    System.out.println("Original probabilities");
    System.out.println(new Vector(originalProbs));
    System.out.println();
}