#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <sstream>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <string>
#include <iterator>

#include "mersenne.h"


using namespace std;

typedef double REAL;
typedef unsigned int uint;
typedef char SPINVAL_T;
typedef vector<SPINVAL_T>::iterator LATTICEITER;

#define ARRAYELEMCOUNT(x) sizeof(x) / sizeof(x[0])
#define TAB "\t"
#define SeedRand(x) seed_mersenne(x)

void ReadParameters(vector<REAL>& betas, uint& Lx, uint& Ly, uint& sweepcount, string& outputfilename, uint& meafreq)
//-----------------------------------
{
	//Paramfile:
	//1. row: beta values
	//2. row: Lx
	//3. row: Ly
	//4. row: sweepcount
	//5. row: outputfilename template(for example with "data", outputfiles will then be like data2.txt).
	//6. row: Measurement frequency.
	ifstream instrmParams("params.txt");
	if(instrmParams.fail())
	{
		cerr << "Opening parameterfile params.txt failed" << endl;
		exit(0);
	}
	string temp;
	getline(instrmParams, temp);
	istringstream strstrm(temp);
	betas.clear();
	copy(istream_iterator<REAL>(strstrm), istream_iterator<REAL>(), back_inserter(betas));
	sort(betas.begin(), betas.end());
	instrmParams >> Lx >> Ly >> sweepcount >> outputfilename >> meafreq;
}

void PrintFileinfo(ostream& infostrm, REAL beta, const string& filename)
//----------------------------------------------------------------------
{
	time_t t;
	time(&t);
	infostrm << "#Time     : " << ctime(&t);
	infostrm << "#beta     : " << beta << endl;
	infostrm << "#filename : " << filename << endl << endl;
}

inline REAL CalculateLocalS(SPINVAL_T val, SPINVAL_T nbstates[4], const REAL& beta)
//---------------------------------------------
{
	char sum = 0;
	for(char i = 0; i<4; i++) {if(val == nbstates[i]) sum++;}
	
	return beta*(4-sum);
}

inline SPINVAL_T GenerateTrial(SPINVAL_T val)
//---------------------------------
{
	//Choose new value with even probability from the two 'free' choices. For example when spinvalue is 0,
	//probability to choose value 1 is 0.5 and for value 2 the same.
	return (mersenne() < 0.5) ? ((val+2) % 3) : ((val+1) % 3);
}


inline REAL GetAbsM(const REAL M[3], const REAL& V)
//--------------------------------
{
	return 1.5 * (*max_element(M, M+3)/V - 1.0/3.0);
}

//Create lattice and generate neighbour array for single index array with periodic boundaries. 
//For element i in lattice, nbarray elements 4*i, 4*i+1, 4*i+2, 4*i+3 tell right, down, left, up neighbours.
void InitArrays(vector<SPINVAL_T>& lattice, vector<SPINVAL_T*>& nbarray, const int Lx, const int Ly)
//------------------------------------------------------------------------------------------------------
{
	lattice.resize(Lx*Ly, 0); //Initialize all elements to state 0.
	nbarray.resize(4 * lattice.size(), 0);
	vector<SPINVAL_T*>::iterator iter = nbarray.begin();
	const uint N = lattice.size();
	for(int i = 0; i<N; i++)
	{
		//Right neighbour
		*iter = ((i+1) % Lx > 0) ? &lattice[i+1] : &lattice[i-Lx+1];
		iter++;
		//Lower neighbour
		*iter = (i-Lx >= 0) ? &lattice[i-Lx] : &lattice[Lx*(Ly-1) + i];
		iter++;
		//Left neighbour
		*iter = (i%Lx - 1 >= 0) ? &lattice[i-1] : &lattice[i+Lx-1];
		iter++;
		//Upper neighbour
		*iter = (i+Lx < N) ? &lattice[i+Lx] : &lattice[i%Lx];
		iter++;
	}
}

//Retrieves neighbour states from neighbour array.
inline void GetNeighbourStates(SPINVAL_T statevals[4], uint elemIndex, const vector<SPINVAL_T*>& nbarray)
//-----------------------------------------------------------------------------------------------
{
	vector<SPINVAL_T*>::const_iterator iter = nbarray.begin() + elemIndex*4;
	statevals[0] = **iter; iter++;
	statevals[1] = **iter; iter++;
	statevals[2] = **iter; iter++;
	statevals[3] = **iter;
}

//Helper for printing vector.
ostream& operator<<(ostream& ostrm, vector<REAL>& v)
//--------------------------------------------------
{
	copy(v.begin(), v.end(), ostream_iterator<REAL>(ostrm, " "));
	return ostrm;
}

//Helper for converting number to string.
template<class T> inline std::string ToString(T& x)
//--------------------------------------------------
{
   std::ostringstream o;
   if(!(o << x)) return "FAILURE";
   else return o.str();
}

int main()
//-------
{
	vector<REAL> betas;
	uint Lx, Ly, sweepcount, measurementfreq;
	string outputfilenametemplate;

	ReadParameters(betas, Lx, Ly, sweepcount, outputfilenametemplate, measurementfreq);

	const REAL V = Lx*Ly;

	vector<SPINVAL_T> lattice;
	vector<SPINVAL_T*> nbarray;
	InitArrays(lattice, nbarray, Lx, Ly);

	const int seed = time(0);
	SeedRand(seed);

	ofstream infostrm((outputfilenametemplate+"Detail.txt").c_str());
	if(infostrm.fail())
	{
		cerr << "Could not open datadetail file" << endl;
		exit(0);
	}
	infostrm << "#Seed : " << seed << endl;
	infostrm << "#Betas: " << betas << endl;
	infostrm << "#Lx   : " << Lx << endl;
	infostrm << "#Ly   : " << Ly << endl;
	infostrm << "#Sweeps: " << sweepcount << endl;

	
	//For different beta values
	for(uint i = 0; i<betas.size(); i++)
	{
		//Printing data to different files for each beta.
		const string filename = outputfilenametemplate + ToString(i) + ".txt";
		ofstream datastrm(filename.c_str(), ios::trunc);
		datastrm.precision(12);

		PrintFileinfo(infostrm, betas[i], filename);

		int t0 = time(0);
		cout << "Beta " << betas[i] << " (" << i+1 << "/" << betas.size() << "). ";
		cout << "Output file: " << filename << endl;
		for(uint s = 0; s < sweepcount; s++)
		{
			SPINVAL_T nbStates[4] = {0,0,0,0};
			REAL M[3] = {0,0,0};
			uint index = 0;
			const REAL beta = betas[i];
			for(LATTICEITER iter = lattice.begin(); iter != lattice.end(); iter++, index++)
			{
				GetNeighbourStates(nbStates, index, nbarray);
				REAL H1 = CalculateLocalS(*iter, nbStates, beta);
				SPINVAL_T trialval = GenerateTrial(*iter);
				REAL H2 = CalculateLocalS(trialval, nbStates, beta);
				REAL dH = H2 - H1;
				if(dH <= 0 || exp(-dH) > mersenne())
				{
					*iter = trialval;
				}
				M[*iter]++;
			}
			if((s % measurementfreq) == 0) datastrm << GetAbsM(M, V) << endl;
		}
		t0 = time(0) - t0;
		cout << "Loop lasted " << t0 << " seconds; " << double(t0) / sweepcount << " seconds per sweep" << endl;

		datastrm.close();
	}
	
	return 0;
}