Python programs
From Statistical Mechanics: Algorithms and Computations
Python programs Statistical Mechanics: Algorithms and Computations discusses more than 100 computer programs. In the book, these programs are presented in pseudocode. In this page, I collect some of these programs in Python, a modern computer language that, as the reader can see below, resembles pseudocode in Statistical Mechanics: Algorithms and Computations. In the classroom, I now mostly use Python.
[edit] Python programs in Chapter 1: Monte Carlo methods
##========+=========+=========+=========+=========+=========+=========+= ## PROGRAM : direct_pi.py ## PURPOSE : implement SMAC algorithm 1.1 (direct-pi) (cf Program:Direct_pi ) ## COMMENT : note that (x,y) is a tuple ## LANGUAGE: Python 2.5 ##========+=========+=========+=========+=========+=========+=========+ from random import uniform def direct_pi(N): N_hits=0 for i in range(N): x,y=uniform(-1.,1.),uniform(-1.,1.) if x**2 + y**2 < 1: N_hits += 1 return N_hits # for k in range(10): print k, direct_pi(10000)
##========+=========+=========+=========+=========+=========+=========+= ## PROGRAM : markov_pi.py ## PURPOSE : implement SMAC algorithm 1.2 (markov-pi) (cf Program:Markov_pi) ## COMMENT : note that (x,y) and (del_x, del_y) are tuples ## LANGUAGE: Python 2.5 ##========+=========+=========+=========+=========+=========+=========+ from random import uniform def markov_pi(delta,N): x,y=1.,1. # (Tuple assignment) N_hits=0 for i in range(N): del_x,del_y=uniform(-delta,delta),uniform(-delta,delta) if abs(x+del_x) < 1 and abs(y+del_y) < 1: x,y=x+del_x, y+del_y if x**2 + y**2 < 1: N_hits += 1 return N_hits ## for k in range(10): print k, markov_pi(0.3,10000)
###========+=========+=========+=========+=========+=========+=========+=
## PROGRAM : ran_perm.py
## PURPOSE : illustrate the built-in function "shuffle",
## implement SMAC algorithm 1.11 (ran-perm)
## in two different ways.
## LANGUAGE: Python 2.5
##========+=========+=========+=========+=========+=========+=========+
import random
#
# python script for ran_perm
#
def ran_perm_1(L):
random.shuffle(L)
return L
#
# explicit program using the bucket-and-shelf algorithm
#
def ran_perm_2(L):
k=len(L)-1
while k > 0:
j=random.randint(0,k)
L[k],L[j] = L[j],L[k]
k -= 1
return L
#
# bucket-and-shelf algorithm with remove
#
def ran_perm_3(L):
M=[]
while L != []:
k=random.choice(L)
L.remove(k)
M.append(k)
return M
#
L= [k for k in range(10)] # List comprehension, see LP p. 365
print ran_perm_1(L)
L= [k for k in range(10)]
print ran_perm_2(L)
L= [k for k in range(10)]
print ran_perm_3(L)
###========+=========+=========+=========+=========+=========+=========+=
## PROGRAM : naive_gauss.py
## PURPOSE : implement SMAC algorithm 1.17 (naive-gauss),
## with plotting by matplotlib (pylab)
## LANGUAGE: Python 2.5
##========+=========+=========+=========+=========+=========+=========+
from random import uniform, random
import math, pylab, numpy
def naive_gauss(K):
sigma = math.sqrt(K/12.)
sum=0
for k in range(K):
sum += uniform(-0.5,0.5)
return sum/sigma
#
data=[]
for k in range(10000):
x=naive_gauss(2)
data.append(x)
#
pylab.hist(data,50,normed=1) # see "http://matplotlib.sourceforge.net/" for more options of "hist"
x=numpy.linspace(-3,3,100)# equally spaced vector between -3 and 3
pylab.plot(x,numpy.exp(-x**2/2)/math.sqrt(2*numpy.pi))
pylab.show()
###========+=========+=========+=========+=========+=========+=========+= ## PROGRAM : gauss.py ## PURPOSE : implement SMAC algorithm 1.18 (gauss) ## COMMENT : note that the module random has a function gauss! ## LANGUAGE: Python 2.5 ##========+=========+=========+=========+=========+=========+=========+ from random import uniform as ran import math def gauss(sigma): phi = ran(0,2*pi) U= -math.log(ran(0.,1.)) r= sigma* math.sqrt(2*U) x= r* math.cos(phi) y= r* math.sin(phi) return x,y # for k in range(10): x,y=gauss(1.) print x,y
###========+=========+=========+=========+=========+=========+=========+=
## PROGRAM : gauss_plot.py
## PURPOSE : implement SMAC algorithm 1.18 (gauss)
## with plotting by matplotlib (pylab)
## COMMENT : note that the module random has a function gauss!
## LANGUAGE: Python 2.5
##========+=========+=========+=========+=========+=========+=========+
from random import uniform, random
import math
from pylab import *
def gauss(sigma):
phi = uniform(0,2*pi)
U= -math.log(random())
r= sigma* math.sqrt(2*U)
x= r* math.cos(phi)
y= r* math.sin(phi)
return x,y
#
figure(1)
axis('scaled')
axis([-4,4,-2,2])
for k in range(1000):
x,y=gauss(1.)
plot([x],[y],'ro-')
show()
###========+=========+=========+=========+=========+=========+=========+= ## PROGRAM : direct_surface_plot.py ## PURPOSE : implement SMAC algorithm 1.22 (direct_surface) (cf Program:Direct_surface) ## with plotting by matplotlib (pylab) ## LANGUAGE: Python 2.5 ##========+=========+=========+=========+=========+=========+=========+ from random import uniform, gauss, random import math from pylab import * def direct_sphere(d): x_vec = [gauss(0.,1.) for k in range(d)] sum=0 for x in x_vec: sum += x**2 for k in range(d): x_vec[k]=x_vec[k]/math.sqrt(sum) return x_vec # figure(1) axis('equal') x_vec=[] y_vec=[] for k in range(100): x,y=direct_sphere(2) plot([x],[y],'ro-') savefig('direct_surface.ps') show()
###========+=========+=========+=========+=========+=========+=========+=
## PROGRAM : data_bunch.py
## PURPOSE : implement SMAC algorithm 1.28 (data_bunch)
## then perform a test.
## LANGUAGE: Python 2.5
##========+=========+=========+=========+=========+=========+=========+
import math
def data_bunch(data):
N=0
New_data=[]
mean=0
meansq=0
while (len(data) >=2):
N +=2
x=data.pop() # note that we are reversing the order of data
y=data.pop()
New_data.append((x+y)/2.0)
mean += x+y
meansq += x**2+y**2
error=math.sqrt((meansq/float(N) - (mean/N)**2)/float(N))
mean=mean/N
return(mean, error, New_data)
##
## test (for independent data)
##
from random import uniform
power=24
data=[uniform(0,1) for k in range(2**power)]
for iter in range(power):
N_data=len(data)
mean,error,data=data_bunch(data)
print iter, N_data,mean, error
[edit] Python programs in Chapter 2: Hard disks and spheres
###========+=========+=========+=========+=========+=========+=========+= ## PROGRAM : direct_disks.py ## PURPOSE : implement SMAC algorithm 2.7 (direct-disks) for 4 disks in a ## square (with periodic boundary conditions) (cf Program:Direct_disks) ## LANGUAGE: Python 2.5 ##========+=========+=========+=========+=========+=========+=========+ from random import random import math def dist(x,y): d_x= abs(x[0]-y[0])%1 d_x = min(d_x,1-d_x) d_y= abs(x[1]-y[1])%1 d_y = min(d_y,1-d_y) return d_x**2 + d_y**2 N=16 sigma=0.1 sigma_sq=sigma**2 condition = False while condition == False: L=[(random(),random())] for k in range(1,N): b=(random(),random()) min_dist = min(dist(b,x) for x in L) if min_dist < 4*sigma_sq: condition = False break else: L.append(b) condition = True print L
###========+=========+=========+=========+=========+=========+=========+= ## PROGRAM : direct_disks_any.py ## PURPOSE : implement SMAC algorithm 2.8 (direct-disks_any) ## plot a curve, just as in figure 2.21. ## LANGUAGE: Python 2.5 ##========+=========+=========+=========+=========+=========+=========+ from random import random import math, pylab def dist(x,y): # box of size 1x1 d_x= abs(x[0]-y[0])%1 d_x = min(d_x,1-d_x) d_y= abs(x[1]-y[1])%1 d_y = min(d_y,1-d_y) return d_x**2 + d_y**2 N=16 N_samp=1000000 eta_list=[] for iter in range(N_samp): L=[(random(),random()) for k in range(N)] sigma = min(dist(L[k],L[l]) for k in range(N) for l in range(k+1,N))/4. eta=math.pi*sigma*N eta_list.append(eta) eta_list.sort() # the integrated histogram is related to the sorted vector y_vec=[1-k/float(N_samp) for k in range(N_samp)] pylab.semilogy(eta_list,y_vec) x_vec=[k/float(100) for k in range(100)] y_vec=[math.exp(-2.*(N-1)*x) for x in x_vec] pylab.semilogy(x_vec,y_vec) pylab.axis([0,0.4,0.000001,1]) pylab.show()
###========+=========+=========+=========+=========+=========+=========+= ## PROGRAM : markov_disks_box.py ## PURPOSE : implement SMAC algorithm 2.9 (markov-disks) for 4 disks in a square (without ## periodic boundary conditions) (cf Program:Markov_disks) ## LANGUAGE: Python 2.5 ##========+=========+=========+=========+=========+=========+=========+ from random import uniform as ran, choice L=[(0.25,0.25),(0.75,0.25),(0.25,0.75),(0.75,0.75)] sigma=0.20 sigma_sq=sigma**2 delta=0.15 for iter in range(1000000): a = choice(L) L.remove(a) b = (a[0] + ran(-delta,delta),a[1] + ran(-delta,delta)) min_dist = min((b[0]-x[0])**2 + (b[1]-x[1])**2 for x in L) box_cond = min(b[0],b[1]) < sigma or max(b[0],b[1]) >1-sigma if box_cond or min_dist < 4*sigma**2: L.append(a) else: L.append(b) print L
###========+=========+=========+=========+=========+=========+=========+=
## PROGRAM : pocket_disks.py
## PURPOSE : implement SMAC algorithm 2.18 (hard-sphere-cluster)
## for four disks in a square of sides 1
## LANGUAGE: Python 2.5
##========+=========+=========+=========+=========+=========+=========+
from random import uniform as ran, choice
def box_it(x):
x[0]=x[0]%1.
if x[0] < 0 : x[0]=x[0]+1
x[1]=x[1]%1.
if x[1] < 0 : x[1]=x[1]+1
return x
def dist(x,y):
d_x= abs(x[0]-y[0])%1
d_x = min(d_x,1-d_x)
d_y= abs(x[1]-y[1])%1
d_y = min(d_y,1-d_y)
return d_x**2 + d_y**2
def T(x,Pivot):
x=(2*Pivot[0]-x[0],2*Pivot[1]-x[1])
return x
#
# Program starts here
#
Others=[(0.25,0.25),(0.25,0.75),(0.75,0.25),(0.75,0.75)]
sigma_sq=0.15**2
for iter in range(10000):
a = choice(Others)
Others.remove(a)
Pocket = [a]
Pivot=(ran(0,1),ran(0,1))
while Pocket != []:
a = choice(Pocket)
Pocket.remove(a)
a = T(a,Pivot)
for b in Others[:]: # "Others[:]" is a copy of "Others"
if dist(a,b) < 4*sigma_sq:
Others.remove(b)
Pocket.append(b)
Others.append(a)
print Others, ' ending config '
[edit] Python programs in Chapter 3: Density matrices and path integrals
###========+=========+=========+=========+=========+=========+=========+=
## PROGRAM : harmonic-wavefunction.py
## PURPOSE : implement SMAC algorithm 3.1 (harmonic=wavefunction)
## using a dictionary
## The 'value' of each 'key' x is the list
## [psi_0(x), psi_1(x), ...]
## LANGUAGE: Python 2.5
##========+=========+=========+=========+=========+=========+=========+
import numpy, pylab
from math import pi, exp, sqrt
grid=numpy.linspace(-5,5,100)
psi={}
pylab.figure(1,figsize=(6,10))
for x in grid:
psi[x]=[pi**(-0.25)*exp(-x**2/2)]
psi[x].append(sqrt(2.)*x*psi[x][0])
for n in range(2,50):
psi[x].append(sqrt(2./n)*x*psi[x][n -1] - sqrt((n-1.)/n)*psi[x][n-2])
for k in range(40):
y=[psi[x][k] + k for x in grid]
pylab.plot(grid,y,'r-')
pylab.xlabel('x')
pylab.ylabel('psi_n(x) (offset)')
pylab.xlim(-5,5)
pylab.savefig('harmonic_wavefunction.eps')
pylab.show()
###========+=========+=========+=========+=========+=========+=========+= ## PROGRAM : naive-harmonic-path.py ## PURPOSE : implement SMAC algorithm 3.4 (naive-harmonic-path) ## pylab histogram output as in SMAC figure 3.10 (p. 151) ## LANGUAGE: Python 2.5 ##========+=========+=========+=========+=========+=========+=========+ import numpy, pylab from math import pi, exp, sqrt from random import random, uniform, randint N=20 beta=4. delta=0.1 del_tau=beta/N def rho_free(x,xp,beta): rho_free=exp(-(x-xp)**2/(2*beta)) return rho_free x =[0. for k in range(N)] data=[] y =[k/float(N)*beta for k in range(N)] for iter in range(1000000): k =randint(0,N-1) x_old=x[k] x_new=x_old + uniform(-delta,delta) xp=x[(k+1)%N] xm=x[k-1] pi_old = rho_free(xm,x_old,del_tau)*rho_free(x_old,xp,del_tau) * exp(-del_tau*x_old**2/2.) pi_new = rho_free(xm,x_new,del_tau)*rho_free(x_new,xp,del_tau) * exp(-del_tau*x_new**2/2.) Upsilon=pi_new/pi_old if (random() < Upsilon): x[k]=x_new if iter%10 ==0: data.append(x[0]) pylab.hist(data,bins=20,range=[-2,2],normed=True) pylab.show()
##========+=========+=========+=========+=========+=========+=========+= ## PROGRAM : levy.py ## PURPOSE : implement SMAC algorithm 3.5 (levy-free-path) ## LANGUAGE: Python 2.5 ##========+=========+=========+=========+=========+=========+=========+ import numpy, pylab,time from math import pi, exp, sqrt from random import gauss N=10000 beta=1. Del_tau=beta/N x=[0 for k in range(N+1)] y=range(N+1) for k in range(1,N): Del_p=(N-k)*Del_tau x_mean=(Del_p*x[k-1] + Del_tau*x[N])/(Del_tau + Del_p) sigma=sqrt(1./(1./Del_tau + 1./Del_p)) x[k]=gauss(x_mean,sigma) pylab.figure(1,figsize=(8.,12.)) pylab.plot(x,y,'b-') pylab.show()
##========+=========+=========+=========+=========+=========+=========+= ## PROGRAM : levy_harmonic_path.py ## PURPOSE : implement SMAC algorithm 3.6 (levy-harmonic-path) ## (with plotting by pylab) ## LANGUAGE: Python 2.5 ##========+=========+=========+=========+=========+=========+=========+ import numpy, pylab from math import pi, exp, sqrt,sinh,tanh from random import gauss N=10000 beta=100. Del_tau=beta/N x=[0 for k in range(N+1)] y=range(N+1) for k in range(1,N): Upsilon_1 = 1./tanh(Del_tau)+1./tanh((N-k)*Del_tau) Upsilon_2 = x[k-1]/sinh(Del_tau)+x[N]/sinh((N-k)*Del_tau) x_mean=Upsilon_2/Upsilon_1 sigma=1./sqrt(Upsilon_1) x[k]=gauss(x_mean,sigma) pylab.figure(1,figsize=(8.,12.)) pylab.plot(x,y,'b-') pylab.show()
[edit] Python programs in Chapter 4: Bosons
###========+=========+=========+=========+=========+=========+=========+=
## PROGRAM : naive-bosons.py
## PURPOSE : implement algorithm 4.3 (naive-bosons)
## with plot of condensate number (using pylab)
## LANGUAGE: Python 2.5
##========+=========+=========+=========+=========+=========+=========+
from pylab import *
Temp=list(linspace(0,4,17))
del Temp[0]
E=[0,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4]
N_0_vec=[]
for T in Temp:
N_0=0
Z=0
beta=1./T
for s1 in range(35):
for s2 in range(s1,35):
for s3 in range(s2,35):
for s4 in range(s3,35):
for s5 in range(s4,35):
dummy=0
E_tot=E[s1] + E[s2] + E[s3] + E[s4] + E[s5]
if s1 ==0: dummy+=1
if s2 ==0: dummy+=1
if s3 ==0: dummy+=1
if s4 ==0: dummy+=1
if s5 ==0: dummy+=1
N_0 += dummy*math.exp(-beta*E_tot)
Z += math.exp(-beta*E_tot)
print T, Z,' T, Z'
N_0_vec.append(N_0/Z/5.)
figure(1)
plot(Temp,N_0_vec,'r-')
show()
###========+=========+=========+=========+=========+=========+=========+=
## PROGRAM : canonic_recursion.py
## PURPOSE : implement SMAC algorithm 4.6 (canonic-recursion)
## with plot of cycle length distribution in pylab (SMAC fig. 4.18)
## LANGUAGE: Python 2.5
##========+=========+=========+=========+=========+=========+=========+
from pylab import *
import math
#--single-particle partition function in a three-dimensional harmonic trap
def z(beta,k):
sum = 1/(1- math.exp(-k*beta))**3
return sum
#--begin program
N=1000
T_star=0.5
T=T_star*N**0.3333
beta=1./T
Z=[1.]
for M in range(1,N+1):
Z.append(sum(Z[k] * z(beta,M-k) for k in range(M))/M)
pi_list=[z(beta,k)*Z[N-k]/Z[N]/N for k in range(1,N+1)]
N_list=linspace(1,N,N)
#--begin figure output
figure(1)
plot(N_list,pi_list)
ylim(0,.003)
xlabel('$k$')
ylabel('$\pi(k)$')
title('$T^*$='+str(T_star)+' $N$= '+str(N))
show()
##========+=========+=========+=========+=========+=========+=========+
## PROGRAM : direct_harmonic_bosons.py
## PURPOSE : implement SMAC algorithm 4.6 (direct-harmonic-bosons)
## with plot of two-dimensional snapshot in pylab (SMAC fig. 4.23)
## LANGUAGE: Python 2.5
##========+=========+=========+=========+=========+=========+=========+
from pylab import *
from math import sqrt, sinh, tanh,exp
from random import uniform as ran,gauss
def z(beta,k):
sum = 1/(1- exp(-k*beta))**3
return sum
def canonic_recursion(beta,N):
Z=[1.]
for M in range(1,N+1):
Z.append(sum(Z[k] * z(beta,M-k) for k in range(M))/M)
return Z
def pi_list_make(Z,M):
pi_list=[0]+[z(beta,k)*Z[M-k]/Z[M]/M for k in range(1,M+1)]
pi_sum=[0]
for k in range(1,M+1):
pi_sum.append(pi_sum[k-1]+pi_list[k])
return pi_sum
def tower_sample(data,Upsilon): #naive tower sampling, cf. SMAC Sect. 1.2.3
for k in range(len(data)):
if Upsilon<data[k]: break
return k
def levy_harmonic_path(Del_tau,N): #
beta=N*Del_tau
xN=gauss(0.,1./sqrt(2*tanh(beta/2.)))
x=[xN]
for k in range(1,N):
Upsilon_1 = 1./tanh(Del_tau)+1./tanh((N-k)*Del_tau)
Upsilon_2 = x[k-1]/sinh(Del_tau)+xN/sinh((N-k)*Del_tau)
x_mean=Upsilon_2/Upsilon_1
sigma=1./sqrt(Upsilon_1)
x.append(gauss(x_mean,sigma))
return x
#-------------------------------------------------------
# -- main program starts here, don't choose N too large
#-------------------------------------------------------
N=512
T_star=.3
T=T_star*N**(1./3.)
beta=1./T
Z=canonic_recursion(beta,N)
M=N
#Dic_cycles={}
x_config=[]
y_config=[]
while M > 0:
pi_sum=pi_list_make(Z,M)
Upsilon=ran(0,1)
k=tower_sample(pi_sum,Upsilon)
x_config+=levy_harmonic_path(beta,k)
y_config+=levy_harmonic_path(beta,k)
M -= k
# if Dic_cycles.has_key(k): Dic_cycles[k]+=1
# else: Dic_cycles[k]=1
#ll = Dic_cycles.keys()
#ll.sort()
# print k, Dic_cycles[k]
##--begin figure output
figure(1)
axis('scaled')
xlim(-5.,5.)
ylim(-5.,5.)
plot(x_config,y_config,'ro')
xlabel('$x$')
ylabel('$y$')
title('3d bosons in trap (2d projection):'+' $N$= '+str(N)+' T* =' + str(T_star))
show()
[edit] Python programs in Chapter 5: Order and disorder in spin systems
##---------------------------------------------------------------------
## PROGRAM : markov_ising.py
## PURPOSE : implement SMAC algorithm 5.7 (markov-ising).
## (Metropolis algorithm). With plotting of configurations
## (using pylab)
## LANGUAGE: Python 2.5
##---------------------------------------------------------------------
from random import uniform as ran, randint,choice
from math import exp
from pylab import hist, figure,savefig,show,plot,axis
##---------------------------------------------------------------------
## Sample geometry
##---------------------------------------------------------------------
def square_neighbors(L):
N = L*L
site_dic = {}
x_y_dic = {}
for j in range(N):
row = j//L
column = j-row*L
site_dic[(row,column)] = j
x_y_dic[j] = (row,column)
nbr=[]
for j in range(N):
row,column = x_y_dic[j]
right_nbr = site_dic[row,(column+1)%L]
up_nbr = site_dic[(row+1)%L,column]
left_nbr = site_dic[row,(column-1+L)%L]
down_nbr = site_dic[(row-1+L)%L,column]
nbr.append((right_nbr,up_nbr,left_nbr,down_nbr))
nbr = tuple(nbr)
return nbr,site_dic,x_y_dic
#
# Program starts here
#
L=32
N=L*L
S=[choice([-1,1]) for k in range(N)]
beta=0.42
nbr,site_dic,x_y_dic=square_neighbors(L)
for i_sweep in range(100):
for iter in range(N):
k=randint(0,N-1)
h=sum(S[nbr[k][j]] for j in range(4))
Delta_E=2.*h*S[k]
Upsilon=exp(-beta*Delta_E)
if ran(0.,1.) < Upsilon: S[k] = -S[k]
figure()
x_plus=[]
y_plus=[]
x_minus=[]
y_minus=[]
for i in range(N):
x,y=x_y_dic[i]
if S[i] ==1:
x_plus.append(x)
y_plus.append(y)
else:
x_minus.append(x)
y_minus.append(y)
axis('scaled')
axis([-0.5,L-0.5,-0.5,L-.5])
plot(x_plus,y_plus,'r+',markersize=10)
plot(x_minus,y_minus,'bx',markersize =10)
savefig('test2.eps')
show()
##---------------------------------------------------------------------
## PROGRAM : cluster_ising.py
## PURPOSE : implement SMAC algorithm 5.9 (cluster-ising).
## complete program with neighbor definitions and graphics
## LANGUAGE: Python 2.5
##---------------------------------------------------------------------
from random import uniform as ran,randint,choice
from math import exp
from pylab import hist, figure,savefig,show,plot,axis
##---------------------------------------------------------------------
## Sample geometry
##---------------------------------------------------------------------
def square_neighbors(L):
N = L*L
site_dic = {}
x_y_dic = {}
for j in range(N):
row = j//L
column = j-row*L
site_dic[(row,column)] = j
x_y_dic[j] = (row,column)
nbr=[]
for j in range(N):
row,column = x_y_dic[j]
right_nbr = site_dic[row,(column+1)%L]
up_nbr = site_dic[(row+1)%L,column]
left_nbr = site_dic[row,(column-1+L)%L]
down_nbr = site_dic[(row-1+L)%L,column]
nbr.append((right_nbr,up_nbr,left_nbr,down_nbr))
nbr = tuple(nbr)
return nbr,site_dic,x_y_dic
#
# Program starts here
#
L=32
N=L*L
S=[choice([-1,1]) for k in range(N)]
beta=0.4407
p=1 - exp(-2*beta)
nbr,site_dic,x_y_dic=square_neighbors(L)
for iter in range(100):
Pocket = [k]
Cluster = [k]
N_cluster = 1
while Pocket != []:
k =choice(Pocket)
for l in nbr[k]:
if S[l] == S[k] and l not in Cluster and ran(0,1) < p:
N_cluster += 1
Pocket.append(l)
Cluster.append(l)
Pocket.remove(k)
for k in Cluster: S[k] = - S[k]
print iter, N_cluster
figure(1)
x_plus=[]
y_plus=[]
x_minus=[]
y_minus=[]
for i in range(N):
x,y=x_y_dic[i]
if S[i] ==1:
x_plus.append(x)
y_plus.append(y)
else:
x_minus.append(x)
y_minus.append(y)
axis('scaled')
axis([-0.5,L-0.5,-0.5,L-.5])
plot(x_plus,y_plus,'r+',markersize=10)
plot(x_minus,y_minus,'bx',markersize =10)
savefig('test2.eps')
show(1)
[edit] Python programs in Chapter 6: Entropic forces
###========+=========+=========+=========+=========+=========+=========+=
## PROGRAM : naive_pin.py
## PURPOSE : implement SMAC algorithm (6.1): naive-pin (a bit modified)
## with plotting by matplotlib (pylab)
## LANGUAGE: Python 2.5
##========+=========+=========+=========+=========+=========+=========+
from random import uniform as ran
import pylab
L=2
N=10
sigma=0.05
data=[]
Nsamp=0
while Nsamp < 1000:
Y = [ran(sigma,L-sigma) for k in range(N)]
Y.sort()
dmin=min(Y[k]-Y[k-1] for k in range(1,N))
if dmin > 2*sigma:
Nsamp +=1
data += Y
pylab.hist(data,range=[0,2],bins=200,normed=1) # see "http://matplotlib.sourceforge.net/" for more options of "hist"
pylab.axis([0.,2.,0.,2.])
pylab.savefig('clothespin.jpg')
pylab.show()
###========+=========+=========+=========+=========+=========+=========+= ## PROGRAM : direct_pin.py ## PURPOSE : implement SMAC algorithm 6.2 (direct-pin) ## with histogram plotting by matplotlib (pylab) ## LANGUAGE: Python 2.5 ##========+=========+=========+=========+=========+=========+=========+ from random import uniform as ran import pylab L=2 N=15 sigma=0.05 data=[] for iter in range(100000): Y = [ran(0,L-2*N*sigma) for k in range(N)] Y.sort() X = [Y[k]+(2*k+1)*sigma for k in range(N)] data += X pylab.hist(data,range=[0,2],bins=200,normed=1) # see "http://matplotlib.sourceforge.net/" for more options of "hist" pylab.axis([0.,2.,0.,2.]) pylab.savefig('clothespin.jpg') pylab.show()
###========+=========+=========+=========+=========+=========+=========+=
## PROGRAM : analytic_pin.py
## PURPOSE : implement SMAC equation (6.7): analytic cloth-pin density
## with plotting by matplotlib (pylab) (yields Fig 6.3)
## LANGUAGE: Python 2.5
##========+=========+=========+=========+=========+=========+=========+
from random import uniform as ran
import pylab
def Z_fun(k,L,sigma):
if k*2*sigma < L:
res=(L-k*2*sigma)**k
else:
res=0
return res
def binom_calc(N): # use a dictionary for the binomial coefficients
binom={}
binom[(0,0)]=1
for k in range(1,N+1):
binom[(k,0)]=1
binom[(k,k)]=1
for l in range (1,k):
binom[(k,l)]=binom[(k-1,l-1)] + binom[(k-1,l)]
return binom
L=2
sigma=0.05
N=15
x_vec=pylab.linspace(0,L,1000)
y_vec=[]
binom=binom_calc(N)
for x in x_vec:
proba=0
for k in range(N):
proba += binom[(N-1,k)]*Z_fun(k,x-sigma,sigma)*Z_fun(N-1-k,L-x-sigma,sigma)
proba=proba/Z_fun(N,L,sigma)
y_vec.append(proba)
pylab.plot(x_vec,y_vec)
pylab.show()
