/*
* LogisticGrowth.java
*
* Copyright (c) 2002-2015 Alexei Drummond, Andrew Rambaut and Marc Suchard
*
* This file is part of BEAST.
* See the NOTICE file distributed with this work for additional
* information regarding copyright ownership and licensing.
*
* BEAST is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* BEAST is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with BEAST; if not, write to the
* Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
* Boston, MA 02110-1301 USA
*/
package dr.evolution.coalescent;
/**
* This class models logistic growth.
*
* @author Alexei Drummond
* @author Andrew Rambaut
* @version $Id: LogisticGrowth.java,v 1.15 2005/05/24 20:25:56 rambaut Exp $
*/
public class LogisticGrowth extends ExponentialGrowth {
/**
* Construct demographic model with default settings
*/
public LogisticGrowth(Type units) {
super(units);
}
public void setShape(double value) {
c = value;
}
public double getShape() {
return c;
}
double lowLimit = 0; // 1e-6;
/**
* An alternative parameterization of this model. This
* function sets the time at which there is a 0.5 proportion
* of N0.
*
* The general form for any k where tk is the time at which Nt = N0/k:
* c = (k - 1.0) / (exp(r * tk) - k);
*/
public void setTime50(double time50) {
c = 1.0 / (Math.exp(getGrowthRate() * time50) - 2.0);
}
public void setShapeFromTimeAtAlpha(double time, double alpha) {
// New parameterization of logistic shape to be the time at which the
// population reached some proportion alpha:
double ert = Math.exp(- getGrowthRate() * time);
c = ((1.0 - alpha) * ert) / (ert - alpha);
}
// Implementation of abstract methods
/**
* Gets the value of the demographic function N(t) at time t.
*
* @param t the time
* @return the value of the demographic function N(t) at time t.
*/
public double getDemographic(double t) {
double nZero = getN0();
double r = getGrowthRate();
double c = getShape();
// return nZero * (1 + c) / (1 + (c * Math.exp(r*t)));
// AER rearranging this to use exp(-rt) may help
// with some overflow situations...
double expOfMRT = Math.exp(-r * t);
return lowLimit + (nZero * (1 + c) * expOfMRT) / (c + expOfMRT);
}
public double getLogDemographic(double t) {
final double d = getDemographic(t);
if( d == 0.0 && lowLimit == 0.0 ) {
double nZero = getN0();
double r = getGrowthRate();
double c = getShape();
int sign = c > 0 ? 1 : -1;
final double v1 = Math.log(c * sign) + r * t;
double ld = Math.log(nZero);
if( v1 < 600 ) {
double v = sign * Math.exp(v1);
if( c > -1 ) {
ld += Math.log1p(c) - Math.log1p(v);
} else {
ld += Math.log((1+c)/(1+v));
}
} else {
ld += Math.log1p(c) - sign * v1;
}
return ld;
}
// if( ! (Math.abs(Math.log(d) - ld) < 1e-12) ) {
// return Math.log(d);
// }
return Math.log(d);
}
/**
* Returns value of demographic intensity function at time t
* (= integral 1/N(x) dx from 0 to t).
*/
public double getIntensity(double t) {
double nZero = getN0();
double r = getGrowthRate();
double c = getShape();
double ert = Math.exp(r * t);
if( lowLimit == 0 ) {
// double emrt = Math.exp(-r * t);
return (c * (ert - 1)/r + t) / ((1+c) * nZero);
}
double z = lowLimit;
return (r*t*z + (1 + c)*nZero*Math.log(nZero + c*nZero + z + c*ert*z))/(r*z*(nZero + c*nZero + z));
}
/**
* Returns value of demographic intensity function at time t
* (= integral 1/N(x) dx from 0 to t).
*/
public double getInverseIntensity(double x) {
throw new RuntimeException("Not implemented!");
}
public double getIntegral(double start, double finish) {
if( lowLimit > 0 ) {
double v1 = getNumericalIntegral(start, finish);
final double v2 = getIntensity(finish) - getIntensity(start);
return v2;
}
double intervalLength = finish - start;
double nZero = getN0();
double r = getGrowthRate();
double c = getShape();
double expOfMinusRT = Math.exp(-r * start);
double expOfMinusRG = Math.exp(-r * intervalLength);
double term1 = nZero * (1.0 + c);
if (term1 == 0.0) {
return Double.POSITIVE_INFINITY;
}
double term2 = c * (1.0 - expOfMinusRG);
double term3 = (term1 * expOfMinusRT) * r * expOfMinusRG;
double term2over3;
if (term3 == 0.0) {
double l1 = expOfMinusRG < 1e-8 ? -r * intervalLength : Math.log1p(expOfMinusRG);
final int sign = c > 0 ? 1 : -1;
term2over3 = (sign/term1) * Math.exp(l1 + r * start + Math.log(c*sign) - Math.log(r));
// throw new RuntimeException("Infinite integral!");
} else {
// if (term3 != 0.0 && term2 == 0.0) {
// term2over3 = 0.0;
// } else if (term3 == 0.0 && term2 == 0.0) {
// throw new RuntimeException("term3 and term2 are both zeros. N0=" + getN0() + " growthRate=" + getGrowthRate() + "c=" + c);
// } else {
term2over3 = term2 / term3;
// }
}
final double term5 = intervalLength / term1;
// double v0 = 1/term1 * (finish + c * (Math.exp(r*finish) - 1) /r);
// double v1 = 1/term1 * (start + c * (Math.exp(r*start) - 1) /r);
// double v = 1/term1 * ((finish + c * (Math.exp(r*finish) - 1) /r) - (start + c * (Math.exp(r*start) - 1) /r));
// double v2 = 1/term1 * ((finish-start) + (c * (Math.exp(r*finish) - 1) /r) - (c * (Math.exp(r*start) - 1) /r));
// double v3 = 1/term1 * ((finish-start) + (c/r)* ( Math.exp(r*finish) - 1) - (c/r) * (Math.exp(r*start) - 1) );
// double v4 = 1/term1 * ((finish-start) + (c/r) * (Math.exp(r*finish) - Math.exp(r*start) ) );
// double v = ( (c * (Math.exp(r*finish) - Math.exp(r*start)) / r) + (start - finish)) / term1;
return term5 + term2over3;
}
public int getNumArguments() {
return 3;
}
public String getArgumentName(int n) {
switch (n) {
case 0:
return "N0";
case 1:
return "r";
case 2:
return "c";
}
throw new IllegalArgumentException("Argument " + n + " does not exist");
}
public double getArgument(int n) {
switch (n) {
case 0:
return getN0();
case 1:
return getGrowthRate();
case 2:
return getShape();
}
throw new IllegalArgumentException("Argument " + n + " does not exist");
}
public void setArgument(int n, double value) {
switch (n) {
case 0:
setN0(value);
break;
case 1:
setGrowthRate(value);
break;
case 2:
setShape(value);
break;
default:
throw new IllegalArgumentException("Argument " + n + " does not exist");
}
}
public double getLowerBound(int n) {
return 0.0;
}
public double getUpperBound(int n) {
return Double.POSITIVE_INFINITY;
}
public DemographicFunction getCopy() {
LogisticGrowth df = new LogisticGrowth(getUnits());
df.setN0(getN0());
df.setGrowthRate(getGrowthRate());
df.c = c;
return df;
}
//
// private stuff
//
private double c;
}