Java tutorial
/* * DemographicFunction.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; import dr.evolution.util.Units; import dr.math.Binomial; import dr.math.MathUtils; import org.apache.commons.math.FunctionEvaluationException; import org.apache.commons.math.MaxIterationsExceededException; import org.apache.commons.math.analysis.UnivariateRealFunction; import org.apache.commons.math.analysis.integration.RombergIntegrator; /** * This interface provides methods that describe a demographic function. * * Parts of this class were derived from C++ code provided by Oliver Pybus. * * @version $Id: DemographicFunction.java,v 1.12 2005/05/24 20:25:55 rambaut Exp $ * * @author Andrew Rambaut * @author Alexei Drummond * @author Korbinian Strimmer */ public interface DemographicFunction extends UnivariateRealFunction, Units { /** * @param t time * @return value of the demographic function N(t) at time t */ double getDemographic(double t); double getLogDemographic(double t); /** * @return value of demographic intensity function at time t (= integral 1/N(x) dx from 0 to t). * @param t time */ double getIntensity(double t); /** * @return value of inverse demographic intensity function * (returns time, needed for simulation of coalescent intervals). */ double getInverseIntensity(double x); /** * Calculates the integral 1/N(x) dx between start and finish. * @param start point * @param finish point * @return integral value */ double getIntegral(double start, double finish); /** * @return the number of arguments for this function. */ int getNumArguments(); /** * @return the name of the n'th argument of this function. */ String getArgumentName(int n); /** * @return the value of the n'th argument of this function. */ double getArgument(int n); /** * Sets the value of the nth argument of this function. */ void setArgument(int n, double value); /** * @return the lower bound of the nth argument of this function. */ double getLowerBound(int n); /** * Returns the upper bound of the nth argument of this function. */ double getUpperBound(int n); /** * Returns a copy of this function. */ // DemographicFunction getCopy(); /** * A threshold for underflow on calculation of likelihood of internode intervals. * Most demo functions could probably return 0.0 but (e.g.,) the Extended Skyline * needs a non zero value to prevent a numerical problem. * @return */ double getThreshold(); public abstract class Abstract implements DemographicFunction { // private static final double LARGE_POSITIVE_NUMBER = 1.0e50; // private static final double LARGE_NEGATIVE_NUMBER = -1.0e50; // private static final double INTEGRATION_PRECISION = 1.0e-5; // private static final double INTEGRATION_MAX_ITERATIONS = 50; RombergIntegrator numericalIntegrator = null; /** * Construct demographic model with default settings */ public Abstract(Type units) { setUnits(units); } // general functions /** * Default implementation * @param t * @return log(demographic at t) */ public double getLogDemographic(double t) { return Math.log(getDemographic(t)); } public double getThreshold() { return 0; } /** * Calculates the integral 1/N(x) dx between start and finish. */ public double getIntegral(double start, double finish) { return getIntensity(finish) - getIntensity(start); } /** * Returns the integral of 1/N(x) between start and finish, calling either the getAnalyticalIntegral or * getNumericalIntegral function as appropriate. */ public double getNumericalIntegral(double start, double finish) { // AER 19th March 2008: I switched this to use the RombergIntegrator from // commons-math v1.2. if (start > finish) { throw new RuntimeException("NumericalIntegration start > finish"); } if (start == finish) { return 0.0; } if (numericalIntegrator == null) { numericalIntegrator = new RombergIntegrator(this); } try { return numericalIntegrator.integrate(start, finish); } catch (MaxIterationsExceededException e) { throw new RuntimeException(e); } catch (FunctionEvaluationException e) { throw new RuntimeException(e); } // double lastST = LARGE_NEGATIVE_NUMBER; // double lastS = LARGE_NEGATIVE_NUMBER; // // assert(finish > start); // // for (int j = 1; j <= INTEGRATION_MAX_ITERATIONS; j++) { // // iterate doTrapezoid() until answer obtained // // double st = doTrapezoid(j, start, finish, lastST); // double s = (4.0 * st - lastST) / 3.0; // // // If answer is within desired accuracy then return // if (Math.abs(s - lastS) < INTEGRATION_PRECISION * Math.abs(lastS)) { // return s; // } // lastS = s; // lastST = st; // } // // throw new RuntimeException("Too many iterations in getNumericalIntegral"); } /** * Performs the trapezoid rule. */ // private double doTrapezoid(int n, double low, double high, double lastS) { // // double s; // // if (n == 1) { // // On the first iteration s is reset // double demoLow = getDemographic(low); // Value of N(x) obtained here // assert(demoLow > 0.0); // // double demoHigh = getDemographic(high); // assert(demoHigh > 0.0); // // s = 0.5 * (high - low) * ( (1.0 / demoLow) + (1.0 / demoHigh) ); // } else { // int it=1; // for (int j = 1; j < n - 1; j++) { // it *= 2; // } // // double tnm = it; // number of points // double del = (high - low) / tnm; // width of spacing between points // // double x = low + 0.5 * del; // // double sum = 0.0; // for (int j = 1; j <= it; j++) { // double demoX = getDemographic(x); // Value of N(x) obtained here // assert(demoX > 0.0); // // sum += (1.0 / demoX); // x += del; // } // s = 0.5 * (lastS + (high - low) * sum / tnm); // New s uses previous s value // } // // return s; // } // ************************************************************** // UnivariateRealFunction IMPLEMENTATION // ************************************************************** /** * Return the intensity at a given time for numerical integration * @param x the time * @return the intensity */ public double value(double x) { return 1.0 / getDemographic(x); } // ************************************************************** // Units IMPLEMENTATION // ************************************************************** /** * Units in which population size is measured. */ private Type units; /** * sets units of measurement. * * @param u units */ public void setUnits(Type u) { units = u; } /** * returns units of measurement. */ public Type getUnits() { return units; } } public static class Utils { private static double getInterval(double U, DemographicFunction demographicFunction, int lineageCount, double timeOfLastCoalescent) { final double intensity = demographicFunction.getIntensity(timeOfLastCoalescent); final double tmp = -Math.log(U) / Binomial.choose2(lineageCount) + intensity; return demographicFunction.getInverseIntensity(tmp) - timeOfLastCoalescent; } private static double getInterval(double U, DemographicFunction demographicFunction, int lineageCount, double timeOfLastCoalescent, double earliestTimeOfFinalCoalescent) { if (timeOfLastCoalescent > earliestTimeOfFinalCoalescent) { throw new IllegalArgumentException( "Given maximum height is smaller than given final coalescent time"); } final double fullIntegral = demographicFunction.getIntegral(timeOfLastCoalescent, earliestTimeOfFinalCoalescent); final double normalisation = 1 - Math.exp(-Binomial.choose2(lineageCount) * fullIntegral); final double intensity = demographicFunction.getIntensity(timeOfLastCoalescent); double tmp = -Math.log(1 - U * normalisation) / Binomial.choose2(lineageCount) + intensity; return demographicFunction.getInverseIntensity(tmp) - timeOfLastCoalescent; } /** * @return a random interval size selected from the Kingman prior of the demographic model. */ public static double getSimulatedInterval(DemographicFunction demographicFunction, int lineageCount, double timeOfLastCoalescent) { final double U = MathUtils.nextDouble(); // create unit uniform random variate return getInterval(U, demographicFunction, lineageCount, timeOfLastCoalescent); } public static double getSimulatedInterval(DemographicFunction demographicFunction, int lineageCount, double timeOfLastCoalescent, double earliestTimeOfFirstCoalescent) { final double U = MathUtils.nextDouble(); return getInterval(U, demographicFunction, lineageCount, timeOfLastCoalescent, earliestTimeOfFirstCoalescent); } public static double getMedianInterval(DemographicFunction demographicFunction, int lineageCount, double timeOfLastCoalescent) { return getInterval(0.5, demographicFunction, lineageCount, timeOfLastCoalescent); } /** * This function tests the consistency of the * getIntensity and getInverseIntensity methods * of this demographic model. If the model is * inconsistent then a RuntimeException will be thrown. * @param demographicFunction the demographic model to test. * @param steps the number of steps between 0.0 and maxTime to test. * @param maxTime the maximum time to test. */ public static void testConsistency(DemographicFunction demographicFunction, int steps, double maxTime) { double delta = maxTime / (double) steps; for (int i = 0; i <= steps; i++) { double time = (double) i * delta; double intensity = demographicFunction.getIntensity(time); double newTime = demographicFunction.getInverseIntensity(intensity); if (Math.abs(time - newTime) > 1e-12) { throw new RuntimeException( "Demographic model not consistent! error size = " + Math.abs(time - newTime)); } } } } }