use of gdsc.smlm.function.gaussian.Gaussian2DFunction in project GDSC-SMLM by aherbert.
the class ErfGaussian2DFunctionTest method factoryDefaultsToErfGaussian2DFunction.
@Test
public void factoryDefaultsToErfGaussian2DFunction() {
Gaussian2DFunction f, f2;
int flags2 = BitFlags.unset(flags, GaussianFunctionFactory.FIT_ERF);
if (this.f2 != null) {
f = GaussianFunctionFactory.create2D(2, maxx, maxy, flags, zModel);
f2 = GaussianFunctionFactory.create2D(2, maxx, maxy, flags2, zModel);
Assert.assertTrue("Incorrect function2", f.getClass() == f2.getClass());
} else {
f = GaussianFunctionFactory.create2D(1, maxx, maxy, flags, zModel);
f2 = GaussianFunctionFactory.create2D(1, maxx, maxy, flags2, zModel);
Assert.assertTrue("Incorrect function1", f.getClass() == f2.getClass());
}
}
use of gdsc.smlm.function.gaussian.Gaussian2DFunction in project GDSC-SMLM by aherbert.
the class FitConfiguration method createGaussianFunction.
/**
* Creates the appropriate 2D Gaussian function for the configuration
*
* @param npeaks
* The number of peaks to fit
* @param maxx
* The width of the XY data
* @param maxx
* The height of the XY data
* @param params
* The Gaussian parameters
* @return The function
*/
public Gaussian2DFunction createGaussianFunction(int npeaks, int maxx, int maxy, double[] params) {
final int flags = getFunctionFlags();
Gaussian2DFunction f = GaussianFunctionFactory.create2D(npeaks, maxx, maxy, flags, null);
//f.initialise(params);
return f;
}
use of gdsc.smlm.function.gaussian.Gaussian2DFunction 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;
}
use of gdsc.smlm.function.gaussian.Gaussian2DFunction 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();
}
use of gdsc.smlm.function.gaussian.Gaussian2DFunction 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);
}
Aggregations