use of uk.ac.sussex.gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator in project GDSC-SMLM by aherbert.
the class SolverSpeedTest method createSolverData.
private static boolean createSolverData(UniformRandomProvider rand, float[][] alpha, float[] beta, boolean positiveDifinite) {
// Generate a 2D Gaussian
final SingleFreeCircularGaussian2DFunction func = new SingleFreeCircularGaussian2DFunction(10, 10);
final double[] params = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
params[Gaussian2DFunction.BACKGROUND] = 2 + rand.nextDouble() * 2;
params[Gaussian2DFunction.SIGNAL] = 100 + rand.nextDouble() * 5;
params[Gaussian2DFunction.X_POSITION] = 4.5 + rand.nextDouble();
params[Gaussian2DFunction.Y_POSITION] = 4.5 + rand.nextDouble();
params[Gaussian2DFunction.X_SD] = 1 + rand.nextDouble();
params[Gaussian2DFunction.Y_SD] = 1 + rand.nextDouble();
params[Gaussian2DFunction.ANGLE] = rand.nextDouble();
final int[] x = new int[100];
final double[] y = new double[100];
func.initialise(params);
for (int i = 0; i < x.length; i++) {
// Add random noise
y[i] = func.eval(i) + ((rand.nextDouble() < 0.5) ? -rand.nextDouble() * 5 : rand.nextDouble() * 5);
}
// Randomise parameters
for (int i = 0; i < params.length; i++) {
params[i] += (rand.nextDouble() < 0.5) ? -rand.nextDouble() : rand.nextDouble();
}
// Compute the Hessian and parameter gradient vector
final GradientCalculator calc = new GradientCalculator(6);
final double[][] alpha2 = new double[6][6];
final double[] beta2 = new double[6];
calc.findLinearised(y.length, y, params, alpha2, beta2, func);
// Update the Hessian using a lambda shift
final double lambda = 1.001;
for (int i = 0; i < alpha2.length; i++) {
alpha2[i][i] *= lambda;
}
// Copy back
for (int i = 0; i < beta.length; i++) {
beta[i] = (float) beta2[i];
for (int j = 0; j < beta.length; j++) {
alpha[i][j] = (float) alpha2[i][j];
}
}
// Check for a positive definite matrix
if (positiveDifinite) {
final EjmlLinearSolver solver = new EjmlLinearSolver();
return solver.solveCholeskyLdlT(copydouble(alpha), copydouble(beta));
}
return true;
}
use of uk.ac.sussex.gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator in project GDSC-SMLM by aherbert.
the class ApacheLvmFitter method computeValue.
@Override
public boolean computeValue(double[] y, double[] fx, double[] a) {
final GradientCalculator calculator = GradientCalculatorUtils.newCalculator(function.getNumberOfGradients(), false);
// Since we know the function is a Gaussian2DFunction from the constructor
value = calculator.findLinearised(y.length, y, fx, a, (NonLinearFunction) function);
return true;
}
use of uk.ac.sussex.gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator in project GDSC-SMLM by aherbert.
the class ApacheLvmFitter method computeFisherInformationMatrix.
@Override
protected FisherInformationMatrix computeFisherInformationMatrix(double[] y, double[] a) {
final GradientCalculator c = GradientCalculatorUtils.newCalculator(function.getNumberOfGradients(), false);
// Since we know the function is a Gaussian2DFunction from the constructor
final double[][] I = c.fisherInformationMatrix(y.length, a, (NonLinearFunction) function);
if (c.isNaNGradients()) {
throw new FunctionSolverException(FitStatus.INVALID_GRADIENTS);
}
return new FisherInformationMatrix(I);
}
use of uk.ac.sussex.gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator in project GDSC-SMLM by aherbert.
the class EjmlLinearSolverTest method runSolverSpeedTest.
private void runSolverSpeedTest(RandomSeed seed, int flags) {
Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
final Gaussian2DFunction f0 = GaussianFunctionFactory.create2D(1, 10, 10, flags, null);
final int n = f0.size();
final double[] y = new double[n];
final LocalList<DenseMatrix64F> aList = new LocalList<>();
final LocalList<DenseMatrix64F> bList = new LocalList<>();
final double[] testbackground = new double[] { 0.2, 0.7 };
final double[] testsignal1 = new double[] { 30, 100, 300 };
final double[] testcx1 = new double[] { 4.9, 5.3 };
final double[] testcy1 = new double[] { 4.8, 5.2 };
final double[] testw1 = new double[] { 1.1, 1.2, 1.5 };
final int np = f0.getNumberOfGradients();
final GradientCalculator calc = GradientCalculatorUtils.newCalculator(np);
final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
// double lambda = 10;
for (final double background : testbackground) {
// Peak 1
for (final double signal1 : testsignal1) {
for (final double cx1 : testcx1) {
for (final double cy1 : testcy1) {
for (final double w1 : testw1) {
final double[] p = new double[] { background, signal1, 0, cx1, cy1, w1, w1 };
f0.initialise(p);
f0.forEach(new ValueProcedure() {
int index = 0;
@Override
public void execute(double value) {
// Poisson data
y[index++] = GdscSmlmTestUtils.createPoissonSampler(rng, value).sample();
}
});
final double[][] alpha = new double[np][np];
final double[] beta = new double[np];
// double ss =
calc.findLinearised(n, y, p, alpha, beta, f0);
// TestLog.fine(logger,"SS = %f", 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));
}
}
}
}
}
final DenseMatrix64F[] a = aList.toArray(new DenseMatrix64F[0]);
final DenseMatrix64F[] b = bList.toArray(new DenseMatrix64F[0]);
final int runs = 100000 / a.length;
final TimingService ts = new TimingService(runs);
final LocalList<SolverTimingTask> tasks = new LocalList<>();
// Added in descending speed order
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 (final SolverTimingTask task : tasks) {
if (!task.badSolver) {
ts.execute(task);
}
}
final int size = ts.getSize();
ts.repeat();
if (logger.isLoggable(Level.INFO)) {
logger.info(ts.getReport(size));
}
// Just check the PseudoInverse is slowest
for (int i = 1; i < size; i++) {
logger.log(TestLogUtils.getTimingRecord(ts.get(-(size)), ts.get(-i)));
}
if (np > 2) {
// The Direct solver may not be faster at size=5
int i = (np == 5) ? 2 : 1;
final int size_1 = size - 1;
for (; i < size_1; i++) {
logger.log(TestLogUtils.getTimingRecord(ts.get(-(size_1)), ts.get(-i)));
}
}
}
use of uk.ac.sussex.gdsc.smlm.fitting.nonlinear.gradient.GradientCalculator in project GDSC-SMLM by aherbert.
the class EjmlLinearSolverTest method runInversionSpeedTest.
private void runInversionSpeedTest(RandomSeed seed, int flags) {
Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
final Gaussian2DFunction f0 = GaussianFunctionFactory.create2D(1, 10, 10, flags, null);
final int n = f0.size();
final double[] y = new double[n];
final LocalList<DenseMatrix64F> aList = new LocalList<>();
final double[] testbackground = new double[] { 0.2, 0.7 };
final double[] testsignal1 = new double[] { 30, 100, 300 };
final double[] testcx1 = new double[] { 4.9, 5.3 };
final double[] testcy1 = new double[] { 4.8, 5.2 };
final double[] testw1 = new double[] { 1.1, 1.2, 1.5 };
final int np = f0.getNumberOfGradients();
final GradientCalculator calc = GradientCalculatorUtils.newCalculator(np);
final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
// double lambda = 10;
for (final double background : testbackground) {
// Peak 1
for (final double signal1 : testsignal1) {
for (final double cx1 : testcx1) {
for (final double cy1 : testcy1) {
for (final double w1 : testw1) {
final double[] p = new double[] { background, signal1, 0, cx1, cy1, w1, w1 };
f0.initialise(p);
f0.forEach(new ValueProcedure() {
int index = 0;
@Override
public void execute(double value) {
// Poisson data
y[index++] = GdscSmlmTestUtils.createPoissonSampler(rng, value).sample();
}
});
final double[][] alpha = new double[np][np];
final double[] beta = new double[np];
// double ss =
calc.findLinearised(n, y, p, alpha, beta, f0);
// TestLog.fine(logger,"SS = %f", ss);
// As per the LVM algorithm
// for (int i = 0; i < np; i++)
// alpha[i][i] *= lambda;
aList.add(EjmlLinearSolver.toA(alpha));
}
}
}
}
}
final DenseMatrix64F[] a = aList.toArray(new DenseMatrix64F[0]);
final boolean[] ignore = new boolean[a.length];
final double[][] answer = new double[a.length][];
final int runs = 100000 / a.length;
final TimingService ts = new TimingService(runs);
final LocalList<InversionTimingTask> tasks = new LocalList<>();
// Added in descending speed order
tasks.add(new PseudoInverseInversionTimingTask(a, ignore, answer));
tasks.add(new LinearInversionTimingTask(a, ignore, answer));
tasks.add(new CholeskyLdltInversionTimingTask(a, ignore, answer));
tasks.add(new CholeskyInversionTimingTask(a, ignore, answer));
tasks.add(new DirectInversionInversionTimingTask(a, ignore, answer));
tasks.add(new DiagonalDirectInversionInversionTimingTask(a, ignore, answer));
for (final InversionTimingTask task : tasks) {
if (!task.badSolver) {
ts.execute(task);
}
}
final int size = ts.getSize();
ts.repeat();
if (logger.isLoggable(Level.INFO)) {
logger.info(ts.getReport(size));
}
// When it is present the DiagonalDirect is fastest (n<=5)
if (np <= 5) {
for (int i = 2; i <= size; i++) {
logger.log(TestLogUtils.getTimingRecord(ts.get(-i), ts.get(-1)));
}
if (np < 5) {
// n < 5 Direct is fastest
for (int i = 3; i <= size; i++) {
logger.log(TestLogUtils.getTimingRecord(ts.get(-i), ts.get(-2)));
}
} else {
// and may not be faster than Direct at n=5 so that comparison is ignored.
for (int i = 4; i <= size; i++) {
logger.log(TestLogUtils.getTimingRecord(ts.get(-i), ts.get(-3)));
}
}
} else {
// Cholesky should be fastest.
for (int i = 2; i <= size; i++) {
logger.log(TestLogUtils.getTimingRecord(ts.get(-i), ts.get(-1)));
}
}
}
Aggregations