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

Example 11 with Gradient1Function

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

the class LvmGradientProcedureTest method gradientProcedureUnrolledComputesSameAsGradientProcedure.

private void gradientProcedureUnrolledComputesSameAsGradientProcedure(RandomSeed seed, int nparams, Type type, boolean precomputed) {
    final int iter = 10;
    final ArrayList<double[]> paramsList = new ArrayList<>(iter);
    final ArrayList<double[]> yList = new ArrayList<>(iter);
    createFakeData(RngUtils.create(seed.getSeed()), nparams, iter, paramsList, yList);
    Gradient1Function func = new FakeGradientFunction(blockWidth, nparams);
    if (precomputed) {
        final double[] b = SimpleArrayUtils.newArray(func.size(), 0.1, 1.3);
        func = OffsetGradient1Function.wrapGradient1Function(func, b);
    }
    final FastLog fastLog = type == Type.FAST_LOG_MLE ? getFastLog() : null;
    final String name = String.format("[%d] %b", nparams, type);
    // Create messages
    final IndexSupplier msgR = new IndexSupplier(1, name + "Result: Not same ", null);
    final IndexSupplier msgOb = new IndexSupplier(1, name + "Observations: Not same beta ", null);
    final IndexSupplier msgOal = new IndexSupplier(1, name + "Observations: Not same alpha linear ", null);
    final IndexSupplier msgOam = new IndexSupplier(1, name + "Observations: Not same alpha matrix ", null);
    for (int i = 0; i < paramsList.size(); i++) {
        final LvmGradientProcedure p1 = createProcedure(type, yList.get(i), func, fastLog);
        p1.gradient(paramsList.get(i));
        final LvmGradientProcedure p2 = LvmGradientProcedureUtils.create(yList.get(i), func, type, fastLog);
        p2.gradient(paramsList.get(i));
        // Exactly the same ...
        Assertions.assertEquals(p1.value, p2.value, msgR.set(0, i));
        Assertions.assertArrayEquals(p1.beta, p2.beta, msgOb.set(0, i));
        Assertions.assertArrayEquals(p1.getAlphaLinear(), p2.getAlphaLinear(), msgOal.set(0, i));
        final double[][] am1 = p1.getAlphaMatrix();
        final double[][] am2 = p2.getAlphaMatrix();
        Assertions.assertArrayEquals(am1, am2, msgOam.set(0, i));
    }
}

Example 12 with Gradient1Function

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

the class PoissonGradientProcedureTest method gradientProcedureIsFasterUnrolledThanGradientProcedure.

private void gradientProcedureIsFasterUnrolledThanGradientProcedure(RandomSeed seed, final int nparams, final boolean precomputed) {
    Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
    final int iter = 100;
    final ArrayList<double[]> paramsList = new ArrayList<>(iter);
    createFakeParams(RngUtils.create(seed.getSeed()), nparams, iter, paramsList);
    // Remove the timing of the function call by creating a dummy function
    final FakeGradientFunction f = new FakeGradientFunction(blockWidth, nparams);
    final Gradient1Function func = (precomputed) ? OffsetGradient1Function.wrapGradient1Function(f, SimpleArrayUtils.newArray(f.size(), 0.1, 1.3)) : f;
    final IntArrayFormatSupplier msg = new IntArrayFormatSupplier("M [%d]", 1);
    for (int i = 0; i < paramsList.size(); i++) {
        final PoissonGradientProcedure p1 = new PoissonGradientProcedure(func);
        p1.computeFisherInformation(paramsList.get(i));
        p1.computeFisherInformation(paramsList.get(i));
        final PoissonGradientProcedure p2 = PoissonGradientProcedureUtils.create(func);
        p2.computeFisherInformation(paramsList.get(i));
        p2.computeFisherInformation(paramsList.get(i));
        // Check they are the same
        Assertions.assertArrayEquals(p1.getLinear(), p2.getLinear(), msg.set(0, i));
    }
    // Realistic loops for an optimisation
    final int loops = 15;
    // Run till stable timing
    final Timer t1 = new Timer() {

        @Override
        void run() {
            for (int i = 0, k = 0; i < paramsList.size(); i++) {
                final PoissonGradientProcedure p1 = new PoissonGradientProcedure(func);
                for (int j = loops; j-- > 0; ) {
                    p1.computeFisherInformation(paramsList.get(k++ % iter));
                }
            }
        }
    };
    final long time1 = t1.getTime();
    final Timer t2 = new Timer(t1.loops) {

        @Override
        void run() {
            for (int i = 0, k = 0; i < paramsList.size(); i++) {
                final PoissonGradientProcedure p2 = PoissonGradientProcedureUtils.create(func);
                for (int j = loops; j-- > 0; ) {
                    p2.computeFisherInformation(paramsList.get(k++ % iter));
                }
            }
        }
    };
    final long time2 = t2.getTime();
    logger.log(TestLogUtils.getTimingRecord("precomputed=" + precomputed + " Standard " + nparams, time1, "Unrolled", time2));
// Assertions.assertTrue(time2 < time1);
}

