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

Example 1 with ValueProcedure

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

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

the class ApacheLvmFitter method computeFit.

@Override
public FitStatus computeFit(double[] y, final double[] fx, double[] a, double[] parametersVariance) {
    final int n = y.length;
    try {
        // Different convergence thresholds seem to have no effect on the resulting fit, only the
        // number of
        // iterations for convergence
        final double initialStepBoundFactor = 100;
        final double costRelativeTolerance = 1e-10;
        final double parRelativeTolerance = 1e-10;
        final double orthoTolerance = 1e-10;
        final double threshold = Precision.SAFE_MIN;
        // Extract the parameters to be fitted
        final double[] initialSolution = getInitialSolution(a);
        // TODO - Pass in more advanced stopping criteria.
        // Create the target and weight arrays
        final double[] yd = new double[n];
        // final double[] w = new double[n];
        for (int i = 0; i < n; i++) {
            yd[i] = y[i];
        // w[i] = 1;
        }
        final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer(initialStepBoundFactor, costRelativeTolerance, parRelativeTolerance, orthoTolerance, threshold);
        // @formatter:off
        final LeastSquaresBuilder builder = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(getMaxEvaluations()).start(initialSolution).target(yd);
        if (function instanceof ExtendedNonLinearFunction && ((ExtendedNonLinearFunction) function).canComputeValuesAndJacobian()) {
            // Compute together, or each individually
            builder.model(new ValueAndJacobianFunction() {

                final ExtendedNonLinearFunction fun = (ExtendedNonLinearFunction) function;

                @Override
                public Pair<RealVector, RealMatrix> value(RealVector point) {
                    final double[] p = point.toArray();
                    final org.apache.commons.lang3.tuple.Pair<double[], double[][]> result = fun.computeValuesAndJacobian(p);
                    return new Pair<>(new ArrayRealVector(result.getKey(), false), new Array2DRowRealMatrix(result.getValue(), false));
                }

                @Override
                public RealVector computeValue(double[] params) {
                    return new ArrayRealVector(fun.computeValues(params), false);
                }

                @Override
                public RealMatrix computeJacobian(double[] params) {
                    return new Array2DRowRealMatrix(fun.computeJacobian(params), false);
                }
            });
        } else {
            // Compute separately
            builder.model(new MultivariateVectorFunctionWrapper((NonLinearFunction) function, a, n), new MultivariateMatrixFunctionWrapper((NonLinearFunction) function, a, n));
        }
        final LeastSquaresProblem problem = builder.build();
        final Optimum optimum = optimizer.optimize(problem);
        final double[] parameters = optimum.getPoint().toArray();
        setSolution(a, parameters);
        iterations = optimum.getIterations();
        evaluations = optimum.getEvaluations();
        if (parametersVariance != null) {
            // Set up the Jacobian.
            final RealMatrix j = optimum.getJacobian();
            // Compute transpose(J)J.
            final RealMatrix jTj = j.transpose().multiply(j);
            final double[][] data = (jTj instanceof Array2DRowRealMatrix) ? ((Array2DRowRealMatrix) jTj).getDataRef() : jTj.getData();
            final FisherInformationMatrix m = new FisherInformationMatrix(data);
            setDeviations(parametersVariance, m);
        }
        // Compute function value
        if (fx != null) {
            final ValueFunction function = (ValueFunction) this.function;
            function.initialise0(a);
            function.forEach(new ValueProcedure() {

                int index;

                @Override
                public void execute(double value) {
                    fx[index++] = value;
                }
            });
        }
        // As this is unweighted then we can do this to get the sum of squared residuals
        // This is the same as optimum.getCost() * optimum.getCost(); The getCost() function
        // just computes the dot product anyway.
        value = optimum.getResiduals().dotProduct(optimum.getResiduals());
    } catch (final TooManyEvaluationsException ex) {
        return FitStatus.TOO_MANY_EVALUATIONS;
    } catch (final TooManyIterationsException ex) {
        return FitStatus.TOO_MANY_ITERATIONS;
    } catch (final ConvergenceException ex) {
        // Occurs when QR decomposition fails - mark as a singular non-linear model (no solution)
        return FitStatus.SINGULAR_NON_LINEAR_MODEL;
    } catch (final Exception ex) {
        // TODO - Find out the other exceptions from the fitter and add return values to match.
        return FitStatus.UNKNOWN;
    }
    return FitStatus.OK;
}

Example 3 with ValueProcedure

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

