public interface Calculateable {
    /**
     * A user-supplied function
     * @param params    When minimizing, array of variables.
     *                  For curve fit array of fit parameters.
     *                  The array contents should not be modified.
     *                  Note that the function can get an array with more
     *                  elements then needed to specify the parameters.
     *                  Ignore the rest (and don't modify them).
     * @param x         For a fit function, the independent variable of the function.
     *                  Ignore it when using the minimizer.
     * @return          The result of the function.
     */
    public double calculate(double[] params, double x);
} public class CurveFitter extends Logger implements Calculateable {
	/** 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,
			EXP_RECOVERY_NOOFFSET = 22, CHAPMAN = 23;

	/** 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,
			EXP_RECOVERY_NOOFFSET, LOG, LOG2, GAUSSIAN, GAUSSIAN_NOOFFSET, RODBARD, RODBARD2, INV_RODBARD, 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)", "Exponential Recovery (no offset)", "Chapman-Richards" }; // 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) [y = c*((x-a)/(d-x))^(1/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))", //EXP_RECOVERY_NOOFFSET
			"y = a*(1-exp(-b*x))^c" //CHAPMAN
	};

	public double[] p;
	public double x;
	public String sPatNumOut;

	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 Calculateable userFunction; // plugin with custom fit function
	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;
	}

	public double p(int i) {
		return p[i];
	}

	public double p(String sI) {
		int i = Integer.parseInt(sI);
		return p[i];
	}

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

	public void logSettings() {
		log("Simplex Fitting Options");
		log("Function name: " + getName() + "\n" + "Formula: " + getFormula());
		for (int i = 0; i < numParams; i++)
			log("p" + i + ":" + U.valToStr(initialParams[i], sPatNumOut));
		log("Maximum iterations:" + minimizer.getMaxIterations());
		log("Number of restarts:" + minimizer.getMaxRestarts());
		log("Error tolerance [1*10^(-x)]:" + (-(Math.log(maxRelError) / Math.log(10))));
	}

	/**
	 * 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 bLogSettings) {
		if (!(fitType >= STRAIGHT_LINE && fitType < fitList.length || fitType == CUSTOM))
			throw new IllegalArgumentException("Invalid fit type");
		if (fitType == CUSTOM && customFormula == 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 (customFormula != null)
				minimizer.setMaximumThreads(1); // macro interpreter does not allow multithreading
			if (!makeInitialParamsAndVariations(fitType)) // also includes some data checking
				return; // initialization failure
			if (bLogSettings) {
				logSettings();
			}
			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)) //params of actual fit to user params
			postProcessModifiedFitType(fitType);

		switch (fitType) { //postprocessing for nicer output:
		case GAUSSIAN: //Gaussians: width (std deviation) should be >0
			finalParams[3] = Math.abs(finalParams[3]);
			break;
		case GAUSSIAN_NOOFFSET:
			finalParams[2] = Math.abs(finalParams[2]);
			break;
		}
		time = System.currentTimeMillis() - startTime;
		p = getParams();
	}

	/**
	 * Fit a function defined as a evaltion String like
	 * "#$p(0)$+$p(1)$*$x$+$p(2)$*LOG($x$)#". Each $p(i)$ should only occur
	 * once. *
	 * 
	 * 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 void doCustomFit(String formula, int numParams, double[] initialParams) {
		doCustomFit(null, formula, numParams, initialParams, null, false);
	}

	/**
	 * Fit a function defined in a user plugin implementing the Calculateable
	 * 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
	 *            calculate(params, x) method. This function must allow
	 *            simultaneous calls in multiple threads.
	 * @param numParams
	 *            Number of parameters of the fit function.
	 * @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 showSettings
	 *            whether log initial parameters and a few numbers controlling
	 *            the minimizer.
	 */
	public void doCustomFit(Calculateable userFunction, int numParams, double[] initialParams) {
		doCustomFit(userFunction, null, numParams, initialParams, null, false);
	}

	/**
	 * Fit a function defined in a user plugin implementing the Calculateable
	 * 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
	 *            calculate(params, x) method.
	 *            This function must allow simultaneous calls in multiple
	 *            threads.
	 * @param numParams
	 *            Number of parameters of the fit function.
	 * @param 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
	 *            whether log initial parameters and a few numbers controlling
	 *            the minimizer.
	 */
	public void doCustomFit(Calculateable userFunction, String formula, int numParams, double[] initialParams, double[] initialParamVariations, boolean showSettings) {
		this.userFunction = userFunction;
		this.customParamCount = numParams;
		this.customFormula = formula;
		this.initialParams = initialParams;
		this.initialParamVariations = initialParamVariations;
		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 EXP_RECOVERY_NOOFFSET:
			return 2;
		case LOG:
			return 2;
		case LOG2:
			return 3;
		case GAUSSIAN_NOOFFSET:
			return 3;
		case EXP_RECOVERY:
			return 3;
		case CHAPMAN:
			return 3;
		case EXP_WITH_OFFSET:
			return 3;
		case RODBARD:
		case RODBARD2:
		case INV_RODBARD:
		case RODBARD_INTERNAL:
			return 4;
		case GAMMA_VARIATE:
			return 4;
		case GAUSSIAN:
		case GAUSSIAN_INTERNAL:
			return 4;
		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) {
		double rtn = Double.NaN;
		this.p = p;
		this.x = x;
		if (fitType != CUSTOM)
			rtn = f(fitType, p, x);
		else {
			if (userFunction != null) { // function defined in plugin
				rtn = userFunction.calculate(p, x);
			} else if (customFormula != null) { // function defined in macro
				rtn = Double.parseDouble(U.eval("#" + customFormula + "#", this));
			}
		}
		return rtn;
	}

	/**
	 * 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 EXP_RECOVERY_NOOFFSET:
			return p[0] * (1 - Math.exp(-p[1] * x));
		case CHAPMAN: // a*(1-exp(-b*x))^c
			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: " + U.valToStr(getSumResidualsSqr(), sPatNumOut) + "\nStandard deviation: " + U.valToStr(getSD(), sPatNumOut) + "\nR^2: "
				+ U.valToStr(getRSquared(), sPatNumOut) + "\nParameters:";
		char pChar = 'a';
		double[] pVal = getParams();
		for (int i = 0; i < numParams; i++) {
			resultS += "\n	" + pChar + " = " + U.valToStr(pVal[i], sPatNumOut);
			pChar++;
		}
		return resultS;
	}

	/**
	 * Return the result formula when the fitting is done using customized
	 * formula (via doCustomFit(String formula, int numParams, double[]
	 * initialParams) or doCustomFit(Calculateable userFunction, String formula,
	 * int numParams, double[] initialParams, double[] initialParamVariations,
	 * boolean showSettings) with a not null formula). When customFormula==null,
	 * the function returns null.
	 */
	public String getResultFormula() {
		String rtn = customFormula;
		if (rtn != null) {
			rtn = rtn.replace("$x$", "x");
			rtn = U.eval(rtn, this, sPatNumOut);
		}
		return rtn;
	}

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

	/**
	 * Create output on the number of iterations in the ImageJ Status line, e.g.
	 * "<ijStatusString> 50 (max 750); ESC to stop"
	 * 
	 * @param ijStatusString
	 *            Displayed in the beginning of the status message. No display
	 *            if null.
	 *            E.g. "Curve Fit: 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) {
		minimizer.setStatusAndEsc(ijStatusString, checkEscape);
	}

	/**
	 * 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 calculate(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]);
			}
			//IJ.log(IJ.d2s(params[0],3,5)+","+IJ.d2s(params[1],3,5)+": r="+IJ.d2s(sumResidualsSqr,3,5)+Thread.currentThread().getName() );

		} 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 (fitType == CUSTOM) {
				for (int i = 0; i < numParams; i++)
					initialParams[i] = 1.0;
			}
		}
		if (fitType == CUSTOM)
			return true; // no way to guess initial parameters or initialParamVariations
		if (!hasInitialParamVariations)
			initialParamVariations = new double[numParams];

		// 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; //with this initialization, the loop starts at 1, not 0
		double yMin = firsty, yMax = firsty;
		double xMean = firstx, yMean = firsty;
		double xOfMax = firstx;
		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 || fitType == CHAPMAN)) {
			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.0 / (xMax - xMin + 1e-100) * Math.signum(yMean) * Math.signum(slope);
				initialParams[0] = yMean / Math.exp(initialParams[1] * xMean); //don't care, done by regression
				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 EXP_RECOVERY_NOOFFSET: // a*(1-exp(-bx))
				initialParams[1] = 1.0 / (xMax - xMin + 1e-100) * Math.signum(yMean) * Math.signum(slope);
				initialParams[0] = yMean / Math.exp(initialParams[1] * xMean); //don't care, done by 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 CHAPMAN: // a*(1-exp(-b*x))^c
				initialParams[0] = yMax;
				initialParams[2] = 1.5; // just assuming any reasonable value
				for (int i = 1; i < numPoints; i++)
					//find when we reach 50% of maximum
					if (yData[i] > 0.5 * yMax) { //approximately (1-exp(-1))^1.5 = 0.5
						initialParams[1] = 1. / xData[i];
						break;
					}
				if (Double.isNaN(initialParams[1]) || initialParams[1] > 1000. / xMax) //in case an outlier at the beginning has fooled us
					initialParams[1] = 10. / xMax;
				break;
			//no case CUSTOM: here, was done above
			}
		}
		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 CHAPMAN:            // a*(1-exp(-b*x))^c use default (10% of value) for b, c
			// 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 EXP_RECOVERY_NOOFFSET: // a*(1-exp(-bx))
			factorParam = 0;
			break;
		case CHAPMAN: // a*(1-exp(-b*x))^c
			factorParam = 0;
			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;
	}

	/**
	 * 
	 * @param initialParams
	 * @param maxIterations
	 *            Maximum iterations
	 * @param maxRestarts
	 *            Number of restarts
	 * @param maxErrorCoefficient
	 *            Error tolerance [1*10^(-x)]
	 */
	private void settInitPars(double initialParams[], int maxIterations, int maxRestarts, int maxErrorCoefficient) {
		if (initialParamVariations == null)
			initialParamVariations = new double[numParams];
		if (initialParams != null) {
			for (int i = 0; i < numParams; i++) {
				double p = initialParams[i];
				if (!Double.isNaN(p)) {
					this.initialParams[i] = p;
					initialParamVariations[i] = Math.max(0.01 * p, 0.001 * initialParamVariations[i]); //assume user-set params are accurate
				}
			}
		}
		if (maxIterations > 0)
			minimizer.setMaxIterations(maxIterations);
		if (maxRestarts >= 0)
			minimizer.setMaxRestarts(maxRestarts);
		setMaxError(Math.pow(10.0, -maxErrorCoefficient));
	}

	/**
	 * Gets index of highest value in an array.
	 * 
	 * @param array
	 *            the 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;
	}

} public class Binomial extends AbstractDiscreteDistribution {
	protected int n;
	protected double p;

	// cache vars for method generateBinomial(...)
	private int n_last = -1, n_prev = -1;
	private double par, np, p0, q, p_last = -1.0, p_prev = -1.0;
	private int b, m, nm;
	private double pq, rc, ss, xm, xl, xr, ll, lr, c, p1, p2, p3, p4, ch;

	// cache vars for method pdf(...)
	private double log_p, log_q, log_n;

	// The uniform random number generated shared by all <b>static</b> methods.
	protected static Binomial shared = new Binomial(1, 0.5, makeDefaultGenerator());
	
	public Binomial(RandomEngine randomGenerator) {
		setRandomGenerator(randomGenerator);
	}

	/**
	 * Constructs a binomial distribution.
	 * Example: n=1, p=0.5.
	 * 
	 * @param n
	 *            the number of trials (also known as <i>sample size</i>).
	 * @param p
	 *            the probability of success.
	 * @param randomGenerator
	 *            a uniform random number generator.
	 * @throws IllegalArgumentException
	 *             if <tt>n*Math.min(p,1-p) &lt;= 0.0</tt>
	 */
	public Binomial(int n, double p, RandomEngine randomGenerator) {
		setRandomGenerator(randomGenerator);
		setNandP(n, p);
	}

	/**
	 * Returns the cumulative distribution function.
	 */
	public double cdf(int k) {
		return Probability.binomial(k, n, p);
	}

	/**
	 * Returns the cumulative distribution function.
	 */
	private double cdfSlow(int k) {
		if (k < 0)
			throw new IllegalArgumentException();

		double sum = 0.0;
		for (int r = 0; r <= k; r++)
			sum += pdf(r);

		return sum;
	}

	/******************************************************************
	 * *
	 * Binomial-Distribution - Acceptance Rejection/Inversion *
	 * *
	 ****************************************************************** 
	 * *
	 * Acceptance Rejection method combined with Inversion for *
	 * generating Binomial random numbers with parameters *
	 * n (number of trials) and p (probability of success). *
	 * For min(n*p,n*(1-p)) < 10 the Inversion method is applied: *
	 * The random numbers are generated via sequential search, *
	 * starting at the lowest index k=0. The cumulative probabilities *
	 * are avoided by using the technique of chop-down. *
	 * For min(n*p,n*(1-p)) >= 10 Acceptance Rejection is used: *
	 * The algorithm is based on a hat-function which is uniform in *
	 * the centre region and exponential in the tails. *
	 * A triangular immediate acceptance region in the centre speeds *
	 * up the generation of binomial variates. *
	 * If candidate k is near the mode, f(k) is computed recursively *
	 * starting at the mode m. *
	 * The acceptance test by Stirling's formula is modified *
	 * according to W. Hoermann (1992): The generation of binomial *
	 * random variates, to appear in J. Statist. Comput. Simul. *
	 * If p < .5 the algorithm is applied to parameters n, p. *
	 * Otherwise p is replaced by 1-p, and k is replaced by n - k. *
	 * *
	 ****************************************************************** 
	 * *
	 * FUNCTION: - samples a random number from the binomial *
	 * distribution with parameters n and p and is *
	 * valid for n*min(p,1-p) > 0. *
	 * REFERENCE: - V. Kachitvichyanukul, B.W. Schmeiser (1988): *
	 * Binomial random variate generation, *
	 * Communications of the ACM 31, 216-222. *
	 * SUBPROGRAMS: - StirlingCorrection() *
	 * ... Correction term of the Stirling *
	 * approximation for log(k!) *
	 * (series in 1/k or table values *
	 * for small k) with long int k *
	 * - randomGenerator ... (0,1)-Uniform engine *
	 * *
	 ******************************************************************/
	public int generateBinomial(int n, double p) {
		final double C1_3 = 0.33333333333333333;
		final double C5_8 = 0.62500000000000000;
		final double C1_6 = 0.16666666666666667;
		final int DMAX_KM = 20;

		int bh, i, K, Km, nK;
		double f, rm, U, V, X, T, E;

		if (n != n_last || p != p_last) { // set-up 
			n_last = n;
			p_last = p;
			par = Math.min(p, 1.0 - p);
			q = 1.0 - par;
			np = n * par;

			// Check for invalid input values

			if (np <= 0.0)
				return -1;

			rm = np + par;
			m = (int) rm; // mode, integer 
			if (np < 10) {
				p0 = Math.exp(n * Math.log(q)); // Chop-down
				bh = (int) (np + 10.0 * Math.sqrt(np * q));
				b = Math.min(n, bh);
			} else {
				rc = (n + 1.0) * (pq = par / q); // recurr. relat.
				ss = np * q; // variance  
				i = (int) (2.195 * Math.sqrt(ss) - 4.6 * q); // i = p1 - 0.5
				xm = m + 0.5;
				xl = (double) (m - i); // limit left
				xr = (double) (m + i + 1L); // limit right
				f = (rm - xl) / (rm - xl * par);
				ll = f * (1.0 + 0.5 * f);
				f = (xr - rm) / (xr * q);
				lr = f * (1.0 + 0.5 * f);
				c = 0.134 + 20.5 / (15.3 + (double) m); // parallelogram
														// height
				p1 = i + 0.5;
				p2 = p1 * (1.0 + c + c); // probabilities
				p3 = p2 + c / ll; // of regions 1-4
				p4 = p3 + c / lr;
			}
		}

		if (np < 10) { //Inversion Chop-down
			double pk;

			K = 0;
			pk = p0;
			U = randomGenerator.raw();
			while (U > pk) {
				++K;
				if (K > b) {
					U = randomGenerator.raw();
					K = 0;
					pk = p0;
				} else {
					U -= pk;
					pk = (double) (((n - K + 1) * par * pk) / (K * q));
				}
			}
			return ((p > 0.5) ? (n - K) : K);
		}

		for (;;) {
			V = randomGenerator.raw();
			if ((U = randomGenerator.raw() * p4) <= p1) { // triangular region
				K = (int) (xm - U + p1 * V);
				return (p > 0.5) ? (n - K) : K; // immediate accept
			}
			if (U <= p2) { // parallelogram
				X = xl + (U - p1) / c;
				if ((V = V * c + 1.0 - Math.abs(xm - X) / p1) >= 1.0)
					continue;
				K = (int) X;
			} else if (U <= p3) { // left tail
				if ((X = xl + Math.log(V) / ll) < 0.0)
					continue;
				K = (int) X;
				V *= (U - p2) * ll;
			} else { // right tail
				if ((K = (int) (xr - Math.log(V) / lr)) > n)
					continue;
				V *= (U - p3) * lr;
			}

			// acceptance test :  two cases, depending on |K - m|
			if ((Km = Math.abs(K - m)) <= DMAX_KM || Km + Km + 2L >= ss) {

				// computation of p(K) via recurrence relationship from the mode
				f = 1.0; // f(m)
				if (m < K) {
					for (i = m; i < K;) {
						if ((f *= (rc / ++i - pq)) < V)
							break; // multiply  f
					}
				} else {
					for (i = K; i < m;) {
						if ((V *= (rc / ++i - pq)) > f)
							break; // multiply  V
					}
				}
				if (V <= f)
					break; // acceptance test
			} else {

				// lower and upper squeeze tests, based on lower bounds for log p(K)
				V = Math.log(V);
				T = -Km * Km / (ss + ss);
				E = (Km / ss) * ((Km * (Km * C1_3 + C5_8) + C1_6) / ss + 0.5);
				if (V <= T - E)
					break;
				if (V <= T + E) {
					if (n != n_prev || par != p_prev) {
						n_prev = n;
						p_prev = par;

						nm = n - m + 1;
						ch = xm * Math.log((m + 1.0) / (pq * nm)) + Arithmetic.stirlingCorrection(m + 1) + Arithmetic.stirlingCorrection(nm);
					}
					nK = n - K + 1;

					// computation of log f(K) via Stirling's formula
					// final acceptance-rejection test
					if (V <= ch + (n + 1.0) * Math.log((double) nm / (double) nK) + (K + 0.5) * Math.log(nK * pq / (K + 1.0)) - Arithmetic.stirlingCorrection(K + 1)
							- Arithmetic.stirlingCorrection(nK))
						break;
				}
			}
		}
		return (p > 0.5) ? (n - K) : K;
	}
	
	/**
	 * Luc Devroye's "Second Waiting Time Method" on page 522 of his text
	 * "Non-Uniform Random Variate Generation." This function works fast
	 * when n*p is small. It also returns random binomial numbers with correct
	 * distribution even when n*p is small. When n*p<1E-2 Acceptance
	 * Rejection/Inversion in
	 * generateBinomial(int n, double p) fails to generate numbers with right
	 * distribution!
	 * 
	 * @author leoking
	 * @param n
	 * @param p
	 * @return
	 */
	public int generateBinomialS(int n, double p) {
		if (p > 0.5) {
			return n - generateBinomialS(n, 1 - p);
		}
		double log_q = Math.log(1.0 - p);
		int x = 0;
		double sum = 0;
		while (true) {
			sum += Math.log(randomGenerator.nextDouble()) / (n - x);
			if (sum < log_q) {
				break;
			}
			x++;
		}
		return x;
	}

	/**
	 * Returns a random number from the distribution using Acceptance Rejection/Inversion method.
	 * This function is fast but can not generate number with right distribution when n*p<1E-2.
	 */
	public int nextInt() {
		return generateBinomial(n, p);
	}
	
	/**
	 * Returns a random number from the distribution using Luc Devroye's
	 * "Second Waiting Time Method".
	 * This function works rigth and fast when n*p is small.
	 * @author leoking
	 * @return
	 */
	public int nextIntS() {
		return generateBinomialS(n, p);
	}

	/**
	 * Returns a random number from the distribution with the given parameters n
	 * and p; bypasses the internal state.
	 * 
	 * @param n
	 *            the number of trials
	 * @param p
	 *            the probability of success.
	 * @throws IllegalArgumentException
	 *             if <tt>n*Math.min(p,1-p) &lt;= 0.0</tt>
	 */
	public int nextInt(int n, double p) {
		if (n * Math.min(p, 1 - p) <= 0.0)
			throw new IllegalArgumentException();
		return generateBinomial(n, p);
	}

	/**
	 * Returns the probability distribution function.
	 */
	public double pdf(int k) {
		if (k < 0)
			throw new IllegalArgumentException();
		int r = this.n - k;
		return Math.exp(this.log_n - Arithmetic.logFactorial(k) - Arithmetic.logFactorial(r) + this.log_p * k + this.log_q * r);
	}

	/**
	 * Sets the parameters number of trials and the probability of success.
	 * 
	 * @param n
	 *            the number of trials
	 * @param p
	 *            the probability of success.
	 * @throws IllegalArgumentException
	 *             if <tt>n*Math.min(p,1-p) &lt;= 0.0</tt>
	 */
	public void setNandP(int n, double p) {
		if (n * Math.min(p, 1 - p) <= 0.0)
			throw new IllegalArgumentException();
		this.n = n;
		this.p = p;

		this.log_p = Math.log(p);
		this.log_q = Math.log(1.0 - p);
		this.log_n = Arithmetic.logFactorial(n);
	}

	/**
	 * Returns a random number from the distribution with the given parameters n
	 * and p.
	 * 
	 * @param n
	 *            the number of trials
	 * @param p
	 *            the probability of success.
	 * @throws IllegalArgumentException
	 *             if <tt>n*Math.min(p,1-p) &lt;= 0.0</tt>
	 */
	public static int staticNextInt(int n, double p) {
		synchronized (shared) {
			return shared.nextInt(n, p);
		}
	}

	/**
	 * Returns a String representation of the receiver.
	 */
	public String toString() {
		return this.getClass().getName() + "(" + n + "," + p + ")";
	}

	/**
	 * Sets the uniform random number generated shared by all <b>static</b>
	 * methods.
	 * 
	 * @param randomGenerator
	 *            the new uniform random number generator to be shared.
	 */
	private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
		synchronized (shared) {
			shared.setRandomGenerator(randomGenerator);
		}
	}
}
 public class DRand extends RandomEngine {
	private int current;
	public static final int DEFAULT_SEED = 1;

	/**
	 * Constructs and returns a random number generator with a default seed,
	 * which is a <b>constant</b>.
	 */
	public DRand() {
		this(DEFAULT_SEED);
	}

	/**
	 * Constructs and returns a random number generator with the given seed.
	 * 
	 * @param seed
	 *            should not be 0, in such a case <tt>DRand.DEFAULT_SEED</tt> is
	 *            substituted.
	 */
	public DRand(long seed) {
		setSeed(seed);
	}

	/**
	 * Constructs and returns a random number generator seeded with the given
	 * date.
	 * 
	 * @param d
	 *            typically <tt>new java.util.Date()</tt>
	 */
	public DRand(Date d) {
		this(d.getTime());
	}

	/**
	 * Returns a 32 bit uniformly distributed random number in the closed
	 * interval <tt>[Integer.MIN_VALUE,Integer.MAX_VALUE]</tt> (including
	 * <tt>Integer.MIN_VALUE</tt> and <tt>Integer.MAX_VALUE</tt>).
	 */
	public int nextInt() {
		current *= 0x278DDE6D; /* z(i+1)=a*z(i) (mod 2**32) */
		// a == 0x278DDE6D == 663608941

		return current;
	}

	/**
	 * Sets the receiver's seed.
	 * This method resets the receiver's entire internal state.
	 * The following condition must hold:
	 * <tt>seed &gt;= 0 && seed &lt; (2<sup>32</sup>-1) / 4</tt>.
	 * 
	 * @param seed
	 *            if the above condition does not hold, a modified seed that
	 *            meets the condition is silently substituted.
	 */
	public void setSeed(long seed) {
		if (seed < 0)
			seed = -seed;
		int limit = (int) ((Math.pow(2, 32) - 1) / 4); // --> 536870911
		if (seed >= limit)
			seed = seed >> 3;

		this.current = (int) (4 * seed + 1);
	}
}
 public class Exponential extends AbstractContinousDistribution { 
	protected double lambda;

 	// The uniform random number generated shared by all <b>static</b> methods.
	protected static Exponential shared = new Exponential(1.0,makeDefaultGenerator());
/**
 * Constructs a Negative Exponential distribution.
 */
public Exponential(double lambda, RandomEngine randomGenerator) {
	setRandomGenerator(randomGenerator);
	setState(lambda);
}
/**
 * Returns the cumulative distribution function.
 */
public double cdf(double x) {
	if (x <= 0.0) return 0.0;
	return 1.0 - Math.exp(-x * lambda);
}
/**
 * Returns a random number from the distribution.
 */
public double nextDouble() {
	return nextDouble(lambda);
}
/**
 * Returns a random number from the distribution; bypasses the internal state.
 */
public double nextDouble(double lambda) {
	return - Math.log(randomGenerator.raw()) / lambda;
}
/**
 * Returns the probability distribution function.
 */
public double pdf(double x) {
	if (x < 0.0) return 0.0;
	return lambda*Math.exp(-x*lambda);
}
/**
 * Sets the mean.
 */
public void setState(double lambda) {
	this.lambda = lambda;
}
/**
 * Returns a random number from the distribution with the given lambda.
 */
public static double staticNextDouble(double lambda) {
	synchronized (shared) {
		return shared.nextDouble(lambda);
	}
}
/**
 * Returns a String representation of the receiver.
 */
public String toString() {
	return this.getClass().getName()+"("+lambda+")";
}
/**
 * Sets the uniform random number generated shared by all <b>static</b> methods.
 * @param randomGenerator the new uniform random number generator to be shared.
 */
private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
	synchronized (shared) {
		shared.setRandomGenerator(randomGenerator);
	}
}
}
 public class Gamma extends AbstractContinousDistribution { 
	protected double alpha;
	protected double lambda;

 	// The uniform random number generated shared by all <b>static</b> methods.
	protected static Gamma shared = new Gamma(1.0,1.0,makeDefaultGenerator());
/**
 * Constructs a Gamma distribution.
 * Example: alpha=1.0, lambda=1.0.
 * @throws IllegalArgumentException if <tt>alpha &lt;= 0.0 || lambda &lt;= 0.0</tt>.
 */
public Gamma(double alpha, double lambda, RandomEngine randomGenerator) {
	setRandomGenerator(randomGenerator);
	setState(alpha,lambda);
}
/**
 * Returns the cumulative distribution function.
 */
public double cdf(double x) {
	return Probability.gamma(alpha,lambda,x);
}
/**
 * Returns a random number from the distribution.
 */
public double nextDouble() {
	return nextDouble(alpha, lambda);
}
/**
 * Returns a random number from the distribution; bypasses the internal state.
 */
public double nextDouble(double alpha, double lambda) {
/******************************************************************
 *                                                                *
 *    Gamma Distribution - Acceptance Rejection combined with     *
 *                         Acceptance Complement                  *
 *                                                                *
 ******************************************************************
 *                                                                *
 * FUNCTION:    - gds samples a random number from the standard   *
 *                gamma distribution with parameter  a > 0.       *
 *                Acceptance Rejection  gs  for  a < 1 ,          *
 *                Acceptance Complement gd  for  a >= 1 .         *
 * REFERENCES:  - J.H. Ahrens, U. Dieter (1974): Computer methods *
 *                for sampling from gamma, beta, Poisson and      *
 *                binomial distributions, Computing 12, 223-246.  *
 *              - J.H. Ahrens, U. Dieter (1982): Generating gamma *
 *                variates by a modified rejection technique,     *
 *                Communications of the ACM 25, 47-54.            *
 * SUBPROGRAMS: - drand(seed) ... (0,1)-Uniform generator with    *
 *                unsigned long integer *seed                     *
 *              - NORMAL(seed) ... Normal generator N(0,1).       *
 *                                                                *
 ******************************************************************/
 	double a = alpha;
	double aa = -1.0, aaa = -1.0, 
		b=0.0, c=0.0, d=0.0, e, r, s=0.0, si=0.0, ss=0.0, q0=0.0,
		q1 = 0.0416666664, q2 =  0.0208333723, q3 = 0.0079849875,
		q4 = 0.0015746717, q5 = -0.0003349403, q6 = 0.0003340332,
		q7 = 0.0006053049, q8 = -0.0004701849, q9 = 0.0001710320,
		a1 = 0.333333333,  a2 = -0.249999949,  a3 = 0.199999867,
		a4 =-0.166677482,  a5 =  0.142873973,  a6 =-0.124385581,
		a7 = 0.110368310,  a8 = -0.112750886,  a9 = 0.104089866,
		e1 = 1.000000000,  e2 =  0.499999994,  e3 = 0.166666848,
		e4 = 0.041664508,  e5 =  0.008345522,  e6 = 0.001353826,
		e7 = 0.000247453;

	double gds,p,q,t,sign_u,u,v,w,x;
	double v1,v2,v12;

	// Check for invalid input values

	if (a <= 0.0) throw new IllegalArgumentException(); 
	if (lambda <= 0.0) new IllegalArgumentException(); 

	if (a < 1.0) { // CASE A: Acceptance rejection algorithm gs
		b = 1.0 + 0.36788794412 * a;              // Step 1
		for(;;) {
			p = b * randomGenerator.raw();
			if (p <= 1.0) {                       // Step 2. Case gds <= 1
				gds = Math.exp(Math.log(p) / a);
				if (Math.log(randomGenerator.raw()) <= -gds) return(gds/lambda);
			}
			else {                                // Step 3. Case gds > 1
				gds = - Math.log ((b - p) / a);
				if (Math.log(randomGenerator.raw()) <= ((a - 1.0) * Math.log(gds))) return(gds/lambda);
			}
		}
	}

	else {        // CASE B: Acceptance complement algorithm gd (gaussian distribution, box muller transformation)
		if (a != aa) {                        // Step 1. Preparations
			aa = a;
			ss = a - 0.5;
			s = Math.sqrt(ss);
			d = 5.656854249 - 12.0 * s;
		}
												  // Step 2. Normal deviate
		do {
			v1 = 2.0 * randomGenerator.raw() - 1.0;
			v2 = 2.0 * randomGenerator.raw() - 1.0;
			v12 = v1*v1 + v2*v2;
		} while ( v12 > 1.0 );
		t = v1*Math.sqrt(-2.0*Math.log(v12)/v12);
		x = s + 0.5 * t;
		gds = x * x;
		if (t >= 0.0) return(gds/lambda);         // Immediate acceptance

		u = randomGenerator.raw();                // Step 3. Uniform random number
		if (d * u <= t * t * t) return(gds/lambda); // Squeeze acceptance

		if (a != aaa) {                           // Step 4. Set-up for hat case
			aaa = a;
			r = 1.0 / a;
			q0 = ((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) *
					  r + q3) * r + q2) * r + q1) * r;
			if (a > 3.686) {
				if (a > 13.022) {
					b = 1.77;
					si = 0.75;
					c = 0.1515 / s;
				}
			    else {
					b = 1.654 + 0.0076 * ss;
					si = 1.68 / s + 0.275;
					c = 0.062 / s + 0.024;
				}
			}
			else {
				b = 0.463 + s - 0.178 * ss;
				si = 1.235;
				c = 0.195 / s - 0.079 + 0.016 * s;
			}
		}
		if (x > 0.0) {                        // Step 5. Calculation of q
			v = t / (s + s);                  // Step 6.
			if (Math.abs(v) > 0.25) {
				q = q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log(1.0 + v);
			}
			else {
				q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) *
				    v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
			}								  // Step 7. Quotient acceptance
			if (Math.log(1.0 - u) <= q) return(gds/lambda);
		}

		for(;;) {              			      // Step 8. Double exponential deviate t
			do {
				e = -Math.log(randomGenerator.raw());
				u = randomGenerator.raw();
				u = u + u - 1.0;
				sign_u = (u > 0)? 1.0 : -1.0;
				t = b + (e * si) * sign_u;
			} while (t <= -0.71874483771719); // Step 9. Rejection of t
			v = t / (s + s);                  // Step 10. New q(t)
			if (Math.abs(v) > 0.25) {
				q = q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log(1.0 + v);
			}
			else {
				q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) *
				    v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
			}
			if (q <= 0.0) continue;           // Step 11.
			if (q > 0.5) {
				w = Math.exp(q) - 1.0;
			}
			else {
				w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) *
					     q + e1) * q;
			}                    			  // Step 12. Hat acceptance
			if ( c * u * sign_u <= w * Math.exp(e - 0.5 * t * t)) {
				x = s + 0.5 * t;
				return(x*x/lambda);
			}
		}
	}
}
/**
 * Returns the probability distribution function.
 */
