Skip to content
Snippets Groups Projects
Commit cf1d5b33 authored by Doak, Peter W.'s avatar Doak, Peter W.
Browse files

working plot miniSofQW.py and CIINT sample script improvements

parent be51515a
No related branches found
No related tags found
No related merge requests found
......@@ -16,16 +16,24 @@ import numpy as np
import re
import matplotlib as mpl
import matplotlib.pyplot as plt
import argparse
if __name__ == '__main__':
file_name = "single_meas_G_100B_10k.bp"
parser = argparse.ArgumentParser(description="read and compare G_k_w (and M_k_w) top miniapp produced values.")
parser.add_argument('--file', dest="file_name", action="store", default="single_meas_G_100B_10k.bp")
parser.add_argument('--force-reload', dest="force_reload", action="store_true", default=False)
input_args = parser.parse_args()
file_name = input_args.file_name
fig = plt.figure(constrained_layout=True, figsize=(20,10))
subfigs = fig.subfigures(1,3, width_ratios=[.3,1,1])
num_samples = 10
rank = 0
adios_fh = pAr.fhFromFilename(file_name)
dca_globals = pAr.readDCARun(adios_fh)
pickle_name = input_args.file_name + "_pars.pickle"
dca_globals = pAr.tryPickleUnlessForce(input_args.force_reload, pickle_name, pAr.readDCARun, adios_fh)
parameters = [(key, item) for key,item in dca_globals.items() if re.search("parameters", key)]
for param in parameters:
print("param", param)
......@@ -42,17 +50,18 @@ if __name__ == '__main__':
axsGreal = subfigs[1].subplots(2,1,squeeze=True)
axsGimag = subfigs[2].subplots(2,1,squeeze=True)
samp_seq = [index * 100 + 1 for index in range(num_samples)]
list_meas = pAr.readSampleConfigG_k_w(adios_fh,rank,0,samp_seq)
for index in range(num_samples):
samp = index * 100 + 1
rank = 0
axsv = axs_vertexes[index]
# axsv.set_xlim(0,100)
# axsv.set_xticks([0,50,100])
# axsv.axis["y=0"] = axsv.new_floating_axis(nth_coord=0, value=0, axis_direction="bottom")
# axsv.axis["y=0"].toggle(all=True)
# axsv.axis["bottom","top","right"].set_visible(False)
G_k_w, log_weight, aux_spins, vertex_times = pAr.readSampleConfigG_k_w(adios_fh,rank,0,samp)
samp = samp_seq[index]
G_k_w, log_weight, aux_spins, vertex_times = list_meas[index]
upTimes = [ t for [s,t] in zip(aux_spins, vertex_times) if s == 1]
downTimes = [ t for [s,t] in zip(aux_spins, vertex_times) if s == 0]
......
......@@ -31,77 +31,84 @@ def trim_axs(axs, N):
ax.remove()
return axs[:N]
if __name__ == '__main__':
parser = argparse.ArgumentParser(description="read and compare G_k_w (and M_k_w) top miniapp produced values.")
parser.add_argument('--file', dest="file_name", action="store", default="single_meas_G_100B_10k.bp")
parser.add_argument('--M_k_w', dest="include_M", action="store_true", default=False)
input_args = parser.parse_args()
parser = argparse.ArgumentParser(description="read and compare G_k_w (and M_k_w) top miniapp produced values.")
parser.add_argument('--file', dest="file_name", action="store", default="single_meas_G_100B_10k.bp")
parser.add_argument('--M_k_w', dest="include_M", action="store_true", default=False)
parser.add_argument('--force-reload', dest="force_reload", action="store_true", default=False)
input_args = parser.parse_args()
if input_args.include_M:
fig = plt.figure(figsize=(20,10), constrained_layout=True)
samp_axs = trim_axs(fig.subplots(4,5),20)
else:
fig = plt.figure(figsize=(20,5), constrained_layout=True)
samp_axs = trim_axs(fig.subplots(2,5),10)
if input_args.include_M:
fig = plt.figure(figsize=(20,10), constrained_layout=True)
samp_axs = trim_axs(fig.subplots(4,5),20)
else:
fig = plt.figure(figsize=(20,5), constrained_layout=True)
samp_axs = trim_axs(fig.subplots(2,5),10)
adios_fh = pAr.fhFromFilename(input_args.file_name)
dca_globals = pAr.readDCARun(adios_fh)
parameters = [(key, item) for key,item in dca_globals.items() if re.search("parameters", key)]
for param in parameters:
print("param", param)
g0_r_t = dca_globals['cluster_excluded_greens_function_G0_r_t']
g0_k_w = dca_globals['cluster_excluded_greens_function_G0_k_w']
g0_r_w = dca_globals['cluster_excluded_greens_function_G0_r_w']
wn = dca_globals['frequency_domain_elements']
beta = dca_globals['physics_beta'][0]
density = dca_globals['physics_density'][0]
U = dca_globals['single_band_Hubbard_model_U'][0]
samp_seq = [ x * 90 + 1 for x in range(10) ]
for samp, ax in zip(samp_seq, samp_axs):
if input_args.include_M:
G_k_w, M_k_w, log_weight, aux_spins, vertex_times, sign = pAr.readSampleConfigM_k_w(adios_fh,0,0,samp)
else:
G_k_w, log_weight, aux_spins, vertex_times = pAr.readSampleConfigG_k_w(adios_fh,0,0,samp)
adios_fh = pAr.fhFromFilename(input_args.file_name)
pickle_name = input_args.file_name + "_pars.pickle"
dca_globals = pAr.tryPickleUnlessForce(input_args.force_reload, pickle_name, pAr.readDCARun, adios_fh)
parameters = [(key, item) for key,item in dca_globals.items() if re.search("parameters", key)]
g0_r_t = dca_globals['cluster_excluded_greens_function_G0_r_t']
g0_k_w = dca_globals['cluster_excluded_greens_function_G0_k_w']
g0_r_w = dca_globals['cluster_excluded_greens_function_G0_r_w']
wn = dca_globals['frequency_domain_elements']
beta = dca_globals['physics_beta'][0]
density = dca_globals['physics_density'][0]
U = dca_globals['single_band_Hubbard_model_U'][0]
samp_seq = [ x * 90 + 1 for x in range(len(samp_axs))]
if input_args.include_M:
list_meas = pAr.readSampleConfigM_k_w(adios_fh,0,0,samp_seq)
else:
list_meas = pAr.readSampleConfigG_k_w(adios_fh,0,0,samp_seq)
mini_GwnUp, mini_GwnDn, mini_MwnUp, mini_MwnDn = CTINT.calculateG_single(CTINT.Sample(vertex_times, aux_spins, log_weight, G_k_w), g0_r_t, g0_k_w, U, beta, wn)
GdcaUp = G_k_w[:,0,1,0,1,0]
GdcaDn = G_k_w[:,0,0,0,0,0]
upTimes = [ t for [s,t] in zip(aux_spins, vertex_times) if s == 1]
downTimes = [ t for [s,t] in zip(aux_spins, vertex_times) if s == 0]
for i, ax in enumerate(samp_axs):
if input_args.include_M:
G_k_w, M_k_w, log_weight, aux_spins, vertex_times, sign = list_meas[i]
else:
G_k_w, log_weight, aux_spins, vertex_times = list_meas[i]
print ("max time: ",max(upTimes))
# plot M's
# mmUp = M_k_w[:,0,1,0,1,0]
# mmDn = M_k_w[:,0,0,0,0,0]
#mini_GwnUp, mini_GwnDn, mini_MwnUp, mini_MwnDn = CTINT.calculateG_single(CTINT.Sample(vertex_times, aux_spins, log_weight, G_k_w), g0_r_t, g0_k_w, U, beta, wn)
# DCAMGwnUp, DCAMGwnDn = CTINT.computeGFromG0kwAndMkw(g0_k_w, M_k_w, beta)
# #ax.set_ylim((-0.5, 0.3))
fig.suptitle("CTINT hubbard model U={} density={} beta={}".format(U, density, beta))
ax.plot(wn,GdcaUp.real, linewidth="0.5", color='maroon', label="up")
ax.plot(wn,GdcaDn.real, linewidth="0.5", color='orange', label="down")
ax.set_title("{} log wgt: {:01.3f}".format(samp,log_weight))
ax.legend(loc=3)
# local_up_line, = ax.plot(wn, mmUp.real , 'black')
# local_up_line.set_linewidth(1)
# local_up_line.set_dashes([1,1,5,1])
# ax.plot(wn, MwnDn.real, 'red')
# local_down_line, = ax.plot(wn, mmDn.real, 'grey')
GdcaUp = G_k_w[:,0,1,0,1,0]
GdcaDn = G_k_w[:,0,0,0,0,0]
upTimes = [ t for [s,t] in zip(aux_spins, vertex_times) if s == 1]
downTimes = [ t for [s,t] in zip(aux_spins, vertex_times) if s == 0]
ax.set_ylim((-0.5, 0.3))
fig.suptitle("CTINT hubbard model U={} density={} beta={}".format(U, density, beta))
ax.plot(wn,GdcaUp.real, linewidth="0.5", color='maroon', label="G_k_w up")
ax.plot(wn,GdcaDn.real, linewidth="0.5", color='orange', label="G_k_w down")
ax.set_title("{} log wgt: {:01.3f}".format(samp_seq[i],log_weight[0]))
ax.legend(loc=3)
# local_up_line.set_linewidth(1)
# local_up_line.set_dashes([1,1,5,1])
ax.set_xlim((0,10))
if input_args.include_M:
ax2 = ax.twinx()
mmUp = M_k_w[:,0,1,0,1,0]
mmDn = M_k_w[:,0,0,0,0,0]
local_up_line, = ax2.plot(wn, mmUp.real , 'black', linewidth="0.5", label="M_k_w real up")
local_down_line = ax2.plot(wn, mmDn.real, 'grey', linewidth="0.5", label="M_k_w real dn")
# local_down_line.set_linewidth(2)
# local_down_line.set_dashes([1,1,5,1])
# ax.set_xlim((-10,10))
# ax.set_title("Sample {}".format(samp))
# ax2 = ax.twinx()
# #ax2.set_ylim((-0.7, 0.7))
# ax2.plot(wn, MwnUp.imag, 'steelblue')
#ax2.set_ylim((-0.7, 0.7))
ax2.plot(wn, mmUp.imag, 'steelblue', linewidth="0.5", label="M_k_w imag up" )
ax2.plot(wn, mmDn.imag, 'dodgerblue', linewidth="0.5", label="M_k_w imag dn")
# local_up_line2, = ax2.plot(wn, mmUp.imag, 'black')
# local_up_line2.set_linewidth(1)
# local_up_line2.set_dashes([1,1,5,1])
# ax2.plot(wn, MwnDn.imag, 'dodgerblue')
# local_down_line2, = ax2.plot(wn, mmDn.imag, 'grey' )
# local_down_line2.set_linewidth(2)
# local_down_line2.set_dashes([1,1,5,1])
plt.show()
plt.show()
......@@ -57,7 +57,6 @@ if __name__ == '__main__':
print("bse.chiL.shape", bse.chiL.shape)
print("bse.chiC.shape", bse.chiL.shape)
chi0_cluster = bse.chic0
fig = plt.figure(figsize=(20,5), constrained_layout=True)
axs = trim_axs(fig.subplots(1,4), 4)
......
......@@ -48,14 +48,15 @@ def logplot_and_inset(resample_count, bin_sizes, ax, direct_bins, start_freq, en
if __name__ == '__main__':
parser = argparse.ArgumentParser(description="read and compare G_k_w (and M_k_w) top miniapp produced values.")
parser = argparse.ArgumentParser(description="read calculate via binning analysis the G_k_w or M_k_w error bars.")
parser.add_argument('--file-stem', dest="filestem", action="store", default="single_meas_G_100B_U4_10k_batch")
parser.add_argument('--start-freq', dest="start_freq", action="store", default=512, type=int)
parser.add_argument('--end-freq', dest="end_freq", action="store", default=1024, type=int)
parser.add_argument('--num-runs', dest="num_runs", action="store", default=20, type=int)
parser.add_argument('--force', dest="force_restrap", action="store", type=bool, default=False)
parser.add_argument('--do-direct', dest="do_direct", action="store_true", default=False)
parser.add_argument('--calculate-M-error', dest="m_error", action="store_true", default=False)
input_args = parser.parse_args()
filename = "{}_{}.bp".format(input_args.filestem,1)
......@@ -74,12 +75,25 @@ if __name__ == '__main__':
num_runs = input_args.num_runs
numpy_file_name = input_args.filestem + ".npy"
if input_args.m_error:
numpy_file_name = input_args.filestem + "_mfunc.npy"
pickle_name = input_args.filestem + "_mfunc.pickle"
else:
numpy_file_name = input_args.filestem + ".npy"
pickle_name = input_args.filestem + ".pickle"
procs = max(8, num_runs)
print("Preparing to read {} measuements.".format(total_measurements))
if os.path.isfile(numpy_file_name) and not input_args.force_restrap:
freq_hist_pool = np.load(numpy_file_name)
if os.path.isfile(pickle_name):
with open(pickle_name, 'rb') as the_pick:
freq_hist_info = pickle.load(the_pick)
else:
raise ValueError("pickle file {} with metadata for {} not found!".format(numpy_file_name, pickle_name))
else:
file_names = []
for id in range(num_runs):
......@@ -89,21 +103,34 @@ if __name__ == '__main__':
#freq_hist_pool_up_trans = freq_hist_pool_trans[0,:,:]
#freq_hist_pool = np.transpose(freq_hist_pool_up_trans)
# if the files are reasonably new data the 10, [10000]*10 are read from the file
freq_hist_pool_trans = ea.mergeSamplesSpinAvg(file_names, fq_start, fq_end, dca_globals, procs=procs)
# if the files are reasonably new data the num ranks and measurements per rank are read from the file
# otherwise they must be supplied
if input_args.m_error:
freq_hist_pool_trans, freq_hist_info = ea.mergeSamplesSpinAvg(file_names, fq_start, fq_end, dca_globals,
procs=procs, ml_function=ea.MLFunctions.MKW)
else:
freq_hist_pool_trans, freq_hist_info = ea.mergeSamplesSpinAvg(file_names, fq_start, fq_end, dca_globals, procs=procs)
freq_hist_pool = np.transpose(freq_hist_pool_trans)
np.save(numpy_file_name, freq_hist_pool)
print("Samples read and merged")
with open(pickle_name, 'wb') as the_pick:
pickle.dump(freq_hist_info, the_pick, pickle.HIGHEST_PROTOCOL)
print("Samples read and merged and saved to numpy archive")
fig = plt.figure(figsize=(10,5), constrained_layout=True)
total_mean = np.mean(freq_hist_pool, 0)
print("Data loaded hist pool shape", freq_hist_pool.shape)
std_Gls, results = ea.recursiveMakePow2Bins(freq_hist_pool, numpy_file_name)
std_Gls, results = ea.recursiveMakePow2Bins(freq_hist_pool, freq_hist_info, numpy_file_name)
# if freq_hist_info['runs'] * freq_hist_info['ranks'] > 1:
# std_Gls = np.mean(std_Gls,0)
# results = np.mean(results,0)
print("results.shape",results.shape)
if input_args.do_direct:
pickle_name = input_args.filestem + "nor_ordered.pickle"
if os.path.isfile(pickle_name) and not input_args.force_restrap:
......@@ -132,9 +159,9 @@ if __name__ == '__main__':
# now plot the recursive method versus the average of the variance for same size bins
for i_freq in range(10):
bins_at_freq = [ results[i,i_freq] for i in range(results.shape[0] - 1) ]
bins_at_freq = [ results[0,i,i_freq] for i in range(results.shape[1] - 1) ]
# mean_std_at_freq = [ direct_bin_results[i_freq,i,3] for i in range(direct_bin_results.shape[1]) ]
powers_of_2 = [i+1 for i in range(results.shape[0] - 1)]
powers_of_2 = [i+1 for i in range(results.shape[1] - 1)]
axs[0].plot(powers_of_2, bins_at_freq, label="$i\omega_0$[{}]".format(i_freq))
#axs[0,1].plot(powers_of_2, mean_std_at_freq, label="$i\omega_0$[{}]".format(i_freq))
......@@ -149,12 +176,10 @@ if __name__ == '__main__':
for i_freq in range(10):
# mean_var_at_freq = [ direct_bin_results[i_freq,i,2] for i in range(direct_bin_results.shape[1]) ]
bins_at_freq = [ results[i,i_freq] for i in range(results.shape[0] - 1) ]
powers_of_2 = [i+1 for i in range(results.shape[0] - 1)]
#axs[1,0].plot(powers_of_2, mean_var_at_freq, label="$i\omega_0$[{}]".format(i_freq))
#axs[1,0].plot(powers_of_2, mean_var_at_freq, label="$i\omega_0$[{}]".format(i_freq))
axs[1].plot(powers_of_2, bins_at_freq, label="$i\omega_0$[{}]".format(i_freq))
#axs[1,0].set_title("Mean variance of non recursive ordered bins")
axs[1].set_title("Reordered Samples Binning Analysis")
#axs[1,0].legend()
......
......@@ -13,6 +13,9 @@ import adios2
from mpi4py import MPI
import numpy as np
import re
import functools
import os.path
import pickle
comm = MPI.COMM_SELF
adios = adios2.ADIOS(comm)
......@@ -40,37 +43,69 @@ def fhFromFilename(file_name):
This keeps the adios_fh based API but removes need to import adios2 from clients"""
return adios2.open(file_name, "r", comm)
def readPickleOrFunc(pickle_name, func, *args, **kwargs):
"""The purpose of this function is to cache objects that are particularly expensive to create
if a file is present at pickle_name if will load it, otherwise it will run func with *args
and **kwargs and dump the results to the pickle file and return the result"""
if os.path.isfile(pickle_name):
try:
with open(pickle_name, 'rb') as the_pick:
return pickle.load(the_pick)
except:
raise Exception("pickle file is present but wrong some how, please delete it unless you know it to be valid")
else:
result = func(*args, **kwargs)
with open(pickle_name, 'wb') as the_pick:
pickle.dump(result, the_pick, pickle.HIGHEST_PROTOCOL)
return result
def tryPickleUnlessForce(force, pickle_name, func, *args, **kwargs):
"""this is convenience function replace frequent if bool then func(*args, **kwargs)...pickle.dump
else readPickleOrFunc(pickle_name, func, *args, **kwargs)"""
if force:
pickle_name = ''
return readPickleOrFunc(pickle_name, func, *args, **kwargs)
def readDCARun(adios_fh):
dca_parameters = {}
#in large files adios_hf.available_variables is terribly expensive.
available_variables = adios_fh.available_variables()
# make the dca run parameter names short and python friendly.
MCVars = [var for var in adios_fh.available_variables() if re.search(r"^/actual/", var)]
#MCVars = [var for var in available_variables if re.search(r"^/actual/", var)]
# the large number os substitutions was do to running a big run of data with an issue
# in the main code with the actual group getting shut.
for mcvar in MCVars:
mc_var_pretty = re.sub(r"^/","", mcvar)
mc_var_pretty = re.sub(r"[-/\+]",r"_", mc_var_pretty)
mc_var_pretty = re.sub(r"actual_domains_","", mc_var_pretty)
mc_var_pretty = re.sub(r"actual_functions_","", mc_var_pretty)
mc_var_pretty = re.sub(r"actual_Configurations","Configurations", mc_var_pretty)
mc_var_pretty = re.sub(r"actual_parameters_","", mc_var_pretty)
mc_var_pretty = re.sub(r"actual_phys","phys", mc_var_pretty)
mc_var_pretty = re.sub(r"actual_(DCA_loop_functions)",r"\1", mc_var_pretty)
dca_parameters[mc_var_pretty] = adios_fh.read(mcvar)
MCVars = [var for var in adios_fh.available_variables() if re.search(r"^/functions.*cluster", var)]
# for mcvar in MCVars:
# mc_var_pretty = re.sub(r"^/","", mcvar)
# mc_var_pretty = re.sub(r"[-/\+]",r"_", mc_var_pretty)
# mc_var_pretty = re.sub(r"actual_domains_","", mc_var_pretty)
# mc_var_pretty = re.sub(r"actual_functions_","", mc_var_pretty)
# mc_var_pretty = re.sub(r"actual_Configurations","Configurations", mc_var_pretty)
# mc_var_pretty = re.sub(r"actual_parameters_","", mc_var_pretty)
# mc_var_pretty = re.sub(r"actual_phys","phys", mc_var_pretty)
# mc_var_pretty = re.sub(r"actual_(DCA_loop_functions)",r"\1", mc_var_pretty)
# dca_parameters[mc_var_pretty] = adios_fh.read(mcvar)
#first lets drop the measurments as that's a huge number
non_meas_vars = [var for var in available_variables if "STQW_Configurations" not in var]
MCVars = [var for var in non_meas_vars if re.search(r"^/functions.*cluster", var)]
for mcvar in MCVars:
mc_var_pretty = re.sub(r"^/functions/","",mcvar)
dca_parameters[mc_var_pretty] = adios_fh.read(mcvar)
MCVars = [var for var in adios_fh.available_variables() if re.search(r"^/functions.*(?!cluster)", var)]
MCVars = [var for var in non_meas_vars if re.search(r"^/functions.*(?!cluster)", var)]
for mcvar in MCVars:
mc_var_pretty = re.sub(r"^/functions/","",mcvar)
dca_parameters[mc_var_pretty] = adios_fh.read(mcvar)
MCVars = [var for var in adios_fh.available_variables() if re.search("^/domains/", var)]
MCVars = [var for var in non_meas_vars if re.search("^/domains/", var)]
for mcvar in MCVars:
mc_var_pretty = re.sub(r"^/domains/","",mcvar)
mc_var_pretty = re.sub("-","_", mc_var_pretty)
mc_var_pretty = re.sub("/","_", mc_var_pretty)
dca_parameters[mc_var_pretty] = adios_fh.read(mcvar)
MCVars = [var for var in adios_fh.available_variables() if re.search("^/parameters/", var)]
MCVars = [var for var in non_meas_vars if re.search("^/parameters/", var)]
for mcvar in MCVars:
mc_var_pretty = re.sub(r"^/parameters/","",mcvar)
mc_var_pretty = re.sub("[-/\+]","_", mc_var_pretty)
......@@ -95,37 +130,67 @@ def readClusterG_k_w(adios_fh):
def readSample(adios_fh, rank, walker, samp):
readSampleConfigM_k_w(adios_fh, rank, walker, samp)
def readSampleConfigG_k_w(adios_fh, rank, walker, samp):
sample = "r_{}_meas_{}_w_{}".format(rank, samp, walker)
MCVars = [var for var in adios_fh.available_variables() if re.search("STQW_Configurations/{}/".format(sample), var)]
#print (rank, samp, walker, MCVars)
G_k_w = adios_fh.read(findVar(MCVars, "G_k_w"))
log_weight = adios_fh.read(findVar(MCVars, "log-weight"))[0]
aux_spins = np.array(adios_fh.read(findVar(MCVars, "spins")))
#print("shape", aux_spins.shape)
#print (aux_spins)
vertex_times = adios_fh.read(findVar(MCVars, "times"))
#sign = adios_fh.read(findVar(MCVars, "sign"))
return G_k_w, log_weight, aux_spins, vertex_times #, sign
def readSampleConfigG_k_w(adios_fh, rank, walker, samp_seq):
list_meas = []
for samp in samp_seq:
print('.', end='', flush=True)
#print (rank, samp, walker, MCVars)
G_k_w = readSampleG_k_w(adios_fh, rank, walker, samp)
log_weight = readSampleLogWeight(adios_fh, rank, walker, samp)
aux_spins = readSampleAuxSpins(adios_fh, rank, walker, samp)
vertex_times = readSampleVertexTimes(adios_fh, rank, walker, samp)
sign = readSampleSign(adios_fh, rank, walker, samp)
print("G_k_w.shape", G_k_w.shape)
list_meas.append([G_k_w, log_weight,aux_spins,vertex_times])
return list_meas
def readSampleConfigM_k_w(adios_fh, rank, walker, samp_seq):
list_meas = []
for samp in samp_seq:
print('.', end='', flush=True)
#print (rank, samp, walker, MCVars)
G_k_w = readSampleG_k_w(adios_fh, rank, walker, samp)
M_k_w = readSampleM_k_w(adios_fh, rank, walker, samp)
log_weight = readSampleLogWeight(adios_fh, rank, walker, samp)
aux_spins = readSampleAuxSpins(adios_fh, rank, walker, samp)
vertex_times = readSampleVertexTimes(adios_fh, rank, walker, samp)
sign = readSampleSign(adios_fh, rank, walker, samp)
list_meas.append([G_k_w, M_k_w, log_weight,aux_spins,vertex_times, sign])
return list_meas
def readSampleG_k_w(adios_fh, rank, walker, samp):
sample = "/STQW_Configurations/r_{}_meas_{}_w_{}/G_k_w".format(rank, samp, walker)
G_k_w = adios_fh.read(sample)
return G_k_w
def readSampleConfigM_k_w(adios_fh, rank, walker, samp):
sample = "r_{}_meas_{}_w_{}".format(rank, samp, walker)
MCVars = [var for var in adios_fh.available_variables() if re.search("STQW_Configurations/{}/".format(sample), var)]
print (rank, samp, walker, MCVars)
G_k_w = adios_fh.read(findVar(MCVars, "G_k_w"))
M_k_w = adios_fh.read(findVar(MCVars, "MFunction"))
log_weight = adios_fh.read(findVar(MCVars, "log-weight"))[0]
aux_spins = np.array(adios_fh.read(findVar(MCVars, "spins")))
print("shape", aux_spins.shape)
print (aux_spins)
vertex_times = adios_fh.read(findVar(MCVars, "times"))
sign = adios_fh.read(findVar(MCVars, "sign"))
return G_k_w, M_k_w, log_weight, aux_spins, vertex_times, sign
def readSampleM_k_w(adios_fh, rank, walker, samp):
sample = "/STQW_Configurations/r_{}_meas_{}_w_{}/MFunction".format(rank, samp, walker)
M_k_w = adios_fh.read(sample)
return M_k_w
def readVarDirect(adios_fh, rank, walker, samp, var_tag):
sample = "/STQW_Configurations/r_{}_meas_{}_w_{}/{}".format(rank, samp, walker, var_tag)
return adios_fh.read(sample)
def makeVarFunc(var_tag):
def var_func(adios_fh, rank, walker, samp):
return readVarDirect(adios_fh,rank,walker,samp, var_tag)
return var_func
readSampleLogWeight = makeVarFunc("log-weight")
# def readSampleLogWeight(fadios_fh, rank, walker, samp):
# return readVarDirect(fadios_fh, rank, walker, samp,"log-weight")
def readSampleAuxSpins(adios_fh, rank, walker, samp):
return readVarDirect(adios_fh, rank, walker, samp,"spins")
def readSampleVertexTimes(adios_fh, rank, walker, samp):
return readVarDirect(adios_fh, rank, walker, samp,"times")
def readSampleSign(adios_fh, rank, walker, samp):
return readVarDirect(adios_fh, rank, walker, samp,"sign")
def readG4Cluster(adios_fh):
# hangs
......
......@@ -12,9 +12,19 @@
import pyDCA.Adios2.read as pAr
import numpy as np
import re
from enum import Enum
from random import choices
from multiprocessing import Pool
class MLFunctions(Enum):
GKW = 1
MKW = 2
ml_reader_funcs = {
MLFunctions.GKW : pAr.readSampleG_k_w,
MLFunctions.MKW : pAr.readSampleM_k_w
}
def accumulating_read_and_std_spin_avg(file_names, fq_start=0, fq_end=-1):
"""Over a vector of DCA output filename produce per run std of measurements and std over all measurements.
You must run the same number of measurements in each run."""
......@@ -57,12 +67,13 @@ def accumulating_read_and_std_spin_avg(file_names, fq_start=0, fq_end=-1):
class SampleReader:
"""Functor class to parallelize reading many Adios2 bps faster when combining samples"""
def __init__(self, fq_start, fq_end, ranks, samples_per_rank):
def __init__(self, fq_start, fq_end, ranks, samples_per_rank, ml_function=MLFunctions.GKW):
"""These are the captures for the concurrent work"""
self.fq_start = fq_start
self.fq_end = fq_end
self.ranks = ranks
self.samples_per_rank = samples_per_rank
self.read_function = ml_reader_funcs[ml_function]
def readRunSamples(self, file_name):
"""This is the parallel work"""
......@@ -73,10 +84,11 @@ class SampleReader:
frequency_histagram_data_up = np.zeros([num_freq, total_samples])
frequency_histagram_data_dn = np.zeros([num_freq, total_samples])
run_fh = pAr.fhFromFilename(file_name)
# the order of these loops is important because we need the measuments for each rank as a contiguous block
for rank in range(self.ranks):
#print (rank, self.samples_per_rank[rank])
for samp in range(self.samples_per_rank[rank]):
G_k_w = pAr.readSampleG_k_w(run_fh,rank,0,samp+1)
G_k_w = self.read_function(run_fh,rank,0,samp+1)
frequency_histagram_data_up[:, meas_index] = np.real(G_k_w[self.fq_start:self.fq_end,0,0,0,0,0])
frequency_histagram_data_dn[:, meas_index] = np.real(G_k_w[self.fq_start:self.fq_end,0,1,0,1,0])
meas_index += 1
......@@ -146,7 +158,7 @@ def mergeSamplesNoSpinAvg(file_names, fq_start=0, fq_end=-1, dca_globals = {}, r
return frequency_histagram_running
def mergeSamplesSpinAvg(file_names, fq_start=0, fq_end=-1, dca_globals = {}, ranks=-1, samples_per_rank=[], procs=8 ):
def mergeSamplesSpinAvg(file_names, fq_start=0, fq_end=-1, dca_globals = {}, ranks=-1, samples_per_rank=[], procs=8, ml_function=MLFunctions.GKW ):
"""Merge samples from file_names into single [frequency, samples] array"""
if len(dca_globals) == 0:
dca_globals = pAr.withFhFromFilname(file_names[0],pAr.readDCARun)
......@@ -163,7 +175,7 @@ def mergeSamplesSpinAvg(file_names, fq_start=0, fq_end=-1, dca_globals = {}, ran
num_files = len(file_names)
total_samples = num_files * samples
sample_reader = SampleReader(fq_start, fq_end, ranks,samples_per_rank)
sample_reader = SampleReader(fq_start, fq_end, ranks,samples_per_rank, ml_function=ml_function)
#Pool.map returns ordered list of results, we could probably us an unordered dispatch here
with Pool(procs) as pool:
frequency_histagram_running_spin_avg_list = pool.map(sample_reader.readRunSamples, file_names)
......@@ -175,7 +187,11 @@ def mergeSamplesSpinAvg(file_names, fq_start=0, fq_end=-1, dca_globals = {}, ran
frequency_histagram_running_spin_avg[:,meas_index:last_index] = run_samp
meas_index = last_index
return frequency_histagram_running_spin_avg
freq_hist_info = { "ranks" : ranks,
"samples_per_rank" : samples_per_rank,
"runs" : num_files}
return frequency_histagram_running_spin_avg, freq_hist_info
class Resampler:
......@@ -435,7 +451,7 @@ def makeDirectBins(freq_hist_data, resamples, sample_size, sample_managers=None,
bin_result[:,3] = np.mean(bin_results[:,1,:],1)
return bin_result, bin_results
def recursiveMakePow2Bins(freq_hist_data, file_name = None):
def recursiveMakePow2Bins(freq_hist_data, freq_hist_info, file_name = None):
num_freq = freq_hist_data.shape[1]
total_meas = freq_hist_data.shape[0]
......@@ -444,41 +460,51 @@ def recursiveMakePow2Bins(freq_hist_data, file_name = None):
raise ValueError("you must have a numpy archive of the freq_hist_data")
data = np.load(file_name)
N = data.shape[0]
lmax = int(np.log2(N))
# this prevents an issue that results in a dropped bin when the Nbin's previous are odd.
N = 2**lmax
Gils = [data]
for l in range(1,lmax+1):
# print("l=",l)
Nbins = N//(2**l)
GG = np.zeros([Nbins, num_freq])
for i in range(Nbins):
# numpy here broadcasts and does this over all freq
GG[i,:] = 0.5*(Gils[l-1][2*i,:] + Gils[l-1][2*i+1,:])
Gils.append(GG)
# Bins for each order. Skipping order 0
Gil_lens = [Gils[l+1].shape[0] for l in range(lmax)]
print("Gil lens", Gil_lens)
delta_Gl = np.zeros([lmax, num_freq])
std_Gl = np.zeros([lmax, num_freq])
for l in range(lmax):
Nl = Gils[l+1].shape[0]
Gil = Gils[l+1]
Gmean = np.mean(Gil,0)
diff_freqs = np.zeros([num_freq])
for frq in range(num_freq):
diff_freqs[frq] = np.sum((Gil[:,frq] - Gmean[frq])**2)
# I find the Troyer argument for moving Nl outside to be dodgy
DGl = 1/Nl * np.sqrt(diff_freqs)
std_Gl[l,:] = np.std(Gil,0)
delta_Gl[l,:] = DGl
print ("For l:", l, DGl[0], std_Gl[l,0])
return std_Gl, delta_Gl
num_walkers = freq_hist_info['runs'] * freq_hist_info['ranks']
samples_per_rank_per_run = freq_hist_info['samples_per_rank'][0]
delta_Gls = []
std_Gls = []
for iw in range(num_walkers):
print(iw*samples_per_rank_per_run,iw*samples_per_rank_per_run+samples_per_rank_per_run)
walker_data = data[iw*samples_per_rank_per_run:iw*samples_per_rank_per_run+samples_per_rank_per_run, :]
N = walker_data.shape[0]
lmax = int(np.log2(N))
# this prevents an issue that results in a dropped bin when the Nbin's previous are odd.
N = 2**lmax
Gils = [walker_data]
for l in range(1,lmax+1):
# print("l=",l)
Nbins = N//(2**l)
GG = np.zeros([Nbins, num_freq])
for i in range(Nbins):
# numpy here broadcasts and does this over all freq
GG[i,:] = 0.5*(Gils[l-1][2*i,:] + Gils[l-1][2*i+1,:])
Gils.append(GG)
# Bins for each order. Skipping order 0
Gil_lens = [Gils[l+1].shape[0] for l in range(lmax)]
print("Gil lens", Gil_lens)
delta_Gl = np.zeros([lmax, num_freq])
std_Gl = np.zeros([lmax, num_freq])
for l in range(lmax):
Nl = Gils[l+1].shape[0]
Gil = Gils[l+1]
Gmean = np.mean(Gil,0)
diff_freqs = np.zeros([num_freq])
for frq in range(num_freq):
diff_freqs[frq] = np.sum((Gil[:,frq] - Gmean[frq])**2)
# I find the Troyer argument for moving Nl outside to be dodgy
DGl = 1/Nl * np.sqrt(diff_freqs)
std_Gl[l,:] = np.std(Gil,0)
delta_Gl[l,:] = DGl
print ("For l:", l, DGl[0], std_Gl[l,0])
std_Gls.append(std_Gl)
delta_Gls.append(delta_Gl)
return np.stack(std_Gls,axis=0), np.stack(delta_Gls,axis=0)
def recursiveMakePow2BinsResample(freq_hist_data, file_name = None):
# total_samples = freq_hist_data.shape[1]
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment