Build in fit list CurveFitter

[Wim Wakelkamp]

Dear ImageJ,

I am working with some bio engineers on a testing program. They're really
interested in the fitted data based on the 3-parametric Chapman-Richards
equation F=a(1-e-bt)c , that's how i get to ImageJ. Right now i am using
the Curvefitter.java class and trying to extends the build in fit list with
this new formula (Chapman-Richards). When i am running doFit() and
acquiring parameters with getparams() based on this fitType, i do not seem
to receive to correct a,b,c values. Attached the edited Curvefitter class.
In this class i edited everything under the CHAPMAN variable.

Hope you guys can put me in the right direction.

--
Kind regards

Wim

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Wim Wakelkamp]

Now with the file actually attached. I also noticed that the formula is not correctly formatted in the previuous message, but it's also within the CurveFitter.java file.

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

** Curve fitting class based on the Simplex method in the Minimizer class
 *
 *
 * Notes on fitting polynomial functions:
 * (i) The range of x values should not be too far from 0, especially for higher-order polynomials.
 * For polynomials of 4th order, the average x value should fulfill |xMean| < 2*(xMax-xMin).
 * For polynomials of 5th order or higher, the x range should encompass x=0; and for
 * 7th and 8th order it is desirable to have x=0 near the center of the x range.
 * (ii) In contrast to some fitting algorithms using the normal equations and matrix inversion, the
 * simplex algorithm used here can cope with parameters having very different orders of magnitude,
 * as long as the coefficients of the polynomial are within a reasonable range (say, 1e-80 to 1e80).
 * Thus, it is usually not needed to scale the x values, even for high-order polynomials.
 *
 * Version history:
 *
 * 2008-01-21: Modified to do Gaussian fitting by Stefan Woerz (s.woerz at dkfz.de).
 * 2012-01-30: Modified for external Minimizer class and UserFunction fit by Michael Schmid.
 *  - Also fixed incorrect equation string for 'Gamma Variate' & 'Rodbard (NIH Image)',
 *  - Added 'Inverse Rodbard', 'Exponential (linear regression)', 'Power (linear regression)'
 * functions and polynomials of order 5-8. Added 'nicely' sorted list of types.
 *  - Added absolute error for minimizer to avoid endless minimization if exact fit is possible.
 *  - Added 'initialParamVariations' (order of magnitude if parameter variations) -
 * this is important for safer and better convergence.
 *  - Linear regression for up to 2 linear parameters, reduces the number of parameters handled
 * by the simplex Minimizer and improves convergence. These parameters can be an offset and
 * either a linear slope or a factor that the full function is multiplied with.
 * 2012-10-07: added GAUSSIAN_NOOFFSET fit type
 * 2012-11-20: Bugfix: exception on Gaussian&Rodbard with initial params, bad initial params for Gaussian
 */

public class CurveFitter implements UserFunction{
        /** Constants for the built-in fit functions */
        public static final int STRAIGHT_LINE=0, POLY2=1, POLY3=2, POLY4=3,
        EXPONENTIAL=4, POWER=5, LOG=6, RODBARD=7, GAMMA_VARIATE=8, LOG2=9,
        RODBARD2=10, EXP_WITH_OFFSET=11, GAUSSIAN=12, EXP_RECOVERY=13,
        INV_RODBARD=14, EXP_REGRESSION=15, POWER_REGRESSION=16,
        POLY5=17, POLY6=18, POLY7=19, POLY8=20,
        GAUSSIAN_NOOFFSET=21, CHAPMAN=22;

        /** Nicer sequence of the built-in function types */
        public static final int[] sortedTypes = { STRAIGHT_LINE, POLY2, POLY3, POLY4,
                        POLY5, POLY6, POLY7, POLY8,
                        POWER, POWER_REGRESSION,
                        EXPONENTIAL, EXP_REGRESSION, EXP_WITH_OFFSET, EXP_RECOVERY,
                        LOG, LOG2,
                        RODBARD, RODBARD2, INV_RODBARD,
                        GAUSSIAN, GAUSSIAN_NOOFFSET,
                        GAMMA_VARIATE, CHAPMAN
        };

        /** Names of the built-in fit functions */
        public static final String[] fitList = {"Straight Line","2nd Degree Polynomial",
        "3rd Degree Polynomial", "4th Degree Polynomial","Exponential","Power",
        "Log","Rodbard", "Gamma Variate", "y = a+b*ln(x-c)","Rodbard (NIH Image)",
        "Exponential with Offset","Gaussian", "Exponential Recovery",
        "Inverse Rodbard", "Exponential (linear regression)", "Power (linear regression)",
        "5th Degree Polynomial","6th Degree Polynomial","7th Degree Polynomial","8th Degree Polynomial",
        "Gaussian (no offset), Chapman"
        }; // fList, doFit(), getNumParams() and makeInitialParamsAndVariations() must also be updated

        /** Equations of the built-in fit functions */
        public static final String[] fList = {
        "y = a+bx","y = a+bx+cx^2", //STRAIGHT_LINE,POLY2
        "y = a+bx+cx^2+dx^3","y = a+bx+cx^2+dx^3+ex^4",
        "y = a*exp(bx)","y = a*x^b", "y = a*ln(bx)", //EXPONENTIAL,POWER,LOG
        "y = d+(a-d)/(1+(x/c)^b)", "y = b*(x-a)^c*exp(-(x-a)/d)", //RODBARD,GAMMA_VARIATE
        "y = a+b*ln(x-c)", "x = d+(a-d)/(1+(y/c)^b)", //LOG2,RODBARD2
        "y = a*exp(-bx) + c", "y = a + (b-a)*exp(-(x-c)*(x-c)/(2*d*d))", //EXP_WITH_OFFSET,GAUSSIAN
        "y = a*(1-exp(-b*x)) + c", "y = c*((x-a)/(d-x))^(1/b)", //EXP_RECOVERY, INV_RODBARD
        "y = a*exp(bx)", "y = a*x^b", //EXP_REGRESSION, POWER_REGRESSION
        "y = a+bx+cx^2+dx^3+ex^4+fx^5", "y = a+bx+cx^2+dx^3+ex^4+fx^5+gx^6",
        "y = a+bx+cx^2+dx^3+ex^4+fx^5+gx^6+hx^7", "y = a+bx+cx^2+dx^3+ex^4+fx^5+gx^6+hx^7+ix^8",
        "y = a*exp(-(x-b)*(x-b)/(2*c*c)))", //GAUSSIAN_NOOFFSET
        "y = a*(1-exp(-b*x))^c" //chapman
        };

        /** @deprecated, now in the Minimizer class (since ImageJ 1.46f).
         * (probably of not much value for anyone anyhow?) */
        public static final int IterFactor = 500;

        private static final int CUSTOM = 100; // user function defined in macro or plugin
        private static final int GAUSSIAN_INTERNAL = 101; // Gaussian with separate offset & multiplier
        private static final int RODBARD_INTERNAL = 102; // Rodbard with separate offset & multiplier

