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

import fr.ensea.chart.*;
import fr.ensea.math.*;
import fr.ensea.montecarlo.canonical.*;
import fr.ensea.montecarlo.multicanonical.*;
import fr.ensea.montecarlo.data.*;
import java.util.*;

/**
 * This class represents a lattice spin model with nearest-neighbor interactions, and helical boundary
 * conditions. Each spin can take on Q integer values b/w 0 and Q-1 (Potts model). The special value
 * Q=2 gives the Ising model.
 *
 * This class holds the state of each spin, the lattice energy, the number of states of the model,
 * and the number of spins per phase (this is used to compute the order parameter).
 *
 * Helical conditions are implemented by storing spins in a one-dimensional array, from left to right
 * and from top to bottom. Thus the "right" neighbor (resp. "left") of a spin belonging to the last
 * column is the first (resp. last) one of the next row. However, vertical boundary conditions are 
 * implemented in the standard way (PBC's, aka Born-Von-Karman). Helical boundary conditions yield
 * a lighter implementation, and are known to reduce sensitivity to finite lattice size.
 
 * [FR]
 * Cette classe représente un modèle de spins sur réseau avec interactions entre plus proches voisins,
 * et des conditions aux limites hélicoïdales. Chaque spin peut prendre Q valeurs entre 0 et Q-1 (modèle de Potts).
 * On retrouve le modèle d'Ising en imposant Q=2.
 * 
 * Sont stockés dans cette classe : l'état de chaque spin et l'énergie totale du réseau, le nombre d'état Q
 * du modèle, et le nombre de spins par phase, cette dernière variable étant nécessaire au calcul du paramètre d'ordre.
 *
 * Les conditions aux limites hélicoïdales sont implémentées en rangeant l'ensemble des spins du réseau, 
 * de gauche à droite et de haut en bas, dans un tableau unidimensionnel. Ainsi, le voisin "de droite" 
 * (respectivement "de gauche") d'un spin situé sur la dernière colonne est le premier (resp. le dernier) 
 * spin de la ligne suivante (resp. précédente). En revanche, les conditions aux limites verticales sont les
 * conditions périodiques habituelles (à la Born-Van-Karman). L'avantage des conditions aux limites hélicoïdales
 * est simplement d'alléger l'implémentation de l'algorithme.
 */
public class PottsLattice implements ILattice2D {

	// package access for use inside PottsWolffCluster etc...
	public int pottsQ; // nb of Q states

	public int size; // lattice size
	public int size2; // L*L

	public int[] spins; // array holding the state of each spin; helical boundary conditions are used 
	public int[] nbSpinsQ; // nbSpins[s]=nb of spins in phase "s" (from 0 to Q-1)

	public double energy; // lattice energy
	public double externalField; // external magnetic field

	/**
	 * Init a 2D square lattice with L=10 (i.e. 100 spins), with each spin inited at random b/w 0 and 2 (Q=3).
	 */
	public PottsLattice() {
		init(3,10);
	}
	
	/**
	 * Init a 2D square lattice with the given size L, and each spin inited at random b/w 0 and Q-1
	 */
	public PottsLattice(int Q, int L) {
		init(Q,L);
	}

	/**
	 * Reinit lattice properties.
	 * @param Q nb of states of the model
	 * @param L lattice size
	 */
	public synchronized void init(int Q, int L){
		pottsQ = Q;
		size = L;
		size2 = L*L;
		
		spins = new int[size2];
		createRandomConfiguration(); 
		updateEnergy();
	}
	
	/**
	 * Reinit  lattice properties and spin values from the given source.
	 * @param Q nb of states of the model
	 * @param L lattice size
	 */
	public synchronized void init(PottsLattice src){
		pottsQ = src.pottsQ;
		size = src.size;
		size2 = src.size2;
		
		spins = new int[src.size2];
		System.arraycopy(src.spins,0, spins,0, spins.length);
	}

	/**
	 * Creates a random configuration
	 */
	public synchronized void createRandomConfiguration(){
		for(int i=0; i<size2; i++) spins[i] = (int)(Math.random() * pottsQ);
		updateEnergy();
		updateNbSpinsPerPhase();
	}
	
