Examples with GradientCalculator gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator used on opensource projects

Example 1 with GradientCalculator

use of gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator in project GDSC-SMLM by aherbert.

the class ApacheLVMFitter method computeFit.

public FitStatus computeFit(double[] y, final double[] y_fit, double[] a, double[] a_dev) {
    int n = y.length;
    try {
        // Different convergence thresholds seem to have no effect on the resulting fit, only the number of
        // iterations for convergence
        final double initialStepBoundFactor = 100;
        final double costRelativeTolerance = 1e-10;
        final double parRelativeTolerance = 1e-10;
        final double orthoTolerance = 1e-10;
        final double threshold = Precision.SAFE_MIN;
        // Extract the parameters to be fitted
        final double[] initialSolution = getInitialSolution(a);
        // TODO - Pass in more advanced stopping criteria.
        // Create the target and weight arrays
        final double[] yd = new double[n];
        final double[] w = new double[n];
        for (int i = 0; i < n; i++) {
            yd[i] = y[i];
            w[i] = 1;
        }
        LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer(initialStepBoundFactor, costRelativeTolerance, parRelativeTolerance, orthoTolerance, threshold);
        //@formatter:off
        LeastSquaresBuilder builder = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(getMaxEvaluations()).start(initialSolution).target(yd).weight(new DiagonalMatrix(w));
        if (f instanceof ExtendedNonLinearFunction && ((ExtendedNonLinearFunction) f).canComputeValuesAndJacobian()) {
            // Compute together, or each individually
            builder.model(new ValueAndJacobianFunction() {

                final ExtendedNonLinearFunction fun = (ExtendedNonLinearFunction) f;

                public Pair<RealVector, RealMatrix> value(RealVector point) {
                    final double[] p = point.toArray();
                    final Pair<double[], double[][]> result = fun.computeValuesAndJacobian(p);
                    return new Pair<RealVector, RealMatrix>(new ArrayRealVector(result.getFirst(), false), new Array2DRowRealMatrix(result.getSecond(), false));
                }

                public RealVector computeValue(double[] params) {
                    return new ArrayRealVector(fun.computeValues(params), false);
                }

                public RealMatrix computeJacobian(double[] params) {
                    return new Array2DRowRealMatrix(fun.computeJacobian(params), false);
                }
            });
        } else {
            // Compute separately
            builder.model(new MultivariateVectorFunctionWrapper((NonLinearFunction) f, a, n), new MultivariateMatrixFunctionWrapper((NonLinearFunction) f, a, n));
        }
        LeastSquaresProblem problem = builder.build();
        Optimum optimum = optimizer.optimize(problem);
        final double[] parameters = optimum.getPoint().toArray();
        setSolution(a, parameters);
        iterations = optimum.getIterations();
        evaluations = optimum.getEvaluations();
        if (a_dev != null) {
            try {
                double[][] covar = optimum.getCovariances(threshold).getData();
                setDeviationsFromMatrix(a_dev, covar);
            } catch (SingularMatrixException e) {
                // Matrix inversion failed. In order to return a solution 
                // return the reciprocal of the diagonal of the Fisher information 
                // for a loose bound on the limit 
                final int[] gradientIndices = f.gradientIndices();
                final int nparams = gradientIndices.length;
                GradientCalculator calculator = GradientCalculatorFactory.newCalculator(nparams);
                double[][] alpha = new double[nparams][nparams];
                double[] beta = new double[nparams];
                calculator.findLinearised(nparams, y, a, alpha, beta, (NonLinearFunction) f);
                FisherInformationMatrix m = new FisherInformationMatrix(alpha);
                setDeviations(a_dev, m.crlb(true));
            }
        }
        // Compute function value
        if (y_fit != null) {
            Gaussian2DFunction f = (Gaussian2DFunction) this.f;
            f.initialise0(a);
            f.forEach(new ValueProcedure() {

                int i = 0;

                public void execute(double value) {
                    y_fit[i] = value;
                }
            });
        }
        // As this is unweighted then we can do this to get the sum of squared residuals
        // This is the same as optimum.getCost() * optimum.getCost(); The getCost() function
        // just computes the dot product anyway.
        value = optimum.getResiduals().dotProduct(optimum.getResiduals());
    } catch (TooManyEvaluationsException e) {
        return FitStatus.TOO_MANY_EVALUATIONS;
    } catch (TooManyIterationsException e) {
        return FitStatus.TOO_MANY_ITERATIONS;
    } catch (ConvergenceException e) {
        // Occurs when QR decomposition fails - mark as a singular non-linear model (no solution)
        return FitStatus.SINGULAR_NON_LINEAR_MODEL;
    } catch (Exception e) {
        // TODO - Find out the other exceptions from the fitter and add return values to match. 
        return FitStatus.UNKNOWN;
    }
    return FitStatus.OK;
}

