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

Example 1 with Gaussian2DFunction

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

the class DoubleDht2DTest method createData.

private static DoubleDht2D createData(double cx, double cy) {
    final Gaussian2DFunction f = GaussianFunctionFactory.create2D(1, size, size, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final double[] a = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    a[Gaussian2DFunction.SIGNAL] = 1;
    a[Gaussian2DFunction.X_POSITION] = cx;
    a[Gaussian2DFunction.Y_POSITION] = cy;
    a[Gaussian2DFunction.X_SD] = 1.2;
    a[Gaussian2DFunction.Y_SD] = 1.1;
    final StandardValueProcedure p = new StandardValueProcedure();
    p.getValues(f, a);
    return new DoubleDht2D(size, size, p.values, false);
}

Example 2 with Gaussian2DFunction

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

the class FloatDht2DTest method createData.

private static FloatDht2D createData(double cx, double cy) {
    final Gaussian2DFunction f = GaussianFunctionFactory.create2D(1, size, size, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final double[] a = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    a[Gaussian2DFunction.SIGNAL] = 1;
    a[Gaussian2DFunction.X_POSITION] = cx;
    a[Gaussian2DFunction.Y_POSITION] = cy;
    a[Gaussian2DFunction.X_SD] = 1.2;
    a[Gaussian2DFunction.Y_SD] = 1.1;
    final StandardFloatValueProcedure p = new StandardFloatValueProcedure();
    p.getValues(f, a);
    return new FloatDht2D(size, size, p.values, false);
}

Example 3 with Gaussian2DFunction

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

the class LvmGradientProcedureTest method gradientProcedureSupportsPrecomputed.

private void gradientProcedureSupportsPrecomputed(RandomSeed seed, final Type type, boolean checkGradients) {
    final int iter = 10;
    final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
    final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(rng, 0, noise);
    final ArrayList<double[]> paramsList = new ArrayList<>(iter);
    final ArrayList<double[]> yList = new ArrayList<>(iter);
    // 3 peaks
    createData(rng, 3, iter, paramsList, yList, true);
    for (int i = 0; i < paramsList.size(); i++) {
        final double[] y = yList.get(i);
        // Add Gaussian read noise so we have negatives
        final double min = MathUtils.min(y);
        for (int j = 0; j < y.length; j++) {
            y[j] = y[i] - min + gs.sample();
        }
    }
    // We want to know that:
    // y|peak1+peak2+peak3 == y|peak1+peak2+peak3(precomputed)
    // We want to know when:
    // y|peak1+peak2+peak3 != y-peak3|peak1+peak2
    // i.e. we cannot subtract a precomputed peak from the data, it must be included in the fit
    // E.G. LSQ - subtraction is OK, MLE/WLSQ - subtraction is not allowed
    final Gaussian2DFunction f123 = GaussianFunctionFactory.create2D(3, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final Gaussian2DFunction f12 = GaussianFunctionFactory.create2D(2, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final Gaussian2DFunction f3 = GaussianFunctionFactory.create2D(1, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final FastLog fastLog = type == Type.FAST_LOG_MLE ? getFastLog() : null;
    final int nparams = f12.getNumberOfGradients();
    final int[] indices = f12.gradientIndices();
    final double[] b = new double[f12.size()];
    // for checking strict equivalence
    final DoubleEquality eq = new DoubleEquality(1e-8, 1e-16);
    // for the gradients
    final double delta = 1e-4;
    final DoubleEquality eq2 = new DoubleEquality(5e-2, 1e-16);
    final double[] a1peaks = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    final double[] y_b = new double[b.length];
    // Count the number of failures for each gradient
    final int failureLimit = TestCounter.computeFailureLimit(iter, 0.1);
    final TestCounter failCounter = new TestCounter(failureLimit, nparams * 2);
    for (int i = 0; i < paramsList.size(); i++) {
        final int ii = i;
        final double[] y = yList.get(i);
        final double[] a3peaks = paramsList.get(i);
        // logger.fine(FunctionUtils.getSupplier("[%d] a=%s", i, Arrays.toString(a3peaks));
        final double[] a2peaks = Arrays.copyOf(a3peaks, 1 + 2 * Gaussian2DFunction.PARAMETERS_PER_PEAK);
        final double[] a2peaks2 = a2peaks.clone();
        for (int j = 1; j < a1peaks.length; j++) {
            a1peaks[j] = a3peaks[j + 2 * Gaussian2DFunction.PARAMETERS_PER_PEAK];
        }
        // Evaluate peak 3 to get the background and subtract it from the data to get the new data
        f3.initialise0(a1peaks);
        f3.forEach(new ValueProcedure() {

            int index = 0;

            @Override
            public void execute(double value) {
                b[index] = value;
                // Remove negatives for MLE
                if (type.isMle()) {
                    y[index] = Math.max(0, y[index]);
                    y_b[index] = Math.max(0, y[index] - value);
                } else {
                    y_b[index] = y[index] - value;
                }
                index++;
            }
        });
        final LvmGradientProcedure p123 = LvmGradientProcedureUtils.create(y, f123, type, fastLog);
        // ///////////////////////////////////
        // These should be the same
        // ///////////////////////////////////
        final LvmGradientProcedure p12b3 = LvmGradientProcedureUtils.create(y, OffsetGradient1Function.wrapGradient1Function(f12, b), type, fastLog);
        // Check they are the same
        p123.gradient(a3peaks);
        final double[][] m123 = p123.getAlphaMatrix();
        p12b3.gradient(a2peaks);
        double value = p12b3.value;
        final double[] beta = p12b3.beta.clone();
        double[][] alpha = p12b3.getAlphaMatrix();
        if (!eq.almostEqualRelativeOrAbsolute(p123.value, value)) {
            Assertions.fail(FunctionUtils.getSupplier("p12b3 Not same value @ %d (error=%s) : %s == %s", i, DoubleEquality.relativeError(p123.value, value), p123.value, value));
        }
        if (!almostEqualRelativeOrAbsolute(eq, beta, p123.beta)) {
            Assertions.fail(FunctionUtils.getSupplier("p12b3 Not same gradient @ %d (error=%s) : %s vs %s", i, relativeError(beta, p123.beta), Arrays.toString(beta), Arrays.toString(p123.beta)));
        }
        for (int j = 0; j < alpha.length; j++) {
            // Arrays.toString(m123[j]));
            if (!almostEqualRelativeOrAbsolute(eq, alpha[j], m123[j])) {
                Assertions.fail(FunctionUtils.getSupplier("p12b3 Not same alpha @ %d,%d (error=%s) : %s vs %s", i, j, relativeError(alpha[j], m123[j]), Arrays.toString(alpha[j]), Arrays.toString(m123[j])));
            }
        }
        // Check actual gradients are correct
        if (checkGradients) {
            for (int j = 0; j < nparams; j++) {
                final int jj = j;
                final int k = indices[j];
                // double d = Precision.representableDelta(a2peaks[k], (a2peaks[k] == 0) ? 1e-3 :
                // a2peaks[k] * delta);
                final double d = Precision.representableDelta(a2peaks[k], delta);
                a2peaks2[k] = a2peaks[k] + d;
                p12b3.value(a2peaks2);
                final double s1 = p12b3.value;
                a2peaks2[k] = a2peaks[k] - d;
                p12b3.value(a2peaks2);
                final double s2 = p12b3.value;
                a2peaks2[k] = a2peaks[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 (%f)", i, k, s,
                // Gaussian2DFunction.getName(k), a2peaks[k], d, beta[j], gradient,
                // DoubleEquality.relativeError(gradient, beta[j]));
                failCounter.run(j, () -> eq2.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)));
                });
            }
        }
        // ///////////////////////////////////
        // This may be different
        // ///////////////////////////////////
        final LvmGradientProcedure p12m3 = LvmGradientProcedureUtils.create(y_b, f12, type, fastLog);
        // Check these may be different.
        // Sometimes they are not different.
        p12m3.gradient(a2peaks);
        value = p12m3.value;
        System.arraycopy(p12m3.beta, 0, beta, 0, p12m3.beta.length);
        alpha = p12m3.getAlphaMatrix();
        if (type != Type.LSQ) {
            if (eq.almostEqualRelativeOrAbsolute(p123.value, value)) {
                logger.log(TestLogUtils.getFailRecord("p12b3 Same value @ %d (error=%s) : %s == %s", i, DoubleEquality.relativeError(p123.value, value), p123.value, value));
            }
            if (almostEqualRelativeOrAbsolute(eq, beta, p123.beta)) {
                logger.log(TestLogUtils.getFailRecord("p12b3 Same gradient @ %d (error=%s) : %s vs %s", i, relativeError(beta, p123.beta), Arrays.toString(beta), Arrays.toString(p123.beta)));
            }
            // Note: Test the matrix is different by finding 1 different column
            int dj = -1;
            for (int j = 0; j < alpha.length; j++) {
                // Arrays.toString(m123[j]));
                if (!almostEqualRelativeOrAbsolute(eq, alpha[j], m123[j])) {
                    // Different column
                    dj = j;
                    break;
                }
            }
            if (dj == -1) {
                // Find biggest error for reporting. This helps set the test tolerance.
                double error = 0;
                dj = -1;
                for (int j = 0; j < alpha.length; j++) {
                    final double e = relativeError(alpha[j], m123[j]);
                    if (error <= e) {
                        error = e;
                        dj = j;
                    }
                }
                logger.log(TestLogUtils.getFailRecord("p12b3 Same alpha @ %d,%d (error=%s) : %s vs %s", i, dj, error, Arrays.toString(alpha[dj]), Arrays.toString(m123[dj])));
            }
        } else {
            if (!eq.almostEqualRelativeOrAbsolute(p123.value, value)) {
                logger.log(TestLogUtils.getFailRecord("p12b3 Not same value @ %d (error=%s) : %s == %s", i, DoubleEquality.relativeError(p123.value, value), p123.value, value));
            }
            if (!almostEqualRelativeOrAbsolute(eq, beta, p123.beta)) {
                logger.log(TestLogUtils.getFailRecord("p12b3 Not same gradient @ %d (error=%s) : %s vs %s", i, relativeError(beta, p123.beta), Arrays.toString(beta), Arrays.toString(p123.beta)));
            }
            for (int j = 0; j < alpha.length; j++) {
                // Arrays.toString(m123[j]));
                if (!almostEqualRelativeOrAbsolute(eq, alpha[j], m123[j])) {
                    logger.log(TestLogUtils.getFailRecord("p12b3 Not same alpha @ %d,%d (error=%s) : %s vs %s", i, j, relativeError(alpha[j], m123[j]), Arrays.toString(alpha[j]), Arrays.toString(m123[j])));
                }
            }
        }
        // Check actual gradients are correct
        if (!checkGradients) {
            continue;
        }
        for (int j = 0; j < nparams; j++) {
            final int jj = j;
            final int k = indices[j];
            // double d = Precision.representableDelta(a2peaks[k], (a2peaks[k] == 0) ? 1e-3 : a2peaks[k]
            // * delta);
            final double d = Precision.representableDelta(a2peaks[k], delta);
            a2peaks2[k] = a2peaks[k] + d;
            p12m3.value(a2peaks2);
            final double s1 = p12m3.value;
            a2peaks2[k] = a2peaks[k] - d;
            p12m3.value(a2peaks2);
            final double s2 = p12m3.value;
            a2peaks2[k] = a2peaks[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 (%f)", i, k, s,
            // Gaussian2DFunction.getName(k), a2peaks[k], d, beta[j], gradient,
            // DoubleEquality.relativeError(gradient, beta[j]));
            failCounter.run(nparams + j, () -> eq2.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)));
            });
        }
    }
}

