Search in sources :

Example 6 with DoubleEquality

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

Example 7 with DoubleEquality

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

Example 8 with DoubleEquality

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

Example 9 with DoubleEquality

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

Example 10 with DoubleEquality

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]]));
                }
            }
        }
    }
}
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