/*  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.canonical;

import fr.ensea.montecarlo.data.*;
import fr.ensea.montecarlo.model.*;
import fr.ensea.montecarlo.canonical.*;
import fr.ensea.montecarlo.misc.*;

import java.io.*;
import java.awt.Color;
import java.text.*;
import java.util.Locale;


/**
 * Single-cluster Wolff algorithm for the 2D Potts model.
 */
public class PottsWolffCluster extends Montecarlo {

	private PottsLattice lattice;
	
	// --- RG ---
	private PottsLattice[] latticeRG; // non-null if RG is active
	private int factorRG;

	private Stack fifo;
	private static final int BUFFER_EMPTY=-1; // used by Stack
	private double bondProba; // 1 - exp(-beta)
	private int bonds; // nb of active bonds for the last flipped cluster
	
	private long totActiveBondsCumul; // incremented by "bonds" at each invokation of sweep()
	private int mcstepCounter; // incremented at each invokation of sweep()
	
	private int newQState, oldQState; // used by "tryNewBond()"
	
	private NumberFormat nfPercent;

	/**
	 * Initialize a new cluster algorithm.
	 */
	public PottsWolffCluster(PottsLattice lattice, SamplesBag bag, double kT) {
		super(bag, kT);
		this.lattice = lattice;
		this.fifo = new Stack();
		nfPercent = NumberFormat.getPercentInstance(Locale.US);
		nfPercent.setGroupingUsed(false);
		nfPercent.setMaximumFractionDigits(0);
	}

	public void activateRG(){
		setFactorRG(1);
	}

	public void setFactorRG(int f){
		if (f<1) f=1;
		this.factorRG=f;
		latticeRG = new PottsLattice[factorRG];
		latticeRG[factorRG-1] = lattice.renormalize(latticeRG[factorRG-1]);
		for (int s=factorRG-2; s >= 0; s--){
			latticeRG[s] = latticeRG[s+1].renormalize(latticeRG[s]);
		}
	}
	
	public PottsLattice getRGLattice(){
		return latticeRG[0];
	}
	
	public synchronized void setTemperature(double kT){
		this.kT = kT;
		this.bondProba = 1.0 - Math.exp(-1.0/kT);
		reinitCounters();
	}
	
	/**
	 * Return a freshly generated random number b/w 0 and 1.
	 */
	public double createRandomNumber(){
		return Math.random(); // [pending] use MersenneTwister
	}

	/**
	 * Build then "flips" a cluster of spins
	 */
	public void sweep(){
		final int qmax = lattice.pottsQ-1;

		mcstepCounter++;
		int spinIdx;

		// *) pick a spin at random
		spinIdx = (int)(lattice.size2 * createRandomNumber()); // b/w 0 and L^2-1 (raw() => slow)

		// *) save old q-state
		this.oldQState = lattice.getSpinValue(spinIdx);

		// *) pick a new spin state
		this.newQState = lattice.getSpinValue(spinIdx) + 1 + (int)(createRandomNumber() * qmax);
		if (this.newQState > qmax) this.newQState -= lattice.pottsQ; // assurons-nous que la nouvelle valeur est dans [0,Q-1]

		// *) flip the cluster initial seed:
		lattice.spins[spinIdx] = newQState;
		this.bonds=0; 
		fifo.reset();
		fifo.push(spinIdx);

		while((spinIdx=fifo.pop()) != BUFFER_EMPTY){ // -1
			
			tryNewBond(lattice.neighborLeft(spinIdx));
			tryNewBond(lattice.neighborRight(spinIdx));
			tryNewBond(lattice.neighborAbove(spinIdx));
			tryNewBond(lattice.neighborBelow(spinIdx));
		}
		int flippedSpins = bonds+1;
		lattice.nbSpinsQ[oldQState] -= flippedSpins;
		lattice.nbSpinsQ[newQState] += flippedSpins;
		lattice.updateEnergy();
		this.totActiveBondsCumul += bonds;

		samplesBag.addEnergySample(lattice.getEnergy()/lattice.size2);
		samplesBag.addMagnetizationSample(lattice.getMagnetization());
		if (latticeRG != null){ // means RG has been activated
			latticeRG[factorRG-1] = lattice.renormalize(latticeRG[factorRG-1]);
			for (int s=factorRG-2; s >= 0; s--){
				latticeRG[s] = latticeRG[s+1].renormalize(latticeRG[s]);
			}
			samplesBag.addEnergyRGSample(latticeRG[0].getEnergy()/latticeRG[0].size2);
		}
	}
	

