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

Example 1 with StandardValueProcedure

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

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

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

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

the class CubicSplineFunctionTest method functionComputesTargetGradient1With2Peaks.

private void functionComputesTargetGradient1With2Peaks(int targetParameter) {
    final int gradientIndex = findGradientIndex(f2, targetParameter);
    final Statistics s = new Statistics();
    final StandardValueProcedure p1a = new StandardValueProcedure();
    final StandardValueProcedure p1b = new StandardValueProcedure();
    final StandardGradient1Procedure p2 = new StandardGradient1Procedure();
    for (final double background : testbackground) {
        // Peak 1
        for (final double signal1 : testsignal1) {
            for (final double cx1 : testcx1) {
                for (final double cy1 : testcy1) {
                    for (final double cz1 : testcz1) {
                        // Peak 2
                        for (final double signal2 : testsignal2) {
                            for (final double cx2 : testcx2) {
                                for (final double cy2 : testcy2) {
                                    for (final double cz2 : testcz2) {
                                        final double[] a = createParameters(background, signal1, cx1, cy1, cz1, signal2, cx2, cy2, cz2);
                                        // System.out.println(java.util.Arrays.toString(a));
                                        // Evaluate all gradients
                                        p2.getValues(f2, a);
                                        // Numerically solve gradient.
                                        // Calculate the step size h to be an exact numerical representation
                                        final double xx = a[targetParameter];
                                        // Get h to minimise roundoff error
                                        final double h = Precision.representableDelta(xx, stepH);
                                        // Evaluate at (x+h) and (x-h)
                                        a[targetParameter] = xx + h;
                                        p1a.getValues(f2, a);
                                        a[targetParameter] = xx - h;
                                        p1b.getValues(f2, a);
                                        // Only test close to the XY centre
                                        for (final int x : testx) {
                                            for (final int y : testy) {
                                                final int i = y * maxx + x;
                                                final double high = p1a.values[i];
                                                final double low = p1b.values[i];
                                                final double gradient = (high - low) / (2 * h);
                                                final double dyda = p2.gradients[i][gradientIndex];
                                                final double error = DoubleEquality.relativeError(gradient, dyda);
                                                s.add(error);
                                                Assertions.assertTrue((gradient * dyda) >= 0, () -> String.format("%s sign != %s", gradient, dyda));
                                                // logger.fine(FunctionUtils.getSupplier("[%d,%d] %f == [%d] %f? (%g)", x,
                                                // y, gradient, gradientIndex, dyda, error);
                                                Assertions.assertTrue(eq.almostEqualRelativeOrAbsolute(gradient, dyda), () -> String.format("%s != %s", gradient, dyda));
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    logger.info(() -> {
        return String.format("functionComputesTargetGradient1With2Peaks %s %s (error %s +/- %s)", f1.getClass().getSimpleName(), CubicSplineFunction.getName(targetParameter), MathUtils.rounded(s.getMean()), MathUtils.rounded(s.getStandardDeviation()));
    });
}

Example 5 with StandardValueProcedure

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

the class CubicSplineFunctionTest method functionComputesTargetGradient1.

private void functionComputesTargetGradient1(int targetParameter) {
    final int gradientIndex = findGradientIndex(f1, targetParameter);
    final Statistics s = new Statistics();
    final StandardValueProcedure p1a = new StandardValueProcedure();
    final StandardValueProcedure p1b = new StandardValueProcedure();
    final StandardGradient1Procedure p2 = new StandardGradient1Procedure();
    for (final double background : testbackground) {
        // Peak 1
        for (final double signal1 : testsignal1) {
            for (final double cx1 : testcx1) {
                for (final double cy1 : testcy1) {
                    for (final double cz1 : testcz1) {
                        final double[] a = createParameters(background, signal1, cx1, cy1, cz1);
                        // System.out.println(java.util.Arrays.toString(a));
                        // Evaluate all gradients
                        p2.getValues(f1, a);
                        // Numerically solve gradient.
                        // Calculate the step size h to be an exact numerical representation
                        final double xx = a[targetParameter];
                        // Get h to minimise roundoff error
                        final double h = Precision.representableDelta(xx, stepH);
                        // Evaluate at (x+h) and (x-h)
                        a[targetParameter] = xx + h;
                        p1a.getValues(f1, a);
                        a[targetParameter] = xx - h;
                        p1b.getValues(f1, a);
                        // Only test close to the XY centre
                        for (final int x : testx) {
                            for (final int y : testy) {
                                final int i = y * maxx + x;
                                final double high = p1a.values[i];
                                final double low = p1b.values[i];
                                final double gradient = (high - low) / (2 * h);
                                final double dyda = p2.gradients[i][gradientIndex];
                                final double error = DoubleEquality.relativeError(gradient, dyda);
                                s.add(error);
                                if ((gradient * dyda) < 0) {
                                    Assertions.fail(String.format("%s sign != %s", gradient, dyda));
                                }
                                // gradient, gradientIndex, dyda, error);
                                if (!eq.almostEqualRelativeOrAbsolute(gradient, dyda)) {
                                    Assertions.fail(String.format("%s != %s", gradient, dyda));
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    logger.info(() -> {
        return String.format("functionComputesTargetGradient1 %s %s (error %s +/- %s)", f1.getClass().getSimpleName(), CubicSplineFunction.getName(targetParameter), MathUtils.rounded(s.getMean()), MathUtils.rounded(s.getStandardDeviation()));
    });
}