Update population function #458

LinguaPhylo · Apr 18, 2024 · 54e72ca · 54e72ca
1 parent 28fa45d
commit 54e72ca
Show file tree

Hide file tree

Showing 5 changed files with 41 additions and 3 deletions.
diff --git a/examples/coalescent/gompertzCoal.lphy → examples/coalescent/gompertzCoalescent.lphy b/examples/coalescent/gompertzCoal.lphy → examples/coalescent/gompertzCoalescent.lphy
diff --git a/lphy-base/src/main/java/lphy/base/evolution/coalescent/PopulationFunctionCoalescent.java b/lphy-base/src/main/java/lphy/base/evolution/coalescent/PopulationFunctionCoalescent.java
@@ -19,9 +19,9 @@
 import java.util.TreeMap;
 
 public class PopulationFunctionCoalescent extends TaxaConditionedTreeGenerator {
-    private Value<PopulationFunction> popFunc;
+    public Value<PopulationFunction> popFunc;
 
-    private final String popFuncParamName = "popFunc";
+    public final String popFuncParamName = "popFunc";
 
     /**
      * Constructs a coalescent model with specified population function, number of taxa, taxa object, and leaf node ages.

diff --git a/...in/java/lphy/base/evolution/coalescent/populationmodel/ExponentialPopulationFunction.java b/...in/java/lphy/base/evolution/coalescent/populationmodel/ExponentialPopulationFunction.java
@@ -32,5 +32,14 @@ public Value<PopulationFunction> apply() {
         return new Value<>(exponentialPopulation, this);
     }
 
+    public Value<Double> getGrowthRate() {
+        return (Value<Double>) getParams().get(GROWTH_RATE_PARAM_NAME);
+    }
+
+
+    public Value<Double> getN0() {
+        return (Value<Double>) getParams().get(N0_PARAM_NAME);
+    }
+
 
 }
diff --git a/...base/src/main/java/lphy/base/evolution/coalescent/populationmodel/GompertzPopulation.java b/...base/src/main/java/lphy/base/evolution/coalescent/populationmodel/GompertzPopulation.java
@@ -4,6 +4,10 @@
 import org.apache.commons.math3.analysis.UnivariateFunction;
 import org.apache.commons.math3.analysis.integration.IterativeLegendreGaussIntegrator;
 import org.apache.commons.math3.analysis.integration.RombergIntegrator;
+import org.apache.commons.math3.analysis.solvers.BrentSolver;
+import org.apache.commons.math3.analysis.solvers.UnivariateSolver;
+import org.apache.commons.math3.exception.NoBracketingException;
+import org.apache.commons.math3.exception.TooManyEvaluationsException;
 
 import java.io.FileWriter;
 import java.io.IOException;
@@ -233,11 +237,22 @@ public double getInverseIntensity(double x) {
         while (intensity < targetIntensity) {
             time += deltaTime;
             intensity = getIntensity(time);
+        }
 
+        double lowerBound = Math.max(0, time - deltaTime);
+        double upperBound = time;
+        UnivariateFunction function = t -> getIntensity(t) - x;
+        UnivariateSolver solver = new BrentSolver(1e-9, 1e-9);
+        try {
+            return solver.solve(100, function, lowerBound, upperBound);
+        } catch (NoBracketingException | TooManyEvaluationsException e) {
+            System.err.println("Solver failed: " + e.getMessage());
+            return Double.NaN;
         }
 
+
         // return time;
-        return Math.max(time, 0);
+        //return Math.max(time, 0);
     }
 
 

diff --git a/.../main/java/lphy/base/evolution/coalescent/populationmodel/LogisticPopulationFunction.java b/.../main/java/lphy/base/evolution/coalescent/populationmodel/LogisticPopulationFunction.java
@@ -35,4 +35,18 @@ public Value<PopulationFunction> apply() {
 
         return new Value<>(logisticPopulation, this);
     }
+    public Value<Double> getT50() {
+        return (Value<Double>) getParams().get(T50_PARAM_NAME);
+    }
+
+
+    public Value<Double> getNCarryingCapacity() {
+        return (Value<Double>) getParams().get(N_CARRYING_CAPACITY_PARAM_NAME);
+    }
+
+
+    public Value<Double> getB() {
+        return (Value<Double>) getParams().get(B_PARAM_NAME);
+    }
 }
+