Examples with PoissonGradientProcedure uk.ac.sussex.gdsc.smlm.fitting.nonlinear.gradient.PoissonGradientProcedure used on opensource projects

Example 1 with PoissonGradientProcedure

use of uk.ac.sussex.gdsc.smlm.fitting.nonlinear.gradient.PoissonGradientProcedure 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 2 with PoissonGradientProcedure

use of uk.ac.sussex.gdsc.smlm.fitting.nonlinear.gradient.PoissonGradientProcedure in project GDSC-SMLM by aherbert.

the class FastMleSteppingFunctionSolver method computeLastFisherInformationMatrix.

@Override
protected FisherInformationMatrix computeLastFisherInformationMatrix(double[] fx) {
    Gradient2Function f2 = (Gradient2Function) function;
    // Capture the y-values if necessary
    if (fx != null && fx.length == f2.size()) {
        f2 = new Gradient2FunctionValueStore(f2, fx);
    }
    // Add the weights if necessary
    if (obsVariances != null) {
        f2 = OffsetGradient2Function.wrapGradient2Function(f2, obsVariances);
    }
    // The fisher information is that for a Poisson process
    final PoissonGradientProcedure p = PoissonGradientProcedureUtils.create(f2);
    initialiseAndRun(p);
    if (p.isNaNGradients()) {
        throw new FunctionSolverException(FitStatus.INVALID_GRADIENTS);
    }
    return new FisherInformationMatrix(p.getLinear(), p.numberOfGradients);
}

Example 3 with PoissonGradientProcedure

use of uk.ac.sussex.gdsc.smlm.fitting.nonlinear.gradient.PoissonGradientProcedure in project GDSC-SMLM by aherbert.

the class FastMleSteppingFunctionSolver method computeFunctionFisherInformationMatrix.

@Override
protected FisherInformationMatrix computeFunctionFisherInformationMatrix(double[] y, double[] a) {
    // The fisher information is that for a Poisson process.
    // We must wrap the gradient function if weights are present.
    Gradient1Function f1 = (Gradient1Function) function;
    if (obsVariances != null) {
        f1 = OffsetGradient1Function.wrapGradient1Function(f1, obsVariances);
    }
    final PoissonGradientProcedure p = PoissonGradientProcedureUtils.create(f1);
    p.computeFisherInformation(a);
    if (p.isNaNGradients()) {
        throw new FunctionSolverException(FitStatus.INVALID_GRADIENTS);
    }
    return new FisherInformationMatrix(p.getLinear(), p.numberOfGradients);
}

Example 4 with PoissonGradientProcedure

use of uk.ac.sussex.gdsc.smlm.fitting.nonlinear.gradient.PoissonGradientProcedure in project GDSC-SMLM by aherbert.

the class BaseFunctionSolverTest method canFitSingleGaussian.

void canFitSingleGaussian(RandomSeed seed, FunctionSolver solver, boolean applyBounds, NoiseModel noiseModel) {
    // Allow reporting the fit deviations
    final boolean report = false;
    double[] crlb = null;
    SimpleArrayMoment moment = null;
    final double[] noise = getNoise(seed, noiseModel);
    if (solver.isWeighted()) {
        solver.setWeights(getWeights(seed, noiseModel));
    }
    final UniformRandomProvider rg = RngUtils.create(seed.getSeed());
    for (final double s : signal) {
        final double[] expected = createParams(1, s, 0, 0, 1);
        final double[] lower = createParams(0, s * 0.5, -0.3, -0.3, 0.8);
        final double[] upper = createParams(3, s * 2, 0.3, 0.3, 1.2);
        if (applyBounds) {
            solver.setBounds(lower, upper);
        }
        if (report) {
            // Compute the CRLB for a Poisson process
            final PoissonGradientProcedure gp = PoissonGradientProcedureUtils.create((Gradient1Function) ((BaseFunctionSolver) solver).getGradientFunction());
            gp.computeFisherInformation(expected);
            final FisherInformationMatrix f = new FisherInformationMatrix(gp.getLinear(), gp.numberOfGradients);
            crlb = f.crlbSqrt();
            // Compute the deviations.
            // Note this is not the same as the CRLB as the fit is repeated
            // against the same data.
            // It should be repeated against different data generated with constant
            // parameters and variable noise.
            moment = new SimpleArrayMoment();
        }
        final double[] data = drawGaussian(expected, noise, noiseModel, rg);
        for (final double db : base) {
            for (final double dx : shift) {
                for (final double dy : shift) {
                    for (final double dsx : factor) {
                        final double[] p = createParams(db, s, dx, dy, dsx);
                        final double[] fp = fitGaussian(solver, data, p, expected);
                        for (int i = 0; i < expected.length; i++) {
                            if (fp[i] < lower[i]) {
                                Assertions.fail(FunctionUtils.getSupplier("Fit Failed: [%d] %.2f < %.2f: %s != %s", i, fp[i], lower[i], Arrays.toString(fp), Arrays.toString(expected)));
                            }
                            if (fp[i] > upper[i]) {
                                Assertions.fail(FunctionUtils.getSupplier("Fit Failed: [%d] %.2f > %.2f: %s != %s", i, fp[i], upper[i], Arrays.toString(fp), Arrays.toString(expected)));
                            }
                            if (report) {
                                fp[i] = expected[i] - fp[i];
                            }
                        }
                        // Store the deviations
                        if (report) {
                            moment.add(fp);
                        }
                    }
                }
            }
        }
        // Report
        if (report) {
            logger.info(FunctionUtils.getSupplier("%s %s %f : CRLB = %s, Deviations = %s", solver.getClass().getSimpleName(), noiseModel, s, Arrays.toString(crlb), Arrays.toString(moment.getStandardDeviation())));
        }
    }
}

Example 5 with PoissonGradientProcedure

use of uk.ac.sussex.gdsc.smlm.fitting.nonlinear.gradient.PoissonGradientProcedure in project GDSC-SMLM by aherbert.

the class MleLvmSteppingFunctionSolver method computeLastFisherInformationMatrix.

@Override
protected FisherInformationMatrix computeLastFisherInformationMatrix(double[] fx) {
    // The Hessian matrix refers to the log-likelihood ratio.
    // Compute and invert a matrix related to the Poisson log-likelihood.
    // This assumes this does achieve the maximum likelihood estimate for a
    // Poisson process.
    Gradient1Function localF1 = (Gradient1Function) function;
    // Capture the y-values if necessary
    if (fx != null && fx.length == localF1.size()) {
        localF1 = new Gradient2FunctionValueStore(localF1, fx);
    }
    // Add the weights if necessary
    if (weights != null) {
        localF1 = OffsetGradient1Function.wrapGradient1Function(localF1, weights);
    }
    final PoissonGradientProcedure p = PoissonGradientProcedureUtils.create(localF1);
    p.computeFisherInformation(lastA);
    if (p.isNaNGradients()) {
        throw new FunctionSolverException(FitStatus.INVALID_GRADIENTS);
    }
    // Re-use space
    p.getLinear(walpha);
    return new FisherInformationMatrix(walpha, beta.length);
}