use of org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.Optimum in project GDSC-SMLM by aherbert.
the class PCPALMFitting method fitEmulsionModel.
/**
* Fits the correlation curve with r>0 to the clustered model using the estimated density and precision. Parameters
* must be fit within a tolerance of the starting values.
*
* @param gr
* @param sigmaS
* The estimated precision
* @param proteinDensity
* The estimated protein density
* @return The fitted parameters [precision, density, clusterRadius, clusterDensity]
*/
private double[] fitEmulsionModel(double[][] gr, double sigmaS, double proteinDensity, String resultColour) {
final EmulsionModelFunctionGradient function = new EmulsionModelFunctionGradient();
emulsionModel = function;
log("Fitting %s: Estimated precision = %f nm, estimated protein density = %g um^-2", emulsionModel.getName(), sigmaS, proteinDensity * 1e6);
emulsionModel.setLogging(true);
for (int i = offset; i < gr[0].length; i++) {
// Only fit the curve above the estimated resolution (points below it will be subject to error)
if (gr[0][i] > sigmaS * fitAboveEstimatedPrecision)
emulsionModel.addPoint(gr[0][i], gr[1][i]);
}
double[] parameters;
// The model is: sigma, density, range, amplitude, alpha
double[] initialSolution = new double[] { sigmaS, proteinDensity, sigmaS * 5, 1, sigmaS * 5 };
int evaluations = 0;
// Constrain the fitting to be close to the estimated precision (sigmaS) and protein density.
// LVM fitting does not support constrained fitting so use a bounded optimiser.
SumOfSquaresModelFunction emulsionModelMulti = new SumOfSquaresModelFunction(emulsionModel);
double[] x = emulsionModelMulti.x;
double[] y = emulsionModelMulti.y;
// Range should be equal to the first time the g(r) curve crosses 1
for (int i = 0; i < x.length; i++) if (y[i] < 1) {
initialSolution[4] = initialSolution[2] = (i > 0) ? (x[i - 1] + x[i]) * 0.5 : x[i];
break;
}
// Put some bounds around the initial guess. Use the fitting tolerance (in %) if provided.
double limit = (fittingTolerance > 0) ? 1 + fittingTolerance / 100 : 2;
double[] lB = new double[] { initialSolution[0] / limit, initialSolution[1] / limit, 0, 0, 0 };
// The amplitude and range should not extend beyond the limits of the g(r) curve.
// TODO - Find out the expected range for the alpha parameter.
double[] uB = new double[] { initialSolution[0] * limit, initialSolution[1] * limit, Maths.max(x), Maths.max(gr[1]), Maths.max(x) * 2 };
log("Fitting %s using a bounded search: %s < precision < %s & %s < density < %s", emulsionModel.getName(), Utils.rounded(lB[0], 4), Utils.rounded(uB[0], 4), Utils.rounded(lB[1] * 1e6, 4), Utils.rounded(uB[1] * 1e6, 4));
PointValuePair constrainedSolution = runBoundedOptimiser(gr, initialSolution, lB, uB, emulsionModelMulti);
if (constrainedSolution == null)
return null;
parameters = constrainedSolution.getPointRef();
evaluations = boundedEvaluations;
// Refit using a LVM
if (useLSE) {
log("Re-fitting %s using a gradient optimisation", emulsionModel.getName());
LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
Optimum lvmSolution;
try {
//@formatter:off
LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(parameters).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, new MultivariateMatrixFunction() {
public double[][] value(double[] point) throws IllegalArgumentException {
return function.jacobian(point);
}
}).build();
//@formatter:on
lvmSolution = optimizer.optimize(problem);
evaluations += lvmSolution.getEvaluations();
double ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
if (ss < constrainedSolution.getValue()) {
log("Re-fitting %s improved the SS from %s to %s (-%s%%)", emulsionModel.getName(), Utils.rounded(constrainedSolution.getValue(), 4), Utils.rounded(ss, 4), Utils.rounded(100 * (constrainedSolution.getValue() - ss) / constrainedSolution.getValue(), 4));
parameters = lvmSolution.getPoint().toArray();
}
} catch (TooManyIterationsException e) {
log("Failed to re-fit %s: Too many iterations (%s)", emulsionModel.getName(), e.getMessage());
} catch (ConvergenceException e) {
log("Failed to re-fit %s: %s", emulsionModel.getName(), e.getMessage());
}
}
emulsionModel.setLogging(false);
// Ensure the width is positive
parameters[0] = Math.abs(parameters[0]);
//parameters[2] = Math.abs(parameters[2]);
double ss = 0;
double[] obs = emulsionModel.getY();
double[] exp = emulsionModel.value(parameters);
for (int i = 0; i < obs.length; i++) ss += (obs[i] - exp[i]) * (obs[i] - exp[i]);
ic3 = Maths.getAkaikeInformationCriterionFromResiduals(ss, emulsionModel.size(), parameters.length);
final double fitSigmaS = parameters[0];
final double fitProteinDensity = parameters[1];
//The radius of the cluster domain
final double domainRadius = parameters[2];
//The density of the cluster domain
final double domainDensity = parameters[3];
//The coherence length between circles
final double coherence = parameters[4];
// This is from the PC-PALM paper. It may not be correct for the emulsion model.
final double nCluster = 2 * domainDensity * Math.PI * domainRadius * domainRadius * fitProteinDensity;
double e1 = parameterDrift(sigmaS, fitSigmaS);
double e2 = parameterDrift(proteinDensity, fitProteinDensity);
log(" %s fit: SS = %f. cAIC = %f. %d evaluations", emulsionModel.getName(), ss, ic3, evaluations);
log(" %s parameters:", emulsionModel.getName());
log(" Average precision = %s nm (%s%%)", Utils.rounded(fitSigmaS, 4), Utils.rounded(e1, 4));
log(" Average protein density = %s um^-2 (%s%%)", Utils.rounded(fitProteinDensity * 1e6, 4), Utils.rounded(e2, 4));
log(" Domain radius = %s nm", Utils.rounded(domainRadius, 4));
log(" Domain density = %s", Utils.rounded(domainDensity, 4));
log(" Domain coherence = %s", Utils.rounded(coherence, 4));
log(" nCluster = %s", Utils.rounded(nCluster, 4));
// Check the fitted parameters are within tolerance of the initial estimates
valid2 = true;
if (fittingTolerance > 0 && (Math.abs(e1) > fittingTolerance || Math.abs(e2) > fittingTolerance)) {
log(" Failed to fit %s within tolerance (%s%%): Average precision = %f nm (%s%%), average protein density = %g um^-2 (%s%%)", emulsionModel.getName(), Utils.rounded(fittingTolerance, 4), fitSigmaS, Utils.rounded(e1, 4), fitProteinDensity * 1e6, Utils.rounded(e2, 4));
valid2 = false;
}
// Check extra parameters. Domain radius should be higher than the precision. Density should be positive
if (domainRadius < fitSigmaS) {
log(" Failed to fit %s: Domain radius is smaller than the average precision (%s < %s)", emulsionModel.getName(), Utils.rounded(domainRadius, 4), Utils.rounded(fitSigmaS, 4));
valid2 = false;
}
if (domainDensity < 0) {
log(" Failed to fit %s: Domain density is negative (%s)", emulsionModel.getName(), Utils.rounded(domainDensity, 4));
valid2 = false;
}
if (ic3 > ic1) {
log(" Failed to fit %s - Information Criterion has increased %s%%", emulsionModel.getName(), Utils.rounded((100 * (ic3 - ic1) / ic1), 4));
valid2 = false;
}
addResult(emulsionModel.getName(), resultColour, valid2, fitSigmaS, fitProteinDensity, domainRadius, domainDensity, nCluster, coherence, ic3);
return parameters;
}
use of org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.Optimum in project GDSC-SMLM by aherbert.
the class BlinkEstimator method fit.
/**
* Fit the dark time to counts of molecules curve. Only use the first n fitted points.
* <p>
* Calculates:<br/>
* N = The number of photoblinking molecules in the sample<br/>
* nBlink = The average number of blinks per flourophore<br/>
* tOff = The off-time
*
* @param td
* The dark time
* @param ntd
* The counts of molecules
* @param nFittedPoints
* @param log
* Write the fitting results to the ImageJ log window
* @return The fitted parameters [N, nBlink, tOff], or null if no fit was possible
*/
public double[] fit(double[] td, double[] ntd, int nFittedPoints, boolean log) {
blinkingModel = new BlinkingFunction();
blinkingModel.setLogging(true);
for (int i = 0; i < nFittedPoints; i++) blinkingModel.addPoint(td[i], ntd[i]);
// Different convergence thresholds seem to have no effect on the resulting fit, only the number of
// iterations for convergence
double initialStepBoundFactor = 100;
double costRelativeTolerance = 1e-6;
double parRelativeTolerance = 1e-6;
double orthoTolerance = 1e-6;
double threshold = Precision.SAFE_MIN;
LevenbergMarquardtOptimizer optimiser = new LevenbergMarquardtOptimizer(initialStepBoundFactor, costRelativeTolerance, parRelativeTolerance, orthoTolerance, threshold);
try {
double[] obs = blinkingModel.getY();
//@formatter:off
LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(1000).start(new double[] { ntd[0], 0.1, td[1] }).target(obs).weight(new DiagonalMatrix(blinkingModel.getWeights())).model(blinkingModel, new MultivariateMatrixFunction() {
public double[][] value(double[] point) throws IllegalArgumentException {
return blinkingModel.jacobian(point);
}
}).build();
//@formatter:on
blinkingModel.setLogging(false);
Optimum optimum = optimiser.optimize(problem);
double[] parameters = optimum.getPoint().toArray();
//double[] exp = blinkingModel.value(parameters);
double mean = 0;
for (double d : obs) mean += d;
mean /= obs.length;
double ssResiduals = 0, ssTotal = 0;
for (int i = 0; i < obs.length; i++) {
//ssResiduals += (obs[i] - exp[i]) * (obs[i] - exp[i]);
ssTotal += (obs[i] - mean) * (obs[i] - mean);
}
// This is true if the weights are 1
ssResiduals = optimum.getResiduals().dotProduct(optimum.getResiduals());
r2 = 1 - ssResiduals / ssTotal;
adjustedR2 = getAdjustedCoefficientOfDetermination(ssResiduals, ssTotal, obs.length, parameters.length);
if (log) {
Utils.log(" Fit %d points. R^2 = %s. Adjusted R^2 = %s", obs.length, Utils.rounded(r2, 4), Utils.rounded(adjustedR2, 4));
Utils.log(" N=%s, nBlink=%s, tOff=%s (%s frames)", Utils.rounded(parameters[0], 4), Utils.rounded(parameters[1], 4), Utils.rounded(parameters[2], 4), Utils.rounded(parameters[2] / msPerFrame, 4));
}
return parameters;
} catch (TooManyIterationsException e) {
if (log)
Utils.log(" Failed to fit %d points: Too many iterations: (%s)", blinkingModel.size(), e.getMessage());
return null;
} catch (ConvergenceException e) {
if (log)
Utils.log(" Failed to fit %d points", blinkingModel.size());
return null;
}
}
use of org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.Optimum in project GDSC-SMLM by aherbert.
the class PCPALMFitting method fitRandomModel.
/**
* Fits the correlation curve with r>0 to the random model using the estimated density and precision. Parameters
* must be fit within a tolerance of the starting values.
*
* @param gr
* @param sigmaS
* The estimated precision
* @param proteinDensity
* The estimate protein density
* @return The fitted parameters [precision, density]
*/
private double[] fitRandomModel(double[][] gr, double sigmaS, double proteinDensity, String resultColour) {
final RandomModelFunction function = new RandomModelFunction();
randomModel = function;
log("Fitting %s: Estimated precision = %f nm, estimated protein density = %g um^-2", randomModel.getName(), sigmaS, proteinDensity * 1e6);
randomModel.setLogging(true);
for (int i = offset; i < gr[0].length; i++) {
// Only fit the curve above the estimated resolution (points below it will be subject to error)
if (gr[0][i] > sigmaS * fitAboveEstimatedPrecision)
randomModel.addPoint(gr[0][i], gr[1][i]);
}
LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
Optimum optimum;
try {
//@formatter:off
LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(new double[] { sigmaS, proteinDensity }).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, new MultivariateMatrixFunction() {
public double[][] value(double[] point) throws IllegalArgumentException {
return function.jacobian(point);
}
}).build();
//@formatter:on
optimum = optimizer.optimize(problem);
} catch (TooManyIterationsException e) {
log("Failed to fit %s: Too many iterations (%s)", randomModel.getName(), e.getMessage());
return null;
} catch (ConvergenceException e) {
log("Failed to fit %s: %s", randomModel.getName(), e.getMessage());
return null;
}
randomModel.setLogging(false);
double[] parameters = optimum.getPoint().toArray();
// Ensure the width is positive
parameters[0] = Math.abs(parameters[0]);
double ss = optimum.getResiduals().dotProduct(optimum.getResiduals());
ic1 = Maths.getAkaikeInformationCriterionFromResiduals(ss, randomModel.size(), parameters.length);
final double fitSigmaS = parameters[0];
final double fitProteinDensity = parameters[1];
// Check the fitted parameters are within tolerance of the initial estimates
double e1 = parameterDrift(sigmaS, fitSigmaS);
double e2 = parameterDrift(proteinDensity, fitProteinDensity);
log(" %s fit: SS = %f. cAIC = %f. %d evaluations", randomModel.getName(), ss, ic1, optimum.getEvaluations());
log(" %s parameters:", randomModel.getName());
log(" Average precision = %s nm (%s%%)", Utils.rounded(fitSigmaS, 4), Utils.rounded(e1, 4));
log(" Average protein density = %s um^-2 (%s%%)", Utils.rounded(fitProteinDensity * 1e6, 4), Utils.rounded(e2, 4));
valid1 = true;
if (fittingTolerance > 0 && (Math.abs(e1) > fittingTolerance || Math.abs(e2) > fittingTolerance)) {
log(" Failed to fit %s within tolerance (%s%%): Average precision = %f nm (%s%%), average protein density = %g um^-2 (%s%%)", randomModel.getName(), Utils.rounded(fittingTolerance, 4), fitSigmaS, Utils.rounded(e1, 4), fitProteinDensity * 1e6, Utils.rounded(e2, 4));
valid1 = false;
}
if (valid1) {
// ---------
// TODO - My data does not comply with this criteria.
// This could be due to the PC-PALM Molecule code limiting the nmPerPixel to fit the images in memory
// thus removing correlations at small r.
// It could also be due to the nature of the random simulations being 3D not 2D membranes
// as per the PC-PALM paper.
// ---------
// Evaluate g(r)protein where:
// g(r)peaks = g(r)protein + g(r)stoch
// g(r)peaks ~ 1 + g(r)stoch
// Verify g(r)protein should be <1.5 for all r>0
double[] gr_stoch = randomModel.value(parameters);
double[] gr_peaks = randomModel.getY();
double[] gr_ = randomModel.getX();
//SummaryStatistics stats = new SummaryStatistics();
for (int i = 0; i < gr_peaks.length; i++) {
// Only evaluate above the fitted average precision
if (gr_[i] < fitSigmaS)
continue;
// Note the RandomModelFunction evaluates g(r)stoch + 1;
double gr_protein_i = gr_peaks[i] - (gr_stoch[i] - 1);
if (gr_protein_i > gr_protein_threshold) {
// Failed fit
log(" Failed to fit %s: g(r)protein %s > %s @ r=%s", randomModel.getName(), Utils.rounded(gr_protein_i, 4), Utils.rounded(gr_protein_threshold, 4), Utils.rounded(gr_[i], 4));
valid1 = false;
}
//stats.addValue(gr_i);
//System.out.printf("g(r)protein @ %f = %f\n", gr[0][i], gr_protein_i);
}
}
addResult(randomModel.getName(), resultColour, valid1, fitSigmaS, fitProteinDensity, 0, 0, 0, 0, ic1);
return parameters;
}
use of org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.Optimum in project GDSC-SMLM by aherbert.
the class PCPALMFitting method fitClusteredModel.
/**
* Fits the correlation curve with r>0 to the clustered model using the estimated density and precision. Parameters
* must be fit within a tolerance of the starting values.
*
* @param gr
* @param sigmaS
* The estimated precision
* @param proteinDensity
* The estimated protein density
* @return The fitted parameters [precision, density, clusterRadius, clusterDensity]
*/
private double[] fitClusteredModel(double[][] gr, double sigmaS, double proteinDensity, String resultColour) {
final ClusteredModelFunctionGradient function = new ClusteredModelFunctionGradient();
clusteredModel = function;
log("Fitting %s: Estimated precision = %f nm, estimated protein density = %g um^-2", clusteredModel.getName(), sigmaS, proteinDensity * 1e6);
clusteredModel.setLogging(true);
for (int i = offset; i < gr[0].length; i++) {
// Only fit the curve above the estimated resolution (points below it will be subject to error)
if (gr[0][i] > sigmaS * fitAboveEstimatedPrecision)
clusteredModel.addPoint(gr[0][i], gr[1][i]);
}
double[] parameters;
// The model is: sigma, density, range, amplitude
double[] initialSolution = new double[] { sigmaS, proteinDensity, sigmaS * 5, 1 };
int evaluations = 0;
// Constrain the fitting to be close to the estimated precision (sigmaS) and protein density.
// LVM fitting does not support constrained fitting so use a bounded optimiser.
SumOfSquaresModelFunction clusteredModelMulti = new SumOfSquaresModelFunction(clusteredModel);
double[] x = clusteredModelMulti.x;
// Put some bounds around the initial guess. Use the fitting tolerance (in %) if provided.
double limit = (fittingTolerance > 0) ? 1 + fittingTolerance / 100 : 2;
double[] lB = new double[] { initialSolution[0] / limit, initialSolution[1] / limit, 0, 0 };
// The amplitude and range should not extend beyond the limits of the g(r) curve.
double[] uB = new double[] { initialSolution[0] * limit, initialSolution[1] * limit, Maths.max(x), Maths.max(gr[1]) };
log("Fitting %s using a bounded search: %s < precision < %s & %s < density < %s", clusteredModel.getName(), Utils.rounded(lB[0], 4), Utils.rounded(uB[0], 4), Utils.rounded(lB[1] * 1e6, 4), Utils.rounded(uB[1] * 1e6, 4));
PointValuePair constrainedSolution = runBoundedOptimiser(gr, initialSolution, lB, uB, clusteredModelMulti);
if (constrainedSolution == null)
return null;
parameters = constrainedSolution.getPointRef();
evaluations = boundedEvaluations;
// Refit using a LVM
if (useLSE) {
log("Re-fitting %s using a gradient optimisation", clusteredModel.getName());
LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
Optimum lvmSolution;
try {
//@formatter:off
LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(3000).start(parameters).target(function.getY()).weight(new DiagonalMatrix(function.getWeights())).model(function, new MultivariateMatrixFunction() {
public double[][] value(double[] point) throws IllegalArgumentException {
return function.jacobian(point);
}
}).build();
//@formatter:on
lvmSolution = optimizer.optimize(problem);
evaluations += lvmSolution.getEvaluations();
double ss = lvmSolution.getResiduals().dotProduct(lvmSolution.getResiduals());
if (ss < constrainedSolution.getValue()) {
log("Re-fitting %s improved the SS from %s to %s (-%s%%)", clusteredModel.getName(), Utils.rounded(constrainedSolution.getValue(), 4), Utils.rounded(ss, 4), Utils.rounded(100 * (constrainedSolution.getValue() - ss) / constrainedSolution.getValue(), 4));
parameters = lvmSolution.getPoint().toArray();
}
} catch (TooManyIterationsException e) {
log("Failed to re-fit %s: Too many iterations (%s)", clusteredModel.getName(), e.getMessage());
} catch (ConvergenceException e) {
log("Failed to re-fit %s: %s", clusteredModel.getName(), e.getMessage());
}
}
clusteredModel.setLogging(false);
// Ensure the width is positive
parameters[0] = Math.abs(parameters[0]);
//parameters[2] = Math.abs(parameters[2]);
double ss = 0;
double[] obs = clusteredModel.getY();
double[] exp = clusteredModel.value(parameters);
for (int i = 0; i < obs.length; i++) ss += (obs[i] - exp[i]) * (obs[i] - exp[i]);
ic2 = Maths.getAkaikeInformationCriterionFromResiduals(ss, clusteredModel.size(), parameters.length);
final double fitSigmaS = parameters[0];
final double fitProteinDensity = parameters[1];
//The radius of the cluster domain
final double domainRadius = parameters[2];
//The density of the cluster domain
final double domainDensity = parameters[3];
// This is from the PC-PALM paper. However that paper fits the g(r)protein exponential convolved in 2D
// with the g(r)PSF. In this method we have just fit the exponential
final double nCluster = 2 * domainDensity * Math.PI * domainRadius * domainRadius * fitProteinDensity;
double e1 = parameterDrift(sigmaS, fitSigmaS);
double e2 = parameterDrift(proteinDensity, fitProteinDensity);
log(" %s fit: SS = %f. cAIC = %f. %d evaluations", clusteredModel.getName(), ss, ic2, evaluations);
log(" %s parameters:", clusteredModel.getName());
log(" Average precision = %s nm (%s%%)", Utils.rounded(fitSigmaS, 4), Utils.rounded(e1, 4));
log(" Average protein density = %s um^-2 (%s%%)", Utils.rounded(fitProteinDensity * 1e6, 4), Utils.rounded(e2, 4));
log(" Domain radius = %s nm", Utils.rounded(domainRadius, 4));
log(" Domain density = %s", Utils.rounded(domainDensity, 4));
log(" nCluster = %s", Utils.rounded(nCluster, 4));
// Check the fitted parameters are within tolerance of the initial estimates
valid2 = true;
if (fittingTolerance > 0 && (Math.abs(e1) > fittingTolerance || Math.abs(e2) > fittingTolerance)) {
log(" Failed to fit %s within tolerance (%s%%): Average precision = %f nm (%s%%), average protein density = %g um^-2 (%s%%)", clusteredModel.getName(), Utils.rounded(fittingTolerance, 4), fitSigmaS, Utils.rounded(e1, 4), fitProteinDensity * 1e6, Utils.rounded(e2, 4));
valid2 = false;
}
// Check extra parameters. Domain radius should be higher than the precision. Density should be positive
if (domainRadius < fitSigmaS) {
log(" Failed to fit %s: Domain radius is smaller than the average precision (%s < %s)", clusteredModel.getName(), Utils.rounded(domainRadius, 4), Utils.rounded(fitSigmaS, 4));
valid2 = false;
}
if (domainDensity < 0) {
log(" Failed to fit %s: Domain density is negative (%s)", clusteredModel.getName(), Utils.rounded(domainDensity, 4));
valid2 = false;
}
if (ic2 > ic1) {
log(" Failed to fit %s - Information Criterion has increased %s%%", clusteredModel.getName(), Utils.rounded((100 * (ic2 - ic1) / ic1), 4));
valid2 = false;
}
addResult(clusteredModel.getName(), resultColour, valid2, fitSigmaS, fitProteinDensity, domainRadius, domainDensity, nCluster, 0, ic2);
return parameters;
}
use of org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.Optimum in project GDSC-SMLM by aherbert.
the class PCPALMFitting method runBoundedOptimiser.
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;
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) {
// No need to restart
break;
} catch (RuntimeException e) {
// No need to restart
break;
} 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.NEGATIVE_INFINITY;
double stopFitness = 0;
boolean isActiveCMA = true;
int diagonalOnly = 0;
int checkFeasableCount = 1;
//Well19937c();
RandomGenerator random = new Well44497b();
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;
}
Aggregations