Skip to content

Commit

Permalink
change gompertz f0 version
Browse files Browse the repository at this point in the history
  • Loading branch information
@yxu927
yxu927 committed May 6, 2024
1 parent cea7f03 commit 564259f
Show file tree
Hide file tree
Showing 4 changed files with 209 additions and 224 deletions.
21 changes: 12 additions & 9 deletions examples/coalescent/gompertzCoalescent.lphy
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,17 @@

data {
L = 500;
L = 200;
}

model {

N0 ~ LogNormal(meanlog=3,sdlog=0.5);
b ~ LogNormal(meanlog=-0.9, sdlog=0.05);
NInfinity ~ LogNormal(meanlog=28, sdlog=0.5);
gompertzPopFunc = gompertz(N0=N0, b=b, NInfinity=NInfinity);
f0 ~ Beta(alpha=20, beta=7);

tree ~ CoalescentPopFunc(popFunc=gompertzPopFunc, n=16);
// rootAge = tree.rootAge();
D ~ PhyloCTMC(tree=tree, L=L, Q=jukesCantor(), mu=1e-9);
}
b ~ LogNormal(meanlog=-0.3, sdlog=0.05);
NInfinity ~ LogNormal(meanlog=35, sdlog=0.5);
gompertzPopFunc = gompertz(f0=f0, b=b, NInfinity=NInfinity);

