ising.cpp

/*

ising.cpp ... Ising Model
Copyright (C) 2010  KennyTM~ <kennytm@gmail.com>

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.

*/

/*

Please compile with
g++-mp-4.5 -O3 -march=native -Wall -std=c++0x -I/opt/local/include -lgsl -fno-unswitch-loops -fno-tree-vectorize -funsafe-loop-optimizations -fno-math-errno ising.cpp

To sample, please use
shark -1 -i ./a.out

*/

#include <vector>
#include <utility>
#include <gsl/gsl_rng.h>
#include <cassert>
#include <cstdio>
#include <cstdint>
#include <algorithm>
#include <ctime>
#include <iterator>
#include <cmath>
#include <numeric>

using namespace std;

static gsl_rng* rng;

typedef int64_t Energy;
typedef int Node;
typedef int64_t EnergySquared;
typedef double Beta; // Beta -> inverse temperature
enum Spin : char {
	spinUp = 1,
	spinDown = -1
};
typedef int TotalSpin;

struct Neighbor {
	Node node;
	Energy energy;
};

struct Edge {
	Node source, target;
	Energy energy;
};

struct Ising;

struct ExpCache {
	Energy e;
	unsigned long expval;
};
static vector<ExpCache> cache;

__attribute__((hot))
bool cmp_memoized_exp(Beta beta, Energy x) {
	unsigned long res;

	for (auto cit = cache.cbegin(); cit != cache.cend(); ++ cit)
		if (cit->e == x) {
			res = cit->expval;
			goto cmp;
		}
	
	res = (unsigned long)(exp(-beta * x) * 0x100000000L);
	cache.push_back(ExpCache{x, res});
	
cmp:
	return gsl_rng_get(rng) < res;
}
void flush_memoized_exp() {
	cache.clear();
}


struct Graph {
	// Construct a graph of given size and edges.
	Graph(const int size, const vector<Edge>& edges) : adjs(size) {
		for (auto ait = adjs.begin(); ait != adjs.end(); ++ ait)
			ait->reserve(4);
		for (auto cit = edges.cbegin(); cit != edges.cend(); ++ cit) {
			adjs[cit->source].push_back(Neighbor{cit->target, cit->energy});
			adjs[cit->target].push_back(Neighbor{cit->source, cit->energy});
		}
	}
	
	// Return a square lattice.
	static Graph square_lattice (const int width, const int height) {
		vector<Edge> edges;
		edges.reserve(2*width*height);
		Node i = 0;
		for (int y = 0; y < height; ++ y) {
			for (int x = 0; x < width; ++ x, ++ i) {
				Node right = (x == width-1) ? i-width+1 : i+1;
				Node down = (y == height-1) ? x : i+width;
				edges.push_back(Edge{i, right, 1});
				edges.push_back(Edge{i, down, 1});
			}
		}
		return Graph(width*height, edges);
	}
	
	// Return a hexagonal lattice with impurity.
	//  - impurity_level must be between 0 and 2^32-1.
	//  - width & height must be even numbers.
	static Graph hexagonal_lattice (const int width, const int height, const unsigned long impurity_level, const Energy pure_energy, const Energy impure_energy) {
		assert(width % 2 == 0);
		assert(height % 2 == 0);
		
		vector<Edge> edges;
		edges.reserve(3*width*height/2);
		Node i = 0;
		for (int y = 0; y < height; ++ y) {
			for (int x = 0; x < width; ++ x, ++ i) {
				Node right = (x == width-1) ? i-width+1 : i+1;
				edges.push_back(Edge{i, right, gsl_rng_get(rng) < impurity_level ? impure_energy : pure_energy});
				
				if ((x + y) % 2 == 0) {
					Node down = (y == height-1) ? x : i+width;
					edges.push_back(Edge{i, down, gsl_rng_get(rng) < impurity_level ? impure_energy : pure_energy});
				}
			}
		}
		return Graph(width*height, edges);
	}
	
	int size() const { return adjs.size(); }
	
private:
	vector<vector<Neighbor>> adjs;
	
public:
	void debug() const {
		Node i = 0;
		for (auto ait = adjs.cbegin(); ait != adjs.cend(); ++ ait) {
			printf("%d: ", i++);
			for (auto eit = ait->cbegin(); eit != ait->cend(); ++ eit) {
				printf("(%d, %lld), ", eit->node, eit->energy);
			}
			printf("\n");
		}
	}
	
	friend struct Ising;
//	friend int main();
};

struct Ising {	
	Ising(const Graph& g) : graph(g), spins(g.size(), spinUp) {}
	
	int size() const { return spins.size(); }
	
	Energy energy() const {
		Energy total_energy = 0;
		Node i = 0;
		auto sit = spins.cbegin();
		auto ait = graph.adjs.cbegin();
		for (; ait != graph.adjs.cend(); ++ ait, ++ sit, ++ i) {
			Spin spin = *sit;
			for (auto eit = ait->cbegin(); eit != ait->cend(); ++ eit) {
				if (spin == spins[eit->node]) 
					total_energy -= eit->energy;
				else
					total_energy += eit->energy;
			}
		}
		return total_energy / 2;
	}

