Refactor PhyloDiscrete framework and add Negative Binomial CNV error model

phylonco-beast/version.xml

@@ -36,6 +36,8 @@
         <provider classname="phylonco.beast.evolution.readcountmodel.GibbsAlignmentOperator"/>
         <provider classname="phylonco.beast.evolution.datatype.ReadCount"/>
         <provider classname="NestedBD.evolution.substitutionmodel.BD"/>
+        <provider classname="NestedBD.evolution.errormodel.NegativeBinomialErrorModel"/>
+
     </service>
 
 </package>

phylonco-lphy/examples/CopyNumber_Yule.lphy

@@ -0,0 +1,16 @@
+//data
+  D = readCopyProfile(file = "examples/data/copyNumberSim_ginkgo.txt");
+  nsites = nchar(D);
+  taxa_list = D.taxa();
+
+//tree prior (Yule)
+λ ~ LogNormal(meanlog=1.0, sdlog=1.5);
+ψ ~ Yule(lambda=λ, taxa = taxa_list);
+
+// Relaxed clock
+clockSD ~ LogNormal(meanlog=-1.5, sdlog=0.75);
+branchRates ~ LogNormal(meanlog=0.0, sdlog=clockSD, replicates=ψ.branchCount());
+
+//model
+cnvModel = CopyNumberBD(lambda=1.0, mu=1.0, rootState=2);
+D ~ PhyloDiscrete(evolutionModel=cnvModel, tree=ψ, L=nsites, branchRates=branchRates);

phylonco-lphy/src/main/java/phylonco/lphy/evolution/copynumbermodel/CopyNumberBD.java

@@ -5,32 +5,77 @@
 import lphy.core.model.annotation.ParameterInfo;
 import lphy.core.simulator.RandomUtils;
 import org.apache.commons.math3.random.RandomGenerator;
-import java.util.Map;
 
+import java.util.HashMap;
+import java.util.Map;
 
