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

Example 6 with ValueProcedure

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

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

the class EJMLLinearSolverTest method runSolverSpeedTest.

private void runSolverSpeedTest(int flags) {
    final Gaussian2DFunction f0 = GaussianFunctionFactory.create2D(1, 10, 10, flags, null);
    int n = f0.size();
    final double[] y = new double[n];
    final TurboList<DenseMatrix64F> aList = new TurboList<DenseMatrix64F>();
    final TurboList<DenseMatrix64F> bList = new TurboList<DenseMatrix64F>();
    double[] testbackground = new double[] { 0.2, 0.7 };
    double[] testsignal1 = new double[] { 30, 100, 300 };
    double[] testcx1 = new double[] { 4.9, 5.3 };
    double[] testcy1 = new double[] { 4.8, 5.2 };
    double[] testw1 = new double[] { 1.1, 1.2, 1.5 };
    int np = f0.getNumberOfGradients();
    GradientCalculator calc = GradientCalculatorFactory.newCalculator(np);
    final RandomDataGenerator rdg = new RandomDataGenerator(new Well19937c(30051977));
    //double lambda = 10;
    for (double background : testbackground) // Peak 1
    for (double signal1 : testsignal1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double w1 : testw1) {
        double[] p = new double[] { background, signal1, 0, cx1, cy1, w1, w1 };
        f0.initialise(p);
        f0.forEach(new ValueProcedure() {

            int i = 0;

            public void execute(double value) {
                // Poisson data 
                y[i++] = rdg.nextPoisson(value);
            }
        });
        double[][] alpha = new double[np][np];
        double[] beta = new double[np];
        //double ss = 
        calc.findLinearised(n, y, p, alpha, beta, f0);
        //System.out.printf("SS = %f\n", ss);
        // As per the LVM algorithm
        //for (int i = 0; i < np; i++)
        //	alpha[i][i] *= lambda;
        aList.add(EJMLLinearSolver.toA(alpha));
        bList.add(EJMLLinearSolver.toB(beta));
    }
    DenseMatrix64F[] a = aList.toArray(new DenseMatrix64F[aList.size()]);
    DenseMatrix64F[] b = bList.toArray(new DenseMatrix64F[bList.size()]);
    int runs = 100000 / a.length;
    TimingService ts = new TimingService(runs);
    TurboList<SolverTimingTask> tasks = new TurboList<SolverTimingTask>();
    tasks.add(new PseudoInverseSolverTimingTask(a, b));
    tasks.add(new LinearSolverTimingTask(a, b));
    tasks.add(new CholeskySolverTimingTask(a, b));
    tasks.add(new CholeskyLDLTSolverTimingTask(a, b));
    tasks.add(new DirectInversionSolverTimingTask(a, b));
    for (SolverTimingTask task : tasks) if (!task.badSolver)
        ts.execute(task);
    ts.repeat();
    ts.report();
}