Examples with ExtendedGradient2Procedure gdsc.smlm.function.ExtendedGradient2Procedure used on opensource projects

Example 1 with ExtendedGradient2Procedure

use of gdsc.smlm.function.ExtendedGradient2Procedure in project GDSC-SMLM by aherbert.

the class ErfGaussian2DFunctionTest method functionComputesExtendedGradientForEach.

@Test
public void functionComputesExtendedGradientForEach() {
    final ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) this.f1;
    final int nparams = f1.getNumberOfGradients();
    final int[] gradientIndices = f1.gradientIndices();
    final double[] du_da = new double[f1.getNumberOfGradients()];
    final double[] du_db = new double[f1.getNumberOfGradients()];
    final ErfGaussian2DFunction[] fHigh = new ErfGaussian2DFunction[nparams];
    final ErfGaussian2DFunction[] fLow = new ErfGaussian2DFunction[nparams];
    final double[] delta = new double[nparams];
    for (int j = 0; j < nparams; j++) {
        fHigh[j] = f1.copy();
        fLow[j] = f1.copy();
    }
    for (double background : testbackground) // Peak 1
    for (double amplitude1 : testamplitude1) for (double shape1 : testshape1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double[] w1 : testw1) {
        double[] a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1]);
        f1.initialiseExtended2(a);
        // Create a set of functions initialised +/- delta in each parameter
        for (int j = 0; j < nparams; j++) {
            int targetParameter = gradientIndices[j];
            // 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
            double h = Precision.representableDelta(xx, h_);
            // Evaluate at (x+h) and (x-h)
            a[targetParameter] = xx + h;
            fHigh[j].initialise1(a);
            a[targetParameter] = xx - h;
            fLow[j].initialise1(a);
            a[targetParameter] = xx;
            delta[j] = 2 * h;
        }
        f1.forEach(new ExtendedGradient2Procedure() {

            int i = -1;

            public void executeExtended(double value, double[] dy_da, double[] d2y_dadb) {
                i++;
                DenseMatrix64F m = DenseMatrix64F.wrap(nparams, nparams, d2y_dadb);
                for (int j = 0; j < nparams; j++) {
                    // Evaluate the function +/- delta for parameter j
                    fHigh[j].eval(i, du_da);
                    fLow[j].eval(i, du_db);
                    // Check the gradient with respect to parameter k
                    for (int k = 0; k < nparams; k++) {
                        double gradient = (du_da[k] - du_db[k]) / delta[j];
                        boolean ok = eq.almostEqualRelativeOrAbsolute(gradient, m.get(j, k));
                        if (!ok) {
                            System.out.printf("%d [%d,%d] %f ?= %f\n", i, j, k, gradient, m.get(j, k));
                            Assert.fail(String.format("%d [%d,%d] %f != %f", i, j, k, gradient, m.get(j, k)));
                        }
                    }
                }
            }
        });
    }
}

Example 2 with ExtendedGradient2Procedure

use of gdsc.smlm.function.ExtendedGradient2Procedure in project GDSC-SMLM by aherbert.

the class ErfGaussian2DFunctionTest method functionComputesGradientForEachWith2Peaks.

@Test
public void functionComputesGradientForEachWith2Peaks() {
    org.junit.Assume.assumeNotNull(f2);
    final ErfGaussian2DFunction f2 = (ErfGaussian2DFunction) this.f2;
    final int n = f2.size();
    double[] du_da = new double[f2.getNumberOfGradients()];
    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 (double background : testbackground) // Peak 1
    for (double amplitude1 : testamplitude1) for (double shape1 : testshape1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double[] w1 : testw1) // Peak 2
    for (double amplitude2 : testamplitude2) for (double shape2 : testshape2) for (double cx2 : testcx2) for (double cy2 : testcy2) for (double[] w2 : testw2) {
        double[] a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1], amplitude2, shape2, cx2, cy2, w2[0], w2[1]);
        f2.initialiseExtended2(a);
        // Compute single
        for (int i = 0; i < n; i++) {
            double o1 = f2.eval(i, du_da);
            double o2 = f2.eval(i, du_db, d2u_da2);
            Assert.assertEquals("Value", o1, o2, 1e-10);
            Assert.assertArrayEquals("Jacobian!=Jacobian", du_da, du_db, 1e-10);
            values[i] = o1;
            jacobian[i] = du_da.clone();
            jacobian2[i] = d2u_da2.clone();
        }
        // Use procedures
        f2.forEach(new ValueProcedure() {

            int i = 0;

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

            int i = 0;

            public void execute(double value, double[] dy_da) {
                Assert.assertEquals("Value Gradient1Procedure", values[i], value, 1e-10);
                Assert.assertArrayEquals("du_da Gradient1Procedure", jacobian[i], dy_da, 1e-10);
                i++;
            }
        });
        f2.forEach(new Gradient2Procedure() {

            int i = 0;

            public void execute(double value, double[] dy_da, double[] d2y_da2) {
                Assert.assertEquals("Value Gradient2Procedure", values[i], value, 1e-10);
                Assert.assertArrayEquals("du_da Gradient2Procedure", jacobian[i], dy_da, 1e-10);
                Assert.assertArrayEquals("d2u_da2 Gradient2Procedure", jacobian2[i], d2y_da2, 1e-10);
                i++;
            }
        });
        f2.forEach(new ExtendedGradient2Procedure() {

            int i = 0;

            public void executeExtended(double value, double[] dy_da, double[] d2y_dadb) {
                Assert.assertEquals("Value ExtendedGradient2Procedure", values[i], value, 1e-10);
                Assert.assertArrayEquals("du_da ExtendedGradient2Procedure", jacobian[i], dy_da, 1e-10);
                for (int j = 0, k = 0; j < d2u_da2.length; j++, k += d2u_da2.length + 1) d2u_da2[j] = d2y_dadb[k];
                Assert.assertArrayEquals("d2u_da2 Gradient2Procedure", jacobian2[i], d2u_da2, 1e-10);
                i++;
            }
        });
    }
}

