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

Example 1 with ExtendedGradient2Procedure

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

the class ErfGaussian2DFunctionTest method functionComputesExtendedGradientForEachWith2Peaks.

@Test
void functionComputesExtendedGradientForEachWith2Peaks() {
    Assumptions.assumeTrue(null != f2);
    final ErfGaussian2DFunction f2 = (ErfGaussian2DFunction) this.f2;
    final int nparams = f2.getNumberOfGradients();
    final int[] gradientIndices = f2.gradientIndices();
    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 (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);
                                                        // Create a set of functions initialised +/- delta in each parameter
                                                        for (int j = 0; j < nparams; j++) {
                                                            final 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
                                                            final double h = Precision.representableDelta(xx, stepH);
                                                            // Evaluate at (x+h) and (x-h)
                                                            a[targetParameter] = xx + h;
                                                            fHigh[j].initialise1(a.clone());
                                                            a[targetParameter] = xx - h;
                                                            fLow[j].initialise1(a.clone());
                                                            a[targetParameter] = xx;
                                                            delta[j] = 2 * h;
                                                        }
                                                        f2.forEach(new ExtendedGradient2Procedure() {

                                                            int index = -1;

                                                            final double[] duDa = new double[f2.getNumberOfGradients()];

                                                            final double[] duDb = new double[f2.getNumberOfGradients()];

                                                            @Override
                                                            public void executeExtended(double value, double[] dyDa, double[] d2yDaDb) {
                                                                index++;
                                                                // if (i!=f2.size()/2) return;
                                                                final DenseMatrix64F m = DenseMatrix64F.wrap(nparams, nparams, d2yDaDb);
                                                                for (int j = 0; j < nparams; j++) {
                                                                    // Evaluate the function +/- delta for parameter j
                                                                    fHigh[j].eval(index, duDa);
                                                                    fLow[j].eval(index, duDb);
                                                                    // Check the gradient with respect to parameter k
                                                                    for (int k = 0; k < nparams; k++) {
                                                                        final double gradient = (duDa[k] - duDb[k]) / delta[j];
                                                                        final boolean ok = eq.almostEqualRelativeOrAbsolute(gradient, m.get(j, k));
                                                                        // gradient, m.get(j, k)));
                                                                        if (!ok) {
                                                                            Assertions.fail(String.format("%d [%d,%d] %f != %f", index, j, k, gradient, m.get(j, k)));
                                                                        }
                                                                    }
                                                                }
                                                            }
                                                        });
                                                    // return;
                                                    }
                                                }
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
}

Example 2 with ExtendedGradient2Procedure

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

Example 3 with ExtendedGradient2Procedure

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

the class ErfGaussian2DFunctionTest method functionComputesGradientForEach.

@Test
void functionComputesGradientForEach() {
    final ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) this.f1;
    final int n = f1.size();
    final double[] du_da = new double[f1.getNumberOfGradients()];
    final double[] du_db = new double[f1.getNumberOfGradients()];
    final double[] d2u_da2 = new double[f1.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) {
                                final double[] a = createParameters(background, signal1, cx1, cy1, cz1, w1[0], w1[1], angle1);
                                f1.initialiseExtended2(a);
                                // Compute single
                                for (int i = 0; i < n; i++) {
                                    final double o1 = f1.eval(i, du_da);
                                    final double o2 = f1.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
                                f1.forEach(new ValueProcedure() {

                                    int index = 0;

                                    @Override
                                    public void execute(double value) {
                                        Assertions.assertEquals(values[index], value, 1e-10, "Value ValueProcedure");
                                        index++;
                                    }
                                });
                                f1.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++;
                                    }
                                });
                                f1.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++;
                                    }
                                });
                                f1.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++;
                                    }
                                });
                            }
                        }
                    }
                }
            }
        }
    }
}

Example 4 with ExtendedGradient2Procedure

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

the class ErfGaussian2DFunctionTest method functionComputesExtendedGradientForEach.

@Test
void functionComputesExtendedGradientForEach() {
    final ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) this.f1;
    final int nparams = f1.getNumberOfGradients();
    final int[] gradientIndices = f1.gradientIndices();
    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 (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) {
                                final double[] a = createParameters(background, signal1, cx1, cy1, cz1, w1[0], w1[1], angle1);
                                f1.initialiseExtended2(a);
                                // Create a set of functions initialised +/- delta in each parameter
                                for (int j = 0; j < nparams; j++) {
                                    final 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
                                    final double h = Precision.representableDelta(xx, stepH);
                                    // Evaluate at (x+h) and (x-h)
                                    a[targetParameter] = xx + h;
                                    fHigh[j].initialise1(a.clone());
                                    a[targetParameter] = xx - h;
                                    fLow[j].initialise1(a.clone());
                                    a[targetParameter] = xx;
                                    delta[j] = 2 * h;
                                }
                                f1.forEach(new ExtendedGradient2Procedure() {

                                    int index = -1;

                                    final double[] duDa = new double[f1.getNumberOfGradients()];

                                    final double[] duDb = new double[f1.getNumberOfGradients()];

                                    @Override
                                    public void executeExtended(double value, double[] dyDa, double[] d2yDaDb) {
                                        index++;
                                        final DenseMatrix64F m = DenseMatrix64F.wrap(nparams, nparams, d2yDaDb);
                                        for (int j = 0; j < nparams; j++) {
                                            // Evaluate the function +/- delta for parameter j
                                            fHigh[j].eval(index, duDa);
                                            fLow[j].eval(index, duDb);
                                            // Check the gradient with respect to parameter k
                                            for (int k = 0; k < nparams; k++) {
                                                final double gradient = (duDa[k] - duDb[k]) / delta[j];
                                                final boolean ok = eq.almostEqualRelativeOrAbsolute(gradient, m.get(j, k));
                                                if (!ok) {
                                                    logger.log(TestLogUtils.getRecord(Level.INFO, "%d [%d,%d] %f ?= %f", index, j, k, gradient, m.get(j, k)));
                                                    Assertions.fail(String.format("%d [%d,%d] %f != %f", index, j, k, gradient, m.get(j, k)));
                                                }
                                            }
                                        }
                                    }
                                });
                            }
                        }
                    }
                }
            }
        }
    }
}