	/**
	 * Creates a ground-state configuration, i.e. all spins in phase s=0.
	 */
	public synchronized void createGroundstateConfiguration(){
		spins = new int[size2];
		updateEnergy();
		updateNbSpinsPerPhase();
	}

	/**
	 * Compute the nb of spins in each phase and store in cache
	 */
	public synchronized void updateNbSpinsPerPhase(){
		nbSpinsQ = new int[pottsQ];
		for (int i=0; i<size2; i++){
			nbSpinsQ[spins[i]]++;
		}
	}
	
	/**
	 * Compute the lattice energy and store in cache
	 */
	public synchronized void updateEnergy(){
		
		// lattice energy:
		
		// n-n interaction:
		//       * - (*)
		//       |
		//      (*)
		//
		// and helical boundary conditions
		energy = 0.0;

		int i;  
		for (i = 0; i < size2; i++) { 
			if (spins[i]==spins[neighborRight(i)]) energy--;
			if (spins[i]==spins[neighborBelow(i)]) energy--;
			if (spins[i]==0) energy -= externalField;
		}
	}
		
	/////////////////////////////////////////////////////////////////
	///// Helical boundary conditions
	/////////////////////////////////////////////////////////////////
	
	/**
	 * Returns the index of the nearest-neighbor left of the given spin
	 */
	public int neighborLeft(int idx){
		int nn=idx-1;
		if (nn < 0) nn += size2;
		return nn;
	}

	/**
	 * Returns the index of the nearest-neighbor right of the given spin
	 */
	public int neighborRight(int idx){
		int nn=idx+1;
		if (nn >= size2) nn -= size2;
		return nn;
	}

	/**
	 * Returns the index of the nearest-neighbor below the given spin
	 */
	public int neighborBelow(int idx){
		int nn=idx-size;
		if (nn < 0) nn += size2;
		return nn;
	}

	/**
	 * Returns the index of the nearest-neighbor above the given spin
	 */
	public int neighborAbove(int idx){
		int nn=idx+size;
		if (nn >= size2) nn -= size2;
		return nn;
	}


	/////////////////////////////////////////////////////////////////
	//// Properties
	/////////////////////////////////////////////////////////////////

	/**
	 * Return the lattice size (=side length)
	 */
	public int getLatticeSize(){
		 return size;
	 }

	/**
	 * Return the nb of spins (=side length squared)
	 */
	public int getSpinCount(){
		 return size2;
	}
	 
	/**
	 * Return the number of states of the model (2 for the Ising model)
	 */
	public int getPottsQ(){
		 return this.pottsQ;
	}
	
	/**
	 * Returns the exact critical temperature
	 */
	public double getExactTc(){
		return 1.0/Math.log(1+Math.sqrt(pottsQ));
	}
	 
	/**
	 * Returns the lattice energy
	 */
	public double getEnergy(){
		return energy;
	}

	/**
	 * Computes then returns the magnetization
	 */
	public double getMagnetization(){
		int maxQ=nbSpinsQ[0]; 
		for (int q=1; q<pottsQ;q++){ // détermine la phase possédant le plus grand nombre de spins
			if (nbSpinsQ[q] > maxQ) maxQ = nbSpinsQ[q];
		}
		double mag = (pottsQ * (double)maxQ / (double)size2 - 1.0)/(pottsQ - 1.0); 
		return mag;
	}
	
	/**
	 * Returns the value of the ith spin (with helical PBC's)
	 */
	public int getSpinValue(int i){
		return spins[i];
	}

	/**
	 * Returns the value of the spin at (i,j)
	 */
	public int getSpinValue(int i, int j){
		return spins[i*size+j];
	}

	/**
	 * Changes the value of the ith spin to val, then updates the nb of spins per phase.
	 */
	public void setSpinValue(int i, int val){
		if (val >= this.pottsQ) val=pottsQ-1;
		nbSpinsQ[spins[i]]--;
		spins[i] = val;
		nbSpinsQ[val]++;
	}
	
	/**
	 * Changes the value of the external magnetic field (forces condensation in phase s=0)
	 */
	 public void setMagneticField(double fieldValue){
		 externalField = fieldValue;
	 }
	
	/**
	 * Returns the current value of the external magnetic field
	 */
	 public double getMagneticField(){
		 return externalField;
	 }

