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

Example 6 with Gaussian2DFunction

use of uk.ac.sussex.gdsc.smlm.function.gaussian.Gaussian2DFunction 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 7 with Gaussian2DFunction

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

the class BoundedFunctionSolverTest method getLvm.

private static NonLinearFit getLvm(int bounded, int clamping, boolean mle) {
    final Gaussian2DFunction f = GaussianFunctionFactory.create2D(1, size, size, flags, null);
    final StoppingCriteria sc = new ErrorStoppingCriteria();
    sc.setMaximumIterations(100);
    final NonLinearFit solver = (bounded != 0 || clamping != 0) ? new BoundedNonLinearFit(f, sc, null) : new NonLinearFit(f, sc);
    if (clamping != 0) {
        final BoundedNonLinearFit bsolver = (BoundedNonLinearFit) solver;
        final ParameterBounds bounds = new ParameterBounds(f);
        bounds.setClampValues(defaultClampValues);
        bounds.setDynamicClamp(clamping == 2);
        bsolver.setBounds(bounds);
    }
    solver.setMle(mle);
    solver.setInitialLambda(1);
    return solver;
}

Example 8 with Gaussian2DFunction

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

the class FisherInformationMatrixTest method computeWithSubsetReducesTheCrlb.

@SeededTest
void computeWithSubsetReducesTheCrlb(RandomSeed seed) {
    final UniformRandomProvider rg = RngUtils.create(seed.getSeed());
    final Gaussian2DFunction f = createFunction(10, 1);
    final int perPeak = f.getGradientParametersPerPeak();
    // Create a matrix with 2 peaks + background
    final FisherInformationMatrix m = createFisherInformationMatrix(rg, 1 + 2 * perPeak, 0);
    // Subset each peak
    final int[] indices = SimpleArrayUtils.natural(1 + perPeak);
    final FisherInformationMatrix m1 = m.subset(indices);
    for (int i = 1; i < indices.length; i++) {
        indices[i] += perPeak;
    }
    final FisherInformationMatrix m2 = m.subset(indices);
    // TestLog.fine(logger,m.getMatrix());
    // TestLog.fine(logger,m1.getMatrix());
    // TestLog.fine(logger,m2.getMatrix());
    final double[] crlb = m.crlb();
    final double[] crlb1 = m1.crlb();
    final double[] crlb2 = m2.crlb();
    final double[] crlbB = Arrays.copyOf(crlb1, crlb.length);
    System.arraycopy(crlb2, 1, crlbB, crlb1.length, perPeak);
    // Removing the interaction between fit parameters lowers the bounds
    for (int i = 0; i < crlb.length; i++) {
        Assertions.assertTrue(crlbB[i] < crlb[i]);
    }
}

Example 9 with Gaussian2DFunction

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

the class UnivariateLikelihoodFisherInformationCalculatorTest method canComputePerPixelPoissonGaussianApproximationFisherInformation.

private static void canComputePerPixelPoissonGaussianApproximationFisherInformation(UniformRandomProvider rng) {
    // Create function
    final Gaussian2DFunction func = GaussianFunctionFactory.create2D(1, 10, 10, GaussianFunctionFactory.FIT_ERF_CIRCLE, null);
    final double[] params = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    params[Gaussian2DFunction.BACKGROUND] = nextUniform(rng, 0.1, 0.3);
    params[Gaussian2DFunction.SIGNAL] = nextUniform(rng, 100, 300);
    params[Gaussian2DFunction.X_POSITION] = nextUniform(rng, 4, 6);
    params[Gaussian2DFunction.Y_POSITION] = nextUniform(rng, 4, 6);
    params[Gaussian2DFunction.X_SD] = nextUniform(rng, 1, 1.3);
    Gradient1Function f1 = func;
    FisherInformation[] fi;
    // Get a per-pixel variance
    final double[] var = new double[func.size()];
    fi = new FisherInformation[var.length];
    for (int i = var.length; i-- > 0; ) {
        var[i] = 0.9 + 0.2 * rng.nextDouble();
        fi[i] = new PoissonGaussianApproximationFisherInformation(Math.sqrt(var[i]));
    }
    f1 = (Gradient1Function) OffsetFunctionFactory.wrapFunction(func, var);
    // This introduces a dependency on a different package, and relies on that
    // computing the correct answer. However that code predates this and so the
    // test ensures that the FisherInformationCalculator functions correctly.
    final PoissonGradientProcedure p1 = PoissonGradientProcedureUtils.create(f1);
    p1.computeFisherInformation(params);
    final double[] e = p1.getLinear();
    final FisherInformationCalculator calc = new UnivariateLikelihoodFisherInformationCalculator(func, fi);
    final FisherInformationMatrix I = calc.compute(params);
    final double[] o = I.getMatrix().data;
    TestAssertions.assertArrayTest(e, o, TestHelper.doublesAreClose(1e-6, 0));
}

Example 10 with Gaussian2DFunction

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

the class GradientCalculatorSpeedTest method gradientCalculatorComputesSameOutputWithBias.

@SeededTest
void gradientCalculatorComputesSameOutputWithBias(RandomSeed seed) {
    final Gaussian2DFunction func = new SingleEllipticalGaussian2DFunction(blockWidth, blockWidth);
    final int nparams = func.getNumberOfGradients();
    final GradientCalculator calc = new GradientCalculator(nparams);
    final int n = func.size();
    final int iter = 50;
    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;
    final boolean report = logger.isLoggable(Level.INFO);
    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 double[][] alpha = new double[nparams][nparams];
            final double[] beta = new double[nparams];
            calc.findLinearised(n, y, a, alpha, beta, func);
            alphaList.add(alpha);
            betaList.add(beta.clone());
            for (int j = 0; j < nparams; j++) {
                if (Math.abs(beta[j]) < 1e-6) {
                    logger.info(FunctionUtils.getSupplier("[%d] Tiny beta %s %g", i, func.getGradientParameterName(j), beta[j]));
                }
            }
            // Solve
            if (!solver.solve(alpha, beta)) {
                throw new AssertionError();
            }
            xList.add(beta);
        // System.out.println(Arrays.toString(beta));
        }
        final double[][] alpha = new double[nparams][nparams];
        final double[] beta = new double[nparams];
        final Statistics[] rel = new Statistics[nparams];
        final Statistics[] abs = new Statistics[nparams];
        for (int i = 0; i < nparams; i++) {
            rel[i] = new Statistics();
            abs[i] = new Statistics();
        }
        final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-10, 0);
        // for (double b : new double[] { -500, -100, -10, -1, -0.1, 0.1, 1, 10, 100, 500 })
        for (final double b : new double[] { -10, -1, -0.1, 0.1, 1, 10 }) {
            if (report) {
                for (int i = 0; i < nparams; i++) {
                    rel[i].reset();
                    abs[i].reset();
                }
            }
            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;
                calc.findLinearised(n, y, a, alpha, beta, func);
                final double[][] alpha2 = alphaList.get(i);
                final double[] beta2 = betaList.get(i);
                final double[] x2 = xList.get(i);
                TestAssertions.assertArrayTest(beta2, beta, predicate, "Beta");
                TestAssertions.assertArrayTest(alpha2, alpha, predicate, "Alpha");
                // Solve
                solver.solve(alpha, beta);
                Assertions.assertArrayEquals(x2, beta, 1e-10, "X");
                if (report) {
                    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 (report) {
                for (int i = 0; i < nparams; i++) {
                    logger.info(FunctionUtils.getSupplier("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;
    }
}