List of usage examples for org.apache.commons.math3.optim PointVectorValuePair getPointRef
public double[] getPointRef()
From source file:gdsc.smlm.fitting.BinomialFitter.java
/** * Fit the binomial distribution (n,p) to the cumulative histogram. Performs fitting assuming a fixed n value and * attempts to optimise p.//from w w w . j a v a 2s . c o m * * @param histogram * The input histogram * @param mean * The histogram mean (used to estimate p). Calculated if NaN. * @param n * The n to evaluate * @param zeroTruncated * True if the model should ignore n=0 (zero-truncated binomial) * @return The best fit (n, p) * @throws IllegalArgumentException * If any of the input data values are negative * @throws IllegalArgumentException * If any fitting a zero truncated binomial and there are no values above zero */ public PointValuePair fitBinomial(double[] histogram, double mean, int n, boolean zeroTruncated) { if (Double.isNaN(mean)) mean = getMean(histogram); if (zeroTruncated && histogram[0] > 0) { log("Fitting zero-truncated histogram but there are zero values - Renormalising to ignore zero"); double cumul = 0; for (int i = 1; i < histogram.length; i++) cumul += histogram[i]; if (cumul == 0) throw new IllegalArgumentException( "Fitting zero-truncated histogram but there are no non-zero values"); histogram[0] = 0; for (int i = 1; i < histogram.length; i++) histogram[i] /= cumul; } int nFittedPoints = Math.min(histogram.length, n + 1) - ((zeroTruncated) ? 1 : 0); if (nFittedPoints < 1) { log("No points to fit (%d): Histogram.length = %d, n = %d, zero-truncated = %b", nFittedPoints, histogram.length, n, zeroTruncated); return null; } // The model is only fitting the probability p // For a binomial n*p = mean => p = mean/n double[] initialSolution = new double[] { FastMath.min(mean / n, 1) }; // Create the function BinomialModelFunction function = new BinomialModelFunction(histogram, n, zeroTruncated); double[] lB = new double[1]; double[] uB = new double[] { 1 }; SimpleBounds bounds = new SimpleBounds(lB, uB); // Fit // CMAESOptimizer or BOBYQAOptimizer support bounds // CMAESOptimiser based on Matlab code: // https://www.lri.fr/~hansen/cmaes.m // Take the defaults from the Matlab documentation int maxIterations = 2000; double stopFitness = 0; //Double.NEGATIVE_INFINITY; boolean isActiveCMA = true; int diagonalOnly = 0; int checkFeasableCount = 1; RandomGenerator random = new Well19937c(); boolean generateStatistics = false; ConvergenceChecker<PointValuePair> checker = new SimpleValueChecker(1e-6, 1e-10); // The sigma determines the search range for the variables. It should be 1/3 of the initial search region. OptimizationData sigma = new CMAESOptimizer.Sigma(new double[] { (uB[0] - lB[0]) / 3 }); OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(2)))); try { PointValuePair solution = null; boolean noRefit = maximumLikelihood; if (n == 1 && zeroTruncated) { // No need to fit solution = new PointValuePair(new double[] { 1 }, 0); noRefit = true; } else { GoalType goalType = (maximumLikelihood) ? GoalType.MAXIMIZE : GoalType.MINIMIZE; // Iteratively fit CMAESOptimizer opt = new CMAESOptimizer(maxIterations, stopFitness, isActiveCMA, diagonalOnly, checkFeasableCount, random, generateStatistics, checker); for (int iteration = 0; iteration <= fitRestarts; iteration++) { try { // Start from the initial solution PointValuePair result = opt.optimize(new InitialGuess(initialSolution), new ObjectiveFunction(function), goalType, bounds, sigma, popSize, new MaxIter(maxIterations), new MaxEval(maxIterations * 2)); //System.out.printf("CMAES Iter %d initial = %g (%d)\n", iteration, result.getValue(), // opt.getEvaluations()); if (solution == null || result.getValue() < solution.getValue()) { solution = result; } } catch (TooManyEvaluationsException e) { } catch (TooManyIterationsException e) { } if (solution == null) continue; try { // Also restart from the current optimum PointValuePair result = opt.optimize(new InitialGuess(solution.getPointRef()), new ObjectiveFunction(function), goalType, bounds, sigma, popSize, new MaxIter(maxIterations), new MaxEval(maxIterations * 2)); //System.out.printf("CMAES Iter %d restart = %g (%d)\n", iteration, result.getValue(), // opt.getEvaluations()); if (result.getValue() < solution.getValue()) { solution = result; } } catch (TooManyEvaluationsException e) { } catch (TooManyIterationsException e) { } } if (solution == null) return null; } if (noRefit) { // Although we fit the log-likelihood, return the sum-of-squares to allow // comparison across different n double p = solution.getPointRef()[0]; double ss = 0; double[] obs = function.p; double[] exp = function.getP(p); for (int i = 0; i < obs.length; i++) ss += (obs[i] - exp[i]) * (obs[i] - exp[i]); return new PointValuePair(solution.getPointRef(), ss); } // We can do a LVM refit if the number of fitted points is more than 1 else if (nFittedPoints > 1) { // Improve SS fit with a gradient based LVM optimizer LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer(); try { final BinomialModelFunctionGradient gradientFunction = new BinomialModelFunctionGradient( histogram, n, zeroTruncated); PointVectorValuePair lvmSolution = optimizer.optimize(new MaxIter(3000), new MaxEval(Integer.MAX_VALUE), new ModelFunctionJacobian(new MultivariateMatrixFunction() { public double[][] value(double[] point) throws IllegalArgumentException { return gradientFunction.jacobian(point); } }), new ModelFunction(gradientFunction), new Target(gradientFunction.p), new Weight(gradientFunction.getWeights()), new InitialGuess(solution.getPointRef())); double ss = 0; double[] obs = gradientFunction.p; double[] exp = lvmSolution.getValue(); for (int i = 0; i < obs.length; i++) ss += (obs[i] - exp[i]) * (obs[i] - exp[i]); // Check the pValue is valid since the LVM is not bounded. double p = lvmSolution.getPointRef()[0]; if (ss < solution.getValue() && p <= 1 && p >= 0) { //log("Re-fitting improved the SS from %s to %s (-%s%%)", // Utils.rounded(solution.getValue(), 4), Utils.rounded(ss, 4), // Utils.rounded(100 * (solution.getValue() - ss) / solution.getValue(), 4)); return new PointValuePair(lvmSolution.getPoint(), ss); } } catch (TooManyIterationsException e) { log("Failed to re-fit: Too many iterations (%d)", optimizer.getIterations()); } catch (ConvergenceException e) { log("Failed to re-fit: %s", e.getMessage()); } catch (Exception e) { // Ignore this ... } } return solution; } catch (Exception e) { log("Failed to fit Binomial distribution with N=%d : %s", n, e.getMessage()); } return null; }
From source file:gdsc.utils.Cell_Outliner.java
/** * Find an ellipse that optimises the fit to the polygon detected edges. * /*from w w w . j av a 2s . c o m*/ * @param roi * @param params * @param weightMap * @param angle * @return */ private double[] fitPolygonalCell(PolygonRoi roi, double[] params, FloatProcessor weightMap, FloatProcessor angle) { // Get an estimate of the starting parameters using the current polygon double[] startPoint = params; startPoint = estimateStartPoint(roi, weightMap.getWidth(), weightMap.getHeight()); int maxEval = 2000; final DifferentiableEllipticalFitFunction func = new DifferentiableEllipticalFitFunction(roi, weightMap); double relativeThreshold = 100 * Precision.EPSILON; double absoluteThreshold = 100 * Precision.SAFE_MIN; ConvergenceChecker<PointVectorValuePair> checker = new SimplePointChecker<PointVectorValuePair>( relativeThreshold, absoluteThreshold); double initialStepBoundFactor = 10; double costRelativeTolerance = 1e-10; double parRelativeTolerance = 1e-10; double orthoTolerance = 1e-10; double threshold = Precision.SAFE_MIN; LevenbergMarquardtOptimizer optimiser = new LevenbergMarquardtOptimizer(initialStepBoundFactor, checker, costRelativeTolerance, parRelativeTolerance, orthoTolerance, threshold); try { PointVectorValuePair solution = optimiser.optimize(new MaxIter(maxEval), new MaxEval(Integer.MAX_VALUE), new ModelFunctionJacobian(new MultivariateMatrixFunction() { public double[][] value(double[] point) throws IllegalArgumentException { return func.jacobian(point); } }), new ModelFunction(func), new Target(func.calculateTarget()), new Weight(func.calculateWeights()), new InitialGuess(startPoint)); if (debug) IJ.log(String.format("Eval = %d (Iter = %d), RMS = %f", optimiser.getEvaluations(), optimiser.getIterations(), optimiser.getRMS())); return solution.getPointRef(); } catch (Exception e) { IJ.log("Failed to find an elliptical solution, defaulting to polygon"); e.printStackTrace(); } return null; }
From source file:gdsc.smlm.ij.plugins.TraceDiffusion.java
/** * Fit the jump distance histogram using a cumulative sum as detailed in * <p>//from ww w . j ava 2 s .co m * Update the plot by adding the fit line. * * @param jumpDistances * @param jdHistogram * @param title * @param plot * @return */ private double[][] fitJumpDistance(StoredDataStatistics jumpDistances, double[][] jdHistogram, String title, Plot2 plot) { final double meanDistance = Math.sqrt(jumpDistances.getMean()) * 1e3; final double beta = meanDistance / precision; Utils.log( "Jump Distance analysis : N = %d, Time = %d frames (%s seconds). Mean Distance = %s nm, Precision = %s nm, Beta = %s", jumpDistances.getN(), settings.jumpDistance, Utils.rounded(settings.jumpDistance * exposureTime, 4), Utils.rounded(meanDistance, 4), Utils.rounded(precision, 4), Utils.rounded(beta, 4)); int n = 0; int N = 10; double[] SS = new double[N]; Arrays.fill(SS, -1); double[] ic = new double[N]; Arrays.fill(ic, Double.POSITIVE_INFINITY); double[][] coefficients = new double[N][]; double[][] fractions = new double[N][]; double[][] fitParams = new double[N][]; double bestIC = Double.POSITIVE_INFINITY; int best = -1; // Guess the D final double estimatedD = jumpDistances.getMean() / 4; Utils.log("Estimated D = %s um^2/s", Utils.rounded(estimatedD, 4)); // Fit using a single population model LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer(); try { final JumpDistanceFunction function = new JumpDistanceFunction(jdHistogram[0], jdHistogram[1], estimatedD); PointVectorValuePair lvmSolution = optimizer.optimize(new MaxIter(3000), new MaxEval(Integer.MAX_VALUE), new ModelFunctionJacobian(new MultivariateMatrixFunction() { public double[][] value(double[] point) throws IllegalArgumentException { return function.jacobian(point); } }), new ModelFunction(function), new Target(function.getY()), new Weight(function.getWeights()), new InitialGuess(function.guess())); fitParams[n] = lvmSolution.getPointRef(); SS[n] = calculateSumOfSquares(function.getY(), lvmSolution.getValueRef()); ic[n] = Maths.getInformationCriterion(SS[n], function.x.length, 1); coefficients[n] = fitParams[n]; fractions[n] = new double[] { 1 }; Utils.log("Fit Jump distance (N=%d) : D = %s um^2/s, SS = %f, IC = %s (%d evaluations)", n + 1, Utils.rounded(fitParams[n][0], 4), SS[n], Utils.rounded(ic[n], 4), optimizer.getEvaluations()); bestIC = ic[n]; best = 0; addToPlot(function, fitParams[n], jdHistogram, title, plot, Color.magenta); } catch (TooManyIterationsException e) { Utils.log("Failed to fit : Too many iterations (%d)", optimizer.getIterations()); } catch (ConvergenceException e) { Utils.log("Failed to fit : %s", e.getMessage()); } n++; // Fit using a mixed population model. // Vary n from 2 to N. Stop when the fit fails or the fit is worse. int bestMulti = -1; double bestMultiIC = Double.POSITIVE_INFINITY; while (n < N) { // Uses a weighted sum of n exponential functions, each function models a fraction of the particles. // An LVM fit cannot restrict the parameters so the fractions do not go below zero. // Use the CMEASOptimizer which supports bounded fitting. MixedJumpDistanceFunctionMultivariate mixedFunction = new MixedJumpDistanceFunctionMultivariate( jdHistogram[0], jdHistogram[1], estimatedD, n + 1); double[] lB = mixedFunction.getLowerBounds(); double[] uB = mixedFunction.getUpperBounds(); SimpleBounds bounds = new SimpleBounds(lB, uB); int maxIterations = 2000; double stopFitness = 0; //Double.NEGATIVE_INFINITY; boolean isActiveCMA = true; int diagonalOnly = 20; int checkFeasableCount = 1; RandomGenerator random = new Well19937c(); boolean generateStatistics = false; ConvergenceChecker<PointValuePair> checker = new SimpleValueChecker(1e-6, 1e-10); // The sigma determines the search range for the variables. It should be 1/3 of the initial search region. double[] s = new double[lB.length]; for (int i = 0; i < s.length; i++) s[i] = (uB[i] - lB[i]) / 3; OptimizationData sigma = new CMAESOptimizer.Sigma(s); OptimizationData popSize = new CMAESOptimizer.PopulationSize( (int) (4 + Math.floor(3 * Math.log(mixedFunction.x.length)))); // Iterate this for stability in the initial guess CMAESOptimizer opt = new CMAESOptimizer(maxIterations, stopFitness, isActiveCMA, diagonalOnly, checkFeasableCount, random, generateStatistics, checker); int evaluations = 0; PointValuePair constrainedSolution = null; for (int i = 0; i <= settings.fitRestarts; i++) { // Try from the initial guess try { PointValuePair solution = opt.optimize(new InitialGuess(mixedFunction.guess()), new ObjectiveFunction(mixedFunction), GoalType.MINIMIZE, bounds, sigma, popSize, new MaxIter(maxIterations), new MaxEval(maxIterations * 2)); if (constrainedSolution == null || solution.getValue() < constrainedSolution.getValue()) { evaluations = opt.getEvaluations(); constrainedSolution = solution; //Utils.log("[%da] Fit Jump distance (N=%d) : SS = %f (%d evaluations)", i, n + 1, // solution.getValue(), evaluations); } } catch (TooManyEvaluationsException e) { } if (constrainedSolution == null) continue; // Try from the current optimum try { PointValuePair solution = opt.optimize(new InitialGuess(constrainedSolution.getPointRef()), new ObjectiveFunction(mixedFunction), GoalType.MINIMIZE, bounds, sigma, popSize, new MaxIter(maxIterations), new MaxEval(maxIterations * 2)); if (constrainedSolution == null || solution.getValue() < constrainedSolution.getValue()) { evaluations = opt.getEvaluations(); constrainedSolution = solution; //Utils.log("[%db] Fit Jump distance (N=%d) : SS = %f (%d evaluations)", i, n + 1, // solution.getValue(), evaluations); } } catch (TooManyEvaluationsException e) { } } if (constrainedSolution == null) { Utils.log("Failed to fit N=%d", n + 1); break; } fitParams[n] = constrainedSolution.getPointRef(); SS[n] = constrainedSolution.getValue(); // TODO - Try a bounded BFGS optimiser // Try and improve using a LVM fit final MixedJumpDistanceFunctionGradient mixedFunctionGradient = new MixedJumpDistanceFunctionGradient( jdHistogram[0], jdHistogram[1], estimatedD, n + 1); PointVectorValuePair lvmSolution; try { lvmSolution = optimizer.optimize(new MaxIter(3000), new MaxEval(Integer.MAX_VALUE), new ModelFunctionJacobian(new MultivariateMatrixFunction() { public double[][] value(double[] point) throws IllegalArgumentException { return mixedFunctionGradient.jacobian(point); } }), new ModelFunction(mixedFunctionGradient), new Target(mixedFunctionGradient.getY()), new Weight(mixedFunctionGradient.getWeights()), new InitialGuess(fitParams[n])); double ss = calculateSumOfSquares(mixedFunctionGradient.getY(), lvmSolution.getValue()); // All fitted parameters must be above zero if (ss < SS[n] && Maths.min(lvmSolution.getPoint()) > 0) { //Utils.log(" Re-fitting improved the SS from %s to %s (-%s%%)", Utils.rounded(SS[n], 4), // Utils.rounded(ss, 4), Utils.rounded(100 * (SS[n] - ss) / SS[n], 4)); fitParams[n] = lvmSolution.getPoint(); SS[n] = ss; evaluations += optimizer.getEvaluations(); } } catch (TooManyIterationsException e) { //Utils.log("Failed to re-fit : Too many evaluations (%d)", optimizer.getEvaluations()); } catch (ConvergenceException e) { //Utils.log("Failed to re-fit : %s", e.getMessage()); } // Since the fractions must sum to one we subtract 1 degree of freedom from the number of parameters ic[n] = Maths.getInformationCriterion(SS[n], mixedFunction.x.length, fitParams[n].length - 1); double[] d = new double[n + 1]; double[] f = new double[n + 1]; double sum = 0; for (int i = 0; i < d.length; i++) { f[i] = fitParams[n][i * 2]; sum += f[i]; d[i] = fitParams[n][i * 2 + 1]; } for (int i = 0; i < f.length; i++) f[i] /= sum; // Sort by coefficient size sort(d, f); coefficients[n] = d; fractions[n] = f; Utils.log("Fit Jump distance (N=%d) : D = %s um^2/s (%s), SS = %f, IC = %s (%d evaluations)", n + 1, format(d), format(f), SS[n], Utils.rounded(ic[n], 4), evaluations); boolean valid = true; for (int i = 0; i < f.length; i++) { // Check the fit has fractions above the minimum fraction if (f[i] < minFraction) { Utils.log("Fraction is less than the minimum fraction: %s < %s", Utils.rounded(f[i]), Utils.rounded(minFraction)); valid = false; break; } // Check the coefficients are different if (i + 1 < f.length && d[i] / d[i + 1] < minDifference) { Utils.log("Coefficients are not different: %s / %s = %s < %s", Utils.rounded(d[i]), Utils.rounded(d[i + 1]), Utils.rounded(d[i] / d[i + 1]), Utils.rounded(minDifference)); valid = false; break; } } if (!valid) break; // Store the best model if (bestIC > ic[n]) { bestIC = ic[n]; best = n; } // Store the best multi model if (bestMultiIC < ic[n]) { break; } bestMultiIC = ic[n]; bestMulti = n; n++; } // Add the best fit to the plot and return the parameters. if (bestMulti > -1) { Function function = new MixedJumpDistanceFunctionMultivariate(jdHistogram[0], jdHistogram[1], 0, bestMulti + 1); addToPlot(function, fitParams[bestMulti], jdHistogram, title, plot, Color.yellow); } if (best > -1) { Utils.log("Best fit achieved using %d population%s: D = %s um^2/s, Fractions = %s", best + 1, (best == 0) ? "" : "s", format(coefficients[best]), format(fractions[best])); } return (best > -1) ? new double[][] { coefficients[best], fractions[best] } : null; }
From source file:org.eclipse.dawnsci.analysis.dataset.roi.fitting.CircleCoordinatesFunction.java
@Override public void geometricFit(IDataset x, IDataset y, double[] init) { residual = Double.NaN;// w w w.ja va 2 s. c om if (x.getSize() < PARAMETERS || y.getSize() < PARAMETERS) { throw new IllegalArgumentException("Need " + PARAMETERS + " or more points"); } if (x.getSize() == PARAMETERS) { init = quickfit(x, y); for (int i = 0; i < PARAMETERS; i++) parameters[i] = init[i]; residual = 0; return; } if (init == null) init = quickfit(x, y); else if (init.length < PARAMETERS) throw new IllegalArgumentException("Need " + PARAMETERS + " parameters"); CircleCoordinatesFunction f = (CircleCoordinatesFunction) getFitFunction(x, y); LevenbergMarquardtOptimizer opt = new LevenbergMarquardtOptimizer(); try { PointVectorValuePair result = opt.optimize(new ModelFunction(f), new ModelFunctionJacobian(f.jacobian()), f.getTarget(), f.getWeight(), f.calcAllInitValues(init), new MaxEval(MAX_EVALUATIONS)); double[] point = result.getPointRef(); for (int i = 0; i < PARAMETERS; i++) parameters[i] = point[i]; residual = opt.getRMS(); logger.trace("Circle fit: rms = {}, x^2 = {}", opt.getRMS(), opt.getChiSquare()); } catch (DimensionMismatchException e) { // cannot happen } catch (IllegalArgumentException e) { // should not happen! } catch (TooManyEvaluationsException e) { throw new IllegalArgumentException("Problem with optimizer converging"); } }
From source file:org.eclipse.dawnsci.analysis.dataset.roi.fitting.EllipseCoordinatesFunction.java
@Override public void geometricFit(IDataset x, IDataset y, double[] init) { residual = Double.NaN;/* w w w . java2 s .c o m*/ if (x.getSize() < PARAMETERS || y.getSize() < PARAMETERS) { throw new IllegalArgumentException("Need " + PARAMETERS + " or more points"); } if (init == null) init = quickfit(x, y); else if (init.length < PARAMETERS) throw new IllegalArgumentException("Need " + PARAMETERS + " parameters"); EllipseCoordinatesFunction f = (EllipseCoordinatesFunction) getFitFunction(x, y); LevenbergMarquardtOptimizer opt = new LevenbergMarquardtOptimizer(); try { PointVectorValuePair result = opt.optimize(new ModelFunction(f), new ModelFunctionJacobian(f.jacobian()), f.getTarget(), f.getWeight(), f.calcAllInitValues(init), new MaxEval(MAX_EVALUATIONS)); double[] point = result.getPointRef(); for (int i = 0; i < PARAMETERS; i++) parameters[i] = point[i]; residual = opt.getRMS(); logger.trace("Ellipse fit: rms = {}, x^2 = {}", residual, opt.getChiSquare()); } catch (DimensionMismatchException e) { // cannot happen } catch (IllegalArgumentException e) { // should not happen! } catch (TooManyEvaluationsException e) { throw new IllegalArgumentException("Problem with optimizer converging"); } }