Example 13 with Gradient1Function

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

Example 14 with Gradient1Function

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

the class MleLvmSteppingFunctionSolver method computeFunctionFisherInformationMatrix.

@Override
protected FisherInformationMatrix computeFunctionFisherInformationMatrix(double[] y, double[] a) {
    // Compute and invert a matrix related to the Poisson log-likelihood.
    // This assumes this does achieve the maximum likelihood estimate for a
    // Poisson process.
    // We must wrap the gradient function if weights are present.
    Gradient1Function localF1 = (Gradient1Function) function;
    if (weights != null) {
        localF1 = OffsetGradient1Function.wrapGradient1Function(localF1, weights);
    }
    final PoissonGradientProcedure p = PoissonGradientProcedureUtils.create(localF1);
    p.computeFisherInformation(a);
    if (p.isNaNGradients()) {
        throw new FunctionSolverException(FitStatus.INVALID_GRADIENTS);
    }
    return new FisherInformationMatrix(p.getLinear(), function.getNumberOfGradients());
}

Example 15 with Gradient1Function

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

the class UnivariateLikelihoodFisherInformationCalculatorTest method computePoissonFisherInformation.

private static void computePoissonFisherInformation(UniformRandomProvider rng, Model model) {
    // 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;
    switch(model) {
        // Get a variance
        case POISSON_GAUSSIAN:
            final double var = 0.9 + 0.2 * rng.nextDouble();
            fi = new PoissonGaussianApproximationFisherInformation(Math.sqrt(var));
            f1 = (Gradient1Function) OffsetFunctionFactory.wrapFunction(func, SimpleArrayUtils.newDoubleArray(func.size(), var));
            break;
        case POISSON:
            fi = new PoissonFisherInformation();
            break;
        case HALF_POISSON:
            fi = new HalfPoissonFisherInformation();
            break;
        default:
            throw new IllegalStateException();
    }
    // 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;
    final boolean emCcd = model == Model.HALF_POISSON;
    if (emCcd) {
        // Assumes half the poisson fisher information
        SimpleArrayUtils.multiply(e, 0.5);
    }
    Assertions.assertArrayEquals(e, o, 1e-6);
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(5e-2, 0);
    if (model == Model.POISSON || model == Model.HALF_POISSON) {
        // Get the Mortensen approximation for fitting Poisson data with a Gaussian.
        // Set a to 100 for the square pixel adjustment.
        final double a = 100;
        final double s = params[Gaussian2DFunction.X_SD] * a;
        final double N = params[Gaussian2DFunction.SIGNAL];
        final double b2 = params[Gaussian2DFunction.BACKGROUND];
        double var = Gaussian2DPeakResultHelper.getMLVarianceX(a, s, N, b2, emCcd);
        // Convert expected variance to pixels
        var /= (a * a);
        // Get the limits by inverting the Fisher information
        final double[] crlb = I.crlb();
        TestAssertions.assertTest(var, crlb[2], predicate);
        TestAssertions.assertTest(var, crlb[3], predicate);
    }
}