public double pdf(double x) {
	if (x < 0) throw new IllegalArgumentException();
	if (x == 0) {
		if (alpha == 1.0) return 1.0/lambda;
		else return 0.0;
	}
	if (alpha == 1.0) return Math.exp(-x/lambda)/lambda;

	return Math.exp((alpha-1.0) * Math.log(x/lambda) - x/lambda - Fun.logGamma(alpha)) / lambda;
}
/**
 * Sets the mean and variance.
 * @throws IllegalArgumentException if <tt>alpha &lt;= 0.0 || lambda &lt;= 0.0</tt>.
 */
public void setState(double alpha, double lambda) {
	if (alpha <= 0.0) throw new IllegalArgumentException();
	if (lambda <= 0.0) throw new IllegalArgumentException();
	this.alpha = alpha;
	this.lambda = lambda;
}
/**
 * Returns a random number from the distribution.
 * @throws IllegalArgumentException if <tt>alpha &lt;= 0.0 || lambda &lt;= 0.0</tt>.
 */
public static double staticNextDouble(double alpha, double lambda) {
	synchronized (shared) {
		return shared.nextDouble(alpha,lambda);
	}
}
/**
 * Returns a String representation of the receiver.
 */
public String toString() {
	return this.getClass().getName()+"("+alpha+","+lambda+")";
}
/**
 * Sets the uniform random number generated shared by all <b>static</b> methods.
 * @param randomGenerator the new uniform random number generator to be shared.
 */
