use of org.apache.commons.math3.random.Well19937c in project GDSC-SMLM by aherbert.
the class PrecomputedFunctionTest method precomputedGradient2FunctionWrapsPrecomputedValues.
@Test
public void precomputedGradient2FunctionWrapsPrecomputedValues() {
int n = 3;
RandomGenerator r = new Well19937c(30051977);
Gradient2Function f0 = new FakeGradientFunction(3, n);
int size = f0.size();
double[] b1 = new PseudoRandomGenerator(size, r).getSequence();
double[] b2 = new PseudoRandomGenerator(size, r).getSequence();
Gradient2Function f1 = PrecomputedGradient2Function.wrapGradient2Function(f0, b1);
Gradient2Function f2 = PrecomputedGradient2Function.wrapGradient2Function(f1, b2);
double[] p = new double[n];
for (int i = 0; i < n; i++) p[i] = r.nextDouble();
double[] d0 = new double[n];
double[] d1 = new double[n];
double[] d2 = new double[n];
double[] d20 = new double[n];
double[] d21 = new double[n];
double[] d22 = new double[n];
double[] v0 = evaluateGradient2Function(f0, p, d0, d20);
double[] v1 = evaluateGradient2Function(f1, p, d1, d21);
double[] v2 = evaluateGradient2Function(f2, p, d2, d22);
for (int i = 0; i < v0.length; i++) {
double e = v0[i] + b1[i] + b2[i];
double o1 = v1[i] + b2[i];
double o2 = v2[i];
Assert.assertEquals("o1", e, o1, 0);
Assert.assertEquals("o2", e, o2, 1e-6);
}
Assert.assertArrayEquals("d1", d0, d1, 0);
Assert.assertArrayEquals("d2", d0, d2, 0);
Assert.assertArrayEquals("d21", d20, d21, 0);
Assert.assertArrayEquals("d22", d20, d22, 0);
}
use of org.apache.commons.math3.random.Well19937c in project GDSC-SMLM by aherbert.
the class ConfigurationTemplateTest method canLoadTemplateImageFromFile.
@Test
public void canLoadTemplateImageFromFile() throws IOException {
ConfigurationTemplate.clearTemplates();
Assert.assertEquals(0, ConfigurationTemplate.getTemplateNamesWithImage().length);
// Create a dummy image
int size = 20;
float[] pixels = new float[size * size];
RandomGenerator r = new Well19937c();
for (int i = pixels.length; i-- > 0; ) pixels[i] = r.nextFloat();
ImagePlus imp = new ImagePlus("test", new FloatProcessor(size, size, pixels));
File tmpFile = File.createTempFile("tmp", ".tif");
tmpFile.deleteOnExit();
IJ.save(imp, tmpFile.getPath());
String name = "canLoadTemplateImageFromFile";
File file = new File(Utils.replaceExtension(tmpFile.getPath(), ".xml"));
ConfigurationTemplate.saveTemplate(name, new GlobalSettings(), file);
Assert.assertEquals(1, ConfigurationTemplate.getTemplateNamesWithImage().length);
ImagePlus imp2 = ConfigurationTemplate.getTemplateImage(name);
Assert.assertNotNull(imp2);
float[] data = (float[]) imp2.getProcessor().toFloat(0, null).getPixels();
Assert.assertArrayEquals(pixels, data, 0);
}
use of org.apache.commons.math3.random.Well19937c in project GDSC-SMLM by aherbert.
the class PulseActivationAnalysisTest method canLinearlyUnmix2Channels.
private void canLinearlyUnmix2Channels(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(2, r, 1, 100);
// Crosstalk should be below 50%
double[] c = create(2, r, 0, 0.5);
// Enumerate
Iterator<int[]> it = CombinatoricsUtils.combinationsIterator(2, n);
while (it.hasNext()) {
final int[] channels = it.next();
double[] dd = new double[2];
for (int i : channels) dd[i] = d[i];
Iterator<int[]> it2 = CombinatoricsUtils.combinationsIterator(2, m);
while (it2.hasNext()) {
final int[] crosstalk = it2.next();
double[] cc = new double[2];
for (int i : crosstalk) cc[i] = c[i];
canLinearlyUnmix2Channels(dd[0], dd[1], cc[0], cc[1]);
}
}
}
} catch (AssertionError e) {
throw new AssertionError(String.format("channels=%d, crosstalk=%d", n, m), e);
}
}
use of org.apache.commons.math3.random.Well19937c in project GDSC-SMLM by aherbert.
the class FIRE method runQEstimation.
private void runQEstimation() {
IJ.showStatus(TITLE + " ...");
if (!showQEstimationInputDialog())
return;
MemoryPeakResults results = ResultsManager.loadInputResults(inputOption, false);
if (results == null || results.size() == 0) {
IJ.error(TITLE, "No results could be loaded");
return;
}
if (results.getCalibration() == null) {
IJ.error(TITLE, "The results are not calibrated");
return;
}
results = cropToRoi(results);
if (results.size() < 2) {
IJ.error(TITLE, "No results within the crop region");
return;
}
initialise(results, null);
// We need localisation precision.
// Build a histogram of the localisation precision.
// Get the initial mean and SD and plot as a Gaussian.
PrecisionHistogram histogram = calculatePrecisionHistogram();
if (histogram == null) {
IJ.error(TITLE, "No localisation precision available.\n \nPlease choose " + PrecisionMethod.FIXED + " and enter a precision mean and SD.");
return;
}
StoredDataStatistics precision = histogram.precision;
//String name = results.getName();
double fourierImageScale = SCALE_VALUES[imageScaleIndex];
int imageSize = IMAGE_SIZE_VALUES[imageSizeIndex];
// Create the image and compute the numerator of FRC.
// Do not use the signal so results.size() is the number of localisations.
IJ.showStatus("Computing FRC curve ...");
FireImages images = createImages(fourierImageScale, imageSize, false);
// DEBUGGING - Save the two images to disk. Load the images into the Matlab
// code that calculates the Q-estimation and make this plugin match the functionality.
//IJ.save(new ImagePlus("i1", images.ip1), "/scratch/i1.tif");
//IJ.save(new ImagePlus("i2", images.ip2), "/scratch/i2.tif");
FRC frc = new FRC();
frc.progress = progress;
frc.setFourierMethod(fourierMethod);
frc.setSamplingMethod(samplingMethod);
frc.setPerimeterSamplingFactor(perimeterSamplingFactor);
FRCCurve frcCurve = frc.calculateFrcCurve(images.ip1, images.ip2, images.nmPerPixel);
if (frcCurve == null) {
IJ.error(TITLE, "Failed to compute FRC curve");
return;
}
IJ.showStatus("Running Q-estimation ...");
// Note:
// The method implemented here is based on Matlab code provided by Bernd Rieger.
// The idea is to compute the spurious correlation component of the FRC Numerator
// using an initial estimate of distribution of the localisation precision (assumed
// to be Gaussian). This component is the contribution of repeat localisations of
// the same molecule to the numerator and is modelled as an exponential decay
// (exp_decay). The component is scaled by the Q-value which
// is the average number of times a molecule is seen in addition to the first time.
// At large spatial frequencies the scaled component should match the numerator,
// i.e. at high resolution (low FIRE number) the numerator is made up of repeat
// localisations of the same molecule and not actual structure in the image.
// The best fit is where the numerator equals the scaled component, i.e. num / (q*exp_decay) == 1.
// The FRC Numerator is plotted and Q can be determined by
// adjusting Q and the precision mean and SD to maximise the cost function.
// This can be done interactively by the user with the effect on the FRC curve
// dynamically updated and displayed.
// Compute the scaled FRC numerator
double qNorm = (1 / frcCurve.mean1 + 1 / frcCurve.mean2);
double[] frcnum = new double[frcCurve.getSize()];
for (int i = 0; i < frcnum.length; i++) {
FRCCurveResult r = frcCurve.get(i);
frcnum[i] = qNorm * r.getNumerator() / r.getNumberOfSamples();
}
// Compute the spatial frequency and the region for curve fitting
double[] q = FRC.computeQ(frcCurve, false);
int low = 0, high = q.length;
while (high > 0 && q[high - 1] > maxQ) high--;
while (low < q.length && q[low] < minQ) low++;
// Require we fit at least 10% of the curve
if (high - low < q.length * 0.1) {
IJ.error(TITLE, "Not enough points for Q estimation");
return;
}
// Obtain initial estimate of Q plateau height and decay.
// This can be done by fitting the precision histogram and then fixing the mean and sigma.
// Or it can be done by allowing the precision to be sampled and the mean and sigma
// become parameters for fitting.
// Check if we can sample precision values
boolean sampleDecay = precision != null && FIRE.sampleDecay;
double[] exp_decay;
if (sampleDecay) {
// Random sample of precision values from the distribution is used to
// construct the decay curve
int[] sample = Random.sample(10000, precision.getN(), new Well19937c());
final double four_pi2 = 4 * Math.PI * Math.PI;
double[] pre = new double[q.length];
for (int i = 1; i < q.length; i++) pre[i] = -four_pi2 * q[i] * q[i];
// Sample
final int n = sample.length;
double[] hq = new double[n];
for (int j = 0; j < n; j++) {
// Scale to SR pixels
double s2 = precision.getValue(sample[j]) / images.nmPerPixel;
s2 *= s2;
for (int i = 1; i < q.length; i++) hq[i] += FastMath.exp(pre[i] * s2);
}
for (int i = 1; i < q.length; i++) hq[i] /= n;
exp_decay = new double[q.length];
exp_decay[0] = 1;
for (int i = 1; i < q.length; i++) {
double sinc_q = sinc(Math.PI * q[i]);
exp_decay[i] = sinc_q * sinc_q * hq[i];
}
} else {
// Note: The sigma mean and std should be in the units of super-resolution
// pixels so scale to SR pixels
exp_decay = computeExpDecay(histogram.mean / images.nmPerPixel, histogram.sigma / images.nmPerPixel, q);
}
// Smoothing
double[] smooth;
if (loessSmoothing) {
// Note: This computes the log then smooths it
double bandwidth = 0.1;
int robustness = 0;
double[] l = new double[exp_decay.length];
for (int i = 0; i < l.length; i++) {
// Original Matlab code computes the log for each array.
// This is equivalent to a single log on the fraction of the two.
// Perhaps the two log method is more numerically stable.
//l[i] = Math.log(Math.abs(frcnum[i])) - Math.log(exp_decay[i]);
l[i] = Math.log(Math.abs(frcnum[i] / exp_decay[i]));
}
try {
LoessInterpolator loess = new LoessInterpolator(bandwidth, robustness);
smooth = loess.smooth(q, l);
} catch (Exception e) {
IJ.error(TITLE, "LOESS smoothing failed");
return;
}
} else {
// Note: This smooths the curve before computing the log
double[] norm = new double[exp_decay.length];
for (int i = 0; i < norm.length; i++) {
norm[i] = frcnum[i] / exp_decay[i];
}
// Median window of 5 == radius of 2
MedianWindow mw = new MedianWindow(norm, 2);
smooth = new double[exp_decay.length];
for (int i = 0; i < norm.length; i++) {
smooth[i] = Math.log(Math.abs(mw.getMedian()));
mw.increment();
}
}
// Fit with quadratic to find the initial guess.
// Note: example Matlab code frc_Qcorrection7.m identifies regions of the
// smoothed log curve with low derivative and only fits those. The fit is
// used for the final estimate. Fitting a subset with low derivative is not
// implemented here since the initial estimate is subsequently optimised
// to maximise a cost function.
Quadratic curve = new Quadratic();
SimpleCurveFitter fit = SimpleCurveFitter.create(curve, new double[2]);
WeightedObservedPoints points = new WeightedObservedPoints();
for (int i = low; i < high; i++) points.add(q[i], smooth[i]);
double[] estimate = fit.fit(points.toList());
double qValue = FastMath.exp(estimate[0]);
//System.out.printf("Initial q-estimate = %s => %.3f\n", Arrays.toString(estimate), qValue);
// This could be made an option. Just use for debugging
boolean debug = false;
if (debug) {
// Plot the initial fit and the fit curve
double[] qScaled = FRC.computeQ(frcCurve, true);
double[] line = new double[q.length];
for (int i = 0; i < q.length; i++) line[i] = curve.value(q[i], estimate);
String title = TITLE + " Initial fit";
Plot2 plot = new Plot2(title, "Spatial Frequency (nm^-1)", "FRC Numerator");
String label = String.format("Q = %.3f", qValue);
plot.addPoints(qScaled, smooth, Plot.LINE);
plot.setColor(Color.red);
plot.addPoints(qScaled, line, Plot.LINE);
plot.setColor(Color.black);
plot.addLabel(0, 0, label);
Utils.display(title, plot, Utils.NO_TO_FRONT);
}
if (fitPrecision) {
// Q - Should this be optional?
if (sampleDecay) {
// If a sample of the precision was used to construct the data for the initial fit
// then update the estimate using the fit result since it will be a better start point.
histogram.sigma = precision.getStandardDeviation();
// Normalise sum-of-squares to the SR pixel size
double meanSumOfSquares = (precision.getSumOfSquares() / (images.nmPerPixel * images.nmPerPixel)) / precision.getN();
histogram.mean = images.nmPerPixel * Math.sqrt(meanSumOfSquares - estimate[1] / (4 * Math.PI * Math.PI));
}
// Do a multivariate fit ...
SimplexOptimizer opt = new SimplexOptimizer(1e-6, 1e-10);
PointValuePair p = null;
MultiPlateauness f = new MultiPlateauness(frcnum, q, low, high);
double[] initial = new double[] { histogram.mean / images.nmPerPixel, histogram.sigma / images.nmPerPixel, qValue };
p = findMin(p, opt, f, scale(initial, 0.1));
p = findMin(p, opt, f, scale(initial, 0.5));
p = findMin(p, opt, f, initial);
p = findMin(p, opt, f, scale(initial, 2));
p = findMin(p, opt, f, scale(initial, 10));
if (p != null) {
double[] point = p.getPointRef();
histogram.mean = point[0] * images.nmPerPixel;
histogram.sigma = point[1] * images.nmPerPixel;
qValue = point[2];
}
} else {
// If so then this should be optional.
if (sampleDecay) {
if (precisionMethod != PrecisionMethod.FIXED) {
histogram.sigma = precision.getStandardDeviation();
// Normalise sum-of-squares to the SR pixel size
double meanSumOfSquares = (precision.getSumOfSquares() / (images.nmPerPixel * images.nmPerPixel)) / precision.getN();
histogram.mean = images.nmPerPixel * Math.sqrt(meanSumOfSquares - estimate[1] / (4 * Math.PI * Math.PI));
}
exp_decay = computeExpDecay(histogram.mean / images.nmPerPixel, histogram.sigma / images.nmPerPixel, q);
}
// Estimate spurious component by promoting plateauness.
// The Matlab code used random initial points for a Simplex optimiser.
// A Brent line search should be pretty deterministic so do simple repeats.
// However it will proceed downhill so if the initial point is wrong then
// it will find a sub-optimal result.
UnivariateOptimizer o = new BrentOptimizer(1e-3, 1e-6);
Plateauness f = new Plateauness(frcnum, exp_decay, low, high);
UnivariatePointValuePair p = null;
p = findMin(p, o, f, qValue, 0.1);
p = findMin(p, o, f, qValue, 0.2);
p = findMin(p, o, f, qValue, 0.333);
p = findMin(p, o, f, qValue, 0.5);
// Do some Simplex repeats as well
SimplexOptimizer opt = new SimplexOptimizer(1e-6, 1e-10);
p = findMin(p, opt, f, qValue * 0.1);
p = findMin(p, opt, f, qValue * 0.5);
p = findMin(p, opt, f, qValue);
p = findMin(p, opt, f, qValue * 2);
p = findMin(p, opt, f, qValue * 10);
if (p != null)
qValue = p.getPoint();
}
QPlot qplot = new QPlot(frcCurve, qValue, low, high);
// Interactive dialog to estimate Q (blinking events per flourophore) using
// sliders for the mean and standard deviation of the localisation precision.
showQEstimationDialog(histogram, qplot, frcCurve, images.nmPerPixel);
IJ.showStatus(TITLE + " complete");
}
use of org.apache.commons.math3.random.Well19937c in project GDSC-SMLM by aherbert.
the class JumpDistanceAnalysis method createCMAESOptimizer.
private CMAESOptimizer createCMAESOptimizer() {
double rel = 1e-8;
double abs = 1e-10;
int maxIterations = 2000;
//Double.NEGATIVE_INFINITY;
double stopFitness = 0;
boolean isActiveCMA = true;
int diagonalOnly = 20;
int checkFeasableCount = 1;
RandomGenerator random = new Well19937c();
boolean generateStatistics = false;
ConvergenceChecker<PointValuePair> checker = new SimpleValueChecker(rel, abs);
// Iterate this for stability in the initial guess
return new CMAESOptimizer(maxIterations, stopFitness, isActiveCMA, diagonalOnly, checkFeasableCount, random, generateStatistics, checker);
}
Aggregations