<jupyter_start><jupyter_text>cf. [`ComputationalPhysics/doc/Programs/LecturePrograms/programs/StatPhys/python/ising2dim.py`](https://github.com/CompPhysics/ComputationalPhysics/blob/master/doc/Programs/LecturePrograms/programs/StatPhys/python/ising2dim.py)<jupyter_code>import numpy
import numpy as np
import sys
import math
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Periodic boundary conditions<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Set up spin matrix, initialize to ground state<jupyter_code>size = 256 # L_x
temp = 10. # temperature T
spin_matrix = np.zeros( (size,size), np.int8) + 1
spin_matrix<jupyter_output><empty_output><jupyter_text>Create and initialize variables<jupyter_code>E = M = 0 
E_av = E2_av = M_av = M2_av = Mabs_av = 0<jupyter_output><empty_output><jupyter_text>Setup array for possible energy changes<jupyter_code>w = np.zeros(17, np.float64)
for de in xrange(-8,9,4):
    print de 
    w[de+8] = math.exp(-de/temp)
print w<jupyter_output>[ 2.22554093  0.          0.          0.          1.4918247   0.          0.
  0.          1.          0.          0.          0.          0.67032005
  0.          0.          0.          0.44932896]<jupyter_text>Calculate initial magnetization<jupyter_code>M = spin_matrix.sum()
print M<jupyter_output>65536<jupyter_text>Calculate initial energy<jupyter_code># 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)))<jupyter_output><empty_output><jupyter_text>Metropolis MonteCarlo computation, 1 single step or iteration, done explicitly:<jupyter_code>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])<jupyter_output>True<jupyter_text>Accept (if True)!<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Initialize (all spins up), explicitly shown<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>Setup array for possible energy changes<jupyter_code>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)<jupyter_output>[  2.98095801e+03   0.00000000e+00   0.00000000e+00   0.00000000e+00
   5.45981483e+01   0.00000000e+00   0.00000000e+00   0.00000000e+00
   1.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00
   1.83156393e-02   0.00000000e+00   0.00000000e+00   0.00000000e+00
   3.35462624e-04]<jupyter_text>Importing from the script `ising2dim.py`<jupyter_code>import os
print(os.getcwd())
print(os.listdir( os.getcwd() ))
sys.path.append('./')<jupyter_output><empty_output><jupyter_text>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`.<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>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,<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>From drafts<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>From [CLaigit](https://github.com/CLaigit/isingmodel/blob/master/ising2d.cu)<jupyter_code>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()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14688383887642642, 'fraction_numerical': 0.04915634448454111, 'mean_word_length': 5.244823386114494, 'pattern_counts': {'":': 0, '<': 62, '<?xml version=': 0, '>': 61, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 15, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '42040', 'license': 'Apache-2.0', 'n_tokens_mistral': 4545, 'n_tokens_neox': 4250, 'n_words': 725, 'path': 'Cpp/Ising/ising.ipynb', 'repo_name': 'ernestyalumni/CompPhys', 'size': '10253'}