<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 cf. ComputationalPhysics/doc/Programs/LecturePrograms/programs/StatPhys/python/ising2dim.py
Step2:  Periodic boundary conditions
Step3:  Set up spin matrix, initialize to ground state
Step4:  Create and initialize variables
Step5:  Setup array for possible energy changes
Step6:  Calculate initial magnetization
Step7:  Calculate initial energy
Step8:  Metropolis MonteCarlo computation, 1 single step or iteration, done explicitly
Step9:  Accept (if True)!
Step10:  Initialize (all spins up), explicitly shown
Step11:  Setup array for possible energy changes
Step12:  Importing from the script ising2dim.py
Step13:  Reading out data from ./IsingGPU/FileIO/output.h
Data is generated by the parallel Metropolis algorithm in CUDA C++ in the subdirectory ./IsingGPU/data/, which is done by the function process_avgs in ./IsingGPU/FileIO/output.h.  The values are saved as a character array, which then can be read in as a NumPy array of float32's.  Be sure to enforce, declare the dtype to be float32.
Step14:  For  
2^10 x 2^10 or 1024 x 1024 grid; 50000 trials, temperature T = 1.0, 1.005,...3. (temperature step of 0.005), so 400 different temperatures, 32 x 32 thread block,
Step15:  From drafts
Step16:  From CLaigit
<ASSISTANT_TASK:>
Python Code:
import numpy
import numpy as np
import sys
import math
import matplotlib.pyplot as plt
Explanation: cf. ComputationalPhysics/doc/Programs/LecturePrograms/programs/StatPhys/python/ising2dim.py
End of explanation
def periodic(i,limit,add): 
    
    Choose correct matrix index with periodic boundary conditions
    
    Input:
    - i:     Base index
    - limit: Highest \"legal\" index
    - add:   Number to add or subtract from i
    
    return (i + limit + add) % limit
Explanation: Periodic boundary conditions
End of explanation
size = 256 # L_x
temp = 10. # temperature T
spin_matrix = np.zeros( (size,size), np.int8) + 1
spin_matrix
Explanation: Set up spin matrix, initialize to ground state
End of explanation
E = M = 0 
E_av = E2_av = M_av = M2_av = Mabs_av = 0 
Explanation: Create and initialize variables
End of explanation
w = np.zeros(17, np.float64)
for de in xrange(-8,9,4):
    print de 
    w[de+8] = math.exp(-de/temp)
print w
Explanation: Setup array for possible energy changes
End of explanation
M = spin_matrix.sum()
print M
Explanation: Calculate initial magnetization
End of explanation
# for i in xrange(16): print i  r
# range creates a list, so if you do range(1, 10000000) it creates a list in memory with 9999999 elements.
# xrange is a sequence object that evaluates lazily. 
for j in xrange(size):
    for i in xrange(size):
        E -= spin_matrix.item(i,j) * (spin_matrix.item(periodic(i,size,-1),j) + spin_matrix.item(i,periodic(j,size,1)))
Explanation: Calculate initial energy
End of explanation
x = int(np.random.random()*size)
print(x)
y = int(np.random.random()*size)
print(y)
deltaE = 2*spin_matrix.item(i,j) * \
          (spin_matrix.item(periodic(x,size,-1),y) + spin_matrix.item(periodic(x,size,1),y) + \
          spin_matrix.item(x,periodic(y,size,-1))+spin_matrix.item(x,periodic(y,size,1)))  
print(deltaE)
print( w[deltaE + 8] )
np.random.random()
print( np.random.random() <= w[deltaE+8])
Explanation: Metropolis MonteCarlo computation, 1 single step or iteration, done explicitly:
End of explanation
print( spin_matrix[x,y] )
print( spin_matrix.item(x,y) )
spin_matrix[x,y] *= -1
M += 2*spin_matrix[x,y]
E += deltaE
print(spin_matrix.item(x,y))
print(M)
print(E)
import pygame
Explanation: Accept (if True)!
End of explanation
Lx=256; Ly=256
spin_matrix = np.zeros((Lx,Ly),np.int8)
print(spin_matrix.shape)
spin_matrix.fill(1)
spin_matrix
def initialize_allup( spin_matrix, J=1.0 ):
    Lx,Ly = spin_matrix.shape
    
    spin_matrix.fill(1)
    M = spin_matrix.sum()
    # Calculate initial energy  
    E=0
    for j in xrange(Ly):
        for i in xrange(Lx):
            E += (-J)*spin_matrix.item(i,j) * \
                (spin_matrix.item(periodic(i,Lx,+1),j) + spin_matrix.item(i,periodic(j,Ly,1)) )
    print "M: ",M," E: ", E
    return E,M    