the class FastMleGradient2ProcedureTest method gradientProcedureComputesSameWithPrecomputed.

@SeededTest
void gradientProcedureComputesSameWithPrecomputed(RandomSeed seed) {
    final int iter = 10;
    final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
    final ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final ErfGaussian2DFunction f2 = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(2, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final double[] a1 = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    final double[] a2 = new double[1 + 2 * Gaussian2DFunction.PARAMETERS_PER_PEAK];
    final double[] x = new double[f1.size()];
    final double[] b = new double[f1.size()];
    for (int i = 0; i < iter; i++) {
        final int ii = i;
        a2[Gaussian2DFunction.BACKGROUND] = nextUniform(rng, 0.1, 0.3);
        a2[Gaussian2DFunction.SIGNAL] = nextUniform(rng, 100, 300);
        a2[Gaussian2DFunction.X_POSITION] = nextUniform(rng, 3, 5);
        a2[Gaussian2DFunction.Y_POSITION] = nextUniform(rng, 3, 5);
        a2[Gaussian2DFunction.Z_POSITION] = nextUniform(rng, -2, 2);
        a2[Gaussian2DFunction.X_SD] = nextUniform(rng, 1, 1.3);
        a2[Gaussian2DFunction.Y_SD] = nextUniform(rng, 1, 1.3);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.SIGNAL] = nextUniform(rng, 100, 300);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.X_POSITION] = nextUniform(rng, 5, 7);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Y_POSITION] = nextUniform(rng, 5, 7);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Z_POSITION] = nextUniform(rng, -3, 1);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.X_SD] = nextUniform(rng, 1, 1.3);
        a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Y_SD] = nextUniform(rng, 1, 1.3);
        // Simulate Poisson data
        f2.initialise0(a2);
        f1.forEach(new ValueProcedure() {

            int index = 0;

            @Override
            public void execute(double value) {
                if (value > 0) {
                    x[index++] = GdscSmlmTestUtils.createPoissonSampler(rng, value).sample();
                } else {
                    x[index++] = 0;
                }
            }
        });
        // Precompute peak 2 (no background)
        a1[Gaussian2DFunction.BACKGROUND] = 0;
        for (int j = 1; j < 7; j++) {
            a1[j] = a2[Gaussian2DFunction.PARAMETERS_PER_PEAK + j];
        }
        f1.initialise0(a1);
        f1.forEach(new ValueProcedure() {

            int index = 0;

            @Override
            public void execute(double value) {
                b[index++] = value;
            }
        });
        // Reset to peak 1
        for (int j = 0; j < 7; j++) {
            a1[j] = a2[j];
        }
        // Compute peak 1+2
        final FastMleGradient2Procedure p12 = FastMleGradient2ProcedureUtils.create(x, f2);
        p12.computeSecondDerivative(a2);
        final double[] d11 = Arrays.copyOf(p12.d1, f1.getNumberOfGradients());
        final double[] d21 = Arrays.copyOf(p12.d2, f1.getNumberOfGradients());
        // Compute peak 1+(precomputed 2)
        final FastMleGradient2Procedure p1b2 = FastMleGradient2ProcedureUtils.create(x, OffsetGradient2Function.wrapGradient2Function(f1, b));
        p1b2.computeSecondDerivative(a1);
        final double[] d12 = p1b2.d1;
        final double[] d22 = p1b2.d2;
        Assertions.assertArrayEquals(p12.u, p1b2.u, 1e-10, () -> " Result: Not same @ " + ii);
        Assertions.assertArrayEquals(d11, d12, 1e-10, () -> " D1: Not same @ " + ii);
        Assertions.assertArrayEquals(d21, d22, 1e-10, () -> " D2: Not same @ " + ii);
        final double[] v1 = p12.computeValue(a2);
        final double[] v2 = p1b2.computeValue(a1);
        Assertions.assertArrayEquals(v1, v2, 1e-10, () -> " Value: Not same @ " + ii);
        final double[] d1 = Arrays.copyOf(p12.computeFirstDerivative(a2), f1.getNumberOfGradients());
        final double[] d2 = p1b2.computeFirstDerivative(a1);
        Assertions.assertArrayEquals(d1, d2, 1e-10, () -> " 1st derivative: Not same @ " + ii);
    }
}

Example 4 with ValueProcedure

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

the class Gaussian2DFunction method computeValues.

