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

Example 1 with ErfGaussian2DFunction

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

the class LsqLvmGradientProcedureTest method gradientProcedureComputesSameOutputWithBias.

@SeededTest
void gradientProcedureComputesSameOutputWithBias(RandomSeed seed) {
    final ErfGaussian2DFunction func = new SingleFreeCircularErfGaussian2DFunction(blockWidth, blockWidth);
    final int nparams = func.getNumberOfGradients();
    final int iter = 100;
    final Level logLevel = Level.FINER;
    final boolean debug = logger.isLoggable(logLevel);
    final ArrayList<double[]> paramsList = new ArrayList<>(iter);
    final ArrayList<double[]> yList = new ArrayList<>(iter);
    final ArrayList<double[]> alphaList = new ArrayList<>(iter);
    final ArrayList<double[]> betaList = new ArrayList<>(iter);
    final ArrayList<double[]> xList = new ArrayList<>(iter);
    // Manipulate the background
    final double defaultBackground = background;
    try {
        background = 1e-2;
        createData(RngUtils.create(seed.getSeed()), 1, iter, paramsList, yList, true);
        final EjmlLinearSolver solver = new EjmlLinearSolver(1e-5, 1e-6);
        for (int i = 0; i < paramsList.size(); i++) {
            final double[] y = yList.get(i);
            final double[] a = paramsList.get(i);
            final BaseLsqLvmGradientProcedure p = LsqLvmGradientProcedureUtils.create(y, func);
            p.gradient(a);
            final double[] beta = p.beta;
            alphaList.add(p.getAlphaLinear());
            betaList.add(beta.clone());
            for (int j = 0; j < nparams; j++) {
                if (Math.abs(beta[j]) < 1e-6) {
                    logger.log(TestLogUtils.getRecord(Level.INFO, "[%d] Tiny beta %s %g", i, func.getGradientParameterName(j), beta[j]));
                }
            }
            // Solve
            if (!solver.solve(p.getAlphaMatrix(), beta)) {
                throw new AssertionError();
            }
            xList.add(beta);
        // System.out.println(Arrays.toString(beta));
        }
        // for (int b = 1; b < 1000; b *= 2)
        for (final double b : new double[] { -500, -100, -10, -1, -0.1, 0, 0.1, 1, 10, 100, 500 }) {
            final Statistics[] rel = new Statistics[nparams];
            final Statistics[] abs = new Statistics[nparams];
            if (debug) {
                for (int i = 0; i < nparams; i++) {
                    rel[i] = new Statistics();
                    abs[i] = new Statistics();
                }
            }
            for (int i = 0; i < paramsList.size(); i++) {
                final double[] y = add(yList.get(i), b);
                final double[] a = paramsList.get(i).clone();
                a[0] += b;
                final BaseLsqLvmGradientProcedure p = LsqLvmGradientProcedureUtils.create(y, func);
                p.gradient(a);
                final double[] beta = p.beta;
                final double[] alpha2 = alphaList.get(i);
                final double[] beta2 = betaList.get(i);
                final double[] x2 = xList.get(i);
                Assertions.assertArrayEquals(beta2, beta, 1e-10, "Beta");
                Assertions.assertArrayEquals(alpha2, p.getAlphaLinear(), 1e-10, "Alpha");
                // Solve
                solver.solve(p.getAlphaMatrix(), beta);
                Assertions.assertArrayEquals(x2, beta, 1e-10, "X");
                if (debug) {
                    for (int j = 0; j < nparams; j++) {
                        rel[j].add(DoubleEquality.relativeError(x2[j], beta[j]));
                        abs[j].add(Math.abs(x2[j] - beta[j]));
                    }
                }
            }
            if (debug) {
                for (int i = 0; i < nparams; i++) {
                    logger.log(TestLogUtils.getRecord(logLevel, "Bias = %.2f : %s : Rel %g +/- %g: Abs %g +/- %g", b, func.getGradientParameterName(i), rel[i].getMean(), rel[i].getStandardDeviation(), abs[i].getMean(), abs[i].getStandardDeviation()));
                }
            }
        }
    } finally {
        background = defaultBackground;
    }
}

Example 2 with ErfGaussian2DFunction

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

the class BaseFunctionSolverTest method fitAndComputeDeviationsMatch.

/**
 * Check the fit and compute deviations match. The first solver will be used to do the fit. This
 * is initialised from the solution so the convergence criteria can be set to accept the first
 * step. The second solver is used to compute deviations (thus is not initialised for fitting).
 *
 * @param seed the seed
 * @param solver1 the solver 1
 * @param solver2 the solver 2
 * @param noiseModel the noise model
 * @param useWeights the use weights
 */
void fitAndComputeDeviationsMatch(RandomSeed seed, BaseFunctionSolver solver1, BaseFunctionSolver solver2, NoiseModel noiseModel, boolean useWeights) {
    final double[] noise = getNoise(seed, noiseModel);
    if (solver1.isWeighted() && useWeights) {
        solver1.setWeights(getWeights(seed, noiseModel));
        solver2.setWeights(getWeights(seed, noiseModel));
    }
    // Draw target data
    final UniformRandomProvider rg = RngUtils.create(seed.getSeed());
    final double[] data = drawGaussian(p12, noise, noiseModel, rg);
    // fit with 2 peaks using the known params.
    // compare to 2 peak deviation computation.
    final Gaussian2DFunction f2 = GaussianFunctionFactory.create2D(2, size, size, flags, null);
    solver1.setGradientFunction(f2);
    solver2.setGradientFunction(f2);
    double[] params = p12.clone();
    double[] expected = new double[params.length];
    double[] observed = new double[params.length];
    solver1.fit(data, null, params, expected);
    // System.out.TestLog.fine(logger,"a="+Arrays.toString(a));
    solver2.computeDeviations(data, params, observed);
    // System.out.TestLog.fine(logger,"e2="+Arrays.toString(e));
    // System.out.TestLog.fine(logger,"o2="+Arrays.toString(o));
    Assertions.assertArrayEquals(observed, expected, "Fit 2 peaks and deviations 2 peaks do not match");
    // Try again with y-fit values
    params = p12.clone();
    final double[] o1 = new double[f2.size()];
    final double[] o2 = new double[o1.length];
    solver1.fit(data, o1, params, expected);
    // System.out.TestLog.fine(logger,"a="+Arrays.toString(a));
    solver2.computeValue(data, o2, params);
    Assertions.assertArrayEquals(observed, expected, "Fit 2 peaks with yFit and deviations 2 peaks do not match");
    final StandardValueProcedure p = new StandardValueProcedure();
    double[] ev = p.getValues(f2, params);
    Assertions.assertArrayEquals(ev, o1, 1e-8, "Fit 2 peaks yFit");
    Assertions.assertArrayEquals(ev, o2, 1e-8, "computeValue 2 peaks yFit");
    if (solver1 instanceof SteppingFunctionSolver) {
        // fit with 1 peak + 1 precomputed using the known params.
        // compare to 2 peak deviation computation.
        final ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, size, size, flags, null);
        final Gradient2Function pf1 = OffsetGradient2Function.wrapGradient2Function(f1, p2v);
        solver1.setGradientFunction(pf1);
        params = p1.clone();
        expected = new double[params.length];
        solver1.fit(data, null, params, expected);
        // To copy the second peak
        final double[] a2 = p12.clone();
        // Add the same fitted first peak
        System.arraycopy(params, 0, a2, 0, params.length);
        solver2.computeDeviations(data, a2, observed);
        // System.out.TestLog.fine(logger,"e1p1=" + Arrays.toString(e));
        // System.out.TestLog.fine(logger,"o2=" + Arrays.toString(o));
        // Deviation should be lower with only 1 peak.
        // Due to matrix inversion this may not be the case for all parameters so count.
        int ok = 0;
        int fail = 0;
        final StringBuilder sb = new StringBuilder();
        for (int i = 0; i < expected.length; i++) {
            if (expected[i] <= observed[i]) {
                ok++;
                continue;
            }
            fail++;
            TextUtils.formatTo(sb, "Fit 1 peak + 1 precomputed is higher than deviations 2 peaks %s: %s > %s", Gaussian2DFunction.getName(i), expected[i], observed[i]);
        }
        if (fail > ok) {
            Assertions.fail(sb.toString());
        }
        // Try again with y-fit values
        params = p1.clone();
        Arrays.fill(o1, 0);
        Arrays.fill(o2, 0);
        observed = new double[params.length];
        solver1.fit(data, o1, params, observed);
        solver2.computeValue(data, o2, a2);
        Assertions.assertArrayEquals(observed, expected, 1e-8, "Fit 1 peak + 1 precomputed with yFit and deviations 1 peak + " + "1 precomputed do not match");
        ev = p.getValues(pf1, params);
        Assertions.assertArrayEquals(ev, o1, 1e-8, "Fit 1 peak + 1 precomputed yFit");
        Assertions.assertArrayEquals(ev, o2, 1e-8, "computeValue 1 peak + 1 precomputed yFit");
    }
}

Example 3 with ErfGaussian2DFunction

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

the class BaseFunctionSolverTest method fitAndComputeValueMatch.

/**
 * Check the fit and compute values match. The first solver will be used to do the fit. This is
 * initialised from the solution so the convergence criteria can be set to accept the first step.
 * The second solver is used to compute values (thus is not initialised for fitting).
 *
 * @param seed the seed
 * @param solver1 the solver
 * @param solver2 the solver 2
 * @param noiseModel the noise model
 * @param useWeights the use weights
 */