E,M = initialize_allup( spin_matrix)
def initialize_allup1( spin_matrix, J=1.0 ):
    Lx,Ly = spin_matrix.shape
    
    spin_matrix.fill(1)
    M = spin_matrix.sum()
    # Calculate initial energy  
    E=0
    for j in xrange(Ly):
        for i in xrange(Lx):
            E -= J*spin_matrix.item(i,j) * \
                (spin_matrix.item(periodic(i,Lx,-1),j) + spin_matrix.item(i,periodic(j,Ly,1)) )
    print "M: ",M," E: ", E
    return E,M   
E,M = initialize_allup( spin_matrix)
Lx=512; Ly=512
spin_matrix = np.zeros((Lx,Ly),np.int8)
E,M = initialize_allup1( spin_matrix)
E,M = initialize_allup( spin_matrix)
Lx=1024; Ly=1024
print(Lx*Ly)
spin_matrix = np.zeros((Lx,Ly),np.int8)
E,M = initialize_allup1( spin_matrix)
E,M = initialize_allup( spin_matrix)
math.pow(2,31)
Explanation: Initialize (all spins up), explicitly shown
End of explanation
temp = 1.0
w = np.zeros(17,np.float32)
for de in xrange(-8,9,4): # include +8
    w[de+8] = math.exp(-de/temp)  
