use of uk.ac.sussex.gdsc.core.utils.DoubleEquality in project GDSC-SMLM by aherbert.
the class FastMleGradient2ProcedureTest method gradientCalculatorComputesGradient.
private void gradientCalculatorComputesGradient(RandomSeed seed, ErfGaussian2DFunction func) {
// Check the first and second derivatives
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 FastMleGradient2Procedure p = FastMleGradient2ProcedureUtils.create(y, func);
// double ll = p.computeLogLikelihood(a);
p.computeSecondDerivative(a);
final double[] d1 = p.d1.clone();
final double[] d2 = p.d2.clone();
for (int j = 0; j < nparams; j++) {
final int j_ = j;
final int k = indices[j];
final double d = Precision.representableDelta(a[k], (a[k] == 0) ? delta : a[k] * delta);
a2[k] = a[k] + d;
final double llh = p.computeLogLikelihood(a2);
p.computeFirstDerivative(a2);
final double[] d1h = p.d1.clone();
a2[k] = a[k] - d;
final double lll = p.computeLogLikelihood(a2);
p.computeFirstDerivative(a2);
final double[] d1l = p.d1.clone();
a2[k] = a[k];
final double gradient1 = (llh - lll) / (2 * d);
final double gradient2 = (d1h[j] - d1l[j]) / (2 * d);
// logger.fine("[%d,%d] ll - %f (%s %f+/-%f) d1 %f ?= %f : d2 %f ?= %f", i, k, ll,
// func.getName(k), a[k], d,
// gradient1, d1[j], gradient2, d2[j]);
failCounter.run(j, () -> eq.almostEqualRelativeOrAbsolute(gradient1, d1[j_]), () -> {
Assertions.fail(FunctionUtils.getSupplier("Not same gradient1 @ %d,%d: %s != %s (error=%s)", ii, j_, gradient1, d1[j_], DoubleEquality.relativeError(gradient1, d1[j_])));
});
failCounter.run(nparams + j, () -> eq.almostEqualRelativeOrAbsolute(gradient2, d2[j_]), () -> {
Assertions.fail(FunctionUtils.getSupplier("Not same gradient2 @ %d,%d: %s != %s (error=%s)", ii, j_, gradient2, d2[j_], DoubleEquality.relativeError(gradient2, d2[j_])));
});
}
}
}
use of uk.ac.sussex.gdsc.core.utils.DoubleEquality in project GDSC-SMLM by aherbert.
the class ScmosLikelihoodWrapperTest method fitEllipticalComputesGradient.
@SeededTest
void fitEllipticalComputesGradient(RandomSeed seed) {
// The elliptical function gradient evaluation is worse
final DoubleEquality tmp = eq;
eq = eqPerDatum;
functionComputesGradient(seed, GaussianFunctionFactory.FIT_ELLIPTICAL);
eq = tmp;
}
use of uk.ac.sussex.gdsc.core.utils.DoubleEquality in project GDSC-SMLM by aherbert.
the class PoissonLikelihoodWrapperTest method fitNbEllipticalComputesGradient.
@Test
void fitNbEllipticalComputesGradient() {
Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
// The elliptical function gradient evaluation is worse
final DoubleEquality tmp = eq;
eq = eqPerDatum;
functionComputesGradient(GaussianFunctionFactory.FIT_SIMPLE_NB_ELLIPTICAL);
eq = tmp;
}
use of uk.ac.sussex.gdsc.core.utils.DoubleEquality in project GDSC-SMLM by aherbert.
the class PoissonLikelihoodWrapperTest method fitEllipticalComputesGradient.
@Test
void fitEllipticalComputesGradient() {
Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
// The elliptical function gradient evaluation is worse
final DoubleEquality tmp = eq;
eq = eqPerDatum;
functionComputesGradient(GaussianFunctionFactory.FIT_ELLIPTICAL);
eq = tmp;
}
use of uk.ac.sussex.gdsc.core.utils.DoubleEquality in project GDSC-SMLM by aherbert.
the class Gaussian2DFunctionSpeedTest method f1ComputesSameAsf2.
void f1ComputesSameAsf2(RandomSeed seed, int npeaks, int flags1, int flags2) {
final DoubleEquality eq = new DoubleEquality(1e-2, 1e-10);
final int iter = 50;
ArrayList<double[]> paramsList2;
if (npeaks == 1) {
paramsList2 = copyList(ensureDataSingle(seed, iter).paramsListSinglePeak, iter);
} else {
paramsList2 = copyList(ensureDataMulti(seed, iter).paramsListMultiPeak, iter);
}
final Gaussian2DFunction f1 = GaussianFunctionFactory.create2D(1, blockWidth, blockWidth, flags1, null);
final Gaussian2DFunction f2 = GaussianFunctionFactory.create2D(1, blockWidth, blockWidth, flags2, null);
final double[] dyda1 = new double[1 + npeaks * Gaussian2DFunction.PARAMETERS_PER_PEAK];
final double[] dyda2 = new double[1 + npeaks * Gaussian2DFunction.PARAMETERS_PER_PEAK];
final int[] gradientIndices = f1.gradientIndices();
final int[] g1 = new int[gradientIndices.length];
final int[] g2 = new int[gradientIndices.length];
int nparams = 0;
for (int i = 0; i < gradientIndices.length; i++) {
final int index1 = f1.findGradientIndex(g1[i]);
final int index2 = f2.findGradientIndex(g2[i]);
if (index1 >= 0 && index2 >= 0) {
g1[nparams] = index1;
g2[nparams] = index2;
nparams++;
}
}
for (int i = 0; i < paramsList2.size(); i++) {
f1.initialise(paramsList2.get(i));
f2.initialise(paramsList2.get(i));
for (int j = 0; j < x.length; j++) {
final double y1 = f1.eval(x[j], dyda1);
final double y2 = f2.eval(x[j], dyda2);
if (!eq.almostEqualRelativeOrAbsolute(y1, y2)) {
Assertions.fail(String.format("Not same y[%d] @ %d : %g != %g", j, i, y1, y2));
}
for (int ii = 0; ii < nparams; ii++) {
if (!eq.almostEqualRelativeOrAbsolute(dyda1[g1[ii]], dyda2[g2[ii]])) {
Assertions.fail(String.format("Not same dyda[%d] @ %d : %g != %g", j, gradientIndices[g1[ii]], dyda1[g1[ii]], dyda2[g2[ii]]));
}
}
}
}
}
Aggregations