        private int fitType = -1; // Number of curve type to fit
        private double[] xData, yData; // x,y data to fit
        private double[] xDataSave, yDataSave; //saved original data after fitting modified data
        private int numPoints; // number of data points in actual fit
        private double ySign = 0; // remember sign of y data for power-law fit via regression
        private double sumY = Double.NaN, sumY2 = Double.NaN; // sum(y), sum(y^2) of the data used for fitting
        private int numParams; // number of parameters
        private double[] initialParams; // user specified or estimated initial parameters
        private double[] initialParamVariations; // estimate of range of parameters
        private double[] minimizerInitialParams;  // modified initialParams of the minimizer
        private double[] minimizerInitialParamVariations; // modified initialParamVariations of the minimizer
        private double maxRelError = 1e-10;// maximum relative deviation of sum of residuals^2 from minimum
        private long time; // elapsed time in ms
        private int customParamCount; // number of parameters of user-supplied function (macro or plugin)
        private String customFormula; // equation defined in macro or text
        private UserFunction userFunction; // plugin with custom fit function
        //private Interpreter macro; // holds macro with custom equation
        private int macroStartProgramCounter; // equation in macro starts here
        private int numRegressionParams;// number of parameters that can be calculated by linear regression
        private int offsetParam = -1; // parameter number: function has this parameter as offset
        private int factorParam = -1; // index of parameter that is slope or multiplicative factor of the function
        private boolean hasSlopeParam; // whether the 'factorParam' is the slope of the function
        private double[] finalParams; // parameters obtained by fit
        private boolean linearRegressionUsed; // whether linear regression alone was used instead of the minimizer
        private boolean restrictPower; // power via linear regression fit: (0,0) requires positive power
        private Minimizer minimizer = new Minimizer();
        private int minimizerStatus = Minimizer.INITIALIZATION_FAILURE; // status of the minimizer after minimizing
        private String errorString; // in case of error before invoking the minimizer
        private static String[] sortedFitList; // names like fitList, but in more logical sequence
        private static Hashtable<String, Integer> namesTable; // converts fitList String into number


        /** Construct a new CurveFitter. */
        public CurveFitter (double[] xData, double[] yData) {
                this.xData = xData;
                this.yData = yData;
                numPoints = xData.length;
        }
       
        /** Perform curve fitting with one of the built-in functions
         * doFit(fitType) does the fit quietly
         * Use getStatus() and/or getStatusString() to see whether fitting was (probably) successful and
         * getParams() to access the result.
         */
        public void doFit(int fitType) {
                doFit(fitType, false);
        }
       
        /** Perform curve fitting with one of the built-in functions
         * doFit(fitType, true) pops up a dialog allowing the user to set the initial
         * fit parameters and various numbers controlling the Minimizer
         * Use getStatus() and/or getStatusString() to see whether fitting was (probably) successful and
         * getParams() to access the result.
         */
        public void doFit(int fitType, boolean showSettings) {
                Log.e("test", "fit type = " + fitType + " && fitlist length = " + fitList.length);
                if (!(fitType>=STRAIGHT_LINE && fitType<=fitList.length || fitType==CUSTOM))
                        throw new IllegalArgumentException("Invalid fit type");
// if (fitType==CUSTOM && macro==null && userFunction==null)
// throw new IllegalArgumentException("No custom formula!");
                this.fitType = fitType;
                if (isModifiedFitType(fitType)) // these fits don't use the original data and a different fit type (this.fitType)
                        if (!prepareModifiedFitType(fitType)) return;
                numParams = getNumParams();
                if (fitType != CUSTOM)
                        getOffsetAndFactorParams();
                //IJ.log("special params: off="+offsetParam+(hasSlopeParam ? " slo=" : " fac=")+factorParam+" numPar="+numParams+" numRegressPar="+numRegressionParams);
                calculateSumYandY2(); // sumY, sumY2 needed for regression, abs Error; R, goodness of modified fit functions
                long startTime = System.currentTimeMillis();
                if (this.fitType == STRAIGHT_LINE) { // no minimizer needed
                        finalParams = new double[] {0, 0, 0}; // (does not work with 0 initial slope)
                        doRegression(finalParams); // fit by regression; save sum(Residuals^2) as last array element
                        linearRegressionUsed = true;
                } else { // use simplex minimizer
                        minimizer.setFunction(this, numParams-numRegressionParams);
                        minimizer.setExtraArrayElements(numRegressionParams);  // reserve space for extra parameters if minimizer has fewer paramters
        // if (macro != null) minimizer.setMaximumThreads(1); // macro interpreter does not allow multithreading
                        if (!makeInitialParamsAndVariations(fitType)) // also includes some data checking
                                return; // initialization failure
                        //if (showSettings) settingsDialog();
                        if (numRegressionParams >0)
                                modifyInitialParamsAndVariations();
                        else {
                                minimizerInitialParams = initialParams;
                                minimizerInitialParamVariations = initialParamVariations;
                        }
                        startTime = System.currentTimeMillis();
                        // The maximum absolute error of the fit must be specified in case the
                        // fit function fits perfectly, i.e. the sume of residuals approaches 0.
                        // In such a case, the maximum relative error is meaningless and the
                        // minimizer would run until it reaches the maximum iteration count.
                        double maxAbsError = Math.min(1e-6,maxRelError)*Math.sqrt(sumY2);
                        minimizer.setMaxError(maxRelError, maxAbsError);
                        //{String s="initVariations:";for(int ii=0;ii<numParams;ii++)s+=" ["+ii+"]:"+IJ.d2s(initialParamVariations[ii],5,9);IJ.log(s);}
                        //{String s="minInitVariations:";for(int ii=0;ii<numParams;ii++)s+=" ["+ii+"]:"+IJ.d2s(minimizerInitialParamVariations[ii],5,9);IJ.log(s);}
                        //{String s="minInitPars:";for(int ii=0;ii<numParams;ii++)s+=" ["+ii+"]:"+IJ.d2s(minimizerInitialParams[ii],5,9);IJ.log(s);}
                        // m i n i m i z a t i o n of squared residuals
                        minimizerStatus = minimizer.minimize(minimizerInitialParams, minimizerInitialParamVariations);
                        finalParams = minimizer.getParams();
                        if (numRegressionParams > 0)
                                minimizerParamsToFullParams(finalParams, false);
                }
                if (isModifiedFitType(fitType))
                        postProcessModifiedFitType(fitType);
                time = System.currentTimeMillis()-startTime;

        }

        /** Fit a function defined as a macro String like "y = a + b*x + c*x*x".
         * Returns the number of parameters, or 0 in case of a macro syntax error.
         *
         * For good performance, it is advisable to set also the typical variation range
         * of the initial parameters by the
         * getMinimizer().setInitialParamVariations(double[]) method (especially if one or
         * more of the initialParams are zero).
         * Use getStatus() and/or getStatusString() to see whether fitting was (probably) successful and
         * getParams() to access the result.
         */
// public int doCustomFit(String equation, double[] initialParams, boolean showSettings) {
// customFormula = null;
// customParamCount = 0;
// Program pgm = (new Tokenizer()).tokenize(equation);
// if (!pgm.hasWord("y")) return 0;
// if (!pgm.hasWord("x")) return 0;
// String[] params = {"a","b","c","d","e","f"};
// for (int i=0; i<params.length; i++) {
// if (pgm.hasWord(params[i]))
// customParamCount++;
// }
// if (customParamCount==0)
// return 0;
// customFormula = equation;
// String code =
// "var x, a, b, c, d, e, f;\n"+
// "function dummy() {}\n"+
// equation+";\n"; // starts at program counter location 21
// macroStartProgramCounter = 21;
// macro = new Interpreter();
// try {
// macro.run(code, null);
// } catch (Exception e) {
// if (!Macro.MACRO_CANCELED.equals(e.getMessage()))
// IJ.handleException(e);
// }
// if (macro.wasError())
// return 0;
// this.initialParams = initialParams;
// doFit(CUSTOM, showSettings);
// return customParamCount;
// }