void fitAndComputeValueMatch(RandomSeed seed, BaseFunctionSolver solver1, BaseFunctionSolver solver2, NoiseModel noiseModel, boolean useWeights) {
    final double[] noise = getNoise(seed, noiseModel);
    if (solver1.isWeighted() && useWeights) {
        solver1.setWeights(getWeights(seed, noiseModel));
        solver2.setWeights(getWeights(seed, noiseModel));
    }
    // Draw target data
    final UniformRandomProvider rg = RngUtils.create(seed.getSeed());
    final double[] data = drawGaussian(p12, noise, noiseModel, rg);
    // fit with 2 peaks using the known params.
    final Gaussian2DFunction f2 = GaussianFunctionFactory.create2D(2, size, size, flags, null);
    solver1.setGradientFunction(f2);
    solver2.setGradientFunction(f2);
    double[] params = p12.clone();
    solver1.fit(data, null, params, null);
    solver2.computeValue(data, null, params);
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-10, 0);
    double v1 = solver1.getValue();
    double v2 = solver2.getValue();
    TestAssertions.assertTest(v1, v2, predicate, "Fit 2 peaks and computeValue");
    final double[] o1 = new double[f2.size()];
    final double[] o2 = new double[o1.length];
    solver1.fit(data, o1, params, null);
    solver2.computeValue(data, o2, params);
    v1 = solver1.getValue();
    v2 = solver2.getValue();
    TestAssertions.assertTest(v1, v2, predicate, "Fit 2 peaks and computeValue with yFit");
    final StandardValueProcedure p = new StandardValueProcedure();
    double[] expected = p.getValues(f2, params);
    Assertions.assertArrayEquals(expected, o1, 1e-8, "Fit 2 peaks yFit");
    Assertions.assertArrayEquals(expected, o2, 1e-8, "computeValue 2 peaks yFit");
    if (solver1 instanceof SteppingFunctionSolver) {
        // fit with 1 peak + 1 precomputed using the known params.
        // compare to 2 peak computation.
        final ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, size, size, flags, null);
        final Gradient2Function pf1 = OffsetGradient2Function.wrapGradient2Function(f1, p2v);
        solver1.setGradientFunction(pf1);
        solver2.setGradientFunction(pf1);
        params = p1.clone();
        solver1.fit(data, null, params, null);
        solver2.computeValue(data, null, params);
        v1 = solver1.getValue();
        v2 = solver2.getValue();
        TestAssertions.assertTest(v1, v2, predicate, "Fit 1 peak + 1 precomputed and computeValue");
        Arrays.fill(o1, 0);
        Arrays.fill(o2, 0);
        solver1.fit(data, o1, params, null);
        solver2.computeValue(data, o2, params);
        v1 = solver1.getValue();
        v2 = solver2.getValue();
        TestAssertions.assertTest(v1, v2, predicate, "Fit 1 peak + 1 precomputed and computeValue with yFit");
        expected = p.getValues(pf1, params);
        Assertions.assertArrayEquals(expected, o1, 1e-8, "Fit 1 peak + 1 precomputed yFit");
        Assertions.assertArrayEquals(expected, o2, 1e-8, "computeValue 1 peak + 1 precomputed yFit");
    }
}

Example 4 with ErfGaussian2DFunction

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

the class PoissonGradientProcedureTest method crlbIsHigherWithPrecomputed.

@SeededTest
void crlbIsHigherWithPrecomputed(RandomSeed seed) {
    final int iter = 10;
    final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
    final ErfGaussian2DFunction func = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final double[] a = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    final int n = func.getNumberOfGradients();
    // Get a background
    final double[] b = new double[func.size()];
    for (int i = 0; i < b.length; i++) {
        b[i] = nextUniform(rng, 1, 2);
    }
    for (int i = 0; i < iter; i++) {
        a[Gaussian2DFunction.BACKGROUND] = nextUniform(rng, 0.1, 0.3);
        a[Gaussian2DFunction.SIGNAL] = nextUniform(rng, 100, 300);
        a[Gaussian2DFunction.X_POSITION] = nextUniform(rng, 4, 6);
        a[Gaussian2DFunction.Y_POSITION] = nextUniform(rng, 4, 6);
        a[Gaussian2DFunction.X_SD] = nextUniform(rng, 1, 1.3);
        a[Gaussian2DFunction.Y_SD] = nextUniform(rng, 1, 1.3);
        final PoissonGradientProcedure p1 = PoissonGradientProcedureUtils.create(func);
        p1.computeFisherInformation(a);
        final PoissonGradientProcedure p2 = PoissonGradientProcedureUtils.create(OffsetGradient1Function.wrapGradient1Function(func, b));
        p2.computeFisherInformation(a);
        final FisherInformationMatrix m1 = new FisherInformationMatrix(p1.getLinear(), n);
        final FisherInformationMatrix m2 = new FisherInformationMatrix(p2.getLinear(), n);
        final double[] crlb1 = m1.crlb();
        final double[] crlb2 = m2.crlb();
        Assertions.assertNotNull(crlb1);
        Assertions.assertNotNull(crlb2);
        // Arrays.toString(crlb2));
        for (int j = 0; j < n; j++) {
            Assertions.assertTrue(crlb1[j] < crlb2[j]);
        }
    }
}