Example 4 with Gaussian2DFunction

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

the class BaseFunctionSolverTest method drawGaussian.

/**
 * Draw a Gaussian with Poisson shot noise and Gaussian read noise.
 *
 * @param params The Gaussian parameters
 * @param noise The read noise
 * @param noiseModel the noise model
 * @param rg the random generator
 * @return The data
 */
static double[] drawGaussian(double[] params, double[] noise, NoiseModel noiseModel, UniformRandomProvider rg) {
    final int n = params.length / Gaussian2DFunction.PARAMETERS_PER_PEAK;
    final Gaussian2DFunction f = GaussianFunctionFactory.create2D(n, size, size, flags, null);
    final double[] data = f.computeValues(params);
    // Poisson noise
    for (int i = 0; i < data.length; i++) {
        if (data[i] > 0) {
            data[i] = GdscSmlmTestUtils.createPoissonSampler(rg, data[i]).sample();
        }
    }
    // Simulate EM-gain
    if (noiseModel == NoiseModel.EMCCD) {
        // Use a gamma distribution
        // Since the call random.nextGamma(...) creates a Gamma distribution
        // which pre-calculates factors only using the scale parameter we
        // create a custom gamma distribution where the shape can be set as a property.
        final MarsagliaTsangGammaSampler gd = new MarsagliaTsangGammaSampler(rg, 1, emGain);
        for (int i = 0; i < data.length; i++) {
            if (data[i] > 0) {
                gd.setAlpha(data[i]);
                // The sample will amplify the signal so we remap to the original scale
                data[i] = gd.sample() / emGain;
            }
        }
    }
    // Read-noise
    final NormalizedGaussianSampler gs = SamplerUtils.createNormalizedGaussianSampler(rg);
    if (noise != null) {
        for (int i = 0; i < data.length; i++) {
            data[i] += gs.sample() * noise[i];
        }
    }
    // uk.ac.sussex.gdsc.core.ij.Utils.display("Spot", data, size, size);
    return data;
}

Example 5 with Gaussian2DFunction

use of uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunction 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");
    }
}