/*
 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;
import org.apache.commons.math3.util.FastMath;

/**
 * A class for eccentricity calculations using a principle component analysis.
 * It calculates variances alone the principal axis and the eccentricity.
 * 
 * @author S.Chekanov
 * 
 */

public class EEcentricity {

	double varMin = 1;
	double varMax = 1;
	double Theta = 0;
	double w[];
	double x[];
	double y[];

	/**
	 * Initialize 2D PCA analysis
	 * 
	 * @param x
	 *            array in X
	 * @param y
	 *            array in Y (the same size as for X)
	 * @param w
	 *            points weights (can be set to 1 if all points are treated with
	 *            the same weight)
	 */

	public EEcentricity(double x[], double y[], double w[]) {

		this.x = x;
		this.y = y;
		this.w = w;
		eval();

	}

	/**
	 * Perform PCA analysis using P1D object (in 2D). All weights for points are
	 * set to 1. X and Y component of P1D are used for the PCA.
	 * 
	 * @param p1d
	 *            P1D input objects
	 */

	public EEcentricity(P1D p1d) {

		int size = p1d.size();
		w = new double[size];
		for (int i = 0; i < size; i++)
			w[i] = 1.0;
		x = p1d.getArrayX();
		y = p1d.getArrayY();
		eval();

	}

	/**
	 * Main method for calculations
	 */

	private void eval() {

		int num = x.length;
		double etaSum = 0.;
		double phiSum = 0.;
		double eTot = 0.;

		double Dphi[] = new double[num];
		double Deta[] = new double[num];

		for (int j = 0; j < num; j++) {
			etaSum = etaSum + x[j] * w[j];
			phiSum = phiSum + y[j] * w[j];
			eTot = eTot + w[j];
		}

		// get differences from weighted mean
		etaSum = etaSum / eTot;
		phiSum = phiSum / eTot;
		for (int j = 0; j < num; j++) {
			Deta[j] = x[j] - etaSum;
			Dphi[j] = y[j] - phiSum;
		}

		double X1 = 0.;
		double X2 = 0;
		for (int i = 0; i < num; i++) {
			X1 += 2. * w[i] * Deta[i] * Dphi[i];
			X2 += w[i] * (Dphi[i] * Dphi[i] - Deta[i] * Deta[i]);
		}

		// variance calculations
		Theta = .5 * FastMath.atan(X1 / X2);
		double sinTheta = FastMath.sin(Theta);
		double cosTheta = FastMath.cos(Theta);
		double VarX = 0.;
		double VarY = 0.;

		for (int i = 0; i < num; i++) {
			double X = cosTheta * Deta[i] - sinTheta * Dphi[i];
			double Y = sinTheta * Deta[i] + cosTheta * Dphi[i];
			VarX += w[i] * X * X;
			VarY += w[i] * Y * Y;
		}

		varMax = VarX;
		varMin = VarY;
		if (varMax < varMin) {
			varMax = VarY;
			varMin = VarX;
		}

	}

	/**
	 * Returns angle of ellipse rotation (angle of the principial axis)
	 * 
	 * @return rotation angle (in rad)
	 */
	public double getAngle() {
		return Theta;
	}

	/**
	 * Get variance alone the minor (shortest component)
	 * 
	 * @return minimum variance
	 */
	public double getVarianceMin() {
		return varMin;

	}

	/**
	 * Get variance alone the major axis (longest component)
	 * 
	 * @return minimum variance
	 */
	public double getVarianceMax() {
		return varMax;

	}

	/**
	 * Get eccentricity, defined as 1-VarianceMin / VarianceMax. The value 0
	 * corresponds to perfectly circular distribution, while 1 is for max.
	 * allongated
	 * 
	 * @return eccentricity
	 */
	public double getEccentricity() {
		return 1 - (varMin / varMax);

	}

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

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

	}

}