private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
	synchronized (shared) {
		shared.setRandomGenerator(randomGenerator);
	}
}
}
 public class Normal extends AbstractContinousDistribution {
	protected double mean;
	protected double variance;
	protected double standardDeviation;

	protected double cache; // cache for Box-Mueller algorithm 
	protected boolean cacheFilled; // Box-Mueller

	protected double SQRT_INV; // performance cache

 	// The uniform random number generated shared by all <b>static</b> methods.
	protected static Normal shared = new Normal(0.0,1.0,makeDefaultGenerator());
/**
 * Constructs a normal (gauss) distribution.
 * Example: mean=0.0, standardDeviation=1.0.
 */
public Normal(double mean, double standardDeviation, RandomEngine randomGenerator) {
	setRandomGenerator(randomGenerator);
	setState(mean,standardDeviation);
}
/**
 * Returns the cumulative distribution function.
 */
public double cdf(double x) {
	return Probability.normal(mean,variance,x);
}
/**
 * Returns a random number from the distribution.
 */
public double nextDouble() {
	return nextDouble(this.mean,this.standardDeviation);
}
/**
 * Returns a random number from the distribution; bypasses the internal state.
 */
public double nextDouble(double mean, double standardDeviation) {
	// Uses polar Box-Muller transformation.
	if (cacheFilled && this.mean == mean && this.standardDeviation == standardDeviation) {
		cacheFilled = false;
		return cache; 
	};

	double x,y,r,z;
	do {
		x = 2.0*randomGenerator.raw() - 1.0; 
		y = 2.0*randomGenerator.raw() - 1.0;		 
		r = x*x+y*y;
	} while (r >= 1.0);

	z = Math.sqrt(-2.0*Math.log(r)/r);
	cache = mean + standardDeviation*x*z;
	cacheFilled = true;
	return mean + standardDeviation*y*z;
}
/**
 * Returns the probability distribution function.
 */
public double pdf(double x) {
	double diff = x-mean;
	return SQRT_INV * Math.exp(-(diff*diff) / (2.0*variance));
}
/**
 * Sets the uniform random generator internally used.
 */
protected void setRandomGenerator(RandomEngine randomGenerator) {
	super.setRandomGenerator(randomGenerator);
	this.cacheFilled = false;
}
/**
 * Sets the mean and variance.
 */
public void setState(double mean, double standardDeviation) {
	if (mean!=this.mean || standardDeviation!=this.standardDeviation) {
		this.mean = mean;
		this.standardDeviation = standardDeviation;
		this.variance = standardDeviation*standardDeviation;
		this.cacheFilled = false;
		
		this.SQRT_INV = 1.0 / Math.sqrt(2.0*Math.PI*variance);
	}
}
/**
 * Returns a random number from the distribution with the given mean and standard deviation.
 */
public static double staticNextDouble(double mean, double standardDeviation) {
	synchronized (shared) {
		return shared.nextDouble(mean,standardDeviation);
	}
}
/**
 * Returns a String representation of the receiver.
 */
public String toString() {
	return this.getClass().getName()+"("+mean+","+standardDeviation+")";
}
/**
 * Sets the uniform random number generated shared by all <b>static</b> methods.
 * @param randomGenerator the new uniform random number generator to be shared.
 */
private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
	synchronized (shared) {
		shared.setRandomGenerator(randomGenerator);
	}
}
}
 public class Poisson extends AbstractDiscreteDistribution {
	protected double mean;

	// precomputed and cached values (for performance only)
	// cache for < SWITCH_MEAN
	protected double my_old = -1.0;
	protected double p,q,p0;
	protected double[] pp = new double[36];
	protected int llll;

	// cache for >= SWITCH_MEAN
	protected double my_last = -1.0;
	protected double ll;
	protected int k2, k4, k1, k5;
	protected double dl, dr, r1, r2, r4, r5, lr, l_my, c_pm;
	protected double f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;

	// cache for both;
	protected int m;


	protected static final double MEAN_MAX = Integer.MAX_VALUE; // for all means larger than that, we don't try to compute a poisson deviation, but return the mean.
	protected static final double SWITCH_MEAN = 10.0; // switch from method A to method B
	

 	// The uniform random number generated shared by all <b>static</b> methods. 
	protected static Poisson shared = new Poisson(0.0,makeDefaultGenerator());
/**
 * Constructs a poisson distribution.
 * Example: mean=1.0.
 */
public Poisson(double mean, RandomEngine randomGenerator) {
	setRandomGenerator(randomGenerator);
	setMean(mean);
}
/**
 * Returns the cumulative distribution function.
 */
public double cdf(int k) {
	return Probability.poisson(k,this.mean);
}
/**
 * Returns a deep copy of the receiver; the copy will produce identical sequences.
 * After this call has returned, the copy and the receiver have equal but separate state.
 *
 * @return a copy of the receiver.
 */
public Object clone() {
	Poisson copy = (Poisson) super.clone();
	if (this.pp != null) copy.pp = (double[]) this.pp.clone();
	return copy;
}
private static double f(int k, double l_nu, double c_pm) {
	return  Math.exp(k * l_nu - Arithmetic.logFactorial(k) - c_pm);
}
/**
 * Returns a random number from the distribution.
 */
public int nextInt() {
	return nextInt(this.mean);
}
/**
 * Returns a random number from the distribution; bypasses the internal state.
 */
public int nextInt(double theMean) {
/******************************************************************
 *                                                                *
 * Poisson Distribution - Patchwork Rejection/Inversion           *
 *                                                                *
 ******************************************************************
 *                                                                *
 * For parameter  my < 10  Tabulated Inversion is applied.        *
 * For my >= 10  Patchwork Rejection is employed:                 *
 * The area below the histogram function f(x) is rearranged in    *
 * its body by certain point reflections. Within a large center   *
 * interval variates are sampled efficiently by rejection from    *
 * uniform hats. Rectangular immediate acceptance regions speed   *
 * up the generation. The remaining tails are covered by          *
 * exponential functions.                                         *
 *                                                                *
 *****************************************************************/
	RandomEngine gen = this.randomGenerator;
	double my = theMean;
	
	double t,g,my_k;

	double gx,gy,px,py,e,x,xx,delta,v;
	int sign;

	//static double p,q,p0,pp[36];
	//static long ll,m;
	double u;
	int k,i;
 
	if (my < SWITCH_MEAN) { // CASE B: Inversion- start new table and calculate p0
		if (my != my_old) {
			my_old = my;
			llll = 0;
			p = Math.exp(-my);
			q = p;
			p0 = p;
			//for (k=pp.length; --k >=0; ) pp[k] = 0;
		}
		m = (my > 1.0) ? (int)my : 1;
		for(;;) {
			u = gen.raw();           // Step U. Uniform sample 
			k = 0;
			if (u <= p0) return(k);
			if (llll != 0) {              // Step T. Table comparison 
				i = (u > 0.458) ? Math.min(llll,m) : 1;
				for (k = i; k <=llll; k++) if (u <= pp[k]) return(k);
				if (llll == 35) continue;
			}
			for (k = llll +1; k <= 35; k++) { // Step C. Creation of new prob. 
				p *= my/(double)k;
				q += p;
				pp[k] = q;
				if (u <= q) {
					llll = k;
					return(k);
				}
			}
			llll = 35;
		}
	}     // end my < SWITCH_MEAN 
	else if (my < MEAN_MAX ) { // CASE A: acceptance complement
		//static double        my_last = -1.0;
		//static long int      m,  k2, k4, k1, k5;
		//static double        dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm,
		//  					 f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
		int    Dk, X, Y;
		double Ds, U, V, W;

		m  = (int) my;
		if (my != my_last) { //  set-up    
			my_last = my;

			// approximate deviation of reflection points k2, k4 from my - 1/2    
			Ds = Math.sqrt(my + 0.25);

			// mode m, reflection points k2 and k4, and points k1 and k5, which    
			// delimit the centre region of h(x)                                    
			k2 = (int) Math.ceil(my - 0.5 - Ds);
			k4 = (int)     (my - 0.5 + Ds);
			k1 = k2 + k2 - m + 1;
			k5 = k4 + k4 - m;

			// range width of the critical left and right centre region            
			dl = (double) (k2 - k1);
			dr = (double) (k5 - k4);

			// recurrence constants r(k) = p(k)/p(k-1) at k = k1, k2, k4+1, k5+1    
			r1 = my / (double) k1;
			r2 = my / (double) k2;
			r4 = my / (double)(k4 + 1);
			r5 = my / (double)(k5 + 1);

			// reciprocal values of the scale parameters of expon. tail envelopes   
			ll =  Math.log(r1);                     // expon. tail left 
			lr = -Math.log(r5);                     // expon. tail right

			// Poisson constants, necessary for computing function values f(k)      
			l_my = Math.log(my);
			c_pm = m * l_my - Arithmetic.logFactorial(m);

			// function values f(k) = p(k)/p(m) at k = k2, k4, k1, k5               
			f2 = f(k2, l_my, c_pm);
			f4 = f(k4, l_my, c_pm);
			f1 = f(k1, l_my, c_pm);
			f5 = f(k5, l_my, c_pm);

			// area of the two centre and the two exponential tail regions          
			// area of the two immediate acceptance regions between k2, k4         
			p1 = f2 * (dl + 1.0);                    // immed. left    
			p2 = f2 * dl         + p1;               // centre left    
			p3 = f4 * (dr + 1.0) + p2;               // immed. right     
			p4 = f4 * dr         + p3;               // centre right     
			p5 = f1 / ll         + p4;               // expon. tail left 
			p6 = f5 / lr         + p5;               // expon. tail right
		} // end set-up

		for (;;) {
			// generate uniform number U -- U(0, p6)                                
			// case distinction corresponding to U                                  
			if ((U = gen.raw() * p6) < p2) {         // centre left      

				// immediate acceptance region R2 = [k2, m) *[0, f2),  X = k2, ... m -1 
				if ((V = U - p1) < 0.0)  return(k2 + (int)(U/f2));
				// immediate acceptance region R1 = [k1, k2)*[0, f1),  X = k1, ... k2-1 
				if ((W = V / dl) < f1 )  return(k1 + (int)(V/f1));

				// computation of candidate X < k2, and its counterpart Y > k2          
				// either squeeze-acceptance of X or acceptance-rejection of Y          
				Dk = (int)(dl * gen.raw()) + 1;
				if (W <= f2 - Dk * (f2 - f2/r2)) {            // quick accept of  
					return(k2 - Dk);                          // X = k2 - Dk      
				}
				if ((V = f2 + f2 - W) < 1.0) {                // quick reject of Y
					Y = k2 + Dk;
					if (V <= f2 + Dk * (1.0 - f2)/(dl + 1.0)) {// quick accept of  
						return(Y);                             // Y = k2 + Dk      
					}
					if (V <= f(Y, l_my, c_pm))  return(Y);    // final accept of Y
				}
				X = k2 - Dk;
			}
			else if (U < p4) {                                 // centre right     
				// immediate acceptance region R3 = [m, k4+1)*[0, f4), X = m, ... k4    
				if ((V = U - p3) < 0.0)  return(k4 - (int)((U - p2)/f4));
				// immediate acceptance region R4 = [k4+1, k5+1)*[0, f5)                
				if ((W = V / dr) < f5 )  return(k5 - (int)(V/f5));

				// computation of candidate X > k4, and its counterpart Y < k4          
				// either squeeze-acceptance of X or acceptance-rejection of Y          
				Dk = (int)(dr * gen.raw()) + 1;
				if (W <= f4 - Dk * (f4 - f4*r4)) {             // quick accept of  
					return(k4 + Dk);                           // X = k4 + Dk      
				}
				if ((V = f4 + f4 - W) < 1.0) {                 // quick reject of Y
					Y = k4 - Dk;
					if (V <= f4 + Dk * (1.0 - f4)/ dr) {       // quick accept of  
						return(Y);                             // Y = k4 - Dk      
					}
					if (V <= f(Y, l_my, c_pm))  return(Y);    // final accept of Y
				}
				X = k4 + Dk;
			}
			else {
				W = gen.raw();
				if (U < p5)	{                                  // expon. tail left 
					Dk = (int)(1.0 - Math.log(W)/ll);
					if ((X = k1 - Dk) < 0)  continue;          // 0 <= X <= k1 - 1 
					W *= (U - p4) * ll;                        // W -- U(0, h(x))  
					if (W <= f1 - Dk * (f1 - f1/r1))  return(X); // quick accept of X
				}
				else {                                         // expon. tail right
					Dk = (int)(1.0 - Math.log(W)/lr);
					X  = k5 + Dk;                              // X >= k5 + 1      
					W *= (U - p5) * lr;                        // W -- U(0, h(x))  
					if (W <= f5 - Dk * (f5 - f5*r5))  return(X); // quick accept of X
				}
			}

			// acceptance-rejection test of candidate X from the original area   
			// test, whether  W <= f(k),    with  W = U*h(x)  and  U -- U(0, 1)  
			// log f(X) = (X - m)*log(my) - log X! + log m!                       
			if (Math.log(W) <= X * l_my - Arithmetic.logFactorial(X) - c_pm)  return(X);	
		}
	}
	else { // mean is too large
		return (int) my;
	}
}
/**
 * Returns the probability distribution function.
 */
public double pdf(int k) {
	return Math.exp(k*Math.log(this.mean) - Arithmetic.logFactorial(k) - this.mean);
	
	// Overflow sensitive:
	// return (Math.pow(mean,k) / cephes.Arithmetic.factorial(k)) * Math.exp(-this.mean);
}
/**
 * Sets the mean.
 */
public void setMean(double mean) {
	this.mean = mean;
}
/**
 * Returns a random number from the distribution with the given mean.
 */
public static int staticNextInt(double mean) {
	synchronized (shared) {
		shared.setMean(mean);
		return shared.nextInt();
	}
}
/**
 * Returns a String representation of the receiver.
 */
public String toString() {
	return this.getClass().getName()+"("+mean+")";
}
/**
 * Sets the uniform random number generated shared by all <b>static</b> methods.
 * @param randomGenerator the new uniform random number generator to be shared.
 */
private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
	synchronized (shared) {
		shared.setRandomGenerator(randomGenerator);
	}
}
}
 public abstract class RandomEngine extends PersistentObject {
	/**
	 * Makes this class non instantiable, but still let's others inherit from
	 * it.
	 */
	protected RandomEngine() {
	}

	/**
	 * Equivalent to <tt>raw()</tt>.
	 * This has the effect that random engines can now be used as function
	 * objects, returning a random number upon function evaluation.
	 */
	public double apply(double dummy) {
		return raw();
	}

	/**
	 * Equivalent to <tt>nextInt()</tt>.
	 * This has the effect that random engines can now be used as function
	 * objects, returning a random number upon function evaluation.
	 */
	public int apply(int dummy) {
		return nextInt();
	}

	/**
	 * Constructs and returns a new uniform random number engine seeded with the
	 * current time.
	 * Currently this is {@link cern.jet.random.engine.MersenneTwister}.
	 */
	public static RandomEngine makeDefault() {
		return new MersenneTwisterFast(System.currentTimeMillis());
	}

	public static RandomEngine makeDefault(long seed) {
		return new MersenneTwisterFast(seed);
	}

	/**
	 * Returns a random double in the half-open range from [0.0,1.0). Thus 0.0
	 * is a valid result but 1.0 is not. <br>
	 * Note: this is different from the original Colt definition which
	 * excludes 0.0.
	 */
	public double nextDouble() {
		double nextDouble;

		do {
			// -9.223372036854776E18 == (double) Long.MIN_VALUE
			// 5.421010862427522E-20 == 1 / Math.pow(2,64) == 1 / ((double) Long.MAX_VALUE - (double) Long.MIN_VALUE);
			nextDouble = ((double) nextLong() - -9.223372036854776E18) * 5.421010862427522E-20;
		}
		// catch loss of precision of long --> double conversion
		while (!(nextDouble >= 0.0 && nextDouble < 1.0));

		// --> in (0.0,1.0)
		return nextDouble;

		/*
			nextLong == Long.MAX_VALUE         --> 1.0
			nextLong == Long.MIN_VALUE         --> 0.0
			nextLong == Long.MAX_VALUE-1       --> 1.0
			nextLong == Long.MAX_VALUE-100000L --> 0.9999999999999946
			nextLong == Long.MIN_VALUE+1       --> 0.0
			nextLong == Long.MIN_VALUE-100000L --> 0.9999999999999946
			nextLong == 1L                     --> 0.5
			nextLong == -1L                    --> 0.5
			nextLong == 2L                     --> 0.5
			nextLong == -2L                    --> 0.5
			nextLong == 2L+100000L             --> 0.5000000000000054
			nextLong == -2L-100000L            --> 0.49999999999999456
		*/
	}

	/**
	 * Returns a random float in the half-open range from [0.0f,1.0f). Thus 0.0f
	 * is a valid result but 1.0f is not. <br>
	 * Note: this is different from the original Colt definition which
	 * excludes 0.0.
	 */
	public float nextFloat() {
		// catch loss of precision of double --> float conversion
		float nextFloat;
		do {
			nextFloat = (float) ((nextInt() & 0xFFFFFFFFL) * 2.3283064365386963E-10);
		} while (nextFloat >= 1.0f);

		// --> in (0.0f,1.0f)
		return nextFloat;
	}

	/**
	 * Returns a 32 bit uniformly distributed random number in the closed
	 * interval <tt>[Integer.MIN_VALUE,Integer.MAX_VALUE]</tt> (including
	 * <tt>Integer.MIN_VALUE</tt> and <tt>Integer.MAX_VALUE</tt>);
	 */
	public abstract int nextInt();

	/**
	 * Returns an integer drawn uniformly from 0 to n-1. Suffice it to say,
	 * n must be > 0, or an IllegalArgumentException is raised.
	 */
	public int nextInt(int n) {
		return nextInt() % n;
	}

	/**
	 * Returns a 64 bit uniformly distributed random number in the closed
	 * interval <tt>[Long.MIN_VALUE,Long.MAX_VALUE]</tt> (including
	 * <tt>Long.MIN_VALUE</tt> and <tt>Long.MAX_VALUE</tt>).
	 */
	public long nextLong() {
		// concatenate two 32-bit strings into one 64-bit string
		return ((nextInt() & 0xFFFFFFFFL) << 32) | ((nextInt() & 0xFFFFFFFFL));
	}

	public abstract void setSeed(long seed);

	/**
	 * Returns a 32 bit uniformly distributed random number in the open unit
	 * interval <code>(0.0,1.0)</code> (excluding 0.0 and 1.0). raw() is used in
	 * many other Colt classes.
	 */
	public double raw() {
		int nextInt;
		do { // accept anything but zero
			nextInt = nextInt(); // in [Integer.MIN_VALUE,Integer.MAX_VALUE]-interval
		} while (nextInt == 0);

		// transform to (0.0,1.0)-interval
		// 2.3283064365386963E-10 == 1.0 / Math.pow(2,32)
		return (double) (nextInt & 0xFFFFFFFFL) * 2.3283064365386963E-10;

		/*
			nextInt == Integer.MAX_VALUE   --> 0.49999999976716936
			nextInt == Integer.MIN_VALUE   --> 0.5
			nextInt == Integer.MAX_VALUE-1 --> 0.4999999995343387
			nextInt == Integer.MIN_VALUE+1 --> 0.5000000002328306
			nextInt == 1                   --> 2.3283064365386963E-10
			nextInt == -1                  --> 0.9999999997671694
			nextInt == 2                   --> 4.6566128730773926E-10
			nextInt == -2                  --> 0.9999999995343387
		*/
	}
}
 public class JGnuplot extends Logger {
	public Integer id;
	public String style, style2d, style3d, header;
	public String plot2d, plot2dBar, plot3d, plotDensity, plotImage, multiplot, plotx;
	public String beforeStyle, beforeStyleVar, afterStyleVar, afterStyle, beforeHeader, extra, afterHeader, bp;
	public String terminal, output;
	public String sPatNumOut;
	protected String sFTemp;
	public String sFLoad;
	protected static JGnuplot tJGnuplot;
	protected static Integer maxId = 0;

	public static final int JG_InNewThread = 0x00000001;
	public static final int JG_DeleteTempFile = 0x00000002;
	public static final int JG_Pause = 0x00000004;
	public static final int JG_Log = 0x00000008;
	public static final int JG_Default = JG_DeleteTempFile | JG_Pause | JG_Log;
	public static final int JG_All = ~0x00000000;

	public JGnuplot(boolean bCopyLoadFromFile) {
		initialize(bCopyLoadFromFile);
	}

	public JGnuplot() {
		initialize(false);
	}

	public static synchronized JGnuplot tJGnuplot() {
		if (tJGnuplot == null) {
			tJGnuplot = new JGnuplot(true);
		}
		return tJGnuplot;
	}

	public static int delAllTempFiles() {
		return U.deleleFiles("JGnuplotTemp.*plt");
	}

	public synchronized static int nextId() {
		maxId++;
		return maxId;
	}

	public boolean initialize(boolean bCopyLoadFromFile) {
		boolean rtn = false;
		id = nextId();
		bp = "";
		sFLoad = "jgnuplot.xml";
		String sFTempEnd = "_" + id;
		sFTemp = "JGnuplotTemp" + sFTempEnd + ".plt";
		if (bCopyLoadFromFile) {
			rtn = copyLoadFromFile();
		} else {
			JGnuplot jg = tJGnuplot();
			style = jg.style;
			style2d = jg.style2d;
			style3d = jg.style3d;
			header = jg.header;
			extra = jg.extra;
			plot2d = jg.plot2d;
			plot2dBar = jg.plot2dBar;
			plot3d = jg.plot3d;
			plotDensity = jg.plotDensity;
			plotImage = jg.plotImage;
			multiplot = jg.multiplot;
			plotx = jg.plotx;
			sPatNumOut = jg.sPatNumOut;
			rtn = true;
		}
		initialize();
		return rtn;
	}

	public void initialize() {
	}

	public boolean load(String str) {
		return U.loadFromString(this, str);
	}

	public boolean copyLoadFromFile() {
		U.copyFileFromClassPath(JGnuplot.class, sFLoad, sFLoad, false);//this also works in a jar file
		return U.loadFromXML(this, sFLoad, false);
	}

	public void afterLoadObj(String sLoadMethod) {
		U.processFields(this, U.modAll, new Processable2.TrimStringField());
	}

	protected boolean prep() {
		return true;
	}

	public boolean compile(Plot plot, String code, String sFile) {
		boolean rtn = false;

		if (plot != null && sFile != null && prep()) {
			String sFileTU = U.toValidFileName(sFile);
			String codeEvaled = code;
			if (codeEvaled == null) {
				codeEvaled = plot2d;
				int nColumn = plot.nColumn(0, 0);
				if (nColumn == 3) {
					codeEvaled = plot3d;
				}
			}
			for (int i = 0; i < 5; i++) {
				codeEvaled = U.eval(codeEvaled, this, U.EVAL_InvalidIgnore | U.EVAL_NullIgnore);
			}
			rtn = plot.compile(codeEvaled, sFileTU, sPatNumOut);
		}

		return rtn;
	}

	public boolean compile(Plot plot, String code) {
		return compile(plot, code, plot.info + ".plt");
	}

	public boolean compile(Plot plot) {
		return compile(plot, null, plot.info + ".plt");
	}

	protected String sCommand(String sFile, boolean bPause) {
		String rtn = "gnuplot " + sFile;
		if (bPause) {
			rtn += " -e \"pause mouse key\"";
		}
		return rtn;
	}

	/**
	 * flags:<br>
	 * <b>JG_InNewThread</b> <br>
	 * <b>JG_DeleteTempFile:</b> use ~JG_DeleteTempFile flag to keep the temp
	 * file. This is useful when doing multi-thread plotting. The temp files
	 * shold be deleted at the end all together (delAllTempFiles()) to guarantee
	 * all figures can be plotted before temp files are deleted.<br>
	 * <b>JG_Pause</b><br>
	 * <b>JG_Log</b><br>
	 * <b>JG_Default</b> = JG_DeleteTempFile | JG_Pause | JG_Log<br>
	 * <b>JG_All</b>
	 * 
	 * @param plot
	 * @param code
	 * @param flags
	 * @return
	 */
	public Command execute(Plot plot, String code, int flags) {
		Command rtn = null;
		U.deleteFile(sFTemp);
		if (compile(plot, code, sFTemp)) {
			String sCmd = sCommand(sFTemp, (flags & JG_Pause) > 0);
			rtn = SysUtil.execCmd(sCmd, (flags & JG_Log) > 0, (flags & JG_InNewThread) > 0);
			if ((flags & JG_DeleteTempFile) > 0) {
				U.deleleFiles("^" + sFTemp + ".*");
			}
		}
		return rtn;
	}

	public Command execute(Plot plot, String code) {
		return execute(plot, code, JG_Default);
	}

	public Command execute(Plot plot) {
		return execute(plot, null, JG_Default);
	}

	/**
	 * Execute the plot in a new thread so that the current thread can go on
	 * without waiting. The execution results can be obtained through calling
	 * waitForRtn() of the return Command obj (this will cause the current
	 * thread to wait).
	 * 
	 * @param plot
	 * @param code
	 * @return
	 */
	public Command executeA(Plot plot, String code) {
		return execute(plot, code, JG_All);
	}

	public Command executeA(Plot plot) {
		return execute(plot, null, JG_All);
	}

	public static class Plot extends Logger {
		public String info;
		public String beforeStyle2, beforeStyleVar2, afterStyleVar2, afterStyle2, beforeHeader2, extra2, extra3, afterHeader2, bp2;
		public String xlabel, ylabel, zlabel, xrange, yrange, zrange;
		protected List<DataTableSet> lDataTableSet;
		protected String sPatNumOut;
		protected String sFile;
		protected Integer fid;

		protected Integer nextFid() {
			fid++;
			return fid;
		}

		public Plot(String info) {
			this.info = info;
			lDataTableSet = new ArrayList<DataTableSet>();
			initialize();
		}

		public void initialize() {
			fid = 0;
			bp2 = "";
		}

		public boolean autoSetLabels() {
			boolean rtn = false;
			DataTableSet dts = getDataTableSet(0);
			if (dts != null) {
				int ndt = dts.size();
				if (ndt > 0) {
					DataTable dt = dts.get(0);
					String[] cns = dt.getColNames();
					if (cns != null) {
						rtn = true;
						if (ndt == 1) {
							if (cns.length > 0) {//1 or more
								xlabel = cns[0];
							}
							if (cns.length > 1) {//2 or more
								ylabel = cns[1];
							}
							if (cns.length > 2) {//3 or more
								zlabel = cns[2];
							}
						} else {//ndt>1
							if (cns.length > 1) {//2 or more
								xlabel = cns[0];
							}
							if (cns.length > 2) {//3 or more
								ylabel = cns[1];
							}
						}

					}
				}
			}
			return rtn;
		}

		public boolean load(String str) {
			return U.loadFromString(this, str);
		}

		public DataTableSet addNewDataTableSet(String info) {
			DataTableSet dataTableSet = new DataTableSet(info);
			add(dataTableSet);
			return dataTableSet;
		}

		public int add(DataTableSet... dataTableSets) {
			int rtn = 0;
			for (int i = 0; i < dataTableSets.length; i++) {
				if (dataTableSets[i] != null) {
					if (lDataTableSet.add(dataTableSets[i])) {
						rtn++;
					}
				}
			}
			return rtn;
		}

		public boolean addAll(Collection<? extends DataTableSet> dataTableSets) {
			boolean rtn = false;
			if (dataTableSets != null) {
				rtn = lDataTableSet.addAll(dataTableSets);
			}
			return rtn;
		}

		public DataTableSet getDataTableSet(int i) {
			DataTableSet rtn = null;
			if (i >= 0 && i < lDataTableSet.size()) {
				return lDataTableSet.get(i);
			}
			return rtn;
		}

		public DataTable getDataTable(int iDataTableSet, int iDataTable) {
			DataTable rtn = null;
			DataTableSet dts = getDataTableSet(iDataTableSet);
			if (dts != null) {
				rtn = dts.get(iDataTable);
			}
			return rtn;
		}

		public List<DataTableSet> getDataTableSets() {
			return lDataTableSet;
		}

		protected boolean prep() {
			return true;
		}

		public boolean compile(String code, String sFile, String sPatNumOut) {
			boolean rtn = false;
			if (code != null && prep()) {
				this.sFile = sFile;
				this.sPatNumOut = sPatNumOut;
				//U.deleteFile(sFile);
				Scanner scanner = new Scanner(code);
				while (scanner.hasNextLine()) {
					String line = scanner.nextLine();
					try {
						//Those lines with invalid/null fields/methods/expressions will be cause a runtime exception here.
						String lineEvaled = line;
						for (int i = 0; i < 5 && (!U.bNoEval(lineEvaled)); i++) {
							lineEvaled = U.eval(lineEvaled, this, U.EVAL_InvalidException | U.EVAL_NullException);
						}
						if (U.bNoEval(lineEvaled)) {
							lineEvaled = lineEvaled.substring(U.sNoEval.length());
						}
						if (lineEvaled.length() > 0) {
							U.appendFile(sFile, lineEvaled + "\n");
						}
						rtn = true;
					} catch (TRuntimeException exception) {
						//log(exception);
					}
				}
				scanner.close();
			}
			return rtn;
		}

		public int size() {
			return lDataTableSet.size();
		}

		public int nColumn(int iDataTableSet, int iDataTable) {
			int rtn = -1;
			if (iDataTableSet >= 0 && iDataTableSet < lDataTableSet.size()) {
				DataTableSet dts = lDataTableSet.get(iDataTableSet);
				if (iDataTable >= 0 && iDataTable < dts.size()) {
					rtn = dts.get(iDataTable).nColumns();
				}
			}
			return rtn;
		}

		public Integer size(String dataTableSetNum) {
			Integer rtn = null;
			int iDataTableSet = Integer.parseInt(dataTableSetNum) - 1;
			if (iDataTableSet >= 0 && iDataTableSet < lDataTableSet.size()) {
				rtn = lDataTableSet.get(iDataTableSet).size();
			}
			return rtn;
		}

		public Integer nColumn(String dataTableSetNum, String dataTableNum) {
			int iDataTableSet = Integer.parseInt(dataTableSetNum) - 1;
			int iDataTable = Integer.parseInt(dataTableNum) - 1;
			Integer rtn = nColumn(iDataTableSet, iDataTable);
			if (rtn < 0) {
				rtn = null;
			}
			return rtn;
		}

		public String info(String dataTableSetNum) {
			String rtn = null;
			int iDataTableSet = Integer.parseInt(dataTableSetNum) - 1;
			if (iDataTableSet >= 0 && iDataTableSet < lDataTableSet.size()) {
				DataTableSet dts = lDataTableSet.get(iDataTableSet);
				rtn = dts.info;
			}
			return rtn;
		}

		public String info(String dataTableSetNum, String dataTableNum) {
			String rtn = null;
			int iDataTableSet = Integer.parseInt(dataTableSetNum) - 1;
			if (iDataTableSet >= 0 && iDataTableSet < lDataTableSet.size()) {
				DataTableSet dts = lDataTableSet.get(iDataTableSet);
				int iDataTable = Integer.parseInt(dataTableNum) - 1;
				if (iDataTable >= 0 && iDataTable < dts.size()) {
					rtn = dts.get(iDataTable).info;
				}
			}
			return rtn;
		}

		public String infos() {
			StringBuffer sBuffer = new StringBuffer();
			for (int i = 0, mi = lDataTableSet.size(); i < mi; i++) {
				DataTableSet dts = lDataTableSet.get(i);
				sBuffer.append("info" + (i + 1) + "=\"" + dts.info + "\";\n");
				for (int j = 0, mj = dts.size(); j < mj; j++) {
					DataTable dt = dts.get(j);
					sBuffer.append("info" + (i + 1) + "_" + (j + 1) + "=\"" + dt.info + "\";\n");
				}
			}
			String sInfos = sBuffer.toString();
			sInfos = U.sNoEval + sInfos.substring(0, sInfos.length() - 1);//remove the last line break;
			return sInfos;
		}

		protected void saveToFile(String sFile, DataTable dt, String type) {
			String sLastColumn1 = null;
			String[] colNames = dt.getColNames();
			try {
				FileWriter fw = new FileWriter(U.createIfNotExist(sFile), true);
				if (colNames != null) {
					String line = "#" + U.wrap(colNames, "", "", " ");
					fw.append(line + "\n");
				}
				for (int i = 0, mj = dt.nRows(); i < mj; i++) {
					Object[] row = dt.getRow(i);
					String[] sRow = U.toStrArray(Arrays.asList(row), sPatNumOut);
					if ("3d".equals(type)) {
						if (sLastColumn1 == null) {
							sLastColumn1 = sRow[0];
						}
						if (!sLastColumn1.equals(sRow[0])) {
							fw.append("\n");
							sLastColumn1 = sRow[0];
						}
					}
					String line = U.wrap(sRow, "", "", " ");
					fw.append(line + "\n");
				}
				fw.append("e\n");
				fw.flush();
				fw.close();
			} catch (IOException e) {
				log(e);
			}
		}

		public String data(String dataTableSetNum, String dataTableNum, String type) {
			String rtn = null;
			if (data2f(sFile, dataTableSetNum, dataTableNum, type) != null) {
				rtn = "";
			}
			return rtn;
		}

		public String data2f(String fn, String dataTableSetNum, String dataTableNum, String type) {
			String rtn = null;
			int iDataTableSet = Integer.parseInt(dataTableSetNum) - 1;
			if (fn != null && iDataTableSet >= 0 && iDataTableSet < lDataTableSet.size() && type != null) {
				DataTableSet dts = lDataTableSet.get(iDataTableSet);
				if (dataTableNum == null) {
					for (int i = 0, mi = dts.size(); i < mi; i++) {
						DataTable dt = dts.get(i);
						saveToFile(fn, dt, type);
					}
					rtn = fn;
				} else {
					int iDataTable = Integer.parseInt(dataTableNum) - 1;
					if (iDataTable >= 0 && iDataTable < dts.size()) {
						DataTable dt = dts.get(iDataTable);
						saveToFile(fn, dt, type);
						rtn = fn;
					}
				}
			}
			return rtn;
		}

		public String data2f(String dataTableSetNum, String dataTableNum, String type) {
			String rtn = data2f(sFile + "-D" + nextFid(), dataTableSetNum, dataTableNum, type);
			return rtn;
		}

		public String data2f(String dataTableSetNum, String type) {
			String rtn = data2f(dataTableSetNum, null, type);
			return rtn;
		}

		public String data(String dataTableSetNum, String type) {
			return data(dataTableSetNum, null, type);
		}

		public String data(String dataTableSetNum) {
			return data(dataTableSetNum, null, "2d");
		}

		public String setInfo(String dataTableSetNum, String info) {
			int iDataTableSet = Integer.parseInt(dataTableSetNum) - 1;
			if (iDataTableSet >= 0 && iDataTableSet < lDataTableSet.size()) {
				DataTableSet dts = lDataTableSet.get(iDataTableSet);
				dts.info = info;
			}
			return "";
		}
	}

}
 public static class Plot extends Logger {
	public String info;
	public String beforeStyle2, beforeStyleVar2, afterStyleVar2, afterStyle2, beforeHeader2, extra2, extra3, afterHeader2, bp2;
	public String xlabel, ylabel, zlabel, xrange, yrange, zrange;
	protected List<DataTableSet> lDataTableSet;
	protected String sPatNumOut;
	protected String sFile;
	protected Integer fid;

	protected Integer nextFid() {
		fid++;
		return fid;
	}

	public Plot(String info) {
		this.info = info;
		lDataTableSet = new ArrayList<DataTableSet>();
		initialize();
	}

	public void initialize() {
		fid = 0;
		bp2 = "";
	}

	public boolean autoSetLabels() {
		boolean rtn = false;
		DataTableSet dts = getDataTableSet(0);
		if (dts != null) {
			int ndt = dts.size();
			if (ndt > 0) {
				DataTable dt = dts.get(0);
				String[] cns = dt.getColNames();
				if (cns != null) {
					rtn = true;
					if (ndt == 1) {
						if (cns.length > 0) {//1 or more
							xlabel = cns[0];
						}
						if (cns.length > 1) {//2 or more
							ylabel = cns[1];
						}
						if (cns.length > 2) {//3 or more
							zlabel = cns[2];
						}
					} else {//ndt>1
						if (cns.length > 1) {//2 or more
							xlabel = cns[0];
						}
						if (cns.length > 2) {//3 or more
							ylabel = cns[1];
						}
					}

				}
			}
		}
		return rtn;
	}

	public boolean load(String str) {
		return U.loadFromString(this, str);
	}

	public DataTableSet addNewDataTableSet(String info) {
		DataTableSet dataTableSet = new DataTableSet(info);
		add(dataTableSet);
		return dataTableSet;
	}

	public int add(DataTableSet... dataTableSets) {
		int rtn = 0;
		for (int i = 0; i < dataTableSets.length; i++) {
			if (dataTableSets[i] != null) {
				if (lDataTableSet.add(dataTableSets[i])) {
					rtn++;
				}
			}
		}
		return rtn;
	}

	public boolean addAll(Collection<? extends DataTableSet> dataTableSets) {
		boolean rtn = false;
		if (dataTableSets != null) {
			rtn = lDataTableSet.addAll(dataTableSets);
		}
		return rtn;
	}

	public DataTableSet getDataTableSet(int i) {
		DataTableSet rtn = null;
		if (i >= 0 && i < lDataTableSet.size()) {
			return lDataTableSet.get(i);
		}
		return rtn;
	}

	public DataTable getDataTable(int iDataTableSet, int iDataTable) {
		DataTable rtn = null;
		DataTableSet dts = getDataTableSet(iDataTableSet);
		if (dts != null) {
			rtn = dts.get(iDataTable);
		}
		return rtn;
	}

	public List<DataTableSet> getDataTableSets() {
		return lDataTableSet;
	}

	protected boolean prep() {
		return true;
	}

	public boolean compile(String code, String sFile, String sPatNumOut) {
		boolean rtn = false;
		if (code != null && prep()) {
			this.sFile = sFile;
			this.sPatNumOut = sPatNumOut;
			//U.deleteFile(sFile);
			Scanner scanner = new Scanner(code);
			while (scanner.hasNextLine()) {
				String line = scanner.nextLine();
				try {
					//Those lines with invalid/null fields/methods/expressions will be cause a runtime exception here.
					String lineEvaled = line;
					for (int i = 0; i < 5 && (!U.bNoEval(lineEvaled)); i++) {
						lineEvaled = U.eval(lineEvaled, this, U.EVAL_InvalidException | U.EVAL_NullException);
					}
					if (U.bNoEval(lineEvaled)) {
						lineEvaled = lineEvaled.substring(U.sNoEval.length());
					}
					if (lineEvaled.length() > 0) {
						U.appendFile(sFile, lineEvaled + "\n");
					}
					rtn = true;
				} catch (TRuntimeException exception) {
					//log(exception);
				}
			}
			scanner.close();
		}
		return rtn;
	}

	public int size() {
		return lDataTableSet.size();
	}

	public int nColumn(int iDataTableSet, int iDataTable) {
		int rtn = -1;
		if (iDataTableSet >= 0 && iDataTableSet < lDataTableSet.size()) {
			DataTableSet dts = lDataTableSet.get(iDataTableSet);
			if (iDataTable >= 0 && iDataTable < dts.size()) {
				rtn = dts.get(iDataTable).nColumns();
			}
		}
		return rtn;
	}

	public Integer size(String dataTableSetNum) {
		Integer rtn = null;
		int iDataTableSet = Integer.parseInt(dataTableSetNum) - 1;
		if (iDataTableSet >= 0 && iDataTableSet < lDataTableSet.size()) {
			rtn = lDataTableSet.get(iDataTableSet).size();
		}
		return rtn;
	}

	public Integer nColumn(String dataTableSetNum, String dataTableNum) {
		int iDataTableSet = Integer.parseInt(dataTableSetNum) - 1;
		int iDataTable = Integer.parseInt(dataTableNum) - 1;
		Integer rtn = nColumn(iDataTableSet, iDataTable);
		if (rtn < 0) {
			rtn = null;
		}
		return rtn;
	}

	public String info(String dataTableSetNum) {
		String rtn = null;
		int iDataTableSet = Integer.parseInt(dataTableSetNum) - 1;
		if (iDataTableSet >= 0 && iDataTableSet < lDataTableSet.size()) {
			DataTableSet dts = lDataTableSet.get(iDataTableSet);
			rtn = dts.info;
		}
		return rtn;
	}

	public String info(String dataTableSetNum, String dataTableNum) {
		String rtn = null;
		int iDataTableSet = Integer.parseInt(dataTableSetNum) - 1;
		if (iDataTableSet >= 0 && iDataTableSet < lDataTableSet.size()) {
			DataTableSet dts = lDataTableSet.get(iDataTableSet);
			int iDataTable = Integer.parseInt(dataTableNum) - 1;
			if (iDataTable >= 0 && iDataTable < dts.size()) {
				rtn = dts.get(iDataTable).info;
			}
		}
		return rtn;
	}

	public String infos() {
		StringBuffer sBuffer = new StringBuffer();
		for (int i = 0, mi = lDataTableSet.size(); i < mi; i++) {
			DataTableSet dts = lDataTableSet.get(i);
			sBuffer.append("info" + (i + 1) + "=\"" + dts.info + "\";\n");
			for (int j = 0, mj = dts.size(); j < mj; j++) {
				DataTable dt = dts.get(j);
				sBuffer.append("info" + (i + 1) + "_" + (j + 1) + "=\"" + dt.info + "\";\n");
			}
		}
		String sInfos = sBuffer.toString();
		sInfos = U.sNoEval + sInfos.substring(0, sInfos.length() - 1);//remove the last line break;
		return sInfos;
	}

	protected void saveToFile(String sFile, DataTable dt, String type) {
		String sLastColumn1 = null;
		String[] colNames = dt.getColNames();
		try {
			FileWriter fw = new FileWriter(U.createIfNotExist(sFile), true);
			if (colNames != null) {
				String line = "#" + U.wrap(colNames, "", "", " ");
				fw.append(line + "\n");
			}
			for (int i = 0, mj = dt.nRows(); i < mj; i++) {
				Object[] row = dt.getRow(i);
				String[] sRow = U.toStrArray(Arrays.asList(row), sPatNumOut);
				if ("3d".equals(type)) {
					if (sLastColumn1 == null) {
						sLastColumn1 = sRow[0];
					}
					if (!sLastColumn1.equals(sRow[0])) {
						fw.append("\n");
						sLastColumn1 = sRow[0];
					}
				}
				String line = U.wrap(sRow, "", "", " ");
				fw.append(line + "\n");
			}
			fw.append("e\n");
			fw.flush();
			fw.close();
		} catch (IOException e) {
			log(e);
		}
	}

	public String data(String dataTableSetNum, String dataTableNum, String type) {
		String rtn = null;
		if (data2f(sFile, dataTableSetNum, dataTableNum, type) != null) {
			rtn = "";
		}
		return rtn;
	}

	public String data2f(String fn, String dataTableSetNum, String dataTableNum, String type) {
		String rtn = null;
		int iDataTableSet = Integer.parseInt(dataTableSetNum) - 1;
		if (fn != null && iDataTableSet >= 0 && iDataTableSet < lDataTableSet.size() && type != null) {
			DataTableSet dts = lDataTableSet.get(iDataTableSet);
			if (dataTableNum == null) {
				for (int i = 0, mi = dts.size(); i < mi; i++) {
					DataTable dt = dts.get(i);
					saveToFile(fn, dt, type);
				}
				rtn = fn;
			} else {
				int iDataTable = Integer.parseInt(dataTableNum) - 1;
				if (iDataTable >= 0 && iDataTable < dts.size()) {
					DataTable dt = dts.get(iDataTable);
					saveToFile(fn, dt, type);
					rtn = fn;
				}
			}
		}
		return rtn;
	}

	public String data2f(String dataTableSetNum, String dataTableNum, String type) {
		String rtn = data2f(sFile + "-D" + nextFid(), dataTableSetNum, dataTableNum, type);
		return rtn;
	}

	public String data2f(String dataTableSetNum, String type) {
		String rtn = data2f(dataTableSetNum, null, type);
		return rtn;
	}

	public String data(String dataTableSetNum, String type) {
		return data(dataTableSetNum, null, type);
	}

	public String data(String dataTableSetNum) {
		return data(dataTableSetNum, null, "2d");
	}

	public String setInfo(String dataTableSetNum, String info) {
		int iDataTableSet = Integer.parseInt(dataTableSetNum) - 1;
		if (iDataTableSet >= 0 && iDataTableSet < lDataTableSet.size()) {
			DataTableSet dts = lDataTableSet.get(iDataTableSet);
			dts.info = info;
		}
		return "";
	}
}
 public class TestUtil extends LogUtil {
	public static double[] compareRunTime(int nRept, int nComp, Runnable... runs) {
		double[] rtn = new double[runs.length];
		Timer timer = new Timer();
		tLog("TestUtil.compareRunTime : ", U.toStr(runs));
		for (int i = 0; i < nComp; i++) {
			tLog("----- " + i + " -----");
			for (int k = 0; k < runs.length; k++) {
				long tTotal = 0;
				for (int j = 0; j < nRept; j++) {
					timer.start();
					runs[k].run();
					tTotal += timer.stop();
				}
				double tAverage = (double) tTotal / nRept;
				rtn[k] = tAverage;
				tLog("r" + k + " [average, total]ms : ", tAverage, tTotal);
			}
		}
		tLog("----- End -----");
		return rtn;
	}

	public static double[] compareRunTime(int nRept, int nComp, Processable0<Double>... pros) {
		double[] rtn = new double[pros.length];
		Timer timer = new Timer();
		long tDifference = 0;
		Double dummyRtn = 0.0;//do dummy operation on p1 and p2's return results to avoid java's optimisation from cancelling each run.

		tLog("TestUtil.compareRunTime : ", U.toStr(pros));
		for (int i = 0; i < nComp; i++) {
			tLog("----- " + i + " -----");
			for (int k = 0; k < pros.length; k++) {
				long tTotal = 0;
				for (int j = 0; j < nRept; j++) {
					timer.start();
					dummyRtn += pros[k].process();
					tTotal += timer.stop();
				}
				double tAverage = (double) tTotal / nRept;
				rtn[k] = tAverage;
				tLog("p" + k + " [average, total]ms : ", tAverage, tTotal);
			}
		}
		tLog("----- End -----");
		tLog("dummySum : " + dummyRtn);
		return rtn;
	}
}
 public class U extends SerialUtil {

}
 public interface Processable0<R> {
	public R process();
}
 public interface Processable2<R, A, B> {
	public R process(A a, B b);

	public static class TrimStringField implements Processable2<Boolean, Object, Field> {

		public Boolean process(Object obj, Field field) {
			boolean rtn = false;
			if (field != null && U.bAssignable(String.class, field.getType())) {
				String sValue = (String) U.getFieldValue(obj, field);
				if (sValue != null) {
					String sValue2 = sValue.trim();
					if (sValue2.length() < sValue.length()) {
						rtn = true;
					}
					U.setFieldValue(obj, field, sValue2);
				}
			}
			return rtn;
		}

	}

	public static class SampleListByIndexB<B> implements Processable2<Boolean, Integer, B> {
		public Integer step;

		/**
		 * The sample results will with these index: 0, 0+step, 0+2*step, ...
		 * 
		 * @param step
		 */
		public SampleListByIndexB(Integer step) {
			this.step = step;
		}

		public Boolean process(Integer i, B b) {
			return i % step == 0;
		}

	}

	public static class SampleListByIndex<B> implements Processable2<B, Integer, B> {
		public Integer step;

		/**
		 * The sample results will with these index: 0, 0+step, 0+2*step, ...
		 * 
		 * @param step
		 */
		public SampleListByIndex(Integer step) {
			this.step = step;
		}

		public B process(Integer i, B b) {
			B rtn = null;
			if (i % step == 0) {
				rtn = b;
			}
			return rtn;
		}
	}

	public static class RemoveRowWithInvalidNumber<B> implements Processable2<B[], Integer, B[]> {
		public String sInvalid;

		/**
		 * E.g. "$%$<0" : remove rows with negative number.
		 * 
		 * @param sInvalid
		 */
		public RemoveRowWithInvalidNumber(String sInvalid) {
			this.sInvalid = sInvalid;
		}

		public B[] process(Integer a, B[] b) {
			B[] rtn = b;
			for (int i = 0; i < b.length; i++) {
				Double d = U.toDouble(b[i] + "");
				if (d != null) {
					if (U.evalCheck(sInvalid, d)) {
						rtn = null;
						break;
					}
				}
			}
			return rtn;
		}
	};

}
 public class DataTable implements Cloneable {
	public String info;
	protected String[] colNames;
	protected List<Object[]> rows;
	private int nColumns;

	public DataTable(String info, int nColumns) {
		this.info = info;
		this.nColumns = nColumns;
		rows = new ArrayList<Object[]>();
	}

	public String[] getColNames() {
		return colNames;
	}

	public boolean setColNames(String[] colNames) {
		boolean rtn = false;
		if (colNames != null && colNames.length == nColumns) {
			rtn = true;
			this.colNames = colNames;
		}
		return rtn;
	}

	public <R, B> boolean addAll(boolean bAddUnEqualSizedList, Processable2<R, Integer, B> pa2Element, List<B>... lists) {
		boolean rtn = false;
		List[] tLists = lists;
		if (tLists != null && tLists.length == nColumns) {
			rtn = true;
			if (pa2Element != null) {
				for (int i = 0; i < nColumns; i++) {
					tLists[i] = U.processList(lists[i], pa2Element);
				}
			}
			int size = tLists[0].size();
			if (!U.bEqualSize(tLists)) {
				if (bAddUnEqualSizedList) {
					size = U.minSize(tLists);
				} else {
					rtn = false;
					throw new TRuntimeException("addAll: lists are not of the same size for DataTable: " + info);
				}
			}
			for (int i = 0; i < size; i++) {
				Object[] row = new Object[nColumns];
				for (int j = 0; j < nColumns; j++) {
					row[j] = tLists[j].get(i);
				}
				add(row);
			}
		}
		return rtn;
	}

	public <R> boolean addAll(boolean bAddUnEqualSizedList, Processable2<R, Integer, Double> pa2Element, double[]... arrays) {
		boolean rtn = false;
		if (arrays != null && arrays.length == nColumns) {
			List<Double>[] lists = new List[nColumns];
			for (int i = 0; i < arrays.length; i++) {
				lists[i] = U.parseList(arrays[i]);
			}
			rtn = addAll(bAddUnEqualSizedList, pa2Element, lists);
		}
		return rtn;
	}

	public <R, B> boolean addAll(boolean bAddUnEqualSizedList, Processable2<R, Integer, B[]> pa2Element, List<B[]> list) {
		boolean rtn = false;
		if (list != null && list.size() > 0) {
			B[] b0 = list.get(0);
			if (b0 != null && b0.length == nColumns) {
				List tList = list;
				if (pa2Element != null) {
					tList = U.processList(list, pa2Element);
				}
				rtn = rows.addAll(tList);
			}
		}
		return rtn;
	}
	
	public boolean insert(int index, Object... objs) {
		boolean rtn = false;
		if (objs != null && objs.length == nColumns) {
			rtn = true;
			rows.add(index, objs);
		}
		return rtn;
	}

	public boolean add(Object... objs) {
		boolean rtn = false;
		if (objs != null && objs.length == nColumns) {
			rtn = true;
			rows.add(objs);
		}
		return rtn;
	}

	public int nColumns() {
		return nColumns;
	}

	public int nRows() {
		return rows.size();
	}

	public Object[] getRow(int i) {
		return rows.get(i);
	}

	public Object[] setRow(int i, Object[] row) {
		return rows.set(i, row);
	}

	public Object[] removeRow(int i) {
		return rows.remove(i);
	}

	public List<Object[]> getRows() {
		return rows;
	}

	public List<Object> getColumn(int i) {
		List<Object> rtn = null;
		if (i >= 0 && i < nColumns) {
			rtn = new ArrayList<Object>();
			for (int j = 0, mj = rows.size(); j < mj; j++) {
				Object[] row = rows.get(j);
				rtn.add(row[i]);
			}
		}
		return rtn;
	}

	/**
	 * Can be used t crate a sub DataTable with rows from iRowFrom to iRowTo and
	 * columns with iCols. There could be repeated ids in iCols and iCols.length
	 * could be bigger than nColumns so as to generate a DataTable with some
	 * repeated columns when necessary.
	 * 
	 * @param iRowFrom
	 * @param iRowTo
	 * @param iCols
	 * @return
	 */
	public DataTable subTable(int iRowFrom, int iRowTo, int... iCols) {
		DataTable rtn = null;
		if (iRowFrom >= 0 && iRowTo >= iRowFrom && iRowTo < rows.size()) {
			int[] iColsTU = iCols;
			if (iColsTU == null) {//null means all columns in the original order.
				iColsTU = new int[nColumns];
				for (int i = 0; i < iColsTU.length; i++) {
					iColsTU[i] = i;
				}
			}
			rtn = new DataTable(info, iColsTU.length);
			if (colNames != null) {
				String[] subColNames = new String[iColsTU.length];
				for (int i = 0; i < iColsTU.length; i++) {
					subColNames[i] = colNames[iColsTU[i]];
				}
				rtn.setColNames(subColNames);
			}
			for (int i = iRowFrom; i <= iRowTo; i++) {
				Object[] row = rows.get(i);
				Object[] subRow = new Object[iColsTU.length];
				for (int j = 0; j < iColsTU.length; j++) {
					subRow[j] = row[iColsTU[j]];
				}
				rtn.add(subRow);
			}
		}
		return rtn;
	}

	public DataTable subRows(int iRowFrom, int iRowTo) {
		return subTable(iRowFrom, iRowTo, null);
	}

	public DataTable subColumns(int... iCols) {
		return subTable(0, rows.size() - 1, iCols);
	}

	@SuppressWarnings("unchecked")
	public <B> int processRows(Processable2<B[], Integer, B[]> pa2Row) {
		int rtn = -1;
		if (pa2Row != null) {
			int nRows = rows.size();
			Processable2 pa2 = pa2Row;
			rows = U.processList(rows, pa2);
			rtn = nRows - rows.size();
		}
		return rtn;
	}

	public void sort(Comparator<? super Object[]> c) {
		Collections.sort(rows, c);
	}

	public void sort(boolean bAscending) {
		sort(new ObjArrayComparator<Object>(bAscending));
	}

	public void sort() {
		sort(true);
	}

	public DataTable clone() {
		return subTable(0, nRows() - 1, null);
	}
}
 public class DataTableSet implements Cloneable {
	public String info;
	public List<DataTable> members;

	public DataTableSet(String info) {
		this.info = info;
		members = new ArrayList<DataTable>();
	}

	public DataTableSet() {
		this(null);
	}

	public boolean add(DataTable dataTable) {
		boolean rtn = false;
		if (dataTable != null) {
			rtn = true;
			members.add(dataTable);
		}
		return rtn;
	}

	public boolean add(DataTable... dataTables) {
		boolean rtn = false;
		if (dataTables != null && dataTables.length > 0) {
			rtn = true;
			for (int i = 0; i < dataTables.length; i++) {
				if (!add(dataTables[i])) {
					rtn = false;
					break;
				}
			}
		}
		return rtn;
	}

	public DataTable addNewDataTable(String info, int nColumn) {
		DataTable rtn = new DataTable(info, nColumn);
		add(rtn);
		return rtn;
	}

	public <R, B> DataTable addNewDataTable(String info, boolean bAddUnEqualSizedList, Processable2<R, Integer, B> pa2Element, List<B>... lists) {
		DataTable rtn = new DataTable(info, lists.length);
		if (rtn.addAll(bAddUnEqualSizedList, pa2Element, lists)) {
			add(rtn);
		} else {
			rtn = null;
		}
		return rtn;
	}

	public <R> DataTable addNewDataTable(String info, boolean bAddUnEqualSizedList, Processable2<R, Integer, Double> pa2Element, double[]... arrays) {
		DataTable rtn = new DataTable(info, arrays.length);
		if (rtn.addAll(bAddUnEqualSizedList, pa2Element, arrays)) {
			add(rtn);
		} else {
			rtn = null;
		}
		return rtn;
	}

	public <R, B> DataTable addNewDataTable(String info, boolean bAddUnEqualSizedList, Processable2<R, Integer, B[]> pa2Element, List<B[]> list) {
		DataTable rtn = null;
		if (list != null && list.size() > 0) {
			B[] b0 = list.get(0);
			if (b0 != null) {
				rtn = new DataTable(info, b0.length);
				if (rtn.addAll(bAddUnEqualSizedList, pa2Element, list)) {
					add(rtn);
				} else {
					rtn = null;
				}
			}
		}
		return rtn;
	}

	public <R, B> DataTable addNewDataTable(String info, Processable2<R, Integer, B> pa2Element, List<B>... lists) {
		return addNewDataTable(info, false, pa2Element, lists);
	}

	public <B> DataTable addNewDataTable(String info, List<B>... lists) {
		return addNewDataTable(info, false, null, lists);
	}

	public <R> DataTable addNewDataTable(String info, Processable2<R, Integer, Double> pa2Element, double[]... arrays) {
		return addNewDataTable(info, false, pa2Element, arrays);
	}

	public DataTable addNewDataTable(String info, double[]... arrays) {
		return addNewDataTable(info, false, null, arrays);
	}

	public <R, B> DataTable addNewDataTable(String info, Processable2<R, Integer, B[]> pa2Element, List<B[]> list) {
		return addNewDataTable(info, false, pa2Element, list);
	}

	public <B> DataTable addNewDataTable(String info, List<B[]> list) {
		return addNewDataTable(info, false, null, list);
	}

	public int size() {
		return members.size();
	}

	/**
	 * 
	 * @param i
	 *            starts from 0
	 * @return
	 */
	public DataTable get(int i) {
		DataTable rtn = null;
		if (i >= 0 && i < members.size()) {
			rtn = members.get(i);
		}
		return rtn;
	}

	public boolean set(int i, DataTable dt) {
		boolean rtn = false;
		if (i >= 0 && i < members.size()) {
			rtn = true;
			members.set(i, dt);
		}
		return rtn;
	}

	public DataTable remove(int i) {
		return members.remove(i);
	}

	public DataTableSet clone() {
		DataTableSet rtn = new DataTableSet(info);
		for (int i = 0, mi = members.size(); i < mi; i++) {
			DataTable dt = members.get(i);
			rtn.add(dt.clone());
		}
		return rtn;
	}
}
 public class Distribution extends Logger implements Serializable {
	private static final long serialVersionUID = 3260826314035959813L;
	protected static double dDiscrete = 1E-10;
	public Integer nTotal;
	public Integer nCounted;
	public Double dMin;
	public Double dMax;
	public Double dInterval;
	public List<Double> data;
	protected List<DistElement> distElements;
	public String sPatNumOut;
	public Integer maxListOutSize; //the maximum output size of lists. default to be null: no limit.
	public String info;

	public static class DistElement implements Comparable, Serializable {
		private static final long serialVersionUID = 59637195809162373L;
		public Double dFromInclude;
		public Double dTo;
		public Integer n;
		public Double p;
		public List<Double> data;

		public int compareTo(Object o) {
			DistElement de = (DistElement) o;
			return dFromInclude.compareTo(de.dFromInclude);
		}
	}

	public Distribution(String sInfo) {
		nTotal = 0;
		nCounted = 0;
		dInterval = null;
		dMax = null;
		dMin = null;
		data = new ArrayList<Double>();
		distElements = new ArrayList<DistElement>();
		sPatNumOut = null;
		maxListOutSize = null;
		info = sInfo;
	}

	public Distribution() {
		this(null);
	}

	/**
	 * dInterval should not be smaller than or equal to 0!
	 * 
	 * @param dInterval
	 * @param dFrom
	 * @param dTo
	 * @return
	 */

	public boolean prepDistElements(Double dInterval, Double dFrom, Double dTo) {
		boolean rtn = false;
		if (dTo >= dFrom && dInterval > 0) {
			rtn = true;
			distElements.clear();
			this.dInterval = dInterval;
			int nElements = (int) ((dTo - dFrom) / dInterval + 1);
			((ArrayList) distElements).ensureCapacity(nElements);
			BigDecimal bd = new BigDecimal(dFrom);
			BigDecimal bdInterval = new BigDecimal(dInterval);
			BigDecimal bdTo = new BigDecimal(dTo);
			for (; bd.compareTo(bdTo) <= 0; bd = bd.add(bdInterval)) {
				DistElement de = new DistElement();
				de.dFromInclude = bd.doubleValue();
				de.dTo = bd.add(bdInterval).doubleValue();
				de.n = 0;
				de.p = 0.0;
				de.data = new ArrayList<Double>();
				distElements.add(de);
			}
		} else {
			log("Error dTo<dFrom! [dTo, dFrom]: ", dTo, dFrom);
		}
		return rtn;
	}

	public boolean prep(Double dInterval, Double dFrom, Double dTo) {
		nTotal = 0;
		nCounted = 0;
		this.dInterval = dInterval;
		this.dMin = dFrom;
		this.dMax = dTo;
		data = new ArrayList<Double>();
		distElements = new ArrayList<DistElement>();
		return prepDistElements(dInterval, dFrom, dTo);
	}

	public boolean put(Double dVal) {
		boolean rtn = false;
		if (dVal != null) {
			if (dMax == null || dVal > dMax) {
				dMax = dVal;
			}
			if (dMin == null || dVal < dMin) {
				dMin = dVal;
			}
			data.add(dVal);
			rtn = true;
			nTotal++;
		}
		return rtn;
	}

	public boolean statDiscrete(Double dVal) {
		boolean rtn = false;
		if (dVal != null) {
			rtn = true;
			boolean found = false;
			for (int i = 0, mi = distElements.size(); i < mi; i++) {
				DistElement de = distElements.get(i);
				if (dVal >= de.dFromInclude && dVal <= de.dTo) {
					de.n++;
					nCounted++;
					found = true;
					break;
				}
			}
			if (!found) {
				DistElement de = new DistElement();
				de.dFromInclude = dVal - dDiscrete;
				de.dTo = dVal + dDiscrete;
				de.n = 0;
				de.p = 0.0;
				distElements.add(de);
			}
		}
		return rtn;
	}

	/**
	 * This function does not put the columns into the columns list, just do
	 * statistics. The construction function Distribution(Double dInterval,
	 * Double dFrom, Double dTo) has to be used before using this function. This
	 * function will save lots of memory compared to the function of put.
	 * 
	 * Note: only use calDist() to get the distribution results after using
	 * stat!
	 * 
	 * @return
	 */
	public boolean stat(Double dVal) {
		boolean rtn = false;

		if (dVal != null) {
			if (dVal >= dMin && dVal <= dMax) {
				int id = (int) ((dVal - dMin) / dInterval);
				DistElement de = distElements.get(id);
				de.n++;
				nCounted++;
				//de.p = (double) de.n / nCounted;
			}
			nTotal++;
			rtn = true;
		}

		return rtn;
	}

	/**
	 * Only use this function after using stat(Double dVal) or
	 * statDiscrete(Double dVal).
	 * 
	 * @return
	 */

	public List<DistElement> calDist() {
		List<DistElement> rtn = distElements;
		Collections.sort(distElements);

		for (int i = 0, size = distElements.size(); i < size; i++) {
			DistElement de = distElements.get(i);
			de.p = (double) de.n / nCounted;
		}

		return rtn;
	}

	public int clear() {
		int rtn = data.size();
		data.clear();
		distElements.clear();
		nCounted = 0;
		nTotal = 0;
		return rtn;
	}

	public int getNTotal() {
		return nTotal;
	}

	public int getNCounted() {
		return nCounted;
	}

	public List<DistElement> calDist(Double dInterval, Double dFrom, Double dTo) {
		List<DistElement> rtn = null;

		if (dInterval != null && dInterval > 0 && nTotal > 0) {
			if (dFrom == null) {
				dFrom = dMin;
			}
			if (dTo == null) {
				dTo = dMax;
			}
			if (prepDistElements(dInterval, dFrom, dTo)) {
				nCounted = 0;
				for (int i = 0, size = data.size(); i < size; i++) {
					Double d = data.get(i);
					if (d >= dFrom && d <= dTo) {
						int id = (int) ((d - dFrom) / dInterval);
						DistElement de = distElements.get(id);
						de.n++;
						de.data.add(d);
						nCounted++;
					}
				}
				for (int i = 0, size = distElements.size(); i < size; i++) {
					DistElement de = distElements.get(i);
					de.p = (double) de.n / nCounted;
				}
				rtn = distElements;
			}
		} else {
			log("Wrong dInterval value: " + dInterval);
		}

		return rtn;
	}

	public List<DistElement> calDist(Double dInterval, Double dFrom) {
		return calDist(dInterval, dFrom, null);
	}

	public List<DistElement> calDist(Double dInterval) {
		return calDist(dInterval, null, null);
	}

	/**
	 * Automatically find a proper interval to separate the range of columns
	 * values
	 * into n parts
	 * 
	 * @param n
	 * @return
	 */
	public List<DistElement> calDistIntoNparts(Integer n) {
		List<DistElement> rtn = null;
		if (n >= 1 && nTotal > 0) {
			Double dInterval = (dMax - dMin) / n;
			if (dInterval == 0) {
				dInterval = 1.0 / 10000;
			} else {
				dInterval += dInterval / 10000;
			}
			rtn = calDist(dInterval);
		}
		return rtn;
	}

	public DataTable getDataTable(DataTable dataTable, int... iCols) {
		DataTable rtn = dataTable;
		if (distElements.size() > 0) {
			int[] iColsTU = iCols;
			if (iColsTU == null) {
				iColsTU = new int[] { 0, 1, 2, 3 };
			}
			if (rtn == null) {
				rtn = new DataTable(info, iColsTU.length);
			}
			String[] colNames = new String[iColsTU.length];
			String[] sDistCols = new String[] { "dFromInclude", "dTo", "n", "p" };
			for (int i = 0; i < iColsTU.length; i++) {
				colNames[i] = sDistCols[iColsTU[i]];
			}
			rtn.setColNames(colNames);
			for (int i = 0, mi = distElements.size(); i < mi; i++) {
				DistElement dei = distElements.get(i);
				Object[] oDistCols = new Object[] { dei.dFromInclude, dei.dTo, dei.n, dei.p };
				Object[] row = new Object[iColsTU.length];
				for (int j = 0; j < iColsTU.length; j++) {
					row[j] = oDistCols[iColsTU[j]];
				}
				rtn.add(row);
			}
		}
		return rtn;
	}

	public DataTable getDataTable() {
		return getDataTable(null, 0, 3);
	}

	protected void output(FileWriter fw, String str) {
		if (!U.append(fw, str)) {
			log(false, str);
		}
	}

	@SuppressWarnings("all")
	public void output(String fn, String before, String after, String delimiter) {
		Processable1 pa1 = new Processable1.MaxListOutSize(maxListOutSize);
		boolean bOutClass = true;
		String outString = U.toOutString(this, before, after, delimiter, sPatNumOut, bOutClass, pa1);
		FileWriter fw = U.createFileWriter(fn);
		output(fw, outString);
		for (int i = 0, mi = distElements.size(); i < mi; i++) {
			DistElement dei = distElements.get(i);
			outString = U.toOutString(dei, before, after, delimiter, sPatNumOut, bOutClass, pa1);
			output(fw, outString);
			bOutClass = false;
		}
		output(fw, U.wrap(before, after, delimiter, "End") + "\n");
		U.close(fw);
	}

	public void saveToCSVFile(String fn) {
		if (fn != null) {
			output(fn, "\"", "\"", ",");
		}
	}

	public void print() {
		output(null, "", "", "\t");
	}
}

