/*  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 fr.ensea.chart.*;
import fr.ensea.montecarlo.model.*;
import hep.aida.ref.*;
import cern.colt.list.*;

/**
 * A 1D histogram with dynamically extensible x-axis.
 */
public class DynamicHistogram1D {

	// Note: methods related to bin heights as well as RMS have been commented out for performance reasons

	protected FixedAxis xAxis;

	// inherited: OVERFLOW=-1 ; UNDERFLOW=-2
	protected int underflowBin, overflowBin; // UNDERFLOW and OVERFLOW bins
	protected int[] entries; // in-range bins only
	protected int inRangeEntries; // nb of in-range entries (i.e. exluding UNDERFLOW and OVERFLOW)
	protected int highestEntry; // value of highest entry

	protected boolean isAxisOK;
	protected double binWidth;
	
	public static final int UNDERFLOW=-2;
	public static final int OVERFLOW=-1;

	/**
	 * Creates a dynamic histogram with the given axis binning
	 */
	public DynamicHistogram1D(double binWidth){
		this(new FixedAxis(1,-1,1));
		this.binWidth=binWidth;
		this.isAxisOK=false; // i.e. xAxis not set, even though we fed it into the super constructor
	}

	/**
	 * Creates a histogram with the given axis binning
	 */
	public DynamicHistogram1D(FixedAxis axis){
		xAxis = axis;
		int bins = axis.bins();
		entries = new int[bins];
		this.binWidth=axis.binWidth(0);
		isAxisOK=true;
		highestEntry=0;
	}

	/**
	 * Create a fixed-width histogram.
	 */
	public DynamicHistogram1D(int bins, double min, double max){
		this(new FixedAxis(bins,min,max));
	}

	/////////////////////////////////////////////////////////////////////////////////
	// Axis
	/////////////////////////////////////////////////////////////////////////////////

	/**
	 * Return the x-axis.
	 */
	public FixedAxis xAxis(){
		return xAxis;
	}

	/**
	 * Assign a new axis to this histogram, resets content, and sets the isAxisOK flag so
	 * that subsequent calls to fill() may enlarge the axis properly.
	 */
	public void resetAxis(FixedAxis axis){
		this.binWidth = axis.binWidth(0);
		this.xAxis = axis;
		this.entries = new int[xAxis.bins()];
		isAxisOK = true;
	}

	public double getBinWidth(){
		return binWidth; // same as xAxis().binWidth(0) since this is a FixedAxis
	}

	/**
	 * Change the bin width, marks the x-axis as invalid, and reset entries. 
	 */
	public void resetBinWidth(double newBinWidth){
		this.binWidth = newBinWidth;
		this.xAxis = new FixedAxis(1,-1,1);
		this.entries = new int[1];
		highestEntry=0;
		inRangeEntries=0;
		isAxisOK=false;
	}

	/**
	 * Updates the x-axis so that the given value fits in the in-range bins (if it doesn't already do so)
	 * @return the bin index corresponding to the given double, after axis has been updated
	 */
	public int updateAxis(double x){
		if (!isAxisOK){ // axis not set yet (we just fed the constructor with a fake axis containing one bin in order to avoir NullPointerException's here and there)
			initAxis(x);
			return 0; // first and only bin, hence index = 0
		}
		else {
			int xAxisBin = xAxis.coordToIndex(x);
			if (xAxisBin == UNDERFLOW){ // enlarge xAxis to the left, prepending as many bins as necessary
				prependBins(x);
				return xAxis.coordToIndex(x); // not the same as before since bins've been preprended !!!
			}
			else if (xAxisBin == OVERFLOW){ // enlarge xAxis to the right: (thus no need to shift array here)
				appendBins(x);
				return xAxis.coordToIndex(x); // ibid.
			}
			return xAxisBin;
		}
	}

	/**
	 * Init the x-axis from the given bin center, resetting bin entries to their default value.
	 */
	protected void initAxis(double x){
		double min = x - 0.5*binWidth; // "x" will be the center of the initial bin
		double max = x + 0.5*binWidth;
		int bins=1; // unchanged wrt constructor, that's fine
		this.xAxis = new FixedAxis(bins, min, max); // change just min and max
		this.entries = new int[1];
		highestEntry=0;
		inRangeEntries=0;
		isAxisOK=true;
	}

