Skip to content

Commit

Permalink
Use Commons Numbers Erf implementation
Browse files Browse the repository at this point in the history
  • Loading branch information
aherbert committed Nov 14, 2024
1 parent 69a55bf commit a1d6a9c
Show file tree
Hide file tree
Showing 6 changed files with 18 additions and 16 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,7 @@ public static double erf0(double x1, double x2) {
* @return the error function erf(x)
*/
public static double erf(double x) {
// return org.apache.commons.math3.special.Erf.erf(x);
// return org.apache.commons.numbers.gamma.Erf.value(x);

final boolean negative = (x < 0);
if (negative) {
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -484,7 +484,7 @@ public double logLikelihood(final double obs, final double exp) {
* @return the cumulative density
*/
double gaussianCdf(final double x) {
// return 0.5 * (1 + org.apache.commons.math3.special.Erf.erf(x / sqrt2sigma2));
// return 0.5 * org.apache.commons.numbers.gamma.Erfc.value(-x / sqrt2sigma2);
// This may not be precise enough.
// Absolute error is <3e-7. Not sure what relative error is.
// The standard Erf is much slower.
Expand All @@ -499,7 +499,7 @@ public double logLikelihood(final double obs, final double exp) {
* @return the cumulative density
*/
double gaussianCdf(final double x, final double x2) {
// return 0.5 * (org.apache.commons.math3.special.Erf.erf(x / sqrt2sigma2, x2 / sqrt2sigma2));
// return 0.5 * org.apache.commons.numbers.gamma.ErfDifference.value(x / sqrt2sigma2, x2 / sqrt2sigma2);
// This may not be precise enough.
// Absolute error is <3e-7. Not sure what relative error is.
// The standard Erf is much slower.
Expand All @@ -513,7 +513,7 @@ public double logLikelihood(final double obs, final double exp) {
* @return the cumulative density
*/
double gaussianErf(final double x) {
// return org.apache.commons.math3.special.Erf.erf(x / sqrt2sigma2);
// return org.apache.commons.numbers.gamma.Erf.value(x / sqrt2sigma2);
// This may not be precise enough.
// Absolute error is <3e-7. Not sure what relative error is.
// The standard Erf is much slower.
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -215,7 +215,7 @@ private double getX(final double D, int q) {
* @return the cumulative density
*/
double gaussianCdf(final double x) {
// return org.apache.commons.math3.special.Erf.erf(x / sqrt_var_by_2)
// return org.apache.commons.numbers.gamma.Erf.value(x / sqrt_var_by_2)
// This may not be precise enough.
// Absolute error is <3e-7. Not sure what relative error is.
// The standard CDF is much slower.
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,9 @@ public enum ErfFunction {
/** Use a fast approximation for the Error function. */
FAST,

/** Use the Apache commons math library Error function. */
/** Use the Apache Commons library Error function.
* <p>As of version 1.2 this uses the implementation in Commons Numbers
* (replacing the Commons Math implementation). */
COMMONS_MATH;
}

Expand Down Expand Up @@ -130,14 +132,14 @@ public double erf(double x) {
}

/**
* Compute the error function using the Commons Math implementation.
* Compute the error function using the Commons Numbers implementation.
*/
private static class CommontsMathErrorFunction implements ErrorFunction {
static final CommontsMathErrorFunction INSTANCE = new CommontsMathErrorFunction();
private static class CommontsNumbersErrorFunction implements ErrorFunction {
static final CommontsNumbersErrorFunction INSTANCE = new CommontsNumbersErrorFunction();

@Override
public double erf(double x) {
return org.apache.commons.math3.special.Erf.erf(x);
return org.apache.commons.numbers.gamma.Erf.value(x);
}
}

Expand Down Expand Up @@ -265,7 +267,7 @@ public ErfFunction getErfFunction() {
public void setErfFunction(ErfFunction erfFunction) {
switch (erfFunction) {
case COMMONS_MATH:
errorFunction = CommontsMathErrorFunction.INSTANCE;
errorFunction = CommontsNumbersErrorFunction.INSTANCE;
break;
case FAST:
errorFunction = FastErrorFunction.INSTANCE;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@
package uk.ac.sussex.gdsc.smlm.model;

import org.apache.commons.lang3.tuple.Pair;
import org.apache.commons.math3.special.Erf;
import org.apache.commons.numbers.gamma.Erf;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler;
import uk.ac.sussex.gdsc.core.utils.MathUtils;
Expand Down Expand Up @@ -289,11 +289,11 @@ public double[] gaussian2D(int x0range, int x1range, double sum, double x0, doub
// Note: The 0.5 factors are moved to reduce computations
for (int x = 0; x <= x0range; x++) {
// erf0[x] = 0.5 * Erf.erf((x - x0) * denom0);
erf0[x] = Erf.erf((x - x0) * denom0);
erf0[x] = Erf.value((x - x0) * denom0);
}
for (int y = 0; y <= x1range; y++) {
// erf1[y] = 0.5 * Erf.erf((y - x1) * denom1);
erf1[y] = Erf.erf((y - x1) * denom1);
erf1[y] = Erf.value((y - x1) * denom1);
}
sum *= 0.5; // Incorporate the 0.5 factor for Y into the sum

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -413,10 +413,10 @@ public static double[] makeErfGaussianKernel(double sigma, double range) {
final double sqrtTwoVar = Math.sqrt(sigma * sigma * 2);

// Use erfc to reach into the long tail at sigma <= 38
double lower = org.apache.commons.math3.special.Erf.erfc(-0.5 / sqrtTwoVar);
double lower = org.apache.commons.numbers.gamma.Erfc.value(-0.5 / sqrtTwoVar);
for (int i = 0; i < kradius; i++) {
final double upper = lower;
lower = org.apache.commons.math3.special.Erf.erfc((i + 0.5) / sqrtTwoVar);
lower = org.apache.commons.numbers.gamma.Erfc.value((i + 0.5) / sqrtTwoVar);
kernel[i] = (upper - lower) * 0.5;
if (kernel[i] == 0) {
break;
Expand Down

0 comments on commit a1d6a9c

Please sign in to comment.