Examples with Gaussian2DFunction uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunction used on opensource projects

Example 11 with Gaussian2DFunction

use of uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunction in project GDSC-SMLM by aherbert.

the class GradientCalculatorSpeedTest method gradientCalculatorComputesGradient.

private void gradientCalculatorComputesGradient(RandomSeed seed, GradientCalculator calc) {
    final int nparams = calc.nparams;
    final Gaussian2DFunction func = new SingleEllipticalGaussian2DFunction(blockWidth, blockWidth);
    // Check the function is the correct size
    final int[] indices = func.gradientIndices();
    Assertions.assertEquals(nparams, indices.length);
    final int iter = 50;
    final double[] beta = new double[nparams];
    final double[] beta2 = new double[nparams];
    final ArrayList<double[]> paramsList = new ArrayList<>(iter);
    final ArrayList<double[]> yList = new ArrayList<>(iter);
    final int[] x = createData(RngUtils.create(seed.getSeed()), 1, iter, paramsList, yList, true);
    final double delta = 1e-3;
    final DoubleEquality eq = new DoubleEquality(1e-3, 1e-3);
    final IntArrayFormatSupplier msg = new IntArrayFormatSupplier("[%d] Not same gradient @ %d", 2);
    for (int i = 0; i < paramsList.size(); i++) {
        msg.set(0, i);
        final double[] y = yList.get(i);
        final double[] a = paramsList.get(i);
        final double[] a2 = a.clone();
        // double s =
        calc.evaluate(x, y, a, beta, func);
        for (int k = 0; k < nparams; k++) {
            final int j = indices[k];
            final double d = Precision.representableDelta(a[j], (a[j] == 0) ? 1e-3 : a[j] * delta);
            a2[j] = a[j] + d;
            final double s1 = calc.evaluate(x, y, a2, beta2, func);
            a2[j] = a[j] - d;
            final double s2 = calc.evaluate(x, y, a2, beta2, func);
            a2[j] = a[j];
            final double gradient = (s1 - s2) / (2 * d);
            // logger.fine(FunctionUtils.getSupplier("[%d,%d] %f (%s %f+/-%f) %f ?= %f", i, j, s,
            // func.getName(j), a[j], d, beta[k],
            // gradient));
            Assertions.assertTrue(eq.almostEqualRelativeOrAbsolute(beta[k], gradient), msg.set(1, j));
        }
    }
}

Example 12 with Gaussian2DFunction

use of uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunction in project GDSC-SMLM by aherbert.

the class LsqVarianceGradientProcedureTest method varianceMatchesFormula.

@SeededTest
void varianceMatchesFormula() {
    // Assumptions.assumeTrue(false);
    final double[] n_ = new double[] { 20, 50, 100, 500 };
    final double[] b2_ = new double[] { 0, 1, 2, 4 };
    final double[] s_ = new double[] { 1, 1.2, 1.5 };
    final double[] x_ = new double[] { 4.8, 5, 5.5 };
    final double a = 100;
    final int size = 10;
    final Gaussian2DFunction f = GaussianFunctionFactory.create2D(1, size, size, GaussianFunctionFactory.FIT_ERF_CIRCLE, null);
    final LsqVarianceGradientProcedure p = LsqVarianceGradientProcedureUtils.create(f);
    final int ix = f.findGradientIndex(Gaussian2DFunction.X_POSITION);
    final int iy = f.findGradientIndex(Gaussian2DFunction.Y_POSITION);
    final double[] params = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    for (final double n : n_) {
        params[Gaussian2DFunction.SIGNAL] = n;
        for (final double b2 : b2_) {
            params[Gaussian2DFunction.BACKGROUND] = b2;
            for (final double s : s_) {
                final double ss = s * a;
                params[Gaussian2DFunction.X_SD] = s;
                for (final double x : x_) {
                    params[Gaussian2DFunction.X_POSITION] = x;
                    for (final double y : x_) {
                        params[Gaussian2DFunction.Y_POSITION] = y;
                        if (p.variance(params) != LsqVarianceGradientProcedure.STATUS_OK) {
                            Assertions.fail("No variance");
                        }
                        final double o1 = Math.sqrt(p.variance[ix]) * a;
                        final double o2 = Math.sqrt(p.variance[iy]) * a;
                        final double e = Gaussian2DPeakResultHelper.getPrecisionX(a, ss, n, b2, false);
                        // logger.fine(FunctionUtils.getSupplier("e = %f : o = %f %f", e, o1, o2);
                        Assertions.assertEquals(e, o1, e * 5e-2);
                        Assertions.assertEquals(e, o2, e * 5e-2);
                    }
                }
            }
        }
    }
}

Example 13 with Gaussian2DFunction

use of uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunction in project GDSC-SMLM by aherbert.

the class CubicSplineFunctionTest method speedTest.

