give warning if branchRates are given and tree has branchRates #491

LinguaPhylo · Jun 6, 2024 · 07ef100 · 07ef100
1 parent c84f649
commit 07ef100
Show file tree

Hide file tree

Showing 2 changed files with 49 additions and 5 deletions.
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;
@@ -66,11 +67,15 @@ public AbstractPhyloCTMC(Value<TimeTree> tree, Value<Number> clockRate, Value<Do
         this.random = RandomUtils.getRandom();
 
         Double[] treeBranchRates = tree.value().getBranchRates();
-        // if tree has branch rates, then use them
-        if (branchRates == null && treeBranchRates != null && treeBranchRates.length > 0) {
-            this.branchRates = new Value<>("branchRates", treeBranchRates);
-        }
 
+        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();
     }
 

diff --git a/lphy-base/src/test/java/lphy/base/evolution/likelihood/PhyloCTMCTest.java b/lphy-base/src/test/java/lphy/base/evolution/likelihood/PhyloCTMCTest.java
@@ -22,7 +22,6 @@ void setUp() {
         newickTree = "((1:2.0, (2:1.0, 3:1.0):1.0):2.0, 4:4.0)";
     }
 
-
     @Test
     void apply() {
         // generate a tree with local clock
@@ -97,4 +96,44 @@ void apply() {
         // index 6 : node (((2,3),1),4)
         assertEquals(allNodes.get(6).getBranchRate(), branchRates[6]);
     }
+
+    @Test
+    void testWarning() {
+        TimeTree tree = Newick.parseNewick(newickTree);
+        for (TimeTreeNode node: tree.getNodes()){
+            node.setBranchRate(0.5);
+        }
+
+        Double[][] Q = {
+                { -1.0,  0.5,  0.3,  0.2 },
+                {  0.4, -1.0,  0.1,  0.5 },
+                {  0.3,  0.2, -1.0,  0.5 },
+                {  0.2,  0.3,  0.5, -1.0 }
+        };
+
+        Value<TimeTree> newTreeValue = new Value<>("newTree", tree);
+        Value<Double[][]> QValue = new Value<>("Q", Q);
+        Value<Double[]> siteRatesValue = new Value<Double[]>("siteRates", new Double[]{0.1, 0.1, 0.1, 0.1, 0.1});
+        Value<Integer> LValue = new Value<Integer>("L", 5);
+        Value<Double[]> branchRatesValue = new Value<>("branchRates", new Double[]{0.2,0.2,0.2,0.2,0.2,0.2,0.2});
+
+        PhyloCTMC phyloCTMCInstance = new PhyloCTMC(
+                newTreeValue,null,null,
+                QValue,siteRatesValue,branchRatesValue,
+                LValue,null, null);
+
+        phyloCTMCInstance.sample();
+        // should have a warning msg
+
+        Double[] branchRates = phyloCTMCInstance.getBranchRates().value();
+
+        List<TimeTreeNode> allNodes = tree.getNodes();
+        assertEquals(allNodes.size(), branchRates.length);
+
+        // should use given branch rates and tree branch rates should not be changed
+        for (int i = 0; i<allNodes.size(); i++){
+            assertEquals(0.2, branchRates[i]);
+            assertEquals(0.5, allNodes.get(i).getBranchRate());
+        }
+    }
 }