Java tutorial
/* * MathUtils.java * * BEAST: Bayesian Evolutionary Analysis by Sampling Trees * Copyright (C) 2014 BEAST Developers * * BEAST is free software: you can redistribute it and/or modify it * under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 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 General Public License for more details. * * You should have received a copy of the GNU General Public License * along with BEAST. If not, see <http://www.gnu.org/licenses/>. */ package beast.math; import beast.util.FileHelpers; import beast.util.NumberFormatter; import beast.util.Serializer; import org.apache.commons.math3.util.CombinatoricsUtils; import java.io.File; import java.io.IOException; import java.text.NumberFormat; import java.text.ParseException; /** * Handy utility functions which have some Mathematical relavance. * * @author Matthew Goode * @author Alexei Drummond * @author Gerton Lunter * @version $Id: MathUtils.java,v 1.13 2006/08/31 14:57:24 rambaut Exp $ */ public class MathUtils { private MathUtils() { } /** * A random number generator that is initialized with the clock when this * class is loaded into the JVM. Use this for all random numbers. * Note: This method or getting random numbers in not thread-safe. Since * MersenneTwisterFast is currently (as of 9/01) not synchronized using * this function may cause concurrency issues. Use the static get methods of the * MersenneTwisterFast class for access to a single instance of the class, that * has synchronization. */ private static final MersenneTwisterFast random; private static final Serializer<MersenneTwisterFast> serializer; private static final File defaultStateFile = FileHelpers.getFile("random.state"); static { final boolean resume = Boolean.valueOf(System.getProperty("resume", "false")); if (resume) { try { final String fileName = System.getProperty("random.state", defaultStateFile.getCanonicalPath()); final File stateFile = new File(fileName); serializer = new Serializer<>(stateFile, MersenneTwisterFast.class); random = serializer.getObject(); } catch (Serializer.SerializationException | IOException e) { throw new RuntimeException("Problem resuming random number generator!", e); } } else { random = MersenneTwisterFast.DEFAULT_INSTANCE; try { final String fileName = System.getProperty("random.state", defaultStateFile.getCanonicalPath()); serializer = new Serializer<>(new File(fileName), random); } catch (Serializer.SerializationException | IOException e) { throw new RuntimeException(e); } } } public static void saveState() throws Serializer.SerializationException { serializer.serialize(); } // Chooses one category if a cumulative probability distribution is given public static int randomChoice(double[] cf) { double U = MathUtils.nextDouble(); int s; if (U <= cf[0]) { s = 0; } else { for (s = 1; s < cf.length; s++) { if (U <= cf[s] && U > cf[s - 1]) { break; } } } return s; } /** * @param pdf array of unnormalized probabilities * @return a sample according to an unnormalized probability distribution */ public static int randomChoicePDF(double[] pdf) { double U = MathUtils.nextDouble() * getTotal(pdf); for (int i = 0; i < pdf.length; i++) { U -= pdf[i]; if (U < 0.0) { return i; } } for (int i = 0; i < pdf.length; i++) { System.out.println(i + "\t" + pdf[i]); } throw new Error("randomChoicePDF falls through -- negative, infinite or NaN components in input " + "distribution, or all zeroes?"); } /** * @param logpdf array of unnormalised log probabilities * @return a sample according to an unnormalised probability distribution * * Use this if probabilities are rounding to zero when converted to real space */ public static int randomChoiceLogPDF(double[] logpdf) { double scalingFactor = Double.NEGATIVE_INFINITY; for (double aLogpdf : logpdf) { if (aLogpdf > scalingFactor) { scalingFactor = aLogpdf; } } if (scalingFactor == Double.NEGATIVE_INFINITY) { throw new Error("randomChoiceLogPDF falls through -- all -INF components in input distribution"); } for (int j = 0; j < logpdf.length; j++) { logpdf[j] = logpdf[j] - scalingFactor; } double[] pdf = new double[logpdf.length]; for (int j = 0; j < logpdf.length; j++) { pdf[j] = Math.exp(logpdf[j]); } return randomChoicePDF(pdf); } /** * @param array to normalize * @return a new double array where all the values sum to 1. * Relative ratios are preserved. */ public static double[] getNormalized(double[] array) { double[] newArray = new double[array.length]; double total = getTotal(array); for (int i = 0; i < array.length; i++) { newArray[i] = array[i] / total; } return newArray; } /** * @param array entries to be summed * @param start start position * @param end the index of the element after the last one to be included * @return the total of a the values in a range of an array */ public static double getTotal(double[] array, int start, int end) { double total = 0.0; for (int i = start; i < end; i++) { total += array[i]; } return total; } /** * @param array to sum over * @return the total of the values in an array */ public static double getTotal(double[] array) { return getTotal(array, 0, array.length); } // ===================== (Synchronized) Static access methods to the private random instance =========== /** * Access a default instance of this class, access is synchronized */ public static long getSeed() { synchronized (random) { return random.getSeed(); } } /** * Access a default instance of this class, access is synchronized */ public static void setSeed(long seed) { synchronized (random) { random.setSeed(seed); } } /** * Access a default instance of this class, access is synchronized */ public static byte nextByte() { synchronized (random) { return random.nextByte(); } } /** * Access a default instance of this class, access is synchronized */ public static boolean nextBoolean() { synchronized (random) { return random.nextBoolean(); } } /** * Access a default instance of this class, access is synchronized */ public static void nextBytes(byte[] bs) { synchronized (random) { random.nextBytes(bs); } } /** * Access a default instance of this class, access is synchronized */ public static char nextChar() { synchronized (random) { return random.nextChar(); } } /** * Access a default instance of this class, access is synchronized */ public static double nextGaussian() { synchronized (random) { return random.nextGaussian(); } } //Mean = alpha / lambda //Variance = alpha / (lambda*lambda) public static double nextGamma(double alpha, double lambda) { synchronized (random) { return random.nextGamma(alpha, lambda); } } //Mean = alpha/(alpha+beta) //Variance = (alpha*beta)/(alpha+beta)^2*(alpha+beta+1) public static double nextBeta(double alpha, double beta) { double x = nextGamma(alpha, 1); double y = nextGamma(beta, 1); return x / (x + y); } /** * Access a default instance of this class, access is synchronized * * @return a pseudo random double precision floating point number in [01) */ public static double nextDouble() { synchronized (random) { return random.nextDouble(); } } /** * @return log of random variable in [0,1] */ public static double randomLogDouble() { return Math.log(nextDouble()); } /** * Access a default instance of this class, access is synchronized */ public static double nextExponential(double lambda) { synchronized (random) { return -1.0 * Math.log(1 - random.nextDouble()) / lambda; } } /** * Access a default instance of this class, access is synchronized */ public static double nextInverseGaussian(double mu, double lambda) { synchronized (random) { /* CODE TAKEN FROM WIKIPEDIA. TESTING DONE WITH RESULTS GENERATED IN R AND LOOK COMPARABLE */ double v = random.nextGaussian(); // sample from a normal distribution with a mean of 0 and 1 standard deviation double y = v * v; double x = mu + (mu * mu * y) / (2 * lambda) - (mu / (2 * lambda)) * Math.sqrt(4 * mu * lambda * y + mu * mu * y * y); double test = MathUtils.nextDouble(); // sample from a uniform distribution between 0 and 1 if (test <= (mu) / (mu + x)) { return x; } else { return (mu * mu) / x; } } } /** * Access a default instance of this class, access is synchronized */ public static float nextFloat() { synchronized (random) { return random.nextFloat(); } } /** * Access a default instance of this class, access is synchronized */ public static long nextLong() { synchronized (random) { return random.nextLong(); } } /** * Access a default instance of this class, access is synchronized */ public static short nextShort() { synchronized (random) { return random.nextShort(); } } /** * Access a default instance of this class, access is synchronized */ public static int nextInt() { synchronized (random) { return random.nextInt(); } } /** * Access a default instance of this class, access is synchronized */ public static int nextInt(int n) { synchronized (random) { return random.nextInt(n); } } /** * * @param low * @param high * @return uniform between low and high */ public static double uniform(double low, double high) { return low + nextDouble() * (high - low); } /** * Shuffles an array. */ public static void shuffle(int[] array) { synchronized (random) { random.shuffle(array); } } /** * Shuffles an array. Shuffles numberOfShuffles times */ public static void shuffle(int[] array, int numberOfShuffles) { synchronized (random) { random.shuffle(array, numberOfShuffles); } } /** * Returns an array of shuffled indices of length l. * * @param l length of the array optional. */ public static int[] shuffled(int l) { synchronized (random) { return random.shuffled(l); } } public static int[] sampleIndicesWithReplacement(int length) { synchronized (random) { int[] result = new int[length]; for (int i = 0; i < length; i++) result[i] = random.nextInt(length); return result; } } /** * Permutes an array. */ public static void permute(int[] array) { synchronized (random) { random.permute(array); } } /** * Returns a uniform random permutation of 0,...,l-1 * * @param l length of the array optional. */ public static int[] permuted(int l) { synchronized (random) { return random.permuted(l); } } /** * Returns sqrt(a^2 + b^2) without under/overflow. */ public static double hypot(double a, double b) { double r; if (Math.abs(a) > Math.abs(b)) { r = b / a; r = Math.abs(a) * Math.sqrt(1 + r * r); } else if (b != 0) { r = a / b; r = Math.abs(b) * Math.sqrt(1 + r * r); } else { r = 0.0; } return r; } /** * return double *.???? * @param value * @param sf * @return */ public static double round(double value, int sf) { NumberFormatter formatter = new NumberFormatter(sf); try { return NumberFormat.getInstance().parse(formatter.format(value)).doubleValue(); } catch (ParseException e) { return value; } } public static double nChoose2(int n) { if (n < 2) return 0.0; else return CombinatoricsUtils.binomialCoefficientDouble(n, 2); } }