        /** Fit a function defined in a user plugin implementing the UserFunction interface
         *
         * Use getStatus() and/or getStatusString() to see whether fitting was (probably) successful and
         * getParams() to access the result.
         *
         * @param userFunction A plugin where the fit function is defined by the
         * userFunction(params, x) method.
         * This function must allow simultaneous calls in multiple threads.
         * @param numParams Number of parameters of the fit function.
         * @params formula A String describing the fit formula, may be null.
         * @param initialParams Starting point for the parameters; may be null (than values
         * of 0 are used). The fit function with these parameters must
         * not return NaN for any of the data points given in the
         * constructor (xData).
         * @param initialParamVariations Each parameter is 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).
         * See Minimizer.minimize for details.
         * @param showSettings Displays a popup dialog for modifying the initial parameters and
         * a few numbers controlling the minimizer.
         */
        public void doCustomFit(UserFunction userFunction, int numParams, String formula,
                        double[] initialParams, double[] initialParamVariations, boolean showSettings) {
                this.userFunction = userFunction;
                this.customParamCount = numParams;
                this.initialParams = initialParams;
                this.initialParamVariations = initialParamVariations;
                customFormula = formula==null ? "(defined in plugin)" : formula;
                doFit(CUSTOM, showSettings);
        }

        /** Sets the initial parameters, which override the default initial parameters. */
        public void setInitialParameters(double[] initialParams) {
                this.initialParams = initialParams;
        }

        /** Returns a reference to the Minimizer used, for accessing Minimizer methods directly.
         * Note that no Minimizer is used if fitType is any of STRAIGHT_LINE, EXP_REGRESSION,
         * and POWER_REGRESSION. */
        public Minimizer getMinimizer() {
                return minimizer;
        }

        /** For improved fitting performance when using a custom fit formula, one may
         * specify parameters that can be calculated directly by linear regression.
         * For values not used, set the index to -1
         *
         * @param offsetParam  Index of a parameter that is a pure offset:
         *   E.g. '0' if f(p0, p1, p2...) = p0 + function(p1, p2, ...).
         * @param multiplyParam Index of a parameter that is purely multiplicative.
         *   E.g. multiplyParams=1 if f(p0, p1, p2, p3...) can be expressed as
         *   p1*func(p0, p2, p3, ...) or p0 +p1*func(p0, p2, p3, ...) with '0' being
         *   the offsetparam.
         * @param slopeParam   Index of a parameter that is multiplied with x and then summed to the function.
         *   E.g. '1' for f(p0, p1, p2, p3...) = p1*x + func(p0, p2, p3, ...)
         *   Only one, multiplyParam and slopeParam can be used (ie.e, the other
         *   should be set to -1)
         */
        public void setOffsetMultiplySlopeParams(int offsetParam, int multiplyParam, int slopeParam) {
                this.offsetParam = offsetParam;
                hasSlopeParam = slopeParam >= 0;
                factorParam = hasSlopeParam ? slopeParam : multiplyParam;
                numRegressionParams = 0;
                if (offsetParam >=0) numRegressionParams++;
                if (factorParam >=0) numRegressionParams++;
        }

        /** Get number of parameters for current fit formula
         * Do not use before 'doFit', because the fit function would be undefined. */
        public int getNumParams() {
                switch (fitType) {
                        case STRAIGHT_LINE: return 2;
                        case POLY2: return 3;
                        case POLY3: return 4;
                        case POLY4: return 5;
                        case POLY5: return 6;
                        case POLY6: return 7;
                        case POLY7: return 8;
                        case POLY8: return 9;
                        case EXPONENTIAL: case EXP_REGRESSION: return 2;
                        case POWER: case POWER_REGRESSION:   return 2;
                        case LOG: return 2;
                        case RODBARD: case RODBARD2: case INV_RODBARD: case RODBARD_INTERNAL: return 4;
                        case GAMMA_VARIATE: return 4;
                        case LOG2: return 3;
                        case EXP_WITH_OFFSET: return 3;
                        case GAUSSIAN: case GAUSSIAN_INTERNAL: return 4;
                        case GAUSSIAN_NOOFFSET: return 3;
                        case EXP_RECOVERY: return 3;
                        case CHAPMAN: return 3;
                        case CUSTOM: return customParamCount;
                }
                return 0;
        }

        /** Returns the formula value for parameters 'p' at 'x'.
         * Do not use before 'doFit', because the fit function would be undefined. */
        public final double f(double[] p, double x) {
                //if (fitType!=CUSTOM)
                        return f(fitType, p, x);
// else {
// if (macro == null) { // function defined in plugin
// return userFunction.userFunction(p, x);
// } else { // function defined in macro
// macro.setVariable("x", x);
// macro.setVariable("a", p[0]);
// if (customParamCount>1) macro.setVariable("b", p[1]);
// if (customParamCount>2) macro.setVariable("c", p[2]);
// if (customParamCount>3) macro.setVariable("d", p[3]);
// if (customParamCount>4) macro.setVariable("e", p[4]);
// if (customParamCount>5) macro.setVariable("f", p[5]);
// macro.run(macroStartProgramCounter);
// return macro.getVariable("y");
// }
// }
        }

        /** Returns value of built-in 'fitType' formula value for parameters "p" at "x" */
        public static double f(int fitType, double[] p, double x) {
                switch (fitType) {
                        case STRAIGHT_LINE:
                                return p[0] + x*p[1];
                        case POLY2:
                                return p[0] + x*(p[1] + x*p[2]);
                        case POLY3:
                                return p[0] + x*(p[1] + x*(p[2] + x*p[3]));
                        case POLY4:
                                return p[0] + x*(p[1] + x*(p[2] + x*(p[3] + x*p[4])));
                        case POLY5:
                                return p[0] + x*(p[1] + x*(p[2] + x*(p[3] + x*(p[4] + x*p[5]))));
                        case POLY6:
                                return p[0] + x*(p[1] + x*(p[2] + x*(p[3] + x*(p[4] + x*(p[5] + x*p[6])))));
                        case POLY7:
                                return p[0] + x*(p[1] + x*(p[2] + x*(p[3] + x*(p[4] + x*(p[5] + x*(p[6] + x*p[7]))))));
                        case POLY8:
                                return p[0] + x*(p[1] + x*(p[2] + x*(p[3] + x*(p[4] + x*(p[5] + x*(p[6] + x*(p[7]+x*p[8])))))));
                        case EXPONENTIAL:
                        case EXP_REGRESSION:
                                return p[0]*Math.exp(p[1]*x);
                        case EXP_WITH_OFFSET:
                                return p[0]*Math.exp(-p[1]*x)+p[2];
                        case EXP_RECOVERY:
                                return p[0]*(1-Math.exp(-p[1]*x))+p[2];
                        case CHAPMAN:
                               
                                double value =  p[0]* (Math.pow((1-(Math.exp(-p[1]*x))), p[2]));
                                Log.e("test", "values = " + value);
                                return value;
                        case GAUSSIAN:
                                return p[0]+(p[1]-p[0])*Math.exp(-(x-p[2])*(x-p[2])/(2.0*p[3]*p[3]));
                        case GAUSSIAN_INTERNAL: // replaces GAUSSIAN for the fitting process
                                return p[0]+p[1]*Math.exp(-(x-p[2])*(x-p[2])/(2.0*p[3]*p[3]));
                        case GAUSSIAN_NOOFFSET:
                                return p[0]*Math.exp(-(x-p[1])*(x-p[1])/(2.0*p[2]*p[2]));
                        case POWER: // ax^b
                        case POWER_REGRESSION:
                                return p[0]*Math.pow(x,p[1]);
                        case LOG:
                                if (x == 0.0)
                                        return -1000*p[0];
                                return p[0]*Math.log(p[1]*x);
                        case RODBARD: { // d+(a-d)/(1+(x/c)^b)
                                double ex = Math.pow(x/p[2], p[1]); //(x/c)^b
                                return p[3]+(p[0]-p[3])/(1.0+ex); }
                        case RODBARD_INTERNAL: { // d+a/(1+(x/c)^b) , replaces RODBARD of the fitting process
                                double ex = Math.pow(x/p[2], p[1]); //(x/c)^b
                                return p[3]+p[0]/(1.0+ex); }
                        case GAMMA_VARIATE: // b*(x-a)^c*exp(-(x-a)/d)
                                if (p[0] >= x) return 0.0;
                                if (p[1] <= 0) return Double.NaN;
                                if (p[2] <= 0) return Double.NaN;
                                if (p[3] <= 0) return Double.NaN;
                               
                                double pw = Math.pow((x - p[0]), p[2]);
                                double e = Math.exp((-(x - p[0]))/p[3]);
                                return p[1]*pw*e;
                        case LOG2:
                                double tmp = x-p[2];
                                if (tmp<=0)
                                        return Double.NaN;
                                return p[0]+p[1]*Math.log(tmp);
                        case INV_RODBARD: // c*((x-a)/(d-x))^(1/b), the inverse Rodbard function
                        case RODBARD2: // also used after the 'Rodbard NIH Image' fit
                                double y;
                                if (p[3]-x < 2*Double.MIN_VALUE || x<p[0]) // avoid x>=d (singularity) and x<a (negative exponent)
                                        y = fitType==INV_RODBARD ? Double.NaN : 0.0;
                                else {
                                        y = (x-p[0])/(p[3]-x); //(x-a)/(d-x) = ( (a-d)/(x-d) -1 )
                                        y = Math.pow(y,1.0/p[1]);  //y=y**(1/b)
                                        y = y*p[2];
                                }
                                return y;
                        default:
                                return 0.0;
                }
        }
       