Example 5 with ErfGaussian2DFunction

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

the class FastMleGradient2ProcedureTest method gradientProcedureComputesSameWithPrecomputed.

@SeededTest
void gradientProcedureComputesSameWithPrecomputed(RandomSeed seed) {
    final int iter = 10;
    final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
    final ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final ErfGaussian2DFunction f2 = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(2, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final double[] a1 = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    final double[] a2 = new double[1 + 2 * Gaussian2DFunction.PARAMETERS_PER_PEAK];
    final double[] x = new double[f1.size()];
    final double[] b = new double[f1.size()];
    for (int i = 0; i < iter; i++) {
        final int ii = i;
        a2[Gaussian2DFunction.BACKGROUND] = nextUniform(rng, 0.1, 0.3);
        a2[Gaussian2DFunction.SIGNAL] = nextUniform(rng, 100, 300);
        a2[Gaussian2DFunction.X_POSITION] = nextUniform(rng, 3, 5);
        a2[Gaussian2DFunction.Y_POSITION] = nextUniform(rng, 3, 5);
        a2[Gaussian2DFunction.Z_POSITION] = nextUniform(rng, -2, 2);
        a2[Gaussian2DFunction.X_SD] = nextUniform(rng, 1, 1.3);
        a2[Gaussian2DFunction.Y_SD] = nextUniform(rng, 1, 1.3);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.SIGNAL] = nextUniform(rng, 100, 300);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.X_POSITION] = nextUniform(rng, 5, 7);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Y_POSITION] = nextUniform(rng, 5, 7);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Z_POSITION] = nextUniform(rng, -3, 1);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.X_SD] = nextUniform(rng, 1, 1.3);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Y_SD] = nextUniform(rng, 1, 1.3);
        // Simulate Poisson data
        f2.initialise0(a2);
        f1.forEach(new ValueProcedure() {

            int index = 0;

            @Override
            public void execute(double value) {
                if (value > 0) {
                    x[index++] = GdscSmlmTestUtils.createPoissonSampler(rng, value).sample();
                } else {
                    x[index++] = 0;
                }
            }
        });
        // Precompute peak 2 (no background)
        a1[Gaussian2DFunction.BACKGROUND] = 0;
        for (int j = 1; j < 7; j++) {
            a1[j] = a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + j];
        }
        f1.initialise0(a1);
        f1.forEach(new ValueProcedure() {

            int index = 0;

            @Override
            public void execute(double value) {
                b[index++] = value;
            }
        });
        // Reset to peak 1
        for (int j = 0; j < 7; j++) {
            a1[j] = a2[j];
        }
        // Compute peak 1+2
        final FastMleGradient2Procedure p12 = FastMleGradient2ProcedureUtils.create(x, f2);
        p12.computeSecondDerivative(a2);
        final double[] d11 = Arrays.copyOf(p12.d1, f1.getNumberOfGradients());
        final double[] d21 = Arrays.copyOf(p12.d2, f1.getNumberOfGradients());
        // Compute peak 1+(precomputed 2)
        final FastMleGradient2Procedure p1b2 = FastMleGradient2ProcedureUtils.create(x, OffsetGradient2Function.wrapGradient2Function(f1, b));
        p1b2.computeSecondDerivative(a1);
        final double[] d12 = p1b2.d1;
        final double[] d22 = p1b2.d2;
        Assertions.assertArrayEquals(p12.u, p1b2.u, 1e-10, () -> " Result: Not same @ " + ii);
        Assertions.assertArrayEquals(d11, d12, 1e-10, () -> " D1: Not same @ " + ii);
        Assertions.assertArrayEquals(d21, d22, 1e-10, () -> " D2: Not same @ " + ii);
        final double[] v1 = p12.computeValue(a2);
        final double[] v2 = p1b2.computeValue(a1);
        Assertions.assertArrayEquals(v1, v2, 1e-10, () -> " Value: Not same @ " + ii);
        final double[] d1 = Arrays.copyOf(p12.computeFirstDerivative(a2), f1.getNumberOfGradients());
        final double[] d2 = p1b2.computeFirstDerivative(a1);
        Assertions.assertArrayEquals(d1, d2, 1e-10, () -> " 1st derivative: Not same @ " + ii);
    }
}