use of gdsc.smlm.function.ValueProcedure 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.ValueProcedure 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.ValueProcedure in project GDSC-SMLM by aherbert.
the class LVMGradientProcedureTest method gradientProcedureSupportsPrecomputed.
private void gradientProcedureSupportsPrecomputed(final Type type) {
int iter = 10;
rdg = new RandomDataGenerator(new Well19937c(30051977));
ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
ArrayList<double[]> yList = new ArrayList<double[]>(iter);
// 3 peaks
createData(3, iter, paramsList, yList, true);
for (int i = 0; i < paramsList.size(); i++) {
final double[] y = yList.get(i);
// Add Gaussian read noise so we have negatives
double min = Maths.min(y);
for (int j = 0; j < y.length; j++) y[j] = y[i] - min + rdg.nextGaussian(0, Noise);
}
// We want to know that:
// y|peak1+peak2+peak3 == y|peak1+peak2+peak3(precomputed)
// We want to know when:
// y|peak1+peak2+peak3 != y-peak3|peak1+peak2
// i.e. we cannot subtract a precomputed peak from the data, it must be included in the fit
// E.G. LSQ - subtraction is OK, MLE/WLSQ - subtraction is not allowed
Gaussian2DFunction f123 = GaussianFunctionFactory.create2D(3, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
Gaussian2DFunction f12 = GaussianFunctionFactory.create2D(2, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
Gaussian2DFunction f3 = GaussianFunctionFactory.create2D(1, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
int nparams = f12.getNumberOfGradients();
int[] indices = f12.gradientIndices();
final double[] b = new double[f12.size()];
double delta = 1e-3;
DoubleEquality eq = new DoubleEquality(5e-3, 1e-6);
double[] a1peaks = new double[7];
final double[] y_b = new double[b.length];
for (int i = 0; i < paramsList.size(); i++) {
final double[] y = yList.get(i);
double[] a3peaks = paramsList.get(i);
double[] a2peaks = Arrays.copyOf(a3peaks, 1 + 2 * 6);
double[] a2peaks2 = a2peaks.clone();
for (int j = 1; j < 7; j++) a1peaks[j] = a3peaks[j + 2 * 6];
// Evaluate peak 3 to get the background and subtract it from the data to get the new data
f3.initialise0(a1peaks);
f3.forEach(new ValueProcedure() {
int k = 0;
public void execute(double value) {
b[k] = value;
// Remove negatives for MLE
if (type == Type.MLE) {
y[k] = Math.max(0, y[k]);
y_b[k] = Math.max(0, y[k] - value);
} else {
y_b[k] = y[k] - value;
}
k++;
}
});
// These should be the same
LVMGradientProcedure p123 = LVMGradientProcedureFactory.create(y, f123, type);
LVMGradientProcedure p12b3 = LVMGradientProcedureFactory.create(y, PrecomputedGradient1Function.wrapGradient1Function(f12, b), type);
// This may be different
LVMGradientProcedure p12m3 = LVMGradientProcedureFactory.create(y_b, f12, type);
// Check they are the same
p123.gradient(a3peaks);
double[][] m123 = p123.getAlphaMatrix();
p12b3.gradient(a2peaks);
double s = p12b3.value;
double[] beta = p12b3.beta.clone();
double[][] alpha = p12b3.getAlphaMatrix();
System.out.printf("MLE=%b [%d] p12b3 %f %f\n", type, i, p123.value, s);
Assert.assertTrue("p12b3 Not same value @ " + i, eq.almostEqualRelativeOrAbsolute(p123.value, s));
Assert.assertTrue("p12b3 Not same gradient @ " + i, eq.almostEqualRelativeOrAbsolute(beta, p123.beta));
for (int j = 0; j < alpha.length; j++) Assert.assertTrue("p12b3 Not same alpha @ " + j, eq.almostEqualRelativeOrAbsolute(alpha[j], m123[j]));
// Check actual gradients are correct
for (int j = 0; j < nparams; j++) {
int k = indices[j];
double d = Precision.representableDelta(a2peaks[k], (a2peaks[k] == 0) ? 1e-3 : a2peaks[k] * delta);
a2peaks2[k] = a2peaks[k] + d;
p12b3.value(a2peaks2);
double s1 = p12b3.value;
a2peaks2[k] = a2peaks[k] - d;
p12b3.value(a2peaks2);
double s2 = p12b3.value;
a2peaks2[k] = a2peaks[k];
// Apply a factor of -2 to compute the actual gradients:
// See Numerical Recipes in C++, 2nd Ed. Equation 15.5.6 for Nonlinear Models
beta[j] *= -2;
double gradient = (s1 - s2) / (2 * d);
System.out.printf("[%d,%d] %f (%s %f+/-%f) %f ?= %f (%f)\n", i, k, s, f12.getName(k), a2peaks[k], d, beta[j], gradient, DoubleEquality.relativeError(gradient, beta[j]));
Assert.assertTrue("Not same gradient @ " + j, eq.almostEqualRelativeOrAbsolute(beta[j], gradient));
}
// Check these may be different
p12m3.gradient(a2peaks);
s = p12m3.value;
beta = p12m3.beta.clone();
alpha = p12m3.getAlphaMatrix();
System.out.printf("%s [%d] p12m3 %f %f\n", type, i, p123.value, s);
if (type != Type.LSQ) {
Assert.assertFalse("p12b3 Same value @ " + i, eq.almostEqualRelativeOrAbsolute(p123.value, s));
Assert.assertFalse("p12b3 Same gradient @ " + i, eq.almostEqualRelativeOrAbsolute(beta, p123.beta));
for (int j = 0; j < alpha.length; j++) {
//System.out.printf("%s != %s\n", Arrays.toString(alpha[j]), Arrays.toString(m123[j]));
Assert.assertFalse("p12b3 Same alpha @ " + j, eq.almostEqualRelativeOrAbsolute(alpha[j], m123[j]));
}
} else {
Assert.assertTrue("p12b3 Not same value @ " + i, eq.almostEqualRelativeOrAbsolute(p123.value, s));
Assert.assertTrue("p12b3 Not same gradient @ " + i, eq.almostEqualRelativeOrAbsolute(beta, p123.beta));
for (int j = 0; j < alpha.length; j++) Assert.assertTrue("p12b3 Not same alpha @ " + j, eq.almostEqualRelativeOrAbsolute(alpha[j], m123[j]));
}
// Check actual gradients are correct
for (int j = 0; j < nparams; j++) {
int k = indices[j];
double d = Precision.representableDelta(a2peaks[k], (a2peaks[k] == 0) ? 1e-3 : a2peaks[k] * delta);
a2peaks2[k] = a2peaks[k] + d;
p12m3.value(a2peaks2);
double s1 = p12m3.value;
a2peaks2[k] = a2peaks[k] - d;
p12m3.value(a2peaks2);
double s2 = p12m3.value;
a2peaks2[k] = a2peaks[k];
// Apply a factor of -2 to compute the actual gradients:
// See Numerical Recipes in C++, 2nd Ed. Equation 15.5.6 for Nonlinear Models
beta[j] *= -2;
double gradient = (s1 - s2) / (2 * d);
System.out.printf("[%d,%d] %f (%s %f+/-%f) %f ?= %f (%f)\n", i, k, s, f12.getName(k), a2peaks[k], d, beta[j], gradient, DoubleEquality.relativeError(gradient, beta[j]));
Assert.assertTrue("Not same gradient @ " + j, eq.almostEqualRelativeOrAbsolute(beta[j], gradient));
}
}
}
use of gdsc.smlm.function.ValueProcedure in project GDSC-SMLM by aherbert.
the class Gaussian2DFunction method computeValues.
/*
* (non-Javadoc)
*
* @see gdsc.smlm.function.ExtendedNonLinearFunction#computeValues(double[])
*/
public double[] computeValues(double[] variables) {
initialise0(variables);
final double[] values = new double[size()];
forEach(new ValueProcedure() {
int i = 0;
public void execute(double value) {
values[i++] = value;
}
});
return values;
}
use of gdsc.smlm.function.ValueProcedure in project GDSC-SMLM by aherbert.
the class FastMLEGradient2ProcedureTest method gradientProcedureComputesSameWithPrecomputed.
@Test
public void gradientProcedureComputesSameWithPrecomputed() {
int iter = 10;
rdg = new RandomDataGenerator(new Well19937c(30051977));
ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
ErfGaussian2DFunction f2 = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(2, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
double[] a1 = new double[7];
double[] a2 = new double[13];
final double[] x = new double[f1.size()];
final double[] b = new double[f1.size()];
for (int i = 0; i < iter; i++) {
a2[Gaussian2DFunction.BACKGROUND] = rdg.nextUniform(0.1, 0.3);
a2[Gaussian2DFunction.SIGNAL] = rdg.nextUniform(100, 300);
a2[Gaussian2DFunction.X_POSITION] = rdg.nextUniform(3, 5);
a2[Gaussian2DFunction.Y_POSITION] = rdg.nextUniform(3, 5);
a2[Gaussian2DFunction.X_SD] = rdg.nextUniform(1, 1.3);
a2[Gaussian2DFunction.Y_SD] = rdg.nextUniform(1, 1.3);
a2[6 + Gaussian2DFunction.SIGNAL] = rdg.nextUniform(100, 300);
a2[6 + Gaussian2DFunction.X_POSITION] = rdg.nextUniform(5, 7);
a2[6 + Gaussian2DFunction.Y_POSITION] = rdg.nextUniform(5, 7);
a2[6 + Gaussian2DFunction.X_SD] = rdg.nextUniform(1, 1.3);
a2[6 + Gaussian2DFunction.Y_SD] = rdg.nextUniform(1, 1.3);
// Simulate Poisson data
f2.initialise0(a2);
f1.forEach(new ValueProcedure() {
int k = 0;
public void execute(double value) {
x[k++] = (value > 0) ? rdg.nextPoisson(value) : 0;
}
});
// Precompute peak 2 (no background)
a1[Gaussian2DFunction.BACKGROUND] = 0;
for (int j = 1; j < 7; j++) a1[j] = a2[6 + j];
f1.initialise0(a1);
f1.forEach(new ValueProcedure() {
int k = 0;
public void execute(double value) {
b[k++] = value;
}
});
// Reset to peak 1
for (int j = 0; j < 7; j++) a1[j] = a2[j];
// Compute peak 1+2
FastMLEGradient2Procedure p12 = FastMLEGradient2ProcedureFactory.create(x, f2);
p12.computeSecondDerivative(a2);
double[] d11 = Arrays.copyOf(p12.d1, f1.getNumberOfGradients());
double[] d21 = Arrays.copyOf(p12.d2, f1.getNumberOfGradients());
// Compute peak 1+(precomputed 2)
FastMLEGradient2Procedure p1b2 = FastMLEGradient2ProcedureFactory.create(x, PrecomputedGradient2Function.wrapGradient2Function(f1, b));
p1b2.computeSecondDerivative(a1);
double[] d12 = p1b2.d1;
double[] d22 = p1b2.d2;
Assert.assertArrayEquals(" Result: Not same @ " + i, p12.u, p1b2.u, 1e-10);
Assert.assertArrayEquals(" D1: Not same @ " + i, d11, d12, 1e-10);
Assert.assertArrayEquals(" D2: Not same @ " + i, d21, d22, 1e-10);
double[] v1 = p12.computeValue(a2);
double[] v2 = p1b2.computeValue(a1);
Assert.assertArrayEquals(" Value: Not same @ " + i, v1, v2, 1e-10);
double[] d1 = Arrays.copyOf(p12.computeFirstDerivative(a2), f1.getNumberOfGradients());
double[] d2 = p1b2.computeFirstDerivative(a1);
Assert.assertArrayEquals(" 1st derivative: Not same @ " + i, d1, d2, 1e-10);
}
}
Aggregations