diff --git a/lphy-base/src/main/java/lphy/base/evolution/coalescent/populationmodel/ExponentialPopulation.java b/lphy-base/src/main/java/lphy/base/evolution/coalescent/populationmodel/ExponentialPopulation.java index ef8df327..b7f20894 100644 --- a/lphy-base/src/main/java/lphy/base/evolution/coalescent/populationmodel/ExponentialPopulation.java +++ b/lphy-base/src/main/java/lphy/base/evolution/coalescent/populationmodel/ExponentialPopulation.java @@ -1,27 +1,20 @@ package lphy.base.evolution.coalescent.populationmodel; import lphy.base.evolution.coalescent.PopulationFunction; +import org.apache.commons.math3.analysis.UnivariateFunction; import org.apache.commons.math3.analysis.integration.IterativeLegendreGaussIntegrator; +import org.apache.commons.math3.analysis.solvers.BrentSolver; +import org.apache.commons.math3.analysis.solvers.UnivariateSolver; -import java.io.FileWriter; -import java.io.IOException; -import java.io.PrintWriter; import java.util.Locale; public class ExponentialPopulation implements PopulationFunction { - private double GrowthRate; + private double growthRate; private double N0; + private double NA; // Ancestral population size + private boolean useAncestralPopulation; // Flag to indicate whether to use NA - /** - * Initializes an IterativeLegendreGaussIntegrator with predefined settings for numerical integration. - * This setup is optimized for accuracy and efficiency in logistic population model computations. - * - * @return Configured IterativeLegendreGaussIntegrator with: - * - 5 Legendre-Gauss points for quadrature precision. - * - Relative accuracy of 1.0e-12 and absolute accuracy of 1.0e-8. - * - A minimum of 2 iterations and a maximum of 10,000 iterations. - */ private IterativeLegendreGaussIntegrator createIntegrator() { int numberOfPoints = 5; // Legendre-Gauss points double relativeAccuracy = 1.0e-12; // relative precision @@ -32,139 +25,233 @@ private IterativeLegendreGaussIntegrator createIntegrator() { } - public ExponentialPopulation(double GrowthRate, double N0) { - this.GrowthRate = GrowthRate; + /** + * Constructor without ancestral population. + * + * @param growthRate Growth rate of the population. + * @param N0 Initial population size at time t=0. + */ + public ExponentialPopulation(double growthRate, double N0) { + if (N0 <= 0) { + throw new IllegalArgumentException("Initial population size N0 must be positive."); + } + this.growthRate = growthRate; this.N0 = N0; + this.useAncestralPopulation = false; + this.NA = 0.0; // Default value when not using NA } /** - * Calculate the growth value at a given time t.。 + * Constructor with ancestral population. + * + * @param growthRate Growth rate of the population. + * @param N0 Initial population size at time t=0. + * @param ancestralPopulationSize Ancestral population size NA. + */ + public ExponentialPopulation(double growthRate, double N0, double ancestralPopulationSize) { + if (N0 <= 0) { + throw new IllegalArgumentException("Initial population size N0 must be positive."); + } + if (ancestralPopulationSize <= 0) { + throw new IllegalArgumentException("Ancestral population size NA must be positive."); + } + if (ancestralPopulationSize > N0) { + throw new IllegalArgumentException("Ancestral population size NA cannot exceed initial population size N0."); + } + this.growthRate = growthRate; + this.N0 = N0; + this.useAncestralPopulation = true; + this.NA = ancestralPopulationSize; + } + + /** + * Calculates the population size at a given time t. + * + * @param t Time at which to calculate the population size. + * @return Population size N(t) at time t. */ @Override public double getTheta(double t) { - - double r = GrowthRate; + double r = growthRate; if (r == 0) { - return N0; + return useAncestralPopulation ? NA : N0; } else { - return N0 * Math.exp(-t * r); + if (useAncestralPopulation) { + return (N0 - NA) * Math.exp(-r * t) + NA; + } else { + return N0 * Math.exp(-r * t); + } } } /** - * Calculate the cumulative intensity over time period from 0 to t. + * Calculates the cumulative intensity from time 0 to t. + * + * @param t Time up to which to calculate the cumulative intensity. + * @return Cumulative intensity up to time t. */ - @Override public double getIntensity(double t) { + if (t == 0) return 0.0; - double r = GrowthRate; - if (r == 0.0) { - return t / N0; + if (!useAncestralPopulation) { + // Without NA, use analytical solution + double r = growthRate; + if (r == 0.0) { + return t / N0; + } + return (Math.exp(r * t) - 1.0) / (r * N0); } else { - return (Math.exp(t * r) - 1.0) / N0 / r; + // Numerical integration when using Ancestral Population (NA) + UnivariateFunction function = time -> 1.0 / getTheta(time); + IterativeLegendreGaussIntegrator integrator = createIntegrator(); + + return integrator.integrate(Integer.MAX_VALUE, function, 0, t); } } -// -// if (t == 0) return 0; -// -// -// UnivariateFunction function = time -> 1 / getTheta(time); -// IterativeLegendreGaussIntegrator integrator = createIntegrator(); -// return integrator.integrate(Integer.MAX_VALUE, function, 0, t); -// } - - -// @Override -// public double getInverseIntensity(double x) { -//// If min is set to 0 in BrentSolver, an error will be reported: -//// org.apache.commons.math3.exception.NumberIsTooLargeException: endpoints do not specify an interval: [0, 0] -//// so I change min equals to 0.1 -// UnivariateFunction function = time -> getIntensity(time) - x; -// UnivariateSolver solver = new BrentSolver(); -// // The range [0, 100] might need to be adjusted depending on the growth model and expected time range. -//// return solver.solve(100, function, 0, 100); -// return solver.solve(100, function, 0.1, 50); -// } + /** + * Calculates the time t corresponding to a given cumulative intensity x. + * + * @param x Cumulative intensity. + * @return Time t such that Intensity(t) = x. + */ @Override public double getInverseIntensity(double x) { - double r = GrowthRate; - if (r == 0.0) { - return N0 * x; + if (x == 0.0) return 0.0; + + if (!useAncestralPopulation) { + // Without NA, use analytical solution + double r = growthRate; + if (r == 0.0) { + return N0 > 0 ? x * N0 : 0.0; + } + return Math.log(1.0 + x * r * N0) / r; } else { - return Math.log(1.0 + N0 * x * r) / r; + // With NA, use numerical solver + UnivariateFunction function = time -> getIntensity(time) - x; + UnivariateSolver solver = new BrentSolver(); + double tMin = 0.0; + + // Estimate initial tMax based on the given intensity x + double tMax = estimateInitialTMax(x); + + // Ensure tMax does not exceed half of Double.MAX_VALUE + tMax = Math.min(tMax, Double.MAX_VALUE / 2); + + // Attempt to solve within [tMin, tMax] + try { + return solver.solve(100, function, tMin, tMax); + } catch (Exception e) { + // If not found, double tMax and retry + tMax *= 2.0; + while (tMax < Double.MAX_VALUE / 2) { + try { + return solver.solve(100, function, tMin, tMax); + } catch (Exception ex) { + tMax *= 2.0; + } + } + throw new RuntimeException("Failed to find a valid time for given intensity: " + x, e); + } } } -// UnivariateFunction function = time -> getIntensity(time) - x; -// UnivariateSolver solver = new BrentSolver(); -// double tMin = 0; -// double tMax = 10; -// -// double lastIntensity = getIntensity(tMin); -// System.out.println("Starting search: tMin=" + tMin + ", Intensity=" + lastIntensity); -// -// -// -// while (lastIntensity < x) { -// tMax *= 2; -// lastIntensity = getIntensity(tMax); -// System.out.println("Expanding search: tMax=" + tMax + ", Intensity=" + lastIntensity); -// -// if (tMax >= Double.MAX_VALUE / 2) { -// System.out.println("Reached maximum double value for tMax."); -// tMax = Double.MAX_VALUE / 2; -// break; -// } -// } -// -// try { -// double result = solver.solve(100, function, tMin, tMax); -// System.out.println("Found time for given intensity: " + result); -// System.out.println("Found time for given intensity: " + result + ", x=" + x); -// return result; -// } catch (Exception e) { -// System.out.println("Failed at tMax=" + tMax + " with lastIntensity=" + lastIntensity); -// throw new RuntimeException("Failed to find a valid time for given intensity: " + x, e); -// } -// } + /** + * Estimates an initial tMax value to improve solver efficiency. + * + * @param x Cumulative intensity. + * @return Estimated tMax. + */ + private double estimateInitialTMax(double x) { + // Based on cumulative intensity x and model parameters, estimate tMax + // Assume when N(t) approaches NA, Intensity(t) ≈ t / (NA * (N0 - NA)) + // Hence, estimate tMax ≈ x * NA * (N0 - NA) * 1.5 + double estimatedTMax = x * NA * (N0 - NA) * 1.5; + + // Prevent tMax from being too small + estimatedTMax = Math.max(estimatedTMax, 1.0); + + return estimatedTMax; + } + + /** + * Checks whether the model is using an ancestral population. + * + * @return True if using NA, false otherwise. + */ + public boolean isUsingAncestralPopulation() { + return useAncestralPopulation; + } + + /** + * Retrieves the ancestral population size NA. + * + * @return Ancestral population size NA. + */ + public double getAncestralPopulation() { + return NA; + } + /** + * Indicates whether the model uses an analytical solution. + * + * @return True if using analytical solution (without NA), false otherwise. + */ @Override public boolean isAnalytical() { - return false; //use numerical method here + return !useAncestralPopulation; // Analytical only when not using NA } + /** + * Provides a string representation of the ExponentialPopulation model. + * + * @return String describing the model parameters. + */ @Override public String toString() { - return "Exponential growth rate: " + GrowthRate + ", initial size: " + N0; + if (useAncestralPopulation) { + return "Exponential Population: Growth Rate = " + growthRate + + ", Initial Size = " + N0 + ", Ancestral Population = " + NA; + } else { + return "Exponential Population: Growth Rate = " + growthRate + + ", Initial Size = " + N0; + } } + /** + * Main method for testing the ExponentialPopulation class functionality. + * + * @param args Command-line arguments. + */ + public static void main(String[] args) { + double growthRate = 0.1; + double N0 = 100.0; + double ancestralPopulation = 50.0; + + // Define specific test time points + double[] testTimes = {0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0}; + // Create an instance without NA + ExponentialPopulation expPop = new ExponentialPopulation(growthRate, N0); - public static void main(String[] args) { - double GrowthRate = 0.1; - double N0 = 3; - double tStart = 0; - double tEnd = 100; - int nPoints = 100; + // Create an instance with NA + ExponentialPopulation expPopWithNA = new ExponentialPopulation(growthRate, N0, ancestralPopulation); - ExponentialPopulation exponentialPopulation = new ExponentialPopulation(GrowthRate, N0); + // Print header for test points + System.out.println("time,theta_noNA,theta_withNA,intensity_noNA,intensity_withNA"); - try (PrintWriter writer = new PrintWriter(new FileWriter("exponential_data.csv"))) { - writer.println("time,theta"); - for (int i = 0; i < nPoints; i++) { - double t = tStart + (i / (double) (nPoints - 1)) * (tEnd - tStart); - double theta = exponentialPopulation.getTheta(t); + // Iterate over each test time point + for (double t : testTimes) { + double theta_noNA = expPop.getTheta(t); + double theta_withNA = expPopWithNA.getTheta(t); + double intensity_noNA = expPop.getIntensity(t); + double intensity_withNA = expPopWithNA.getIntensity(t); - writer.printf(Locale.US, "%.4f,%.4f%n", t, theta); - } - } catch (IOException e) { - e.printStackTrace(); + // Print the results in CSV format + System.out.printf(Locale.US, "%.2f,%.4f,%.4f,%.6f,%.6f%n", + t, theta_noNA, theta_withNA, intensity_noNA, intensity_withNA); } } - } - - diff --git a/lphy-base/src/main/java/lphy/base/evolution/coalescent/populationmodel/ExponentialPopulationFunction.java b/lphy-base/src/main/java/lphy/base/evolution/coalescent/populationmodel/ExponentialPopulationFunction.java index 2bad5022..2e65fd82 100644 --- a/lphy-base/src/main/java/lphy/base/evolution/coalescent/populationmodel/ExponentialPopulationFunction.java +++ b/lphy-base/src/main/java/lphy/base/evolution/coalescent/populationmodel/ExponentialPopulationFunction.java @@ -11,12 +11,37 @@ public class ExponentialPopulationFunction extends DeterministicFunction GrowthRate, - @ParameterInfo(name = N0_PARAM_NAME, description = "The initial population size.") Value N0) { + public static final String NA_PARAM_NAME = "NA"; + + /** + * Constructor for ExponentialPopulationFunction without Ancestral Population (NA). + * + * @param GrowthRate Growth rate of the population. + * @param N0 Initial population size at time t=0. + */ + public ExponentialPopulationFunction( + @ParameterInfo(name = GROWTH_RATE_PARAM_NAME, description = "The growth rate of the population.") Value GrowthRate, + @ParameterInfo(name = N0_PARAM_NAME, description = "The initial population size.") Value N0) { setParam(GROWTH_RATE_PARAM_NAME, GrowthRate); setParam(N0_PARAM_NAME, N0); } + /** + * Constructor for ExponentialPopulationFunction with Ancestral Population (NA). + * + * @param GrowthRate Growth rate of the population. + * @param N0 Initial population size at time t=0. + * @param NA Ancestral population size. + */ + public ExponentialPopulationFunction( + @ParameterInfo(name = GROWTH_RATE_PARAM_NAME, description = "The growth rate of the population.") Value GrowthRate, + @ParameterInfo(name = N0_PARAM_NAME, description = "The initial population size.") Value N0, + @ParameterInfo(name = NA_PARAM_NAME, description = "The ancestral population size.") Value NA) { + setParam(GROWTH_RATE_PARAM_NAME, GrowthRate); + setParam(N0_PARAM_NAME, N0); + setParam(NA_PARAM_NAME, NA); + } + @GeneratorInfo(name="exponentialPopFunc", narrativeName = "Exponential growth function", category = GeneratorCategory.COAL_TREE, examples = {" expCoal.lphy, expCoalJC.lphy" }, description = "Models population growth using an exponential growth function.") @@ -25,8 +50,18 @@ public Value apply() { double GrowthRate = ((Number) getParams().get(GROWTH_RATE_PARAM_NAME).value()).doubleValue(); double N0 = ((Number) getParams().get(N0_PARAM_NAME).value()).doubleValue(); + // Check if NA parameter is provided + Value naValue = getParams().get(NA_PARAM_NAME); + PopulationFunction exponentialPopulation; - PopulationFunction exponentialPopulation = new ExponentialPopulation(GrowthRate, N0); + if (naValue != null && naValue.value() != null && naValue.value() > 0.0) { + double NA = naValue.value().doubleValue(); + // Create ExponentialPopulation with NA + exponentialPopulation = new ExponentialPopulation(GrowthRate, N0, NA); + } else { + // Create ExponentialPopulation without NA + exponentialPopulation = new ExponentialPopulation(GrowthRate, N0); + } return new Value<>(exponentialPopulation, this); } @@ -40,5 +75,9 @@ public Value getN0() { return (Value) getParams().get(N0_PARAM_NAME); } + public Value getNA() { + return (Value) getParams().get(NA_PARAM_NAME); + } + } diff --git a/lphy-base/src/test/java/lphy/base/evolution/coalescent/populationmodel/ExponentialPopulationTest.java b/lphy-base/src/test/java/lphy/base/evolution/coalescent/populationmodel/ExponentialPopulationTest.java index 6e174a02..0cd93c95 100644 --- a/lphy-base/src/test/java/lphy/base/evolution/coalescent/populationmodel/ExponentialPopulationTest.java +++ b/lphy-base/src/test/java/lphy/base/evolution/coalescent/populationmodel/ExponentialPopulationTest.java @@ -1,6 +1,7 @@ package lphy.base.evolution.coalescent.populationmodel; - +import org.junit.jupiter.api.BeforeEach; +import org.junit.jupiter.api.DisplayName; import org.junit.jupiter.api.Test; import static org.junit.jupiter.api.Assertions.assertEquals; @@ -8,57 +9,136 @@ public class ExponentialPopulationTest { private static final double DELTA = 1e-6; + private double growthRate; + private double N0; + private double NA; + private ExponentialPopulation expPop; + private ExponentialPopulation expPopWithNA; /** - * Test the getTheta method of the ExponentialPopulation class. - * This test verifies whether the getTheta method can correctly calculate theta value based on the exponential growth model. + * Setup method to initialize common test parameters and instances before each test. */ - @Test - public void testGetTheta() { - // Set the parameters for the exponential growth model - double growthRate = 0.1; // Growth rate (r) - double N0 = 100; // Initial population size + @BeforeEach + public void setUp() { + // Initialize parameters + growthRate = 0.1; // Growth rate (r) + N0 = 100.0; // Initial population size + NA = 50.0; // Ancestral population size + + // Create instances + expPop = new ExponentialPopulation(growthRate, N0); + expPopWithNA = new ExponentialPopulation(growthRate, N0, NA); + } - ExponentialPopulation exponentialPopulation = new ExponentialPopulation(growthRate, N0); + /** + * Test the getTheta method without using Ancestral Population (NA). + * This test verifies whether the getTheta method correctly calculates theta based on the exponential growth model without NA. + */ + @Test + @DisplayName("Test getTheta without Ancestral Population (NA)") + public void testGetThetaWithoutNA() { + double t = 10.0; // Time at which to calculate theta - // Test getTheta method for a given time t - double t = 80000; // Time at which to calculate theta - // Calculate expected theta using the exponential growth formula: N0 * exp(growthRate * t) - double expectedTheta = N0 * Math.exp(growthRate * -t); - double actualTheta = exponentialPopulation.getTheta(t); + // Calculate expected theta using the exponential growth formula: N0 * exp(-r * t) + double expectedTheta = N0 * Math.exp(-growthRate * t); + double actualTheta = expPop.getTheta(t); - // Assert that the expected value is equal to the actual value, allowing a certain error range - assertEquals(expectedTheta, actualTheta, DELTA, "The theta value at time t should match the expected value calculated using the exponential growth formula."); + // Assert that the expected value is equal to the actual value within the specified delta + assertEquals(expectedTheta, actualTheta, DELTA, "Theta without NA should match the expected exponential decay value."); } - - + /** + * Test the getTheta method with using Ancestral Population (NA). + * This test verifies whether the getTheta method correctly calculates theta based on the exponential growth model with NA. + */ @Test - public void testIntensityAndInverseIntensity() { - // Parameters for the exponential growth model - double growthRate = 0.1; // Growth rate - double N0 = 100; // Initial population size + @DisplayName("Test getTheta with Ancestral Population (NA)") + public void testGetThetaWithNA() { + double t = 10.0; // Time at which to calculate theta - ExponentialPopulation model = new ExponentialPopulation(growthRate, N0); + // Calculate expected theta using the exponential growth formula with NA: + // N(t) = (N0 - NA) * exp(-r * t) + NA + double expectedTheta = (N0 - NA) * Math.exp(-growthRate * t) + NA; + double actualTheta = expPopWithNA.getTheta(t); + + // Assert that the expected value is equal to the actual value within the specified delta + assertEquals(expectedTheta, actualTheta, DELTA, "Theta with NA should match the expected value considering ancestral population."); + } - double t = 2; // Time point to calculate intensity and then retrieve it through inverse intensity + /** + * Test the getIntensity and getInverseIntensity methods without using Ancestral Population (NA). + * This test verifies the consistency between intensity and its inverse without NA. + */ + @Test + @DisplayName("Test Intensity and InverseIntensity without Ancestral Population (NA)") + public void testIntensityAndInverseIntensityWithoutNA() { + double t = 2.0; // Time point to calculate intensity and then retrieve it through inverse intensity - // Calculate cumulative intensity at a given point in time - double intensityAtTimeT = model.getIntensity(t); + // Calculate cumulative intensity at time t + double intensityAtTimeT = expPop.getIntensity(t); // Use the accumulated intensity value to calculate its corresponding time - double inverseIntensityResult = model.getInverseIntensity(intensityAtTimeT); + double inverseIntensityResult = expPop.getInverseIntensity(intensityAtTimeT); // Output for debugging + System.out.println("Without NA:"); System.out.println("t = " + t); System.out.println("intensity = " + intensityAtTimeT); System.out.println("inverse intensity = " + inverseIntensityResult); - // Verify whether the inverse intensity calculation result is close to the original time point - assertEquals(t, inverseIntensityResult, DELTA, "Inverse intensity calculation should return the original time point within an acceptable error margin."); + // Assert that the inverse intensity result is equal to the original time within the specified delta + assertEquals(t, inverseIntensityResult, DELTA, "Inverse intensity should return the original time without NA."); } + /** + * Test the getIntensity and getInverseIntensity methods with using Ancestral Population (NA). + * This test verifies the consistency between intensity and its inverse with NA. + */ + @Test + @DisplayName("Test Intensity and InverseIntensity with Ancestral Population (NA)") + public void testIntensityAndInverseIntensityWithNA() { + double t = 122.000007; // Time point to calculate intensity and then retrieve it through inverse intensity + + // Calculate cumulative intensity at time t with NA + double intensityAtTimeT = expPopWithNA.getIntensity(t); -} + // Use the accumulated intensity value to calculate its corresponding time + double inverseIntensityResult = expPopWithNA.getInverseIntensity(intensityAtTimeT); + + // Output for debugging + System.out.println("With NA:"); + System.out.println("t = " + t); + System.out.println("intensity = " + intensityAtTimeT); + System.out.println("inverse intensity = " + inverseIntensityResult); + // Assert that the inverse intensity result is equal to the original time within the specified delta + assertEquals(t, inverseIntensityResult, DELTA, "Inverse intensity should return the original time with NA."); + } + /** + * Test the getInverseIntensity method with a range of intensity values with NA. + * This test ensures that inverseIntensity correctly retrieves the original time for various intensities. + */ + @Test + @DisplayName("Test InverseIntensity with multiple intensity values with NA") + public void testInverseIntensityMultipleWithNA() { + double[] testTimes = {0.0, 5.0, 10.0, 20.0, 50.0, 100.0}; // Various time points + + for (double t : testTimes) { + // Calculate cumulative intensity at time t with NA + double intensity = expPopWithNA.getIntensity(t); + + // Use the accumulated intensity value to calculate its corresponding time + double inverseIntensityResult = expPopWithNA.getInverseIntensity(intensity); + + // Output for debugging + System.out.println("With NA:"); + System.out.println("t = " + t); + System.out.println("intensity = " + intensity); + System.out.println("inverse intensity = " + inverseIntensityResult); + + // Assert that the inverse intensity result is equal to the original time within the specified delta + assertEquals(t, inverseIntensityResult, DELTA, "Inverse intensity should return the original time for multiple points with NA."); + } + } +}