use of uk.ac.sussex.gdsc.core.utils.DoubleEquality in project GDSC-SMLM by aherbert.
the class ScmosLikelihoodWrapperTest method fitNbEllipticalComputesGradient.
@SeededTest
void fitNbEllipticalComputesGradient(RandomSeed seed) {
// The elliptical function gradient evaluation is worse
final DoubleEquality tmp = eq;
eq = eqPerDatum;
functionComputesGradient(seed, GaussianFunctionFactory.FIT_SIMPLE_NB_ELLIPTICAL);
eq = tmp;
}
use of uk.ac.sussex.gdsc.core.utils.DoubleEquality in project GDSC-SMLM by aherbert.
the class LsqLvmGradientProcedureTest method gradientProcedureComputesGradient.
private void gradientProcedureComputesGradient(RandomSeed seed, ErfGaussian2DFunction func) {
final int nparams = func.getNumberOfGradients();
final int[] indices = func.gradientIndices();
final int iter = 100;
final ArrayList<double[]> paramsList = new ArrayList<>(iter);
final ArrayList<double[]> yList = new ArrayList<>(iter);
createData(RngUtils.create(seed.getSeed()), 1, iter, paramsList, yList, true);
// for the gradients
final double delta = 1e-4;
final DoubleEquality eq = new DoubleEquality(5e-2, 1e-16);
// Must compute most of the time
final int failureLimit = TestCounter.computeFailureLimit(iter, 0.1);
final TestCounter failCounter = new TestCounter(failureLimit, nparams);
for (int i = 0; i < paramsList.size(); i++) {
final int ii = i;
final double[] y = yList.get(i);
final double[] a = paramsList.get(i);
final double[] a2 = a.clone();
final BaseLsqLvmGradientProcedure p = LsqLvmGradientProcedureUtils.create(y, func);
p.gradient(a);
// double s = p.ssx;
final double[] beta = p.beta.clone();
for (int j = 0; j < nparams; j++) {
final int jj = j;
final int k = indices[j];
final double d = Precision.representableDelta(a[k], (a[k] == 0) ? 1e-3 : a[k] * delta);
a2[k] = a[k] + d;
p.value(a2);
final double s1 = p.value;
a2[k] = a[k] - d;
p.value(a2);
final double s2 = p.value;
a2[k] = a[k];
// Apply a factor of -2 to compute the actual gradients:
// See Numerical Recipes in C++, 2nd Ed. Equation 15.5.6 for Nonlinear Models
beta[j] *= -2;
final double gradient = (s1 - s2) / (2 * d);
// logger.fine(FunctionUtils.getSupplier("[%d,%d] %f (%s %f+/-%f) %f ?= %f", i, k, s,
// func.getName(k), a[k], d, beta[j],
// gradient);
failCounter.run(j, () -> eq.almostEqualRelativeOrAbsolute(beta[jj], gradient), () -> {
Assertions.fail(() -> String.format("Not same gradient @ %d,%d: %s != %s (error=%s)", ii, jj, beta[jj], gradient, DoubleEquality.relativeError(beta[jj], gradient)));
});
}
}
}
use of uk.ac.sussex.gdsc.core.utils.DoubleEquality in project GDSC-SMLM by aherbert.
the class LvmGradientProcedureTest method gradientProcedureComputesGradient.
@SuppressWarnings("null")
private void gradientProcedureComputesGradient(RandomSeed seed, ErfGaussian2DFunction func, Type type, boolean precomputed) {
final int nparams = func.getNumberOfGradients();
final int[] indices = func.gradientIndices();
final int iter = 100;
final ArrayList<double[]> paramsList = new ArrayList<>(iter);
final ArrayList<double[]> yList = new ArrayList<>(iter);
createData(RngUtils.create(seed.getSeed()), 1, iter, paramsList, yList, true);
// for the gradients
final double delta = 1e-4;
final DoubleEquality eq = new DoubleEquality(5e-2, 1e-16);
final double[] b = (precomputed) ? new double[func.size()] : null;
final FastLog fastLog = type == Type.FAST_LOG_MLE ? getFastLog() : null;
// Must compute most of the time
final int failureLimit = TestCounter.computeFailureLimit(iter, 0.1);
final TestCounter failCounter = new TestCounter(failureLimit, nparams);
for (int i = 0; i < paramsList.size(); i++) {
final int ii = i;
final double[] y = yList.get(i);
final double[] a = paramsList.get(i);
final double[] a2 = a.clone();
LvmGradientProcedure gp;
if (precomputed) {
// Mock fitting part of the function already
for (int j = 0; j < b.length; j++) {
b[j] = y[j] * 0.5;
}
gp = LvmGradientProcedureUtils.create(y, OffsetGradient1Function.wrapGradient1Function(func, b), type, fastLog);
} else {
gp = LvmGradientProcedureUtils.create(y, func, type, fastLog);
}
gp.gradient(a);
// double s = p.value;
final double[] beta = gp.beta.clone();
for (int j = 0; j < nparams; j++) {
final int jj = j;
final int k = indices[j];
// double d = Precision.representableDelta(a[k], (a[k] == 0) ? 1e-3 : a[k] * delta);
final double d = Precision.representableDelta(a[k], delta);
a2[k] = a[k] + d;
gp.value(a2);
final double s1 = gp.value;
a2[k] = a[k] - d;
gp.value(a2);
final double s2 = gp.value;
a2[k] = a[k];
// Apply a factor of -2 to compute the actual gradients:
// See Numerical Recipes in C++, 2nd Ed. Equation 15.5.6 for Nonlinear Models
beta[j] *= -2;
final double gradient = (s1 - s2) / (2 * d);
// logger.fine(FunctionUtils.getSupplier("[%d,%d] %f (%s %f+/-%f) %f ?= %f", i, k, s,
// Gaussian2DFunction.getName(k),
// a[k], d, beta[j], gradient);
failCounter.run(j, () -> eq.almostEqualRelativeOrAbsolute(beta[jj], gradient), () -> {
Assertions.fail(() -> String.format("Not same gradient @ %d,%d: %s != %s (error=%s)", ii, jj, beta[jj], gradient, DoubleEquality.relativeError(beta[jj], gradient)));
});
}
}
}
use of uk.ac.sussex.gdsc.core.utils.DoubleEquality in project GDSC-SMLM by aherbert.
the class ErfTest method analyticErfGradientCorrectForErfApproximation.
@Test
void analyticErfGradientCorrectForErfApproximation() {
final BaseErf erf = new Erf();
final int range = 7;
final int steps = 10000;
final double step = (double) range / steps;
final double delta = 1e-3;
final DoubleEquality eq = new DoubleEquality(5e-4, 1e-6);
for (int i = 0; i < steps; i++) {
final double x = i * step;
final double x1 = x + Precision.representableDelta(x, delta);
final double x2 = x - Precision.representableDelta(x, delta);
final double o1 = erf.erf(x1);
final double o2 = erf.erf(x2);
final double delta2 = x1 - x2;
final double g = (o1 - o2) / delta2;
final double e = uk.ac.sussex.gdsc.smlm.function.Erf.erfDerivative(x);
if (!eq.almostEqualRelativeOrAbsolute(e, g)) {
Assertions.fail(x + " : " + e + " != " + g);
}
}
}
use of uk.ac.sussex.gdsc.core.utils.DoubleEquality in project GDSC-SMLM by aherbert.
the class NonLinearFit method init.
private void init(StoppingCriteria sc, double maxRelativeError, double maxAbsoluteError) {
setStoppingCriteria(sc);
solver.setEqual(new DoubleEquality(maxRelativeError, maxAbsoluteError));
}
Aggregations