prior problem in birthDeathOnRootAgeAndTaxa.lphy #492 #496

LinguaPhylo · Jun 14, 2024 · c58147c · c58147c
1 parent 2b1fd6e
commit c58147c
Show file tree

Hide file tree

Showing 2 changed files with 22 additions and 7 deletions.
diff --git a/examples/birth-death/birthDeathOnRootAgeAndTaxa.lphy b/examples/birth-death/birthDeathOnRootAgeAndTaxa.lphy
@@ -17,8 +17,8 @@ model{
   // Conditioned on root age and on number of taxa
   ψ ~ BirthDeath(lambda=λ, mu=death_rate, rootAge=root_time, taxa=taxa);
 
-  ucln_mean ~ Exp(mean=2.0);
-  ucln_sigma ~ Exp(mean=3.0);
+  ucln_mean ~ Normal(mean=-4.5, sd=0.5);
+  ucln_sigma ~ Exp(mean=0.5);
   branch_rates ~ LogNormal(meanlog=ucln_mean, sdlog=ucln_sigma, replicates=2*n - 2);
 
   D ~ PhyloCTMC(tree=ψ, Q=Q, branchRates=branch_rates, L=L);

diff --git a/lphy-base/src/main/java/lphy/base/evolution/likelihood/AbstractPhyloCTMC.java b/lphy-base/src/main/java/lphy/base/evolution/likelihood/AbstractPhyloCTMC.java
@@ -5,6 +5,7 @@
 import lphy.base.evolution.alignment.SimpleAlignment;
 import lphy.base.evolution.tree.TimeTree;
 import lphy.base.evolution.tree.TimeTreeNode;
+import lphy.core.logger.LoggerUtils;
 import lphy.core.model.GenerativeDistribution;
 import lphy.core.model.Value;
 import lphy.core.simulator.RandomUtils;
@@ -47,11 +48,15 @@ public abstract class AbstractPhyloCTMC implements GenerativeDistribution<Alignm
     protected Value<Double[]> rootFreqs;
     protected SortedMap<String, Integer> idMap = new TreeMap<>();
     protected double[][] transProb;
+    /**
+     * <code>e^{Qt} = Ee^{At}E^-1</code>, where A is a diagonal matrix of eigenvalues (Eval),
+     * E is the matrix of right eigenvectors (Evec), and E^-1 is the matrix of left eigenvectors (Ievc).
+     */
     private EigenDecomposition decomposition;
-    private double[][] Ievc;
-    private double[][] Evec;
-    private double[][] iexp;
-    private double[] Eval;
+    private double[][] Ievc; // inverse Eigen vectors
+    private double[][] Evec; // Eigen vectors
+    private double[][] iexp; // intermediate matrix
+    private double[] Eval; // Eigenvalues
 
 
     public AbstractPhyloCTMC(Value<TimeTree> tree, Value<Number> clockRate, Value<Double[]> freq,
@@ -65,7 +70,17 @@ public AbstractPhyloCTMC(Value<TimeTree> tree, Value<Number> clockRate, Value<Do
 
         this.random = RandomUtils.getRandom();
 
-//        checkCompatibilities();
+        Double[] treeBranchRates = tree.value().getBranchRates();
+
+        if (treeBranchRates != null && treeBranchRates.length > 0) {
+            if (this.branchRates != null) { // have branchRates from input but tree also has branch rates
+                LoggerUtils.log.warning("PhyloCTMC has branchRates from input parameter and tree has branch rates, " +
+                        "default to using input parameter branchRates.");
+            } else { // if tree has branch rates, then use them
+                this.branchRates = new Value<>("branchRates", treeBranchRates);
+            }
+        }
+        //        checkCompatibilities();
     }
 
     //+++ protected methods +++//