<commit_before>/*
 * LogisticGrowthN0.java
 *
 * <NAME> 4th October 2011
 *
 */

package dr.evomodel.epidemiology;

import dr.evolution.coalescent.*;

/**
 * This class models logistic growth.
 *
 * @author <NAME>
 * @author <NAME>
 * @author <NAME>
 * @version $Id: LogisticGrowth.java,v 1.15 2008/03/21 20:25:56 rambaut Exp $
 */
public class LogisticGrowthN0 extends ExponentialGrowth {

    /**
     * Construct demographic model with default settings
     */
    public LogisticGrowthN0(Type units) {
        super(units);
    }

	public void setT50(double value) {
		t50 = value;
	}
	
	public double getT50() {
		return t50;
	}

    // 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 N0 = getN0();
        double r = getGrowthRate();
        double T50 = getT50();
		
		return N0 * (1 + Math.exp(-r * T50)) / (1 + Math.exp(r * (t-T50)));
    }

    public double getLogDemographic(double t) {
		return Math.log(getDemographic(t));
	}

    /**
     * Returns value of demographic intensity function at time t
     * (= integral 1/N(x) dx from 0 to t).
     */
    public double getIntensity(double t) {
        double N0 = getN0();
        double r = getGrowthRate();
        double T50 = getT50();
		double exp_rT50 = Math.exp(-r*T50);
		
		return (t + exp_rT50 * (Math.exp(r * t) - 1)/r) / (N0 * (1 + exp_rT50));
    }

    /**
     * Returns the inverse function of getIntensity
     */
    public double getInverseIntensity(double x) {
        throw new RuntimeException("Not implemented!");
    }

    public double getIntegral(double start, double finish) {
		return getIntensity(finish) - getIntensity(start);
    }

    public int getNumArguments() {
        return 3;
    }

    public String getArgumentName(int n) {
        switch (n) {
            case 0:
                return "N0";
            case 1:
                return "r";
            case 2:
                return "t50";
        }<commit_msg>Implemented getInverseIntensity in LogisticGrowthN0 - should be able to simulate intervals using a logistic model now, but needs testing.<commit_after>/*
 * LogisticGrowthN0.java
 *
 * <NAME> 4th October 2011
 *
 */

package dr.evomodel.epidemiology;

import dr.evolution.coalescent.*;

/**
 * This class models logistic growth.
 *
 * @author <NAME>
 * @author <NAME>
 * @author <NAME>
 * @version $Id: LogisticGrowth.java,v 1.15 2008/03/21 20:25:56 rambaut Exp $
 */
public class LogisticGrowthN0 extends ExponentialGrowth {

    /**
     * Construct demographic model with default settings
     */
    public LogisticGrowthN0(Type units) {
        super(units);
    }

	public void setT50(double value) {
		t50 = value;
	}
	
	public double getT50() {
		return t50;
	}

    // 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 N0 = getN0();
        double r = getGrowthRate();
        double T50 = getT50();
		
		return N0 * (1 + Math.exp(-r * T50)) / (1 + Math.exp(r * (t-T50)));
    }

    public double getLogDemographic(double t) {
		return Math.log(getDemographic(t));
	}

    /**
     * Returns value of demographic intensity function at time t
     * (= integral 1/N(x) dx from 0 to t).
     */
    public double getIntensity(double t) {
        double N0 = getN0();
        double r = getGrowthRate();
        double T50 = getT50();
		double exp_rT50 = Math.exp(-r*T50);
		
		return (t + exp_rT50 * (Math.exp(r * t) - 1)/r) / (N0 * (1 + exp_rT50));
    }

    /**
     * Returns the inverse function of getIntensity
     *
     * If exp(-qt) = a(x-k) then t = k + (1/q) * W(q*exp(-q*k)/a) where W(x) is the Lambert W function.
     *
     * for our purposes:
     *
     * q = -r
     * a = (q/exp(q*T50))
     * k = N0*(1+exp(q*T50))*x - (1/q)exp(q*T50)
     *
     * For large x, W0(x) is approximately equal to ln(x) - ln(ln(x)); if q*exp(-q*k)/a rounds to infinity, we log it
     * and use this instead
     */
    public double getInverseIntensity(double x) {

        double q = -getGrowthRate();
        double T50 = getT50();
        double N0 = getN0();
        double a = (q/Math.exp(q*T50));
        double k = N0*(1+Math.exp(q*T50))*x - (1/q)*Math.exp(q*T50);

        double lambertInput = q*Math.exp(-q*k)/a;

        double lambertResult;

        if(lambertInput==Double.POSITIVE_INFINITY){

            //use the asymptote

            double logInput = Math.log(-q)-Math.log(-a)-q*k;
            lambertResult = logInput - Math.log(logInput);

        } else {
            lambertResult = LambertW.branch0(lambertInput);
        }

        return k + (1/q)*lambertResult;

    }

    public double getIntegral(double start, double finish) {
		return getIntensity(finish) - getIntensity(start);
    }

    public int getNumArguments() {
        return 3;
    }

    public String getArgumentName(int n) {
        switch (n) {
            case 0:
                return "N0";
            case 1:
                return "r";
            case 2:
                return "t50";
        }
{'fraction_non_alphanumeric': 0.12523191094619665, 'fraction_numerical': 0.03209647495361781, 'mean_word_length': 2.1416083916083917, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5894602', 'n_tokens_mistral': 1960, 'n_tokens_neox': 1758, 'n_words': 495}