        /** Get the result of fitting, i.e. the set of parameter values for the best fit.
         * Note that the array returned may have more elements than numParams; ignore the rest.
         * May return an array with only NaN values if the minimizer could not start properly,
         * i.e., if getStatus() returns Minimizer.INITILIZATION_FAILURE.
         * See Minimizer.getParams() for details.
         */
        public double[] getParams() {
                return finalParams==null ? minimizer.getParams() : finalParams; //if we have no result, take all_NaN result from the Minimizer
        }
       
        /** Returns residuals array, i.e., differences between data and curve.
         * The residuals are with respect to the real data, also for fit types where the data are
         * modified before fitting (power&exp fit by linear regression, 'Rodbard NIH Image' ).
         * This is in contrast to sum of squared residuals, which is for the fit that was actually done.
         */
        public double[] getResiduals() {
                double[] params = getParams();
                double[] residuals = new double[xData.length];
                for (int i=0; i<xData.length; i++)
                        residuals[i] = yData[i] - f(params, xData[i]);
                return residuals;
        }

        /* Get the sum of the residuals (may be NaN if the minimizer could not start properly
         * i.e., if getStatus() returns Minimizer.INITILIZATION_FAILURE).
         */
        public double getSumResidualsSqr() {
                return getParams()[numParams]; // value is stored as last element by the minimizer
        }

        /** Returns the standard deviation of the residuals.
         * Here, the standard deviation is defined here as the root-mean-square of the residuals
         * times sqrt(n/(n-1)); where n is the number of points.
         */
        public double getSD() {
                double[] residuals = getResiduals();
                int n = residuals.length;
                double sum=0.0, sum2=0.0;
                for (int i=0; i<n; i++) {
                        sum += residuals[i];
                        sum2 += residuals[i]*residuals[i];
                }
                double stdDev = (sum2-sum*sum/n); //sum of squared residuals
                return Math.sqrt(stdDev/(n-1.0));
        }
       
        /** Returns R^2, where 1.0 is best.
        <pre>
         r^2 = 1 - SSE/SSD
         
         where: SSE = sum of the squared errors
                                 SSD = sum of the squared deviations about the mean.
        </pre>
         * For power, exp by linear regression and 'Rodbard NIH Image', this is calculated for the
         * fit actually done, not for the residuals of the original data.
        */
        public double getRSquared() {
                if (Double.isNaN(sumY)) calculateSumYandY2();
                double sumMeanDiffSqr = sumY2 - sumY*sumY/numPoints;
                double rSquared = 0.0;
                if (sumMeanDiffSqr > 0.0)
                        rSquared = 1.0 - getSumResidualsSqr()/sumMeanDiffSqr;
                return rSquared;
        }

        /** Get a measure of "goodness of fit" where 1.0 is best.
         * Approaches R^2 if the number of points is much larger than the number of fit parameters.
         * For power, exp by linear regression and 'Rodbard NIH Image', this is calculated for the
         * fit actually done, not for the residuals of the original data.
         */
        public double getFitGoodness() {
                if (Double.isNaN(sumY)) calculateSumYandY2();
                double sumMeanDiffSqr = sumY2 - sumY*sumY/numPoints;
                double fitGoodness = 0.0;
                int degreesOfFreedom = numPoints - getNumParams();
                if (sumMeanDiffSqr > 0.0 && degreesOfFreedom > 0)
                        fitGoodness = 1.0 - (getSumResidualsSqr()/ sumMeanDiffSqr) * numPoints / (double)degreesOfFreedom;

                return fitGoodness;
        }

        public int getStatus() {
                return linearRegressionUsed ? Minimizer.SUCCESS : minimizerStatus;
        }

        /** Get a short text with a short description of the status. Should be preferred over
         * Minimizer.STATUS_STRING[getMinimizer().getStatus()] because getStatusString()
         * better explains the problem in some cases of initialization failure (data not
         * compatible with the fit function chosen) */
        public String getStatusString() {
                return errorString != null ? errorString : minimizer.STATUS_STRING[getStatus()];
        }
       
        /** Get a string with detailed description of the curve fitting results (several lines,
         * including the fit parameters).
         */
// public String getResultString() {
// String resultS =  "\nFormula: " + getFormula() +
// "\nStatus: "+getStatusString();
// if (!linearRegressionUsed) resultS += "\nNumber of completed minimizations: " + minimizer.getCompletedMinimizations();
// resultS += "\nNumber of iterations: " + getIterations();
// if (!linearRegressionUsed) resultS += " (max: " + minimizer.getMaxIterations() + ")";
// resultS += "\nTime: "+time+" ms" +
// "\nSum of residuals squared: " + IJ.d2s(getSumResidualsSqr(),5,9) +
// "\nStandard deviation: " + IJ.d2s(getSD(),5,9) +
// "\nR^2: " + IJ.d2s(getRSquared(),5) +
// "\nParameters:";
// char pChar = 'a';
// double[] pVal = getParams();
// for (int i = 0; i < numParams; i++) {
// resultS += "\n " + pChar + " = " + IJ.d2s(pVal[i],5,9);
// pChar++;
// }
// return resultS;
// }

        /** Set maximum number of simplex restarts to do. See Minimizer.setMaxRestarts for details. */
        public void setRestarts(int maxRestarts) {
                minimizer.setMaxRestarts(maxRestarts);
        }

        /** Set the maximum error. by which the sum of residuals may deviate from the true value
         * (relative w.r.t. full sum of rediduals). Possible range: 0.1 ... 10^-16 */
        public void setMaxError(double maxRelError) {
                if (Double.isNaN(maxRelError)) return;
                if (maxRelError > 0.1) maxRelError = 0.1;
                if (maxRelError < 1e-16) maxRelError = 1e-16; // can't be less than numerical accuracy
                this.maxRelError = maxRelError;
        }

        /** Get number of iterations performed. Returns 1 in case the fit was done by linear regression only. */
        public int getIterations() {
                return linearRegressionUsed ? 1 : minimizer.getIterations();
        }
       
