```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
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
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
*                      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
if (maxRestarts>0 && !useSingleThread) {  // set up 2nd thread to minimize
final int seed = randomSeed+1000000+i;
new Runnable() {
final public void run() {
minimizeOnce(initialParams, initialParamVariations, seed);
}
}, "Minimizer-1"
);
}
minimizeOnce(initialParams, initialParamVariations, randomSeed+i); //minimize in main thread
if (secondThread != null) try {
} 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() {
}

/** 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;
}

*  target function does not allow moultithreading). Currently a maximum of 2 thread is used
*  irrespective of any higher value. */
}

/** 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() {
} */

/** 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) {
} 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;          // 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];

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, Math.min(initialParams.length, numParams));
}
evaluate(simp);
if (Double.isNaN(value(simp))) {
if (IJ.debugMode) showVertex(simp, "Warning: Initial Parameters yield NaN:");
findValidInitalParams(simp, initialParamVariations, random);
}
if (Double.isNaN(value(simp))) {
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;
simp = 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[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[i])<1e-10)
range = Math.abs(simp[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[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[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[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) {