All Downloads are FREE. Search and download functionalities are using the official Maven repository.

uk.ac.sussex.gdsc.smlm.function.Erf Maven / Gradle / Ivy

Go to download

Genome Damage and Stability Centre SMLM Package Software for single molecule localisation microscopy (SMLM)

The newest version!
/*-
 * #%L
 * Genome Damage and Stability Centre SMLM Package
 *
 * Software for single molecule localisation microscopy (SMLM)
 * %%
 * Copyright (C) 2011 - 2023 Alex Herbert
 * %%
 * This program 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.
 *
 * This program 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 this program.  If not, see
 * .
 * #L%
 */

package uk.ac.sussex.gdsc.smlm.function;

import uk.ac.sussex.gdsc.smlm.utils.StdMath;

/**
 * Class for computing the error function using approximations.
 *
 * 

Methods in this class are ordered by speed, not accuracy. The default call to * {@link #erf(double)} is a good compromise between both. * * @see https://en.wikipedia.org/wiki/Error_function#Approximation_with_elementary_functions * @see * Winitzki, Sergei (6 February 2008). * "A handy approximation for the error function and its inverse" */ public final class Erf { private static final double FOUR_OVER_PI = 4.0 / Math.PI; private static double DERIVATIVE_FACTOR = 2 / Math.sqrt(Math.PI); /** No public constructor. */ private Erf() {} /** * Returns the error function. * *

erf(x) = 2/√π 0x e-t*tdt

* *

This implementation computes erf(x) using the approximation by Abramowitz and Stegun, (1972) * Handbook of Mathematical Functions. New York: Dover. Formulas 7.1.27, p299. The maximum * absolute error is about 5e-4 for all x.

* *

The value returned is always between -1 and 1 (inclusive). If {@code abs(x) > 40}, then * {@code erf(x)} is indistinguishable from either 1 or -1 as a double, so the appropriate extreme * value is returned.

* * @param x the value. * @return the error function erf(x) */ public static double erf0(double x) { final boolean negative = (x < 0); if (negative) { x = -x; } // Set this to 40 when computing the limit in the JUnit test if (x > 19.581294086464855) { return negative ? -1 : 1; } final double x2 = x * x; final double ret = 1 - 1 / pow4(1.0 + 0.278393 * x + 0.230389 * x2 + 0.000972 * x2 * x + 0.078108 * x2 * x2); return (negative) ? -ret : ret; } /** * Returns the difference between erf(x1) and erf(x2). Uses the fast approximation * {@link #erf0(double)}. * * @param x1 the first value * @param x2 the second value * @return erf(x2) - erf(x1) */ public static double erf0(double x1, double x2) { return erf0(x2) - erf0(x1); } /** * Returns the error function. * *

erf(x) = 2/√π 0x e-t*tdt

* *

This implementation computes erf(x) using the approximation by Abramowitz and Stegun, (1972) * Handbook of Mathematical Functions. New York: Dover. Formulas 7.1.28, p299. The maximum * absolute error is about 3e-7 for all x.

* *

The value returned is always between -1 and 1 (inclusive). If {@code abs(x) > 40}, then * {@code erf(x)} is indistinguishable from either 1 or -1 as a double, so the appropriate extreme * value is returned.

* * @param x the value. * @return the error function erf(x) */ public static double erf(double x) { // return org.apache.commons.math3.special.Erf.erf(x); final boolean negative = (x < 0); if (negative) { x = -x; } // Set this to 40 when computing the limit in the JUnit test if (x > 6.183574750897915) { return negative ? -1 : 1; } final double x2 = x * x; final double x3 = x2 * x; final double ret = 1 - 1 / pow16(1.0 + 0.0705230784 * x + 0.0422820123 * x2 + 0.0092705272 * x3 + 0.0001520143 * x2 * x2 + 0.0002765672 * x2 * x3 + 0.0000430638 * x3 * x3); return (negative) ? -ret : ret; } /** * Returns the difference between erf(x1) and erf(x2). Uses the fast approximation * {@link #erf(double)}. * * @param x1 the first value * @param x2 the second value * @return erf(x2) - erf(x1) */ public static double erf(double x1, double x2) { return erf(x2) - erf(x1); } /** * Returns the error function. * *

erf(x) = 2/√π 0x e-t*tdt

* *

This implementation computes erf(x) using the approximation by Sergei Winitzki. This * involves a sqrt() and exp() function call and so is slower than {@link #erf(double)} and * {@link #erf0(double)}. The maximum absolute error is about 0.00012 for all x.

* *

The value returned is always between -1 and 1 (inclusive). If {@code abs(x) > 40}, then * {@code erf(x)} is indistinguishable from either 1 or -1 as a double, so the appropriate extreme * value is returned.

* * @param x the value. * @return the error function erf(x) */ public static double erf2(double x) { final boolean negative = (x < 0); if (negative) { x = -x; } // Set this to 40 when computing the limit in the JUnit test if (x > 5.9889490707148445) { return negative ? -1 : 1; } final double x2 = x * x; final double ax2 = 0.147 * x2; final double ret = Math.sqrt(-Math.expm1(-x2 * (FOUR_OVER_PI + ax2) / (1 + ax2))); return negative ? -ret : ret; } /** * Returns the difference between erf(x1) and erf(x2). Uses the fast approximation * {@link #erf2(double)}. * * @param x1 the first value * @param x2 the second value * @return erf(x2) - erf(x1) */ public static double erf2(double x1, double x2) { return erf2(x2) - erf2(x1); } /** * Compute the first derivative of the Error function = (2 / sqrt(pi)) * exp(-x*x). * * @see http://mathworld.wolfram.com/Erf.html * * @param x the x * @return the first derivative */ public static double erfDerivative(double x) { return DERIVATIVE_FACTOR * StdMath.exp(-x * x); } /** * Compute the value to the power of 4. * * @param value the value * @return value^4 */ static double pow4(double value) { value = value * value; // power 2 return value * value; } /** * Compute the value to the power of 16. * * @param value the value * @return value^16 */ static double pow16(double value) { value = value * value; // power2 value = value * value; // power4 value = value * value; // power8 return value * value; } }




© 2015 - 2024 Weber Informatics LLC | Privacy Policy