tree ~ CoalescentPopFunc(popFunc=gompertzPopFunc, n=16);
rootAge = tree.rootAge();
D ~ PhyloCTMC(tree=tree, L=L, Q=jukesCantor(), mu=1e-2);
}
161 changes: 76 additions & 85 deletions ...base/src/main/java/lphy/base/evolution/coalescent/populationmodel/GompertzPopulation.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,10 +22,6 @@ public class GompertzPopulation implements PopulationFunction {
public static final String T50ParamName = "t50";
public static final String F0ParamName = "f0";

public static final String N0ParamName = "N0";



private double N0; // Initial population size
private double b; // Initial growth rate of tumor growth
private double NInfinity; // Carrying capacity
Expand All @@ -41,9 +37,9 @@ public static double computeT50(double NInfinity, double N0, double b) {
}
double ratio = NInfinity / N0;
double proportion = 0.5;
// if (proportion <= 0 || proportion >= 1) {
// throw new IllegalArgumentException("Proportion must be between 0 and 1.");
// }
// if (proportion <= 0 || proportion >= 1) {
// throw new IllegalArgumentException("Proportion must be between 0 and 1.");
// }
double t50 = Math.log(1 - Math.log(proportion) / Math.log(ratio)) / b;
return t50;
}
Expand Down Expand Up @@ -89,38 +85,37 @@ private IterativeLegendreGaussIntegrator createIntegrator() {
* @param b
* @param NInfinity
*/
// public GompertzPopulation(double t50, double b, double NInfinity) { //(this is for t50 )
// this.b = b;
// this.t50 = t50;
// this.NInfinity = NInfinity;
// // Calculate N0 based on t50, b, and NInfinity
// // N(t50) = NInfinity / 2
// // t50 is a time location given by the user, t50 < 0 means it is in the early exponential phase
// this.N0 = NInfinity * Math.pow(2, -Math.exp(-b * this.t50));
// }

// public GompertzPopulation(double f0, double b, double NInfinity) { //(this is for f0 )
// this.f0 = f0;
// this.b = b;
// this.NInfinity = NInfinity;
// this.N0 = NInfinity * this .f0;
// }

public GompertzPopulation(double N0, double b, double NInfinity) {
// this.f0 = f0;
// public GompertzPopulation(double t50, double b, double NInfinity) { //(this is for t50 )
// this.b = b;
// this.t50 = t50;
// this.NInfinity = NInfinity;
// // Calculate N0 based on t50, b, and NInfinity
// // N(t50) = NInfinity / 2
// // t50 is a time location given by the user, t50 < 0 means it is in the early exponential phase
// this.N0 = NInfinity * Math.pow(2, -Math.exp(-b * this.t50));
// }

public GompertzPopulation(double f0, double b, double NInfinity) { //(this is for f0 )
this.f0 = f0;
this.b = b;
this.NInfinity = NInfinity;
this.N0 = N0;
//this.N0 = NInfinity * this.f0;
this.N0 = NInfinity * this .f0;
}



// @Override //(this is for t50 )
// public double getTheta(double t) {
// // the sign of b * t is such that t = 0 is present time and t > 0 is time in the past
// return N0 * Math.exp(Math.log(NInfinity / N0) * (1 - Math.exp(b * t)));
// }
/**
* Implement the Gompertz function to calculate theta at time t
* Assuming theta is proportional to population size for simplicity
*
* @param t time, where t > 0 is time in the past
* @return N0 * Math.exp(Math.log(NInfinity / N0) * (1 - Math.exp(b * t)))
*/
// @Override //(this is for t50 )
// public double getTheta(double t) {
// // the sign of b * t is such that t = 0 is present time and t > 0 is time in the past
// return N0 * Math.exp(Math.log(NInfinity / N0) * (1 - Math.exp(b * t)));
// }

@Override //(this is for f0 )
public double getTheta(double t) {
Expand All @@ -137,13 +132,13 @@ public double getIntensity(double t) {
UnivariateFunction function = time -> 1 / getTheta(time);

// Use the separate method to create the integrator
// return legrandeIntegrator(function, t);
// return legrandeIntegrator(function, t);
IterativeLegendreGaussIntegrator integrator = createIntegrator();
return integrator.integrate(Integer.MAX_VALUE, function, 0, t);
}

// return rombergIntegrator(function, t);
//}
// return rombergIntegrator(function, t);
//}

private double legrandeIntegrator(UnivariateFunction function, double t) {
IterativeLegendreGaussIntegrator integrator = createIntegrator();
Expand All @@ -153,68 +148,68 @@ private double legrandeIntegrator(UnivariateFunction function, double t) {
private double rombergIntegrator(UnivariateFunction function, double t) {

int maxBound = Integer.MAX_VALUE;
// int maxBound = 100000;
// int maxBound = 100000;

RombergIntegrator integrator = new RombergIntegrator();
double absoluteAccuracy = 1e-6;
integrator = new RombergIntegrator(RombergIntegrator.DEFAULT_RELATIVE_ACCURACY,
absoluteAccuracy, RombergIntegrator.DEFAULT_MIN_ITERATIONS_COUNT,
RombergIntegrator.ROMBERG_MAX_ITERATIONS_COUNT);

// try {
// System.out.println("Absolute accuracy = " + integrator.getAbsoluteAccuracy());
// System.out.println("Relative accuracy = " + integrator.getRelativeAccuracy());
// try {
// System.out.println("Absolute accuracy = " + integrator.getAbsoluteAccuracy());
// System.out.println("Relative accuracy = " + integrator.getRelativeAccuracy());
return integrator.integrate(maxBound, function, 0, t);
// }
// }
}






// @Override
// public double getIntensity(double t) {
//
// if (t == 0) return 0;
//
// if (getTheta(t) < NInfinity/resolution_magic_number) {
// throw new RuntimeException("Theta too small to calculate intensity!");
// }
//
// UnivariateFunction function = time -> 1.0 / getTheta(time);
// UnivariateIntegrator integrator = new TrapezoidIntegrator();
// // The number 10000 here represents a very high number of iterations for accuracy.
// return integrator.integrate(100000, function, 0, t);
// }
//
// @Override
// public double getIntensity(double t) {
//
// if (t == 0) return 0;
//
// if (getTheta(t) < NInfinity/resolution_magic_number) {
// throw new RuntimeException("Theta too small to calculate intensity!");
// }
//
// UnivariateFunction function = time -> 1.0 / getTheta(time);
// UnivariateIntegrator integrator = new TrapezoidIntegrator();
// // The number 10000 here represents a very high number of iterations for accuracy.
// return integrator.integrate(100000, function, 0, t);
// }
//


// passes unit test without magic number
// Error when running LPhy script that signs at interval endpoints are not different
// @Override
// public double getInverseIntensity(double x) {
////
////
//// UnivariateFunction thetaFunction = time -> getTheta(time) - NInfinity / resolution_magic_number;
//// UnivariateSolver thetaSolver = new BrentSolver();
//// double startValue = thetaFunction.value(0);
//// double endValue = thetaFunction.value(1000);
////
//// System.out.println("Function value at start of the interval (0): " + startValue);
//// System.out.println("Function value at end of the interval (1000): " + endValue);
//// double maxTime = thetaSolver.solve(100, thetaFunction, 1e-6, t50 * 10);
//// System.out.println("maxTime = " + maxTime);
////
//// 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.001, maxTime);
////
////}
//
//
// Error when running LPhy script that signs at interval endpoints are not different
// @Override
// public double getInverseIntensity(double x) {
////
////
//// UnivariateFunction thetaFunction = time -> getTheta(time) - NInfinity / resolution_magic_number;
//// UnivariateSolver thetaSolver = new BrentSolver();
//// double startValue = thetaFunction.value(0);
//// double endValue = thetaFunction.value(1000);
////
//// System.out.println("Function value at start of the interval (0): " + startValue);
//// System.out.println("Function value at end of the interval (1000): " + endValue);
//// double maxTime = thetaSolver.solve(100, thetaFunction, 1e-6, t50 * 10);
//// System.out.println("maxTime = " + maxTime);
////
//// 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.001, maxTime);
////
////}
//
//



Expand Down Expand Up @@ -295,8 +290,4 @@ public static void main(String[] args) {



}




}
32 changes: 10 additions & 22 deletions .../main/java/lphy/base/evolution/coalescent/populationmodel/GompertzPopulationFunction.java
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@

package lphy.base.evolution.coalescent.populationmodel;

import lphy.base.evolution.coalescent.PopulationFunction;
Expand All @@ -11,20 +12,15 @@

public class GompertzPopulationFunction extends DeterministicFunction<PopulationFunction> {
public GompertzPopulationFunction(//@ParameterInfo(name = T50ParamName, description = "Time when population is half of carrying capacity.") Value<Double> t50,
@ParameterInfo(name = N0ParamName, description = "Time when population is half of carrying capacity.") Value<Double> N0,
// @ParameterInfo(name = F0ParamName, description = "Time when population is half of carrying capacity.") Value<Double> f0,
@ParameterInfo(name = BParamName, description = "Initial growth rate of tumor growth.") Value<Double> b,
@ParameterInfo(name = NINFINITYParamName, description = "Limiting population size (carrying capacity).") Value<Double> NInfinity) {
@ParameterInfo(name = F0ParamName, description = "Time when population is half of carrying capacity.") Value<Double> f0,
@ParameterInfo(name = BParamName, description = "Initial growth rate of tumor growth.") Value<Number> b,
@ParameterInfo(name = NINFINITYParamName, description = "Limiting population size (carrying capacity).") Value<Number> NInfinity) {
//setParam(T50ParamName, t50);
// setParam(F0ParamName, f0);
setParam(F0ParamName, f0);
setParam(BParamName, b);
setParam(NINFINITYParamName, NInfinity);

setParam(N0ParamName, N0);

}


@GeneratorInfo(name="gompertz", narrativeName = "Gompertz growth function",
category = GeneratorCategory.COAL_TREE, examples = {" .lphy" },
description = "Models population growth using the Gompertz growth function.")
Expand All @@ -33,27 +29,19 @@ public GompertzPopulationFunction(//@ParameterInfo(name = T50ParamName, descript
public Value<PopulationFunction> apply() {

//double t50 = ((Number) getParams().get(T50ParamName).value()).doubleValue();
//double f0 = ((Number) getParams().get(F0ParamName).value()).doubleValue();
double f0 = ((Number) getParams().get(F0ParamName).value()).doubleValue();
double b = ((Number) getParams().get(BParamName).value()).doubleValue();
double NInfinity = ((Number) getParams().get(NINFINITYParamName).value()).doubleValue();
double N0 = ((Number) getParams().get(N0ParamName).value()).doubleValue();

//PopulationFunction gompertzPopulation = new GompertzPopulation(t50, b, NInfinity);

//PopulationFunction gompertzPopulation = new GompertzPopulation(f0, b, NInfinity);
PopulationFunction gompertzPopulation = new GompertzPopulation(N0, b, NInfinity);
PopulationFunction gompertzPopulation = new GompertzPopulation(f0, b, NInfinity);

return new Value<>(gompertzPopulation, this);
return new Value<>( gompertzPopulation, this);
}



// public Value<Double> getF0() {
// return getParams().get(F0ParamName);
// }

public Value<Number> getN0() {
return getParams().get(N0ParamName);
public Value<Double> getF0() {
return getParams().get(F0ParamName);
}

public Value<Number> getB() {
Expand Down
Loading

0 comments on commit 564259f

Please sign in to comment.