| Minimizer.java |
package ij.measure;
import ij.*;
import ij.gui.*;
import ij.macro.*;
import ij.gui.Roi;
import ij.gui.PolygonRoi;
import ij.gui.Line;
import ij.util.Tools;
import ij.plugin.frame.RoiManager;
import java.util.Random;
import java.util.Arrays;
import java.util.Vector;
import java.util.Hashtable;
/** Minimizer based on Nelder-Mead simplex method (also known as polytope method),
* including the 'outside contraction' as described in:
* J.C. Lagarias, J.A. Reeds, M.H. Wright, P. Wright:
* Convergence properties of the Nelder-Mead simplex algorithm in low dimensions.
* SIAM J. Optim. 9, 112-147 (1998).
* Differences w.r.t. this publication:
* - If outside contraction is rejected, instead of shrinking the whole simplex, an inside
* contraction is tried first. Own experiments show that this results in slightly better
* performance for some test functions (Perm, Mc Kinnon's function with tau=2, theta=6,
* Osborne1 curve-fitting problem). In most cases, there is no difference, however.
* - This implementation does not include any 'ordering rules' in case of equal function values.
* - When checking for convergence, a special iteration step may be performed, improving
* the best vertex of the simplex.
*
* Re-initialization within a minimization run:
* In some cases, the simplex algorithm may fail to find the minimum or the convergence
* check might stop it prematurely. Therefore, the search is initialized again, keeping the
* best vertex of the simplex and setting the other vertices to random values, but keeping
* variations of the parameters w.r.t. the best vertex at a similar value. If re-initializing
* the simplex does not lead to a significant improvement, the value is accepted as true
* (local) minimum.
*
* Multiple minimization runs:
* In spite of re-initializing (see above), there are rare cases where minimization is stopped
* too early. Also, minimization may result in a local minimum. Therefore, unless determined
* otherwise by setting 'setRestarts', two minimization runs with different initialization
* of the simplex are started in parallel threads. If the results don't agree within the
* error limit, two more minimization runs are started. This is repeated until the two best
* results agree within the error limits or the number of restarts (determined by 'setRestarts';
* default 2, i.e., up to 3 runs with two threads each) is exceeded.
* This does not guarantee that the minimum is a global minimum, however: A local minimum
* will be accepted if the minimizer finds a local minimum twice (or two different local
* minima with the same function value within the error bounds), but no better minimum has
* been found at that time.
*
* The user-supplied target function should return NaN for out-of-bounds parameters instead
* of a high (penalty) value (minimization is faster and more reliable with NaNs).
* The region where the function is defined (e.g. not returning NaN) must be convex.
* Sharp corners of the region where the function value is defined (especially in higher dimensions)
* may cause a problem with finding suitable test points when (re-)initializing the simplex.
* If all attempts to find initial points result in NaN, the status returned is
* INITIALIZATION_FAILURE.
*
* Versions:
* Michael Schmid 2012-01-30: first version, based on previous CurveFitter
* 2012-11-20: mor tries to find initial params not leading to NaN
* 2013-09-24: 50% higher maximum iteration count, and never uses more than 0.4*maxIter
* iterations per minimization to avoid trying too few sets of initial params
* 2013-10-13: setStatusAndEscape to show iterations and enable abort by ESC
* 2014-09-16: normalization bug fixed. makeNewParamVariations checks for paramResolutions
*
*/
public class Minimizer {
/** Status returned: successful completion */
public final static int SUCCESS = 0;
/** Status returned: Could not initialize the simplex because either the initialParams
* resulted in the target function returning NaN or all attempts to find starting
* parameters for the other simplex points resulted in the target function returning NaN.
* No minimization was possible. */
public final static int INITIALIZATION_FAILURE = 1;
/** Status returned: aborted by call to abort method. */
public final static int ABORTED = 2;
/** Status returned: Could not reinitialize the simplex because all attempts to find restarting
* parameters resulted in the target function returning NaN. Reinitialization is
* needed to obtain a reliable result; so the result may be inaccurate or wrong. */
public final static int REINITIALIZATION_FAILURE = 3;
/** Status returned: no convergence detected after maximum iteration count */
public final static int MAX_ITERATIONS_EXCEEDED = 4;
/** Status returned: not two equal solutions after maximum number of restarts */
public final static int MAX_RESTARTS_EXCEEDED = 5;
/** Strings describing the status codes */
public final static String[] STATUS_STRING = { "Success",
"Initialization failure; no result",
"Aborted",
"Re-initialization failure (inaccurate result?)",
"Max. no. of iterations reached (inaccurate result?)",
"Max. no. of restarts reached (inaccurate result?)"};
private final static double C_REFLECTION = 1.0; // reflection coefficient
private final static double C_CONTRACTION = 0.5; // contraction coefficient
private final static double C_EXPANSION = 2.0; // expansion coefficient
private final static double C_SHRINK = 0.5; // shrink coefficient
private final static int ITER_FACTOR = 750; // maximum number of iterations per numParams^2; twice that value for 2 threads
private final static int WORST=0, NEXT_WORST=1, BEST=2;//indices in array to pass the numbers of the respective vertices
private int numParams; // number of independent variables (parameters)
private int numVertices; // numParams+1, number of vertices of the simplex
private int numExtraArrayElements; // number of extra array elements for private use in user function
private UserFunction userFunction; // target function to minimize
private double maxRelError = 1e-10; // max relative error
private double maxAbsError = 1e-100; // max absolute error
private double[] paramResolutions; // differences of parameters less than these values are considered 0
private int maxIter; // stops after this number of iterations
private int totalNumIter; // number of iterations performed in all tries so far
private int numCompletedMinimizations; // number of minimizations completed
private int maxRestarts = 2; // number of times to try finding local minima; each uses 2 threads if restarts>0
private int randomSeed; // for starting the random number generator
private boolean useSingleThread = // single thread if one processor or set by useSingleThread(true)
Runtime.getRuntime().availableProcessors() <= 1;
private int status; // SUCCESS or error code
private boolean wasInitialized; // initialization was successful at least once
private double[] result; // result data+function value
private Vector<double[]> resultsVector; // holds several results if multiple tries; the best one is kept.
private String ijStatusString = null; // shown together with iteration count in status, no display if null
private boolean checkEscape; // whether to stop when Escape is pressed
private int nextIterationForStatus = 10;// next iteration when we should display the status
private long startTime; // of the whole minimization process
/*private Hashtable<Thread, double[][]> simpTable =
new Hashtable<Thread, double[][]>(); //for each thread, holds a reference to its simplex */
//-----------------------------------------------------------------------------
// We can't set the function in the constructor because the CurveFitter
// allows specifying the fit function after other variables
/**
* Set the the target function, i.e. function whose value should be minimized.
* @param userFunction The class having a function to minimize (implementing
* the UserFunction interface).
* This function must allow simultaneous calls in multiple threads unless
* setMaximumThreads(1); has been called.
* Note that the function will be called with at least numParams+1 array
* elements; the last one is needed to store the function value. Further
* array elements for private use in the user function may be added by
* calling setExtraArrayElements.
* @param numParams Number of independent variables (also called parameters)
* of the function.
*/
public void setFunction(UserFunction userFunction, int numParams) {
if (maxIter<=0) {
maxIter = ITER_FACTOR*numParams*numParams;
if (maxRestarts > 0)
maxIter *= 2;
}
this.userFunction = userFunction;
this.numParams = numParams;
this.numVertices = numParams+1;
}
/** Perform minimization with the gradient-enhanced simplex method once or a few
* times, depending on the value of 'restarts'. Running it several times helps
* to reduce the probability of finding local minima or accepting one of the rare
* results where the minimizer has got stuck before finding the true minimum.
* We are using two threads and terminate after two equal results. Thus, apart
* from the overhead of starting a new thread (typically < 1 ms), for unproblematic
* minimization problems on a dual-core machine this is almost as fast as running
* it once.
*
* Use 'setFunction' first to define the function and number of parameters.
* Afterwards, use the 'getParams' method to access the result.
*
* @param initialParams Array with starting values of the parameters (variables).
* When null, the starting values are assumed to be 0.
* The target function must be defined (not returning NaN) for
* the values specified as initialParams.
* @param initialParamVariations Parameters (variables) are initially varied by up to +/-
* this value. If not given (null), initial variations are taken as
* 10% of initial parameter value or 0.01 for parameters that are zero.
* When this array is given, all elements must be positive (nonzero).
* If one or several initial parameters are zero, is advisable to set the initialParamVariations
* array to useful values indicating the typical order of magnitude of the parameters.
* For target functions with only one minimum, convergence is fastest with large values of
* initialParamVariations, so that the expected value is within initialParam+/-initialParamVariations.
* If local minima can occur, it is best to use a value close to the expected global minimum,
* and rather small initialParamVariations, much lower than the distance to the nearest local
* minimum.
*
* @return status code; SUCCESS if two attempts have found minima with the
* same value (within the error bounds); so a minimum has been found
* with very high probability.
*/
public int minimize(final double[] initialParams, final double[] initialParamVariations) {
status = SUCCESS;
resultsVector = new Vector<double[]>();
int maxLoopCount = maxRestarts+1;
if (useSingleThread) maxLoopCount*=2; // if we have only one thread, loop twice as many times
for (int i=0; i<maxLoopCount; i++) { // try several times, until we have twice the same result
Thread secondThread = null;
if (maxRestarts>0 && !useSingleThread) { // set up 2nd thread to minimize
final int seed = randomSeed+1000000+i;
final Thread thread = new Thread(
new Runnable() {
final public void run() {
minimizeOnce(initialParams, initialParamVariations, seed);
}
}, "Minimizer-1"
);
thread.setPriority( Thread.currentThread().getPriority() );
thread.start();
secondThread = thread;
}
minimizeOnce(initialParams, initialParamVariations, randomSeed+i); //minimize in main thread
if (secondThread != null) try {
secondThread.join(); // wait until send thread is done
} catch (InterruptedException e) {}
if (resultsVector.size() == 0 && result==null)
return status;
if (result==null)
result = (double[])resultsVector.get(0);
for (double[] r : resultsVector) // find best result so far
if (value(r) < value(result))
result = r;
if (status != SUCCESS && status != REINITIALIZATION_FAILURE && status != MAX_ITERATIONS_EXCEEDED)
return status; // no more tries if permanent error or aborted
if (totalNumIter >= maxIter)
return MAX_ITERATIONS_EXCEEDED; // no more tries if too many iterations
for (int ir=0; ir<resultsVector.size(); ir++)
if (!belowErrorLimit(value((double[])resultsVector.get(ir)), value(result), 1.0)) {
resultsVector.remove(ir); // discard results that are significantly worse
ir --;
}
if (resultsVector.size() >= 2) return SUCCESS; // if we have two (almost) equal results, it's enough
} //for i <= maxRestarts
if (ijStatusString != null)
IJ.showStatus(""); // reset status display
return maxRestarts>0 ?
MAX_RESTARTS_EXCEEDED : // number of restarts exceeded without two equal results
status; // if only one run was required, we can't have 2 equal results
}
/** Perform minimization with the simplex method once, including re-initialization until
* we have a stable solution.
* Use the 'getParams' method to access the result.
*
* @param initialParams Array with starting values of the parameters (variables).
* When null, the starting values are assumed to be 0.
* The target function must be defined (not returning NaN) for
* the values specified as initialParams.
* @param initialParamVariations Parameters (variables) are initially varied by up to +/-
* this value. If not given (null), iniital variations are taken as
* 10% of inial parameter value or 0.01 for parameters that are zero.
* When this array is given, all elements must be positive (nonzero).
* If one or several initial parameters are zero, is advisable to set the initialParamVariations
* array to useful values indicating the typical order of magnitude of the parameters.
* For target functions with only one minimum, convergence is fastest with large values of
* initialParamVariations, so that the expected value is within initialParam+/-initialParamVariations.
* If local minima can occur, it is best to use a value close to the expected global minimum,
* and rather small initialParamVariations, much lower than the distance to the nearest local
* minimum.
*
* @return status code; SUCCESS if it is considered likely that a minimum of the
* target function has been found.
*/
public int minimizeOnce(double[] initialParams, double[] initialParamVariations) {
status = SUCCESS;
minimizeOnce(initialParams, initialParamVariations, randomSeed);
return status;
}
/** Get the result, i.e., the set of parameter values (i.e., variable values)
* from the best corner of the simplex. Note that the array returned may have more
* elements than numParams; ignore the rest.
* May return an array with only NaN values in case the minimize call has returned
* an INITIALIZATION_FAILURE status or that abort() has been called at the very
* beginning of the minimization.
* Do not call this method before minimization. */
public double[] getParams() {
if (result == null) {
result = new double[numParams+1+numExtraArrayElements];
Arrays.fill(result, Double.NaN);
}
return result;
}
/** Get the value of the minimum, i.e. the value associated with the resulting parameters
* as obtained by getParams(). May return NaN in case the minimize call has returned
* an INITIALIZATION_FAILURE status or that abort() has been called at the very
* beginning of the minimization.
* Do not call this method before minimization. */
public double getFunctionValue() {
if (result == null) {
result = new double[numParams+1];
Arrays.fill(result, Double.NaN);
}
return value(result);
}
/** Get number of iterations performed (includes all restarts). One iteration needs
* between one and numParams+3 calls of the target function (typically two calls
* per iteration) */
public int getIterations() {
return totalNumIter;
}
/** Set maximum number of iterations allowed (including all restarts and all threads).
* The number of function calls will be higher, up to about twice the number of
* iterations.
* Note that the number of iterations needed typically scales with the square of
* the dimensions (i.e., numParams^2).
* Default value is 1000 * numParams^2 (half that value if maxRestarts=0), which is
* enough for >99% of all cases (if the maximum number of restarts is set to 2);
* typical number of iterations are below 10 and 20% of this value.
*/
public void setMaxIterations(int x) {
maxIter = x;
}
/** Get maximum number of iterations allowed. Unless given by 'setMaxIterations',
* this value is defined only after running 'setFunction' */
public int getMaxIterations() {
return maxIter;
}
/** Set maximum number of minimization restarts to do.
* With n=0, the minimizer is run once in a single thread.
* With n>0, two threads are used, and if the two results do not agree within
* the error bounds, additional optimizations are done up to n times, each
* with two threads. In any case, if the two best results are within the error
* bounds, the best result is accepted.
* Thus, on dual-core machines running no other jobs, values of n=1 or n=2 (default)
* do not cause a notable increase of computing time for 'easy' optimization problems,
* while greatly reducing the risk of running into spurious local minima or non-
* optimal results due to the minimizer getting stuck. In problematic cases, the
* improved
* The 'n' value does not refer to the restarts within one minimization run
* (there, at least one restart is performed, and restart is repeated until the result
* does not change within the error bounds).
* This value does not affect the 'minimizeOnce' function call.
* When setting the maximum number of restarts to a value much higher than 3, remember
* to adjust the maximum number of iterations (see setMaxIterations).
*/
public void setMaxRestarts(int n) {
maxRestarts = n;
}
/** Get maximum number of minimization restarts to do */
public int getMaxRestarts() {
return maxRestarts;
}
/** Get number of minimizations completed (i.e. not aborted or stopped because the
* number of minimization was exceeded). After a minimize(..) call, typically 2
* for unproblematic cases. Higher number indicate a functin that is difficult to
* minimize or the existence of more than one minimum.
*/
public int getCompletedMinimizations() {
return numCompletedMinimizations;
}
/** Set a seed to initialize the random number generator, which is used for initialization
* of the simplex. */
public void setRandomSeed(int n) {
randomSeed = n;
}
/** Sets the accuracy of convergence. Minimizing is done as long as the
* relative error of the function value is more than this number (Default: 1e-10).
*/
public void setMaxError(double maxRelError) {
this.maxRelError = maxRelError;
}
/** Sets the accuracy of convergence. Minimizing is done as long as the
* relative error of the function value is more than maxRelError (Default: 1e-10)
* and the maximum absolute error is more than maxAbsError
* (i.e. it is enough to fulfil one of these two criteria)
*/
public void setMaxError(double maxRelError, double maxAbsError) {
this.maxRelError = maxRelError;
this.maxAbsError = maxAbsError;
}
/** Set the resolution of the parameters, for problems where the target function is not smooth
* but suffers from numerical noise. If all parameters of all vertices are closer to the
* best value than the respective resolution value, minimization is finished, irrespective
* of the difference of the target function values at the vertices */
public void setParamResolutions(double[] paramResolutions) {
this.paramResolutions = paramResolutions;
}
/** Call setMaximuThreads(1) to avoid multi-threaded execution (in case the user-provided
* target function does not allow moultithreading). Currently a maximum of 2 thread is used
* irrespective of any higher value. */
public void setMaximumThreads (int numThreads) {
useSingleThread = numThreads <= 1;
}
/** Aborts minimization. Calls to getParams() will return the best solution found so far.
* This method may be called from the user-supplied target function.
* If displayStatusAndCheckEsc has been called before, the Minimizer itself checks for the ESC key.
*/
public void abort() {
status = ABORTED;
}
/** Create output on the number of iterations in the ImageJ Status line, e.g.
* "<ijStatusString> 50 (max 750); ESC to stop"
* @param ijStatusString Displayed in the beginning of the status message. No display if null.
* E.g. "Optimization: Iteration "
* @param checkEscape When true, the Minimizer stops if escape is pressed and the status
* becomes ABORTED. Note that checking for ESC does not work in the Event Queue thread. */
public void setStatusAndEsc(String ijStatusString, boolean checkEscape) {
this.ijStatusString = ijStatusString;
this.checkEscape = checkEscape;
}
/** Add a given number of extra elements to array of parameters (independent vaiables)
* for private use in the user function. Note that the first numParams+1 elements
* should not be touched.*/
public void setExtraArrayElements(int numExtraArrayElements) {
this.numExtraArrayElements = numExtraArrayElements;
}
/** Get the full simplex of the current thread. This may be useful if the target function
* wants to modify the simplex */
/* public double[][] getSimplex() {
return simpTable.get(Thread.currentThread());
} */
/** One minimization run (including reinitializations of the simplex until the result is stable) */
private void minimizeOnce(double[] initialParams, double[] initialParamVariations, int seed) {
Random random = new Random(seed);
double[][] simp = makeSimplex(initialParams, initialParamVariations, random);
if (simp == null) {
status = wasInitialized ? REINITIALIZATION_FAILURE : INITIALIZATION_FAILURE;
return;
}
wasInitialized = true;
if (startTime == 0)
startTime = System.currentTimeMillis();
//if (IJ.debugMode) showSimplex(simp, seed+" Initialized:");
int bestVertexNumber = minimize(simp); // first minimization
double bestValueSoFar = value(simp[bestVertexNumber]);
//reinitialize until converged or error/aborted (don't care about reinitialization failure in other thread)
boolean reinitialisationFailure = false;
while (status == SUCCESS || status == REINITIALIZATION_FAILURE) {
double[] paramVariations =
makeNewParamVariations(simp, bestVertexNumber, initialParams, initialParamVariations);
if (!reInitializeSimplex(simp, bestVertexNumber, paramVariations, random)) {
reinitialisationFailure = true;
break;
}
//if (IJ.debugMode) showSimplex(simp, seed+" Reinitialized:");
bestVertexNumber = minimize(simp); // minimize with reinitialized simplex
if (belowErrorLimit(value(simp[bestVertexNumber]), bestValueSoFar, 2.0)) break;
bestValueSoFar = value(simp[bestVertexNumber]);
}
if (reinitialisationFailure)
status = REINITIALIZATION_FAILURE;
else if (status == SUCCESS || status == REINITIALIZATION_FAILURE) //i.e. not aborted, not max iterations exceeded
numCompletedMinimizations++; // it was a complete minimization
//if (IJ.debugMode) showSimplex(simp, seed+" Final:");
if (resultsVector != null) synchronized(resultsVector) {
resultsVector.add(simp[bestVertexNumber]);
} else
result = simp[bestVertexNumber];
}
/** Minimizes the target function by variation of the simplex.
* Note that one call to this function never does more than 0.4*maxIter iterations.
* @return index of the best value in simp
*/
private int minimize(double[][] simp) {
int[] worstNextBestArray = new int[3]; // used to pass these indices from 'order' function
double[] center = new double[numParams+1+numExtraArrayElements]; // center of all vertices except worst
double[] reflected = new double[numParams+1+numExtraArrayElements]; // the 1st new vertex to try
double[] secondTry = new double[numParams+1+numExtraArrayElements]; // expanded or contracted vertex
order(simp, worstNextBestArray);
int worst = worstNextBestArray[WORST];
int nextWorst = worstNextBestArray[NEXT_WORST];
int best = worstNextBestArray[BEST];
//showSimplex(simp, "before minimization, value="+value(simp[best]));
//String operation="ini";
int thisNumIter=0;
while (true) {
totalNumIter++; // global count over all threads
thisNumIter++; // local count for this minimize call
// THE MINIMIZAION ALGORITHM IS HERE
iteration: {
getCenter(simp, worst, center); // centroid of vertices except worst
// Reflect worst vertex through centriod of not-worst
getVertexAndEvaluate(center, simp[worst], -C_REFLECTION, reflected);
if (value(reflected) <= value(simp[best])) { // if it's better than the best...
// try expanding it
getVertexAndEvaluate(center, simp[worst], -C_EXPANSION, secondTry);
if (value(secondTry) <= value(reflected)) {
copyVertex(secondTry, simp[worst]); // if expanded is better than reflected, keep it
//operation="expa";
break iteration;
}
}
if (value(reflected) < value(simp[nextWorst])) {
copyVertex(reflected, simp[worst]); // keep reflected if better than 2nd worst
//operation="refl";
break iteration;
} else if (value(reflected) < value(simp[worst])) {
// try outer contraction
getVertexAndEvaluate(center, simp[worst], -C_CONTRACTION, secondTry);
if (value(secondTry) <= value(reflected)) {
copyVertex(secondTry, simp[worst]); // keep outer contraction
//operation="outC";
break iteration;
}
} else if (value(reflected) > value(simp[worst]) || Double.isNaN(value(reflected))) {
// else inner contraction
getVertexAndEvaluate(center, simp[worst], C_CONTRACTION, secondTry);
if (value(secondTry) < value(simp[worst])) {
copyVertex(secondTry, simp[worst]); // keep contracted if better than 2nd worst
//operation="innC";
break iteration;
}
}
// if everything else has failed, contract simplex in on best
shrinkSimplexAndEvaluate(simp, best);
//operation="shri";
break iteration;
} // iteration:
boolean checkParamResolution = // if new 'worst' is not close to 'best', don't check any further
paramResolutions!=null && belowResolutionLimit(simp[worst], simp[best]);
order(simp, worstNextBestArray);
worst = worstNextBestArray[WORST];
nextWorst = worstNextBestArray[NEXT_WORST];
best = worstNextBestArray[BEST];
if (checkParamResolution)
if (belowResolutionLimit(simp, best)) // check whether all parameters are within the resolution limit
break;
if (belowErrorLimit(value(simp[best]), value(simp[worst]), 4.0)) { // make sure we are at the minimum:
getCenter(simp, -1, secondTry); // check center of the simplex
evaluate(secondTry);
if (value(secondTry) < value(simp[best]))
copyVertex(secondTry, simp[best]); // better than best: keep
}
if (belowErrorLimit(value(simp[best]), value(simp[worst]), 4.0)) // no large spread of values
break; // looks like we are at the minimum
if (totalNumIter > maxIter || thisNumIter>4*(maxIter/10))
status = MAX_ITERATIONS_EXCEEDED;
if (status != SUCCESS)
break;
if ((ijStatusString != null || checkEscape) && totalNumIter > nextIterationForStatus) {
long time = System.currentTimeMillis();
nextIterationForStatus = totalNumIter + (int)(totalNumIter*500L/(time-startTime+1)); //next display 0.5 sec later
if (time - startTime > 1000L) { // display status and check for ESC after the first second
if (checkEscape && IJ.escapePressed()) {
status = ABORTED;
IJ.resetEscape();
IJ.showStatus(ijStatusString+" ABORTED");
break;
}
if (ijStatusString != null) {
String statusString = ijStatusString+totalNumIter+" ("+maxIter+" max)";
if (checkEscape) statusString += " ESC to stop";
IJ.showStatus(statusString);
}
}
}
}
//showSimplex(simp, "after "+totalNumIter+" iterations: value="+value(simp[best]));
return best;
}
/** Move along the line between center and worst, where +1 corresponds to worstVertex
* The new point is written to'newVertex' and target function is evaluated there
*/
private void getVertexAndEvaluate(double[] center, double[] worstVertex, double howFar, double[] newVertex) {
for (int i = 0; i < numParams; i++)
newVertex[i] = (1.-howFar)*center[i] + howFar*worstVertex[i];
evaluate(newVertex);
}
/** Get center of all vertices except one to exclude
* (may be -1 to exclude none)
* Does not care about function values (i.e., last array elements) */
private void getCenter(double[][]simp, int excludeVertex, double[] center) {
Arrays.fill(center, 0.);
int nV = 0;
for (int v=0; v<numVertices; v++)
if (v != excludeVertex) {
for (int i=0; i<numParams; i++)
center[i] += simp[v][i];
nV++;
}
double norm = 1.0/nV;
for (int i=0; i<numParams; i++)
center[i] *= norm;
}
private void shrinkSimplexAndEvaluate(double[][] simp, int best) {
for (int v=0; v<numVertices; v++)
if (v != best) {
for (int i=0; i<numParams; i++)
simp[v][i] = C_SHRINK*simp[v][i] + (1-C_SHRINK)*simp[best][i];
evaluate(simp[v]);
}
}
/** Whether the highest and lowest values fulfill the error limit criterion.
* Error limits are reduced by a factor of 'sensitivity' */
private boolean belowErrorLimit(double highest, double lowest, double sensitivity) {
double absError = sensitivity*Math.abs(highest - lowest);
double relError = absError/(Math.max(Math.abs(highest),Math.abs(lowest))+1e-100);
return relError < maxRelError || absError < maxAbsError;
}
/** Initialise the simplex and evaluate it. Returns null on failure. */
private double[][] makeSimplex(double[] initialParams, double[] initialParamVariations, Random random) {
double[][] simp = new double[numVertices][numParams+1+numExtraArrayElements];
/* simpTable.put(Thread.currentThread(), simp); */
if (initialParams!=null) {
for (int i=0; i<numParams; i++)
if (Double.isNaN(initialParams[i]))
if (IJ.debugMode) IJ.log("Warning: Initial Parameter["+i+"] is NaN");
System.arraycopy(initialParams, 0, simp[0], 0, Math.min(initialParams.length, numParams));
}
evaluate(simp[0]);
if (Double.isNaN(value(simp[0]))) {
if (IJ.debugMode) showVertex(simp[0], "Warning: Initial Parameters yield NaN:");
findValidInitalParams(simp[0], initialParamVariations, random);
}
if (Double.isNaN(value(simp[0]))) {
if (IJ.debugMode) IJ.log("Error: Could not find initial parameters not yielding NaN:");
return null;
}
if (initializeSimplex(simp, initialParamVariations, random))
return simp;
else {
if (IJ.debugMode) showSimplex(simp, "Error: Could not make simplex vertices not yielding NaN");
return null;
}
}
/** Whether the distance between the best and any other simplex vertex
* is less than the resolution of the parameters */
private boolean belowResolutionLimit(double[][] simp, int best) {
for (int v=0; v<numVertices; v++)
if (v!= best && !belowResolutionLimit(simp[v], simp[best]))
return false;
return true;
}
/** Whether the distance between two vertices is less than the resolution of the parameters */
private boolean belowResolutionLimit(double[] vertex1, double[] vertex2) {
for (int i=0; i<numParams; i++)
if (Math.abs(vertex1[i]-vertex2[i]) >= paramResolutions[i])
return false;
return true;
}
/** Find initial parameters not yielding a result of NaN in case those given yield NaN
* Called with params containing the initial parameters that have been tried previously */
private void findValidInitalParams(double[] params, double[] initialParamVariations, Random random) {
final int maxAttempts = 50*numParams*numParams; //max number of attempts to find params that do not lead to NaN
double rangeFactor = 1; // will gradually become larger to handle different orders of magnitude
final double rangeMultiplyLog = Math.log(1e20)/(maxAttempts-1); //will try up to 1e-20 to 1e20*initialParamVariations
double[] firstParams = new double[numParams]; // remember starting params (which may be modified)
double[] variations = new double[numParams]; // new values of parameter variations
for (int i=0; i<numParams; i++) {
firstParams[i] = Double.isNaN(params[i]) ? 0 : params[i];
variations[i] = initialParamVariations!=null ? initialParamVariations[i] : 0.1*firstParams[i];
if (Double.isNaN(variations[i]) || Math.abs(variations[i])<1e-10 || Math.abs(variations[i])>1e10)
variations[i] = 0.1;
}
for (int attempt=0; attempt<maxAttempts; attempt++) {
for (int i=0; i<numParams; i++) {
double multiplier = attempt<maxAttempts/10 ? 1 : //after a while try different orders of magnitude
Math.exp(rangeMultiplyLog*attempt*2*(random.nextDouble()-0.5));
params[i] = multiplier*(firstParams[i]+2*(random.nextDouble()-0.5)*variations[i]);
}
evaluate(params);
//showVertex(params,"findValidInitalParams attempt "+attempt);
if (!Double.isNaN(value(params))) return; //found a valid parameter set
}
}
/** Reinitialize an existing simplex: Create new vertices around the best one, keeping the rough size of
* the simplex. This helps to avoid premature termination or cases where the simplex has become degenerate.
* All points of the new simplex are evaluated already.
* Returns false on error (if it could not create a simplex not resulting in NaN). */
private boolean reInitializeSimplex(double[][] simp, int bestVertexNumber, double[] paramVariations, Random random) {
if (bestVertexNumber != 0) {
double[] swap = simp[0];
simp[0] = simp[bestVertexNumber];
simp[bestVertexNumber] = swap;
}
return initializeSimplex(simp, paramVariations, random);
}
/** Initialize simplex vertices except the number 0, which will be taken as starting point.
* All points of the new simplex are evaluated already.
* Returns false on error (if it could not create a simplex not resulting in NaN).
*/
private boolean initializeSimplex(double[][] simp, double[] paramVariations, Random random) {
double[] variations = new double[numParams]; // improved values of parameter variations
for (int i=0; i<numParams; i++) {
double range = paramVariations!=null && i<paramVariations.length ?
paramVariations[i] : 0.1 * Math.abs(simp[0][i]); // if nothing else given, based on initial parameter
if (range < 1e-100 || Double.isNaN(range))
range = 0.01; //range must be nonzero
if (Math.abs(range/simp[0][i])<1e-10)
range = Math.abs(simp[0][i]*1e-10); //range must be more than very last digits
variations[i] = range;
//IJ.log("param#"+i+" variation="+(float)variations[i]);
}
final int maxAttempts = 100*numParams; //max number of attempts to find params that do not lead to NaN
for (int v=1; v<numVertices; v++) {
int numTries = 0;
do { //try finding a vertex that does not yield NaN
if (numTries++ > maxAttempts)
return false;
// Create a vector orthogonal to all others in a coordinate system where every value is
// normalized to its respective variations value. We first use this coordinate system where
// every parameter is in the [-1, +1] range.
// Don't orthogonalize after getting only NaN function values in many attempts; it might be
// easier to find viable parameters without normalization.
for (int i=0; i<numParams; i++)
simp[v][i] = 2*random.nextFloat() - 1.0;
if (numTries < maxAttempts/3) // (don't orthogonalize if finding valid params was very difficult
//IJ.log("orthogonalize vertex "+v);
for (int v2=1; v2<v; v2++) { // to avoid a degenerate simplex, make orthogonal to all others
double lengthSqr = 0;
double innerProduct = 0;
for (int i=0; i<numParams; i++) {
double x2 = (simp[v2][i] - simp[0][i])/variations[i]; // convert to [-1, +1] range
lengthSqr += x2*x2;
innerProduct += x2*simp[v][i];
}
for (int i=0; i<numParams; i++) //subtract component parallel to the other vector
simp[v][i] -= ((simp[v2][i] - simp[0][i])/variations[i])*(innerProduct/lengthSqr);
}
double sumSqr = 0; // normalize the 'excursion vector' to length = 1
for (int i=0; i<numParams; i++)
sumSqr += simp[v][i]*simp[v][i];
if (sumSqr < 1e-14) continue; // bad luck with random numbers, excursion vector almost zero
if (numTries > 2 && sumSqr > 1e-6)
sumSqr = 1; // don't normalize if attempts with normalized were unsuccessful
for (int i=0; i<numParams; i++)
simp[v][i] = simp[0][i] + simp[v][i]/Math.sqrt(sumSqr)*variations[i];
evaluate(simp[v]);
} while (Double.isNaN(value(simp[v])) && (status == SUCCESS || status == REINITIALIZATION_FAILURE));
}
//showSimplex(simp, 0);
return true;
}
/** Calculate ranges for parameter variations for re-initializing:
* Note that the simplex should have adapted in size to the useful range of each
* parameter, but in special cases it might have become degenerate [i.e. the
* (hyper-)volume spanned by the vertices might be zero; in such a case it may
* happen that one parameter has the same value for all vertices of the simplex].
*
* For each parameter, as a variation value we use the typical difference w.r.t. the best
* vertex of the simplex.
* To avoid a degenerate simplex, we have to know what the the relative orders of magnitude of
* the parameter variations should be. If initialParamVariations are given, we use them as
* indication of the typoical relative values; otherwise we use the parameter values themselves
* and assume that the variation of each parameter should be roughly the same fraction of the value.
* If the variation of a parameter in the simplex w.r.t. this estimated variation is much
* smaller (<10^3) than typical for the parameters, we increase the variation for this parameter.
*/
private double[] makeNewParamVariations(double[][] simp, int bestVertexNumber,
double[] initialParams, double[] initialParamVariations) {
double[] paramVariations = new double[numParams]; // will be the return value
double[] relatedTo = new double[numParams]; // parameter variation should be related to (in size) these values
double logTypicalRelativeVariation = 0;
for (int i=0; i<numParams; i++) { // for all parameters
double variation = 0; // roughly size of simplex in direction of that parameter
for (int v=0; v<numVertices; v++)
if (v != bestVertexNumber) {
double delta = Math.abs(simp[v][i] - simp[bestVertexNumber][i]);
paramVariations[i] += delta*delta;
}
paramVariations[i] = 10*Math.sqrt(paramVariations[i]); // make the new simplex larger than the old one
relatedTo[i] = initialParamVariations!=null && initialParamVariations.length>= numParams ?
initialParamVariations[i] : // relate to initialParamVariations (if given),
Math.max(Math.abs(initialParams[i]), Math.abs(simp[bestVertexNumber][i]));
double logRelativeVariation = paramVariations[i]>relatedTo[i] ? 0 :
Math.log(paramVariations[i]/relatedTo[i]);
logTypicalRelativeVariation += logRelativeVariation;
}
logTypicalRelativeVariation /= numParams;
double typicalRelativeVariation = Math.exp(logTypicalRelativeVariation);
final double WORST_RATIO = 1e-3; //parameter variation should not be lower than typical value by more than this
for (int i=0; i<numParams; i++) {
if (paramVariations[i]<relatedTo[i] && paramVariations[i]/relatedTo[i] < typicalRelativeVariation*WORST_RATIO)
paramVariations[i] = relatedTo[i]*typicalRelativeVariation*WORST_RATIO;
if (paramResolutions != null && paramVariations[i] < 10*paramResolutions[i])
paramVariations[i] = 10*paramResolutions[i];// parameter variation must be well above numeric noise limit
}
return paramVariations;
}
/** Evaluate the target function for the parameters given by 'vertex'
* The function value is stored into the last (extra) array element.
*/
private void evaluate(double[] vertex) {
vertex[numParams] = userFunction.userFunction(vertex, 0);
}
/** Function value of a vertex after call to 'evaluate'
* (stored it as last array element) */
private double value(double[] vertex) {
return vertex[numParams];
}
/** Keep the newVertex: replace a given Vertex with it */
void copyVertex(double[] newVertex, double[] vertex) {
System.arraycopy(newVertex, 0, vertex, 0, newVertex.length);
}
/** Find the simplex vertices with the worst, nextWorst and best values
* of the target function */
private void order(double[][] simp, int[] worstNextBestArray) {
int worst = 0, nextWorst = 0, best = 0;
for (int i = 0; i < numVertices; i++) {
if (value(simp[i]) < value(simp[best])) best = i;
if (value(simp[i]) > value(simp[worst])) worst = i;
}
nextWorst = best;
for (int i = 0; i < numVertices; i++)
if (i != worst && value(simp[i]) > value(simp[nextWorst]))
nextWorst = i;
worstNextBestArray[WORST] = worst;
worstNextBestArray[NEXT_WORST] = nextWorst;
worstNextBestArray[BEST] = best;
}
// Display simplex [Iteration: s0(p1, p2....), s1(),....] in Log window
private synchronized void showSimplex(double[][] simp, String heading) {
IJ.log("Minimizer: "+heading);
for (int i = 0; i < numVertices; i++)
showVertex(simp[i], null);
}
private synchronized void showVertex(double[] vertex, String heading) {
if (heading != null)
IJ.log(heading);
String s = "";
for (int j=0; j < numParams; j++)
s += " " + IJ.d2s(vertex[j], 8,12);
s += " -> " + IJ.d2s(value(vertex), 8,12);
IJ.log(s);
}
}