List of usage examples for org.apache.commons.math3.optim.nonlinear.scalar MultivariateFunctionMappingAdapter unboundedToBounded
public double[] unboundedToBounded(double[] point)
From source file:gdsc.smlm.ij.plugins.EMGainAnalysis.java
/** * Fit the EM-gain distribution (Gaussian * Gamma) * /*from w ww . j ava 2s .c o m*/ * @param h * The distribution */ private void fit(int[] h) { final int[] limits = limits(h); final double[] x = getX(limits); final double[] y = getY(h, limits); Plot2 plot = new Plot2(TITLE, "ADU", "Frequency"); double yMax = Maths.max(y); plot.setLimits(limits[0], limits[1], 0, yMax); plot.setColor(Color.black); plot.addPoints(x, y, Plot2.DOT); Utils.display(TITLE, plot); // Estimate remaining parameters. // Assuming a gamma_distribution(shape,scale) then mean = shape * scale // scale = gain // shape = Photons = mean / gain double mean = getMean(h) - bias; // Note: if the bias is too high then the mean will be negative. Just move the bias. while (mean < 0) { bias -= 1; mean += 1; } double photons = mean / gain; if (simulate) Utils.log("Simulated bias=%d, gain=%s, noise=%s, photons=%s", (int) _bias, Utils.rounded(_gain), Utils.rounded(_noise), Utils.rounded(_photons)); Utils.log("Estimate bias=%d, gain=%s, noise=%s, photons=%s", (int) bias, Utils.rounded(gain), Utils.rounded(noise), Utils.rounded(photons)); final int max = (int) x[x.length - 1]; double[] g = pdf(max, photons, gain, noise, (int) bias); plot.setColor(Color.blue); plot.addPoints(x, g, Plot2.LINE); Utils.display(TITLE, plot); // Perform a fit CustomPowellOptimizer o = new CustomPowellOptimizer(1e-6, 1e-16, 1e-6, 1e-16); double[] startPoint = new double[] { photons, gain, noise, bias }; final int maxEval = 3000; // Set bounds double[] lower = new double[] { 0, 0.5 * gain, 0, bias - noise }; double[] upper = new double[] { (limits[1] - limits[0]) / gain, 2 * gain, gain, bias + noise }; // Restart until converged. // TODO - Maybe fix this with a better optimiser. This needs to be tested on real data. PointValuePair solution = null; for (int iter = 0; iter < 3; iter++) { IJ.showStatus("Fitting histogram ... Iteration " + iter); try { // Basic Powell optimiser MultivariateFunction fun = getFunction(limits, y, max, maxEval); PointValuePair optimum = o.optimize(new MaxEval(maxEval), new ObjectiveFunction(fun), GoalType.MINIMIZE, new InitialGuess((solution == null) ? startPoint : solution.getPointRef())); if (solution == null || optimum.getValue() < solution.getValue()) { solution = optimum; } } catch (Exception e) { } try { // Bounded Powell optimiser MultivariateFunction fun = getFunction(limits, y, max, maxEval); MultivariateFunctionMappingAdapter adapter = new MultivariateFunctionMappingAdapter(fun, lower, upper); PointValuePair optimum = o.optimize(new MaxEval(maxEval), new ObjectiveFunction(adapter), GoalType.MINIMIZE, new InitialGuess(adapter .boundedToUnbounded((solution == null) ? startPoint : solution.getPointRef()))); double[] point = adapter.unboundedToBounded(optimum.getPointRef()); optimum = new PointValuePair(point, optimum.getValue()); if (solution == null || optimum.getValue() < solution.getValue()) { solution = optimum; } } catch (Exception e) { } } IJ.showStatus(""); IJ.showProgress(1); if (solution == null) { Utils.log("Failed to fit the distribution"); return; } double[] point = solution.getPointRef(); photons = point[0]; gain = point[1]; noise = point[2]; bias = (int) Math.round(point[3]); String label = String.format("Fitted bias=%d, gain=%s, noise=%s, photons=%s", (int) bias, Utils.rounded(gain), Utils.rounded(noise), Utils.rounded(photons)); Utils.log(label); if (simulate) { Utils.log("Relative Error bias=%s, gain=%s, noise=%s, photons=%s", Utils.rounded(relativeError(bias, _bias)), Utils.rounded(relativeError(gain, _gain)), Utils.rounded(relativeError(noise, _noise)), Utils.rounded(relativeError(photons, _photons))); } // Show the PoissonGammaGaussian approximation double[] f = null; if (showApproximation) { f = new double[x.length]; PoissonGammaGaussianFunction fun = new PoissonGammaGaussianFunction(1.0 / gain, noise); final double expected = photons * gain; for (int i = 0; i < f.length; i++) { f[i] = fun.likelihood(x[i] - bias, expected); //System.out.printf("x=%d, g=%f, f=%f, error=%f\n", (int) x[i], g[i], f[i], // gdsc.smlm.fitting.utils.DoubleEquality.relativeError(g[i], f[i])); } yMax = Maths.maxDefault(yMax, f); } // Replot g = pdf(max, photons, gain, noise, (int) bias); plot = new Plot2(TITLE, "ADU", "Frequency"); plot.setLimits(limits[0], limits[1], 0, yMax * 1.05); plot.setColor(Color.black); plot.addPoints(x, y, Plot2.DOT); plot.setColor(Color.red); plot.addPoints(x, g, Plot2.LINE); plot.addLabel(0, 0, label); if (showApproximation) { plot.setColor(Color.blue); plot.addPoints(x, f, Plot2.LINE); } Utils.display(TITLE, plot); }
From source file:org.hawkular.datamining.forecast.models.AbstractModelOptimizer.java
protected void optimize(double[] initialGuess, MultivariateFunctionMappingAdapter costFunction) { // Nelder-Mead Simplex SimplexOptimizer optimizer = new SimplexOptimizer(0.0001, 0.0001); PointValuePair unBoundedResult = optimizer.optimize(GoalType.MINIMIZE, new MaxIter(MAX_ITER), new MaxEval(MAX_EVAL), new InitialGuess(initialGuess), new ObjectiveFunction(costFunction), new NelderMeadSimplex(initialGuess.length)); result = costFunction.unboundedToBounded(unBoundedResult.getPoint()); }
From source file:org.hawkular.datamining.forecast.models.performance.OptimizationAlgorithmsTests.java
private void printOptimizationResult(MultivariateFunctionMappingAdapter adapter, double[] unbounded, ModelData modelData) {/*from www. j av a2s. com*/ double[] boundedResult = adapter.unboundedToBounded(unbounded); DoubleExponentialSmoothing bestModel = DoubleExponentialSmoothing.createCustom(boundedResult[0], boundedResult[1], ImmutableMetricContext.getDefault()); bestModel.init(modelData.getData()); System.out.println(bestModel.initStatistics()); }
From source file:smlm.fitting.FittingClassical.java
@SuppressWarnings("unused") @Override/*from w w w . j a v a2 s .c o m*/ public boolean FitThis() { // Initial estimates (your initial x) // Calculating the centroid ImageStatistics stat = roi.getStatistics(); double x0 = stat.xCenterOfMass; double y0 = stat.yCenterOfMass; double[] start = { x0, y0, param.psfSigma, param.psfSigma, i0Max, backgroundLevel }; PointValuePair solutionMult = null; MultivariateFunctionMappingAdapter fitFunc = null; if (param.fitting != Params.Fitting.CentroidFit) { // initial step sizes (take a good guess) double[] step = { 1, 1, 0.1, 0.1, stdBackground, stdBackground }; // convergence tolerance double ftol = 0.0001;// 0.000001; pixelPrecision = 3;// 5; if (param.fitting == Params.Fitting.FastGaussianFit) { ftol = 0.1; pixelPrecision = 3; } SimplexOptimizer Fit = new SimplexOptimizer(ftol, ftol * ftol); double[] low = new double[] { 0, 0, 0, 0, 0, 0 }; double[] up = new double[] { roi.getWidth(), roi.getHeight(), Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY }; fitFunc = new MultivariateFunctionMappingAdapter(new LLH(), low, up); // maximal number of iterations int maxIter = 5000; // Nelder and Mead maximisation procedure // x0 e [0, xmax] // Fit.addConstraint(0, -1, 0); // Fit.addConstraint(0, 1, roi.getWidth()); // y0 e [0, ymax] // Fit.addConstraint(1, -1, 0); // Fit.addConstraint(1, 1, roi.getHeight()); /* * // sigmax e [PSFSigma/3, 3*PSFSigma] Fit.addConstraint(2, -1, * PSFSigmaInt-PSFSigmaInt/2); Fit.addConstraint(2, 1, * 2*PSFSigmaInt); // sigmay e [PSFSigma/3, 3*PSFSigma] * Fit.addConstraint(3, -1, PSFSigmaInt-PSFSigmaInt/2); * Fit.addConstraint(3, 1, 2*PSFSigmaInt); // I0 e [StdBackground, * 50*Intensities[i]] Fit.addConstraint(4, -1, StdBackground); * Fit.addConstraint(4, 1, 50*XY[3][i]); // PoissonNoise e * [BackgroundLevel/3, 3*BackgroundLevel] Fit.addConstraint(5, -1, * BackgroundLevel/3); Fit.addConstraint(5, 1, 3*BackgroundLevel); */ solutionMult = Fit.optimize(new MaxEval(maxIter), new ObjectiveFunction(fitFunc), GoalType.MAXIMIZE, new InitialGuess(fitFunc.boundedToUnbounded(start)), new MultiDirectionalSimplex(step)); } // Result of minimisation // Save the fit results this.fit.incrementCounter(); if (param.fitting == Params.Fitting.CentroidFit) { results = start; } else { results = fitFunc.unboundedToBounded(solutionMult.getPoint()); } for (int i = 0; i < isArgumentFixed.length; i++) { if (isArgumentFixed[i]) { results[i] = fixedArguments[i]; } } if (cycle != -1) this.fit.addValue(ResultsTableMt.CYCLE, cycle); this.fit.addValue(ResultsTableMt.FRAME, realFrame); this.fit.addValue(ResultsTableMt.X0, results[0] + xMax - ((double) roiWidth * param.psfSigmaInt)); this.fit.addValue(ResultsTableMt.Y0, results[1] + yMax - ((double) roiWidth * param.psfSigmaInt)); this.fit.addValue(ResultsTableMt.SIGMAX, results[2]); this.fit.addValue(ResultsTableMt.SIGMAY, results[3]); this.fit.addValue(ResultsTableMt.I0, results[4]); this.fit.addValue(ResultsTableMt.NOISE, results[5]); if (param.fitting == Params.Fitting.CentroidFit || true) { this.fit.addValue(ResultsTableMt.IS_FITTED, 1); if (param.fitting != Params.Fitting.CentroidFit) { this.fit.addValue("MinFit", solutionMult.getValue()); } } else this.fit.addValue(ResultsTableMt.IS_FITTED, 0); // Save results if (param.debug) { DrawInitialRoi(true); DrawFit(true); } return true; }