14 from scipy.interpolate
import interp1d
15 import matplotlib.pyplot
as plt
20 p = optparse.OptionParser()
21 p.add_option(
"-i",
"--input" , action=
"store", dest=
"rundir", default=
'.', \
22 help=
"Simulation directory to summarize (default: current directory)")
23 p.add_option(
"-o",
"--output", action=
"store", dest=
"outdir", default=
'./summary.hdf5', \
24 help=
"Output path (default: summary.hdf5)")
25 p.add_option(
"-b",
"--buf" , action=
"store", dest=
"buffersize", default=
'500', \
26 help=
"Buffer size before writing to the file (default: 500)")
27 p.add_option(
"-f",
"--force", action=
"store_true", dest=
"force", default=
False, \
28 help=
"Force overwrite of the output file (default: False)")
29 p.add_option(
"-v",
"--verbose", action=
"store_true", dest=
"verbose", default=
False, \
30 help=
"Enable verbose output. If enabled, output is written to 'debug.log' (default: False)")
31 p.add_option(
"--time_file", action=
"store", dest=
"time_file", default=
None, \
32 help=
"File containing a time grid in seconds to map the data to (default: None)")
33 p.add_option(
"--time_final", action=
"store", dest=
"time_final", default=
None, \
34 help=
"Final time for the time grid. Only used if '--time_file' "+\
35 "is not given (default: Read from template file)")
36 p.add_option(
"--time_initial", action=
"store", dest=
"time_initial", default=
None, \
37 help=
"Initial time for the time grid. Only used if '--time_file' "+\
38 "is not given (default: 1e-5)")
39 p.add_option(
"--time_number", action=
"store", dest=
"time_number", default=
None, \
40 help=
"Number of time steps for the time grid. Only used if '--time_file' "+\
41 "is not given (default: 200)")
42 p.add_option(
"--sunet_path", action=
"store", dest=
"sunet_path", default=
None, \
43 help=
"Path to the sunet file (default: Read from template)")
44 p.add_option(
"--disable_mainout", action=
"store_true", dest=
"disable_mainout", default=
False, \
45 help=
"Disable the summary of the mainout output (default: False)")
46 p.add_option(
"--disable_energy", action=
"store_true", dest=
"disable_energy", default=
False, \
47 help=
"Disable the summary of the energy output (default: False)")
48 p.add_option(
"--disable_timescales", action=
"store_true", dest=
"disable_timescales", default=
False, \
49 help=
"Disable the summary of the timescales output (default: False)")
50 p.add_option(
"--disable_tracked_nuclei", action=
"store_true", dest=
"disable_tracked_nuclei", default=
False, \
51 help=
"Disable the summary of the tracked_nuclei output (default: False)")
52 p.add_option(
"--disable_nuloss", action=
"store_true", dest=
"disable_nuloss", default=
False, \
53 help=
"Disable the summary of the nuloss output (default: False)")
54 p.add_option(
"--disable_snapshots", action=
"store_true", dest=
"disable_snapshots", default=
False, \
55 help=
"Disable the summary of the snapshots output (default: False)")
56 p.add_option(
"-r",
"--recursive", action=
"store_true", dest=
"recursive", default=
False, \
57 help=
"Give a folder that contains subfolders with runs inside (default: False)")
60 Usage: ./summarize.py -i <rundir>
61 Example: ./summarize.py -i runs/test""")
65 (options,args) = p.parse_args()
66 run_path = options.rundir
70 buffsize = int(options.buffersize)
77 format=
"%(asctime)s - %(levelname)-10s - %(message)s",
79 datefmt=
"%Y-%m-%d %H:%M",
82 file_handler = logging.FileHandler(
"debug.log", mode=
"w", encoding=
"utf-8")
83 logger = logging.getLogger(__name__)
85 file_handler.setFormatter(logging.Formatter(
86 fmt=
"%(asctime)s - %(levelname)s - %(message)s",
87 datefmt=
"%Y-%m-%d %H:%M",
90 logger.addHandler(file_handler)
93 logging.basicConfig(level=logging.ERROR)
94 logger = logging.getLogger(__name__)
97 logger.propagate =
False
99 logger.info(f
"Started summarizing run at {run_path}.")
102 if options.recursive:
104 folders = [f
for f
in os.listdir(run_path)
if os.path.isdir(os.path.join(run_path, f))]
106 logger.info(f
"Recursively going through {len(folders)} folders.")
108 fold = options.outdir.replace(
".hdf5",
"")
110 if not os.path.exists(fold):
114 if not options.force:
115 logger.error(f
"Folder {fold} already exists. Exiting.")
116 raise ValueError(f
"Folder {fold} already exists. Exiting.")
118 logger.warning(f
"Folder {fold} already exists. Overwriting it.")
119 os.system(f
"rm -r {fold}")
127 if options.recursive:
130 output_file = os.path.join(fold, folders[0]+
".hdf5")
131 run_path = os.path.join(basepath, folders.pop(0))
137 output_file = options.outdir
140 dirs = [d
for d
in os.listdir(run_path)
if os.path.isdir(os.path.join(run_path, d))
and d !=
"network_data"]
143 logger.info(f
"Found {len(dirs)} directories in {run_path}.")
147 if os.path.exists(output_file)
and not options.force:
148 logger.error(f
"Output file {output_file} already exists. Exiting.")
149 raise ValueError(f
"Output file {output_file} already exists. Either delete or use -f option to overwrite.")
150 elif os.path.exists(output_file)
and options.force:
151 logger.warning(f
"Output file {output_file} already exists. Overwriting it.")
152 f_hdf = h5py.File(output_file,
'w')
160 if not data.is_crashed:
161 logger.info(f
"Found {d} to look up the output.")
166 if options.recursive:
167 logger.warning(
"Detected folder where all runs have crashed!")
169 logger.error(
"All runs are crashed!")
170 raise ValueError(
"All runs are crashed!")
173 template_name = [f
for f
in os.listdir(os.path.join(run_path, d))
if f.endswith(
'.par')][0]
174 t =
template(os.path.join(run_path, d, template_name))
177 if options.sunet_path
is not None:
178 net_source = options.sunet_path
180 if (
not "net_source" in t.entries):
182 raise ValueError(
"net_source not given in the template file.")
185 net_source = t[
"net_source"]
188 logger.info(f
"Using sunet file from {net_source}.")
189 nuclei = np.loadtxt(net_source,dtype=str)
194 if options.time_file
is not None:
195 logger.info(f
"Using time grid from {options.time_file}.")
196 mainout_time = np.loadtxt(options.time_file, dtype=float, unpack=
True)
198 if options.time_final
is not None:
199 final_time = float(options.time_final)
204 if "termination_criterion" in t.entries:
205 if t[
"termination_criterion"] ==
"1":
206 if "final_time" in t.entries:
207 final_time = float(t[
"final_time"])
209 if options.time_initial
is not None:
210 initial_time = float(options.time_initial)
214 if options.time_number
is not None:
215 time_number = int(options.time_number)
219 mainout_time = np.logspace(np.log10(initial_time), np.log10(final_time), time_number)
223 possible_entries = []
224 if not options.disable_mainout:
225 possible_entries.append(
"mainout")
227 logger.info(
"Ignoring mainout output.")
228 if not options.disable_energy:
229 possible_entries.append(
"energy")
231 logger.info(
"Ignoring energy output.")
232 if not options.disable_timescales:
233 possible_entries.append(
"timescales")
235 logger.info(
"Ignoring timescales output.")
236 if not options.disable_tracked_nuclei:
237 possible_entries.append(
"tracked_nuclei")
239 logger.info(
"Ignoring tracked_nuclei output.")
240 if not options.disable_nuloss:
241 possible_entries.append(
"nuloss")
243 logger.info(
"Ignoring nuloss output.")
247 for entry
in possible_entries:
249 if data.check_existence(entry) != 0:
250 entry_dict[entry] = {}
251 for key
in data[entry].keys():
253 if ((key ==
"iteration")
or (key ==
"time")
or (key ==
"A")
or (key ==
"Z")
or (key ==
"N")
254 or (key ==
"names")
or (key ==
"latex_names")):
257 if (entry ==
"nuloss")
and ((key ==
"temp")
or (key ==
"dens")
or (key ==
"rad")):
260 if data[entry][key].ndim == 1:
261 entry_dict[entry][key] = np.zeros((len(mainout_time),buffsize))
262 elif data[entry][key].ndim == 2:
263 entry_dict[entry][key] = np.zeros((len(mainout_time),data[entry][key].shape[1],buffsize))
265 raise ValueError(f
"Invalid shape of {entry}/{key} in the data.")
268 f_hdf[entry+
"/time"] = mainout_time
271 logger.info(f
"Found {entry} in the data.")
282 if (data.check_existence(
"snapshot") != 0)
and (
not options.disable_snapshots):
283 if (
"custom_snapshots" in t.entries)
or (
"h_custom_snapshots" in t.entries):
285 summarize_snapshots =
False
286 if (
"custom_snapshots" in t.entries):
287 summarize_snapshots = (t[
"custom_snapshots"].strip().lower() ==
"yes")
289 if (t[
"custom_snapshots"].strip().lower() ==
"yes"):
290 logger.info(
"Found custom snapshots in ascii format.")
291 if (
"h_custom_snapshots" in t.entries):
292 summarize_snapshots = (summarize_snapshots
or (t[
"h_custom_snapshots"].strip().lower() ==
"yes"))
294 if (t[
"h_custom_snapshots"].strip().lower() ==
"yes"):
295 logger.info(
"Found custom snapshots in hdf5 format.")
298 if summarize_snapshots:
299 if "snapshot_file" not in t.entries:
300 raise ValueError(
"Invalid template file. snapshot_file not given in the template file.")
301 snapshot_time = np.loadtxt(os.path.join(t[
"snapshot_file"]),dtype=float)
303 snapshot_time *= 24*3600
305 f_hdf[
"snapshots/time"] = snapshot_time
307 f_hdf[
"snapshots/A"] = nuclei_data.A
308 f_hdf[
"snapshots/Z"] = nuclei_data.Z
309 f_hdf[
"snapshots/N"] = nuclei_data.N
311 snapshot_data = np.zeros((len(nuclei),len(snapshot_time),buffsize))
312 logger.info(f
"Summarize custom snapshots as well.")
314 summarize_snapshots =
False
316 summarize_snapshots =
False
317 if not options.disable_snapshots:
318 logger.info(
"Ignoring snapshots output.")
322 finab_data_Y = np.zeros((len(nuclei),buffsize))
323 finab_data_X = np.zeros((len(nuclei),buffsize))
326 f_hdf[
"finab/A"] = nuclei_data.A
327 f_hdf[
"finab/Z"] = nuclei_data.Z
328 f_hdf[
"finab/N"] = nuclei_data.N
331 dtype_nuclei = np.dtype([(
'A', int), (
'Z', int)])
332 nuclei_struct = np.array(list(zip(nuclei_data.A.astype(int), nuclei_data.Z.astype(int))), dtype=dtype_nuclei)
334 nuclei_sorted_idx = np.argsort(nuclei_struct)
335 sorted_nuclei_struct = nuclei_struct[nuclei_sorted_idx]
339 run_names = np.zeros(buffsize,dtype=
"S100")
340 run_ids = np.zeros(buffsize,dtype=int)
344 for counter, d
in enumerate(tqdm(dirs)):
346 data =
wreader(os.path.join(run_path, d),silent=
True)
349 logger.warning(f
"Run {d} is crashed. Skipping it.")
365 if ind % buffsize == 0
and ind != 0:
367 if "finab/Y" not in f_hdf:
368 f_hdf.create_dataset(
"finab/Y", (len(nuclei),ind+1), maxshape=(len(nuclei),
None))
369 f_hdf.create_dataset(
"finab/X", (len(nuclei),ind+1), maxshape=(len(nuclei),
None))
372 f_hdf[
"finab/Y"].resize((len(nuclei),ind+1))
373 f_hdf[
"finab/X"].resize((len(nuclei),ind+1))
375 f_hdf[
"finab/Y"][:,ind-buffsize:ind] = finab_data_Y
376 f_hdf[
"finab/X"][:,ind-buffsize:ind] = finab_data_X
379 finab_struct = np.array(list(zip(data.finab[
"A"].astype(int), data.finab[
"Z"].astype(int))), dtype=dtype_nuclei)
382 matching_idx = np.searchsorted(sorted_nuclei_struct, finab_struct)
384 indices = nuclei_sorted_idx[matching_idx]
387 finab_data_Y[:,ind % buffsize] = 0
388 finab_data_X[:,ind % buffsize] = 0
389 finab_data_Y[indices,ind % buffsize] = data.finab[
"Y"][:]
390 finab_data_X[indices,ind % buffsize] = data.finab[
"X"][:]
396 if ind % buffsize == 0
and ind != 0:
398 if "run_names" not in f_hdf:
399 f_hdf.create_dataset(
"run_names", (ind+1,), maxshape=(
None,),dtype=
"S100")
402 f_hdf[
"run_names"].resize((ind+1,))
404 f_hdf[
"run_names"][ind-buffsize:ind] = run_names
407 run_names[ind % buffsize] = d
412 if ind % buffsize == 0
and ind != 0:
414 if "run_ids" not in f_hdf:
415 f_hdf.create_dataset(
"run_ids", (ind+1,), maxshape=(
None,),dtype=int)
418 f_hdf[
"run_ids"].resize((ind+1,))
420 f_hdf[
"run_ids"][ind-buffsize:ind] = run_ids
423 run_ids[ind % buffsize] = int(re.findall(
r'\d+', d)[-1])
425 run_ids[ind % buffsize] = -1
430 if summarize_snapshots:
432 snapstime = data.snapshot_time
435 indexes = np.searchsorted(snapstime, snapshot_time)
438 if (len(snapstime) < len(snapshot_time)):
439 mask = np.zeros(len(snapshot_time),dtype=bool)
441 mask[np.min(np.abs(snapstime - snapshot_time[:,np.newaxis]),axis=1) > 1e-5] =
True
447 finab_struct = np.array(list(zip(data.A.astype(int), data.Z.astype(int))), dtype=dtype_nuclei)
450 matching_idx = np.searchsorted(sorted_nuclei_struct, finab_struct)
452 indices_nuclei = nuclei_sorted_idx[matching_idx]
455 if ind % buffsize == 0
and ind != 0:
457 if "snapshots/Y" not in f_hdf:
458 f_hdf.create_dataset(
"snapshots/Y", (len(nuclei),len(snapshot_time),ind+1), maxshape=(len(nuclei),len(snapshot_time),
None))
459 f_hdf.create_dataset(
"snapshots/X", (len(nuclei),len(snapshot_time),ind+1), maxshape=(len(nuclei),len(snapshot_time),
None))
462 f_hdf[
"snapshots/Y"].resize((len(nuclei),len(snapshot_time),ind+1))
463 f_hdf[
"snapshots/X"].resize((len(nuclei),len(snapshot_time),ind+1))
465 f_hdf[
"snapshots/Y"][:,:,ind-buffsize:ind] = snapshot_data
466 f_hdf[
"snapshots/X"][:,:,ind-buffsize:ind] = snapshot_data*nuclei_data.A[:,np.newaxis,np.newaxis]
469 snapshot_data[:,:,ind % buffsize] = 0
470 snapshot_data[indices_nuclei,:,ind % buffsize] = data.Y[indexes][:, :].T
474 snapshot_data[:,mask,ind % buffsize] = np.nan
480 for entry
in entry_dict.keys():
482 if ind % buffsize == 0
and ind != 0:
483 for key
in entry_dict[entry].keys():
485 dim = entry_dict[entry][key].ndim
488 if entry+
"/"+key
not in f_hdf:
489 f_hdf.create_dataset(entry+
"/"+key, (len(mainout_time),ind+1), maxshape=(len(mainout_time),
None))
492 f_hdf[entry+
"/"+key].resize((len(mainout_time),ind+1))
494 f_hdf[entry+
"/"+key][:,ind-buffsize:ind] = entry_dict[entry][key]
496 if entry+
"/"+key
not in f_hdf:
497 f_hdf.create_dataset(entry+
"/"+key, (len(mainout_time),data[entry][key].shape[1],ind+1), maxshape=(len(mainout_time),data[entry][key].shape[1],
None))
500 f_hdf[entry+
"/"+key].resize((len(mainout_time),data[entry][key].shape[1],ind+1))
502 f_hdf[entry+
"/"+key][:,:,ind-buffsize:ind] = entry_dict[entry][key]
505 for key
in entry_dict[entry].keys():
506 value =
interp1d(data[entry][
"time"],data[entry][key],
507 bounds_error =
False, fill_value = np.nan, axis = 0)(mainout_time)
508 if entry_dict[entry][key].ndim == 2:
509 entry_dict[entry][key][:,ind % buffsize] = value
510 elif entry_dict[entry][key].ndim == 3:
511 entry_dict[entry][key][:,:,ind % buffsize] = value
516 logger.info(f
"Finished looping over all directories, writing the last data to the hdf5 file.")
524 if "finab/Y" not in f_hdf:
525 f_hdf.create_dataset(
"finab/Y", (len(nuclei),ind+1), maxshape=(len(nuclei),
None))
526 f_hdf.create_dataset(
"finab/X", (len(nuclei),ind+1), maxshape=(len(nuclei),
None))
528 f_hdf[
"finab/Y"].resize((len(nuclei),ind+1))
529 f_hdf[
"finab/X"].resize((len(nuclei),ind+1))
532 f_hdf[
"finab/Y"][:,ind-buffsize+1:ind+1] = finab_data_Y[:,:buffsize]
533 f_hdf[
"finab/X"][:,ind-buffsize+1:ind+1] = finab_data_X[:,:buffsize]
535 f_hdf[
"finab/Y"][:,:ind+1] = finab_data_Y[:,:ind+1]
536 f_hdf[
"finab/X"][:,:ind+1] = finab_data_X[:,:ind+1]
542 if "run_names" not in f_hdf:
543 f_hdf.create_dataset(
"run_names", (ind+1,), maxshape=(
None,),dtype=
"S100")
545 f_hdf[
"run_names"].resize((ind+1,))
548 f_hdf[
"run_names"][ind-buffsize+1:ind+1] = run_names[:buffsize]
550 f_hdf[
"run_names"][:ind+1] = run_names[:ind+1]
553 if "run_ids" not in f_hdf:
554 f_hdf.create_dataset(
"run_ids", (ind+1,), maxshape=(
None,),dtype=int)
556 f_hdf[
"run_ids"].resize((ind+1,))
559 f_hdf[
"run_ids"][ind-buffsize+1:ind+1] = run_ids[:buffsize]
561 f_hdf[
"run_ids"][:ind+1] = run_ids[:ind+1]
567 if summarize_snapshots:
568 if "snapshots/Y" not in f_hdf:
569 f_hdf.create_dataset(
"snapshots/Y", (len(nuclei),len(snapshot_time),ind+1), maxshape=(len(nuclei),len(snapshot_time),
None))
570 f_hdf.create_dataset(
"snapshots/X", (len(nuclei),len(snapshot_time),ind+1), maxshape=(len(nuclei),len(snapshot_time),
None))
572 f_hdf[
"snapshots/Y"].resize((len(nuclei),len(snapshot_time),ind+1))
573 f_hdf[
"snapshots/X"].resize((len(nuclei),len(snapshot_time),ind+1))
576 f_hdf[
"snapshots/Y"][:,:,ind-buffsize+1:ind+1] = snapshot_data[:,:,:buffsize]
577 f_hdf[
"snapshots/X"][:,:,ind-buffsize+1:ind+1] = snapshot_data[:,:,:buffsize]*nuclei_data.A[:,np.newaxis,np.newaxis]
579 f_hdf[
"snapshots/Y"][:,:,:ind+1] = snapshot_data[:,:,:ind+1]
580 f_hdf[
"snapshots/X"][:,:,:ind+1] = snapshot_data[:,:,:ind+1]*nuclei_data.A[:,np.newaxis,np.newaxis]
586 for entry
in entry_dict.keys():
587 for key
in entry_dict[entry].keys():
589 dim = entry_dict[entry][key].ndim
592 if entry+
"/"+key
not in f_hdf:
593 f_hdf.create_dataset(entry+
"/"+key, (len(mainout_time),ind+1), maxshape=(len(mainout_time),
None))
595 f_hdf[entry+
"/"+key].resize((len(mainout_time),ind+1))
598 f_hdf[entry+
"/"+key][:,ind-buffsize+1:ind+1] = entry_dict[entry][key][:,:buffsize]
600 f_hdf[entry+
"/"+key][:,:ind+1] = entry_dict[entry][key][:,:ind+1]
602 if entry+
"/"+key
not in f_hdf:
603 f_hdf.create_dataset(entry+
"/"+key, (len(mainout_time),data[entry][key].shape[1],ind+1), maxshape=(len(mainout_time),data[entry][key].shape[1],
None))
605 f_hdf[entry+
"/"+key].resize((len(mainout_time),data[entry][key].shape[1],ind+1))
608 f_hdf[entry+
"/"+key][:,:,ind-buffsize+1:ind+1] = entry_dict[entry][key][:,:,:buffsize]
610 f_hdf[entry+
"/"+key][:,:,:ind+1] = entry_dict[entry][key][:,:,:ind+1]
615 logger.info(f
"Finished summarizing run at {run_path}.")
620 if not options.recursive: