/*  Monte-Carlo simulation code for statistical physics
    Copyright (C) 2001-2005  Sylvain Reynal
    
    Département de Physique
    Ecole Nationale Supérieure de l'Electronique et de ses Applications (ENSEA)
    6, avenue du Ponceau,    
    F-95014 CERGY CEDEX

    et
    
    Laboratoire de Physique Théorique et Modélisation (LPTM)
    Université de Cergy-Pontoise - Site de Neuville
    F-95031 CERGY CEDEX
    
	
    Tel : 00 +33 130 736 245
    Fax : 00 +33 130 736 667
    e-mail : reynal@ensea.fr
    web page : http://www.ensea.fr/staff/reynal/
*/

package fr.ensea.montecarlo.data;

import com.imsl.math.Complex;
import fr.ensea.chart.Curve2D;
import fr.ensea.chart.Point2D;
import fr.ensea.math.FFT;
import java.util.Vector;
import java.awt.Color;

/**
 * Cette classe permet de stocker une série d'échantillons d'une grandeur thermodynamique obtenus au cours d'une simulation,
 * puis de calculer quelques estimateurs statistiques de cette série.
 *
 */
public class Samples extends Vector implements Curve2D {
	
	private double minX, minY, maxX, maxY;
	private Histogram histo;
	private boolean autocomputeCorrelationTime=false;
	private Double correlationTime;
	private int mcsWait=1, mcsTherm=0;
	private final static double LOG2 = Math.log(2);
	
	/**
	 * Construit un nouveau conteneur d'échantillon, initialement vide.
	 */
	public Samples(){
		minX= minY= maxX= maxY=0; // pas encore initialisé
	}
	
	/**
	 * Ajoute un nouvel échantillon
	 */
	public synchronized void addSample(double sample){
		addElement(new Double(sample));
		if (size()==1){ // premier élément ajouté
			minY = maxY = sample;
		}
		else {
			if (sample < minY) minY = sample;
			else if (sample > maxY) maxY = sample;
		}
		maxX = size();
		if (histo != null) histo.updateIncremental(sample);
		// recompute correlation time ?
		if (autocomputeCorrelationTime){
			double n = log2(size());
			if (n==Math.floor(n)) computeCorrelationTime();
		}
	}
	
	/**
	 * Retourne le ième échantillon
	 */
	public synchronized double getSample(int i){
		if (i>=size()) return 0.0;
		return ((Double)elementAt(i)).doubleValue();
	}
	
	public void clear(){
		super.removeAllElements();
		if (histo != null) histo.update();
	}
	
	///////////////////////////////////////
	//// Statistiques
	///////////////////////////////////////
	
	/**
	 * Calcule et retourne <x>
	 */
	public synchronized double computeMean(){
		if (isEmpty()) return 0.0;
		int inc=0;
		if (mcsWait<1) inc = (int)(2.0*getCorrelationTime());
		else inc = mcsWait;
		double mean = 0;
		int n=0;
		for (int i=mcsTherm; i<size(); i+=inc){
			mean += getSample(i);
			n++;
		}
		if (n==0) return 0.0; // it's the case if size() < mcsTherm
		return mean/n;
	}

	/**
	 * Calcule et retourne <x^2> - <x><x> (i.e. les fluctuations thermodynamiques de la grandeur thermodynamique)
	 * @param mean la moyenne de la série d'échantillons, préalablement obtenue en invoquant "computeMean()".
	 */
	public synchronized double computeVariance(double mean){
		if (isEmpty()) return 0.0;
		int inc=0;
		if (mcsWait<1) inc = (int)(2.0*getCorrelationTime());
		else inc = mcsWait;
		double sumOfSquaress = 0;
		int n=0;
		for (int i=mcsTherm; i<size(); i+=inc){
			double ech = getSample(i);
			sumOfSquaress += ech * ech;
			n++;
		}
		if (n==0) return 0.0; // it's the case if size() < mcsTherm
		return sumOfSquaress/n - mean * mean;
	}
	
	/**
	 * Créé puis retourne un histogramme de cette série d'échantillons.
	 * L'histogramme est mis-à-jour à chaque nouvel ajout d'échantillons.
	 */
	public Histogram createHistogram(int nbBins){
		histo = new Histogram(this, nbBins);
		return histo;
	}
	
	///////////////////////////////////////
	//// Autocorrelation time
	///////////////////////////////////////
	
	/**
	 * If b is true, the cached autocorrelation-time will be computed anew each time size() crosses a 2^n border.
	 * False by default, in which case computeCorrelationTime() must be called explicitely when needed.
	 */
	public void setAutocomputeCorrelationTime(boolean b){
		this.autocomputeCorrelationTime = b;
	}
	
	/**
	 * Sets the number of steps between independent samples. This is used for computing mean values.
	 * @param mcsWait set to -1 to compute mcsWait automatically from the autocorrelation time. 
	 */
	public void setMCSWait(int mcsWait){
		this.mcsWait = mcsWait;
	}
	