import java.util.ArrayList;
import java.util.List;
import java.util.Random;
import org.leores.math.Calculateable;
import org.leores.math.CurveFitter;
import org.leores.math.rand.Binomial;
import org.leores.math.rand.DRand;
import org.leores.math.rand.Exponential;
import org.leores.math.rand.Gamma;
import org.leores.math.rand.Normal;
import org.leores.math.rand.Poisson;
import org.leores.math.rand.RandomEngine;
import org.leores.plot.JGnuplot;
import org.leores.plot.JGnuplot.Plot;
import org.leores.util.TestUtil;
import org.leores.util.U;
import org.leores.util.able.Processable0;
import org.leores.util.able.Processable2;
import org.leores.util.data.DataTable;
import org.leores.util.data.DataTableSet;
import org.leores.util.data.Distribution;

package org.leores.demo;



public class MathDemo extends Demo {

	public void randDist() {
		Distribution dist = new Distribution();
		Distribution distSys = new Distribution();
		long seed = System.currentTimeMillis();
		RandomEngine rand = RandomEngine.makeDefault(seed);
		Random randSys = new Random(seed);
		for (int i = 0; i < 1000; i++) {
			dist.put(rand.nextDouble());
			distSys.put(randSys.nextDouble());
		}
		dist.calDist(0.001);
		dist.saveToCSVFile("distRand.csv");
		distSys.calDist(0.001);
		distSys.saveToCSVFile("distRandSys.csv");
	}