@Override
public double[] computeValues(double[] variables) {
    initialise0(variables);
    final double[] values = new double[size()];
    forEach(new ValueProcedure() {

        int index;

        @Override
        public void execute(double value) {
            values[index++] = value;
        }
    });
    return values;
}

Example 5 with ValueProcedure

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

the class ErfGaussian2DFunctionTest method functionComputesGradientForEachWith2Peaks.

@Test
void functionComputesGradientForEachWith2Peaks() {
    Assumptions.assumeTrue(null != f2);
    final ErfGaussian2DFunction f2 = (ErfGaussian2DFunction) this.f2;
    final int n = f2.size();
    final double[] du_da = new double[f2.getNumberOfGradients()];
    final double[] du_db = new double[f2.getNumberOfGradients()];
    final double[] d2u_da2 = new double[f2.getNumberOfGradients()];
    final double[] values = new double[n];
    final double[][] jacobian = new double[n][];
    final double[][] jacobian2 = new double[n][];
    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) {
                        for (final double[] w1 : testw1) {
                            for (final double angle1 : testangle1) {
                                // Peak 2
                                for (final double signal2 : testsignal2) {
                                    for (final double cx2 : testcx2) {
                                        for (final double cy2 : testcy2) {
                                            for (final double cz2 : testcz2) {
                                                for (final double[] w2 : testw2) {
                                                    for (final double angle2 : testangle2) {
                                                        final double[] a = createParameters(background, signal1, cx1, cy1, cz1, w1[0], w1[1], angle1, signal2, cx2, cy2, cz2, w2[0], w2[1], angle2);
                                                        f2.initialiseExtended2(a);
                                                        // Compute single
                                                        for (int i = 0; i < n; i++) {
                                                            final double o1 = f2.eval(i, du_da);
                                                            final double o2 = f2.eval2(i, du_db, d2u_da2);
                                                            Assertions.assertEquals(o1, o2, 1e-10, "Value");
                                                            Assertions.assertArrayEquals(du_da, du_db, 1e-10, "Jacobian!=Jacobian");
                                                            values[i] = o1;
                                                            jacobian[i] = du_da.clone();
                                                            jacobian2[i] = d2u_da2.clone();
                                                        }
                                                        // Use procedures
                                                        f2.forEach(new ValueProcedure() {

                                                            int index = 0;

                                                            @Override
                                                            public void execute(double value) {
                                                                Assertions.assertEquals(values[index], value, 1e-10, "Value ValueProcedure");
                                                                index++;
                                                            }
                                                        });
                                                        f2.forEach(new Gradient1Procedure() {

                                                            int index = 0;

                                                            @Override
                                                            public void execute(double value, double[] dyDa) {
                                                                Assertions.assertEquals(values[index], value, 1e-10, "Value Gradient1Procedure");
                                                                Assertions.assertArrayEquals(jacobian[index], dyDa, 1e-10, "du_da Gradient1Procedure");
                                                                index++;
                                                            }
                                                        });
                                                        f2.forEach(new Gradient2Procedure() {

                                                            int index = 0;

                                                            @Override
                                                            public void execute(double value, double[] dyDa, double[] d2yDa2) {
                                                                Assertions.assertEquals(values[index], value, 1e-10, "Value Gradient2Procedure");
                                                                Assertions.assertArrayEquals(jacobian[index], dyDa, 1e-10, "du_da Gradient2Procedure");
                                                                Assertions.assertArrayEquals(jacobian2[index], d2yDa2, 1e-10, "d2u_da2 Gradient2Procedure");
                                                                index++;
                                                            }
                                                        });
                                                        f2.forEach(new ExtendedGradient2Procedure() {

                                                            int index = 0;

                                                            @Override
                                                            public void executeExtended(double value, double[] dyDa, double[] d2yDaDb) {
                                                                Assertions.assertEquals(values[index], value, 1e-10, "Value ExtendedGradient2Procedure");
                                                                Assertions.assertArrayEquals(jacobian[index], dyDa, 1e-10, "du_da ExtendedGradient2Procedure");
                                                                for (int j = 0, k = 0; j < d2u_da2.length; j++, k += d2u_da2.length + 1) {
                                                                    d2u_da2[j] = d2yDaDb[k];
                                                                }
                                                                Assertions.assertArrayEquals(jacobian2[index], d2u_da2, 1e-10, "d2u_da2 Gradient2Procedure");
                                                                index++;
                                                            }
                                                        });
                                                    }
                                                }
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
}