@SuppressWarnings("null")
private void speedTest(int n, int order) {
    // No assertions, this is just a report
    Assumptions.assumeTrue(logger.isLoggable(Level.INFO));
    // Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
    final CubicSplineFunction cf = (n == 2) ? f2 : f1;
    Assumptions.assumeTrue(null != cf);
    final CubicSplineFunction cff = (n == 2) ? f2f : f1f;
    final ErfGaussian2DFunction gf = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(n, maxx, maxy, GaussianFunctionFactory.FIT_ASTIGMATISM, zModel);
    final Gaussian2DFunction gf2 = (order < 2) ? GaussianFunctionFactory.create2D(n, maxx, maxy, GaussianFunctionFactory.FIT_SIMPLE_FREE_CIRCLE, zModel) : null;
    final LocalList<double[]> l1 = new LocalList<>();
    final LocalList<double[]> l2 = new LocalList<>();
    final LocalList<double[]> l3 = new LocalList<>();
    final double[] a = new double[1 + n * CubicSplineFunction.PARAMETERS_PER_PEAK];
    final double[] b = new double[1 + n * Gaussian2DFunction.PARAMETERS_PER_PEAK];
    double[] bb = null;
    a[CubicSplineFunction.BACKGROUND] = 0.1;
    b[Gaussian2DFunction.BACKGROUND] = 0.1;
    for (int i = 0; i < n; i++) {
        a[i * CubicSplineFunction.PARAMETERS_PER_PEAK + CubicSplineFunction.SIGNAL] = 10;
        b[i * Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.SIGNAL] = 10;
    }
    if (n == 2) {
        // Fix second peak parameters
        a[CubicSplineFunction.PARAMETERS_PER_PEAK + CubicSplineFunction.X_POSITION] = testcx1[0];
        a[CubicSplineFunction.PARAMETERS_PER_PEAK + CubicSplineFunction.Y_POSITION] = testcy1[0];
        a[CubicSplineFunction.PARAMETERS_PER_PEAK + CubicSplineFunction.Z_POSITION] = testcz1[0];
        b[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.X_POSITION] = testcx1[0];
        b[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Y_POSITION] = testcy1[0];
        b[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Z_POSITION] = testcz1[0];
    }
    if (gf2 != null) {
        bb = b.clone();
        if (n == 2) {
            // Fix second peak parameters
            bb[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.X_SD] = zModel.getSx(testcz1[0]);
            bb[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Y_SD] = zModel.getSy(testcz1[0]);
        }
    }
    for (int x = 0; x <= maxx; x++) {
        a[CubicSplineFunction.X_POSITION] = x;
        b[Gaussian2DFunction.X_POSITION] = x;
        for (int y = 0; y <= maxy; y++) {
            a[CubicSplineFunction.Y_POSITION] = y;
            b[Gaussian2DFunction.Y_POSITION] = y;
            for (int z = -zDepth; z <= zDepth; z++) {
                a[CubicSplineFunction.Z_POSITION] = z;
                b[Gaussian2DFunction.Z_POSITION] = z;
                l1.add(a.clone());
                l2.add(b.clone());
                if (gf2 != null) {
                    bb[Gaussian2DFunction.X_SD] = zModel.getSx(z);
                    bb[Gaussian2DFunction.Y_SD] = zModel.getSy(z);
                    l3.add(bb.clone());
                }
            }
        }
    }
    final double[][] x1 = l1.toArray(new double[0][]);
    final double[][] x2 = l2.toArray(new double[0][]);
    final double[][] x3 = l3.toArray(new double[0][]);
    final TimingService ts = new TimingService(5);
    ts.execute(new FunctionTimingTask(gf, x2, order));
    if (gf2 != null) {
        ts.execute(new FunctionTimingTask(gf2, x3, order));
    }
    ts.execute(new FunctionTimingTask(cf, x1, order));
    ts.execute(new FunctionTimingTask(cff, x1, order, " single-precision"));
    final int size = ts.getSize();
    ts.repeat(size);
    logger.info(ts.getReport(size));
}

Example 14 with Gaussian2DFunction

use of uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunction in project GDSC-SMLM by aherbert.

the class ScmosLikelihoodWrapperTest method functionComputesGradient.

