/* 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 { ILattice lattice; public HistoCurve2D(ILattice lattice){ this.lattice = lattice; } public synchronized Point2D getPoint(int i){ if (i>=entries.length) return new Point2D(); double x = xAxis.binCentre(i)/lattice.getSpinCount(); 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)/lattice.getSpinCount(); 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)/lattice.getSpinCount(); case Y_AXIS: return highestEntry; default: throw new IllegalArgumentException("Wrong axis index : "+axis); } } } }