        /** Get maximum number of iterations allowed (sum of iteration count for all restarts) */
        public int getMaxIterations() {
                return minimizer.getMaxIterations();
        }
       
        /** Set maximum number of iterations allowed (sum of iteration count for all restarts) */
        public void setMaxIterations(int maxIter) {
                minimizer.setMaxIterations(maxIter);
        }
       
        /** Get maximum number of simplex restarts to do. See Minimizer.setMaxRestarts for details. */
        public int getRestarts() {
                return minimizer.getMaxRestarts();
        }

        /** Returns the status of the Minimizer after doFit.  Minimizer.SUCCESS indicates a
         * successful completion. In case of Minimizer.INITIALIZATION_FAILURE, fitting could
         * not be performed because the data and/or initial parameters are not compatible
         * with the function value.  In that case, getStatusString may explain the problem.
         * For further status codes indicating problems during fitting, see the status codes
         * of the Minimzer class. */

        /** returns the array with the x data */
        public double[] getXPoints() {
                return xData;
        }

        /** returns the array with the y data */
        public double[] getYPoints() {
                return yData;
        }

        /** returns the code of the fit type of the fit performed */
        public int getFit() {
                return fitType;
        }

        /** returns the name of the fit function of the fit performed */
        public String getName() {
                if (fitType==CUSTOM)
                        return "User-defined";
                if (fitType==GAUSSIAN_INTERNAL)
                        fitType = GAUSSIAN;
                else if (fitType==RODBARD_INTERNAL)
                        fitType = RODBARD;
                return fitList[fitType];
        }

        /** returns a String with the formula of the fit function used */
        public String getFormula() {
                if (fitType==CUSTOM)
                        return customFormula;
                else
                        return fList[fitType];
        }

        /** Returns an array of fit names with nicer sorting */
        public static String[] getSortedFitList() {
                if (sortedFitList == null) {
                        String[] l = new String[fitList.length];
                        for (int i=0; i<fitList.length; i++)
                                sortedFitList[i] = fitList[sortedTypes[i]];
                        sortedFitList = l;
                }
                return sortedFitList;
        }

        /** Returns the code for a fit with given name as defined in fitList, or -1 if not found */
        public static int getFitCode(String fitName) {
                if (namesTable == null) {
                        Hashtable<String,Integer> h = new Hashtable<String,Integer>();
                        for (int i=0; i<fitList.length; i++)
                                h.put(fitList[i], new Integer(i));
                        namesTable = h;
                }
                Integer i = (Integer)namesTable.get(fitName);
                return i!=null? i.intValue() : -1;
        }

        /** This function is called by the Minimizer and calculates the sum of squared
         * residuals for given parameters.
         * To improve the efficiency, simple linear dependencies are solved directly
         * by linear regression; in that case the corresponding parameters are modified.
         * This effectively reduces the number of free parameters by one or two and thereby
         * significantly improves the performance of minimization.
         */
        public final double userFunction(double[] params, double dummy) {
                double sumResidualsSqr = 0.0;
                if (numRegressionParams == 0) { // simply calculate sum of residuals
                        for (int i=0; i<numPoints; i++) {
                                double fValue = f(params,xData[i]);
                                sumResidualsSqr += sqr(fValue-yData[i]);
                        }
                } else { // handle simple linear dependencies by linear regression:
                        //if(getIterations()<1){String s="minimizerPar:";for(int ii=0;ii<=numParams;ii++)s+=" ["+ii+"]:"+IJ.d2s(params[ii],5,9);IJ.log(s);}
                        minimizerParamsToFullParams(params, true);
                        doRegression(params);
                        sumResidualsSqr = fullParamsToMinimizerParams(params);
                }
                return sumResidualsSqr;
        }

        /** For fits where linear regression is used to reduce the number of parameters handled
         * by the Minimizer, convert Minimizer parameters to the complete set of parameters.
         * When not for calculating regression, we use the sum of squared residuals,
         * offset and multiplication factor stored in the extra array elements:
         * The Minimizer stores the sum of squared residuals directly behind its last parameter.
         * The next element is the value of the offsetParam (if any).
         * The final element is the value of the factorParam (if any).
         */
        private void minimizerParamsToFullParams(double[] params, boolean forRegression) {
                boolean shouldTransformToSmallerParams = false;
                double offset = 0;
                double factor = hasSlopeParam ? 0 : 1; //for regression, calculate function value without x*slope, but multiplied with 1
                double sumResidualsSqr = 0;
                if (!forRegression) { // recover regression-calculated parameters from extra array elements
                        int i = params.length - 1;
                        if (factorParam >= 0)
                                factor = params[i--];
                        if (offsetParam >= 0)
                                offset = params[i];
                        sumResidualsSqr = params[numParams-numRegressionParams]; // sum of squared residuals has been calculated already
                        params[numParams] = sumResidualsSqr; // ... and is now stored in its new (proper) place
                }
                // move array elements to position in array with full number of parameters
                for (int i=numParams-1, iM=numParams-numRegressionParams-1; i>=0; i--) {
                        if (i == offsetParam)
                                params[i] = offset;
                        else if (i == factorParam)
                                params[i] = factor;
                        else
                                params[i] = params[iM--];
                }
                params[numParams] = sumResidualsSqr;
        }

        /** Determine sum of squared residuals with linear regression.
         * The sum of squared residuals is written to the array element with index 'numParams',
         * the offset and factor params (if any) are written to their proper positions in the
         * params array */
        private void doRegression(double[] params) {
                double sumX=0, sumX2=0, sumXY=0; //sums for regression; here 'x' are function values
                double sumY=0, sumY2=0; //only calculated for 'slope', otherwise we use the values calculated already
                for (int i=0; i<numPoints; i++) {
                        double fValue = fitType == STRAIGHT_LINE ? 0 : f(params, xData[i]); // function value
                        if (Double.isNaN(fValue)) { //check for NaN now; later we need NaN checking for division-by-zero check.
                                params[numParams] = Double.NaN;
                                return; //sum of squared residuals is NaN if any value is NaN
                        }
                        //if(getIterations()==0)IJ.log(xData[i]+"\t"+yData[i]+"\t"+fValue); //x,y,function
                        if (hasSlopeParam) { // fit y = offset + slope*x + function(of other params)
                                double x = xData[i];
                                double y = yData[i] - fValue;
                                sumX += x;
                                sumX2 += x*x;
                                sumXY += x*y;
                                sumY2 += y*y;
                                sumY += y;
                        } else { // fit y = offset + factor * function(of other params)
                                double x = fValue;
                                double y = yData[i];
                                sumX += fValue;
                                sumX2 += fValue*fValue;
                                sumXY += fValue*yData[i];
                        }
                }
                if (!hasSlopeParam) {
                        sumY = this.sumY;
                        sumY2 = this.sumY2;
                }
                double factor = 0; // factor or slope
                double sumResidualsSqr = 0;
                if (offsetParam<0) { // only 'factor' parameter, no offset:
                        factor = sumXY/sumX2;
                        if (Double.isNaN(factor) || Double.isInfinite(factor))
                                factor = 0; // all 'x' values are 0, any factor (slope) will fit
                        sumResidualsSqr = sumY2 + factor*(factor*sumX2 - 2*sumXY);
                        if (sumResidualsSqr < 2e-15*sumY2)
                                sumResidualsSqr = 2e-15*sumY2;
                } else { // full linear regression or offset only. Slope is named 'factor' here
                        if (factorParam >= 0) {
                                factor = (sumXY-sumX*sumY/numPoints)/(sumX2-sumX*sumX/numPoints);
                                if (restrictPower & factor<=0) // power-law fit with (0,0) point: power must be >0
                                        factor = 1e-100;
                                else if (Double.isNaN(factor) || Double.isInfinite(factor))
                                        factor = 0; // all 'x' values are equal, any factor (slope) will fit
                        }
                        double offset = (sumY-factor*sumX)/numPoints;
                        params[offsetParam] = offset;
                        sumResidualsSqr = sqr(factor)*sumX2 + numPoints*sqr(offset) + sumY2 +
                                        2*factor*offset*sumX - 2*factor*sumXY - 2*offset*sumY;
                        // check for accuracy problem: large difference of small numbers?
                        // Don't report unrealistic or even negative values, otherwise minimization could lead
                        // into parameters where we have a numeric problem
                        if (sumResidualsSqr < 2e-15*(sqr(factor)*sumX2 + numPoints*sqr(offset) + sumY2))
                                sumResidualsSqr = 2e-15*(sqr(factor)*sumX2 + numPoints*sqr(offset) + sumY2);
                        //if(){IJ.log("sumX="+sumX+" sumX2="+sumX2+" sumXY="+sumXY+" factor="+factor+" offset=="+offset);}
                }
                params[numParams] = sumResidualsSqr;
                if (factorParam >= 0)
                        params[factorParam] = factor;
        }
        /** Convert full set of parameters to minimizer parameters and returns the sum of residuals squared.
         * The last array elements, not used by the minimizer, are the offset and factor parameters (if any)
         */
        private double fullParamsToMinimizerParams(double[] params) {
                double offset = offsetParam >=0 ? params[offsetParam] : 0;
                double factor = factorParam >=0 ? params[factorParam] : 0;
                double sumResidualsSqr = params[numParams];

                for (int i=0, iNew=0; i<numParams; i++) // leave only the parameters for the minimizer in the beginning of the array
                        if (i != factorParam && i != offsetParam)
                                params[iNew++] = params[i];
                int i = params.length - 1;
                if (factorParam >= 0)
                        params[i--] = factor;
                if (offsetParam >= 0)
                        params[i--] = offset;
                params[i--] = sumResidualsSqr;
                return sumResidualsSqr;
        }

