Search in sources :

Example 11 with DoubleEquality

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;
}
Also used : DoubleEquality(uk.ac.sussex.gdsc.core.utils.DoubleEquality) SeededTest(uk.ac.sussex.gdsc.test.junit5.SeededTest)

Example 12 with DoubleEquality

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)));
            });
        }
    }
}
Also used : TestCounter(uk.ac.sussex.gdsc.test.utils.TestCounter) ArrayList(java.util.ArrayList) DoubleEquality(uk.ac.sussex.gdsc.core.utils.DoubleEquality)

Example 13 with DoubleEquality

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)));
            });
        }
    }
}
Also used : TestCounter(uk.ac.sussex.gdsc.test.utils.TestCounter) ArrayList(java.util.ArrayList) DoubleEquality(uk.ac.sussex.gdsc.core.utils.DoubleEquality) FastLog(uk.ac.sussex.gdsc.smlm.function.FastLog)

Example 14 with DoubleEquality

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);
        }
    }
}
Also used : DoubleEquality(uk.ac.sussex.gdsc.core.utils.DoubleEquality) SeededTest(uk.ac.sussex.gdsc.test.junit5.SeededTest) Test(org.junit.jupiter.api.Test)

Example 15 with DoubleEquality

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));
}
Also used : DoubleEquality(uk.ac.sussex.gdsc.core.utils.DoubleEquality)

Aggregations

DoubleEquality (uk.ac.sussex.gdsc.core.utils.DoubleEquality)15 ArrayList (java.util.ArrayList)6 SeededTest (uk.ac.sussex.gdsc.test.junit5.SeededTest)5 TestCounter (uk.ac.sussex.gdsc.test.utils.TestCounter)5 Test (org.junit.jupiter.api.Test)4 UniformRandomProvider (org.apache.commons.rng.UniformRandomProvider)2 ParameterizedTest (org.junit.jupiter.params.ParameterizedTest)2 FastLog (uk.ac.sussex.gdsc.smlm.function.FastLog)2 Gaussian2DFunction (uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunction)2 SharedStateContinuousSampler (org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler)1 DenseMatrix64F (org.ejml.data.DenseMatrix64F)1 ValueProcedure (uk.ac.sussex.gdsc.smlm.function.ValueProcedure)1 EllipticalGaussian2DFunction (uk.ac.sussex.gdsc.smlm.function.gaussian.EllipticalGaussian2DFunction)1 SingleCircularGaussian2DFunction (uk.ac.sussex.gdsc.smlm.function.gaussian.SingleCircularGaussian2DFunction)1 SingleEllipticalGaussian2DFunction (uk.ac.sussex.gdsc.smlm.function.gaussian.SingleEllipticalGaussian2DFunction)1 SingleFixedGaussian2DFunction (uk.ac.sussex.gdsc.smlm.function.gaussian.SingleFixedGaussian2DFunction)1 SingleFreeCircularGaussian2DFunction (uk.ac.sussex.gdsc.smlm.function.gaussian.SingleFreeCircularGaussian2DFunction)1 SingleNbFixedGaussian2DFunction (uk.ac.sussex.gdsc.smlm.function.gaussian.SingleNbFixedGaussian2DFunction)1 ErfGaussian2DFunction (uk.ac.sussex.gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction)1 SingleFreeCircularErfGaussian2DFunction (uk.ac.sussex.gdsc.smlm.function.gaussian.erf.SingleFreeCircularErfGaussian2DFunction)1