	TotalSpin total_spin() const {
		TotalSpin total_spin_1 = 0;
		for (auto sit = spins.cbegin(); sit != spins.cend(); ++ sit)
			total_spin_1 += *sit;
		return total_spin_1;
	}
		
	Energy energy_change_by_flipping(const Node node) const {
		Spin spin = spins[node];
		const vector<Neighbor>& adj = graph.adjs[node];
		Energy total_energy_change = 0;
		for (auto eit = adj.cbegin(); eit != adj.cend(); ++ eit) {
			if (spin == spins[eit->node])
				total_energy_change += eit->energy;
			else
				total_energy_change -= eit->energy;
		}
		return total_energy_change * 2;
	}
	
	void flip(const Node node) {
		spins[node] = spins[node] == spinUp ? spinDown : spinUp;
	}
	
	__attribute__((hot))
	void monte_carlo_step(const Beta beta, Energy* delta_energy, TotalSpin* delta_total_spin) {
		Node node = gsl_rng_get(rng) % size();	// warning! requires size() is a power of 2.
		Energy dE = energy_change_by_flipping(node);
		
		if (dE < 0 || cmp_memoized_exp(beta, dE)) {
			flip(node);
			if (delta_energy)
				*delta_energy = dE;
			if (delta_total_spin)
				*delta_total_spin = spins[node]*2;
		} else {
			if (delta_energy)
				*delta_energy = 0;
			if (delta_total_spin)
				*delta_total_spin = 0;
		}
	}
	
	void run_monte_carlo(const Beta beta, int stats_steps, Energy* energy_history, TotalSpin* spin_history) {
		int steps = 20000 * size();
		
		// dry run for 20000N - 180000 steps.
		for (int i = stats_steps; i < steps; ++ i)
			monte_carlo_step(beta, NULL, NULL);
		
		// collect stats for the last 180000 steps;
		Energy cur_energy = energy();
		TotalSpin cur_spin = total_spin();
		for (int i = 0; i < stats_steps; ++ i) {
			Energy dE;
			TotalSpin dM;
			monte_carlo_step(beta, &dE, &dM);
			cur_energy += dE;
			cur_spin += dM;
			energy_history[i] = cur_energy;
			spin_history[i] = abs(cur_spin);
		}
	}
	
	void reset_spins() { fill(spins.begin(), spins.end(), spinUp); }
	
private:
	Graph graph;
	vector<Spin> spins;
};


template <typename T>
T sqsum(const T res, const T val) { return res + val * val; }

template <typename Iter>
double mean(const Iter beg, const Iter end) {
	auto sum = accumulate(beg, end, 0);
	auto size = distance(beg, end);
	return (double)sum / size;
}

template <typename Iter>
pair<double, double> stats(const Iter beg, const Iter end) {
	typedef typename iterator_traits<Iter>::value_type T;
	T sum = accumulate(beg, end, (T)0);
	T sq_sum = accumulate(beg, end, (T)0, sqsum<T>);
	auto size = distance(beg, end);
	double xmean = (double)sum / size;
	double var = (double)(size * sq_sum - sum * sum) / (size * (size-1));
	return make_pair(xmean, var);
}

int main (int argc, char* argv[]) {
	rng = gsl_rng_alloc(gsl_rng_taus2);
	gsl_rng_set(rng, time(NULL));
	
	auto level = atol(argv[1]) * (0x100000000L / 10);
	
	Graph g = Graph::hexagonal_lattice(8, 16, level, 1, -1000);
	
#if 1
		
	static const long stats_steps = 180000;
	auto energy_history = new Energy[stats_steps];
	auto spin_history = new TotalSpin[stats_steps];
	
	static const long rep_steps = 48;
	auto heat_caps = new double[rep_steps];
	auto mags = new double[rep_steps];
		
	printf("T\tCv\tdCv\tM\tdM\n");
	
	for (int T = 1; T <= 100; ++ T) {
		double beta = 20.0 / T;		

		int j;
		for (j = 0; j < rep_steps; ++ j)  {
			Graph g = Graph::hexagonal_lattice(8, 16, level, -1, 1000);
//			Graph g = Graph::square_lattice(5,5);
			Ising ising (g);

			ising.run_monte_carlo(beta, stats_steps, energy_history, spin_history);
			
			auto energy_p = stats(energy_history, energy_history+stats_steps);
			double heat_capacity = energy_p.second * beta * beta / ising.size();
			double magnetization = mean(spin_history, spin_history+stats_steps) / ising.size();
			
			heat_caps[j] = heat_capacity;
			mags[j] = magnetization;
			
//			if (j == 2 && heat_caps[j] == 0) {
//				j ++;
//				break;
//			}
		}
		
//		auto heat_cap_p = stats(heat_caps, heat_caps+j);
//		auto mags_p = stats(mags, mags+j);

		for (--j; j > 0; --j)
			printf("%g\t%g\t%g\n", 1/beta, heat_caps[j], mags[j]);
					
//		printf("%g\t%g\t%g\t%g\t%g\n", 1/beta, heat_cap_p.first, sqrt(heat_cap_p.second/j), mags_p.first, sqrt(mags_p.second/j));

		flush_memoized_exp();
	}
		
	delete[] spin_history;
	delete[] energy_history;
	delete[] heat_caps;
	delete[] mags;
	
#endif
	
	gsl_rng_free(rng);
	
	return 0;
}