	/**
	 * Prepends new bins to the left of the current axis, and shift array accordingly.
	 *
	 * We allocate a new entries[] array with size = new_nb_bins, and copy old entries[] array into it:
	 * - new_entries[0:new_nb_bins-old_nb_bins-1]=inited to 0 (new bins)
	 * - new_entries[new_nb_bins-old_nb_bins:new_nb_bins-1]=old_entries[0:old_nb_bins-1] (copy of old values at shifted position)
	 *
	 * Examples: old_nb_bins=2; new_nb_bins=4
	 *
	 * old_entries: 							[0]=first_bin	[1]=last_bin
	 *                   v (prepend.) v          v  (copied)     v
	 * new_entries:  	[0]=0 		 [1]=0 		[2]=old[0] 		[3]=old[1]
	 *
	 * underflowBin and overflowBin left unchanged (should be empty if everything went fine till now)
	 *
	 * @param x the new entry the new (enlarged) axis should include; must be lower than xAxis.lowerEdge() 
	 */
	protected void prependBins(double x){
		double max = xAxis.upperEdge(); // unchanged
		int newNbBins = 1+(int)Math.floor((max - x)/binWidth); // new nb of bins
		int oldNbBins = xAxis.bins();
		if (newNbBins == oldNbBins) return;
		if (newNbBins < oldNbBins) throw new IllegalArgumentException("when prepending bins, new nb of bins must be greater than old one !");

		double newMin = max - newNbBins * binWidth; // new axis's min value
		this.xAxis = new FixedAxis(newNbBins, newMin, max); // new axis

		int[] _entries = new int[newNbBins];

		int destPos = newNbBins-oldNbBins;
		System.arraycopy(entries,0,_entries, destPos,oldNbBins); // don't copy neither UNDERFLOW nor OVERFLOW, they should be 0 !
		entries=_entries;
	}

	/**
	 * Appends new bins to the right of the current axis.
	 *
	 * We allocate a new entries[] array with size = new_nb_bins+2, and copy old entries[] array into it:
	 * - new_entries[0:old_nb_bins-1]=old_entries[0:old_nb_bins-1] (unchanged, copy of old values at the same position)
	 * - new_entries[old_nb_bins:new_nb_bins-1]=inited to 0 (new bins)
	 *
	 * @param x the entry which the new (enlarged) axis should include  
	 */
	protected void appendBins(double x){
		double min = xAxis.lowerEdge(); // unchanged
		int newNbBins = 1+(int)Math.floor((x - min)/binWidth); // new nb of bins
		int oldNbBins = entries.length;
		if (newNbBins == oldNbBins) return;
		if (newNbBins < oldNbBins) throw new IllegalArgumentException("when appending bins, new nb of bins must be greater than old one !");

		double newMax = min + newNbBins * binWidth; // new axis's max value
		this.xAxis = new FixedAxis(newNbBins, min, newMax); // new axis

		int[] _entries = new int[newNbBins];

		System.arraycopy(entries,0,_entries,0,oldNbBins); // don't copy neither UNDERFLOW nor OVERFLOW, they should be 0 !
		entries=_entries;
	}

	/////////////////////////////////////////////////////////////////////////////////
	// Filling and related
	/////////////////////////////////////////////////////////////////////////////////

	/**
	 * Reset contents; as if just constructed.
	 */
	public void reset(){
		for (int i=entries.length; --i >= 0;) entries[i] = 0;
		underflowBin=0;
		overflowBin=0;
		inRangeEntries = 0;
		highestEntry=0;
	}

	/**
	 * Fill the given bin index
	 */
	protected synchronized void fillFromBin(int bin){
		switch (bin){
		case UNDERFLOW: underflowBin++; break; // -1
		case OVERFLOW: overflowBin++; break; // -2
		default: 
			entries[bin]++;
			inRangeEntries++; 
			if (entries[bin]>highestEntry) highestEntry = entries[bin];
			break;
		}
	}

	/**
	 * Fill histogram with weight 1
	 */
	public synchronized void fill(double x){
		int bin = updateAxis(x);
		fillFromBin(bin);
	}

	/**
	 * Number of in-range entries in the histogram.
	 */
	public int entries(){
		return inRangeEntries;
	}

	/**
	 * Number of all entries in all (both in-range and under/overflow) bins in the histogram. 
	 */
	public int allEntries(){
		return underflowBin+inRangeEntries+overflowBin;
	}

	/**
	 * Number of under and overflow entries in the histogram.
	 */
	public int extraEntries(){
		return underflowBin + overflowBin;
	}

	/////////////////////////////////////////////////////////////////////////////////
	// Bin entries and related bin properties
	/////////////////////////////////////////////////////////////////////////////////

	/**
	 * Returns the centre of the bin specified.
	 */
	public double binCenter(int binIndex){
		return xAxis.binCentre(binIndex);
	}

	/**
	 * Returns the nb of bins (excluding UNDERFLOW and OVERFLOW)
	 */
	public int bins(){
		return xAxis.bins();
	}

	/**
	 * Number of entries in the corresponding bin (ie the number of times fill was called for this bin)
	 */
	public int binEntries(int index){
		switch (index){
		case UNDERFLOW: return underflowBin;
		case OVERFLOW: return overflowBin;
		default: return entries[index];
		}
	}

	/**
	 * Indexes of the in-range bins containing the smallest and largest binEntries(), respectively.
	 */
	public int[] minMaxBins(){
		return minMaxBins(entries.length);
	}

