use of org.apache.commons.math3.random.Well19937c 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 org.apache.commons.math3.random.Well19937c in project GDSC-SMLM by aherbert.
the class PoissonGammaGaussianFunctionTest method fasterThan.
private void fasterThan(PoissonGammaGaussianFunction f1, PoissonGammaGaussianFunction f2) {
// Generate realistic data from the probability mass function
double[][] samples = new double[photons.length][];
for (int j = 0; j < photons.length; j++) {
int start = (int) (4 * -s);
int u = start;
StoredDataStatistics stats = new StoredDataStatistics();
while (stats.getSum() < 0.995) {
stats.add(f1.likelihood(u, photons[j]));
u++;
}
// Generate cumulative probability
double[] data = stats.getValues();
for (int i = 1; i < data.length; i++) data[i] += data[i - 1];
// Sample
RandomGenerator rand = new Well19937c();
double[] sample = new double[1000];
for (int i = 0; i < sample.length; i++) {
final double p = rand.nextDouble();
int x = 0;
while (x < data.length && data[x] < p) x++;
sample[i] = x;
}
samples[j] = sample;
}
// Warm-up
run(f1, samples, photons);
run(f2, samples, photons);
long t1 = 0;
for (int i = 0; i < 5; i++) t1 += run(f1, samples, photons);
long t2 = 0;
for (int i = 0; i < 5; i++) t2 += run(f2, samples, photons);
System.out.printf("%s %d -> %s %d = %f x\n", getName(f1), t1, getName(f2), t2, (double) t1 / t2);
}
use of org.apache.commons.math3.random.Well19937c in project GDSC-SMLM by aherbert.
the class FastMLEJacobianGradient2ProcedureTest method gradientProcedureComputesSameAsBaseGradientProcedure.
private void gradientProcedureComputesSameAsBaseGradientProcedure(int nparams) {
int iter = 10;
rdg = new RandomDataGenerator(new Well19937c(30051977));
ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
ArrayList<double[]> yList = new ArrayList<double[]>(iter);
createFakeData(nparams, iter, paramsList, yList);
FakeGradientFunction func = new FakeGradientFunction(blockWidth, nparams);
for (int i = 0; i < paramsList.size(); i++) {
FastMLEGradient2Procedure p = FastMLEGradient2ProcedureFactory.createUnrolled(yList.get(i), func);
FastMLEJacobianGradient2Procedure p2 = new FastMLEJacobianGradient2Procedure(yList.get(i), func);
p.computeSecondDerivative(paramsList.get(i));
p2.computeSecondDerivative(paramsList.get(i));
// Virtually the same ...
Assert.assertArrayEquals(p.d1, p2.d1, 1e-5);
Assert.assertArrayEquals(p.d2, p2.d2, 1e-5);
}
}
use of org.apache.commons.math3.random.Well19937c in project GDSC-SMLM by aherbert.
the class GradientCalculatorSpeedTest method gradientCalculatorComputesSameOutputWithBias.
@Test
public void gradientCalculatorComputesSameOutputWithBias() {
Gaussian2DFunction func = new SingleEllipticalGaussian2DFunction(blockWidth, blockWidth);
int nparams = func.getNumberOfGradients();
GradientCalculator calc = new GradientCalculator(nparams);
int n = func.size();
int iter = 100;
rdg = new RandomDataGenerator(new Well19937c(30051977));
ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
ArrayList<double[]> yList = new ArrayList<double[]>(iter);
ArrayList<double[][]> alphaList = new ArrayList<double[][]>(iter);
ArrayList<double[]> betaList = new ArrayList<double[]>(iter);
ArrayList<double[]> xList = new ArrayList<double[]>(iter);
// Manipulate the background
double defaultBackground = Background;
try {
Background = 1e-2;
createData(1, iter, paramsList, yList, true);
EJMLLinearSolver solver = new EJMLLinearSolver(1e-5, 1e-6);
for (int i = 0; i < paramsList.size(); i++) {
double[] y = yList.get(i);
double[] a = paramsList.get(i);
double[][] alpha = new double[nparams][nparams];
double[] beta = new double[nparams];
calc.findLinearised(n, y, a, alpha, beta, func);
alphaList.add(alpha);
betaList.add(beta.clone());
for (int j = 0; j < nparams; j++) {
if (Math.abs(beta[j]) < 1e-6)
System.out.printf("[%d] Tiny beta %s %g\n", i, func.getName(j), beta[j]);
}
// Solve
if (!solver.solve(alpha, beta))
throw new AssertionError();
xList.add(beta);
//System.out.println(Arrays.toString(beta));
}
double[][] alpha = new double[nparams][nparams];
double[] beta = new double[nparams];
//for (int b = 1; b < 1000; b *= 2)
for (double b : new double[] { -500, -100, -10, -1, -0.1, 0, 0.1, 1, 10, 100, 500 }) {
Statistics[] rel = new Statistics[nparams];
Statistics[] abs = new Statistics[nparams];
for (int i = 0; i < nparams; i++) {
rel[i] = new Statistics();
abs[i] = new Statistics();
}
for (int i = 0; i < paramsList.size(); i++) {
double[] y = add(yList.get(i), b);
double[] a = paramsList.get(i).clone();
a[0] += b;
calc.findLinearised(n, y, a, alpha, beta, func);
double[][] alpha2 = alphaList.get(i);
double[] beta2 = betaList.get(i);
double[] x2 = xList.get(i);
Assert.assertArrayEquals("Beta", beta2, beta, 1e-10);
for (int j = 0; j < nparams; j++) {
Assert.assertArrayEquals("Alpha", alpha2[j], alpha[j], 1e-10);
}
// Solve
solver.solve(alpha, beta);
Assert.assertArrayEquals("X", x2, beta, 1e-10);
for (int j = 0; j < nparams; j++) {
rel[j].add(DoubleEquality.relativeError(x2[j], beta[j]));
abs[j].add(Math.abs(x2[j] - beta[j]));
}
}
for (int i = 0; i < nparams; i++) System.out.printf("Bias = %.2f : %s : Rel %g +/- %g: Abs %g +/- %g\n", b, func.getName(i), rel[i].getMean(), rel[i].getStandardDeviation(), abs[i].getMean(), abs[i].getStandardDeviation());
}
} finally {
Background = defaultBackground;
}
}
use of org.apache.commons.math3.random.Well19937c in project GDSC-SMLM by aherbert.
the class PulseActivationAnalysisTest method canLinearlyUnmix3Channels.
private void canLinearlyUnmix3Channels(int n, int m) {
RandomGenerator r = new Well19937c(30051977);
try {
for (int loop = 0; loop < 10; loop++) {
// A rough mix of each channel
double[] d = create(6, r, 100, 100);
// Total crosstalk per channel should be below 50%
double[] c = create(6, r, 0, 0.25);
// Enumerate
Iterator<int[]> it = CombinatoricsUtils.combinationsIterator(3, n);
while (it.hasNext()) {
final int[] channels = it.next();
double[] dd = new double[3];
for (int i : channels) dd[i] = d[i];
Iterator<int[]> it2 = CombinatoricsUtils.combinationsIterator(6, m);
while (it2.hasNext()) {
final int[] crosstalk = it2.next();
double[] cc = new double[6];
for (int i : crosstalk) cc[i] = c[i];
canLinearlyUnmix3Channels(dd[0], dd[1], dd[2], cc[0], cc[1], cc[2], cc[3], cc[4], cc[5]);
}
}
}
} catch (AssertionError e) {
throw new AssertionError(String.format("channels=%d, crosstalk=%d", n, m), e);
}
}
Aggregations