/*
 This is a proprietary part of DataMelt package.
 License: the Creative Commons Attribution-Share Alike License; either version 3.0 of the License, 
 or any later version.
 http://creativecommons.org/licenses/by-sa/3.0/
 

 Permission to use, copy and modify this software and its documentation for 
 NON-COMMERCIAL purposes is granted, 
 without fee, provided that an acknowledgment to the author, Dr Sergei Chekanov
 appears in all copies and associated documentation or publications.
 

 Redistributions of the source code, or parts of the source codes, 
 must retain the copyright notice, this list of conditions and the following disclaimer 
 (all at the top of each source code) and requires the written permission of Sergei Chekanov
 

 Redistribution in binary form of all or parts of a class must reproduce the 
 copyright notice, this list of conditions and the following disclaimer in the documentation and/or other 
 materials provided with the distribution and requires the written permission of Sergei Chekanov
 

 Details about license agreement  or commercial use are:  http://jwork.org/scavis
 */

package jhpro.stat;

import jhplot.P1D;
import jhplot.gui.HelpBrowser;

/**
 * Calculations of normalized factorial moments of a distribution. The
 * normalised factorial moments are calculated as up to an orbitrary order n,
 * where n is the oder of the moment. The factorial moments are defined as 

 * F_1=1 (1st order), 

 * F_2=/^2 (second order),

 * F_2=/N^3 (3rd order) 

 * and so on. The general form is F_q= /^q, where
 * N^[q]=N*(N-1)*..(N-q+1). The code works only for integer random variables.
 * Statistical errors are also evaluated. If they are set to 0, this means that
 * there is no enough events for their calculation.
 * 
 * @author S.Chekanov
 * 
 */
public class MomentsFacNorm {

	private int oder = 0;
	private double[] IL, IE, IEE, VAL, AVAL;
	private int nev;

	/**
	 * Initialization of the normalised factorial moments.
	 * 
	 * @param order
	 *            Max oder of the normalised factorial moment
	 */
	public MomentsFacNorm(int order) {

		nev = 0;
		this.oder = order;
		IL = new double[oder];
		IE = new double[oder];
		IEE = new double[oder];
		VAL = new double[oder];
		AVAL = new double[oder];

		for (int i = 0; i < oder; i++) {
			IL[i] = 0;
			IE[i] = 0;
			IEE[i] = 0;
			VAL[i] = 0;
			AVAL[i] = 0;
		}
	}

	/**
	 * Process current event for an integer random variable
	 * 
	 * @param n
	 *            input variable for the moment calculations
	 */
	public void process(int n) {

		nev++;
		for (int i = 0; i < oder; i++) {

			int h = n;
			for (int j = 0; j < i; j++)
				h = h * (n - 1 - j);
			IL[i] = IL[i] + h;
			IE[i] = IE[i] + h * h;
			IEE[i] = IEE[i] + h * n; // for covariance matrix

		}
		;

	}

	/**
	 * Get values of the normalised factorial moments (NFM).
	 * 
	 * @return aeeay with values of the NFM.
	 */

	public double[] getValues() {

		double events = (double) (nev - 1);
		// everage over
		for (int i = 0; i < oder; i++)
			AVAL[i] = IL[i] / events;

		// take the ratio
		for (int i = 0; i < oder; i++)
			VAL[i] = AVAL[i] / Math.pow(AVAL[0], i + 1);

		return VAL;

	}

	/**
	 * Get errors on the factorial moment calculations. Call first getValues()
	 * to process the errors correctly. If the call returns 0, this means there
	 * is no enough evens for the calculation.
	 * 
	 * @return array with statistical errors
	 */
	public double[] getErrors() {
		// use HEN-455 W.Metzger way

		double events = (double) (nev - 1);

		// calculate variances
		double[] IER = new double[oder];
		for (int i = 0; i < oder; i++) {
			double X1 = IL[i] / events;
			double X2 = IE[i] / events;
			double VAR = X2 - X1 * X1;
			IER[i] = 0;
			if (VAR > 0)
				IER[i] = VAR / events;

		}

		// calculate covariances
		double[] COV = new double[oder];
		for (int i = 0; i < oder; i++) {
			double x1 = IEE[i] / events;
			double x2 = (IL[i] / events) * (IL[0] / events);
			COV[i] = (x1 - x2) / events;
		}

		// error propogation
		double[] ERROR = new double[oder];
		for (int i = 0; i < oder; i++) {
			int n = i + 1;
			double x1 = IER[i] / Math.pow(AVAL[0], 2 * n);
			double x2 = n * n
					* ((AVAL[i] * AVAL[i]) / Math.pow(AVAL[0], 2 * n + 2))
					* IER[0];
			double x3 = 2 * n * (AVAL[i] / Math.pow(AVAL[0], 2 * n + 1))
					* COV[i];
			double VARTOT = x1 + x2 - x3;
			ERROR[i] = 0;
			if (VARTOT > 0)
				ERROR[i] = Math.sqrt(VARTOT);

		}

		return ERROR;

	}

	/**
	 * Get results in form of P1D.
	 * 
	 * 
	 * @return p1d with values and errors
	 */
	public P1D getResults() {

		P1D tmp = new P1D("Factorial Moments");
		double[] val = getValues();
		double[] err = getErrors();
		for (int i = 0; i < oder; i++) {
			tmp.add(i + 1, val[i], err[i]);
		}

		return tmp;

	}

	/**
	 * Show online documentation.
	 */
	public void doc() {

		String a = this.getClass().getName();
		a = a.replace(".", "/") + ".html";
		new HelpBrowser(HelpBrowser.JHPLOT_HTTP + a);

	}

}