	/** try and add a new spin to the cluster with probability "bondProba" if spins are identical */
	private void tryNewBond(int nn){
		if (oldQState==lattice.getSpinValue(nn) && createRandomNumber() < bondProba) { // oldQState = initial value of cluster seed
			lattice.spins[nn]=newQState;  // flip spin
			bonds++; // increment B = number of active bonds
			fifo.push(nn);
		}
	}
	
	
	/**
	 * Implementation of a FIFO buffer used to stores indices of spins belonging to the cluster 
	 */
	class Stack {
		private int[] fifoBuffer; // FIFO buffer for cluster construction
		private int fifoPointer; // index of last added spin index ; FIFO grows upwards

		Stack(){
			fifoBuffer = new int[400]; // 20x20
			reset();
		}

		void reset(){
			fifoPointer=-1;
		}

		void push(int i){
			if (fifoPointer==fifoBuffer.length-1){ // stack overflow !
				int[] _fifoBuffer = new int[fifoBuffer.length*2];
				System.arraycopy(fifoBuffer,0,_fifoBuffer,0,fifoBuffer.length);
				fifoBuffer = _fifoBuffer;
			}
			fifoBuffer[++fifoPointer]=i;
		}

		/** returns BUFFER_EMPTY whenever buffer is empty */
		int pop(){
			if (fifoPointer < 0) return BUFFER_EMPTY; // -1
			else {
				int i = fifoBuffer[fifoPointer];
				fifoPointer--;
				return i;
			}
		}
		
		public String toString(){
			StringBuffer s = new StringBuffer();
			s.append('[');
			for (int i=0; i<=fifoPointer; i++){
				s.append(fifoBuffer[i]);
				s.append(' ');
			}
			s.append(']');
			return s.toString();
		}
	}

	///////////////////////////////////////////////////////////
	/// Statistics
	///////////////////////////////////////////////////////////

	public double getFractionActiveBonds(){
		return (double)totActiveBondsCumul / (double)(mcstepCounter * lattice.size2 * 2);
	}
	
	public synchronized void reinitCounters(){
		this.totActiveBondsCumul = this.mcstepCounter = 0;
	}
	
	public String getCountersRateLabel(){
		return "Active bonds="+nfPercent.format(getFractionActiveBonds());
	}		
	
	/////////////////////////////////////////////////////////////////////////////////
	// TEST
	/////////////////////////////////////////////////////////////////////////////////

	public String toString(){
		String s = "MC Wolff algorithm :\n";
		s += "\tmean acceptance rate = " + nfPercent.format(getFractionActiveBonds()) + "\n";
		return s;
	}

	public static void main(String args[]){

		int L=10, Q=3;
		PottsLattice lattice = new PottsLattice(L,Q);

		SamplesBag bag = new SamplesBag(0);
		PottsWolffCluster algo = new PottsWolffCluster(lattice, bag, 0.5);
		for (double kT=0.5; kT<2.0; kT+=0.2){
			System.out.println("---------- kT="+kT+" ----------");
			algo.setTemperature(kT);

			final int nbPas = 10000;

			for (int i=0; i<1000; i++) algo.sweep();
			for (int i=0; i<nbPas; i++) {
				algo.sweep();
			}

			double meanE = bag.computeEnergyMean();
			double varE = bag.computeEnergyVariance(meanE);
			double meanM = bag.computeMagnetizationMean();
			double varM = bag.computeMagnetizationVariance(meanM);
			System.out.print("<E>="+meanE);
			System.out.println("\t\tCv="+(varE/algo.kT/algo.kT/L));
			System.out.print("<M>="+meanM);
			System.out.println("\tChi="+(L*varM/algo.kT));
			System.out.println("Fraction of active bonds =" + (100.0*algo.getFractionActiveBonds())+" %");
		}
	}
}

