use of org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler in project GDSC-SMLM by aherbert.
the class FrcTest method canComputeMirrored.
@SeededTest
void canComputeMirrored(RandomSeed seed) {
// Sample lines through an image to create a structure.
final int size = 1024;
final double[][] data = new double[size * 2][];
final UniformRandomProvider r = RngUtils.create(seed.getSeed());
final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(r, 0, 5);
for (int x = 0, y = 0, y2 = size, i = 0; x < size; x++, y++, y2--) {
data[i++] = new double[] { x + gs.sample(), y + gs.sample() };
data[i++] = new double[] { x + gs.sample(), y2 + gs.sample() };
}
// Create 2 images
final Rectangle bounds = new Rectangle(0, 0, size, size);
ImageJImagePeakResults i1 = createImage(bounds);
ImageJImagePeakResults i2 = createImage(bounds);
final int[] indices = SimpleArrayUtils.natural(data.length);
PermutationSampler.shuffle(r, indices);
for (final int i : indices) {
final ImageJImagePeakResults image = i1;
i1 = i2;
i2 = image;
image.add((float) data[i][0], (float) data[i][1], 1);
}
i1.end();
i2.end();
final ImageProcessor ip1 = i1.getImagePlus().getProcessor();
final ImageProcessor ip2 = i2.getImagePlus().getProcessor();
// Test
final Frc frc = new Frc();
FloatProcessor[] fft1;
FloatProcessor[] fft2;
fft1 = frc.getComplexFft(ip1);
fft2 = frc.getComplexFft(ip2);
final float[] dataA1 = (float[]) fft1[0].getPixels();
final float[] dataB1 = (float[]) fft1[1].getPixels();
final float[] dataA2 = (float[]) fft2[0].getPixels();
final float[] dataB2 = (float[]) fft2[1].getPixels();
final float[] numeratorE = new float[dataA1.length];
final float[] absFft1E = new float[dataA1.length];
final float[] absFft2E = new float[dataA1.length];
Frc.compute(numeratorE, absFft1E, absFft2E, dataA1, dataB1, dataA2, dataB2);
Assertions.assertTrue(checkSymmetry(numeratorE, size), "numeratorE");
Assertions.assertTrue(checkSymmetry(absFft1E, size), "absFft1E");
Assertions.assertTrue(checkSymmetry(absFft2E, size), "absFft2E");
final float[] numeratorA = new float[dataA1.length];
final float[] absFft1A = new float[dataA1.length];
final float[] absFft2A = new float[dataA1.length];
Frc.computeMirrored(size, numeratorA, absFft1A, absFft2A, dataA1, dataB1, dataA2, dataB2);
// for (int y=0, i=0; y<size; y++)
// for (int x=0; x<size; x++, i++)
// {
// logger.fine(FunctionUtils.getSupplier("[%d,%d = %d] %f ?= %f", x, y, i, numeratorE[i],
// numeratorA[i]);
// }
Assertions.assertArrayEquals(numeratorE, numeratorA, "numerator");
Assertions.assertArrayEquals(absFft1E, absFft1A, "absFft1");
Assertions.assertArrayEquals(absFft2E, absFft2A, "absFft2");
Frc.computeMirroredFast(size, numeratorA, absFft1A, absFft2A, dataA1, dataB1, dataA2, dataB2);
// Check this.
for (int y = 1; y < size; y++) {
for (int x = 1, i = y * size + 1; x < size; x++, i++) {
Assertions.assertEquals(numeratorE[i], numeratorA[i], "numerator");
Assertions.assertEquals(absFft1E[i], absFft1A[i], "absFft1");
Assertions.assertEquals(absFft2E[i], absFft2A[i], "absFft2");
}
}
}
use of org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler in project GDSC-SMLM by aherbert.
the class FrcTest method computeMirroredIsFaster.
@SeededTest
void computeMirroredIsFaster(RandomSeed seed) {
Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
// Sample lines through an image to create a structure.
final int N = 2048;
final double[][] data = new double[N * 2][];
final UniformRandomProvider r = RngUtils.create(seed.getSeed());
final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(r, 0, 5);
for (int x = 0, y = 0, y2 = N, i = 0; x < N; x++, y++, y2--) {
data[i++] = new double[] { x + gs.sample(), y + gs.sample() };
data[i++] = new double[] { x + gs.sample(), y2 + gs.sample() };
}
// Create 2 images
final Rectangle bounds = new Rectangle(0, 0, N, N);
ImageJImagePeakResults i1 = createImage(bounds);
ImageJImagePeakResults i2 = createImage(bounds);
final int[] indices = SimpleArrayUtils.natural(data.length);
PermutationSampler.shuffle(r, indices);
for (final int i : indices) {
final ImageJImagePeakResults image = i1;
i1 = i2;
i2 = image;
image.add((float) data[i][0], (float) data[i][1], 1);
}
i1.end();
i2.end();
final ImageProcessor ip1 = i1.getImagePlus().getProcessor();
final ImageProcessor ip2 = i2.getImagePlus().getProcessor();
// Test
final Frc frc = new Frc();
FloatProcessor[] fft1;
FloatProcessor[] fft2;
fft1 = frc.getComplexFft(ip1);
fft2 = frc.getComplexFft(ip2);
final float[] dataA1 = (float[]) fft1[0].getPixels();
final float[] dataB1 = (float[]) fft1[1].getPixels();
final float[] dataA2 = (float[]) fft2[0].getPixels();
final float[] dataB2 = (float[]) fft2[1].getPixels();
final float[] numerator = new float[dataA1.length];
final float[] absFft1 = new float[dataA1.length];
final float[] absFft2 = new float[dataA1.length];
final TimingService ts = new TimingService(10);
ts.execute(new MyTimingTask("compute") {
@Override
public Object run(Object data) {
Frc.compute(numerator, absFft1, absFft2, dataA1, dataB1, dataA2, dataB2);
return null;
}
});
ts.execute(new MyTimingTask("computeMirrored") {
@Override
public Object run(Object data) {
Frc.computeMirrored(N, numerator, absFft1, absFft2, dataA1, dataB1, dataA2, dataB2);
return null;
}
});
ts.execute(new MyTimingTask("computeMirroredFast") {
@Override
public Object run(Object data) {
Frc.computeMirroredFast(N, numerator, absFft1, absFft2, dataA1, dataB1, dataA2, dataB2);
return null;
}
});
final int size = ts.getSize();
ts.repeat(size);
if (logger.isLoggable(Level.INFO)) {
logger.info(ts.getReport(size));
}
// The 'Fast' method may not always be faster so just log results
final TimingResult slow = ts.get(-3);
final TimingResult fast = ts.get(-2);
final TimingResult fastest = ts.get(-1);
logger.log(TestLogUtils.getTimingRecord(slow, fastest));
logger.log(TestLogUtils.getTimingRecord(fast, fastest));
// It should be faster than the non mirrored version
Assertions.assertTrue(ts.get(-1).getMean() <= ts.get(-3).getMean());
}
use of org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler in project GDSC-SMLM by aherbert.
the class LvmGradientProcedureTest method gradientProcedureSupportsPrecomputed.
private void gradientProcedureSupportsPrecomputed(RandomSeed seed, final Type type, boolean checkGradients) {
final int iter = 10;
final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(rng, 0, noise);
final ArrayList<double[]> paramsList = new ArrayList<>(iter);
final ArrayList<double[]> yList = new ArrayList<>(iter);
// 3 peaks
createData(rng, 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
final double min = MathUtils.min(y);
for (int j = 0; j < y.length; j++) {
y[j] = y[i] - min + gs.sample();
}
}
// 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
final Gaussian2DFunction f123 = GaussianFunctionFactory.create2D(3, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
final Gaussian2DFunction f12 = GaussianFunctionFactory.create2D(2, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
final Gaussian2DFunction f3 = GaussianFunctionFactory.create2D(1, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
final FastLog fastLog = type == Type.FAST_LOG_MLE ? getFastLog() : null;
final int nparams = f12.getNumberOfGradients();
final int[] indices = f12.gradientIndices();
final double[] b = new double[f12.size()];
// for checking strict equivalence
final DoubleEquality eq = new DoubleEquality(1e-8, 1e-16);
// for the gradients
final double delta = 1e-4;
final DoubleEquality eq2 = new DoubleEquality(5e-2, 1e-16);
final double[] a1peaks = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
final double[] y_b = new double[b.length];
// Count the number of failures for each gradient
final int failureLimit = TestCounter.computeFailureLimit(iter, 0.1);
final TestCounter failCounter = new TestCounter(failureLimit, nparams * 2);
for (int i = 0; i < paramsList.size(); i++) {
final int ii = i;
final double[] y = yList.get(i);
final double[] a3peaks = paramsList.get(i);
// logger.fine(FunctionUtils.getSupplier("[%d] a=%s", i, Arrays.toString(a3peaks));
final double[] a2peaks = Arrays.copyOf(a3peaks, 1 + 2 * Gaussian2DFunction.PARAMETERS_PER_PEAK);
final double[] a2peaks2 = a2peaks.clone();
for (int j = 1; j < a1peaks.length; j++) {
a1peaks[j] = a3peaks[j + 2 * Gaussian2DFunction.PARAMETERS_PER_PEAK];
}
// 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 index = 0;
@Override
public void execute(double value) {
b[index] = value;
// Remove negatives for MLE
if (type.isMle()) {
y[index] = Math.max(0, y[index]);
y_b[index] = Math.max(0, y[index] - value);
} else {
y_b[index] = y[index] - value;
}
index++;
}
});
final LvmGradientProcedure p123 = LvmGradientProcedureUtils.create(y, f123, type, fastLog);
// ///////////////////////////////////
// These should be the same
// ///////////////////////////////////
final LvmGradientProcedure p12b3 = LvmGradientProcedureUtils.create(y, OffsetGradient1Function.wrapGradient1Function(f12, b), type, fastLog);
// Check they are the same
p123.gradient(a3peaks);
final double[][] m123 = p123.getAlphaMatrix();
p12b3.gradient(a2peaks);
double value = p12b3.value;
final double[] beta = p12b3.beta.clone();
double[][] alpha = p12b3.getAlphaMatrix();
if (!eq.almostEqualRelativeOrAbsolute(p123.value, value)) {
Assertions.fail(FunctionUtils.getSupplier("p12b3 Not same value @ %d (error=%s) : %s == %s", i, DoubleEquality.relativeError(p123.value, value), p123.value, value));
}
if (!almostEqualRelativeOrAbsolute(eq, beta, p123.beta)) {
Assertions.fail(FunctionUtils.getSupplier("p12b3 Not same gradient @ %d (error=%s) : %s vs %s", i, relativeError(beta, p123.beta), Arrays.toString(beta), Arrays.toString(p123.beta)));
}
for (int j = 0; j < alpha.length; j++) {
// Arrays.toString(m123[j]));
if (!almostEqualRelativeOrAbsolute(eq, alpha[j], m123[j])) {
Assertions.fail(FunctionUtils.getSupplier("p12b3 Not same alpha @ %d,%d (error=%s) : %s vs %s", i, j, relativeError(alpha[j], m123[j]), Arrays.toString(alpha[j]), Arrays.toString(m123[j])));
}
}
// Check actual gradients are correct
if (checkGradients) {
for (int j = 0; j < nparams; j++) {
final int jj = j;
final int k = indices[j];
// double d = Precision.representableDelta(a2peaks[k], (a2peaks[k] == 0) ? 1e-3 :
// a2peaks[k] * delta);
final double d = Precision.representableDelta(a2peaks[k], delta);
a2peaks2[k] = a2peaks[k] + d;
p12b3.value(a2peaks2);
final double s1 = p12b3.value;
a2peaks2[k] = a2peaks[k] - d;
p12b3.value(a2peaks2);
final 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;
final double gradient = (s1 - s2) / (2 * d);
// logger.fine(FunctionUtils.getSupplier("[%d,%d] %f (%s %f+/-%f) %f ?= %f (%f)", i, k, s,
// Gaussian2DFunction.getName(k), a2peaks[k], d, beta[j], gradient,
// DoubleEquality.relativeError(gradient, beta[j]));
failCounter.run(j, () -> eq2.almostEqualRelativeOrAbsolute(beta[jj], gradient), () -> {
Assertions.fail(() -> String.format("Not same gradient @ %d,%d: %s != %s (error=%s)", ii, jj, beta[jj], gradient, DoubleEquality.relativeError(beta[jj], gradient)));
});
}
}
// ///////////////////////////////////
// This may be different
// ///////////////////////////////////
final LvmGradientProcedure p12m3 = LvmGradientProcedureUtils.create(y_b, f12, type, fastLog);
// Check these may be different.
// Sometimes they are not different.
p12m3.gradient(a2peaks);
value = p12m3.value;
System.arraycopy(p12m3.beta, 0, beta, 0, p12m3.beta.length);
alpha = p12m3.getAlphaMatrix();
if (type != Type.LSQ) {
if (eq.almostEqualRelativeOrAbsolute(p123.value, value)) {
logger.log(TestLogUtils.getFailRecord("p12b3 Same value @ %d (error=%s) : %s == %s", i, DoubleEquality.relativeError(p123.value, value), p123.value, value));
}
if (almostEqualRelativeOrAbsolute(eq, beta, p123.beta)) {
logger.log(TestLogUtils.getFailRecord("p12b3 Same gradient @ %d (error=%s) : %s vs %s", i, relativeError(beta, p123.beta), Arrays.toString(beta), Arrays.toString(p123.beta)));
}
// Note: Test the matrix is different by finding 1 different column
int dj = -1;
for (int j = 0; j < alpha.length; j++) {
// Arrays.toString(m123[j]));
if (!almostEqualRelativeOrAbsolute(eq, alpha[j], m123[j])) {
// Different column
dj = j;
break;
}
}
if (dj == -1) {
// Find biggest error for reporting. This helps set the test tolerance.
double error = 0;
dj = -1;
for (int j = 0; j < alpha.length; j++) {
final double e = relativeError(alpha[j], m123[j]);
if (error <= e) {
error = e;
dj = j;
}
}
logger.log(TestLogUtils.getFailRecord("p12b3 Same alpha @ %d,%d (error=%s) : %s vs %s", i, dj, error, Arrays.toString(alpha[dj]), Arrays.toString(m123[dj])));
}
} else {
if (!eq.almostEqualRelativeOrAbsolute(p123.value, value)) {
logger.log(TestLogUtils.getFailRecord("p12b3 Not same value @ %d (error=%s) : %s == %s", i, DoubleEquality.relativeError(p123.value, value), p123.value, value));
}
if (!almostEqualRelativeOrAbsolute(eq, beta, p123.beta)) {
logger.log(TestLogUtils.getFailRecord("p12b3 Not same gradient @ %d (error=%s) : %s vs %s", i, relativeError(beta, p123.beta), Arrays.toString(beta), Arrays.toString(p123.beta)));
}
for (int j = 0; j < alpha.length; j++) {
// Arrays.toString(m123[j]));
if (!almostEqualRelativeOrAbsolute(eq, alpha[j], m123[j])) {
logger.log(TestLogUtils.getFailRecord("p12b3 Not same alpha @ %d,%d (error=%s) : %s vs %s", i, j, relativeError(alpha[j], m123[j]), Arrays.toString(alpha[j]), Arrays.toString(m123[j])));
}
}
}
// Check actual gradients are correct
if (!checkGradients) {
continue;
}
for (int j = 0; j < nparams; j++) {
final int jj = j;
final int k = indices[j];
// double d = Precision.representableDelta(a2peaks[k], (a2peaks[k] == 0) ? 1e-3 : a2peaks[k]
// * delta);
final double d = Precision.representableDelta(a2peaks[k], delta);
a2peaks2[k] = a2peaks[k] + d;
p12m3.value(a2peaks2);
final double s1 = p12m3.value;
a2peaks2[k] = a2peaks[k] - d;
p12m3.value(a2peaks2);
final 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;
final double gradient = (s1 - s2) / (2 * d);
// logger.fine(FunctionUtils.getSupplier("[%d,%d] %f (%s %f+/-%f) %f ?= %f (%f)", i, k, s,
// Gaussian2DFunction.getName(k), a2peaks[k], d, beta[j], gradient,
// DoubleEquality.relativeError(gradient, beta[j]));
failCounter.run(nparams + j, () -> eq2.almostEqualRelativeOrAbsolute(beta[jj], gradient), () -> {
Assertions.fail(() -> String.format("Not same gradient @ %d,%d: %s != %s (error=%s)", ii, jj, beta[jj], gradient, DoubleEquality.relativeError(beta[jj], gradient)));
});
}
}
}
use of org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler in project GDSC-SMLM by aherbert.
the class BaseFunctionSolverTest method createData.
private static double[][] createData(RandomSeed source) {
// Per observation read noise.
// This is generated once so create the randon generator here.
final UniformRandomProvider rg = RngUtils.create(source.getSeed());
final ContinuousSampler ed = SamplerUtils.createExponentialSampler(rg, variance);
final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(rg, gain, gainSD);
final double[] w = new double[size * size];
final double[] n = new double[size * size];
for (int i = 0; i < w.length; i++) {
final double pixelVariance = ed.sample();
final double pixelGain = Math.max(0.1, gs.sample());
// weights = var / g^2
w[i] = pixelVariance / (pixelGain * pixelGain);
// Convert to standard deviation for simulation
n[i] = Math.sqrt(w[i]);
}
return new double[][] { w, n };
}
use of org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler in project GDSC-SMLM by aherbert.
the class WeightedKernelFilterTest method filterPerformsWeightedKernelFiltering.
@SeededTest
void filterPerformsWeightedKernelFiltering(RandomSeed seed) {
final DataFilter filter = createDataFilter();
final UniformRandomProvider rg = RngUtils.create(seed.getSeed());
final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(rg, 2, 0.2);
final float[] offsets = getOffsets(filter);
final int[] boxSizes = getBoxSizes(filter);
final TDoubleArrayList l1 = new TDoubleArrayList();
final FloatFloatBiPredicate equality = TestHelper.floatsAreClose(1e-6, 0);
for (final int width : primes) {
for (final int height : new int[] { 29 }) {
final float[] data = createData(width, height, rg);
l1.reset();
// Uniform weights
final float[] w1 = new float[width * height];
Arrays.fill(w1, 0.5f);
// Weights simulating the variance of sCMOS pixels
final float[] w2 = new float[width * height];
for (int i = 0; i < w2.length; i++) {
w2[i] = (float) (1.0 / Math.max(0.01, gs.sample()));
}
for (final int boxSize : boxSizes) {
for (final float offset : offsets) {
for (final boolean internal : checkInternal) {
// For each pixel over the range around the pixel (vi).
// kernel filter: sum(vi * ki) / sum(ki)
// Weighted kernel filter: sum(vi * wi * ki) / sum(ki * wi)
// Note: The kernel filter is like a weighted filter
// (New kernel = wi * ki)
filter.setWeights(null, width, height);
// Uniform weights
testfilterPerformsWeightedFiltering(filter, width, height, data, w1, boxSize, offset, internal, equality);
// Random weights.
testfilterPerformsWeightedFiltering(filter, width, height, data, w2, boxSize, offset, internal, equality);
}
}
}
}
}
}
Aggregations