Examples with ErfGaussian2DFunction gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction used on opensource projects

Example 1 with ErfGaussian2DFunction

use of gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction in project GDSC-SMLM by aherbert.

the class BaseSteppingFunctionSolverTest method getSolver.

SteppingFunctionSolver getSolver(SteppingFunctionSolverClamp clamp, SteppingFunctionSolverType type) {
    ErfGaussian2DFunction f = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, size, size, flags, null);
    ToleranceChecker tc = new ToleranceChecker(1e-5, 1e-5, 0, 0, 100);
    ParameterBounds bounds = new ParameterBounds(f);
    switch(clamp) {
        case DYNAMIC_CLAMP:
            bounds.setDynamicClamp(true);
        case CLAMP:
            bounds.setClampValues(defaultClampValues);
        case NO_CLAMP:
        default:
            break;
    }
    SteppingFunctionSolver solver;
    switch(type) {
        case LSELVM:
            solver = new LSELVMSteppingFunctionSolver(f, tc, bounds);
            break;
        case MLELVM:
            solver = new MLELVMSteppingFunctionSolver(f, tc, bounds);
            // MLE requires a positive function value so use a lower bound
            solver.setBounds(new double[7], null);
            break;
        case WLSELVM:
            solver = new WLSELVMSteppingFunctionSolver(f, tc, bounds);
            break;
        case FastMLE:
            solver = new FastMLESteppingFunctionSolver(f, tc, bounds);
            // MLE requires a positive function value so use a lower bound
            solver.setBounds(new double[7], null);
            break;
        case BTFastMLE:
            solver = new BacktrackingFastMLESteppingFunctionSolver(f, tc, bounds);
            // MLE requires a positive function value so use a lower bound
            solver.setBounds(new double[7], null);
            break;
        case JFastMLE:
            FastMLESteppingFunctionSolver s = new FastMLESteppingFunctionSolver(f, tc, bounds);
            s.enableJacobianSolution(true);
            // MLE requires a positive function value so use a lower bound
            solver = s;
            solver.setBounds(new double[7], null);
            break;
        default:
            throw new NotImplementedException();
    }
    if (solver instanceof LVMSteppingFunctionSolver)
        ((LVMSteppingFunctionSolver) solver).setInitialLambda(1);
    return solver;
}

Example 2 with ErfGaussian2DFunction

use of gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction in project GDSC-SMLM by aherbert.

the class LSQLVMGradientProcedureTest method gradientProcedureComputesSameOutputWithBias.