	/**
	 * Indexes of the in-range bins containing the smallest and largest binEntries(), respectively, yet
	 * including bins below the given bin index only (inclusive). 
	 */
	public int[] minMaxBins(int maxBinIdx){
		int minValue = Integer.MAX_VALUE;
		int maxValue = Integer.MIN_VALUE;
		int minBinX = -1;
		int maxBinX = -1;
		final int IN_RANGE_BINS = entries.length;
		maxBinIdx = (maxBinIdx > IN_RANGE_BINS ? IN_RANGE_BINS : maxBinIdx);
		for (int i=maxBinIdx; --i >= 0; ) {
			int value = entries[i];
			if (value < minValue) {
				minValue = value;
				minBinX = i;
			}
			if (value > maxValue) {
				maxValue = value;
				maxBinX = i;
			}
		}
		int[] result = {minBinX,maxBinX};
		return result;
	}

	/**
	 * Return indices of non-empty bins, in ascending order.
	 */
	public IntArrayList getIndicesOfNonEmptyBins(){
		return getIndicesOfNonEmptyBins(entries.length);
	}

	/**
	 * Return indices of non-empty bins below the given bin index (exclusive), in ascending order.
	 * Note : used by BergEnergyInterpolatorUpdater.
	 */
	public IntArrayList getIndicesOfNonEmptyBins(int maxBinIdx){
		IntArrayList indices = new IntArrayList();
		final int IN_RANGE_BINS = entries.length;
		maxBinIdx = (maxBinIdx > IN_RANGE_BINS ? IN_RANGE_BINS : maxBinIdx);
		for (int bin=0; bin < maxBinIdx; bin++){
			if (entries[bin] > 0) indices.add(bin);
		}
		//indices.sortFromTo(0,indices.size()-1);
		return indices;
	}

	/**
	 * Return the index of the lowest non-empty bin
	 * @return -1 if all in-range bins are empty
	 */
	public int getLowestNonEmptyBin(){
		final int IN_RANGE_BINS = entries.length;
		for (int bin=0; bin < IN_RANGE_BINS; bin++){
			if (entries[bin] > 0) return bin;
		}
		return -1;
	}

	/**
	 * Return the index of the highest non-empty bin
	 * @return -1 if all in-range bins are empty
	 */
	public int getHighestNonEmptyBin(){
		for (int bin=entries.length; --bin >= 0; ) {
			if (entries[bin] > 0) return bin;
		}
		return -1;
	}

	/////////////////////////////////////////////////////////////////////////////////
	// Statistics and related properties
	/////////////////////////////////////////////////////////////////////////////////



	/**
	 * Compute the flatness value of this histogram
	 * @return 0.0 for a perfectly flat histogram
	 */
	public double flatness(){
		double entriesMean = 0.0;
		double entriesMean2 = 0.0;
		final int maxBinIdx = xAxis().bins();
		for (int bin=0; bin < maxBinIdx; bin++){
			int binHeight = binEntries(bin);
			entriesMean += binHeight;
			entriesMean2 += binHeight * binHeight;
		}
		entriesMean /= maxBinIdx;
		entriesMean2 /= maxBinIdx;
		double flatness= Math.sqrt(entriesMean2 - entriesMean*entriesMean)/entriesMean;

		return flatness;
	}

	/**
	 * Compute the flatness value of this histogram, restricted over non-empty bins
	 * @return 0.0 for a perfectly flat histogram
	 */
	public double flatnessOverNonEmptyBins(){
		final int maxBinIdx = xAxis().bins();
		double entriesMean = 0.0;
		double entriesMean2 = 0.0;
		IntArrayList nonEmptyBins = getIndicesOfNonEmptyBins(maxBinIdx);
		for (int i=0; i < nonEmptyBins.size(); i++){
			int bin = nonEmptyBins.get(i);
			int	binHeight = binEntries(bin);
			entriesMean += binHeight;
			entriesMean2 += binHeight * binHeight;
		}
		entriesMean /= nonEmptyBins.size();
		entriesMean2 /= nonEmptyBins.size();
		double flatness= Math.sqrt(entriesMean2 - entriesMean*entriesMean)/entriesMean;

		return flatness;
	}

	///////////////////////////////////////////////////////////////
	/// UI
	///////////////////////////////////////////////////////////////

	public HistoCurve2D createHistoCurve2D(ILattice lattice){
		return new HistoCurve2D(lattice);
	}

	public class HistoCurve2D implements Curve2D {

		int VOL;

		public HistoCurve2D(ILattice lattice){
			VOL = lattice.getSpinCount();
		}

		public synchronized Point2D getPoint(int i){
			double x = xAxis.binCentre(i)/VOL;
			double y = entries[i];
			return new Point2D(x,y);
		}

		public int getNbPoints(){
			return bins();
		}

		public double getMin(int axis){
			switch (axis){
			case X_AXIS: return xAxis.binCentre(0)/VOL;
			case Y_AXIS: return 0;
			default: throw new IllegalArgumentException("Wrong axis index : "+axis);
			}
		}

		public double getMax(int axis){
			switch (axis){
			case X_AXIS: return xAxis.binCentre(bins()-1)/VOL;
			case Y_AXIS: return highestEntry;
			default: throw new IllegalArgumentException("Wrong axis index : "+axis);
			}
		}
	}

}
