/*  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)
    CNRS/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.multicanonical;

import fr.ensea.montecarlo.data.*;
import fr.ensea.chart.*;

/**
 * A factory that produces reweighted statistical averages from a sample set and a weight w(E) obtained from S(E).
 */
public class MomentsReweighter {

	private WLDOSHistogram SE;
	private int VOL;
	private SamplesBag samplesBag;
	private Samples energySamples; 
	private Samples magnetizationSamples;
	private ThermalAveragesSet set;
	public DefaultCurve2D meanEnergyVsKTCurve; // <E>
	public DefaultCurve2D meanMagVsKTCurve; // <M>
	public DefaultCurve2D susceptibilityVsKTCurve; // chi
	
	/**
	 * reweighting w/o using flat histogram. 
	 */
	public MomentsReweighter(WLDOSHistogram SE, SamplesBag samplesBag, ThermalAveragesSet set, double[] kTarray, int VOL){
		this.SE = SE;
		this.samplesBag = samplesBag;
		this.energySamples = samplesBag.getEnergySamples();
		this.magnetizationSamples = samplesBag.getMagnetizationSamples();
		this.VOL = VOL;
		this.set=set;
		this.meanEnergyVsKTCurve = set.meanEnergyVsKTCurve;
		this.meanEnergyVsKTCurve.clear();
		this.meanMagVsKTCurve = set.meanMagVsKTCurve;
		this.meanMagVsKTCurve.clear();
		this.susceptibilityVsKTCurve = set.susceptibilityVsKTCurve;
		this.susceptibilityVsKTCurve.clear();
		for (int i=0; i<kTarray.length; i++){
			//System.out.print((i+1)+"/"+kTarray.length+" ");
			computeMoments(kTarray[i]);
		}
	}
	

	/**
	 * Compute Z, F, and reweighted thermodynamical averages at the given temperature kT.
	 * @param filterIndices an array of indices indicating which samples will be taken into account ; this may serve either
	 *        for resampling or for bootstrapping.
	 */
	protected void computeMoments(double kT){
		MuCaWeighter mw = new MuCaWeighter(kT, SE);
			
		// 1°) compute partition function (up to a anti-overflow factor) ; Z will be used for averages later on 
		double Z = computeZ(kT,mw);
			
		// 2°) compute free energy F(kT) ; Note : we must compensate for the anti-overflow factor in Z !!!
		double freeEnergy = -kT * Math.log(Z) - kT * mw.getAntiOverflowK(); // yeah ! otherwise, F(T) just makes no sense !!!
		freeEnergy/=VOL; // per spin
		
		// 3°) compute rew'd averages, e.g. <A> = 1/Z * sum_E <A>(E) N(E) w(E),
		double sumE, sumM, sumM2;
		sumE=sumM=sumM2=0;
		// 3a°) accumulate sum_states A(state) w(E(state)) ...
		for (int imeas=0; imeas < samplesBag.size(); imeas++){
			double E = energySamples.getSample(imeas)*VOL; // lattice energy (energySamples stores E/N !)
			double W = mw.getWeight(E);  // w(E(state))

			// Energy Per Spin: (E_ID=0)
			double e = energySamples.getSample(imeas);
			double e2 = e * e;
			sumE += e * W;    
			
			// Magnetization (per spin, as usual): (M_ID=1)
			double M = magnetizationSamples.getSample(imeas);
			double M2 = M * M;
			sumM += M * W;    
			sumM2 += M2 * W;
		}
		
		// 3b°) normalize, then fills moments[i], i.e. for this temperature point :
		double mean,variance;

		//System.out.print("kT="+kT);
		// E:
		mean = sumE/Z;
		//variance = sumE2/Z - sumE/Z * sumE/Z;
		if (!Double.isNaN(mean) && !Double.isInfinite(mean)) meanEnergyVsKTCurve.addXY(kT, mean);
		//System.out.print("\t<E>="+mean);
		
		// M and chi
		mean = sumM/Z;
		if (!Double.isNaN(mean) && !Double.isInfinite(mean)) meanMagVsKTCurve.addXY(kT, mean);
		variance = sumM2/Z - sumM/Z * sumM/Z;
		if (!Double.isNaN(variance) && !Double.isInfinite(variance)) susceptibilityVsKTCurve.addXY(kT, VOL * variance / kT);// chi
		//System.out.print("\t<M>="+mean);
		//System.out.println("\tchi="+(VOL * variance / kT));
		set.updateGUI();
	}
	
	////////////////// PRIVATE METHODS /////////////////////

	/**
	 * Compute Z * exp(-K) for the given temperature w/o flat-histogram, where K is here to avoid overflows.
	 * @param kT temperature point
	 * @param mw Emu(E) reweighter
	 */
	private double computeZ(double kT, MuCaWeighter mw){
		
		// compute Z = sum_states w(E(state)) ; actually what we compute is Z * exp(-K), where K is the normalization
		// constant, since otherwise Z would be too prone to overflowing !
		double Z = 0.0;
		for (int imeas=0; imeas < samplesBag.size(); imeas++){
			Z += mw.getWeight(energySamples.getSample(imeas)*VOL); // Z += w(E) for each state
		}
		return Z;
	}

		

}