        /** In case one or two parameters are calculated by regression and not by the minimizer:
         * Make modified initialParams and initialParamVariations for the Minimizer
         */
        private void modifyInitialParamsAndVariations() {
                minimizerInitialParams = initialParams.clone();
                minimizerInitialParamVariations = initialParamVariations.clone();
                if (numRegressionParams >  0) // convert to shorter arrays with only the parameters used by the minimizer
                        for (int i=0, iNew=0; i<numParams; i++)
                                if (i != factorParam && i != offsetParam) {
                                        minimizerInitialParams[iNew] = minimizerInitialParams[i];
                                        minimizerInitialParamVariations[iNew] = minimizerInitialParamVariations[i];
                                        iNew++;
                                }
        }

        /** Estimate values for initial parameters and their typical range for the built-in
         * function types. For fits with modifications for regression, 'fitType' is still
         * the type of the original (unmodified) fit.
         * Also checks for x values being non-negative for fit types that require this,
         * and returns false if the data cannot be fitted because of negative x.
         * In such a case, 'errorString' contains a message about the problem. */
        private boolean makeInitialParamsAndVariations(int fitType) {
                boolean hasInitialParams = initialParams != null;
                boolean hasInitialParamVariations = initialParamVariations != null;
                if (!hasInitialParams)
                        initialParams = new double[numParams];
                if (!hasInitialParamVariations)
                        initialParamVariations = new double[numParams];
                if (fitType == CUSTOM)
                        return true; //no way to guess initial parameters & variations for an unknown formula

                // Calculate some things that might be useful for predicting parameters
                double firstx = xData[0];
                double firsty = yData[0];
                double lastx = xData[numPoints-1];
                double lasty = yData[numPoints-1];
                double xMin=firstx, xMax=firstx;
                double yMin=firsty, yMax=firsty;
                double xMean=0, yMean=0;
                double xOfMax = xData[0];
                for (int i=1; i<numPoints; i++) {
                        xMean += xData[i];
                        yMean += yData[i];
                        if (xData[i]>xMax) xMax = xData[i];
                        if (xData[i]<xMin) xMin = xData[i];
                        if (yData[i]>yMax) { yMax = yData[i]; xOfMax = xData[i]; }
                        if (yData[i]<yMin) yMin = yData[i];
                }
                xMean /= numPoints;
                yMean /= numPoints;

                double slope = (lasty - firsty)/(lastx - firstx);
                if (Double.isNaN(slope) || Double.isInfinite(slope)) slope = 0;
                double yIntercept = yMean - slope * xMean;

                //We cannot fit the following cases because we would always get NaN function values
                if (xMin < 0 && fitType==POWER) {
                        errorString = "Cannot fit "+fitList[fitType]+" when x<0";
                        return false;
                } else if (xMin < 0 && xMax > 0 && fitType==RODBARD) {
                        errorString = "Cannot fit "+fitList[fitType]+" to mixture of x<0 and x>0";
                        return false;
                } else if (xMin <= 0 && fitType==LOG) {
                        errorString = "Cannot fit "+fitList[fitType]+" when x<=0";
                        return false;
                }

                if (!hasInitialParams) {
                        switch (fitType) {
                                //case POLY2: case POLY3: case POLY4: case POLY5: case POLY6: case POLY7: case POLY8:
                                // offset&slope calculated via regression; leave the others at 0

                                // also for the other cases, some initial parameters are unused; only to show them with 'showSettings'
                                case EXPONENTIAL: // a*exp(bx)   assuming change by factor of e between xMin & xMax
                                        initialParams[1] = 1./(xMax-xMin+1e-100) * Math.signum(yMean) * Math.signum(slope); //actually ignored, done by regression
                                        initialParams[0] = yMean/Math.exp(initialParams[1]*xMean);
                                        break;
                                case EXP_WITH_OFFSET: // a*exp(-bx) + c assuming b>0, change by factor of e between xMin & xMax
                                case EXP_RECOVERY: // a*(1-exp(-bx)) + c
                                        initialParams[1] = 1./(xMax-xMin+1e-100);
                                        initialParams[0] = (yMax-yMin)/Math.exp(initialParams[1]*xMean) * Math.signum(slope) *
                                                        fitType==EXP_RECOVERY ? 1 : -1; // don't care, we will do this via regression
                                        initialParams[2] = 0.5*yMean; // don't care, we will do this via regression
                                        break;
                                case POWER: // ax^b, assume linear for the beginning
                                        initialParams[0] = yMean/(Math.abs(xMean+1e-100)); // don't care, we will do this via regression
                                        initialParams[1] = 1.0;
                                        break;
                                case LOG: // a*ln(bx), assume b=e/xMax
                                        initialParams[0] = yMean; // don't care, we will do this via regression
                                        initialParams[1] = Math.E/(xMax+1e-100);
                                        break;
                                case LOG2: // y = a+b*ln(x-c)
                                        initialParams[0] = yMean; // don't care, we will do this via regression
                                        initialParams[1] = (yMax-yMin)/(xMax-xMin+1e-100); // don't care, we will do this via regression
                                        initialParams[2] = Math.min(0., -xMin-0.1*(xMax-xMin)-1e-100);
                                        break;
                                case RODBARD: // d+(a-d)/(1+(x/c)^b)
                                        initialParams[0] = firsty; // don't care, we will do this via regression
                                        initialParams[1] = 1.0;
                                        initialParams[2] = xMin < 0 ? xMin : xMax; //better than xMean;
                                        initialParams[3] = lasty; // don't care, we will do this via regression
                                        break;
                                case INV_RODBARD: case RODBARD2: // c*((x-a)/(d-x))^(1/b)
                                        initialParams[0] = xMin - 0.1 * (xMax-xMin);
                                        initialParams[1] = 1.0;
                                        initialParams[2] = yMax; // don't care, we will do this via regression
                                        initialParams[3] = xMax + (xMax - xMin);
                                        break;
                                case GAMMA_VARIATE: // // b*(x-a)^c*exp(-(x-a)/d)
                                // First guesses based on following observations (naming the 'x' coordinate 't'):
                                // t0 [a] = time of first rise in gamma curve - so use the user specified first x value
                                // tmax = t0 + c*d where tmX is the time of the peak of the curve
                                // therefore an estimate for c and d is sqrt(tm-t0)
                                // K [a] can now be calculated from these estimates
                                        initialParams[0] = xMin;
                                        double cd = xOfMax - xMin;
                                        if (cd < 0.1*(xMax-xMin)) cd = 0.1*(xMax-xMin);
                                        initialParams[2] = Math.sqrt(cd);
                                        initialParams[3] = Math.sqrt(cd);
                                        initialParams[1] = yMax / (Math.pow(cd, initialParams[2]) * Math.exp(-cd/initialParams[3])); // don't care, we will do this via regression
                                        break;
                                case GAUSSIAN: // a + (b-a)*exp(-(x-c)^2/(2d^2))
                                        initialParams[0] = yMin; //actually don't care, we will do this via regression
                                        initialParams[1] = yMax; //actually don't care, we will do this via regression
                                        initialParams[2] = xOfMax;
                                        initialParams[3] = 0.39894 * (xMax-xMin) * (yMean-yMin)/(yMax-yMin+1e-100);
                                        break;
                                case GAUSSIAN_NOOFFSET: // a*exp(-(x-b)^2/(2c^2))
                                        initialParams[0] = yMax; //actually don't care, we will do this via regression
                                        initialParams[1] = xOfMax;  //actually don't care, we will do this via regression
                                        initialParams[2] = 0.39894 * (xMax-xMin) * yMean/(yMax+1e-100);
                                        break;  
                                //case CUSTOM: // initial params should be specified anyhow, otherwise use minimizer default
                        }
                }
                if (!hasInitialParamVariations) { // estimate initial range for parameters
                        for (int i=0; i<numParams; i++)
                                initialParamVariations[i] = 0.1 * initialParams[i]; //default, should be overridden if it can be zero
                        switch (fitType) {
                                case POLY2: case POLY3: case POLY4: case POLY5: case POLY6: case POLY7: case POLY8:
                                        double xFactor = 0.5*Math.max(Math.abs(xMax+xMin), (xMax-xMin));
                                        initialParamVariations[numParams-1] = (yMax-yMin)/(Math.pow(0.5*(xMax-xMin), numParams-1)+1e-100);
                                        for (int i=numParams-2; i>=0; i--)
                                                initialParamVariations[i] = initialParamVariations[i+1]*xFactor;
                                        break;
                                case EXPONENTIAL: // a*exp(bx)  a (and c) is calculated by regression
                                case EXP_WITH_OFFSET: // a*exp(-bx) + c
                                case EXP_RECOVERY: // a*(1-exp(-bx)) + c
                                        initialParamVariations[1] = 0.1/(xMax-xMin+1e-100);
                                        break;
                                // case POWER: // ax^b; use default for b
                                // case LOG: // a*ln(bx); use default for b
                                // case LOG2: // y = a+b*ln(x-c); use default for c
                                case RODBARD: // d+(a-d)/(1+(x/c)^b); a and d calculated by regression
                                        initialParamVariations[2] = 0.5*Math.max((xMax-xMin), Math.abs(xMean));
                                        initialParamVariations[3] = 0.5*Math.max(yMax-yMin, Math.abs(yMax));
                                        break;
                                case INV_RODBARD: // c*((x-a)/(d-x))^(1/b); c calculated by regression
                                        initialParamVariations[0] = 0.01*Math.max(xMax-xMin, Math.abs(xMax));
                                        initialParamVariations[2] = 0.1*Math.max(yMax-yMin, Math.abs(yMax));
                                        initialParamVariations[3] = 0.1*Math.max((xMax-xMin), Math.abs(xMean));
                                        break;
                                case GAMMA_VARIATE: // // b*(x-a)^c*exp(-(x-a)/d); b calculated by regression
                                // First guesses based on following observations:
                                // t0 [b] = time of first rise in gamma curve - so use the user specified first limit
                                // tm = t0 + a*B [c*d] where tm is the time of the peak of the curve
                                // therefore an estimate for a and B is sqrt(tm-t0)
                                // K [a] can now be calculated from these estimates
                                        initialParamVariations[0] = 0.1*Math.max(yMax-yMin, Math.abs(yMax));
                                        double ab = xOfMax - firstx + 0.1*(xMax-xMin);
                                        initialParamVariations[2] = 0.1*Math.sqrt(ab);
                                        initialParamVariations[3] = 0.1*Math.sqrt(ab);
                                        break;
                                case GAUSSIAN: // a + (b-a)*exp(-(x-c)^2/(2d^2)); a,b calculated by regression
                                        initialParamVariations[2] = 0.2*initialParams[3]; //(and default for d)
                                        break;
                                case GAUSSIAN_NOOFFSET: // a*exp(-(x-b)^2/(2c^2))
                                        initialParamVariations[1] = 0.2*initialParams[2]; //(and default for c)
                                        break;
                        }
                }
                return true;
        }