print(w)
Explanation: Setup array for possible energy changes
End of explanation
import os
print(os.getcwd())
print(os.listdir( os.getcwd() ))
sys.path.append('./')
Explanation: Importing from the script ising2dim.py
End of explanation
avgsresults_GPU = np.fromfile("./IsingGPU/data/IsingMetroGPU.bin",dtype=np.float32)
print(avgsresults_GPU.shape)
print(avgsresults_GPU.size)
avgsresults_GPU = avgsresults_GPU.reshape(201,7) # 7 different averages
print(avgsresults_GPU.shape)
print(avgsresults_GPU.size)
avgsresults_GPU
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
T = avgsresults_GPU[:,0]
E_avg = avgsresults_GPU[:,1]
ax.scatter( T, E_avg)
plt.show()
Evar_avg = avgsresults_GPU[:,2]
plt.scatter( T, Evar_avg)
plt.show()
M_avg = avgsresults_GPU[:,3]
Mvar_avg = avgsresults_GPU[:,4]
absM_avg = avgsresults_GPU[:,5]
M4_avg = avgsresults_GPU[:,6]
#fig = plt.figure()
#ax = fig.add_subplot(4,1,1)
plt.scatter( T, M_avg)
#fig.add_subplot(4,1,2)
#plt.scatter(T,Mvar_avg)
#fig.add_subplot(4,1,3)
#plt.scatter(T,absM_avg)
#fig.add_subplot(4,1,4)
#plt.scatter(T,M4_avg)
plt.show()
plt.scatter(T,Mvar_avg)
plt.show()
plt.scatter(T,absM_avg)
plt.show()
plt.scatter(T,M4_avg)
plt.show()
Explanation: Reading out data from ./IsingGPU/FileIO/output.h
Data is generated by the parallel Metropolis algorithm in CUDA C++ in the subdirectory ./IsingGPU/data/, which is done by the function process_avgs in ./IsingGPU/FileIO/output.h.  The values are saved as a character array, which then can be read in as a NumPy array of float32's.  Be sure to enforce, declare the dtype to be float32.
End of explanation
avgsresults_GPU = np.fromfile("./IsingGPU/data/IsingMetroGPU.bin",dtype=np.float32)
print(avgsresults_GPU.shape)
print(avgsresults_GPU.size)
avgsresults_GPU = avgsresults_GPU.reshape( avgsresults_GPU.size/7 ,7) # 7 different averages
print(avgsresults_GPU.shape)
print(avgsresults_GPU.size)
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
T = avgsresults_GPU[:,0]
E_avg = avgsresults_GPU[:,1]
Evar_avg = avgsresults_GPU[:,2]
M_avg = avgsresults_GPU[:,3]
Mvar_avg = avgsresults_GPU[:,4]
absM_avg = avgsresults_GPU[:,5]
M4_avg = avgsresults_GPU[:,6]
ax.scatter( T, E_avg)
plt.show()
plt.scatter( T, Evar_avg)
plt.show()
plt.scatter( T, M_avg)
plt.show()
plt.scatter(T,Mvar_avg)
plt.show()
plt.scatter(T,absM_avg)
plt.show()
plt.scatter(T,M4_avg)
plt.show()
Explanation: For  
2^10 x 2^10 or 1024 x 1024 grid; 50000 trials, temperature T = 1.0, 1.005,...3. (temperature step of 0.005), so 400 different temperatures, 32 x 32 thread block,
End of explanation
avgsresults_GPU = np.fromfile("./IsingGPU/drafts/data/IsingMetroGPU.bin",dtype=np.float32)
print(avgsresults_GPU.shape)
print(avgsresults_GPU.size)
avgsresults_GPU = avgsresults_GPU.reshape( avgsresults_GPU.size/7 ,7) # 7 different averages
print(avgsresults_GPU.shape)
print(avgsresults_GPU.size)
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
T = avgsresults_GPU[:,0]
E_avg = avgsresults_GPU[:,1]
Evar_avg = avgsresults_GPU[:,2]
M_avg = avgsresults_GPU[:,3]
Mvar_avg = avgsresults_GPU[:,4]
absM_avg = avgsresults_GPU[:,5]
M4_avg = avgsresults_GPU[:,6]
ax.scatter( T, E_avg)
plt.show()
avgsresults_GPU = np.fromfile("./IsingGPU/drafts/IsingGPU/data/IsingMetroGPU_runs10.bin",dtype=np.float32)
print(avgsresults_GPU.shape)
print(avgsresults_GPU.size)
avgsresults_GPU = avgsresults_GPU.reshape( avgsresults_GPU.size/7 ,7) # 7 different averages
print(avgsresults_GPU.shape)
print(avgsresults_GPU.size)
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
T = avgsresults_GPU[:,0]
E_avg = avgsresults_GPU[:,1]
Evar_avg = avgsresults_GPU[:,2]
M_avg = avgsresults_GPU[:,3]
Mvar_avg = avgsresults_GPU[:,4]
absM_avg = avgsresults_GPU[:,5]
M4_avg = avgsresults_GPU[:,6]
ax.scatter( T, E_avg)
plt.show()
plt.scatter( T, Evar_avg)
plt.show()
plt.scatter( T, M_avg)
plt.show()
plt.scatter(T,Mvar_avg)
plt.show()
plt.scatter(T,absM_avg)
plt.show()
plt.scatter(T,M4_avg)
plt.show()
avgsresults_GPU = np.fromfile("./IsingGPU/drafts/IsingGPU/data/IsingMetroGPU.bin",dtype=np.float32)
print(avgsresults_GPU.shape)
print(avgsresults_GPU.size)
avgsresults_GPU = avgsresults_GPU.reshape( avgsresults_GPU.size/7 ,7) # 7 different averages
print(avgsresults_GPU.shape)
print(avgsresults_GPU.size)
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
T = avgsresults_GPU[:,0]
E_avg = avgsresults_GPU[:,1]
Evar_avg = avgsresults_GPU[:,2]
M_avg = avgsresults_GPU[:,3]
Mvar_avg = avgsresults_GPU[:,4]
absM_avg = avgsresults_GPU[:,5]
M4_avg = avgsresults_GPU[:,6]
ax.scatter( T, E_avg)
plt.show()
plt.scatter( T, Evar_avg)
plt.show()
plt.scatter( T, M_avg)
plt.show()
plt.scatter(T,Mvar_avg)
plt.show()
plt.scatter(T,absM_avg)
plt.show()
plt.scatter(T,M4_avg)
plt.show()
Explanation: From drafts
End of explanation
avgsresults_GPU = []
for temp in range(10,31,2): 
    avgsresults_GPU.append( np.fromfile("./data/ising2d_CLaigit" + str(temp) + ".bin",dtype=np.float64) )
avgsresults_GPU = np.array( avgsresults_GPU) 
print( avgsresults_GPU.shape, avgsresults_GPU.size)
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
T = avgsresults_GPU[:,0]
E_avg = avgsresults_GPU[:,1]
M_avg = avgsresults_GPU[:,2]
heat_cap_avg = avgsresults_GPU[:,3]
mag_sus_avg = avgsresults_GPU[:,4]
ax.scatter( T, E_avg)
plt.show()
plt.scatter( T, M_avg)
plt.show()
plt.scatter( T, heat_cap_avg)
plt.show()
plt.scatter( T, mag_sus_avg)
plt.show()
Explanation: From CLaigit
End of explanation
<END_TASK>