package edu.stanford.rsl.conrad.calibration;

import java.util.ArrayList;

import Jama.EigenvalueDecomposition;
import Jama.Matrix;
import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND;

/**
 * Implements the Principal Component Analysis for a set of PointND.java, restricted to 2D. First the centroid S is computed. 
 * Then the set of points is weighted and the covariance matrix C is built. Eigenvalues and eigenvectors are estimated using
 * the Eigenvalue Decomposition. Imports Jama.Matrix and Jama.EigenvalueDecomposition.
 * 
 * @author Philipp Roser
 */
public class PrincipalComponentAnalysis2D {
	
	/**
	 * double array containing the eigenvalues of the covariance matrix C. Should be of length 2.
	 */
	private double[] eigenvalues;
	
	/**
	 * Jama.Matrix array containing the eigenvectors of the covariance matrix C. Should be of length 2.
	 */
	private Matrix[] eigenvectors;
	
	/**
	 * Covariance matrix C of the data set.
	 */
	private Matrix C;
	
	/**
	 * Centroid S of the data set.
	 */
	private PointND S;
	
	/**
	 * PointND array containing the data to be analyzed.
	 */
	private PointND[] set;
	
	/**
	 * PointND array containing the data translated by the centroid.
	 */
	private PointND[] weightedSet;

	/**
	 * Constructor, expecting the data set to be analyzed as array PointND[]. Executes the complete Principal Component Analysis.
	 * 
	 * @param set
	 */
	public PrincipalComponentAnalysis2D(PointND[] set) {
		this.set = set;
		computeS();
		weightSet();
		buildC();
		computeEigenvalues();
		computeEigenvectors();
	}
	
	
	/**
	 * Constructor, expecting the data set to be analyzed as ArrayList<PointND>. Executes the complete Principal Component Analysis.
	 * 
	 * @param set
	 */
	public PrincipalComponentAnalysis2D(ArrayList<PointND> set) {
		this.set = new PointND[set.size()];
		for (int i = 0; i < this.set.length; i++) {
			this.set[i] = new PointND(set.get(i));
		}
		computeS();
		weightSet();
		buildC();
		computeEigenvalues();
		computeEigenvectors();
	}

	/**
	 * 
	 * @return Principal Components of the data set as Jama.Matrix[] of length 2.
	 */
	public Matrix[] getEigenvectors() {
		return this.eigenvectors;
	}

	/**
	 * 
	 * @return Principal Components of the data set as String.
	 */
	public String stringEigenvectors() {
		String result = "";
		for (int i = 0; i < eigenvectors.length; i++) {
			result += "e" + i + " = (" + eigenvectors[i].get(0, 0) + ", "
					+ eigenvectors[i].get(1, 0) + ")\n";
		}
		return result;
	}

	/**
	 * 
	 * @return eigenvalues of the covariance matrix
	 */
	public double[] getEigenvalues() {
		return this.eigenvalues;
	}

	/**
	 * Computes the centroid S of the data set.
	 */
	private void computeS() {

		double xsum = 0.0;
		double ysum = 0.0;

		for (int i = 0; i < set.length; i++) {
			xsum = xsum + set[i].get(0);
			ysum = ysum + set[i].get(1);
		}

		S = new PointND(xsum / set.length, ysum / set.length);
	}

	/**
	 * Translates all points in the data set by the centroid S.
	 */
	private void weightSet() {
		weightedSet = new PointND[set.length];
		for (int i = 0; i < set.length; i++) {
			weightedSet[i] = new PointND(set[i].get(0) - S.get(0),
					set[i].get(1) - S.get(1));
		}
	}

	/**
	 * Builds the covariance matrix.
	 */
	private void buildC() {
		C = new Matrix(2, 2);
		int n = set.length;
		for (int i = 0; i < n; i++) {
			Matrix P = new Matrix(2, 2);
			P.set(0, 0, ((double) 1.0 / (n - 1)) * weightedSet[i].get(0)
					* weightedSet[i].get(0));
			P.set(0, 1, ((double) 1.0 / (n - 1)) * weightedSet[i].get(0)
					* weightedSet[i].get(1));
			P.set(1, 0, ((double) 1.0 / (n - 1)) * weightedSet[i].get(1)
					* weightedSet[i].get(0));
			P.set(1, 1, ((double) 1.0 / (n - 1)) * weightedSet[i].get(1)
					* weightedSet[i].get(1));
			C = C.plus(P);
		}

	}

	/**
	 * Estimates the eigenvalues of the covariance Matrix C using Jama.EigenvalueDecomposition.
	 */
	private void computeEigenvalues() {
		EigenvalueDecomposition eig = C.eig();
		eigenvalues = eig.getRealEigenvalues();
	}

	/**
	 * Estimates the eigenvectors of the covaricance matrix C using Jama.EigenvalueDecomposition.
	 */
	private void computeEigenvectors() {
		// Matrix c1 = C.minus(Matrix.identity(2, 2).times(eigenvalues[0]));
		// Matrix c2 = C.minus(Matrix.identity(2, 2).times(eigenvalues[1]));
		EigenvalueDecomposition m = C.eig();
		Matrix eig = m.getV();
		eigenvectors = new Matrix[2];
		eigenvectors[0] = eig.getMatrix(0, 1, 0, 0);
		eigenvectors[1] = eig.getMatrix(0, 1, 1, 1);
		/*
		 * eigenvectors[0] = new Matrix(2, 1); eigenvectors[0].set(0, 0,
		 * 1.0); eigenvectors[0].set(1, 0, -1.0 * c1.get(0, 0) / c1.get(0,
		 * 1)); eigenvectors[1] = new Matrix(2, 1); eigenvectors[1].set(0,
		 * 0, 1.0); eigenvectors[1].set(1, 0, -1.0 * c2.get(0, 0) /
		 * c2.get(0, 1));
		 */
	}

}