use of org.ejml.data.DenseMatrix64F in project GDSC-SMLM by aherbert.
the class EJMLLinearSolver method solve.
/**
* Solves (one) linear equation, a x = b
* <p>
* On output b replaced by x. Matrix a may be modified.
*
* @return False if the equation is singular (no solution)
*/
private boolean solve(LinearSolver<DenseMatrix64F> solver, DenseMatrix64F A, DenseMatrix64F B) {
boolean copy = (errorChecking || solver.modifiesB());
DenseMatrix64F x = (copy) ? getX() : B;
if (!solve(solver, A, B, x))
return false;
// Copy back result if necessary
if (copy)
System.arraycopy(x.data, 0, B.data, 0, B.numRows);
return true;
}
use of org.ejml.data.DenseMatrix64F 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 org.ejml.data.DenseMatrix64F in project GDSC-SMLM by aherbert.
the class ChiSquaredDistributionTableTest method canComputeChiSquared.
@Test
public void canComputeChiSquared() {
// We have to use the transpose of the table
DenseMatrix64F m = new DenseMatrix64F(chi2);
CommonOps.transpose(m);
int max = m.numCols;
double[] et = m.data;
for (int i = 0, j = 0; i < p.length; i++) {
ChiSquaredDistributionTable upperTable = ChiSquaredDistributionTable.createUpperTailed(p[i], max);
// Use 1-p as the significance level to get the same critical values
ChiSquaredDistributionTable lowerTable = ChiSquaredDistributionTable.createLowerTailed(1 - p[i], max);
for (int df = 1; df <= max; df++) {
double o = upperTable.getCrititalValue(df);
double e = et[j++];
//System.out.printf("p=%.3f,df=%d = %f\n", p[i], df, o);
Assert.assertEquals(e, o, 1e-2);
// The test only stores 2 decimal places so use the computed value to set upper/lower
double upper = o * 1.01;
double lower = o * 0.99;
Assert.assertTrue("Upper did not reject higher", upperTable.reject(upper, df));
Assert.assertFalse("Upper did not reject actual value", upperTable.reject(o, df));
Assert.assertFalse("Upper did not accept lower", upperTable.reject(lower, df));
Assert.assertTrue("Lower did not reject lower", lowerTable.reject(lower, df));
Assert.assertFalse("Lower did not accept actual value", lowerTable.reject(o, df));
Assert.assertFalse("Lower did not accept higher", lowerTable.reject(upper, df));
}
}
}
use of org.ejml.data.DenseMatrix64F in project GDSC-SMLM by aherbert.
the class LSQLVMGradientProcedureTest method gradientProcedureComputesSameAsGradientCalculator.
private void gradientProcedureComputesSameAsGradientCalculator(int nparams, BaseLSQLVMGradientProcedureFactory factory) {
int iter = 10;
rdg = new RandomDataGenerator(new Well19937c(30051977));
double[][] alpha = new double[nparams][nparams];
double[] beta = new double[nparams];
ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
ArrayList<double[]> yList = new ArrayList<double[]>(iter);
int[] x = createFakeData(nparams, iter, paramsList, yList);
int n = x.length;
FakeGradientFunction func = new FakeGradientFunction(blockWidth, nparams);
GradientCalculator calc = GradientCalculatorFactory.newCalculator(nparams, false);
String name = factory.getClass().getSimpleName();
for (int i = 0; i < paramsList.size(); i++) {
BaseLSQLVMGradientProcedure p = factory.createProcedure(yList.get(i), func);
p.gradient(paramsList.get(i));
double s = p.value;
double s2 = calc.findLinearised(n, yList.get(i), paramsList.get(i), alpha, beta, func);
// Exactly the same ...
Assert.assertEquals(name + " Result: Not same @ " + i, s, s2, 0);
Assert.assertArrayEquals(name + " Observations: Not same beta @ " + i, p.beta, beta, 0);
double[] al = p.getAlphaLinear();
Assert.assertArrayEquals(name + " Observations: Not same alpha @ " + i, al, new DenseMatrix64F(alpha).data, 0);
double[][] am = p.getAlphaMatrix();
for (int j = 0; j < nparams; j++) Assert.assertArrayEquals(name + " Observations: Not same alpha @ " + i, am[j], alpha[j], 0);
}
}
use of org.ejml.data.DenseMatrix64F in project GDSC-SMLM by aherbert.
the class LVMGradientProcedureTest method gradientProcedureComputesSameAsGradientCalculator.
private void gradientProcedureComputesSameAsGradientCalculator(int nparams, boolean mle) {
int iter = 10;
rdg = new RandomDataGenerator(new Well19937c(30051977));
double[][] alpha = new double[nparams][nparams];
double[] beta = new double[nparams];
ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
ArrayList<double[]> yList = new ArrayList<double[]>(iter);
int[] x = createFakeData(nparams, iter, paramsList, yList);
int n = x.length;
FakeGradientFunction func = new FakeGradientFunction(blockWidth, nparams);
GradientCalculator calc = GradientCalculatorFactory.newCalculator(nparams, mle);
String name = String.format("[%d] %b", nparams, mle);
for (int i = 0; i < paramsList.size(); i++) {
LVMGradientProcedure p = LVMGradientProcedureFactory.create(yList.get(i), func, (mle) ? Type.MLE : Type.LSQ);
p.gradient(paramsList.get(i));
double s = p.value;
double s2 = calc.findLinearised(n, yList.get(i), paramsList.get(i), alpha, beta, func);
// Exactly the same ...
Assert.assertEquals(name + " Result: Not same @ " + i, s, s2, 0);
Assert.assertArrayEquals(name + " Observations: Not same beta @ " + i, p.beta, beta, 0);
double[] al = p.getAlphaLinear();
Assert.assertArrayEquals(name + " Observations: Not same alpha @ " + i, al, new DenseMatrix64F(alpha).data, 0);
double[][] am = p.getAlphaMatrix();
for (int j = 0; j < nparams; j++) Assert.assertArrayEquals(name + " Observations: Not same alpha @ " + i, am[j], alpha[j], 0);
}
}
Aggregations