        /** Set multiplyParams and offsetParam for built-in functions. This allows us to use linear
         * regression, reducing the number of parameters used by the Minimizer by up to 2, and
         * improving the speed and success rate of the minimization process */
        private void getOffsetAndFactorParams() {
                offsetParam = -1;
                factorParam = -1;
                hasSlopeParam = false;
                switch (fitType) {
                        case STRAIGHT_LINE:
                        case POLY2: case POLY3: case POLY4: case POLY5: case POLY6: case POLY7: case POLY8:
                                offsetParam = 0;
                                factorParam = 1;
                                hasSlopeParam = true;
                                break;
                        case EXPONENTIAL: // a*exp(bx)
                                factorParam = 0;
                                break;
                        case EXP_WITH_OFFSET: // a*exp(-bx) + c
                        case EXP_RECOVERY: // a*(1-exp(-bx)) + c
                                offsetParam = 2;
                                factorParam = 0;
                                break;
                        case CHAPMAN:
                                offsetParam = 2;
                                factorParam = 1;
                                break;
                        case POWER: // ax^b
                                factorParam = 0;
                                break;
                        case LOG: // a*ln(bx)
                                factorParam = 0;
                                break;
                        case LOG2: // y = a+b*ln(x-c)
                                offsetParam = 0;
                                factorParam = 1;
                                break;
                        case RODBARD_INTERNAL: // d+a/(1+(x/c)^b)
                                offsetParam = 3;
                                factorParam = 0;
                                break;
                        case INV_RODBARD: // c*((x-a)/(d-x))^(1/b)
                                factorParam = 2;
                                break;
                        case GAMMA_VARIATE: // b*(x-a)^c*exp(-(x-a)/d)
                                factorParam = 1;
                                break;
                        case GAUSSIAN_INTERNAL: // a + b*exp(-(x-c)^2/(2d^2))
                                offsetParam = 0;
                                factorParam = 1;
                                break;
                        case GAUSSIAN_NOOFFSET: // a*exp(-(x-b)^2/(2c^2))
                                factorParam = 0;
                                break;  
                }
                numRegressionParams = 0;
                if (offsetParam >= 0) numRegressionParams++;
                if (factorParam >= 0) numRegressionParams++;
        }