	/**
	 * Returns the energy of the ground-state
	 */
	 public double getGroundstateEnergy(){
		 return -this.size2*2;
	 }

	/////////////////////////////////////////////////////////////////
	//// Renormalization
	/////////////////////////////////////////////////////////////////
	
	/**
	 * returns the power of 2 that is immediately below (or equal to) the given number 
	 */
	public static int floorToPowerOf2(final int x){
		for(int i=1,y=1<<30;i<32;i++,y>>=1) { // generate 2^31,...,4,2,1
			int z=x&y; // isolate one bit, starting from MSB
			if(z!=0) { // check if this bit is non-null => this yields the power of 2 just below x
				return y;
			}
		}
		return x;
	}

	/**
	 * Generates a renormalized lattice using the majority rule and 2x2 plaquettes.
	 * Note that the lattice size must be a power of 2.
	 * @param result renormalized lattice (may be null, in which case it's allocated here)
	 * @return result, for convenience.
	 */ 
	public PottsLattice renormalize(PottsLattice result){
		if (size==1 || (((size>>1)<<1)!=size)) { // can't renormalize if size=1 or not even
			if (result==null) result = new  PottsLattice(pottsQ,size);
			result.init(this);
			return result; 
		}
		if (result==null || result.size != (size >>1)) result = new PottsLattice(pottsQ,size>>1);
		int spinDest=0;
		int spinSrc;
		int valRenorm;
		switch (pottsQ){
			
			// --- Q=2 ---
			case 2:
			// from bottom to top and from left to right:
			// val3 val4
			// val1 val2
			// with val=0 or 1
			// * valRenorm = sum(val1 à val4) = 0 to 4
			// * if valRenorm==2, draw at random
			// * if valRenorm<2, set it to 0
			// * if valRenorm>2, set it to 1
			for (int row=0; row< size; row+=2){
				for(int col=0; col < size; col+=2){
					spinSrc=row*size+col;
					// majority rule:
					valRenorm = this.spins[spinSrc];
					valRenorm += this.spins[neighborRight(spinSrc)];
					valRenorm += this.spins[neighborAbove(spinSrc)];
					valRenorm += this.spins[neighborRight(neighborAbove(spinSrc))];
					if (valRenorm > 2) valRenorm=1;
					else if (valRenorm < 2) valRenorm=0;
					else valRenorm = (int)(2*Math.random());
					result.spins[spinDest] = valRenorm;
					spinDest++;
				}
			}
			break;
			
			// --- Q>2 ---
			default:
			//System.out.println("========== RG:Q>2 ==========");
			// from bottom to top and from left to right:
			// val3 val4
			// val1 val2
			// with val=0 or Q-1
			// * valRenorm = maj(val1,val2,val3,val4), and we draw at random if undetermined.
			int[] nbSpinsQRG = new int[pottsQ];
			for(int row=0; row < size; row+=2){
				for (int col=0; col< size; col+=2){
					//System.out.println("[col,row]=["+col+","+row+"]");
					// reset array
					for (int i=0; i<pottsQ; i++) nbSpinsQRG[i]=0;
					valRenorm=-1; // used also as a flag, see blow in for() loop
					spinSrc=row*size+col;
					// majority rule:
					nbSpinsQRG[this.spins[spinSrc]]++;
					nbSpinsQRG[this.spins[neighborRight(spinSrc)]]++;
					nbSpinsQRG[this.spins[neighborAbove(spinSrc)]]++;
					nbSpinsQRG[this.spins[neighborRight(neighborAbove(spinSrc))]]++;
					int maxSpinCountPerPhase=0; // max nb of spins in each phase 
					int majorityPhase=0; // spin value appearing most often in the plaquette
					for (int i=0; i<pottsQ; i++) { // for each phase
						//System.out.println("nbSpinsQRG["+i+"]="+nbSpinsQRG[i]);
						if (nbSpinsQRG[i] > maxSpinCountPerPhase) {
							maxSpinCountPerPhase=nbSpinsQRG[i];
							majorityPhase=i;
							//System.out.println("\tmaxSpinCountPerPhase="+maxSpinCountPerPhase);
							//System.out.println("\tmajorityPhase="+majorityPhase);
						}
						// if two phases have the same nb of spins, then there are two spins per phase, and we
						// must draw at random NOW which phase we choose (it's more efficient than doing it later)
						else if (maxSpinCountPerPhase==2 && nbSpinsQRG[i] == maxSpinCountPerPhase){
							// the current phase "i" and the phase "majorityPhase" contain 2 spins each:
							int b = (int)(2.0*Math.random()); // 0 or 1
							if (b==0) valRenorm=i;
							else valRenorm=majorityPhase;
							break;
							// thus valRenorm > -1 !
						}
					}
					if (valRenorm==-1){ // "valRenorm" not set yet
						// {X X , X X} or {X X , X Y} or {X X , Z Y}
						// note that {X X , Y Y} was handled in the for() loop just above
						if (maxSpinCountPerPhase > 1) valRenorm=majorityPhase;
						else { // {X Y, Z T} => draw one spin value at random
							int b = (int)(4.0*Math.random()); // 0, 1, 2 or 3
							switch (b){
								case 0: valRenorm = this.spins[spinSrc]; break;
								case 1: valRenorm = this.spins[neighborRight(spinSrc)]; break;
								case 2: valRenorm = this.spins[neighborAbove(spinSrc)]; break;
								case 3: valRenorm = this.spins[neighborRight(neighborAbove(spinSrc))]; break;
								default:
							}
						}
					}
					//System.out.println("valRenorm="+valRenorm);
					result.spins[spinDest] = valRenorm;
					spinDest++;
				}
			}
			break;
		}
		result.updateEnergy();
		result.updateNbSpinsPerPhase();
		//System.out.println("lattice = " + toString());
		//System.out.println("RG lattice = " + result.toString());
		return result;
	}