Example 2 with GradientCalculator

use of gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator in project GDSC-SMLM by aherbert.

the class ApacheLVMFitter method computeValue.

@Override
public boolean computeValue(double[] y, double[] y_fit, double[] a) {
    final int nparams = f.gradientIndices().length;
    GradientCalculator calculator = GradientCalculatorFactory.newCalculator(nparams, false);
    // Since we know the function is a Gaussian2DFunction
    value = calculator.findLinearised(y.length, y, y_fit, a, (NonLinearFunction) f);
    return true;
}

Example 3 with GradientCalculator

use of gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator in project GDSC-SMLM by aherbert.

the class EJMLLinearSolverTest method runInversionSpeedTest.

private void runInversionSpeedTest(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>();
    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));
    }
    DenseMatrix64F[] a = aList.toArray(new DenseMatrix64F[aList.size()]);
    boolean[] ignore = new boolean[a.length];
    double[][] answer = new double[a.length][];
    int runs = 100000 / a.length;
    TimingService ts = new TimingService(runs);
    TurboList<InversionTimingTask> tasks = new TurboList<InversionTimingTask>();
    tasks.add(new PseudoInverseInversionTimingTask(a, ignore, answer));
    tasks.add(new CholeskyInversionTimingTask(a, ignore, answer));
    tasks.add(new LinearInversionTimingTask(a, ignore, answer));
    tasks.add(new CholeskyLDLTInversionTimingTask(a, ignore, answer));
    tasks.add(new DirectInversionInversionTimingTask(a, ignore, answer));
    tasks.add(new DiagonalDirectInversionInversionTimingTask(a, ignore, answer));
    for (InversionTimingTask task : tasks) if (!task.badSolver)
        ts.execute(task);
    ts.repeat();
    ts.report();
}

Example 4 with GradientCalculator

use of gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator in project GDSC-SMLM by aherbert.

the class FisherInformationMatrixTest method createFisherInformationMatrix.

private FisherInformationMatrix createFisherInformationMatrix(int n, int k) {
    int maxx = 10;
    int size = maxx * maxx;
    RandomGenerator randomGenerator = new Well19937c(30051977);
    RandomDataGenerator rdg = new RandomDataGenerator(randomGenerator);
    // Use a real Gaussian function here to compute the Fisher information.
    // The matrix may be sensitive to the type of equation used.
    int npeaks = 1;
    while (1 + npeaks * 6 < n) npeaks++;
    Gaussian2DFunction f = GaussianFunctionFactory.create2D(npeaks, maxx, maxx, GaussianFunctionFactory.FIT_ELLIPTICAL, null);
    double[] a = new double[1 + npeaks * 6];
    a[Gaussian2DFunction.BACKGROUND] = rdg.nextUniform(1, 5);
    for (int i = 0, j = 0; i < npeaks; i++, j += 6) {
        a[j + Gaussian2DFunction.SIGNAL] = rdg.nextUniform(100, 300);
        a[j + Gaussian2DFunction.SHAPE] = rdg.nextUniform(-Math.PI, Math.PI);
        // Non-overlapping peaks otherwise the CRLB are poor
        a[j + Gaussian2DFunction.X_POSITION] = rdg.nextUniform(2 + i * 2, 4 + i * 2);
        a[j + Gaussian2DFunction.Y_POSITION] = rdg.nextUniform(2 + i * 2, 4 + i * 2);
        a[j + Gaussian2DFunction.X_SD] = rdg.nextUniform(1.5, 2);
        a[j + Gaussian2DFunction.Y_SD] = rdg.nextUniform(1.5, 2);
    }
    f.initialise(a);
    GradientCalculator c = GradientCalculatorFactory.newCalculator(a.length);
    double[][] I = c.fisherInformationMatrix(size, a, f);
    //System.out.printf("n=%d, k=%d, I=\n", n, k);
    //for (int i = 0; i < I.length; i++)
    //	System.out.println(Arrays.toString(I[i]));
    // Reduce to the desired size
    I = Arrays.copyOf(I, n);
    for (int i = 0; i < n; i++) I[i] = Arrays.copyOf(I[i], n);
    // Zero selected columns
    if (k > 0) {
        int[] zero = new RandomDataGenerator(randomGenerator).nextPermutation(n, k);
        for (int i : zero) {
            for (int j = 0; j < n; j++) {
                I[i][j] = I[j][i] = 0;
            }
        }
    }
    // Create matrix
    return new FisherInformationMatrix(I, 1e-3);
}

Example 5 with GradientCalculator

use of gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator 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();
}