	public void randSpeed() {
		final long seed = System.currentTimeMillis();
		final RandomEngine randTWF = RandomEngine.makeDefault(seed);
		final RandomEngine randDR = new DRand(seed);
		final Random randSys = new Random(seed);

		Runnable rTWF = new Runnable() {
			public void run() {
				randTWF.nextDouble();
			}
		};
		Runnable rDR = new Runnable() {
			public void run() {
				randDR.nextDouble();
			}
		};
		Runnable rSys = new Runnable() {
			public void run() {
				randSys.nextDouble();
			}
		};
		TestUtil.compareRunTime(10000, 2, rTWF, rSys, rDR);
	}

	public void genRandNumWithDistribution() {
		RandomEngine rand = RandomEngine.makeDefault();
		Binomial binomial = new Binomial(100, 0.5, rand);
		Poisson poisson = new Poisson(10, rand);
		Normal normal = new Normal(10, 1, rand);
		Exponential exponential = new Exponential(1, rand);
		Gamma gamma = new Gamma(1, 1, rand);
		for (int i = 0; i < 100; i++) {
			log(binomial.nextInt(), poisson.nextInt(), normal.nextDouble(), exponential.nextDouble(), gamma.nextDouble());
		}
	}

	class BinomialSim {
		public int n;
		public double p;
		public RandomEngine rand;

