From b7d3b2b0f3ea52130944e7129fd3a3d28e0e789a Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Fri, 31 Jan 2025 11:03:26 -0800 Subject: [PATCH 01/21] adds MarkovRewardDensityCalculator for comparing implementations --- .../tools/MarkovRewardDensityCalculator.java | 162 ++++++++++++++++++ 1 file changed, 162 insertions(+) create mode 100644 src/dr/app/tools/MarkovRewardDensityCalculator.java diff --git a/src/dr/app/tools/MarkovRewardDensityCalculator.java b/src/dr/app/tools/MarkovRewardDensityCalculator.java new file mode 100644 index 0000000000..d96088bfb9 --- /dev/null +++ b/src/dr/app/tools/MarkovRewardDensityCalculator.java @@ -0,0 +1,162 @@ +package dr.app.tools; + +import java.io.IOException; + +import dr.app.util.Arguments; +import dr.inference.markovjumps.SericolaSeriesMarkovReward; +import dr.inference.markovjumps.TwoStateOccupancyMarkovReward; +import dr.inference.markovjumps.TwoStateSericolaSeriesMarkovReward; + +/** + * The goal is to calculate the pdf and cdf of time spent in a given state + * for a continuous-time Markov chain with rewards. + * + * There are several implementations that can achieve this. + * We will compare them here and output the results in a json + * file format. + * @author JT McCrone + */ +public class MarkovRewardDensityCalculator { + private double[] Q ; + private final int dim = 2; + private final double[] r = new double[]{0.0,1.0}; // rewards + + private TwoStateOccupancyMarkovReward twoStateOccupancyMarkovReward; + private TwoStateSericolaSeriesMarkovReward twoStateSericolaSeriesMarkovReward; + private SericolaSeriesMarkovReward sericolaSeriesMarkovReward; + + private MarkovRewardDensityCalculator(double rate,double bias){ + this.Q = new double[]{ + -rate * bias, rate * bias, + rate * (1.0 - bias), -rate * (1.0 - bias) + }; + + this.twoStateOccupancyMarkovReward = new TwoStateOccupancyMarkovReward(Q); + this.twoStateSericolaSeriesMarkovReward = new TwoStateSericolaSeriesMarkovReward(Q, r, dim); + this.sericolaSeriesMarkovReward = new SericolaSeriesMarkovReward(Q, r, dim); + + } + + + void run(double time){ + double s=0.0; + double[] times = new double[100]; + + double pdf; + double cdf; + double conditional; + + System.out.print("{\"data\":["); + + for(int i=0; i Date: Fri, 31 Jan 2025 11:03:40 -0800 Subject: [PATCH 02/21] updates latent model to use TwoStateSericolaMR --- .../SericolaLatentStateBranchRateModel.java | 17 +++++++++++++++-- .../LatentStateBranchRateModelParser.java | 16 ++++++++-------- 2 files changed, 23 insertions(+), 10 deletions(-) diff --git a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java index 230e7bee9b..f3ab5b9917 100644 --- a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java +++ b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java @@ -34,6 +34,7 @@ import dr.evomodel.tree.TreeParameterModel; import dr.inference.markovjumps.MarkovReward; import dr.inference.markovjumps.TwoStateOccupancyMarkovReward; +import dr.inference.markovjumps.TwoStateSericolaSeriesMarkovReward; import dr.inference.model.AbstractModelLikelihood; import dr.inference.model.Model; import dr.inference.model.Parameter; @@ -172,7 +173,8 @@ private static double[] createReward() { private MarkovReward createSeries() { // MarkovReward series = new SericolaSeriesMarkovReward(createLatentInfinitesimalMatrix(), // createReward(), 2); - MarkovReward series = new TwoStateOccupancyMarkovReward(createLatentInfinitesimalMatrix()); + // MarkovReward series = new TwoStateOccupancyMarkovReward(createLatentInfinitesimalMatrix()); + MarkovReward series = new TwoStateSericolaSeriesMarkovReward(createLatentInfinitesimalMatrix(),createReward(),2); return series; } @@ -359,6 +361,9 @@ private double calculateLogLikelihood() { double reward = branchLength * latentProportion; double density = getBranchRewardDensity(reward, branchLength); + if(latentProportion>0){ + density *=branchLength;// jacobian for sampling on 0,1 + } branchLikelihoods[node.getNumber()] = Math.log(density); } logLike += branchLikelihoods[node.getNumber()]; @@ -390,7 +395,15 @@ public double getBranchRewardDensity(double reward, double branchLength) { // Reward is [0,1], and we want to track time in latent state (= 1). // Therefore all nodes are in state 0 // double joint = series.computePdf(reward, branchLength)[state]; - double joint = series.computePdf(reward, branchLength, 0, 0); + + double joint; + + if(reward==0){ + joint = series.computeCdf(reward, branchLength, 0, 0); + }else{ + joint = series.computePdf(reward, branchLength, 0, 0); + } + double marg = series.computeConditionalProbability(branchLength, 0, 0); // TODO Overhead in creating double[] could be saved by changing signature to computePdf diff --git a/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java b/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java index 95300da423..8f57fd5476 100644 --- a/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java +++ b/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java @@ -61,14 +61,14 @@ public Object parseXMLObject(XMLObject xo) throws XMLParseException { Logger.getLogger("dr.evomodel").info("\nCreating a latent state branch rate model"); - return new LatentStateBranchRateModel(LatentStateBranchRateModel.LATENT_STATE_BRANCH_RATE_MODEL, - tree, nonLatentRateModel, - latentTransitionRateParameter, latentTransitionFrequencyParameter, /* 0/1 CTMC have two parameters */ - latentStateProportionParameter, branchCategoryProvider); -// return new SericolaLatentStateBranchRateModel(SericolaLatentStateBranchRateModel.LATENT_STATE_BRANCH_RATE_MODEL, -// tree, nonLatentRateModel, -// latentTransitionRateParameter, latentTransitionFrequencyParameter, /* 0/1 CTMC have two parameters */ -// latentStateProportionParameter, branchCategoryProvider); + // return new LatentStateBranchRateModel(LatentStateBranchRateModel.LATENT_STATE_BRANCH_RATE_MODEL, + // tree, nonLatentRateModel, + // latentTransitionRateParameter, latentTransitionFrequencyParameter, /* 0/1 CTMC have two parameters */ + // latentStateProportionParameter, branchCategoryProvider); + return new SericolaLatentStateBranchRateModel(SericolaLatentStateBranchRateModel.LATENT_STATE_BRANCH_RATE_MODEL, + tree, nonLatentRateModel, + latentTransitionRateParameter, latentTransitionFrequencyParameter, /* 0/1 CTMC have two parameters */ + latentStateProportionParameter, branchCategoryProvider); } //************************************************************************ From 1e5b3f1a857eaff5e573f701f95bf3b65e9db02e Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Fri, 31 Jan 2025 11:38:02 -0800 Subject: [PATCH 03/21] adds indicatorProductparameter to sample 0 without touching treeparametermodel --- .../app/beast/development_parsers.properties | 5 +- .../model/IndicatorProductParameter.java | 112 ++++++++++++++++++ .../IndicatorProductParameterParser.java | 51 ++++++++ 3 files changed, 167 insertions(+), 1 deletion(-) create mode 100644 src/dr/inference/model/IndicatorProductParameter.java create mode 100644 src/dr/inferencexml/model/IndicatorProductParameterParser.java diff --git a/src/dr/app/beast/development_parsers.properties b/src/dr/app/beast/development_parsers.properties index c0143ceaa9..f5bbcaa08e 100644 --- a/src/dr/app/beast/development_parsers.properties +++ b/src/dr/app/beast/development_parsers.properties @@ -1,7 +1,7 @@ # # development_parsers.properties # -# Copyright © 2002-2024 the BEAST Development Team +# Copyright � 2002-2024 the BEAST Development Team # http://beast.community/about # # This file is part of BEAST. @@ -158,6 +158,9 @@ dr.evomodelxml.coalescent.demographicmodel.PowerLawGrowthModelParser dr.evomodelxml.coalescent.demographicmodel.PeakAndDeclineModelParser dr.evomodelxml.coalescent.demographicmodel.AsymptoticGrowthModelParser +#PARAMETER +dr.inferencexml.model.IndicatorProductParameterParser + # UNIFORM INTERNAL NODE HEIGHT PRIOR dr.evomodelxml.operators.FunkyPriorMixerOperatorParser diff --git a/src/dr/inference/model/IndicatorProductParameter.java b/src/dr/inference/model/IndicatorProductParameter.java new file mode 100644 index 0000000000..e28cb0a8c4 --- /dev/null +++ b/src/dr/inference/model/IndicatorProductParameter.java @@ -0,0 +1,112 @@ +package dr.inference.model; + +import java.util.List; + +/** + * This is parameter class that is intended to allow sampling values on the range + * [0,) inclusive of 0. It represents a specific case of a produce parameter. + + */ +public class IndicatorProductParameter extends Parameter.Abstract implements VariableListener { + Parameter indicatorParameter; + Parameter continuousParameter; + + private Bounds bounds = null; + + public IndicatorProductParameter(Parameter indicatorParameter, Parameter continuousParameter) { + + this.indicatorParameter = indicatorParameter; + this.continuousParameter = continuousParameter; + + indicatorParameter.addVariableListener(this); + Parameter.CONNECTED_PARAMETER_SET.add(indicatorParameter); + + continuousParameter.addVariableListener(this); + Parameter.CONNECTED_PARAMETER_SET.add(continuousParameter); + } + public int getDimension() { + return indicatorParameter.getDimension(); + } + @Override + public boolean isImmutable() { + return true; + } + + protected void storeValues() { + indicatorParameter.storeParameterValues(); + continuousParameter.storeParameterValues(); + } + protected void restoreValues() { + indicatorParameter.restoreParameterValues(); + continuousParameter.restoreParameterValues(); + } + protected void acceptValues() { + indicatorParameter.acceptParameterValues(); + continuousParameter.acceptParameterValues(); + } + + public double getParameterValue(int dim) { + double value = indicatorParameter.getParameterValue(dim) *continuousParameter.getParameterValue(dim); + return value; + } + + public void setParameterValue(int dim, double value) { + if(value==0){ + indicatorParameter.setParameterValue(dim, 0); + }else{ + indicatorParameter.setParameterValue(dim, 1); + continuousParameter.setParameterValue(dim, value); + } + } + + public void setParameterValueQuietly(int dim, double value) { + throw new RuntimeException("Not implemented"); + } + + public void setParameterValueNotifyChangedAll(int dim, double value){ + throw new RuntimeException("Not implemented"); + } + @Override + protected void adoptValues(Parameter source) { + // TODO Auto-generated method stub + throw new UnsupportedOperationException("Unimplemented method 'adoptValues'"); + } + + public String getParameterName() { + if (getId() == null) { + StringBuilder sb = new StringBuilder("IndicatorProduct"); + + sb.append(".").append(indicatorParameter.getId()); + sb.append(".").append(continuousParameter.getId()); + + setId(sb.toString()); + } + return getId(); + } + public void addBounds(Bounds bounds) { + this.bounds = bounds; + } + + public Bounds getBounds() { + if (bounds == null) { + return indicatorParameter.getBounds(); // TODO + } else { + return bounds; + } + } + + public void addDimension(int index, double value) { + throw new RuntimeException("Not yet implemented."); + } + + public double removeDimension(int index) { + throw new RuntimeException("Not yet implemented."); + } + + public void variableChangedEvent(Variable variable, int index, ChangeType type) { + fireParameterChangedEvent(index,type); + } + + + +} diff --git a/src/dr/inferencexml/model/IndicatorProductParameterParser.java b/src/dr/inferencexml/model/IndicatorProductParameterParser.java new file mode 100644 index 0000000000..0fc14f0f5c --- /dev/null +++ b/src/dr/inferencexml/model/IndicatorProductParameterParser.java @@ -0,0 +1,51 @@ +package dr.inferencexml.model; + +import dr.inference.model.IndicatorProductParameter; +import dr.inference.model.Parameter; + +import dr.xml.*; + + +public class IndicatorProductParameterParser extends AbstractXMLObjectParser { + + public static final String INDICATOR_PRODUCT_PARAMETER = "indicatorProductParameter"; + public static final String INDICATOR = "indicator"; + public static final String CONTINUOUS_PARAMETER = "continuousParameter"; + + public Object parseXMLObject(XMLObject xo) throws XMLParseException { + + Parameter continuous = (Parameter) xo.getChild(CONTINUOUS_PARAMETER).getChild(Parameter.class); + Parameter indicator = (Parameter) xo.getChild(INDICATOR).getChild(Parameter.class); + + IndicatorProductParameter prodParam = new IndicatorProductParameter(indicator, continuous); + + return prodParam; + + + } + + public XMLSyntaxRule[] getSyntaxRules() { + return rules; + } + + private final XMLSyntaxRule[] rules = { + new ElementRule(INDICATOR, new XMLSyntaxRule[]{ + new ElementRule(Parameter.class) + }, false), + new ElementRule(CONTINUOUS_PARAMETER, new XMLSyntaxRule[]{ + new ElementRule(Parameter.class) + }, false) + }; + + public String getParserDescription() { + return "A element-wise product of an indicator and continuous paramter."; + } + + public Class getReturnType() { + return Parameter.class; + } + + public String getParserName() { + return INDICATOR_PRODUCT_PARAMETER; + } +} From f2fd62b2e6120d8dc1e5cb0a457150d1407a6d19 Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Wed, 5 Feb 2025 16:02:59 -0600 Subject: [PATCH 04/21] latent model iterates over all nodes --- .../branchratemodel/SericolaLatentStateBranchRateModel.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java index f3ab5b9917..eca3b13fb5 100644 --- a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java +++ b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java @@ -350,7 +350,7 @@ private double calculateLogLikelihood() { double logLike = 0.0; - for (int i = 0; i < tree.getInternalNodeCount(); ++i) { + for (int i = 0; i < tree.getNodeCount(); ++i) { NodeRef node = tree.getNode(i); if (node != tree.getRoot()) { if (updateNeededForNode(tree, node)) { From 22da8c0b024d74dbb46d6993b49f3c60ebfe4141 Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Thu, 6 Feb 2025 15:45:58 -0600 Subject: [PATCH 05/21] adds option to disallow latency on root children --- .../SericolaLatentStateBranchRateModel.java | 18 +++++++++++++++++- .../LatentStateBranchRateModelParser.java | 5 +++-- 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java index eca3b13fb5..3596388fb9 100644 --- a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java +++ b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java @@ -91,6 +91,7 @@ public class SericolaLatentStateBranchRateModel extends AbstractModelLikelihood private boolean[] updateCategory; private boolean[] storedUpdateCategory; + private boolean excludeRoot; public SericolaLatentStateBranchRateModel(String name, TreeModel treeModel, @@ -98,7 +99,8 @@ public SericolaLatentStateBranchRateModel(String name, Parameter latentTransitionRateParameter, Parameter latentTransitionFrequencyParameter, Parameter latentStateProportionParameter, - CountableBranchCategoryProvider branchCategoryProvider) { + CountableBranchCategoryProvider branchCategoryProvider, + boolean excludeRoot) { super(name); this.tree = treeModel; @@ -142,8 +144,18 @@ public SericolaLatentStateBranchRateModel(String name, setUpdateAllBranches(); } + this.excludeRoot=excludeRoot; } + public SericolaLatentStateBranchRateModel(String name, + TreeModel treeModel, + BranchRateModel nonLatentRateModel, + Parameter latentTransitionRateParameter, + Parameter latentTransitionFrequencyParameter, + Parameter latentStateProportionParameter, + CountableBranchCategoryProvider branchCategoryProvider){ + this(name, treeModel, nonLatentRateModel, latentTransitionRateParameter, latentTransitionFrequencyParameter, latentStateProportionParameter, branchCategoryProvider, false); + } public SericolaLatentStateBranchRateModel(Parameter rate, Parameter prop) { super(LATENT_STATE_BRANCH_RATE_MODEL); tree = null; @@ -359,6 +371,10 @@ private double calculateLogLikelihood() { assert(latentProportion < 1.0); + if(excludeRoot && tree.isRoot(tree.getParent(node)) && latentProportion>0){ + return Double.NEGATIVE_INFINITY; + } + double reward = branchLength * latentProportion; double density = getBranchRewardDensity(reward, branchLength); if(latentProportion>0){ diff --git a/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java b/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java index 8f57fd5476..87a927eb1d 100644 --- a/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java +++ b/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java @@ -40,6 +40,7 @@ public class LatentStateBranchRateModelParser extends AbstractXMLObjectParser { public static final String LATENT_TRANSITION_RATE = "latentTransitionRate"; public static final String LATENT_TRANSITION_FREQUENCY = "latentTransitionFrequency"; public static final String LATENT_STATE_PROPORTIONS = "latentStateProportions"; + public static final String EXCLUDE_ROOT = "excludeRoot"; public String getParserName() { @@ -58,7 +59,7 @@ public Object parseXMLObject(XMLObject xo) throws XMLParseException { if (xo.hasChildNamed(LATENT_STATE_PROPORTIONS)) { latentStateProportionParameter = (Parameter) xo.getElementFirstChild(LATENT_STATE_PROPORTIONS); } - + boolean excludeRoot = xo.getAttribute(EXCLUDE_ROOT, false); Logger.getLogger("dr.evomodel").info("\nCreating a latent state branch rate model"); // return new LatentStateBranchRateModel(LatentStateBranchRateModel.LATENT_STATE_BRANCH_RATE_MODEL, @@ -68,7 +69,7 @@ public Object parseXMLObject(XMLObject xo) throws XMLParseException { return new SericolaLatentStateBranchRateModel(SericolaLatentStateBranchRateModel.LATENT_STATE_BRANCH_RATE_MODEL, tree, nonLatentRateModel, latentTransitionRateParameter, latentTransitionFrequencyParameter, /* 0/1 CTMC have two parameters */ - latentStateProportionParameter, branchCategoryProvider); + latentStateProportionParameter, branchCategoryProvider,excludeRoot); } //************************************************************************ From 5ebe5fc414c256a0af589fe5aa6540e9eaa3b367 Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Fri, 7 Mar 2025 20:40:44 -0600 Subject: [PATCH 06/21] fully implements option to disallow latency on root children --- .project | 11 +++++++++++ packaging_tools/windows/launch4j/.project | 11 +++++++++++ packaging_tools/windows/launch4j/maven/.project | 11 +++++++++++ .../SericolaLatentStateBranchRateModel.java | 11 +++++++---- 4 files changed, 40 insertions(+), 4 deletions(-) diff --git a/.project b/.project index 9662e139ba..3adcbc53ff 100644 --- a/.project +++ b/.project @@ -14,4 +14,15 @@ org.eclipse.jdt.core.javanature + + + 1738341918935 + + 30 + + org.eclipse.core.resources.regexFilterMatcher + node_modules|\.git|__CREATED_BY_JAVA_LANGUAGE_SERVER__ + + + diff --git a/packaging_tools/windows/launch4j/.project b/packaging_tools/windows/launch4j/.project index 0607fa0ec4..5261e55db2 100644 --- a/packaging_tools/windows/launch4j/.project +++ b/packaging_tools/windows/launch4j/.project @@ -14,4 +14,15 @@ org.eclipse.jdt.core.javanature + + + 1738341918954 + + 30 + + org.eclipse.core.resources.regexFilterMatcher + node_modules|\.git|__CREATED_BY_JAVA_LANGUAGE_SERVER__ + + + diff --git a/packaging_tools/windows/launch4j/maven/.project b/packaging_tools/windows/launch4j/maven/.project index 08f86897f0..4aaabb4986 100644 --- a/packaging_tools/windows/launch4j/maven/.project +++ b/packaging_tools/windows/launch4j/maven/.project @@ -20,4 +20,15 @@ org.eclipse.m2e.core.maven2Nature org.eclipse.jdt.core.javanature + + + 1738341918944 + + 30 + + org.eclipse.core.resources.regexFilterMatcher + node_modules|\.git|__CREATED_BY_JAVA_LANGUAGE_SERVER__ + + + diff --git a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java index 3596388fb9..277c9ad2bd 100644 --- a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java +++ b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java @@ -200,6 +200,10 @@ public double getBranchRate(Tree tree, NodeRef node) { } public double getLatentProportion(Tree tree, NodeRef node) { + //if excluderoot is true it's children don't have latency. + if(excludeRoot&&tree.isRoot(tree.getParent(node))){ + return 0; + } if (latentStateProportions != null) { return latentStateProportions.getNodeValue(tree, node); @@ -364,6 +368,9 @@ private double calculateLogLikelihood() { for (int i = 0; i < tree.getNodeCount(); ++i) { NodeRef node = tree.getNode(i); + if(excludeRoot && tree.isRoot(tree.getParent(node))){ // don't include these nodes if excludeRoot is true + continue; + } if (node != tree.getRoot()) { if (updateNeededForNode(tree, node)) { double branchLength = tree.getBranchLength(node); @@ -371,10 +378,6 @@ private double calculateLogLikelihood() { assert(latentProportion < 1.0); - if(excludeRoot && tree.isRoot(tree.getParent(node)) && latentProportion>0){ - return Double.NEGATIVE_INFINITY; - } - double reward = branchLength * latentProportion; double density = getBranchRewardDensity(reward, branchLength); if(latentProportion>0){ From 54447bbe174c0ac301ea15f60990db1b00b79b05 Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Tue, 11 Mar 2025 17:58:58 -0700 Subject: [PATCH 07/21] Revert "fully implements option to disallow latency on root children" This reverts commit 5ebe5fc414c256a0af589fe5aa6540e9eaa3b367. --- .project | 11 ----------- packaging_tools/windows/launch4j/.project | 11 ----------- packaging_tools/windows/launch4j/maven/.project | 11 ----------- .../SericolaLatentStateBranchRateModel.java | 11 ++++------- 4 files changed, 4 insertions(+), 40 deletions(-) diff --git a/.project b/.project index 3adcbc53ff..9662e139ba 100644 --- a/.project +++ b/.project @@ -14,15 +14,4 @@ org.eclipse.jdt.core.javanature - - - 1738341918935 - - 30 - - org.eclipse.core.resources.regexFilterMatcher - node_modules|\.git|__CREATED_BY_JAVA_LANGUAGE_SERVER__ - - - diff --git a/packaging_tools/windows/launch4j/.project b/packaging_tools/windows/launch4j/.project index 5261e55db2..0607fa0ec4 100644 --- a/packaging_tools/windows/launch4j/.project +++ b/packaging_tools/windows/launch4j/.project @@ -14,15 +14,4 @@ org.eclipse.jdt.core.javanature - - - 1738341918954 - - 30 - - org.eclipse.core.resources.regexFilterMatcher - node_modules|\.git|__CREATED_BY_JAVA_LANGUAGE_SERVER__ - - - diff --git a/packaging_tools/windows/launch4j/maven/.project b/packaging_tools/windows/launch4j/maven/.project index 4aaabb4986..08f86897f0 100644 --- a/packaging_tools/windows/launch4j/maven/.project +++ b/packaging_tools/windows/launch4j/maven/.project @@ -20,15 +20,4 @@ org.eclipse.m2e.core.maven2Nature org.eclipse.jdt.core.javanature - - - 1738341918944 - - 30 - - org.eclipse.core.resources.regexFilterMatcher - node_modules|\.git|__CREATED_BY_JAVA_LANGUAGE_SERVER__ - - - diff --git a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java index 277c9ad2bd..3596388fb9 100644 --- a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java +++ b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java @@ -200,10 +200,6 @@ public double getBranchRate(Tree tree, NodeRef node) { } public double getLatentProportion(Tree tree, NodeRef node) { - //if excluderoot is true it's children don't have latency. - if(excludeRoot&&tree.isRoot(tree.getParent(node))){ - return 0; - } if (latentStateProportions != null) { return latentStateProportions.getNodeValue(tree, node); @@ -368,9 +364,6 @@ private double calculateLogLikelihood() { for (int i = 0; i < tree.getNodeCount(); ++i) { NodeRef node = tree.getNode(i); - if(excludeRoot && tree.isRoot(tree.getParent(node))){ // don't include these nodes if excludeRoot is true - continue; - } if (node != tree.getRoot()) { if (updateNeededForNode(tree, node)) { double branchLength = tree.getBranchLength(node); @@ -378,6 +371,10 @@ private double calculateLogLikelihood() { assert(latentProportion < 1.0); + if(excludeRoot && tree.isRoot(tree.getParent(node)) && latentProportion>0){ + return Double.NEGATIVE_INFINITY; + } + double reward = branchLength * latentProportion; double density = getBranchRewardDensity(reward, branchLength); if(latentProportion>0){ From af3cf91bf04a33e7cdf27a3fe486b326fefde631 Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Tue, 11 Mar 2025 18:37:19 -0700 Subject: [PATCH 08/21] adds epsilon to xml --- .../SericolaLatentStateBranchRateModel.java | 11 ++++------- .../LatentStateBranchRateModelParser.java | 6 +++--- 2 files changed, 7 insertions(+), 10 deletions(-) diff --git a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java index 3596388fb9..910d24146b 100644 --- a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java +++ b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java @@ -82,6 +82,7 @@ public class SericolaLatentStateBranchRateModel extends AbstractModelLikelihood private boolean storedLikelihoodKnown; private double logLikelihood; private double storedLogLikelihood; + private double epsilon; private double[] branchLikelihoods; private double[] storedbranchLikelihoods; @@ -91,7 +92,6 @@ public class SericolaLatentStateBranchRateModel extends AbstractModelLikelihood private boolean[] updateCategory; private boolean[] storedUpdateCategory; - private boolean excludeRoot; public SericolaLatentStateBranchRateModel(String name, TreeModel treeModel, @@ -100,7 +100,7 @@ public SericolaLatentStateBranchRateModel(String name, Parameter latentTransitionFrequencyParameter, Parameter latentStateProportionParameter, CountableBranchCategoryProvider branchCategoryProvider, - boolean excludeRoot) { + double epsilon) { super(name); this.tree = treeModel; @@ -144,7 +144,7 @@ public SericolaLatentStateBranchRateModel(String name, setUpdateAllBranches(); } - this.excludeRoot=excludeRoot; + this.epsilon=epsilon; } public SericolaLatentStateBranchRateModel(String name, @@ -154,7 +154,7 @@ public SericolaLatentStateBranchRateModel(String name, Parameter latentTransitionFrequencyParameter, Parameter latentStateProportionParameter, CountableBranchCategoryProvider branchCategoryProvider){ - this(name, treeModel, nonLatentRateModel, latentTransitionRateParameter, latentTransitionFrequencyParameter, latentStateProportionParameter, branchCategoryProvider, false); + this(name, treeModel, nonLatentRateModel, latentTransitionRateParameter, latentTransitionFrequencyParameter, latentStateProportionParameter, branchCategoryProvider, 1E-10); } public SericolaLatentStateBranchRateModel(Parameter rate, Parameter prop) { super(LATENT_STATE_BRANCH_RATE_MODEL); @@ -371,9 +371,6 @@ private double calculateLogLikelihood() { assert(latentProportion < 1.0); - if(excludeRoot && tree.isRoot(tree.getParent(node)) && latentProportion>0){ - return Double.NEGATIVE_INFINITY; - } double reward = branchLength * latentProportion; double density = getBranchRewardDensity(reward, branchLength); diff --git a/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java b/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java index 87a927eb1d..932945ddad 100644 --- a/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java +++ b/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java @@ -40,7 +40,7 @@ public class LatentStateBranchRateModelParser extends AbstractXMLObjectParser { public static final String LATENT_TRANSITION_RATE = "latentTransitionRate"; public static final String LATENT_TRANSITION_FREQUENCY = "latentTransitionFrequency"; public static final String LATENT_STATE_PROPORTIONS = "latentStateProportions"; - public static final String EXCLUDE_ROOT = "excludeRoot"; + public static final String EPSILON = "epsilon"; public String getParserName() { @@ -59,7 +59,7 @@ public Object parseXMLObject(XMLObject xo) throws XMLParseException { if (xo.hasChildNamed(LATENT_STATE_PROPORTIONS)) { latentStateProportionParameter = (Parameter) xo.getElementFirstChild(LATENT_STATE_PROPORTIONS); } - boolean excludeRoot = xo.getAttribute(EXCLUDE_ROOT, false); + double epsilon = xo.getAttribute(EPSILON, 1E-10); Logger.getLogger("dr.evomodel").info("\nCreating a latent state branch rate model"); // return new LatentStateBranchRateModel(LatentStateBranchRateModel.LATENT_STATE_BRANCH_RATE_MODEL, @@ -69,7 +69,7 @@ public Object parseXMLObject(XMLObject xo) throws XMLParseException { return new SericolaLatentStateBranchRateModel(SericolaLatentStateBranchRateModel.LATENT_STATE_BRANCH_RATE_MODEL, tree, nonLatentRateModel, latentTransitionRateParameter, latentTransitionFrequencyParameter, /* 0/1 CTMC have two parameters */ - latentStateProportionParameter, branchCategoryProvider,excludeRoot); + latentStateProportionParameter, branchCategoryProvider,epsilon); } //************************************************************************ From 5adaee307ce0ca1883a9066397f27d7b4c032b85 Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Wed, 17 Sep 2025 09:05:46 -0500 Subject: [PATCH 09/21] allow for excluding root children from latency - test --- .../SericolaLatentStateBranchRateModel.java | 18 +++++++++++++++--- .../LatentStateBranchRateModelParser.java | 4 +++- 2 files changed, 18 insertions(+), 4 deletions(-) diff --git a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java index 910d24146b..ecd1aec502 100644 --- a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java +++ b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java @@ -80,6 +80,7 @@ public class SericolaLatentStateBranchRateModel extends AbstractModelLikelihood private MarkovReward storedSeries; private boolean likelihoodKnown = false; private boolean storedLikelihoodKnown; + private boolean excludeRoot; private double logLikelihood; private double storedLogLikelihood; private double epsilon; @@ -100,11 +101,14 @@ public SericolaLatentStateBranchRateModel(String name, Parameter latentTransitionFrequencyParameter, Parameter latentStateProportionParameter, CountableBranchCategoryProvider branchCategoryProvider, - double epsilon) { + double epsilon, + boolean excludeRoot) { super(name); this.tree = treeModel; addModel(tree); + + this.excludeRoot = excludeRoot; this.nonLatentRateModel = nonLatentRateModel; addModel(nonLatentRateModel); @@ -154,7 +158,7 @@ public SericolaLatentStateBranchRateModel(String name, Parameter latentTransitionFrequencyParameter, Parameter latentStateProportionParameter, CountableBranchCategoryProvider branchCategoryProvider){ - this(name, treeModel, nonLatentRateModel, latentTransitionRateParameter, latentTransitionFrequencyParameter, latentStateProportionParameter, branchCategoryProvider, 1E-10); + this(name, treeModel, nonLatentRateModel, latentTransitionRateParameter, latentTransitionFrequencyParameter, latentStateProportionParameter, branchCategoryProvider, 1E-10,false); } public SericolaLatentStateBranchRateModel(Parameter rate, Parameter prop) { super(LATENT_STATE_BRANCH_RATE_MODEL); @@ -366,6 +370,13 @@ private double calculateLogLikelihood() { NodeRef node = tree.getNode(i); if (node != tree.getRoot()) { if (updateNeededForNode(tree, node)) { + if(excludeRoot && tree.getParent(node)== tree.getRoot()){ + if(getLatentProportion(tree,node)>0){ + // branchLikelihoods[node.getNumber()]=Double.NEGATIVE_INFINITY; // reject we don't allow latency on these branches + return Double.NEGATIVE_INFINITY; + } + branchLikelihoods[node.getNumber()]=0.0; // if it is zero we don't get credit for it since it couldn't be otherwise + }else{ double branchLength = tree.getBranchLength(node); double latentProportion = getLatentProportion(tree, node); @@ -379,11 +390,12 @@ private double calculateLogLikelihood() { } branchLikelihoods[node.getNumber()] = Math.log(density); } + } logLike += branchLikelihoods[node.getNumber()]; // TODO More importantly, MH proposals on [0,1] may be missing a Jacobian for which we should adjust. // TODO This is easy to test and we should do it when sampling appears to work. - } } + } clearUpdateAllBranches(); clearAllCategories(); diff --git a/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java b/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java index 932945ddad..8a216e72d7 100644 --- a/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java +++ b/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java @@ -41,6 +41,7 @@ public class LatentStateBranchRateModelParser extends AbstractXMLObjectParser { public static final String LATENT_TRANSITION_FREQUENCY = "latentTransitionFrequency"; public static final String LATENT_STATE_PROPORTIONS = "latentStateProportions"; public static final String EPSILON = "epsilon"; + public static final String EXCLUDE_ROOT = "excludeRoot"; public String getParserName() { @@ -52,6 +53,7 @@ public Object parseXMLObject(XMLObject xo) throws XMLParseException { TreeModel tree = (TreeModel) xo.getChild(TreeModel.class); Parameter latentTransitionRateParameter = (Parameter) xo.getElementFirstChild(LATENT_TRANSITION_RATE); Parameter latentTransitionFrequencyParameter = (Parameter) xo.getElementFirstChild(LATENT_TRANSITION_FREQUENCY); + boolean excludeRoot = xo.getAttribute("excludeRoot", false); CountableBranchCategoryProvider branchCategoryProvider = (CountableBranchCategoryProvider)xo.getChild(CountableBranchCategoryProvider.class); @@ -69,7 +71,7 @@ public Object parseXMLObject(XMLObject xo) throws XMLParseException { return new SericolaLatentStateBranchRateModel(SericolaLatentStateBranchRateModel.LATENT_STATE_BRANCH_RATE_MODEL, tree, nonLatentRateModel, latentTransitionRateParameter, latentTransitionFrequencyParameter, /* 0/1 CTMC have two parameters */ - latentStateProportionParameter, branchCategoryProvider,epsilon); + latentStateProportionParameter, branchCategoryProvider,epsilon,excludeRoot); } //************************************************************************ From 93ab2744c54b99556a8cf28815f7a0c9c254f086 Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Wed, 17 Sep 2025 09:37:00 -0500 Subject: [PATCH 10/21] pass epsilon from xml to model --- .../branchratemodel/SericolaLatentStateBranchRateModel.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java index ecd1aec502..47cfed56c1 100644 --- a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java +++ b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java @@ -190,7 +190,7 @@ private MarkovReward createSeries() { // MarkovReward series = new SericolaSeriesMarkovReward(createLatentInfinitesimalMatrix(), // createReward(), 2); // MarkovReward series = new TwoStateOccupancyMarkovReward(createLatentInfinitesimalMatrix()); - MarkovReward series = new TwoStateSericolaSeriesMarkovReward(createLatentInfinitesimalMatrix(),createReward(),2); + MarkovReward series = new TwoStateSericolaSeriesMarkovReward(createLatentInfinitesimalMatrix(),createReward(),2,epsilon); return series; } From f69571821028f93fa3f6525d0e6aec76fcb46cf7 Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Thu, 18 Sep 2025 20:07:25 -0500 Subject: [PATCH 11/21] speed up cdf for when there is no latency --- .../SericolaLatentStateBranchRateModel.java | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java index 47cfed56c1..94bea79f8a 100644 --- a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java +++ b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java @@ -420,12 +420,24 @@ public double getBranchRewardDensity(double reward, double branchLength) { // Reward is [0,1], and we want to track time in latent state (= 1). // Therefore all nodes are in state 0 // double joint = series.computePdf(reward, branchLength)[state]; + double rate = latentTransitionRateParameter.getParameterValue(0); + double prop = latentTransitionFrequencyParameter.getParameterValue(0); + double noJump = Math.exp(-1 * rate * prop * branchLength); // no jump double joint; - if(reward==0){ - joint = series.computeCdf(reward, branchLength, 0, 0); - }else{ + if (reward == 0) { + // joint = series.computeCdf(reward, branchLength, 0, 0); + // no jump. + joint = noJump; + // series.computeCdf(reward, branchLength, 0, 0)- is no jumps using the series + // machinery + assert Math.abs(series.computeCdf(reward, branchLength, 0, 0) - joint) < epsilon + : "0 jumps not calculated correctly: " + + Math.abs(series.computeCdf(reward, branchLength, 0, 0) - joint) + "epsilon: " + epsilon + + " rate: " + rate + " prop: " + prop + " bl: " + branchLength; + // System.out.println(joint + " = " + noJump); + } else { joint = series.computePdf(reward, branchLength, 0, 0); } From 6c428ddb0e186f623f55fd971f2860d021947513 Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Thu, 18 Sep 2025 20:57:02 -0500 Subject: [PATCH 12/21] caching for some speed ups --- .../TwoStateSericolaSeriesMarkovReward.java | 28 +++++++++++++++++-- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/src/dr/inference/markovjumps/TwoStateSericolaSeriesMarkovReward.java b/src/dr/inference/markovjumps/TwoStateSericolaSeriesMarkovReward.java index b8f8a192e9..7daffa3ee4 100644 --- a/src/dr/inference/markovjumps/TwoStateSericolaSeriesMarkovReward.java +++ b/src/dr/inference/markovjumps/TwoStateSericolaSeriesMarkovReward.java @@ -262,7 +262,7 @@ private double accumulatePdf(double x, int n, double time, int uv) { double w = 0.0; final double premult = Math.exp( - -lambda * time + n * (Math.log(lambda) + Math.log(time)) - GammaFunction.lnGamma(n + 1.0) + -lambda * time + n * (Math.log(lambda) + Math.log(time)) - getLnGamma(n + 1.0) ); // TODO Make factorial/choose static look-up tables @@ -278,7 +278,7 @@ private double accumulatePdf(double x, int n, double time, int uv) { // double[] inc = new double[dim2]; // W^{\epsilon}(x(i),t,n) double inc = 0.0; for (int k = 0; k <= n; k++) { - final double binomialCoef = Binomial.choose(n, k) * Math.pow(xh, k) * Math.pow(1.0 - xh, n - k); + final double binomialCoef = choose(n, k) * Math.pow(xh, k) * Math.pow(1.0 - xh, n - k); // for (int uv = 0; uv < dim2; ++uv) { inc += binomialCoef * (C(h, n + 1, k + 1)[uv] - C(h, n + 1, k)[uv]); // } @@ -471,7 +471,7 @@ private int determineNumberOfSteps(double time, double lambda) { while (Math.abs(sum2 - tolerance2) > epsilon && sum2 < 1.0) { // while (sum2 < tolerance2) { i++; - double logDensity = -lambda * time + i * (Math.log(lambda) + Math.log(time)) - GammaFunction.lnGamma(i + 1); + double logDensity = -lambda * time + i * (Math.log(lambda) + Math.log(time)) -getLnGamma(i + 1); sum2 += Math.exp(logDensity); // sum2 = LogTricks.logSum(sum2, logDensity); if (DEBUG2) { @@ -487,6 +487,25 @@ private int determineNumberOfSteps(double time, double lambda) { return i; } + private double choose(int n, int k){ + if (k < 0 || k > n) return 0; + double lchoose = getLnGamma(n + 1.0) - + getLnGamma(k + 1.0) - getLnGamma(n - k + 1.0); + + return Math.floor(Math.exp(lchoose) + 0.5); + } + // should be the same as calling GammaFunction.lnGamma + private double getLnGamma(double n){ + n = Math.floor(n + 0.5); + if (lnGamma == null || lnGamma.length <= n ) { // fill it up! + lnGamma = new double[(int) (n+50)]; + for (int i = 0; i < (n+50); ++i) { + lnGamma[i] = GammaFunction.lnGamma(i); + } + } + return lnGamma[(int) n]; + } + public String toString() { StringBuilder sb = new StringBuilder(); sb.append("Q: " + new Vector(Q) + "\n"); @@ -528,6 +547,9 @@ public double computeConditionalProbability(double distance, int i, int j) { private final int dim; private final double epsilon; + private double[][] binomialCoefficients; + private double[] lnGamma; + private final EigenSystem eigenSystem; private double maxTime; From 51ed95e8a6e646d1dc35f2931f1d5143e8effc1c Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Thu, 18 Sep 2025 21:33:40 -0500 Subject: [PATCH 13/21] cache powers and loop calcs for more speed --- .../TwoStateSericolaSeriesMarkovReward.java | 45 ++++++++++++++----- 1 file changed, 34 insertions(+), 11 deletions(-) diff --git a/src/dr/inference/markovjumps/TwoStateSericolaSeriesMarkovReward.java b/src/dr/inference/markovjumps/TwoStateSericolaSeriesMarkovReward.java index 7daffa3ee4..31830fc93e 100644 --- a/src/dr/inference/markovjumps/TwoStateSericolaSeriesMarkovReward.java +++ b/src/dr/inference/markovjumps/TwoStateSericolaSeriesMarkovReward.java @@ -145,8 +145,18 @@ public double computePdf(double x, double time, int i, int j) { double w = 0.0; final int N = getNfromC() - 1; + + // moving these caluculations out of the loop + final double premultFactor = (Math.log(lambda) + Math.log(time)); + final int h = 1; + final double factor = lambda / (r[h] - r[h - 1]); + double xh = (x - r[h - 1] * time) / ((r[h] - r[h - 1]) * time); + for (int n = 0; n <= N; ++n) { - w += accumulatePdf(x, n, time, uv); + final double premult = Math.exp( + -lambda * time + n * premultFactor - getLnGamma(n + 1.0) + ); + w += accumulatePdf(x, n, time, uv,premult,xh,factor); } // if (DEBUG2) { @@ -257,33 +267,46 @@ private void accumulateCdf(double[][] W, double[] X, int[] H, int n, double time } } - private double accumulatePdf(double x, int n, double time, int uv) { + + private double accumulatePdf(double x, int n, double time, int uv, double premult, double xh,double factor) { double w = 0.0; - final double premult = Math.exp( - -lambda * time + n * (Math.log(lambda) + Math.log(time)) - getLnGamma(n + 1.0) - ); + // final double premult = Math.exp( + // // -lambda * time + n * (Math.log(lambda) + Math.log(time)) - getLnGamma(n + 1.0) + // -lambda * time + n * (Math.log(lambda) + Math.log(time)) - getLnGamma(n + 1.0) + // ); // TODO Make factorial/choose static look-up tables // for (int t = 0; t < X.length; ++t) { // For each time point int h = 1; - final double factor = lambda / (r[h] - r[h - 1]); + // final double factor = lambda / (r[h] - r[h - 1]); - double xh = (x - r[h - 1] * time) / ((r[h] - r[h - 1]) * time); + // double xh = (x - r[h - 1] * time) / ((r[h] - r[h - 1]) * time); // final int dim2 = dim * dim; // double[] inc = new double[dim2]; // W^{\epsilon}(x(i),t,n) double inc = 0.0; + double xhPow = 1.0; + double oneMinusXhPow = Math.pow(1.0 - xh, n); for (int k = 0; k <= n; k++) { - final double binomialCoef = choose(n, k) * Math.pow(xh, k) * Math.pow(1.0 - xh, n - k); -// for (int uv = 0; uv < dim2; ++uv) { - inc += binomialCoef * (C(h, n + 1, k + 1)[uv] - C(h, n + 1, k)[uv]); -// } + double binomialCoef = choose(n, k) * xhPow * oneMinusXhPow; + inc += binomialCoef * (C(h, n + 1, k + 1)[uv] - C(h, n + 1, k)[uv]); + // update powers for next iteration + xhPow *= xh; + oneMinusXhPow /= (1.0 - xh); } + +// for (int k = 0; k <= n; k++) { +// final double binomialCoef = choose(n, k) * Math.pow(xh, k) * Math.pow(1.0 - xh, n - k); +// // for (int uv = 0; uv < dim2; ++uv) { +// inc += binomialCoef * (C(h, n + 1, k + 1)[uv] - C(h, n + 1, k)[uv]); +// // } +// } + // for (int uv = 0; uv < dim2; ++uv) { w += factor * premult * inc; // } From edf73f9ddb018290a6cf953399f8e041104ce19e Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Mon, 13 Oct 2025 20:44:44 -0500 Subject: [PATCH 14/21] move latent work to package adds latent statistic --- src/dr/app/beast/release_parsers.properties | 8 +- .../LatentStateBranchRateLengthStatistic.java | 127 ++++++++++++++++++ .../LatentStateBranchRateModel.java | 4 +- ...tStateBranchRateLengthStatisticParser.java | 44 ++++++ .../LatentStateBranchRateModelParser.java | 2 +- .../TwoStateOccupancyMarkovRewardsTest.java | 2 +- .../dr/math/GeneralizedIntegerGammaTest.java | 2 +- 7 files changed, 183 insertions(+), 6 deletions(-) create mode 100644 src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatistic.java rename src/dr/evomodel/branchratemodel/{ => latentStateBranchRate}/LatentStateBranchRateModel.java (98%) create mode 100644 src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatisticParser.java rename src/dr/evomodelxml/branchratemodel/{ => latentStateBranchRate}/LatentStateBranchRateModelParser.java (96%) diff --git a/src/dr/app/beast/release_parsers.properties b/src/dr/app/beast/release_parsers.properties index e68fdf2578..7fb3eef235 100644 --- a/src/dr/app/beast/release_parsers.properties +++ b/src/dr/app/beast/release_parsers.properties @@ -1,7 +1,7 @@ # # release_parsers.properties # -# Copyright © 2002-2024 the BEAST Development Team +# Copyright � 2002-2024 the BEAST Development Team # http://beast.community/about # # This file is part of BEAST. @@ -194,7 +194,6 @@ dr.evomodelxml.branchratemodel.RandomLocalClockModelParser dr.evomodelxml.tree.RLTVLoggerOnTreeParser dr.evomodelxml.branchratemodel.BranchCategoriesParser dr.evomodelxml.branchratemodel.CountableMixtureBranchRatesParser -dr.evomodelxml.branchratemodel.LatentStateBranchRateModelParser dr.evomodelxml.branchratemodel.RelaxedDriftModelParser dr.evomodelxml.branchratemodel.FixedDriftModelParser dr.evomodelxml.branchratemodel.BranchSpecificFixedEffectsParser @@ -204,6 +203,11 @@ dr.evomodelxml.branchratemodel.BranchRateGradientWrtIncrementsParser dr.evomodelxml.branchratemodel.AncestralTraitBranchRatesParser dr.evomodelxml.branchratemodel.RandomEffectsTreeTraitProviderParser +## Latent model +dr.evomodelxml.branchratemodel.LatentStateBranchRateModelParser +dr.evomodelxml.branchratemodel.latentStateBranchRate.LatentStateBranchRateLengthStatisticParser + + #COALESCENT dr.evomodelxml.coalescent.demographicmodel.CataclysmicDemographicModelParser dr.evomodelxml.coalescent.demographicmodel.ExpConstExpDemographicModelParser diff --git a/src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatistic.java b/src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatistic.java new file mode 100644 index 0000000000..bc209ec591 --- /dev/null +++ b/src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatistic.java @@ -0,0 +1,127 @@ +package dr.evomodel.branchratemodel.latentStateBranchRate; + +import dr.evolution.tree.Tree; +import dr.evomodel.branchratemodel.SericolaLatentStateBranchRateModel; +import dr.evomodel.tree.TreeModel; +import dr.evomodel.tree.TreeStatistic; +import dr.inference.model.Bounds; +import dr.inference.model.Model; +import dr.inference.model.ModelListener; +import dr.inference.model.Parameter; +/** + * This statistic reports the length a tree spent either replicating or in the latent state. + * It is a parameter so it can be used in distributions like the gamma distribution. + * @author JT McCrone + * + */ + + +public class LatentStateBranchRateLengthStatistic extends Parameter.Abstract implements ModelListener { + public LatentStateBranchRateLengthStatistic( SericolaLatentStateBranchRateModel model, Tree tree, STATE state) { + + this.model = model; + + this.state = state; + this.tree= tree; + model.addModelListener(this); + if(tree instanceof TreeModel){ + ((TreeModel) tree).addModelListener(this); + } + } + + + public int getDimension() { + return 1; + } + + /** + * @return the total length of all the branches in the tree + */ + public double getParameterValue(int dim) { + double length =0.0; + for(int i=0; i< tree.getNodeCount(); i++) { + if(!tree.isRoot(tree.getNode(i))){ + double branchLength =tree.getBranchLength(tree.getNode(i)); + double latentProportion = model.getLatentProportion(tree, tree.getNode(i)); + if(state == STATE.REPLICATING) { + length += branchLength * (1-latentProportion); + } else if(state == STATE.LATENT) { + length += branchLength * latentProportion; + } + } + } + return length; + } + + @Override + public boolean isImmutable() { + return true; + } + + protected void storeValues() { + //do nothing + } + + protected void restoreValues() { + //do nothing + } + + protected void acceptValues() { + //do nothing + } + + protected void adoptValues(Parameter source) { + throw new RuntimeException("Not implemented"); + } + + + public void setParameterValue(int dim, double value) { + throw new RuntimeException("Not implemented"); + } + + public void setParameterValueQuietly(int dim, double value) { + throw new RuntimeException("Not implemented"); + } + + public void setParameterValueNotifyChangedAll(int dim, double value){ + throw new RuntimeException("Not implemented"); + } + + public String getParameterName() { + return getId(); + } + + public void addBounds(Bounds bounds) { + this.bounds = bounds; + } + + public Bounds getBounds() { + return bounds; + } + + public void addDimension(int index, double value) { + throw new RuntimeException("Not yet implemented."); + } + + public double removeDimension(int index) { + throw new RuntimeException("Not yet implemented."); + } + @Override + public void modelChangedEvent(Model model, Object object, int index) { + fireParameterChangedEvent(); + } + + + @Override + public void modelRestored(Model model) { + // do nothing + } + + private SericolaLatentStateBranchRateModel model = null; + private STATE state; + private Tree tree; + public static String LATENT_STATE_BRANCH_RATE_LENGTH_STATISTIC = "latentStateBranchRateLengthStatistic"; + public enum STATE {REPLICATING, LATENT} + private Bounds bounds = null; + +} diff --git a/src/dr/evomodel/branchratemodel/LatentStateBranchRateModel.java b/src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateModel.java similarity index 98% rename from src/dr/evomodel/branchratemodel/LatentStateBranchRateModel.java rename to src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateModel.java index 21a6616baa..df56d87a4c 100644 --- a/src/dr/evomodel/branchratemodel/LatentStateBranchRateModel.java +++ b/src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateModel.java @@ -25,11 +25,13 @@ * */ -package dr.evomodel.branchratemodel; +package dr.evomodel.branchratemodel.latentStateBranchRate; import dr.evolution.tree.NodeRef; import dr.evolution.tree.Tree; import dr.evolution.tree.TreeTrait; +import dr.evomodel.branchratemodel.BranchRateModel; +import dr.evomodel.branchratemodel.CountableBranchCategoryProvider; import dr.evomodel.tree.TreeModel; import dr.evomodel.tree.TreeParameterModel; import dr.inference.markovjumps.TwoStateOccupancyMarkovReward; diff --git a/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatisticParser.java b/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatisticParser.java new file mode 100644 index 0000000000..3a09cf9d3d --- /dev/null +++ b/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatisticParser.java @@ -0,0 +1,44 @@ +package dr.evomodelxml.branchratemodel.latentStateBranchRate; + + +import dr.evomodel.branchratemodel.SericolaLatentStateBranchRateModel; +import dr.evomodel.branchratemodel.latentStateBranchRate.LatentStateBranchRateLengthStatistic; +import dr.evomodel.tree.TreeModel; +import dr.xml.AbstractXMLObjectParser; +import dr.xml.ElementRule; +import dr.xml.XMLObject; +import dr.xml.XMLParseException; +import dr.xml.XMLSyntaxRule; + +public class LatentStateBranchRateLengthStatisticParser extends AbstractXMLObjectParser { + + private static final String STATE="state"; + + public String getParserName() { + return LatentStateBranchRateLengthStatistic.LATENT_STATE_BRANCH_RATE_LENGTH_STATISTIC; + } + + public Object parseXMLObject(XMLObject xo) throws XMLParseException { + SericolaLatentStateBranchRateModel latenBranchRateModel = (SericolaLatentStateBranchRateModel) xo.getChild(SericolaLatentStateBranchRateModel.class); + TreeModel tree = (TreeModel) xo.getChild(TreeModel.class); + LatentStateBranchRateLengthStatistic.STATE state = LatentStateBranchRateLengthStatistic.STATE.valueOf(xo.getAttribute(STATE, LatentStateBranchRateLengthStatistic.STATE.REPLICATING.name()).toUpperCase()); + return new LatentStateBranchRateLengthStatistic(latenBranchRateModel, tree, state); + } + public String getParserDescription() { + return "This element provides a statistic for the length of a tree spent in a replicating or latent state."; + } + + public Class getReturnType() { + return LatentStateBranchRateLengthStatistic.class; + } + + public XMLSyntaxRule[] getSyntaxRules() { + return rules; + } + + private XMLSyntaxRule[] rules = new XMLSyntaxRule[]{ + new ElementRule(SericolaLatentStateBranchRateModel.class, "A branch rate model to provide the rates for the non-latent state"), + new ElementRule(TreeModel.class, "The tree on which this will operate"), + }; + +} diff --git a/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java b/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateModelParser.java similarity index 96% rename from src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java rename to src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateModelParser.java index 8a216e72d7..89322e67cc 100644 --- a/src/dr/evomodelxml/branchratemodel/LatentStateBranchRateModelParser.java +++ b/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateModelParser.java @@ -25,7 +25,7 @@ * */ -package dr.evomodelxml.branchratemodel; +package dr.evomodelxml.branchratemodel.latentStateBranchRate; import dr.evomodel.branchratemodel.*; import dr.evomodel.tree.TreeModel; diff --git a/src/test/dr/evomodel/substmodel/TwoStateOccupancyMarkovRewardsTest.java b/src/test/dr/evomodel/substmodel/TwoStateOccupancyMarkovRewardsTest.java index a225ba317f..1b7b7fef8d 100644 --- a/src/test/dr/evomodel/substmodel/TwoStateOccupancyMarkovRewardsTest.java +++ b/src/test/dr/evomodel/substmodel/TwoStateOccupancyMarkovRewardsTest.java @@ -27,7 +27,7 @@ package test.dr.evomodel.substmodel; -import dr.evomodel.branchratemodel.LatentStateBranchRateModel; +import dr.evomodel.branchratemodel.latentStateBranchRate.LatentStateBranchRateModel; import dr.inference.markovjumps.MarkovReward; import dr.inference.markovjumps.SericolaSeriesMarkovReward; import dr.inference.markovjumps.TwoStateOccupancyMarkovReward; diff --git a/src/test/dr/math/GeneralizedIntegerGammaTest.java b/src/test/dr/math/GeneralizedIntegerGammaTest.java index 12c21d3df9..9dd3319dd0 100644 --- a/src/test/dr/math/GeneralizedIntegerGammaTest.java +++ b/src/test/dr/math/GeneralizedIntegerGammaTest.java @@ -27,7 +27,7 @@ package test.dr.math; -import dr.evomodel.branchratemodel.LatentStateBranchRateModel; +import dr.evomodel.branchratemodel.latentStateBranchRate.LatentStateBranchRateModel; import dr.inference.markovjumps.MarkovReward; import dr.inference.markovjumps.SericolaSeriesMarkovReward; import dr.inference.markovjumps.TwoStateOccupancyMarkovReward; From d1dcc793736454a2dd3f10be70425de19d5f56bb Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Tue, 21 Oct 2025 21:32:14 -0500 Subject: [PATCH 15/21] correcting parser path in properties --- src/dr/app/beast/development_parsers.properties | 3 --- src/dr/app/beast/release_parsers.properties | 3 ++- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/dr/app/beast/development_parsers.properties b/src/dr/app/beast/development_parsers.properties index f5bbcaa08e..f69738fe96 100644 --- a/src/dr/app/beast/development_parsers.properties +++ b/src/dr/app/beast/development_parsers.properties @@ -158,9 +158,6 @@ dr.evomodelxml.coalescent.demographicmodel.PowerLawGrowthModelParser dr.evomodelxml.coalescent.demographicmodel.PeakAndDeclineModelParser dr.evomodelxml.coalescent.demographicmodel.AsymptoticGrowthModelParser -#PARAMETER -dr.inferencexml.model.IndicatorProductParameterParser - # UNIFORM INTERNAL NODE HEIGHT PRIOR dr.evomodelxml.operators.FunkyPriorMixerOperatorParser diff --git a/src/dr/app/beast/release_parsers.properties b/src/dr/app/beast/release_parsers.properties index 7fb3eef235..022f9b565f 100644 --- a/src/dr/app/beast/release_parsers.properties +++ b/src/dr/app/beast/release_parsers.properties @@ -204,8 +204,9 @@ dr.evomodelxml.branchratemodel.AncestralTraitBranchRatesParser dr.evomodelxml.branchratemodel.RandomEffectsTreeTraitProviderParser ## Latent model -dr.evomodelxml.branchratemodel.LatentStateBranchRateModelParser +dr.evomodelxml.branchratemodel.latentStateBranchRate.LatentStateBranchRateModelParser dr.evomodelxml.branchratemodel.latentStateBranchRate.LatentStateBranchRateLengthStatisticParser +dr.inferencexml.model.IndicatorProductParameterParser #COALESCENT From 321dd39d9b419c087b928c94dc1c7396b26ab1a7 Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Thu, 12 Mar 2026 15:05:47 -0500 Subject: [PATCH 16/21] fixes copyright character --- src/dr/app/beast/development_parsers.properties | 2 +- src/dr/app/beast/release_parsers.properties | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dr/app/beast/development_parsers.properties b/src/dr/app/beast/development_parsers.properties index ca82024cb5..cd97e0f51d 100644 --- a/src/dr/app/beast/development_parsers.properties +++ b/src/dr/app/beast/development_parsers.properties @@ -1,7 +1,7 @@ # # development_parsers.properties # -# Copyright � 2002-2024 the BEAST Development Team +# Copyright © 2002-2024 the BEAST Development Team # http://beast.community/about # # This file is part of BEAST. diff --git a/src/dr/app/beast/release_parsers.properties b/src/dr/app/beast/release_parsers.properties index 8c8ba7e8ba..eec6387277 100644 --- a/src/dr/app/beast/release_parsers.properties +++ b/src/dr/app/beast/release_parsers.properties @@ -1,7 +1,7 @@ # # release_parsers.properties # -# Copyright � 2002-2024 the BEAST Development Team +# Copyright © 2002-2024 the BEAST Development Team # http://beast.community/about # # This file is part of BEAST. From 4fe5fba3ba537b2716a904e65c12a40900955d5a Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Thu, 12 Mar 2026 16:06:33 -0500 Subject: [PATCH 17/21] adds short label arguments to MRDC (and builds) --- src/dr/app/tools/MarkovRewardDensityCalculator.java | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/dr/app/tools/MarkovRewardDensityCalculator.java b/src/dr/app/tools/MarkovRewardDensityCalculator.java index d96088bfb9..71c657e948 100644 --- a/src/dr/app/tools/MarkovRewardDensityCalculator.java +++ b/src/dr/app/tools/MarkovRewardDensityCalculator.java @@ -111,10 +111,10 @@ void run(double time){ public static void main(String[] args) throws IOException, Arguments.ArgumentException { Arguments arguments = new Arguments( new Arguments.Option[]{ - new Arguments.RealOption("rate", "Rate parameter"), - new Arguments.RealOption("bias", "Bias parameter"), - new Arguments.RealOption("length", "branchlength to evaluate"), - new Arguments.Option("help", "option to print this message"), + new Arguments.RealOption("rate","r", "Rate parameter"), + new Arguments.RealOption("bias","b", "Bias parameter"), + new Arguments.RealOption("length","l", "branchlength to evaluate"), + new Arguments.Option("help", "h","option to print this message"), } ); From 65f0f510c41d1ede66e8e44c3826f1b3661ba0dd Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Fri, 13 Mar 2026 09:42:46 -0500 Subject: [PATCH 18/21] moves indicators to the branch rate model --- .../SericolaLatentStateBranchRateModel.java | 18 ++++++++++++++---- .../LatentStateBranchRateModelParser.java | 8 +++++++- 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java index 94bea79f8a..bea537c263 100644 --- a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java +++ b/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java @@ -73,6 +73,7 @@ public class SericolaLatentStateBranchRateModel extends AbstractModelLikelihood private final Parameter latentTransitionRateParameter; private final Parameter latentTransitionFrequencyParameter; private final TreeParameterModel latentStateProportions; + private final TreeParameterModel latentStateIndicators; private final Parameter latentStateProportionParameter; private final CountableBranchCategoryProvider branchCategoryProvider; @@ -100,6 +101,7 @@ public SericolaLatentStateBranchRateModel(String name, Parameter latentTransitionRateParameter, Parameter latentTransitionFrequencyParameter, Parameter latentStateProportionParameter, + Parameter latentStateIndicatorParameter, CountableBranchCategoryProvider branchCategoryProvider, double epsilon, boolean excludeRoot) { @@ -121,12 +123,15 @@ public SericolaLatentStateBranchRateModel(String name, if (branchCategoryProvider == null) { this.latentStateProportions = new TreeParameterModel(tree, latentStateProportionParameter, false, Intent.BRANCH); + this.latentStateIndicators = new TreeParameterModel(tree, latentStateIndicatorParameter, false, Intent.BRANCH); addModel(latentStateProportions); + addModel(latentStateIndicators); this.latentStateProportionParameter = null; this.branchCategoryProvider = null; } else { this.latentStateProportions = null; + this.latentStateIndicators = null; this.branchCategoryProvider = branchCategoryProvider; this.latentStateProportionParameter = latentStateProportionParameter; this.latentStateProportionParameter.setDimension(branchCategoryProvider.getCategoryCount()); @@ -157,8 +162,9 @@ public SericolaLatentStateBranchRateModel(String name, Parameter latentTransitionRateParameter, Parameter latentTransitionFrequencyParameter, Parameter latentStateProportionParameter, + Parameter latentStateIndicatorParameter, CountableBranchCategoryProvider branchCategoryProvider){ - this(name, treeModel, nonLatentRateModel, latentTransitionRateParameter, latentTransitionFrequencyParameter, latentStateProportionParameter, branchCategoryProvider, 1E-10,false); + this(name, treeModel, nonLatentRateModel, latentTransitionRateParameter, latentTransitionFrequencyParameter, latentStateProportionParameter, latentStateIndicatorParameter,branchCategoryProvider, 1E-10,false); } public SericolaLatentStateBranchRateModel(Parameter rate, Parameter prop) { super(LATENT_STATE_BRANCH_RATE_MODEL); @@ -168,6 +174,7 @@ public SericolaLatentStateBranchRateModel(Parameter rate, Parameter prop) { latentTransitionFrequencyParameter = prop; latentStateProportions = null; this.latentStateProportionParameter = null; + this.latentStateIndicators = null; this.branchCategoryProvider = null; } @@ -206,7 +213,7 @@ public double getBranchRate(Tree tree, NodeRef node) { public double getLatentProportion(Tree tree, NodeRef node) { if (latentStateProportions != null) { - return latentStateProportions.getNodeValue(tree, node); + return latentStateProportions.getNodeValue(tree, node)*latentStateIndicators.getNodeValue(tree,node); } else { return latentStateProportionParameter.getParameterValue(branchCategoryProvider.getBranchCategory(tree, node)); } @@ -229,7 +236,7 @@ protected void handleModelChangedEvent(Model model, Object object, int index) { } else if (model == nonLatentRateModel) { // rates will change but the latent proportions haven't so the density is unchanged - } else if (model == latentStateProportions) { + } else if (model == latentStateProportions || model ==latentStateIndicators) { likelihoodKnown = false; // argument of density has changed if (index == -1) { @@ -470,7 +477,10 @@ public TreeTrait getTreeTrait(final String key) { return this; } else if (latentStateProportions != null && key.equals(latentStateProportions.getTraitName())) { return latentStateProportions; - } else if (branchCategoryProvider != null && key.equals(branchCategoryProvider.getTraitName())) { + }else if(latentStateIndicators!=null && key.equals(latentStateIndicators.getTraitName())){ + return latentStateIndicators; + } + else if (branchCategoryProvider != null && key.equals(branchCategoryProvider.getTraitName())) { return branchCategoryProvider; } else { throw new IllegalArgumentException("Unrecognised Tree Trait key, " + key); diff --git a/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateModelParser.java b/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateModelParser.java index 89322e67cc..255b3299dc 100644 --- a/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateModelParser.java +++ b/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateModelParser.java @@ -40,6 +40,7 @@ public class LatentStateBranchRateModelParser extends AbstractXMLObjectParser { public static final String LATENT_TRANSITION_RATE = "latentTransitionRate"; public static final String LATENT_TRANSITION_FREQUENCY = "latentTransitionFrequency"; public static final String LATENT_STATE_PROPORTIONS = "latentStateProportions"; + public static final String LATENT_STATE_INDICATORS = "latentStateIndicators"; public static final String EPSILON = "epsilon"; public static final String EXCLUDE_ROOT = "excludeRoot"; @@ -58,9 +59,14 @@ public Object parseXMLObject(XMLObject xo) throws XMLParseException { CountableBranchCategoryProvider branchCategoryProvider = (CountableBranchCategoryProvider)xo.getChild(CountableBranchCategoryProvider.class); Parameter latentStateProportionParameter = null; + if (xo.hasChildNamed(LATENT_STATE_PROPORTIONS)) { latentStateProportionParameter = (Parameter) xo.getElementFirstChild(LATENT_STATE_PROPORTIONS); } + Parameter latentStateIndicatorParameter = null; + if (xo.hasChildNamed(LATENT_STATE_INDICATORS)) { + latentStateIndicatorParameter = (Parameter) xo.getElementFirstChild(LATENT_STATE_INDICATORS); + } double epsilon = xo.getAttribute(EPSILON, 1E-10); Logger.getLogger("dr.evomodel").info("\nCreating a latent state branch rate model"); @@ -71,7 +77,7 @@ public Object parseXMLObject(XMLObject xo) throws XMLParseException { return new SericolaLatentStateBranchRateModel(SericolaLatentStateBranchRateModel.LATENT_STATE_BRANCH_RATE_MODEL, tree, nonLatentRateModel, latentTransitionRateParameter, latentTransitionFrequencyParameter, /* 0/1 CTMC have two parameters */ - latentStateProportionParameter, branchCategoryProvider,epsilon,excludeRoot); + latentStateProportionParameter,latentStateIndicatorParameter, branchCategoryProvider,epsilon,excludeRoot); } //************************************************************************ From 0de248e5ece343540a65c0e511473c56549f0da9 Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Fri, 13 Mar 2026 09:48:17 -0500 Subject: [PATCH 19/21] cleans up packaging for latentBranchRateModel --- .../LatentStateBranchRateLengthStatistic.java | 1 - .../latentStateBranchRate/LatentStateBranchRateModel.java | 2 ++ .../SericolaLatentStateBranchRateModel.java | 5 +++-- .../LatentStateBranchRateLengthStatisticParser.java | 2 +- .../LatentStateBranchRateModelParser.java | 1 + 5 files changed, 7 insertions(+), 4 deletions(-) rename src/dr/evomodel/branchratemodel/{ => latentStateBranchRate}/SericolaLatentStateBranchRateModel.java (99%) diff --git a/src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatistic.java b/src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatistic.java index bc209ec591..5bf6b8e0da 100644 --- a/src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatistic.java +++ b/src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatistic.java @@ -1,7 +1,6 @@ package dr.evomodel.branchratemodel.latentStateBranchRate; import dr.evolution.tree.Tree; -import dr.evomodel.branchratemodel.SericolaLatentStateBranchRateModel; import dr.evomodel.tree.TreeModel; import dr.evomodel.tree.TreeStatistic; import dr.inference.model.Bounds; diff --git a/src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateModel.java b/src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateModel.java index df56d87a4c..31a1cdc5ef 100644 --- a/src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateModel.java +++ b/src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateModel.java @@ -52,6 +52,8 @@ * $LastChangedDate$ * $LastChangedRevision$ */ +// Deprecated in favor of SericolaLatentStateBranchRateModel +@Deprecated public class LatentStateBranchRateModel extends AbstractModelLikelihood implements BranchRateModel { public static final String LATENT_STATE_BRANCH_RATE_MODEL = "latentStateBranchRateModel"; diff --git a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java b/src/dr/evomodel/branchratemodel/latentStateBranchRate/SericolaLatentStateBranchRateModel.java similarity index 99% rename from src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java rename to src/dr/evomodel/branchratemodel/latentStateBranchRate/SericolaLatentStateBranchRateModel.java index bea537c263..33ef95c21a 100644 --- a/src/dr/evomodel/branchratemodel/SericolaLatentStateBranchRateModel.java +++ b/src/dr/evomodel/branchratemodel/latentStateBranchRate/SericolaLatentStateBranchRateModel.java @@ -25,15 +25,16 @@ * */ -package dr.evomodel.branchratemodel; +package dr.evomodel.branchratemodel.latentStateBranchRate; import dr.evolution.tree.NodeRef; import dr.evolution.tree.Tree; import dr.evolution.tree.TreeTrait; +import dr.evomodel.branchratemodel.BranchRateModel; +import dr.evomodel.branchratemodel.CountableBranchCategoryProvider; import dr.evomodel.tree.TreeModel; import dr.evomodel.tree.TreeParameterModel; import dr.inference.markovjumps.MarkovReward; -import dr.inference.markovjumps.TwoStateOccupancyMarkovReward; import dr.inference.markovjumps.TwoStateSericolaSeriesMarkovReward; import dr.inference.model.AbstractModelLikelihood; import dr.inference.model.Model; diff --git a/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatisticParser.java b/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatisticParser.java index 3a09cf9d3d..7da535e11f 100644 --- a/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatisticParser.java +++ b/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatisticParser.java @@ -1,8 +1,8 @@ package dr.evomodelxml.branchratemodel.latentStateBranchRate; -import dr.evomodel.branchratemodel.SericolaLatentStateBranchRateModel; import dr.evomodel.branchratemodel.latentStateBranchRate.LatentStateBranchRateLengthStatistic; +import dr.evomodel.branchratemodel.latentStateBranchRate.SericolaLatentStateBranchRateModel; import dr.evomodel.tree.TreeModel; import dr.xml.AbstractXMLObjectParser; import dr.xml.ElementRule; diff --git a/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateModelParser.java b/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateModelParser.java index 255b3299dc..ee77715145 100644 --- a/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateModelParser.java +++ b/src/dr/evomodelxml/branchratemodel/latentStateBranchRate/LatentStateBranchRateModelParser.java @@ -28,6 +28,7 @@ package dr.evomodelxml.branchratemodel.latentStateBranchRate; import dr.evomodel.branchratemodel.*; +import dr.evomodel.branchratemodel.latentStateBranchRate.SericolaLatentStateBranchRateModel; import dr.evomodel.tree.TreeModel; import dr.inference.model.Parameter; import dr.xml.*; From f719744b46195ce4f0b28e5f11d8629d21c404f5 Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Fri, 13 Mar 2026 09:51:37 -0500 Subject: [PATCH 20/21] Removing indicator product parameter as it is too easy to abuse --- src/dr/app/beast/release_parsers.properties | 2 +- .../model/IndicatorProductParameter.java | 112 ------------------ .../IndicatorProductParameterParser.java | 51 -------- 3 files changed, 1 insertion(+), 164 deletions(-) delete mode 100644 src/dr/inference/model/IndicatorProductParameter.java delete mode 100644 src/dr/inferencexml/model/IndicatorProductParameterParser.java diff --git a/src/dr/app/beast/release_parsers.properties b/src/dr/app/beast/release_parsers.properties index eec6387277..3aab17d598 100644 --- a/src/dr/app/beast/release_parsers.properties +++ b/src/dr/app/beast/release_parsers.properties @@ -206,7 +206,7 @@ dr.evomodelxml.branchratemodel.RandomEffectsTreeTraitProviderParser ## Latent model dr.evomodelxml.branchratemodel.latentStateBranchRate.LatentStateBranchRateModelParser dr.evomodelxml.branchratemodel.latentStateBranchRate.LatentStateBranchRateLengthStatisticParser -dr.inferencexml.model.IndicatorProductParameterParser + #COALESCENT diff --git a/src/dr/inference/model/IndicatorProductParameter.java b/src/dr/inference/model/IndicatorProductParameter.java deleted file mode 100644 index e28cb0a8c4..0000000000 --- a/src/dr/inference/model/IndicatorProductParameter.java +++ /dev/null @@ -1,112 +0,0 @@ -package dr.inference.model; - -import java.util.List; - -/** - * This is parameter class that is intended to allow sampling values on the range - * [0,) inclusive of 0. It represents a specific case of a produce parameter. - - */ -public class IndicatorProductParameter extends Parameter.Abstract implements VariableListener { - Parameter indicatorParameter; - Parameter continuousParameter; - - private Bounds bounds = null; - - public IndicatorProductParameter(Parameter indicatorParameter, Parameter continuousParameter) { - - this.indicatorParameter = indicatorParameter; - this.continuousParameter = continuousParameter; - - indicatorParameter.addVariableListener(this); - Parameter.CONNECTED_PARAMETER_SET.add(indicatorParameter); - - continuousParameter.addVariableListener(this); - Parameter.CONNECTED_PARAMETER_SET.add(continuousParameter); - } - public int getDimension() { - return indicatorParameter.getDimension(); - } - @Override - public boolean isImmutable() { - return true; - } - - protected void storeValues() { - indicatorParameter.storeParameterValues(); - continuousParameter.storeParameterValues(); - } - protected void restoreValues() { - indicatorParameter.restoreParameterValues(); - continuousParameter.restoreParameterValues(); - } - protected void acceptValues() { - indicatorParameter.acceptParameterValues(); - continuousParameter.acceptParameterValues(); - } - - public double getParameterValue(int dim) { - double value = indicatorParameter.getParameterValue(dim) *continuousParameter.getParameterValue(dim); - return value; - } - - public void setParameterValue(int dim, double value) { - if(value==0){ - indicatorParameter.setParameterValue(dim, 0); - }else{ - indicatorParameter.setParameterValue(dim, 1); - continuousParameter.setParameterValue(dim, value); - } - } - - public void setParameterValueQuietly(int dim, double value) { - throw new RuntimeException("Not implemented"); - } - - public void setParameterValueNotifyChangedAll(int dim, double value){ - throw new RuntimeException("Not implemented"); - } - @Override - protected void adoptValues(Parameter source) { - // TODO Auto-generated method stub - throw new UnsupportedOperationException("Unimplemented method 'adoptValues'"); - } - - public String getParameterName() { - if (getId() == null) { - StringBuilder sb = new StringBuilder("IndicatorProduct"); - - sb.append(".").append(indicatorParameter.getId()); - sb.append(".").append(continuousParameter.getId()); - - setId(sb.toString()); - } - return getId(); - } - public void addBounds(Bounds bounds) { - this.bounds = bounds; - } - - public Bounds getBounds() { - if (bounds == null) { - return indicatorParameter.getBounds(); // TODO - } else { - return bounds; - } - } - - public void addDimension(int index, double value) { - throw new RuntimeException("Not yet implemented."); - } - - public double removeDimension(int index) { - throw new RuntimeException("Not yet implemented."); - } - - public void variableChangedEvent(Variable variable, int index, ChangeType type) { - fireParameterChangedEvent(index,type); - } - - - -} diff --git a/src/dr/inferencexml/model/IndicatorProductParameterParser.java b/src/dr/inferencexml/model/IndicatorProductParameterParser.java deleted file mode 100644 index 0fc14f0f5c..0000000000 --- a/src/dr/inferencexml/model/IndicatorProductParameterParser.java +++ /dev/null @@ -1,51 +0,0 @@ -package dr.inferencexml.model; - -import dr.inference.model.IndicatorProductParameter; -import dr.inference.model.Parameter; - -import dr.xml.*; - - -public class IndicatorProductParameterParser extends AbstractXMLObjectParser { - - public static final String INDICATOR_PRODUCT_PARAMETER = "indicatorProductParameter"; - public static final String INDICATOR = "indicator"; - public static final String CONTINUOUS_PARAMETER = "continuousParameter"; - - public Object parseXMLObject(XMLObject xo) throws XMLParseException { - - Parameter continuous = (Parameter) xo.getChild(CONTINUOUS_PARAMETER).getChild(Parameter.class); - Parameter indicator = (Parameter) xo.getChild(INDICATOR).getChild(Parameter.class); - - IndicatorProductParameter prodParam = new IndicatorProductParameter(indicator, continuous); - - return prodParam; - - - } - - public XMLSyntaxRule[] getSyntaxRules() { - return rules; - } - - private final XMLSyntaxRule[] rules = { - new ElementRule(INDICATOR, new XMLSyntaxRule[]{ - new ElementRule(Parameter.class) - }, false), - new ElementRule(CONTINUOUS_PARAMETER, new XMLSyntaxRule[]{ - new ElementRule(Parameter.class) - }, false) - }; - - public String getParserDescription() { - return "A element-wise product of an indicator and continuous paramter."; - } - - public Class getReturnType() { - return Parameter.class; - } - - public String getParserName() { - return INDICATOR_PRODUCT_PARAMETER; - } -} From a5dbe6fdb6c58a68366dd7c5731f4f7332438e49 Mon Sep 17 00:00:00 2001 From: jtmccr1 Date: Fri, 13 Mar 2026 11:15:23 -0500 Subject: [PATCH 21/21] adds reward density test for latent model classes --- .../LatentRewardDensityTest.java | 146 ++++++++++++++++++ 1 file changed, 146 insertions(+) create mode 100644 src/test/dr/evomodel/branchratemodel/latentStateBranchRate/LatentRewardDensityTest.java diff --git a/src/test/dr/evomodel/branchratemodel/latentStateBranchRate/LatentRewardDensityTest.java b/src/test/dr/evomodel/branchratemodel/latentStateBranchRate/LatentRewardDensityTest.java new file mode 100644 index 0000000000..cde5beb195 --- /dev/null +++ b/src/test/dr/evomodel/branchratemodel/latentStateBranchRate/LatentRewardDensityTest.java @@ -0,0 +1,146 @@ +package test.dr.evomodel.branchratemodel.latentStateBranchRate; + +import dr.inference.markovjumps.SericolaSeriesMarkovReward; +import dr.inference.markovjumps.TwoStateOccupancyMarkovReward; +import dr.inference.markovjumps.TwoStateSericolaSeriesMarkovReward; +import junit.framework.TestCase; +/** + * @author JT McCrone + * This class tests the various classes that estimate the reward for a two state + * model. The TwoStateSericolaSeriesMarkovReward is used in the SericolaLatentStateBranchModel. + * TwoStateOccupancyMarkovReward is used in the deprecated LatentStateBranchRateModel + * We compare the cdf, pdf, and conditional probabilities of these classes where applicable. + */ +public class LatentRewardDensityTest extends TestCase{ + + private final int dim = 2; + private final double[] r = new double[]{0.0,1.0}; // rewards + private final double epsilon = 1e-5; + private double[] Q(double rate, double bias){ + return new double[]{ + -rate * bias, rate * bias, + rate * (1.0 - bias), -rate * (1.0 - bias) + }; + } + private final double totalTime = 100.0; + + public void testPdfEqualRates(){ + double rate = 0.05; + double bias = 0.5; + TwoStateOccupancyMarkovReward TSOMR = createTwoStateOccupancyMarkovReward(rate, bias); + TwoStateSericolaSeriesMarkovReward TSSMR = createTwoStateSericolaSeriesMarkovReward(rate, bias); + SericolaSeriesMarkovReward SSMR = createSericolaSeriesMarkovReward(rate, bias); + + // pdf + double[] latentTimes = new double[10]; + + for(int i=0; i