private void functionComputesGradient(RandomSeed seed, int flags) {
    final Gaussian2DFunction f1 = GaussianFunctionFactory.create2D(1, maxx, maxx, flags, null);
    // Setup
    testbackground = testbackgroundOptions;
    testsignal1 = testsignal1Options;
    testcx1 = testcx1Options;
    testcy1 = testcy1Options;
    testcz1 = testcz1Options;
    testw1 = testw1Options;
    testangle1 = testangle1Options;
    if (!f1.evaluatesBackground()) {
        testbackground = new double[] { testbackground[0] };
    }
    if (!f1.evaluatesSignal()) {
        testsignal1 = new double[] { testsignal1[0] };
    }
    if (!f1.evaluatesZ()) {
        testcz1 = new double[] { 0 };
    }
    boolean noSecondWidth = false;
    if (!f1.evaluatesSD0()) {
        // Just use 1 width
        testw1 = new double[][] { testw1[0] };
        // If no width 0 then assume we have no width 1 as well
        noSecondWidth = true;
    } else if (!f1.evaluatesSD1()) {
        // No evaluation of second width needs only variation in width 0 so truncate
        testw1 = Arrays.copyOf(testw1, 2);
        noSecondWidth = true;
    }
    if (noSecondWidth) {
        for (int i = 0; i < testw1.length; i++) {
            testw1[i][1] = testw1[i][0];
        }
    }
    if (!f1.evaluatesAngle()) {
        testangle1 = new double[] { 0 };
    }
    final double fraction = 85;
    if (f1.evaluatesBackground()) {
        functionComputesTargetGradient(seed, f1, Gaussian2DFunction.BACKGROUND, fraction);
    }
    if (f1.evaluatesSignal()) {
        functionComputesTargetGradient(seed, f1, Gaussian2DFunction.SIGNAL, fraction);
    }
    functionComputesTargetGradient(seed, f1, Gaussian2DFunction.X_POSITION, fraction);
    functionComputesTargetGradient(seed, f1, Gaussian2DFunction.Y_POSITION, fraction);
    if (f1.evaluatesZ()) {
        functionComputesTargetGradient(seed, f1, Gaussian2DFunction.Z_POSITION, fraction);
    }
    if (f1.evaluatesSD0()) {
        functionComputesTargetGradient(seed, f1, Gaussian2DFunction.X_SD, fraction);
    }
    if (f1.evaluatesSD1()) {
        functionComputesTargetGradient(seed, f1, Gaussian2DFunction.Y_SD, fraction);
    }
    if (f1.evaluatesAngle()) {
        functionComputesTargetGradient(seed, f1, Gaussian2DFunction.ANGLE, fraction);
    }
}

Example 15 with Gaussian2DFunction

use of uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunction in project GDSC-SMLM by aherbert.

the class PoissonLikelihoodWrapperTest method functionComputesGradientPerDatum.

private void functionComputesGradientPerDatum(int flags) {
    final Gaussian2DFunction f1 = GaussianFunctionFactory.create2D(1, maxx, maxx, flags, null);
    // Setup
    testbackground = testbackgroundOptions;
    testsignal1 = testsignal1Options;
    testcx1 = testcx1Options;
    testcy1 = testcy1Options;
    testcz1 = testcz1Options;
    testw1 = testw1Options;
    testangle1 = testangle1Options;
    if (!f1.evaluatesBackground()) {
        testbackground = new double[] { testbackground[0] };
    }
    if (!f1.evaluatesSignal()) {
        testsignal1 = new double[] { testsignal1[0] };
    }
    if (!f1.evaluatesZ()) {
        testcz1 = new double[] { 0 };
    }
    boolean noSecondWidth = false;
    if (!f1.evaluatesSD0()) {
        // Just use 1 width
        testw1 = new double[][] { testw1[0] };
        // If no width 0 then assume we have no width 1 as well
        noSecondWidth = true;
    } else if (!f1.evaluatesSD1()) {
        // No evaluation of second width needs only variation in width 0 so truncate
        testw1 = Arrays.copyOf(testw1, 2);
        noSecondWidth = true;
    }
    if (noSecondWidth) {
        for (int i = 0; i < testw1.length; i++) {
            testw1[i][1] = testw1[i][0];
        }
    }
    if (!f1.evaluatesAngle()) {
        testangle1 = new double[] { 0 };
    }
    if (f1.evaluatesBackground()) {
        functionComputesTargetGradientPerDatum(f1, Gaussian2DFunction.BACKGROUND);
    }
    if (f1.evaluatesSignal()) {
        functionComputesTargetGradientPerDatum(f1, Gaussian2DFunction.SIGNAL);
    }
    functionComputesTargetGradientPerDatum(f1, Gaussian2DFunction.X_POSITION);
    functionComputesTargetGradientPerDatum(f1, Gaussian2DFunction.Y_POSITION);
    if (f1.evaluatesZ()) {
        functionComputesTargetGradientPerDatum(f1, Gaussian2DFunction.Z_POSITION);
    }
    if (f1.evaluatesSD0()) {
        functionComputesTargetGradientPerDatum(f1, Gaussian2DFunction.X_SD);
    }
    if (f1.evaluatesSD1()) {
        functionComputesTargetGradientPerDatum(f1, Gaussian2DFunction.Y_SD);
    }
    if (f1.evaluatesAngle()) {
        functionComputesTargetGradientPerDatum(f1, Gaussian2DFunction.ANGLE);
    }
}