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 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 simpTable =
new Hashtable(); //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();
int maxLoopCount = maxRestarts+1;
if (useSingleThread) maxLoopCount*=2; // if we have only one thread, loop twice as many times
for (int i=0; i0 && !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= 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.
* " 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= 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; i1e10)
variations[i] = 0.1;
}
for (int attempt=0; attempt 2 && sumSqr > 1e-6)
sumSqr = 1; // don't normalize if attempts with normalized were unsuccessful
for (int i=0; i= 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 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);
}
}
Ads help maintain this website.