        /** calculates the sum of y and y^2 */
        private void calculateSumYandY2() {
                sumY = 0.0; sumY2 = 0.0;
                for (int i=0; i<numPoints; i++) {
                        double y = yData[i];
                        sumY += y;
                        sumY2 += y*y;
                }
        }

        /** returns whether this a fit type that acutally fits modified data with a modified function */
        private boolean isModifiedFitType(int fitType) {
                return fitType == POWER_REGRESSION || fitType == EXP_REGRESSION || fitType == RODBARD2 ||
                                fitType == GAUSSIAN;
        }

        /** For fits don't use the original data, prepare modified data and fit type.
         * Returns false if the data are incompatible with the fit type
         * In that case, 'errorString' is set to a message explaining the problem
         */
        private boolean prepareModifiedFitType(int fitType) {
                if (fitType == GAUSSIAN) {
                        this.fitType = GAUSSIAN_INTERNAL; // different definition of parameters for using regression
                        return true;
                } else if (fitType == RODBARD) {
                        this.fitType = RODBARD_INTERNAL; // different definition of parameters for using regression
                        return true;
                } else if (fitType == POWER_REGRESSION || fitType == EXP_REGRESSION) {
                        if (fitType == POWER_REGRESSION) {
                                xDataSave = xData;
                                xData = new double[numPoints];
                        }
                        yDataSave = yData;
                        yData = new double[numPoints];
                        ySign = 0;
                        numPoints=0;  // we may have lower number of points if there is a (0,0) point that we ignore
                        for (int i=0; i<xData.length; i++) {
                                double y = yDataSave[i];
                                if (fitType == POWER_REGRESSION) {
                                        double x = xDataSave[i];
                                        if (x==0 && y==0) {
                                                restrictPower = true;
                                                continue;  // ignore (0,0) point in power law
                                        }
                                        if (x<=0) {
                                                errorString = "Cannot fit x<=0";
                                                return false;
                                        }
                                        xData[numPoints] = Math.log(x);
                                }
                                if (ySign == 0) ySign = Math.signum(y); //if unknown, determine whether y data are positive or negative
                                if (y*ySign<=0) {
                                        errorString = "Cannot fit y=0 or mixture of y>0, y<0";
                                        return false;
                                }
                                yData[numPoints] = Math.log(y*ySign);
                                numPoints++;
                        }
                        this.fitType = STRAIGHT_LINE;
                } else if (fitType == RODBARD2) { // 'Rodbard (NIH Image)' fit is inverse Rodbard function, by fitting x(y); for compatibility with NIH Image
                        xDataSave = xData;
                        yDataSave = yData;
                        xData = yDataSave; //swap
                        yData = xDataSave;
                        this.fitType = RODBARD_INTERNAL;
                }
                return true;
        }

        /** Get correct params and revert xData, yData if fit has been done via another function */
        private void postProcessModifiedFitType(int fitType) {
                if (fitType == POWER_REGRESSION || fitType == EXP_REGRESSION) // ln y = ln (a*x^b) = ln a + b ln x
                        finalParams[0] = ySign * Math.exp(finalParams[0]); //or: ln (+,-)y = ln ((+,-)a*exp(bx)) = ln (+,-)a + bx  
                if (fitType == GAUSSIAN) // a + b exp(-...) to  a + (b-a)*exp(-...)
                        finalParams[1] += finalParams[0];
                else if (fitType == RODBARD || fitType == RODBARD2) //d+a/(1+(x/c)^b) to d+(a-d)/(1+(x/c)^b)
                        finalParams[0] += finalParams[3];

                if (xDataSave != null) {
                        xData = xDataSave;
                        numPoints = xData.length;
                }
                if (yDataSave != null) yData = yDataSave;
                this.fitType = fitType;
        }

        private final double sqr(double d) { return d * d; }


        /** Pop up a dialog allowing control over simplex starting parameters */
// private void settingsDialog() {
// GenericDialog gd = new GenericDialog("Simplex Fitting Options");
// gd.addMessage("Function name: " + getName() + "\n" +
// "Formula: " + getFormula());
// char pChar = 'a';
// for (int i = 0; i < numParams; i++)
// gd.addNumericField("Initial_"+(char)(pChar+i)+":", initialParams[i], 2);
// gd.addNumericField("Maximum iterations:", minimizer.getMaxIterations(), 0);
// gd.addNumericField("Number of restarts:", minimizer.getMaxRestarts(), 0);
// gd.addNumericField("Error tolerance [1*10^(-x)]:", -(Math.log(maxRelError)/Math.log(10)), 0);
// gd.showDialog();
// if (gd.wasCanceled())
// return;
// // read initial parameters:
// for (int i = 0; i < numParams; i++) {
// double p = gd.getNextNumber();
// if (!Double.isNaN(p)) {
// initialParams[i] = p;
// initialParamVariations[i] = Math.max(0.01*p, 0.001*initialParamVariations[i]); //assume user-set params are accurate
// }
// }
// double n = gd.getNextNumber();
// if (n>0)
// minimizer.setMaxIterations((int)n);
// n = gd.getNextNumber();
// if (n>=0)
// minimizer.setMaxRestarts((int)n);
// n = gd.getNextNumber();
// setMaxError(Math.pow(10.0, -n));
// }
       
         /**
         * Gets index of highest value in an array.
         *
         * @param   Double array.
         * @return   Index of highest value.
         */
        public static int getMax(double[] array) {
                double max = array[0];
                int index = 0;
                for(int i = 1; i < array.length; i++) {
                        if(max < array[i]) {
                                max = array[i];
                                index = i;
                        }
                }
                return index;
        }

}


--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Michael Schmid]

Hi Wim,

there was no attachment.

Anyhow, the new versions of ImageJ allow you to use the curveFitter from a plugin that implements the ij.measure.UserFunction interface. Thus, you need not modify ImageJ.

CurveFitter has a call for this:
doCustomFit(UserFunction userFunction, int numParams, String formula,
  double[] initialParams, double[] initialParamVariations, boolean showSettings)

I assume that your function is: a*((1-exp(-b*t))^c).

Then, for best results, when using the CurveFitter, tell it that 'a' (parameter number 0) is a multiplicative parameter:
  setOffsetMultiplySlopeParams(-1, 0, -1)

Also, make sure that you get a reasonable guess for the 'b' parameter in the initialParams (it should be at least the correct order of magnitude). For the c parameter, you can probably specify 1 as an initial value.

---

Of course, if the Chapman-Richards equation is popular among ImageJ users, of course, it could be added to the built-in fit functions.


Michael
________________________________________________________________


On Apr 17, 2013, at 11:39, Wim Wakelkamp wrote:

> Dear ImageJ,
>
> I am working with some bio engineers on a testing program. They're really
> interested in the fitted data based on the 3-parametric Chapman-Richards
> equation F=a(1-e-bt)c , that's how i get to ImageJ. Right now i am using
> the Curvefitter.java class and trying to extends the build in fit list with
> this new formula (Chapman-Richards). When i am running doFit() and
> acquiring parameters with getparams() based on this fitType, i do not seem
> to receive to correct a,b,c values. Attached the edited Curvefitter class.
> In this class i edited everything under the CHAPMAN variable.
>
> Hope you guys can put me in the right direction.
>
> --
> Kind regards
>
> Wim
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Michael Schmid]

Hi Wim,

looking at your code now, there is one obvious problem: you should not set an offsetParam; parameter 'c' is not additive in the function.

You must set an initial value for the 'b' and 'c' parameters. This should be done in the makeInitialParamsAndVariations. If you describe a typical growth curve, you could try using the time where you reach half of the maximum as 1/b.

Michael

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html