@Test
public void gradientProcedureComputesSameOutputWithBias() {
    ErfGaussian2DFunction func = new SingleFreeCircularErfGaussian2DFunction(blockWidth, blockWidth);
    int nparams = func.getNumberOfGradients();
    int iter = 100;
    rdg = new RandomDataGenerator(new Well19937c(30051977));
    ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
    ArrayList<double[]> yList = new ArrayList<double[]>(iter);
    ArrayList<double[]> alphaList = new ArrayList<double[]>(iter);
    ArrayList<double[]> betaList = new ArrayList<double[]>(iter);
    ArrayList<double[]> xList = new ArrayList<double[]>(iter);
    // Manipulate the background
    double defaultBackground = Background;
    try {
        Background = 1e-2;
        createData(1, iter, paramsList, yList, true);
        EJMLLinearSolver solver = new EJMLLinearSolver(1e-5, 1e-6);
        for (int i = 0; i < paramsList.size(); i++) {
            double[] y = yList.get(i);
            double[] a = paramsList.get(i);
            BaseLSQLVMGradientProcedure p = LSQLVMGradientProcedureFactory.create(y, func);
            p.gradient(a);
            double[] beta = p.beta;
            alphaList.add(p.getAlphaLinear());
            betaList.add(beta.clone());
            for (int j = 0; j < nparams; j++) {
                if (Math.abs(beta[j]) < 1e-6)
                    System.out.printf("[%d] Tiny beta %s %g\n", i, func.getName(j), beta[j]);
            }
            // Solve
            if (!solver.solve(p.getAlphaMatrix(), beta))
                throw new AssertionError();
            xList.add(beta);
        //System.out.println(Arrays.toString(beta));
        }
        //for (int b = 1; b < 1000; b *= 2)
        for (double b : new double[] { -500, -100, -10, -1, -0.1, 0, 0.1, 1, 10, 100, 500 }) {
            Statistics[] rel = new Statistics[nparams];
            Statistics[] abs = new Statistics[nparams];
            for (int i = 0; i < nparams; i++) {
                rel[i] = new Statistics();
                abs[i] = new Statistics();
            }
            for (int i = 0; i < paramsList.size(); i++) {
                double[] y = add(yList.get(i), b);
                double[] a = paramsList.get(i).clone();
                a[0] += b;
                BaseLSQLVMGradientProcedure p = LSQLVMGradientProcedureFactory.create(y, func);
                p.gradient(a);
                double[] beta = p.beta;
                double[] alpha2 = alphaList.get(i);
                double[] beta2 = betaList.get(i);
                double[] x2 = xList.get(i);
                Assert.assertArrayEquals("Beta", beta2, beta, 1e-10);
                Assert.assertArrayEquals("Alpha", alpha2, p.getAlphaLinear(), 1e-10);
                // Solve
                solver.solve(p.getAlphaMatrix(), beta);
                Assert.assertArrayEquals("X", x2, beta, 1e-10);
                for (int j = 0; j < nparams; j++) {
                    rel[j].add(DoubleEquality.relativeError(x2[j], beta[j]));
                    abs[j].add(Math.abs(x2[j] - beta[j]));
                }
            }
            for (int i = 0; i < nparams; i++) System.out.printf("Bias = %.2f : %s : Rel %g +/- %g: Abs %g +/- %g\n", b, func.getName(i), rel[i].getMean(), rel[i].getStandardDeviation(), abs[i].getMean(), abs[i].getStandardDeviation());
        }
    } finally {
        Background = defaultBackground;
    }
}

Example 3 with ErfGaussian2DFunction

use of gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction in project GDSC-SMLM by aherbert.

the class LVMGradientProcedureTest method doubleCreateGaussianData.

/**
	 * Create random elliptical Gaussian data an returns the data plus an estimate of the parameters.
	 * Only the chosen parameters are randomised and returned for a maximum of (background, amplitude, angle, xpos,
	 * ypos, xwidth, ywidth }
	 *
	 * @param npeaks
	 *            the npeaks
	 * @param params
	 *            set on output
	 * @param randomiseParams
	 *            Set to true to randomise the params
	 * @return the double[]
	 */
private double[] doubleCreateGaussianData(int npeaks, double[] params, boolean randomiseParams) {
    int n = blockWidth * blockWidth;
    // Generate a 2D Gaussian
    ErfGaussian2DFunction func = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(npeaks, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    params[0] = random(Background);
    for (int i = 0, j = 1; i < npeaks; i++, j += 6) {
        params[j] = random(Signal);
        params[j + 2] = random(Xpos);
        params[j + 3] = random(Ypos);
        params[j + 4] = random(Xwidth);
        params[j + 5] = random(Ywidth);
    }
    double[] y = new double[n];
    func.initialise(params);
    for (int i = 0; i < y.length; i++) {
        // Add random Poisson noise
        y[i] = rdg.nextPoisson(func.eval(i));
    }
    if (randomiseParams) {
        params[0] = random(params[0]);
        for (int i = 0, j = 1; i < npeaks; i++, j += 6) {
            params[j] = random(params[j]);
            params[j + 2] = random(params[j + 2]);
            params[j + 3] = random(params[j + 3]);
            params[j + 4] = random(params[j + 4]);
            params[j + 5] = random(params[j + 5]);
        }
    }
    return y;
}

Example 4 with ErfGaussian2DFunction

use of gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction in project GDSC-SMLM by aherbert.

the class FastMLEGradient2ProcedureTest method doubleCreateGaussianData.

/**
	 * Create random elliptical Gaussian data an returns the data plus an estimate of the parameters.
	 * Only the chosen parameters are randomised and returned for a maximum of (background, amplitude, angle, xpos,
	 * ypos, xwidth, ywidth }
	 *
	 * @param npeaks
	 *            the npeaks
	 * @param params
	 *            set on output
	 * @param randomiseParams
	 *            Set to true to randomise the params
	 * @return the double[]
	 */
private double[] doubleCreateGaussianData(int npeaks, double[] params, boolean randomiseParams) {
    int n = blockWidth * blockWidth;
    // Generate a 2D Gaussian
    ErfGaussian2DFunction func = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(npeaks, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    params[0] = random(Background);
    for (int i = 0, j = 1; i < npeaks; i++, j += 6) {
        params[j] = random(Signal);
        params[j + 2] = random(Xpos);
        params[j + 3] = random(Ypos);
        params[j + 4] = random(Xwidth);
        params[j + 5] = random(Ywidth);
    }
    double[] y = new double[n];
    func.initialise(params);
    for (int i = 0; i < y.length; i++) {
        // Add random Poisson noise
        y[i] = rdg.nextPoisson(func.eval(i));
    }
    if (randomiseParams) {
        params[0] = random(params[0]);
        for (int i = 0, j = 1; i < npeaks; i++, j += 6) {
            params[j] = random(params[j]);
            params[j + 2] = random(params[j + 2]);
            params[j + 3] = random(params[j + 3]);
            params[j + 4] = random(params[j + 4]);
            params[j + 5] = random(params[j + 5]);
        }
    }
    return y;
}

Example 5 with ErfGaussian2DFunction

use of gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction in project GDSC-SMLM by aherbert.

the class PoissonGradientProcedureTest method crlbIsHigherWithPrecomputed.

@Test
public void crlbIsHigherWithPrecomputed() {
    int iter = 10;
    rdg = new RandomDataGenerator(new Well19937c(30051977));
    ErfGaussian2DFunction func = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    double[] a = new double[7];
    int n = func.getNumberOfGradients();
    // Get a background
    double[] b = new double[func.size()];
    for (int i = 0; i < b.length; i++) b[i] = rdg.nextUniform(1, 2);
    for (int i = 0; i < iter; i++) {
        a[Gaussian2DFunction.BACKGROUND] = rdg.nextUniform(0.1, 0.3);
        a[Gaussian2DFunction.SIGNAL] = rdg.nextUniform(100, 300);
        a[Gaussian2DFunction.X_POSITION] = rdg.nextUniform(4, 6);
        a[Gaussian2DFunction.Y_POSITION] = rdg.nextUniform(4, 6);
        a[Gaussian2DFunction.X_SD] = rdg.nextUniform(1, 1.3);
        a[Gaussian2DFunction.Y_SD] = rdg.nextUniform(1, 1.3);
        PoissonGradientProcedure p1 = PoissonGradientProcedureFactory.create(func);
        p1.computeFisherInformation(a);
        PoissonGradientProcedure p2 = PoissonGradientProcedureFactory.create(PrecomputedGradient1Function.wrapGradient1Function(func, b));
        p2.computeFisherInformation(a);
        FisherInformationMatrix m1 = new FisherInformationMatrix(p1.getLinear(), n);
        FisherInformationMatrix m2 = new FisherInformationMatrix(p2.getLinear(), n);
        double[] crlb1 = m1.crlb();
        double[] crlb2 = m2.crlb();
        Assert.assertNotNull(crlb1);
        Assert.assertNotNull(crlb2);
        //System.out.printf("%s : %s\n", Arrays.toString(crlb1), Arrays.toString(crlb2));
        for (int j = 0; j < n; j++) Assert.assertTrue(crlb1[j] < crlb2[j]);
    }
}