Merge work from Latent model

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.

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,12 @@ dr.evomodelxml.branchratemodel.BranchRateGradientWrtIncrementsParser
 dr.evomodelxml.branchratemodel.AncestralTraitBranchRatesParser
 dr.evomodelxml.branchratemodel.RandomEffectsTreeTraitProviderParser
 
+## Latent model
+dr.evomodelxml.branchratemodel.latentStateBranchRate.LatentStateBranchRateModelParser
+dr.evomodelxml.branchratemodel.latentStateBranchRate.LatentStateBranchRateLengthStatisticParser
+
+
+
 #COALESCENT
 dr.evomodelxml.coalescent.demographicmodel.CataclysmicDemographicModelParser
 dr.evomodelxml.coalescent.demographicmodel.ExpConstExpDemographicModelParser

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<times.length; i++){
+            times[i] = time * i / 100.0;
+
+            System.out.print("{");
+
+//TODO remove duplication
+                    pdf = twoStateOccupancyMarkovReward.computePdf(times[i], time, 0, 0);
+                    // cdf = twoStateOccupancyMarkovReward.computeCdf(times[i], time, 0, 0);
+                    conditional = twoStateOccupancyMarkovReward.computeConditionalProbability( time, 0, 0);
+
+                    System.out.print("\"time\": " + times[i] + ",");
+                    System.out.print("\"totalTime\": " + time + ",");
+                    // System.out.print("\"twoStateOccupancyMarkovReward\": {");
+                    System.out.print("\"pdf\": " + pdf + ",");
+                    // System.out.print("\"cdf\": " + cdf + ",");
+                    System.out.print("\"conditional\": " + conditional + ",");
+                    System.out.print("\"implementation\": \"twoStateOccupancyMarkovReward\"");
+
+            System.out.print("},");
+
+            System.out.print("{");
+
+                    pdf = twoStateSericolaSeriesMarkovReward.computePdf(times[i], time, 0, 0);
+                    cdf = twoStateSericolaSeriesMarkovReward.computeCdf(times[i], time, 0, 0);
+                    conditional = twoStateSericolaSeriesMarkovReward.computeConditionalProbability(time, 0, 0);
+                    System.out.print("\"time\": " + times[i] + ",");
+                    System.out.print("\"totalTime\": " + time + ",");
+                    System.out.print("\"pdf\": " + pdf + ",");
+                    System.out.print("\"cdf\": " + cdf + ",");
+                    System.out.print("\"conditional\": " + conditional + ",");
+                    System.out.print("\"implementation\": \"twoStateSericolaSeriesMarkovReward\"");
+
+            System.out.print("},");
+
+            System.out.print("{");
+                    pdf = sericolaSeriesMarkovReward.computePdf(times[i], time, 0, 0);
+                    cdf = sericolaSeriesMarkovReward.computeCdf(times[i], time, 0, 0);
+                    conditional = sericolaSeriesMarkovReward.computeConditionalProbability(time, 0, 0);
+
+                    System.out.print("\"time\": " + times[i] + ",");
+                    System.out.print("\"totalTime\": " + time + ",");
+                    System.out.print("\"pdf\": " + pdf + ",");
+                    System.out.print("\"cdf\": " + cdf + ",");
+                    System.out.print("\"conditional\": " + conditional + ",");
+                    System.out.print("\"implementation\": \"sericolaSeriesMarkovReward\"");
+
+                if(i < times.length - 1){
+                    System.out.print("},");
+                }else{
+                    System.out.print("}");
+                }
+
+        }
+
+        System.out.println("]}");
+
+
+
+    }
+
+public static void main(String[] args) throws IOException, Arguments.ArgumentException {
+    Arguments arguments = new Arguments(
+        new Arguments.Option[]{
+            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"),
+        }
+    );
+
+
+    try {
+        arguments.parseArguments(args);
+    } catch (Arguments.ArgumentException ae) {
+        System.err.println(ae);
+        System.exit(1);
+    }
+
+    double rate = arguments.getRealOption("rate");
+    double bias = arguments.getRealOption("bias");
+    double length = arguments.getRealOption("length");
+
+    MarkovRewardDensityCalculator calculator = new MarkovRewardDensityCalculator(rate, bias);
+
+    calculator.run(length);
+
+
+}
+
+
+private class Result{
+    double pdf;
+    double cdf;
+    double length;
+    double s;
+    double rate;
+    double bias;
+    String calculator;
+
+
+public Result(double pdf, double cdf, double length, double s, double rate, double bias, String calculator) {
+    this.pdf = pdf;
+    this.cdf = cdf;
+    this.length = length;
+    this.s = s;
+    this.rate = rate;
+    this.bias = bias;
+    this.calculator = calculator;
+}
+}
+
+}

src/dr/evomodel/branchratemodel/latentStateBranchRate/LatentStateBranchRateLengthStatistic.java

@@ -0,0 +1,126 @@
+package dr.evomodel.branchratemodel.latentStateBranchRate;
+
+import dr.evolution.tree.Tree;
+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<Double> 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;
+
+}

src/dr/evomodel/branchratemodel/LatentStateBranchRateModel.java → 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;
@@ -50,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";