		public BinomialSim(int n, double p, RandomEngine rand) {
			this.n = n;
			this.p = p;
			this.rand = rand;
		}

		public int nextInt() {
			int x = 0;
			for (int i = 0; i < n; i++) {
				if (rand.nextDouble() < p)
					x++;
			}
			return x;
		}

		//The implementation is a variant of Luc Devroye's "Second Waiting Time Method" on page 522 of his text "Non-Uniform Random Variate Generation."
		public int nextInt2() {
			double log_q = Math.log(1.0 - p);
			int x = 0;
			double sum = 0;
			while (true) {
				sum += Math.log(rand.nextDouble()) / (n - x);
				if (sum < log_q) {
					return x;
				}
				x++;
			}
		}
	}

	public void binomialSpeedDist() {
		final RandomEngine rand = RandomEngine.makeDefault();
		//		final int n = 5000;
		//		final double p = 9.4E-8;
		final int n = 10000;
		final double p = 100d / n;
		final Binomial biCal = new Binomial(n, p, rand);
		final BinomialSim biSim = new BinomialSim(n, p, rand);

		log(biCal.nextInt(), biSim.nextInt());

		Processable0<Double> pCal = new Processable0<Double>() {
			public Double process() {
				return (double) biCal.nextInt();
			}
		};
		Processable0<Double> pCalS = new Processable0<Double>() {
			public Double process() {
				return (double) biCal.nextIntS();
			}
		};
		Processable0<Double> pSim = new Processable0<Double>() {
			public Double process() {
				return (double) biSim.nextInt();
			}
		};

		TestUtil.compareRunTime(1000, 3, pCal, pCalS, pSim);

		Distribution distCal = new Distribution("distCal");
		distCal.prep(1.0, 0.0, 200.0);
		Distribution distCalS = new Distribution("distCalS");
		distCalS.prep(1.0, 0.0, 200.0);
		Distribution distSim = new Distribution("distSim");
		distSim.prep(1.0, 0.0, 200.0);
		for (int i = 0; i < 10000; i++) {
			distCal.stat((double) biCal.nextInt());
			distSim.stat((double) biSim.nextInt());
			distCalS.stat((double) biCal.nextIntS());
		}
		distCal.calDist();
		distSim.calDist();
		distCalS.calDist();

		JGnuplot jg = new JGnuplot();