	/**
	 * Returns the number of steps between independent samples. This is used for computing mean values.
	 * @return -1 if we compute mcsWait automatically from the autocorrelation time. 
	 */
	public int getMCSWait(){
		return this.mcsWait;
	}

	/**
	 * Sets the number of samples to be discarded at the beginning of the serie. 
	 */
	public void setThermalizationTime(int mcsTherm){
		if (mcsTherm < 0) mcsTherm=0;
		this.mcsTherm = mcsTherm;
	}
	
	public int getThermalizationTime(){
		return this.mcsTherm;
	}

	/**
	 * Return the cached correlation-time for this set of data. Invoke computeCorrelationTime() if you need to update
	 * the cache.
	 * @see #setAutocomputeCorrelationTime
	 */
	public double getCorrelationTime(){
		if (this.correlationTime == null) computeCorrelationTime();
		return this.correlationTime.doubleValue();
	}
	
	private static double log2(double x){
		return Math.log(x)/LOG2;
	}
	
	/**
	 * Re-compute the cached correlation-time (using integration of the normalized time-displaced auto-correlation 
	 * function) for the given Samples
	 * - chi(t) = [sum_i m(i) m(i+t) - <m>²]/chi(0)
	 * The computed time may then be retrieved by calling getCorrelationTime();
	 */
	public synchronized void computeCorrelationTime(){
		if (size()==0) {
			this.correlationTime = new Double(1.0);
			return;
		}

		// ensure length= 2**n :
		int n = (int)Math.floor(log2(size()));
		if (n < 5) {
			this.correlationTime = new Double(1.0);
			return; // if less than 32, not reliable => return 1.0
		}

		int npts=1;
		for (int i=0; i<n ; i++) npts *= 2;

		//System.out.println("Number of FFT pts="+npts);

		// take subpart of samples from 0 to 2**n-1 and compute mean value:
		double[] data = new double[npts];
		double mean=0;
		for (int i=0; i<data.length; i++) {
			data[i] = getSample(i);
			mean += data[i];
		}
		mean /= data.length;
		
		// remove mean value :
		for (int i=0; i<data.length; i++) data[i]-=mean;
		
		Complex[] cArray = FFT.transform(data);
		// compute fft*fft :
		for (int i=0; i<data.length; i++){
			data[i]=Complex.abs(cArray[i]); // data = abs(fft2)
			data[i]*=data[i];
			//System.out.println("fft2="+data[i]);
		}
		cArray = FFT.inverseTransform(data); // ifft(abs(fft)**2)
		data[0]=Complex.real(cArray[0]);
		for (int i=1; i<data.length; i++){
			data[i]=Complex.real(cArray[i])/data[0]; // now data = correlation function
			//System.out.println("chi="+data[i]);
		}
		data[0]=1.0;
		double[] chi = data;

		// correlation time :
		double tau=0.0;
		for (int i=0 ; i<chi.length && chi[i] > 0.1; i++){ // cut at 10%
			tau += chi[i];
		}
		this.correlationTime = new Double(tau);
	}
	
	
	////////////////////////////////////////////////////////////////////
	// test
	////////////////////////////////////////////////////////////////////
	public static void main(String[] args){
		Samples samples = new Samples();
		for (int i=0; i<100; i++){
			samples.addSample(Math.random());
		}
		double mean = samples.computeMean();
		double variance = samples.computeVariance(mean);
		System.out.println("mean="+mean);
		System.out.println("variance="+variance);
	}
	
	/*
	public String toString(){
		String s = "Samples:\n";
		s += "* size:" + size() + "\n";
		s += "* x-axis: min="+minX+" max="+maxX+"\n";
		s += "* y-axis: min="+minY+" max="+maxY+"\n";
		double mean = computeMean();
		s += "* mean: "+mean+"\n";
		s += "* variance: "+ computeVariance(mean)+"\n";
		return s;
	}
	*/
		
	/////////////////////////////////////////////////
	// Implémentation de l'interface "Courbe"
	/////////////////////////////////////////////////
	
	/**
	 * returns the ith point
	 */
	public Point2D getPoint(int i){
		return new Point2D(i, getSample(i));
	}
	
	/**
	 * Return the number of points
	 */
	public int getNbPoints(){
		return size();
	}
	
	/**
	 * Returns the minimum on the given axis
	 */
	public double getMin(int axis){
		switch (axis){
			case X_AXIS: return minX;
			case Y_AXIS: return minY;
			default: throw new IllegalArgumentException("Wrong axis index : "+axis);
		}
	}
	
	/**
	 * Returns the maximum on the given axis
	 */
	public double getMax(int axis){
		switch (axis){
			case X_AXIS: return maxX;
			case Y_AXIS: return maxY;
			default: throw new IllegalArgumentException("Wrong axis index : "+axis);
		}
	}
}
