Merge branch 'master' of https://gitlab.ethz.ch/periv/cqpfs2019

Exercises/ex06/Solutions/ising2D.py

+import numpy as np
+from timeit import default_timer as timer
+import random
+import matplotlib.pyplot as plt
+
+#does the next level of binning
+def binning_step(data):
+	newdata = []
+	for i in range(int(len(data)/2)):
+		newdata.append((data[2*i]+data[2*i+1])/2.)
+	stdev = np.std(newdata)/np.sqrt(len(newdata)-1)
+	return [stdev,newdata]
+
+#perform a binning analysis of the data, up to the level where at least minsamples samples remain
+def binning(data, minsamples=8):
+	levels = []
+	stdevs = []
+	stdevs.append(np.std(data)/np.sqrt(len(data)-1))
+	levels.append(0)
+	level = 1
+	while(len(data) >= 2*minsamples):
+		stdev,data = binning_step(data)
+		stdevs.append(stdev)
+		levels.append(level)
+		level += 1
+	return [levels, stdevs]
+
+
+#perform a Wolff cluster update step on sys of size H*L
+#absolute value of return is the cluster size, twice the return is the change in magnetisation
+def Wolff_step(sys,H,L,probability):
+	flipped = 1
+	#choose random spin
+	spin_v = random.randint(0,H-1)
+	spin_h = random.randint(0,L-1)
+	#flip the spin
+	sys[spin_v][spin_h] *= -1
+	#keep track of what we do to the magnetisation
+	flip_sign = sys[spin_v][spin_h]
+	#check neighbours recursively
+	spins_to_check = [[spin_v,spin_h]]
+	while(len(spins_to_check) > 0):
+		indices = spins_to_check.pop()
+		spin_v = indices[0]
+		spin_h = indices[1]
+		current_spin = sys[spin_v][spin_h]
+		#add neighbours to cluster and schedule them with correct probability
+		if(current_spin * sys[(spin_v+1)%H][spin_h] == -1 and random.random() < probability):
+			spins_to_check.append([(spin_v+1)%H,spin_h])
+			sys[(spin_v+1)%H][spin_h] *= -1
+			flipped += 1
+
+		if(current_spin * sys[(spin_v-1)%H][spin_h] == -1 and random.random() < probability):
+			spins_to_check.append([(spin_v-1)%H,spin_h])
+			sys[(spin_v-1)%H][spin_h] *= -1
+			flipped += 1
+
+		if(current_spin * sys[spin_v][(spin_h+1)%L] == -1 and random.random() < probability):
+			spins_to_check.append([spin_v,(spin_h+1)%L])
+			sys[spin_v][(spin_h+1)%L] *= -1
+			flipped += 1
+
+		if(current_spin * sys[spin_v][(spin_h-1)%L] == -1 and random.random() < probability):
+			spins_to_check.append([spin_v,(spin_h-1)%L])
+			sys[spin_v][(spin_h-1)%L] *= -1
+			flipped += 1
+	return flipped*flip_sign
+
+
+	
+#perform one single spin update on sys of extent H*L. 
+#returns change in magnetisation
+def MCMC_step(sys,H,L,exptable):
+	#pick random spin
+	spin_v = random.randint(0,H-1)
+	spin_h = random.randint(0,L-1)
+	#calculate local field
+	h = sys[spin_v][(spin_h-1)%L] + sys[spin_v][(spin_h+1)%L] + sys[(spin_v+1)%H][spin_h]  + sys[(spin_v-1)%H][spin_h] 
+	
+	#flip with metropolis probability
+	if(random.random() < exptable[sys[spin_v][spin_h]*h]):
+		if(sys[spin_v][spin_h] == 1):
+			sys[spin_v][spin_h] = -1
+			return -2
+		else:
+			sys[spin_v][spin_h] = 1
+			return 2
+	else:
+		return 0
+
+#returns the magnetisation of sys
+def magnetisation(sys):
+	magn = 0
+	L = len(sys[0])
+	H = len(sys)
+	for a in sys:
+		for s in a:
+			magn += s
+	return float(magn)/float(H*L)
+
+
+#run a simulation using cluster updates
+def run_cluster_simulation(dT,Tmin=1.8,Tmax=3.2,meas=8192*2,equil=1000*2):
+
+	#######################
+	###	Initialisation
+	#######################
+
+	#arrays to store info about correlations
+	statistics_magnetisation = []
+	statistics_magnetisationq = []
+
+	#arrays to store the expectation values and their standard deviations for all temperatures
+	magnetisations = []
+	magnetisationsq = []
+	magnetisations_std = []
+	magnetisationsq_std = []
+	temperatures = []
+
+	#define system parameters
+	L = 16
+	H = 16
+	J = 1.
+	#initialise the system
+	sys = []
+	for i in range(H):
+		a = []
+		for j in range(L):
+			a.append(1)
+		sys.append(a)
+	#initialise the magnetisation
+	total_magnetisation = L*H
+
+	T = Tmin
+
+	
+	#######################
+	###	Simulation
+	#######################
+	#simulate for the specified temperatures
+	while(T <= Tmax):
+		#arrays to store all samples
+		magdata = [] #magnetisation
+		magqdata = [] #magnetisation squared
+		T += dT
+		beta = 1./T
+		probability = 1. - np.exp(-2.*beta) #probability of accepting the addition of a site to the cluster
+
+		#equilibrate
+		for i in range(int(equil)):
+			total_magnetisation += 2*Wolff_step(sys,H,L,probability)
+		#measure
+		for i in range(int(meas)):
+			wolff_return = Wolff_step(sys,H,L,probability)
+			m = np.abs(wolff_return)/float(L*H)
+			total_magnetisation += 2*wolff_return
+			magdata.append(m)
+			magqdata.append((total_magnetisation/float(L*H))**2)
+
+		#extract means and their standard deviations
+		magnetisations.append(np.mean(magdata))
+		magnetisationsq.append(np.mean(magqdata))
+		magnetisations_std.append(np.std(magdata))
+		magnetisationsq_std.append(np.std(magqdata))
+
+		#save temperature
+		temperatures.append(T)
+
+		#do a binning analysis of the data and store it
+		statistics_magnetisation.append(binning(magdata))
+		statistics_magnetisationq.append(binning(magqdata))
+
+	#calculate correlation times for all temperatures from the binning results. binning level 5 seems reasonable for most temperatures
+	taus = []
+	tausq = []
+	for i in range(len(temperatures)):
+		taus.append(0.5*(statistics_magnetisation[i][1][5]**2/magnetisations_std[i]**2*(meas-1.)-1.))
+		tausq.append(0.5*(statistics_magnetisationq[i][1][5]**2/magnetisationsq_std[i]**2*(meas-1.)-1.))
+
+
+
+	#######################
+	###	Plotting
+	#######################
+
+	#correlation time magnetisation
+	plt.plot(temperatures,taus)
+	plt.title("Correlation time Wolff, magnetization")
+	plt.xlabel("Temperature")
+	plt.ylabel("Correlation time")
+	#plt.show()
+	plt.savefig("wolff_corr_times_m.pdf")
+	plt.clf()
+
+	#correlation time magnetisation squared
+	plt.plot(temperatures,taus)
+	plt.title("Correlation time Wolff, magnetization$^2$")
+	plt.xlabel("Temperature")
+	plt.ylabel("Correlation time")
+	#plt.show()
+	plt.savefig("wolff_corr_times_mq.pdf")
+	plt.clf()
+
+	#binning analysis results
+	for i in range(len(statistics_magnetisation)):
+		plt.plot(statistics_magnetisation[i][0],statistics_magnetisation[i][1])
+	leg = []
+	for i in range(len(temperatures)):
+		leg.append("T=" + str(round(temperatures[i],3)))
+	plt.legend(leg)
+	plt.title("Binning analysis Wolff, magnetisation")
+	plt.xlabel("Binning Level")
+	plt.ylabel("Error")
+	#plt.show()
+	plt.savefig("wolff_binning_m.pdf")
+	plt.clf()
+
+	for i in range(len(statistics_magnetisationq)):
+		plt.plot(statistics_magnetisationq[i][0],statistics_magnetisationq[i][1])
+	leg = []
+	for i in range(len(temperatures)):
+		leg.append("T=" + str(round(temperatures[i],3)))
+	plt.legend(leg)
+	plt.title("Binning analysis Wolff, magnetisation$^2$")
+	plt.xlabel("Binning Level")
+	plt.ylabel("Error")
+	#plt.show()
+	plt.savefig("wolff_binning_mq.pdf")
+	plt.clf()
+
+	#calculate errorbars
+	#binning level 5 seems reasonable for most temperatures
+	errorbars_m = []
+	errorbars_mq = []
+	for i in range(len(temperatures)):
+		errorbars_m.append(statistics_magnetisation[i][1][5])
+		errorbars_mq.append(statistics_magnetisationq[i][1][5])
+
+	#plot magnetisation with errorbars
+	plt.errorbar(temperatures,magnetisations,yerr=errorbars_m,fmt='o',ms=1)
+	plt.title("Magnetisation, Wolff")
+	plt.xlabel("Temperature")
+	plt.ylabel("Magnetisation")
+	#plt.show()
+	plt.savefig("wolff_m.pdf")
+	plt.clf()
+
+	#plot magnetisation squared with errorbars
+	plt.errorbar(temperatures,magnetisationsq,yerr=errorbars_mq,fmt='o',ms=1)
+	plt.title("Magnetisation$^2$, Wolff")
+	plt.xlabel("Temperature")
+	plt.ylabel("Magnetisation^2")
+	#plt.show()
+	plt.savefig("wolff_mq.pdf")
+	plt.clf()
+
+
+
+#run a simulation using single spin updates
+def run_metropolis_simulation(dT,Tmin=1.8,Tmax=3.2,meas=131072*4,equil=8192*2):
+
+	#######################
+	###	Initialisation
+	#######################
+
+	#arrays to store info about correlations
+	statistics_magnetisation = []
+	statistics_magnetisationq = []
+
+	#arrays to store the expectation values and their standard deviations for all temperatures
+	magnetisations = []
+	magnetisationsq = []
+	magnetisations_std = []
+	magnetisationsq_std = []
+	temperatures = []
+
+	#define system parameters
+	L = 16
+	H = 16
+	J = 1.
+	#initialise the system
+	sys = []
+	for i in range(H):
+		a = []
+		for j in range(L):
+			a.append(1)
+		sys.append(a)
+	#initialise the magnetisation
+	total_magnetisation = L*H
+
+	T = Tmin
+
+	
+	#######################
+	###	Simulation
+	#######################
+	#simulate for the specified temperatures
+	while(T <= Tmax):
+		#arrays to store all samples
+		magdata = [] #magnetisation
+		magqdata = [] #magnetisation squared
+		T += dT
+		beta = 1./T
+		exptable = [1.,-1.,np.exp(-4.*J*beta),-1.,np.exp(-8.*J*beta),1.,-1.,1.,-1.] #lookuptable for the exponents we may need
+
+		#equilibrate
+		for i in range(int(equil)):
+			total_magnetisation += MCMC_step(sys,H,L,exptable)
+		#measure
+		for i in range(int(meas)):
+			total_magnetisation += MCMC_step(sys,H,L,exptable)
+			m = np.abs(total_magnetisation)/float(L*H)
+			magdata.append(m)
+			magqdata.append((m)**2)
+
+		#extract means and their standard deviations
+		magnetisations.append(np.mean(magdata))
+		magnetisationsq.append(np.mean(magqdata))
+		magnetisations_std.append(np.std(magdata))
+		magnetisationsq_std.append(np.std(magqdata))
+
+		#save temperature
+		temperatures.append(T)
+
+		#do a binning analysis of the data and store it
+		statistics_magnetisation.append(binning(magdata))
+		statistics_magnetisationq.append(binning(magqdata))
+
+	#calculate correlation times for all temperatures from the binning results. binning level 14 seems okay, but still not good at certain temperatures (CSD)
+	taus = []
+	tausq = []
+	for i in range(len(temperatures)):
+		taus.append(0.5*(statistics_magnetisation[i][1][14]**2/magnetisations_std[i]**2*(meas-1.)-1.))
+		tausq.append(0.5*(statistics_magnetisationq[i][1][14]**2/magnetisationsq_std[i]**2*(meas-1.)-1.))
+
+
+
+	#######################
+	###	Plotting
+	#######################
+
+	#correlation time magnetisation
+	plt.plot(temperatures,taus)
+	plt.title("Correlation time Metropolis, magnetization")
+	plt.xlabel("Temperature")
+	plt.ylabel("Correlation time")
+	#plt.show()
+	plt.savefig("metropolis_corr_times_m.pdf")
+	plt.clf()
+
+	#correlation time magnetisation squared
+	plt.plot(temperatures,taus)
+	plt.title("Correlation time Metropolis, magnetization$^2$")
+	plt.xlabel("Temperature")
+	plt.ylabel("Correlation time")
+	#plt.show()
+	plt.savefig("metropolis_corr_times_mq.pdf")
+	plt.clf()
+
+	#binning analysis results
+	for i in range(len(statistics_magnetisation)):
+		plt.plot(statistics_magnetisation[i][0],statistics_magnetisation[i][1])
+	leg = []
+	for i in range(len(temperatures)):
+		leg.append("T=" + str(round(temperatures[i],3)))
+	plt.legend(leg)
+	plt.title("Binning analysis Metropolis, magnetisation")
+	plt.xlabel("Binning Level")
+	plt.ylabel("Error")
+	#plt.show()
+	plt.savefig("metropolis_binning_m.pdf")
+	plt.clf()
+
+	for i in range(len(statistics_magnetisationq)):
+		plt.plot(statistics_magnetisationq[i][0],statistics_magnetisationq[i][1])
+	leg = []
+	for i in range(len(temperatures)):
+		leg.append("T=" + str(round(temperatures[i],3)))
+	plt.legend(leg)
+	plt.title("Binning analysis Metropolis, magnetisation$^2$")
+	plt.xlabel("Binning Level")
+	plt.ylabel("Error")
+	#plt.show()
+	plt.savefig("metropolis_binning_mq.pdf")
+	plt.clf()
+
+	#calculate errorbars
+	#binning level 14 seems reasonable for most temperatures
+	errorbars_m = []
+	errorbars_mq = []
+	for i in range(len(temperatures)):
+		errorbars_m.append(statistics_magnetisation[i][1][14])
+		errorbars_mq.append(statistics_magnetisationq[i][1][14])
+
+	#plot magnetisation with errorbars
+	plt.errorbar(temperatures,magnetisations,yerr=errorbars_m,fmt='o',ms=1)
+	plt.title("Magnetisation, Metropolis")
+	plt.xlabel("Temperature")
+	plt.ylabel("Magnetisation")
+	#plt.show()
+	plt.savefig("metropolis_m.pdf")
+	plt.clf()
+
+	#plot magnetisation squared with errorbars
+	plt.errorbar(temperatures,magnetisationsq,yerr=errorbars_mq,fmt='o',ms=1)
+	plt.title("Magnetisation$^2$, Metropolis")
+	plt.xlabel("Temperature")
+	plt.ylabel("Magnetisation^2")
+	#plt.show()
+	plt.savefig("metropolis_mq.pdf")
+	plt.clf()
+
+
+
+
+
+
+
+
+
+#run
+run_metropolis_simulation(0.1)
+run_cluster_simulation(0.1)