Example 3 with ExtendedGradient2Procedure

use of gdsc.smlm.function.ExtendedGradient2Procedure in project GDSC-SMLM by aherbert.

the class ErfGaussian2DFunctionTest method functionComputesExtendedGradientForEachWith2Peaks.

@Test
public void functionComputesExtendedGradientForEachWith2Peaks() {
    org.junit.Assume.assumeNotNull(f2);
    final ErfGaussian2DFunction f2 = (ErfGaussian2DFunction) this.f2;
    final int nparams = f2.getNumberOfGradients();
    final int[] gradientIndices = f2.gradientIndices();
    final double[] du_da = new double[f2.getNumberOfGradients()];
    final double[] du_db = new double[f2.getNumberOfGradients()];
    final ErfGaussian2DFunction[] fHigh = new ErfGaussian2DFunction[nparams];
    final ErfGaussian2DFunction[] fLow = new ErfGaussian2DFunction[nparams];
    final double[] delta = new double[nparams];
    for (int j = 0; j < nparams; j++) {
        fHigh[j] = f2.copy();
        fLow[j] = f2.copy();
    }
    for (double background : testbackground) // Peak 1
    for (double amplitude1 : testamplitude1) for (double shape1 : testshape1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double[] w1 : testw1) // Peak 2
    for (double amplitude2 : testamplitude2) for (double shape2 : testshape2) for (double cx2 : testcx2) for (double cy2 : testcy2) for (double[] w2 : testw2) {
        double[] a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1], amplitude2, shape2, cx2, cy2, w2[0], w2[1]);
        f2.initialiseExtended2(a);
        // Create a set of functions initialised +/- delta in each parameter
        for (int j = 0; j < nparams; j++) {
            int targetParameter = gradientIndices[j];
            // 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
            double h = Precision.representableDelta(xx, h_);
            // Evaluate at (x+h) and (x-h)
            a[targetParameter] = xx + h;
            fHigh[j].initialise1(a);
            a[targetParameter] = xx - h;
            fLow[j].initialise1(a);
            a[targetParameter] = xx;
            delta[j] = 2 * h;
        }
        f2.forEach(new ExtendedGradient2Procedure() {

            int i = -1;

            public void executeExtended(double value, double[] dy_da, double[] d2y_dadb) {
                i++;
                //if (i!=f2.size()/2) return;
                DenseMatrix64F m = DenseMatrix64F.wrap(nparams, nparams, d2y_dadb);
                for (int j = 0; j < nparams; j++) {
                    // Evaluate the function +/- delta for parameter j
                    fHigh[j].eval(i, du_da);
                    fLow[j].eval(i, du_db);
                    // Check the gradient with respect to parameter k
                    for (int k = 0; k < nparams; k++) {
                        double gradient = (du_da[k] - du_db[k]) / delta[j];
                        boolean ok = eq.almostEqualRelativeOrAbsolute(gradient, m.get(j, k));
                        if (!ok) {
                            System.out.printf("%d [%d,%d] %f ?= %f\n", i, j, k, gradient, m.get(j, k));
                            Assert.fail(String.format("%d [%d,%d] %f != %f", i, j, k, gradient, m.get(j, k)));
                        }
                    }
                }
            }
        });
    //return;
    }
}