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

Example 1 with Gradient1Procedure

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

the class UnivariateLikelihoodFisherInformationCalculator method compute.

/**
 * {@inheritDoc}
 *
 * @throws DataException If the Fisher information cannot be computed for a function value
 * @throws DataException If the Fisher information is infinite for a function value
 */
@Override
public FisherInformationMatrix compute(double[] parameters) {
    final int n = gf.getNumberOfGradients();
    final double[] data = new double[n * (n + 1) / 2];
    gf.initialise1(parameters);
    gf.forEach(new Gradient1Procedure() {

        // CHECKSTYLE.OFF: MemberName
        int k = -1;

        // CHECKSTYLE.ON: MemberName
        @Override
        public void execute(double value, double[] dvDt) {
            k++;
            if (!fi[k].isValid(value)) {
                return;
            }
            // Get the Fisher information of the value
            final double f = fi[k].getFisherInformation(value);
            if (f == 0) {
                // No summation
                return;
            }
            if (f == Double.POSITIVE_INFINITY) {
                throw new DataException("Fisher information is infinite at f(" + k + ")");
            }
            // Compute the actual matrix data
            for (int i = 0, c = 0; i < n; i++) {
                final double wgt = f * dvDt[i];
                for (int j = 0; j <= i; j++) {
                    data[c++] += wgt * dvDt[j];
                }
            }
        }
    });
    // Generate symmetric matrix
    final double[] matrix = new double[n * n];
    for (int i = 0, c = 0; i < n; i++) {
        for (int j = 0; j <= i; j++) {
            matrix[i * n + j] = matrix[j * n + i] = data[c++];
        }
    }
    return new FisherInformationMatrix(matrix, n);
}

Example 2 with Gradient1Procedure

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

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

the class Gaussian2DFunction method computeValuesAndJacobian.

@Override
public Pair<double[], double[][]> computeValuesAndJacobian(double[] variables) {
    initialise1(variables);
    final int n = size();
    final double[][] jacobian = new double[n][];
    final double[] values = new double[n];
    forEach(new Gradient1Procedure() {

        int index;

        @Override
        public void execute(double value, double[] derivative) {
            values[index] = value;
            jacobian[index++] = derivative.clone();
        }
    });
    return Pair.of(values, jacobian);
}

Example 4 with Gradient1Procedure

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

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

the class PsfModelGradient1FunctionTest method canComputeValueAndGradient.

@Test
void canComputeValueAndGradient() {
    // Use a reasonable z-depth function from the Smith, et al (2010) paper (page 377)
    final double sx = 1.08;
    final double sy = 1.01;
    final double gamma = 0.389;
    final double d = 0.531;
    final double Ax = -0.0708;
    final double Bx = -0.073;
    final double Ay = 0.164;
    final double By = 0.0417;
    final AstigmatismZModel zModel = HoltzerAstigmatismZModel.create(sx, sy, gamma, d, Ax, Bx, Ay, By);
    // Small size ensure the PSF model covers the entire image
    final int maxx = 11;
    final int maxy = 11;
    final double[] ve = new double[maxx * maxy];
    final double[] vo = new double[maxx * maxy];
    final double[][] ge = new double[maxx * maxy][];
    final double[][] go = new double[maxx * maxy][];
    final PsfModelGradient1Function psf = new PsfModelGradient1Function(new GaussianPsfModel(zModel), maxx, maxy);
    final ErfGaussian2DFunction f = new SingleAstigmatismErfGaussian2DFunction(maxx, maxy, zModel);
    f.setErfFunction(ErfFunction.COMMONS_MATH);
    final double[] a2 = new double[Gaussian2DFunction.PARAMETERS_PER_PEAK + 1];
    final DoubleDoubleBiPredicate equality = TestHelper.doublesAreClose(1e-8, 0);
    final double c = maxx * 0.5;
    for (int i = -1; i <= 1; i++) {
        final double x0 = c + i * 0.33;
        for (int j = -1; j <= 1; j++) {
            final double x1 = c + j * 0.33;
            for (int k = -1; k <= 1; k++) {
                final double x2 = k * 0.33;
                for (final double in : new double[] { 23.2, 405.67 }) {
                    // Background is constant for gradients so just use 1 value
                    final double[] a = new double[] { 2.2, in, x0, x1, x2 };
                    psf.initialise1(a);
                    psf.forEach(new Gradient1Procedure() {

                        int index = 0;

                        @Override
                        public void execute(double value, double[] dyDa) {
                            vo[index] = value;
                            go[index] = dyDa.clone();
                            index++;
                        }
                    });
                    a2[Gaussian2DFunction.BACKGROUND] = a[0];
                    a2[Gaussian2DFunction.SIGNAL] = a[1];
                    a2[Gaussian2DFunction.X_POSITION] = a[2] - 0.5;
                    a2[Gaussian2DFunction.Y_POSITION] = a[3] - 0.5;
                    a2[Gaussian2DFunction.Z_POSITION] = a[4];
                    f.initialise1(a2);
                    f.forEach(new Gradient1Procedure() {

                        int index = 0;

                        @Override
                        public void execute(double value, double[] dyDa) {
                            ve[index] = value;
                            ge[index] = dyDa.clone();
                            index++;
                        }
                    });
                    for (int ii = 0; ii < ve.length; ii++) {
                        TestAssertions.assertTest(ve[ii], vo[ii], equality);
                        TestAssertions.assertArrayTest(ge[ii], go[ii], equality);
                    }
                }
            }
        }
    }
}