forked from yale-phys-382L/ising
-
Notifications
You must be signed in to change notification settings - Fork 0
/
ising.py
408 lines (346 loc) · 14.9 KB
/
ising.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
import sys
import os
import time
import csv
import logging
# import click
import numpy as np
# import logging
import matplotlib.pyplot as plt
import multiprocessing as mp
# from IsingLattice import IsingLattice as IsingLattice_c
from sys import exit, argv
# from pandas import DataFrame
# from tqdm import tqdm #fancy progress bar generator
# from ising_c import run_ising #import run_ising function from ising.py
# conditionally import python and C++ version of IsingLattice
has_cpp = False
from IsingLattice_python import IsingLattice as IsingLattice_py
if os.path.isfile('ising_lattice_lib.so'):
has_cpp = True
from IsingLattice_cpp import IsingLattice as IsingLattice_cpp
## main program
def make_B_generator(inp, t_final=None):
"""Return a generator that makes values of B (magetic field in each step)
note: you *should* play with this function
default implementation: always return 0
"""
for val in range(inp['n_steps']):
yield inp['B']
def make_T_generator(inp, t_final):
"""Return a generator that makes values of T (temperature in each step)
note: you *should* play with this function
default implementation:
start at T=t_top, and linearly decrease temperature to t_final
hold at t_final for n_burnin
hold at t_final for n_analyze
"""
n_slope = inp['n_steps'] - inp['n_burnin'] - inp['n_analyze']
if n_slope < 0:
print('fatal error: n_steps - n_burnin - n_slope < 0')
print('terminating program')
exit(2)
# get linearly decreasing values from t_top to t_final
for val in np.linspace(start=inp['t_top'], stop=t_final, num=n_slope):
yield val
for val in range(inp['n_burnin']):
yield t_final
for val in range(inp['n_analyze']):
yield t_final
def set_input(cmd_line_args):
"""Parse command-line parameters into an input dictionary.
use default values
input: sys.argv
use syntax of keyword:value on command line
return: dict[key] = value
note: Any value which can be turned into a number will be a float
if it has a '.', otherwise it will be a float.
"""
inp = dict()
inp['t_min'] = 1.6 # minimum temperature
inp['t_max'] = 3.6 # maximum temperature
inp['t_step'] = 0.1 # step size from min to max temperature
inp['t_top'] = 4.0 # start temperature (arbitrary; feel free to change)
inp['N'] = 10 # sqrt(lattice size) (i.e. lattice = N^2 points
inp['n_steps'] = 10000 # number of lattice steps in simulation
inp['n_burnin'] = 2000 # optional parameter, used as naive default
inp['n_analyze'] = 5000 # number of lattice steps at end of simulation calculated for averages and std.dev.
# inp['J'] = 1.0 # **great** default value -- spin-spin interaction strength
inp['B'] = 0.0 # magnetic field strength
inp['flip_perc'] = 0.1 # ratio of sites examined to flip in each step
inp['dir_out'] = 'data' # output directory for file output
inp['plots'] = False # whether or not plots are generated
inp['print_inp'] = False # temperature option
inp['use_cpp'] = True # use 1 for True and 0 for False
inp['date_output'] = False
inp['file_prefix'] = ''
inp['multiprocess'] = False
inp['skip_prog_print'] = False
for x in cmd_line_args[1:]:
if ':' in x:
try:
key, val = x.split(':')
try:
if '.' in val:
inp[key] = float(val)
print('%-20s'%('inp["%s"]'%key),'set to float ',inp[key])
elif val.lower() == 'false' or val.lower() == 'f':
inp[key] = False
elif val.lower() == 'true' or val.lower() == 't':
inp[key] = True
else:
inp[key] = int(val)
print('%-20s'%('inp["%s"]'%key),'set to int ',inp[key])
except:
inp[key] = val
print('%-20s'%('inp["%s"]'%key),'set to string ',inp[key])
except:
print('warning: input "%s" not added to arguments'%x)
else:
print('ignoring command line input: %s'%x)
if inp['print_inp']:
print('Printed list of input keys:')
for key in sorted(inp.keys()):
print('%-20s'%key,' ',inp[key])
return inp
class check_progress(object):
"""A class that will print a simple status bar of percentage of work done"""
def __init__(self, inp, T, skip_progress=False):
self.skip_print = skip_progress
if self.skip_print:
return
self.start_time = time.time()
self.n_steps = inp['n_steps']
self.n_check = 1000
if 'check_per_steps' in inp:
self.n_check = inp['check_per_steps']
self.fmt_print = ('%ix%i (T=%.2f) steps: %%7i/%7i, %%5.1f%%%% '
'run time: %%8s est.time-to-go: %%8s'%
(inp['N'], inp['N'], T, self.n_steps))
self.n_called = -1 # will progress untill n_called = n_steps
self.check()
# print(self.fmt_print)
# if False:
# print('%ix%i IsingLattice. Finished %10i / %10i steps (%4.1f%%). '
# 'Started '+str(time.strftime('%m-%d %H:%M:%S')))
# dir_out += str(time.strftime("_%Y%m%d-%H%M%S"))
def check(self, final=False):
if self.skip_print:
return
self.n_called += 1
if not self.n_called % 1000 == 0 and not final:
return
ratio = float(self.n_called) / self.n_steps
if ratio == 0:
time_pass_str = '00:00:00'
est_str = 'n/a'
else:
time_pass = time.time() - self.start_time
time_pass_str = time.strftime('%H:%M:%S', time.gmtime(time_pass));
est_time = (1-ratio)*time_pass/ratio
# print('est_time ', est_time)
est_str = time.strftime('%H:%M:%S', time.gmtime(est_time))
if final:
print(self.fmt_print%(self.n_called, ratio*100., time_pass_str, 'done!'))
else:
print(self.fmt_print%(self.n_called, ratio*100., time_pass_str, est_str), end='\r')
# def check_progress(inp):
# '''Print the progress of the plot'''
# pass
def run_ising_lattice(inp, T_final, skip_print=False):
'''Run a 2-D Ising model on a lattice.
Return three objects (each a numpy.array, which you can treat
as identical to a numpy.array)
M_avg: the average magnetization of each site for each n_analyze step
E_avg: '.........' energy '..................................'
correlation: an (N/2-1) array of the correlation function values at the final data frame
Note that this will use the generators from functions:
(1) make_T_generator
(2) make_B_generator
'''
time_start = time.time()
lattice = None
if inp['use_cpp'] and has_cpp:
lattice = IsingLattice_cpp(inp['N'], inp['flip_perc'])
elif inp['use_cpp'] and not has_cpp:
print('Warning: although use_cpp is set to 1, '
' the shared library IsingLattice.so is not present.\n '
' Therefore, the python implementation of Ising will be used.')
lattice = IsingLattice_py(inp['N'],inp['flip_perc'])
else:
lattice = IsingLattice_py(inp['N'],inp['flip_perc'])
# Make the run loop
try: # try loop that can be interrupted by the user
#first loop through all steps up to n_analyze
T_generator = make_T_generator(inp, T_final)
B_generator = make_B_generator(inp, T_final)
n_prior = inp['n_steps'] - inp['n_analyze']
progress = check_progress(inp, T_final, skip_print)
for T, B, step in zip(T_generator, B_generator, range(n_prior)):
lattice.step(T,B)
progress.check()
# loop through the analyze section of generators
E_avg = []
M_avg = []
for T, B, step in zip(T_generator, B_generator, range(inp['n_analyze'])):
lattice.step(T,B)
E_avg.append(lattice.get_E())
M_avg.append(lattice.get_M())
progress.check()
progress.check(True)
spin_correlation = np.array(lattice.calc_auto_correlation())
lattice.free_memory()
return (
np.array(E_avg),
np.array(M_avg),
np.array(spin_correlation)
)
except KeyboardInterrupt:
try:
lattice.free_memory()
except:
pass
print("\n\nProgram terminated by keyboard. Good Bye!")
sys.exit()
def plot_graphs(data): #T,E_mean,E_std,M_mean,M_std): #plot graphs at end
dat = np.array(data)
# print('data: ', dat)
# print('x: ', dat[:,1])
# x = dat[:,1][0]
# print('xlist: ', x)
plt.figure(1)
# plt.ylim(0,1)
plt.errorbar(dat[:,0], dat[:,1], yerr=dat[:,2], fmt='o')
plt.xlabel('Temperature')
plt.ylabel('Average Site Energy')
plt.figure(2)
plt.errorbar(dat[:,0], np.absolute(dat[:,3]), yerr=dat[:,4], uplims=True, lolims=True,fmt='o')
plt.xlabel('Temperature')
plt.ylabel('Aveage Site Magnetization')
plt.show()
def get_filenames(inp): #make data folder if doesn't exist, then specify filename
'''Generate the output file names for the EM (energy and megnetism) and SC (spin correlation) files'''
try:
dir_out = inp['dir_out']
prefix = inp['file_prefix']
if inp['date_output']:
dir_out += str(time.strftime("_%Y%m%d-%H%M%S"))
if not os.path.isdir(dir_out):
os.makedirs(dir_out)
# file name = [file_prefix]##_EM_v#.csv if only one temperature (example: runA_4.20_EM_v0.csv)
# [file_prefix]##T##_EM_v#.csv if there are two temperatures (example: runA_4.2T5.3_EM_v0.csv)
# the other file name is identical, but with "SC" (for spin correlation)) instead of EM
if inp['t_max'] <= inp['t_min']:
t_name = '%.2f'%inp['t_min']
else:
t_name = '%.2fT%.2f'%(inp['t_min'],inp['t_max'])
# print('%s%s_SC_v%i.csv'%(prefix,t_name,v))
v = 0
while (os.path.isfile( os.path.join(dir_out, '%s%s_EM_v%i.csv'%(prefix,t_name,v))) or
os.path.isfile( os.path.join(dir_out, '%s%s_SC_v%i.csv'%(prefix,t_name,v)))):
v += 1
return ( os.path.join(dir_out, '%s%s_EM_v%i.csv'%(prefix,t_name,v)),
os.path.join(dir_out, '%s%s_SC_v%i.csv'%(prefix,t_name,v)) )
except:
print ('fatal: Failed to make output file names')
sys.exit()
def print_results(inp, data, corr):
data_filename, corr_filename = get_filenames(inp)
with open(data_filename,'w') as f_out:
writer = csv.writer(f_out, delimiter=',', lineterminator='\n')
writer.writerow(['N', 'n_steps', 'n_analyze', 'flip_perc'])
writer.writerow([inp['N'], inp['n_steps'], inp['n_analyze'], inp['flip_perc']])
writer.writerow([])
writer.writerow(['Temp','E_mean','E_std','M_mean','M_std'])
for entry in data:
writer.writerow(entry)
# for t, e_mean, e_std, m_mean, m_std in zip(T, E_mean, E_std, M_mean, M_std):
# writer.writerow([t, e_mean, e_std, m_mean, m_std])
with open(corr_filename,'w') as f_out:
writer = csv.writer(f_out, delimiter=',', lineterminator='\n')
writer.writerow(['N', 'n_steps', 'n_analyze', 'flip_perc'])
writer.writerow([inp['N'], inp['n_steps'], inp['n_analyze'], inp['flip_perc']])
writer.writerow([])
writer.writerow(['Temp']+['d=%i'%i for i in range(1,len(corr[0])+1)])
for entry in corr:
writer.writerow(entry)
def run_indexed_process( inp, T, data_listener):
# def run_simulation(
# temp, n, num_steps, num_burnin, num_analysis, flip_prop, j, b, data_filename, corr_filename, data_listener, corr_listener):
print("Starting Temp {0}".format(round(T,3)))
try:
E, M, C = run_ising_lattice(inp, T, skip_print=True)
data_listener.put(([T,E.mean(),E.std(), M.mean(), M.std()], [T,]+[x[1] for x in C]))
# corr_listener.put([T,]+[x[1] for x in C])
print("Finished Temp {0}".format(round(T,3)))
return True
except KeyboardInterrupt:
print("\n\nProgram Terminated. Good Bye!")
data_listener.put('kill')
# corr_listener.put('kill')
sys.exit()
except:
logging.error("Temp="+str(round(T,3))+": Simulation Failed. No Data Written")
return False
def listener(queue, inp, data):
'''listen for messages on the queue
appends messages to data'''
# f = open(fn, 'a')
# writer = csv.writer(f, delimiter=',', lineterminator='\n')
while True:
message = queue.get()
# print('message: ', message)
if message == 'kill':
data['data'].sort()
data['corr'].sort()
print_results(inp, data['data'], data['corr'])
print ('Closing listener')
# print('killing')
break
data['data'].append(message[0])
data['corr'].append(message[1])
# print('--------\n',data)
def make_T_array(inp):
if inp['t_max'] <= inp['t_min']:
return [inp['t_min'],]
else:
return np.arange(inp['t_min'], inp['t_max'], inp['t_step'])
def run_multi_core(inp):
print("\n2D Ising Model Simulation; multi-core\n")
T_array = make_T_array(inp)
#must use Manager queue here, or will not work
manager = mp.Manager()
data_listener = manager.Queue()
# corr_listener = manager.Queue()
pool = mp.Pool(mp.cpu_count())
# arrays of results:
data = {'data':[], 'corr':[]}
# corr = []
#put listener to work first
data_watcher = pool.apply_async(listener, args=(data_listener, inp, data,))
# corr_watcher = pool.apply_async(listener, args=(corr_listener, inp, corr,))
#fire off workers
jobs = [pool.apply_async(run_indexed_process,args=(inp,T, data_listener)) for T in T_array]
# collect results from the workers through the pool result queue
[job.get() for job in jobs]
data_listener.put('kill')
pool.close()
def run_single_core(inp):
print("\n2D Ising Model Simulation; single core\n")
# sequentially run through the desired temperatures and collect the output for each temperature
data = []
corr = []
for temp in make_T_array(inp):
E, M, C = run_ising_lattice(inp, temp, skip_print=inp['skip_prog_print'])
data.append( (temp, E.mean(), E.std(), M.mean(), M.std() ) )
corr.append([temp,]+[x[1] for x in C])
print_results(inp, data, corr)
if inp['plots']:
plot_graphs(data)
if __name__ == "__main__":
"""Main program: run Ising Lattice here"""
inp = set_input(argv)
if inp['multiprocess']:
run_multi_core(inp)
if not inp['multiprocess']:
run_single_core(inp)