List of usage examples for org.apache.commons.math3.optim.nonlinear.scalar.noderiv CMAESOptimizer CMAESOptimizer
public CMAESOptimizer(int maxIterations, double stopFitness, boolean isActiveCMA, int diagonalOnly, int checkFeasableCount, RandomGenerator random, boolean generateStatistics, ConvergenceChecker<PointValuePair> checker)
From source file:net.sf.tweety.math.opt.solver.ApacheCommonsCMAESOptimizer.java
@Override public Map<Variable, Term> solve(ConstraintSatisfactionProblem problem) throws GeneralMathException { // only optimization problems if (!(problem instanceof OptimizationProblem)) throw new IllegalArgumentException("Only optimization problems allowed for this solver."); OptimizationProblem p = (OptimizationProblem) problem; // only box constraints allowed (so far) if (!p.isEmpty()) throw new IllegalArgumentException( "Only optimization problems with box constraints on variables allowed for this solver (no other constraints."); final List<Variable> vars = new ArrayList<Variable>(p.getTargetFunction().getVariables()); double[] lb = new double[vars.size()]; double[] ub = new double[vars.size()]; double[] s = new double[vars.size()]; double[] sigma = new double[vars.size()]; for (int i = 0; i < vars.size(); i++) { lb[i] = vars.get(i).getLowerBound(); ub[i] = vars.get(i).getUpperBound(); s[i] = (lb[i] + ub[i]) / 2;/*from w w w . java2 s . c o m*/ sigma[i] = ub[i] - lb[i]; } final Term targetFunction = p.getTargetFunction(); MultivariateFunction target = new MultivariateFunction() { @Override public double value(double[] arg0) { return targetFunction.replaceAllTerms(arg0, vars).doubleValue(); } }; // construct solver CMAESOptimizer optimizer = new CMAESOptimizer(this.maxIterations, this.stopFitness, this.isActiveCMA, this.diagonalOnly, this.checkFeasableCount, new JDKRandomGenerator(), false, new SimplePointChecker<PointValuePair>(this.precision, this.precision)); PointValuePair val = optimizer.optimize(new CMAESOptimizer.Sigma(sigma), new ObjectiveFunction(target), new InitialGuess(s), p.getType() == OptimizationProblem.MAXIMIZE ? GoalType.MAXIMIZE : GoalType.MINIMIZE, new MaxEval(this.maxIterations), new SimpleBounds(lb, ub), new CMAESOptimizer.PopulationSize(this.populationSize)); Map<Variable, Term> result = new HashMap<Variable, Term>(); for (int i = 0; i < vars.size(); i++) result.put(vars.get(i), new FloatConstant(val.getPoint()[i])); return result; }
From source file:com.insightml.math.optimization.AbstractOptimizable.java
private PointValuePair cmaes(final double[] initialValues, final int initialSize) { final CMAESOptimizer cmaes = new CMAESOptimizer(30000, 0, true, 1, 0, new MersenneTwister(1334498400 * 1000), false, conv()); return optimize(cmaes, initialValues, new PopulationSize(initialSize), new Sigma(Vectors.fill(0.3, initialValues.length)), bounds); }
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 ww . jav a 2 s .co 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:edu.cmu.tetrad.sem.GeneralizedSemEstimator.java
private double[] optimize(MultivariateFunction function, double[] values, int optimizer) { PointValuePair pair = null;//from w ww.j a v a2s. c o m if (optimizer == 1) { // 0.01, 0.000001 //2.0D * FastMath.ulp(1.0D), 1e-8 MultivariateOptimizer search = new PowellOptimizer(1e-7, 1e-7); pair = search.optimize(new InitialGuess(values), new ObjectiveFunction(function), GoalType.MINIMIZE, new MaxEval(100000)); } else if (optimizer == 2) { MultivariateOptimizer search = new SimplexOptimizer(1e-7, 1e-7); pair = search.optimize(new InitialGuess(values), new ObjectiveFunction(function), GoalType.MINIMIZE, new MaxEval(100000), new NelderMeadSimplex(values.length)); } else if (optimizer == 3) { int dim = values.length; int additionalInterpolationPoints = 0; final int numIterpolationPoints = 2 * dim + 1 + additionalInterpolationPoints; BOBYQAOptimizer search = new BOBYQAOptimizer(numIterpolationPoints); pair = search.optimize(new MaxEval(100000), new ObjectiveFunction(function), GoalType.MINIMIZE, new InitialGuess(values), SimpleBounds.unbounded(dim)); } else if (optimizer == 4) { MultivariateOptimizer search = new CMAESOptimizer(3000000, .05, false, 0, 0, new MersenneTwister(), false, new SimplePointChecker<PointValuePair>(0.5, 0.5)); pair = search.optimize(new MaxEval(30000), new ObjectiveFunction(function), GoalType.MINIMIZE, new InitialGuess(values), new CMAESOptimizer.Sigma(new double[values.length]), new CMAESOptimizer.PopulationSize(1000)); } else if (optimizer == 5) { // 0.01, 0.000001 //2.0D * FastMath.ulp(1.0D), 1e-8 MultivariateOptimizer search = new PowellOptimizer(.05, .05); pair = search.optimize(new InitialGuess(values), new ObjectiveFunction(function), GoalType.MINIMIZE, new MaxEval(100000)); } else if (optimizer == 6) { MultivariateOptimizer search = new PowellOptimizer(1e-7, 1e-7); pair = search.optimize(new InitialGuess(values), new ObjectiveFunction(function), GoalType.MAXIMIZE, new MaxEval(10000)); } return pair.getPoint(); }
From source file:eu.europeana.querylog.learn.Evaluate.java
/** * Performs the CMA learning algorithm./* w ww.j a v a2 s .c om*/ */ public String learningToRankWithCMAES() { startTime = System.currentTimeMillis(); pool = Executors.newFixedThreadPool(properties.getInt("bm25f.learn.concurrency")); try { final int dim = bm25fParams.length; int lambda = 4 + (int) (3. * Math.log(dim)); int maxEvaluations = 1800; double[] inSigma = floatArrayToDouble(getParamsVector(0.5f, 0.5f, 0.25f)); CMAESOptimizer optim = new CMAESOptimizer(maxEvaluations / 10, 1.0, false, 0, 10, new MersenneTwister(), true, new ConvergenceChecker<PointValuePair>() { @Override public boolean converged(int iteration, PointValuePair previous, PointValuePair pv) { maxValue = new Point(doubleArrayToFloat(pv.getFirst()), pv.getSecond().floatValue()); writeLogFile(); logger.info("{}", maxValue); return false; } }); PointValuePair pv = optim.optimize(new MaxEval(maxEvaluations), new ObjectiveFunction(new MultivariateFunction() { @Override public double value(double[] point) { return evaluateAssessments(doubleArrayToFloat(point)); } }), GoalType.MAXIMIZE, new SimpleBounds(floatArrayToDouble(minValues), floatArrayToDouble(maxValues)), new InitialGuess(floatArrayToDouble(bm25fParams)), new CMAESOptimizer.Sigma(inSigma), new CMAESOptimizer.PopulationSize(lambda)); maxValue = new Point(doubleArrayToFloat(pv.getFirst()), pv.getSecond().floatValue()); writeLogFile(); logger.info("final value: \n{}", maxValue); return paramsToXML(); } finally { pool.shutdown(); pool = null; } }
From source file:com.wwidesigner.optimization.ObjectiveFunctionOptimizer.java
protected static PointValuePair runCmaes(BaseObjectiveFunction objective, double[] startPoint) throws TooManyEvaluationsException { // Rely on relative difference to test convergence, // to allow for vast differences in absolute error scale // between objective functions. ConvergenceChecker<PointValuePair> convergenceChecker = new SimpleValueChecker(1.e-6, 1.e-14); MultivariateOptimizer optimizer = new CMAESOptimizer(objective.getMaxEvaluations(), 0.0001 * initialNorm, true, 0, 0, new MersenneTwister(), false, convergenceChecker); return optimizer.optimize(GoalType.MINIMIZE, new ObjectiveFunction(objective), new MaxEval(objective.getMaxEvaluations()), MaxIter.unlimited(), new InitialGuess(startPoint), new SimpleBounds(objective.getLowerBounds(), objective.getUpperBounds()), new CMAESOptimizer.PopulationSize(objective.getNrInterpolations()), new CMAESOptimizer.Sigma(objective.getStdDev())); }
From source file:gdsc.smlm.fitting.nonlinear.MaximumLikelihoodFitter.java
public FitStatus fit(int n, double[] y, double[] y_fit, double[] a, double[] a_dev, double[] error, double noise) { numberOfFittedPoints = n;//from ww w.j a v a2 s.c o m LikelihoodWrapper maximumLikelihoodFunction; // We can use different likelihood wrapper functions: switch (likelihoodFunction) { case POISSON_GAMMA_GAUSSIAN: // Poisson-Gamma-Gaussian - EM-CCD data if (alpha > 0 && sigma > 0) { maximumLikelihoodFunction = new PoissonGammaGaussianLikelihoodWrapper(f, a, y, n, alpha, sigma); break; } case POISSON_GAUSSIAN: // Poisson-Gaussian - CCD data if (sigma > 0) { maximumLikelihoodFunction = new PoissonGaussianLikelihoodWrapper(f, a, y, n, sigma); break; } case POISSON: default: // Poisson - most counting data maximumLikelihoodFunction = new PoissonLikelihoodWrapper(f, a, y, n); } // Check if the method requires the gradient but it cannot be computed if (searchMethod.usesGradient && !maximumLikelihoodFunction.canComputeGradient()) { maximumLikelihoodFunction = new PoissonLikelihoodWrapper(f, a, y, n); } try { double[] startPoint = getInitialSolution(a); PointValuePair optimum = null; if (searchMethod == SearchMethod.POWELL || searchMethod == SearchMethod.POWELL_BOUNDED) { // Non-differentiable version using Powell Optimiser // This is as per the method in Numerical Recipes 10.5 (Direction Set (Powell's) method) // I could extend the optimiser and implement bounds on the directions moved. However the mapping // adapter seems to work OK. final boolean basisConvergence = false; // Perhaps these thresholds should be tighter? // The default is to use the sqrt() of the overall tolerance //final double lineRel = FastMath.sqrt(relativeThreshold); //final double lineAbs = FastMath.sqrt(absoluteThreshold); //final double lineRel = relativeThreshold * 1e2; //final double lineAbs = absoluteThreshold * 1e2; // Since we are fitting only a small number of parameters then just use the same tolerance // for each search direction final double lineRel = relativeThreshold; final double lineAbs = absoluteThreshold; CustomPowellOptimizer o = new CustomPowellOptimizer(relativeThreshold, absoluteThreshold, lineRel, lineAbs, null, basisConvergence); OptimizationData maxIterationData = null; if (getMaxIterations() > 0) maxIterationData = new MaxIter(getMaxIterations()); if (searchMethod == SearchMethod.POWELL) { if (powellFunction == null) { // We must map all the parameters into the same range. This is done in the Mortensen MLE // Python code by using the sqrt of the number of photons and background. if (mapGaussian) { Gaussian2DFunction gf = (Gaussian2DFunction) f; // Re-map signal and background using the sqrt int[] indices = gf.gradientIndices(); int[] map = new int[indices.length]; int count = 0; // Background is always first if (indices[0] == Gaussian2DFunction.BACKGROUND) { map[count++] = 0; } // Look for the Signal in multiple peak 2D Gaussians for (int i = 1; i < indices.length; i++) if (indices[i] % 6 == Gaussian2DFunction.SIGNAL) { map[count++] = i; } if (count > 0) { powellFunction = new MappedMultivariateLikelihood(maximumLikelihoodFunction, Arrays.copyOf(map, count)); } } if (powellFunction == null) { powellFunction = new MultivariateLikelihood(maximumLikelihoodFunction); } } // Update the maximum likelihood function in the Powell function wrapper powellFunction.fun = maximumLikelihoodFunction; OptimizationData positionChecker = null; // new org.apache.commons.math3.optim.PositionChecker(relativeThreshold, absoluteThreshold); if (powellFunction.isMapped()) { MappedMultivariateLikelihood adapter = (MappedMultivariateLikelihood) powellFunction; optimum = o.optimize(maxIterationData, new MaxEval(getMaxEvaluations()), new ObjectiveFunction(powellFunction), GoalType.MINIMIZE, new InitialGuess(adapter.map(startPoint)), positionChecker); double[] solution = adapter.unmap(optimum.getPointRef()); optimum = new PointValuePair(solution, optimum.getValue()); } else { optimum = o.optimize(maxIterationData, new MaxEval(getMaxEvaluations()), new ObjectiveFunction(powellFunction), GoalType.MINIMIZE, new InitialGuess(startPoint), positionChecker); } } else { // Try using the mapping adapter for a bounded Powell search MultivariateFunctionMappingAdapter adapter = new MultivariateFunctionMappingAdapter( new MultivariateLikelihood(maximumLikelihoodFunction), lower, upper); optimum = o.optimize(maxIterationData, new MaxEval(getMaxEvaluations()), new ObjectiveFunction(adapter), GoalType.MINIMIZE, new InitialGuess(adapter.boundedToUnbounded(startPoint))); double[] solution = adapter.unboundedToBounded(optimum.getPointRef()); optimum = new PointValuePair(solution, optimum.getValue()); } iterations = o.getIterations(); evaluations = o.getEvaluations(); } else if (searchMethod == SearchMethod.BOBYQA) { // Differentiable approximation using Powell's BOBYQA algorithm. // This is slower than the Powell optimiser and requires a high number of evaluations. int numberOfInterpolationPoints = this.getNumberOfFittedParameters() + 2; BOBYQAOptimizer o = new BOBYQAOptimizer(numberOfInterpolationPoints); optimum = o.optimize(new MaxEval(getMaxEvaluations()), new ObjectiveFunction(new MultivariateLikelihood(maximumLikelihoodFunction)), GoalType.MINIMIZE, new InitialGuess(startPoint), new SimpleBounds(lower, upper)); iterations = o.getIterations(); evaluations = o.getEvaluations(); } else if (searchMethod == SearchMethod.CMAES) { // TODO - Understand why the CMAES optimiser does not fit very well on test data. It appears // to converge too early and the likelihood scores are not as low as the other optimisers. // CMAESOptimiser based on Matlab code: // https://www.lri.fr/~hansen/cmaes.m // Take the defaults from the Matlab documentation double stopFitness = 0; //Double.NEGATIVE_INFINITY; boolean isActiveCMA = true; int diagonalOnly = 0; int checkFeasableCount = 1; RandomGenerator random = new Well19937c(); boolean generateStatistics = false; // The sigma determines the search range for the variables. It should be 1/3 of the initial search region. double[] sigma = new double[lower.length]; for (int i = 0; i < sigma.length; i++) sigma[i] = (upper[i] - lower[i]) / 3; int popSize = (int) (4 + Math.floor(3 * Math.log(sigma.length))); // The CMAES optimiser is random and restarting can overcome problems with quick convergence. // The Apache commons documentations states that convergence should occur between 30N and 300N^2 // function evaluations final int n30 = FastMath.min(sigma.length * sigma.length * 30, getMaxEvaluations() / 2); evaluations = 0; OptimizationData[] data = new OptimizationData[] { new InitialGuess(startPoint), new CMAESOptimizer.PopulationSize(popSize), new MaxEval(getMaxEvaluations()), new CMAESOptimizer.Sigma(sigma), new ObjectiveFunction(new MultivariateLikelihood(maximumLikelihoodFunction)), GoalType.MINIMIZE, new SimpleBounds(lower, upper) }; // Iterate to prevent early convergence int repeat = 0; while (evaluations < n30) { if (repeat++ > 1) { // Update the start point and population size data[0] = new InitialGuess(optimum.getPointRef()); popSize *= 2; data[1] = new CMAESOptimizer.PopulationSize(popSize); } CMAESOptimizer o = new CMAESOptimizer(getMaxIterations(), stopFitness, isActiveCMA, diagonalOnly, checkFeasableCount, random, generateStatistics, new SimpleValueChecker(relativeThreshold, absoluteThreshold)); PointValuePair result = o.optimize(data); iterations += o.getIterations(); evaluations += o.getEvaluations(); //System.out.printf("CMAES [%d] i=%d [%d], e=%d [%d]\n", repeat, o.getIterations(), iterations, // o.getEvaluations(), totalEvaluations); if (optimum == null || result.getValue() < optimum.getValue()) { optimum = result; } } } else if (searchMethod == SearchMethod.BFGS) { // BFGS can use an approximate line search minimisation where as Powell and conjugate gradient // methods require a more accurate line minimisation. The BFGS search does not do a full // minimisation but takes appropriate steps in the direction of the current gradient. // Do not use the convergence checker on the value of the function. Use the convergence on the // point coordinate and gradient //BFGSOptimizer o = new BFGSOptimizer(new SimpleValueChecker(rel, abs)); BFGSOptimizer o = new BFGSOptimizer(); // Configure maximum step length for each dimension using the bounds double[] stepLength = new double[lower.length]; for (int i = 0; i < stepLength.length; i++) { stepLength[i] = (upper[i] - lower[i]) * 0.3333333; if (stepLength[i] <= 0) stepLength[i] = Double.POSITIVE_INFINITY; } // The GoalType is always minimise so no need to pass this in OptimizationData positionChecker = null; //new org.apache.commons.math3.optim.PositionChecker(relativeThreshold, absoluteThreshold); optimum = o.optimize(new MaxEval(getMaxEvaluations()), new ObjectiveFunctionGradient(new MultivariateVectorLikelihood(maximumLikelihoodFunction)), new ObjectiveFunction(new MultivariateLikelihood(maximumLikelihoodFunction)), new InitialGuess(startPoint), new SimpleBounds(lowerConstraint, upperConstraint), new BFGSOptimizer.GradientTolerance(relativeThreshold), positionChecker, new BFGSOptimizer.StepLength(stepLength)); iterations = o.getIterations(); evaluations = o.getEvaluations(); } else { // The line search algorithm often fails. This is due to searching into a region where the // function evaluates to a negative so has been clipped. This means the upper bound of the line // cannot be found. // Note that running it on an easy problem (200 photons with fixed fitting (no background)) the algorithm // does sometimes produces results better than the Powell algorithm but it is slower. BoundedNonLinearConjugateGradientOptimizer o = new BoundedNonLinearConjugateGradientOptimizer( (searchMethod == SearchMethod.CONJUGATE_GRADIENT_FR) ? Formula.FLETCHER_REEVES : Formula.POLAK_RIBIERE, new SimpleValueChecker(relativeThreshold, absoluteThreshold)); // Note: The gradients may become unstable at the edge of the bounds. Or they will not change // direction if the true solution is on the bounds since the gradient will always continue // towards the bounds. This is key to the conjugate gradient method. It searches along a vector // until the direction of the gradient is in the opposite direction (using dot products, i.e. // cosine of angle between them) // NR 10.7 states there is no advantage of the variable metric DFP or BFGS methods over // conjugate gradient methods. So I will try these first. // Try this: // Adapt the conjugate gradient optimiser to use the gradient to pick the search direction // and then for the line minimisation. However if the function is out of bounds then clip the // variables at the bounds and continue. // If the current point is at the bounds and the gradient is to continue out of bounds then // clip the gradient too. // Or: just use the gradient for the search direction then use the line minimisation/rest // as per the Powell optimiser. The bounds should limit the search. // I tried a Bounded conjugate gradient optimiser with clipped variables: // This sometimes works. However when the variables go a long way out of the expected range the gradients // can have vastly different magnitudes. This results in the algorithm stalling since the gradients // can be close to zero and the some of the parameters are no longer adjusted. // Perhaps this can be looked for and the algorithm then gives up and resorts to a Powell optimiser from // the current point. // Changed the bracketing step to very small (default is 1, changed to 0.001). This improves the // performance. The gradient direction is very sensitive to small changes in the coordinates so a // tighter bracketing of the line search helps. // Tried using a non-gradient method for the line search copied from the Powell optimiser: // This also works when the bracketing step is small but the number of iterations is higher. // 24.10.2014: I have tried to get conjugate gradient to work but the gradient function // must not behave suitably for the optimiser. In the current state both methods of using a // Bounded Conjugate Gradient Optimiser perform poorly relative to other optimisers: // Simulated : n=1000, signal=200, x=0.53, y=0.47 // LVM : n=1000, signal=171, x=0.537, y=0.471 (1.003s) // Powell : n=1000, signal=187, x=0.537, y=0.48 (1.238s) // Gradient based PR (constrained): n=858, signal=161, x=0.533, y=0.474 (2.54s) // Gradient based PR (bounded): n=948, signal=161, x=0.533, y=0.473 (2.67s) // Non-gradient based : n=1000, signal=151.47, x=0.535, y=0.474 (1.626s) // The conjugate optimisers are slower, under predict the signal by the most and in the case of // the gradient based optimiser, fail to converge on some problems. This is worse when constrained // fitting is used and not tightly bounded fitting. // I will leave the code in as an option but would not recommend using it. I may remove it in the // future. // Note: It is strange that the non-gradient based line minimisation is more successful. // It may be that the gradient function is not accurate (due to round off error) or that it is // simply wrong when far from the optimum. My JUnit tests only evaluate the function within the // expected range of the answer. // Note the default step size on the Powell optimiser is 1 but the initial directions are unit vectors. // So our bracketing step should be a minimum of 1 / average length of the first gradient vector to prevent // the first step being too large when bracketing. final double gradient[] = new double[startPoint.length]; maximumLikelihoodFunction.likelihood(startPoint, gradient); double l = 0; for (double d : gradient) l += d * d; final double bracketingStep = FastMath.min(0.001, ((l > 1) ? 1.0 / l : 1)); //System.out.printf("Bracketing step = %f (length=%f)\n", bracketingStep, l); o.setUseGradientLineSearch(gradientLineMinimisation); optimum = o.optimize(new MaxEval(getMaxEvaluations()), new ObjectiveFunctionGradient(new MultivariateVectorLikelihood(maximumLikelihoodFunction)), new ObjectiveFunction(new MultivariateLikelihood(maximumLikelihoodFunction)), GoalType.MINIMIZE, new InitialGuess(startPoint), new SimpleBounds(lowerConstraint, upperConstraint), new BoundedNonLinearConjugateGradientOptimizer.BracketingStep(bracketingStep)); iterations = o.getIterations(); evaluations = o.getEvaluations(); //maximumLikelihoodFunction.value(solution, gradient); //System.out.printf("Iter = %d, %g @ %s : %s\n", iterations, ll, Arrays.toString(solution), // Arrays.toString(gradient)); } final double[] solution = optimum.getPointRef(); setSolution(a, solution); //System.out.printf("Iter = %d, Eval = %d, %g @ %s\n", iterations, evaluations, optimum.getValue(), // java.util.Arrays.toString(solution)); // Compute residuals for the FunctionSolver interface if (y_fit == null || y_fit.length < n) y_fit = new double[n]; f.initialise(a); residualSumOfSquares = 0; for (int i = 0; i < n; i++) { y_fit[i] = f.eval(i); final double residual = y[i] - y_fit[i]; residualSumOfSquares += residual * residual; } if (a_dev != null) { // Assume the Maximum Likelihood estimator returns the optimum fit (achieves the Cramer Roa // lower bounds) and so the covariance can be obtained from the Fisher Information Matrix. final int[] gradientIndices = f.gradientIndices(); final int nparams = gradientIndices.length; GradientCalculator calculator = GradientCalculatorFactory.newCalculator(nparams); final double[] I = calculator.fisherInformationDiagonal(n, a, f); for (int i = 0; i < gradientIndices.length; i++) a_dev[gradientIndices[i]] = 1.0 / Math.sqrt(I[i]); } error[0] = NonLinearFit.getError(residualSumOfSquares, noise, n, f.gradientIndices().length); totalSumOfSquares = getSumOfSquares(n, y); } catch (TooManyIterationsException e) { //System.out.printf("Too many iterations = %d\n", e.getMax()); //e.printStackTrace(); return FitStatus.FAILED_TO_CONVERGE; } catch (TooManyEvaluationsException e) { //System.out.printf("Too many evaluations = %d\n", e.getMax()); //e.printStackTrace(); return FitStatus.FAILED_TO_CONVERGE; } catch (ConvergenceException e) { // Occurs when QR decomposition fails - mark as a singular non-linear model (no solution) //System.out.printf("Singular non linear model = %s\n", e.getMessage()); return FitStatus.SINGULAR_NON_LINEAR_MODEL; } catch (BFGSOptimizer.LineSearchRoundoffException e) { //System.out.println("BFGS error: " + e.getMessage()); //e.printStackTrace(); return FitStatus.FAILED_TO_CONVERGE; } catch (Exception e) { //System.out.printf("Unknown error = %s\n", e.getMessage()); e.printStackTrace(); return FitStatus.UNKNOWN; } return FitStatus.OK; }
From source file:gdsc.smlm.ij.plugins.pcpalm.PCPALMFitting.java
private PointValuePair runBoundedOptimiser(double[][] gr, double[] initialSolution, double[] lB, double[] uB, SumOfSquaresModelFunction function) { // Create the functions to optimise ObjectiveFunction objective = new ObjectiveFunction(new SumOfSquaresMultivariateFunction(function)); ObjectiveFunctionGradient gradient = new ObjectiveFunctionGradient( new SumOfSquaresMultivariateVectorFunction(function)); final boolean debug = false; // Try a BFGS optimiser since this will produce a deterministic solution and can respect bounds. PointValuePair optimum = null;// ww w . j a v a 2 s. c o m boundedEvaluations = 0; final MaxEval maxEvaluations = new MaxEval(2000); MultivariateOptimizer opt = null; for (int iteration = 0; iteration <= fitRestarts; iteration++) { try { opt = new BFGSOptimizer(); final double relativeThreshold = 1e-6; // Configure maximum step length for each dimension using the bounds double[] stepLength = new double[lB.length]; for (int i = 0; i < stepLength.length; i++) stepLength[i] = (uB[i] - lB[i]) * 0.3333333; // The GoalType is always minimise so no need to pass this in optimum = opt.optimize(maxEvaluations, gradient, objective, new InitialGuess((optimum == null) ? initialSolution : optimum.getPointRef()), new SimpleBounds(lB, uB), new BFGSOptimizer.GradientTolerance(relativeThreshold), new BFGSOptimizer.StepLength(stepLength)); if (debug) System.out.printf("BFGS Iter %d = %g (%d)\n", iteration, optimum.getValue(), opt.getEvaluations()); } catch (TooManyEvaluationsException e) { break; // No need to restart } catch (RuntimeException e) { break; // No need to restart } finally { boundedEvaluations += opt.getEvaluations(); } } // Try a CMAES optimiser which is non-deterministic. To overcome this we perform restarts. // CMAESOptimiser based on Matlab code: // https://www.lri.fr/~hansen/cmaes.m // Take the defaults from the Matlab documentation double stopFitness = 0; //Double.NEGATIVE_INFINITY; boolean isActiveCMA = true; int diagonalOnly = 0; int checkFeasableCount = 1; RandomGenerator random = new Well44497b(); //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[] range = new double[lB.length]; for (int i = 0; i < lB.length; i++) range[i] = (uB[i] - lB[i]) / 3; OptimizationData sigma = new CMAESOptimizer.Sigma(range); OptimizationData popSize = new CMAESOptimizer.PopulationSize( (int) (4 + Math.floor(3 * Math.log(initialSolution.length)))); SimpleBounds bounds = new SimpleBounds(lB, uB); opt = new CMAESOptimizer(maxEvaluations.getMaxEval(), stopFitness, isActiveCMA, diagonalOnly, checkFeasableCount, random, generateStatistics, checker); // Restart the optimiser several times and take the best answer. for (int iteration = 0; iteration <= fitRestarts; iteration++) { try { // Start from the initial solution PointValuePair constrainedSolution = opt.optimize(new InitialGuess(initialSolution), objective, GoalType.MINIMIZE, bounds, sigma, popSize, maxEvaluations); if (debug) System.out.printf("CMAES Iter %d initial = %g (%d)\n", iteration, constrainedSolution.getValue(), opt.getEvaluations()); boundedEvaluations += opt.getEvaluations(); if (optimum == null || constrainedSolution.getValue() < optimum.getValue()) { optimum = constrainedSolution; } } catch (TooManyEvaluationsException e) { } catch (TooManyIterationsException e) { } finally { boundedEvaluations += maxEvaluations.getMaxEval(); } if (optimum == null) continue; try { // Also restart from the current optimum PointValuePair constrainedSolution = opt.optimize(new InitialGuess(optimum.getPointRef()), objective, GoalType.MINIMIZE, bounds, sigma, popSize, maxEvaluations); if (debug) System.out.printf("CMAES Iter %d restart = %g (%d)\n", iteration, constrainedSolution.getValue(), opt.getEvaluations()); if (constrainedSolution.getValue() < optimum.getValue()) { optimum = constrainedSolution; } } catch (TooManyEvaluationsException e) { } catch (TooManyIterationsException e) { } finally { boundedEvaluations += maxEvaluations.getMaxEval(); } } return optimum; }
From source file:gdsc.smlm.ij.plugins.TraceDiffusion.java
/** * Fit the jump distance histogram using a cumulative sum as detailed in * <p>/* w w w . jav a2 s. com*/ * 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:norbert.mynemo.core.selection.SvdRecommenderSelector.java
/** * Evaluates the given recommender type for several configurations. The number of iterations and * the number of features for each evaluation is chosen via an optimizer. The number of evaluation * performed depends of the number of steps necessary for the evaluation results to converge. The * number of evaluations performed is also bounded by an internal maximum. * * <p>// w w w. java 2 s . co m * The given recommender type must be part of the SVD based family. * * @return all evaluations done during the selection process */ public Collection<RecommenderEvaluation> select(RecommenderType type, double minimumCoverage) { checkArgument(type.getFamily() == RecommenderFamily.SVD_BASED); // initialize necessary optimization data ConvergenceChecker<PointValuePair> checker = new MaxIterationChecker<PointValuePair>(CMAES_MAX_ITERATIONS); SvdRecommenderEvalFunction function = new SvdRecommenderEvalFunction(configuration, type, minimumCoverage); ObjectiveFunction objectiveFunction = new ObjectiveFunction(function); CMAESOptimizer optimizer = new CMAESOptimizer(CMAES_MAX_ITERATIONS, 1.0, true, 2, 0, new JDKRandomGenerator(), false, checker); MaxEval maxEval = MaxEval.unlimited(); GoalType goalType = GoalType.MINIMIZE; SimpleBounds bounds = new SimpleBounds(new double[] { MIN_SVD_FEATURES, MIN_SVD_ITERATIONS }, new double[] { MAX_SVD_FEATURES, MAX_SVD_ITERATIONS }); InitialGuess initialGuess = new InitialGuess( new double[] { INITIAL_GUESS_FEATURES, INITIAL_GUESS_ITERATIONS }); CMAESOptimizer.PopulationSize populationSize = new CMAESOptimizer.PopulationSize(16); CMAESOptimizer.Sigma sigma = new CMAESOptimizer.Sigma( new double[] { CMAES_SIGMA_FEATURE, CMAES_SIGMA_ITERATION }); // run the optimizer optimizer.optimize(objectiveFunction, goalType, initialGuess, populationSize, sigma, bounds, maxEval); return function.getEvaluations(); }