+/**
+ * Copy Number Birth-Death Model for phylogenetic inference of copy number variation.
+ *
+ * <p>Models copy number evolution as a birth-death process where copy numbers can
+ * gain (birth) or loss (death) by one unit per event.
+ * State 0 is absorbing and cannot be recovered.</p>
+ * <p>PARAMETERIZATION: Can be specified in two ways:
+ * <ul>
+ *   <li>Single rate: bdRate (assumes lambda=mu)</li>
+ *   <li>Separate rates: lambda and mu (for future flexibility)</li>
+ * </ul>
+ * NOTE: Internally, lambda and mu are ALWAYS used for simulation.
+ * When bdRate is provided, lambda and mu are set equal to bdRate.</p>
+ *
+ * <p>Used with {@link PhyloDiscrete} for simulating CNV data along phylogenetic trees.</p>
+ */
 public class CopyNumberBD extends DeterministicFunction<MarkovTraitEvolution<Integer>> implements MarkovTraitEvolution<Integer> {
     private Value<Double> lambda;
     private Value<Double> mu;
+    private Value<Double> bdRate;
     private Value<Integer> rootState;
+    private Value<Integer> nstate;
 
     public static final String birthRateString = "lambda";
     public static final String deathRateString = "mu";
+    public static final String bdRateString = "bdRate";
     public static final String rootStateString = "rootState";
+    public static final String nstateString = "nstate";
 
-    // Constructor with all parameters including optional rootState
+    // Constructor 1: Using single bdRate (lambda=mu)
     public CopyNumberBD(
-            @ParameterInfo(name = birthRateString, description = "The birth rate >= 0 (generation) of copy number variants") Value<Double> lambda,
-            @ParameterInfo(name = deathRateString, description = "The death rate >= 0 (loss) of copy number variants") Value<Double> mu,
-            @ParameterInfo(name = rootStateString, description = "The initial copy number state at the root node") Value<Integer> rootState
+            @ParameterInfo(name = bdRateString, description = "Birth-death rate when we assume lambda=mu")
+            Value<Double> bdRate,
+            @ParameterInfo(name = rootStateString, description = "Initial copy number at root", optional = true)
+            Value<Integer> rootState,
+            @ParameterInfo(name = nstateString, description = "Maximum number of copy number states (default=15)", optional = true)
+            Value<Integer> nstate
     ) {
         super();
-        this.lambda = lambda;
-        this.mu = mu;
+        assert (bdRate.value() >= 0) : "bdRate must be non-negative";
         this.rootState = rootState != null ? rootState : new Value<>("rootState", 2);
+        this.nstate = nstate != null ? nstate : new Value<>("nstate", 15);
+        this.bdRate = bdRate;
+        this.lambda = new Value<>("lambda", bdRate.value());
+        this.mu = new Value<>("mu", bdRate.value());
+    }
 
-        // Birth rate(lambda) and Death rate(mu) are non-negative
-        assert (lambda.value() >= 0);
-        assert (mu.value() >= 0);
+    // Constructor 2: Using separate lambda and mu (allows lambda≠mu for future flexibility)
+    public CopyNumberBD(
+            @ParameterInfo(name = birthRateString, description = "Birth rate (gain)")
+            Value<Double> lambda,
+            @ParameterInfo(name = deathRateString, description = "Death rate (loss)")
+            Value<Double> mu,
+            @ParameterInfo(name = rootStateString, description = "Initial copy number at root", optional = true)
+            Value<Integer> rootState,
+            @ParameterInfo(name = nstateString, description = "Maximum number of copy number states (default=15)", optional = true)
+            Value<Integer> nstate
+    ) {
+        super();
+        // Use separate lambda and mu (allows lambda≠mu for future flexibility)
+        assert (lambda.value() >= 0) : "lambda must be non-negative";
+        assert (mu.value() >= 0) : "mu must be non-negative";
+        this.rootState = rootState != null ? rootState : new Value<>("rootState", 2);
+        this.nstate = nstate != null ? nstate : new Value<>("nstate", 15);
+        this.lambda = lambda;
+        this.mu = mu;
+        this.bdRate = null;
     }
 
     @Override
@@ -82,7 +127,7 @@ public int simulateCopiesOnBranchBin(double startTime, double endTime, int start
             if (u < birthRate / totalRate) {
                 // Birth event (copy gain)
                 m += 1;
-        } else {
+            } else {
                 // Death event (copy loss)
                 m -= 1;
             }
@@ -97,20 +142,41 @@ public Value<MarkovTraitEvolution<Integer>> apply() {
     }
 
     public Map<String, Value> getParams() {
-        return Map.of(
-                birthRateString, lambda,
-                deathRateString, mu,
-                rootStateString, rootState
-        );
+        Map<String, Value> params = new HashMap<>();
+
+        if (bdRate != null) {
+            // User provided bdRate
+            params.put(bdRateString, bdRate);
+        } else {
+            // User provided lambda and mu separately
+            params.put(birthRateString, lambda);
+            params.put(deathRateString, mu);
+        }
+        params.put(rootStateString, rootState);
+        params.put(nstateString, nstate);
+        return params;
     }
 
     public void setParam(String paramName, Value value) {
-        if (paramName.equals(birthRateString)) {
+        if (paramName.equals(bdRateString)) {
+            bdRate = value;
+            // The simulator always uses lambda/mu internally
+            lambda = new Value<>("lambda", bdRate.value());
+            mu = new Value<>("mu", bdRate.value());
+        } else if (paramName.equals(birthRateString)) {
+            if (bdRate != null) {
+                throw new RuntimeException("Cannot set lambda when using bdRate mode.");
+            }
             lambda = value;
         } else if (paramName.equals(deathRateString)) {
+            if (bdRate != null) {
+                throw new RuntimeException("Cannot set mu when using bdRate mode.");
+            }
             mu = value;
         } else if (paramName.equals(rootStateString)) {
             rootState = value;
+        } else if (paramName.equals(nstateString)) {
+            nstate = value;
         } else {
             throw new RuntimeException("Unrecognised parameter name: " + paramName);
         }
@@ -128,4 +194,10 @@ public Value<Integer> getRootState() {
         return rootState;
     }
 
-}
+    public Value<Double> getBdRate() {
+        return bdRate;
+    }
+    public Value<Integer> getNstate() {
+        return nstate;
+    }
+}

phylonco-lphy/src/main/java/phylonco/lphy/evolution/copynumbermodel/IntegerCharacterMatrix.java

@@ -5,6 +5,14 @@
 import lphy.core.model.annotation.GeneratorCategory;
 import lphy.core.model.annotation.MethodInfo;
 
+/**
+ * A character matrix for integer-valued traits across taxa.
+ *
+ * <p>Stores discrete integer states for multiple taxa across multiple sites or bins.</p>
+ *
+ * <p>This is analogous to {@link lphy.base.evolution.alignment.Alignment} for
+ * molecular sequences, but designed for integer-valued data.</p>
+ */
 public class IntegerCharacterMatrix implements TaxaCharacterMatrix<Integer> {
 
     private int ntaxa;
@@ -23,8 +31,7 @@ public IntegerCharacterMatrix(Taxa t, int nBins) {
     }
 
     @MethodInfo(description = "the taxa of the copy number matrix.", narrativeName = "list of taxa",
-            category = GeneratorCategory.TAXA_ALIGNMENT,
-            examples = {"copynumber.lphy"}) //needs to check
+            category = GeneratorCategory.TAXA_ALIGNMENT)
     public Taxa taxa() {
         return getTaxa();
     }
@@ -73,18 +80,18 @@ public int getDimension() {
         return 2;
     }
 
-@Override
-public String toString() {
-    StringBuilder res = new StringBuilder();
-    for (int i = 0; i < ntaxa; i++) {
-        for (int j = 0; j < nBins; j++) {
-            res.append(data[i][j]);
-            // Add separator between values
-            if (j < nBins-1) res.append(" ");
+    @Override
+    public String toString() {
+        StringBuilder res = new StringBuilder();
+        for (int i = 0; i < ntaxa; i++) {
+            for (int j = 0; j < nBins; j++) {
+                res.append(data[i][j]);
+                // Add separator between values
+                if (j < nBins - 1) res.append(" ");
+            }
+            // Add newline between taxa
+            res.append("\n");
         }
-        // Add newline between taxa
-        res.append("\n");
+        return res.toString();
     }
-    return res.toString();
-}
 }