/*  Monte-Carlo simulation code for statistical physics
    Copyright (C) 2001-2004  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;

/**
 * Returns exp(Emu(E)/kT0 - E/kT - antiOverflowK), where "antiOverflowK"  enables us to avoid overflows.
 *
 * This is based on the better estimate of n(E), which is N(E)/PI_mu(E), where:
 * - PI_mu(E) is exp(-Emu(E)/kT0),
 * - N(E) is the energy histogram.
 *  
 * Hence n(E) ~ N(E) exp(Emu(E)/kT0). Once reweighted at temperature kT,
 * we thus get : N(E) exp(Emu(E)/kT0 - E/kT). Hence the reweighting factor w(E) is given by exp(Emu(E)/kT0 - E/kT).
 */		
public class MuCaWeighter { 
		
	private double kT; 
	private double antiOverflowK;
	private WLDOSHistogram hmue; // store Emu(E)

	/**
	 * Construct a new MuCa weighter from the given S(E)
	 */
	public MuCaWeighter(double kT, WLDOSHistogram hmue){
		this.kT = kT;
		this.hmue = hmue;
		this.antiOverflowK=0;
		computeAntiOverflowConstant();
	}
	
	
	/** 
	 * Overflow work-around : max=Math.exp(709.78) and min=Math.exp(-744.44)
	 * The idea is to first determine max{w(E)} (actually max{Emu(E)/kT0 - E/kT}), that is, log w(E),
	 * and possibly renormalize to avoid overflow 
	 */
	private void computeAntiOverflowConstant(){
		double maxEmuCa = -Double.MAX_VALUE; 
		//double dE = muca.getDeltaE(); [pending]
		double minEnergyFlat = hmue.xAxis().binCentre(0);
		double maxEnergyFlat = hmue.xAxis().binCentre(hmue.xAxis().bins()-1);
		double dE = (maxEnergyFlat-minEnergyFlat)/100.0; 
		for (double ene=minEnergyFlat+dE; ene < maxEnergyFlat; ene+=dE){
			if (getExponent(ene) == Double.POSITIVE_INFINITY) continue; // i.e. if S(E)==REJECT_MOVE
			maxEmuCa = Math.max(maxEmuCa, getExponent(ene));
		}
		//System.out.println("Max Emuca = " + maxEmuCa);
		if (maxEmuCa >= Math.log(Double.MAX_VALUE)-40){
			this.antiOverflowK = maxEmuCa - Math.log(Double.MAX_VALUE) + 40; // e.g. if maxEmuCa=800, antiOverflowK=121
		}
		else this.antiOverflowK=0.0;
		// hence, now there can't be any overflow (yet underflows are still allowed!)
	}
		
	/** 
	 * return S(E) - E/kT ; allow computing max value for overflow workaround 
	 */
	private double getExponent(double energy){
		return this.hmue.getSE(energy)-energy/this.kT;
	}


	//////////////////////////////////////////////////////////////////////////////////////////
	/// IWeighter inteface
	//////////////////////////////////////////////////////////////////////////////////////////
	
	/**
	 * Returns w(E,kT) = exp(S(E) - E/kT - antiOverflowK), except if E > Emax or E < Emin (i.e. flatness boundaries)
	 * in which case 0.0 is returned.
	 */
	public double getWeight(double energy){
		double SE = this.hmue.getSE(energy);
		return Math.exp(SE-energy/this.kT-this.antiOverflowK);
	}

	/** 
	 * return the normalization constant 
	 */
	public double getAntiOverflowK(){
		return antiOverflowK;
	}
}


