/* 
Loesungsvorschlag fuer 
              :  Uebung 1 
zur Vorlesung :  Moderne numerische Methoden der Festkoerperphysik

(c) Eberhard Jakobi, 2007 
*/

// includes
#include <iostream>
#include <string>
#include <fstream>
#include <unistd.h>
#include <cmath>

class Potential {
	public:
		double left() { return -3.; }
		double right() { return 3.; }
		double operator() (double fx) {
			if (fx<-3. || fx>3.) {
				return std::numeric_limits<double>::max();
			} else {
				return
				-(std::pow(fx,2)+.3)*std::exp(-std::pow((fx-.1),2));
			}
		}
};

bool inline metropolis(double beta, double E_n, double E_c) {
	if (E_n<E_c) return true;
	double r = (double)(rand())/((double)(RAND_MAX));
	if (std::exp(-beta*(E_n-E_c))>r) return true; else return false;
}

class Move {
	private:
		Potential& pot;

	public:
		enum Observable { MeanPosition, MeanEnergy };

		Move(Potential& _pot,size_t _max_iter, Observable _observable):pot(_pot) {
			max_iter = _max_iter;
			observable = _observable;
			warm_up_sweeps = max_iter/100;
			max_iter += warm_up_sweeps;
		}
		double operator() (double fx, std::ostream& x_out, std::ostream& E_out, double s=-1.) {
			size_t iter = 0;
			size_t N = 0;
			beta = 1./fx;
			
			// average
			double E1 = 0.;
			double x1 = 0.;
			
			double x_akt = 0.;
			double E_akt = pot(x_akt);

			size_t accepted = 0;

			if (s<0.) { s = 1.+(fx-0.005)*10.; }

			while(iter++ < max_iter) {
				double r = rand()/((double)(RAND_MAX));

				double x_neu = x_akt + (r-.5)*s;
				double E_neu = pot(x_neu);

				if (metropolis(beta,E_neu,E_akt)) {
					x_akt = x_neu;
					E_akt = E_neu;					
					if (iter>warm_up_sweeps) accepted++;
				}

				if (iter>warm_up_sweeps) {
					N++;
					E1+=E_akt;
					x1+=x_akt;

					E_out << E_akt << "\n";
					x_out << x_akt << "\n";
				}
			}
			// average
			x1 /= (double)N;
			E1 /= (double)N;

			// return observable
			switch(observable) {
				case MeanPosition:
					return x1;
				case MeanEnergy:
					return E1;
				default:
					std::cerr << "ERROR : unknown selected.\n";
					return 0.;
			}
		}
	private:
		double beta;
		double x_c;
		double s;
		size_t max_iter, warm_up_sweeps;
		Observable observable;
};


int main(int argc, char *argv[]) {
	std::cout << "# Metropolis-Tester.\n";
	Potential pot;

	const size_t max_samples = 1000;
	std::ofstream pot_file("Potential.out");
	for (size_t i = 0; i < max_samples; i++) {
		double x = pot.left()+(pot.right()-pot.left())*i/(double)max_samples;
		pot_file << x << "\t" << pot(x) << "\n";
	}
	pot_file.close();
	
	std::cout << "Starting ... \n";
	Move::Observable ob = Move::MeanPosition;
	Move move(pot,20000,ob);

	std::ofstream x1_file("T1.00_s0.50_20k.x");
	std::ofstream E1_file("T1.00_s0.50_20k.E");
	move(1.0,x1_file,E1_file,0.5);
	x1_file.close(); E1_file.close();

	std::ofstream x2_file("T1.00_s0.30_20k.x");
	std::ofstream E2_file("T1.00_s0.30_20k.E");
	move(1.0,x2_file,E2_file,0.3);
	x2_file.close(); E2_file.close();

	std::ofstream x3_file("T1.00_s0.10_20k.x");
	std::ofstream E3_file("T1.00_s0.10_20k.E");
	move(1.0,x3_file,E3_file,0.1);
	x3_file.close(); E3_file.close();

	std::ofstream x4_file("T0.10_s0.10_20k.x");
	std::ofstream E4_file("T0.10_s0.10_20k.E");
	move(0.1,x4_file,E4_file,0.1);
	x4_file.close(); E4_file.close();

	std::ofstream x5_file("T0.003_s0.10_20k.x");
	std::ofstream E5_file("T0.003_s0.10_20k.E");
	move(0.003,x5_file,E5_file,0.1);
	x5_file.close(); E5_file.close();
	
	std::cout << "single runs finished\n";

	std::ofstream T_file("T_log_1000k.x");
	ob = Move::MeanPosition;
	Move move2(pot,2000000,ob);
	std::cout << "Temperature run ... \n";
	for (size_t i = 200; i > 0; i--) {
		double temperature = std::exp(-(double)i*0.05);
		std::ofstream x_file("/dev/null");
		std::ofstream E_file("/dev/null"); 
 		if (i/30 == 0) std::cout << "." << std::flush;
		T_file << temperature << "\t" << move2(temperature,x_file,E_file,-1.0) << "\n";
		x_file.close(); E_file.close();
	}
	T_file.close();

	std::cout << "\n ... finished\n";
}