	/////////////////////////////////////////////////////////
	//// Metropolis
	/////////////////////////////////////////////////////////
	
	/**
	 * Compute the energy difference b/w this lattice and the same lattice with spin "spinIdx" moved to "newValue"
	 */
	public double getEnergyChange(int spinIdx, int newVal){

			double dE=0;
			int oldVal = getSpinValue(spinIdx);
			if (newVal==oldVal) return 0;

			int neighborVal = spins[neighborLeft(spinIdx)];
			if (oldVal==neighborVal) dE++;
			if (newVal==neighborVal) dE--;

			neighborVal = spins[neighborRight(spinIdx)];
			if (oldVal==neighborVal) dE++;
			if (newVal==neighborVal) dE--;

			neighborVal = spins[neighborAbove(spinIdx)];
			if (oldVal==neighborVal) dE++;
			if (newVal==neighborVal) dE--;

			neighborVal = spins[neighborBelow(spinIdx)];
			if (oldVal==neighborVal) dE++;
			if (newVal==neighborVal) dE--;
			
			// external field
			if (oldVal==0) dE += externalField;
			else if (newVal==0) dE -= externalField;

			return dE;
	}

	//////////////////////////////////////
	//// algorithms
	//////////////////////////////////////
	public Montecarlo createMetropolis(SamplesBag bag, double kT){
		return new PottsMetropolis(this, bag, kT);
	}
	
	public Montecarlo createWolffCluster(SamplesBag bag, double kT){
		 return new PottsWolffCluster(this, bag, kT);
	}
	
	public WangLandau createWangLandau(SamplesBag bag, double weight){
		if (this.getPottsQ()==2) return new PottsWL(this,bag,2.0, weight);
		else return new PottsWL(this,bag,1.0, weight);
	}
	
	/////////////////////////////////////////////////////////
	//// debug
	/////////////////////////////////////////////////////////
	public String toString(){

		String s = "Modèle de Potts bidimensionnel  : \n";
		s += "* size de réseau=" + Integer.toString(size) + "\n";
		s += "* M=" + getMagnetization() + "\n";
		s += "* E=" + getEnergy() + "\n";
		// affiche le contenu de chaque ligne du réseau :
		int i=0;
		for (int row = 0; row < size; row++){
			for(int col = 0; col < size; col++){
				s += Integer.toString(spins[i++]) + " ";
			}
			s += "\n";
		}
		return s;
	}
	
	// for testing purpose
	public static void main(String[] args){
		PottsLattice lat = new PottsLattice(4, 2);
		for (int i=0; i<lat.size2; i++) lat.setSpinValue(i,0);
		lat.updateEnergy();
		System.out.println(lat);
	}
	

}

