From 5a58e140b8d1569d237d9127e8afbba7c5cf29d5 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Mon, 25 Feb 2013 18:51:14 -0500 Subject: [PATCH 01/12] added new function to extract data and made script compatible with revision 7 temporarily --- src/utils/statepoint.py | 83 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 79 insertions(+), 4 deletions(-) diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 1227d89cbe..0a4653aceb 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -195,10 +195,11 @@ class StatePoint(BinaryFile): self.n_inactive, self.gen_per_batch = self._get_int(2) self.k_batch = self._get_double(self.current_batch) self.entropy = self._get_double(self.current_batch*self.gen_per_batch) - self.k_col_abs = self._get_double()[0] - self.k_col_tra = self._get_double()[0] - self.k_abs_tra = self._get_double()[0] - self.k_combined = self._get_double(2) + if self.revision >= 8: + self.k_col_abs = self._get_double()[0] + self.k_col_tra = self._get_double()[0] + self.k_abs_tra = self._get_double()[0] + self.k_combined = self._get_double(2) # Read number of meshes n_meshes = self._get_int()[0] @@ -423,3 +424,77 @@ class StatePoint(BinaryFile): # sum of squares, or it could be mean and stdev if self.generate_stdev() # has been called already. return t.results[filter_index, score_index] + + def extract_results(self,tally_id,score_str): + + # get tally + try: + tally = self.tallies[tally_id-1] + except: + print 'Tally does not exist' + return + + # get the score index if it is present + try: + idx = tally.scores.index(score_str) + except ValueError: + print 'Score does not exist' + return + + # create numpy array for mean and 95% CI + n_bins = len(tally.results) + n_filters = len(tally.filters) + n_scores = len(tally.scores) + meanv = np.zeros(n_bins) + unctv = np.zeros(n_bins) + filters = np.zeros((n_bins,n_filters)) + filtmax = np.zeros(n_filters+1) + meshmax = np.zeros(4) + filtmax[0] = 1 + meshmax[0] = 1 + + # get number of realizations + n = tally.n_realizations + + # get t-value + t_value = scipy.stats.t.ppf(0.975, n - 1) + + # calculate mean + meanv = tally.results[:,idx,0] + meanv = meanv / n + + # calculate 95% two-sided CI + unctv = tally.results[:,idx,1] + unctv = t_value*np.sqrt((unctv/n - meanv*meanv)/(n-1))/meanv + + # create output dictionary + data = {'mean':meanv,'CI95':unctv} + + # get bounds of filter bins + for akey in tally.filters.keys(): + idx = tally.filters.keys().index(akey) + filtmax[n_filters - idx] = tally.filters[akey].length + + # compute bin info + i = 0 + while i < n_filters: + + # compute indices for filter combination + filters[:,n_filters - i - 1] = np.floor((np.arange(n_bins) % np.prod(filtmax[0:i+2]))/(np.prod(filtmax[0:i+1]))) + 1 + + # append in dictionary bin with filter + data.update({tally.filters.keys()[n_filters - i - 1]:filters[:,n_filters - i - 1]}) + + # check for mesh + if tally.filters.keys()[n_filters - i - 1] == 'mesh': + mesh_idx = tally.filters['mesh'].bins + meshmax[1:4] = self.meshes[tally.filters['mesh'].bins[0] - 1].dimension + mesh_bins = np.zeros((n_bins,3)) + mesh_bins[:,0] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:2]))/(np.prod(meshmax[0:1]))) + 1 + mesh_bins[:,1] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:3]))/(np.prod(meshmax[0:2]))) + 1 + mesh_bins[:,2] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:4]))/(np.prod(meshmax[0:3]))) + 1 + data.update({'mesh':zip(mesh_bins[:,0],mesh_bins[:,1],mesh_bins[:,2])}) + i += 1 + + return data + From 39a46f930412c3ec3ecdac08f734545810b39adc Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 26 Feb 2013 10:09:15 -0500 Subject: [PATCH 02/12] bin info is now saved in data dictionary when extracting results --- src/utils/statepoint.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 0a4653aceb..f90b626257 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -496,5 +496,16 @@ class StatePoint(BinaryFile): data.update({'mesh':zip(mesh_bins[:,0],mesh_bins[:,1],mesh_bins[:,2])}) i += 1 - return data + # add in maximum bin filters and order + b = tally.filters.keys() + b.reverse() + filtmax = list(filtmax[1:]) + try: + idx = b.index('mesh') + filtmax[idx] = np.max(mesh_bins[:,2]) + filtmax.insert(idx,np.max(mesh_bins[:,1])) + filtmax.insert(idx,np.max(mesh_bins[:,0])) + except ValueError: pass + data.update({'bin_order':b,'bin_max':filtmax}) + return data From 5406a9cb37089e00384b6b0b22bde01b7df40e1b Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 26 Feb 2013 10:09:45 -0500 Subject: [PATCH 03/12] reads in statepoint and process source distribution for gnuplot --- src/utils/eigenfunction_rms.py | 292 +++++++++++---------------------- 1 file changed, 97 insertions(+), 195 deletions(-) diff --git a/src/utils/eigenfunction_rms.py b/src/utils/eigenfunction_rms.py index 9777a65ec0..16888a7788 100644 --- a/src/utils/eigenfunction_rms.py +++ b/src/utils/eigenfunction_rms.py @@ -1,206 +1,108 @@ #!/usr/bin/python # Filename: eigenfunction_rms.py -import h5py +# import packages +import statepoint import numpy as np -import cPickle -import matplotlib.pyplot as plt -import sys -# -class EigenFunction: -# - '''Represents the reference eigenfunction''' - reference = 0. -# - def __init__(self,data): -# - '''Initializes the eigenfunction''' - self.function = 0. - self.data = data - self.meanfunction = 0. -# - def read_hdf5(self,h5_file,cycle): -# - '''Read data from HDF5 file''' - self.cycle = cycle - f = h5py.File(h5_file,'r') - group = '/cycle'+str(cycle)+'/'+self.data - dataset = f[group] - self.function = np.empty(dataset.shape,dataset.dtype) - dataset.read_direct(self.function) - self.function = self.function - self.iamref = 'F' -# - def set_reference(self): -# - '''Sets instance to be reference calc''' - self.iamref = 'T' - EigenFunction.reference = self.function -# - def compute_rms(self): -# - '''Computes RMS value''' - Np = self.function.size - Np = 41772 - tmp = (self.meanfunction - EigenFunction.reference)**2 - tmp2 = tmp.sum() - self.rms = np.sqrt((1.0/float(Np))*tmp2) -# -def read_runs(runpath,hdfile,cycle_start,cycle_end,run_start,run_end,data): - runlist = [] - tmp = EigenFunction(data) - tmp.read_hdf5(runpath+str(run_start)+'/'+hdfile,cycle_start) - indices = tmp.function.shape # extent of all dimensions - i = cycle_start - while i <= cycle_end: - meantmp = EigenFunction(data) # init mean object - j = run_start - runs = np.zeros((run_end - run_start + 1,indices[0],indices[1],indices[2],indices[3])) # init runs array - while j <= run_end: - tmp.read_hdf5(runpath+str(j)+'/'+hdfile,i) # read hdf5 file - runs[j-run_start] = tmp.function # put function into runs - if i == cycle_start: - print runs[j-run_start,0,150,150,0] - j += 1 - meantmp.meanfunction = np.average(runs, axis=0) # compute the mean - meantmp.function = runs - runlist.append(meantmp) - print 'Read in from path: '+runpath+' Cycle: '+str(i) - i += 1 - return runlist +import os -def create_reference(runpath,hdfile,cycle,run_start,run_end,data): +def main(): -# calculate reference solution - print 'Calculating Reference solution...' - tmp = EigenFunction(data) - tmp.read_hdf5(runpath+str(run_start)+'/'+hdfile,cycle) # load first eigenfunction - print 'Read in: '+runpath+str(1)+hdfile+' '+str(tmp.function[0,150,150,0]) - indices = tmp.function.shape # extent of all dimensions - ref = np.zeros((run_end,indices[0],indices[1],indices[2],indices[3])) # initialize ref array - ref[0] = tmp.function # set the first run in ref - i = run_start + 1 - while i <= run_end: # begin loop around all runs - tmp.read_hdf5(runpath+str(i)+'/'+hdfile,cycle) - print 'Read in: '+runpath+str(i)+hdfile+' '+str(tmp.function[0,150,150,0]) - ref[i - run_start] = tmp.function - i += 1 - meanref = np.average(ref, axis=0) # compute average of all runs - EigenFunction.reference = meanref # set to global space in EigenFunction instances - return meanref + # read in statepoint header data + sp = statepoint.StatePoint('statepoint.ref.binary') -def compute_rms(runlist): + # read in results + sp.read_results() - # calculate rms array - print 'Calculating rms...' - rms = np.zeros(len(runlist)) - i = 0 - while i < len(runlist): - runlist[i].compute_rms() - rms[i] = runlist[i].rms - i += 1 - return rms + # extract results + results = sp.extract_results(2,'nu-fission') -def plot_rms(rms): - print 'Generating plot...' - ax = plt.subplot(111) - size = rms.shape[0] - x = np.linspace(1,size,size)*1e6 - y = (rms[0]/x[0]**(-0.5))*x**(-0.5) - plt.loglog(x,rms*100,'b+') - plt.loglog(x,y*100,'g--') - ax.xaxis.grid(True,'minor') - ax.yaxis.grid(True,'minor') - ax.xaxis.grid(True,'major',linewidth=2) - ax.yaxis.grid(True,'major',linewidth=2) - plt.xlabel('# of Total Neutron Histories (active cycles)') - plt.ylabel('RMS Error [%]') - plt.legend(('1 million (10 runs)','Ideal Error 1 mil')) - return + # extract means and copy + mean = results['mean'].copy() + + # reshape and integrate over energy + mean = mean.reshape(results['bin_max'],order='F') + mean = np.sum(mean,0) + mean = np.sum(mean,0) + mean = mean/mean.sum()*(mean > 1.e-8).sum() + + # write gnuplot file + write_gnuplot('testsrc_pin','Pin Mesh',mean,np.size(mean,0)) + + # extract results + results = sp.extract_results(5,'nu-fission') + + # extract means and copy + mean = results['mean'].copy() + + # reshape and integrate over energy + mean = mean.reshape(results['bin_max'],order='F') + mean = np.sum(mean,0) + mean = np.sum(mean,0) + mean = mean/mean.sum()*(mean > 1.e-8).sum() + + # write gnuplot file + write_gnuplot('testsrc_quart','Quarter Assembly Mesh',mean,np.size(mean,0)) + + # extract results + results = sp.extract_results(8,'nu-fission') + + # extract means and copy + mean = results['mean'].copy() + + # reshape and integrate over energy + mean = mean.reshape(results['bin_max'],order='F') + mean = np.sum(mean,0) + mean = np.sum(mean,0) + mean = mean/mean.sum()*(mean > 1.e-8).sum() + + # write gnuplot file + write_gnuplot('testsrc_assy','Assembly Mesh',mean,np.size(mean,0)) + +def write_gnuplot(path,name,src,size): + + # Header String for GNUPLOT + headerstr = """#!/usr/bin/env gnuplot + +set terminal pdf enhanced +set output '{output}' +set palette defined (0 '#000090', 1 '#000fff', 2 '#0090ff', 3 '#0fffee', 4 '#90ff70', 5 '#ffee00', 6 '#ff7000', 7 '#ee0000', 8 '#7f0000') +set view map +set size ratio -1 +set lmargin at screen 0.10 +set rmargin at screen 0.90 +set bmargin at screen 0.15 +set tmargin at screen 0.90 +unset xtics +unset ytics +set title '{title}'""".format(output=path+'.pdf',title=name) + + # Write out the plot string + pltstr = "splot '-' matrix with image " + + # Write out the data string + i = 0 + datastr = '' + while i < size: + j = 0 + while j < size: + datastr = datastr + '{0} '.format(src[i,j][0]) + j += 1 + datastr = datastr + '\n' + i += 1 + + # replace all nan with zero + datastr = datastr.replace('nan','0.0') + + # Concatenate all + outstr = headerstr + '\n' + pltstr + '\n' + datastr + + # Write File + with open(path+".plot",'w') as f: + f.write(outstr) + + # Run GNUPLOT + os.system('gnuplot ' + path+".plot") -def plot_source(source): - plt.figure() - X = np.linspace(1,272,272) - Y = np.linspace(1,272,272) - Y,X = np.meshgrid(Y,X) - plt.contourf(X,Y,source[0,:,:,0],100) - plt.colorbar() - plt.xlabel('Mesh Cell in x-direction') - plt.ylabel('Mesh Cell in y-direction') - plt.title('OPR Converged Fission Source Distribution') - plt.show() -# if __name__ == "__main__": - - if sys.argv[1] == 'restart': - - # load in data - print 'Loading input...' - filein = open('rms.out','r') - output = cPickle.load(filein) - filein.close() - rms = output['1milrms'] - meanref = output['ref'] - EigenFunction.reference = meanref - - # plot rms - plot_rms(rms) - - # plot mean source distribution - plot_source(meanref) - - elif sys.argv[1] == 'interactive': - - # load in data - print 'Loading input...' - filein = open('rms.out','r') - output = cPickle.load(filein) - filein.close() - rms = output['1milrms'] - meanref = output['ref'] - EigenFunction.reference = meanref - - # pop an interactive python shell - from IPython import embed - embed() - - else: - - # calculate reference solution - runpath = '/media/Backup/opr_runs/64mil/run' - hdfile = 'output.h5' - cycle = 210 - run_start = 1 - run_end = 4 - data = 'openmc_src' - meanref = create_reference(runpath,hdfile,cycle,run_start,run_end,data) - - # calculate rms for 1 million case - runpath = '/media/Backup/opr_runs/1mil/run' - hdfile = 'output.h5' - cycle_start = 201 - cycle_end = 840 - run_start = 1 - run_end = 10 - data = 'openmc_src' - onemil = read_runs(runpath,hdfile,cycle_start,cycle_end,run_start,run_end,data) - - # calculate rms array - rms = compute_rms(onemil) - - # write out numpy array to binary file - print 'Writing output...' - output = {} - output.update({'1milrms':rms}) - output.update({'ref':meanref}) - fileout = open('rms.out','wb') - cPickle.dump(output,fileout) - fileout.close() - - # plot rms - plot_rms(rms) - - # plot mean source distribution - plot_source(meanref) - - + main() From 1d24a867fa8a78c72bc249af501f5dcac12d2400 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 26 Feb 2013 11:33:44 -0500 Subject: [PATCH 04/12] L-2 norms between reference and statepoints are now calculated --- src/utils/eigenfunction_rms.py | 109 ++++++++++++++++++++++++--------- 1 file changed, 81 insertions(+), 28 deletions(-) diff --git a/src/utils/eigenfunction_rms.py b/src/utils/eigenfunction_rms.py index 16888a7788..2487db3e06 100644 --- a/src/utils/eigenfunction_rms.py +++ b/src/utils/eigenfunction_rms.py @@ -5,8 +5,9 @@ import statepoint import numpy as np import os +import sys -def main(): +def main(batch_start, batch_end): # read in statepoint header data sp = statepoint.StatePoint('statepoint.ref.binary') @@ -14,8 +15,53 @@ def main(): # read in results sp.read_results() + # extract reference mean + mean_ref = extract_mean(sp, 2,'nu-fission') + + # write gnuplot file + write_src_gnuplot('testsrc_pin','Pin Mesh',mean_ref,np.size(mean_ref,0)) + + # preallocate arrays + hists = np.zeros(batch_end - batch_start + 1) + norms = np.zeros(batch_end - batch_start + 1) + + i = batch_start + while i <= batch_end: + + # process statepoint + sp = statepoint.StatePoint('statepoint.'+str(i)+'.binary') + sp.read_results() + + # extract mean + mean = extract_mean(sp, 2, 'nu-fission') + + # calculate L2 norm + norm = np.linalg.norm(mean - mean_ref) + + # get history information + n_inactive = sp.n_inactive + current_batch = sp.current_batch + n_particles = sp.n_particles + gen_per_batch = sp.gen_per_batch + n_histories = (current_batch - n_inactive)*n_particles*gen_per_batch + + # batch in vectors + hists[i - n_inactive - 1] = n_histories + norms[i - n_inactive - 1] = norm + + # print + print 'Batch: '+str(i)+' Histories: '+str(n_histories)+' Norm: '+str(norm) + + i += 1 + + # write out gnuplot file + write_norm_gnuplot('norms','NORMS',hists,norms,np.size(hists)) + + +def extract_mean(sp, tally_id,score_id): + # extract results - results = sp.extract_results(2,'nu-fission') + results = sp.extract_results(tally_id,score_id) # extract means and copy mean = results['mean'].copy() @@ -26,40 +72,45 @@ def main(): mean = np.sum(mean,0) mean = mean/mean.sum()*(mean > 1.e-8).sum() - # write gnuplot file - write_gnuplot('testsrc_pin','Pin Mesh',mean,np.size(mean,0)) + return mean - # extract results - results = sp.extract_results(5,'nu-fission') +def write_norm_gnuplot(path,name,xdat,ydat,size): - # extract means and copy - mean = results['mean'].copy() + # Header String for GNUPLOT + headerstr = """#!/usr/bin/env gnuplot - # reshape and integrate over energy - mean = mean.reshape(results['bin_max'],order='F') - mean = np.sum(mean,0) - mean = np.sum(mean,0) - mean = mean/mean.sum()*(mean > 1.e-8).sum() +set terminal pdf enhanced +set output '{output}' +set ylabel 'L-2 norm' +set xlabel 'Histories' +set log x +set log y +set title '{title}'""".format(output=path+'.pdf',title=name) - # write gnuplot file - write_gnuplot('testsrc_quart','Quarter Assembly Mesh',mean,np.size(mean,0)) + # Write out the plot string + pltstr = "plot '-' using 1:2 with lines" - # extract results - results = sp.extract_results(8,'nu-fission') + # Write out the data string + i = 0 + datastr = '' + while i < size: + datastr = datastr + '{0} {1}\n'.format(xdat[i],ydat[i]) + i += 1 - # extract means and copy - mean = results['mean'].copy() + # replace all nan with zero +# datastr = datastr.replace('nan','0.0') - # reshape and integrate over energy - mean = mean.reshape(results['bin_max'],order='F') - mean = np.sum(mean,0) - mean = np.sum(mean,0) - mean = mean/mean.sum()*(mean > 1.e-8).sum() + # Concatenate all + outstr = headerstr + '\n' + pltstr + '\n' + datastr - # write gnuplot file - write_gnuplot('testsrc_assy','Assembly Mesh',mean,np.size(mean,0)) + # Write File + with open(path+".plot",'w') as f: + f.write(outstr) -def write_gnuplot(path,name,src,size): + # Run GNUPLOT +# os.system('gnuplot ' + path+".plot") + +def write_src_gnuplot(path,name,src,size): # Header String for GNUPLOT headerstr = """#!/usr/bin/env gnuplot @@ -105,4 +156,6 @@ set title '{title}'""".format(output=path+'.pdf',title=name) os.system('gnuplot ' + path+".plot") if __name__ == "__main__": - main() + batch_start = int(sys.argv[1]) + batch_end = int(sys.argv[2]) + main(batch_start, batch_end) From d02d1bdce565dc7141fa1de87021008544ff35af Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 6 Mar 2013 11:06:06 -0500 Subject: [PATCH 05/12] removed commented lines of code --- src/utils/eigenfunction_rms.py | 35 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/src/utils/eigenfunction_rms.py b/src/utils/eigenfunction_rms.py index 2487db3e06..5f87389da6 100644 --- a/src/utils/eigenfunction_rms.py +++ b/src/utils/eigenfunction_rms.py @@ -7,7 +7,7 @@ import numpy as np import os import sys -def main(batch_start, batch_end): +def main(tally_id, score_id, batch_start, batch_end, name): # read in statepoint header data sp = statepoint.StatePoint('statepoint.ref.binary') @@ -16,10 +16,10 @@ def main(batch_start, batch_end): sp.read_results() # extract reference mean - mean_ref = extract_mean(sp, 2,'nu-fission') + mean_ref = extract_mean(sp, tally_id, score_id) # write gnuplot file - write_src_gnuplot('testsrc_pin','Pin Mesh',mean_ref,np.size(mean_ref,0)) + write_src_gnuplot('testsrc_pin','Pin mesh',mean_ref,np.size(mean_ref,0)) # preallocate arrays hists = np.zeros(batch_end - batch_start + 1) @@ -33,7 +33,7 @@ def main(batch_start, batch_end): sp.read_results() # extract mean - mean = extract_mean(sp, 2, 'nu-fission') + mean = extract_mean(sp, tally_id, score_id) # calculate L2 norm norm = np.linalg.norm(mean - mean_ref) @@ -46,8 +46,8 @@ def main(batch_start, batch_end): n_histories = (current_batch - n_inactive)*n_particles*gen_per_batch # batch in vectors - hists[i - n_inactive - 1] = n_histories - norms[i - n_inactive - 1] = norm + hists[i - batch_start] = n_histories + norms[i - batch_start] = norm # print print 'Batch: '+str(i)+' Histories: '+str(n_histories)+' Norm: '+str(norm) @@ -55,8 +55,7 @@ def main(batch_start, batch_end): i += 1 # write out gnuplot file - write_norm_gnuplot('norms','NORMS',hists,norms,np.size(hists)) - + write_norm_gnuplot(name,hists,norms,np.size(hists)) def extract_mean(sp, tally_id,score_id): @@ -74,7 +73,7 @@ def extract_mean(sp, tally_id,score_id): return mean -def write_norm_gnuplot(path,name,xdat,ydat,size): +def write_norm_gnuplot(path,xdat,ydat,size): # Header String for GNUPLOT headerstr = """#!/usr/bin/env gnuplot @@ -85,7 +84,7 @@ set ylabel 'L-2 norm' set xlabel 'Histories' set log x set log y -set title '{title}'""".format(output=path+'.pdf',title=name) +""".format(output=path+'.pdf') # Write out the plot string pltstr = "plot '-' using 1:2 with lines" @@ -97,9 +96,6 @@ set title '{title}'""".format(output=path+'.pdf',title=name) datastr = datastr + '{0} {1}\n'.format(xdat[i],ydat[i]) i += 1 - # replace all nan with zero -# datastr = datastr.replace('nan','0.0') - # Concatenate all outstr = headerstr + '\n' + pltstr + '\n' + datastr @@ -108,7 +104,7 @@ set title '{title}'""".format(output=path+'.pdf',title=name) f.write(outstr) # Run GNUPLOT -# os.system('gnuplot ' + path+".plot") + os.system('gnuplot ' + path+".plot") def write_src_gnuplot(path,name,src,size): @@ -126,7 +122,7 @@ set bmargin at screen 0.15 set tmargin at screen 0.90 unset xtics unset ytics -set title '{title}'""".format(output=path+'.pdf',title=name) +set title '{title}'""".format(output=path+'.pdf',title=name) # Write out the plot string pltstr = "splot '-' matrix with image " @@ -156,6 +152,9 @@ set title '{title}'""".format(output=path+'.pdf',title=name) os.system('gnuplot ' + path+".plot") if __name__ == "__main__": - batch_start = int(sys.argv[1]) - batch_end = int(sys.argv[2]) - main(batch_start, batch_end) + tally_id = int(sys.argv[1]) + score_id = sys.argv[2] + batch_start = int(sys.argv[3]) + batch_end = int(sys.argv[4]) + name = sys.argv[5] + main(tally_id, score_id, batch_start, batch_end, name) From e195163727ab585a22910ed0b53c828f66c6562d Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 26 Mar 2013 15:50:37 -0400 Subject: [PATCH 06/12] added current tally handling and reveresed mesh indices to z,y,x ordering --- src/utils/statepoint.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index f90b626257..11a2f4dfe6 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -439,6 +439,7 @@ class StatePoint(BinaryFile): idx = tally.scores.index(score_str) except ValueError: print 'Score does not exist' + print tally.scores return # create numpy array for mean and 95% CI @@ -488,11 +489,16 @@ class StatePoint(BinaryFile): # check for mesh if tally.filters.keys()[n_filters - i - 1] == 'mesh': mesh_idx = tally.filters['mesh'].bins - meshmax[1:4] = self.meshes[tally.filters['mesh'].bins[0] - 1].dimension + dims = self.meshes[tally.filters['mesh'].bins[0] - 1].dimension + dims.reverse() + dims = np.asarray(dims) + if score_str == 'current': + dims += 1 + meshmax[1:4] = dims mesh_bins = np.zeros((n_bins,3)) - mesh_bins[:,0] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:2]))/(np.prod(meshmax[0:1]))) + 1 + mesh_bins[:,2] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:2]))/(np.prod(meshmax[0:1]))) + 1 mesh_bins[:,1] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:3]))/(np.prod(meshmax[0:2]))) + 1 - mesh_bins[:,2] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:4]))/(np.prod(meshmax[0:3]))) + 1 + mesh_bins[:,0] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:4]))/(np.prod(meshmax[0:3]))) + 1 data.update({'mesh':zip(mesh_bins[:,0],mesh_bins[:,1],mesh_bins[:,2])}) i += 1 @@ -505,6 +511,7 @@ class StatePoint(BinaryFile): filtmax[idx] = np.max(mesh_bins[:,2]) filtmax.insert(idx,np.max(mesh_bins[:,1])) filtmax.insert(idx,np.max(mesh_bins[:,0])) + except ValueError: pass data.update({'bin_order':b,'bin_max':filtmax}) From 20ff4a5a6496014fd7df193672b5bd594866f473 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 2 Apr 2013 14:16:37 -0400 Subject: [PATCH 07/12] edited coding style --- src/utils/statepoint.py | 65 +++++++++++++++++++++++------------------ 1 file changed, 36 insertions(+), 29 deletions(-) diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 11a2f4dfe6..76ecf56645 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -429,10 +429,10 @@ class StatePoint(BinaryFile): # get tally try: - tally = self.tallies[tally_id-1] + tally = self.tallies[tally_id-1] except: - print 'Tally does not exist' - return + print 'Tally does not exist' + return # get the score index if it is present try: @@ -473,46 +473,53 @@ class StatePoint(BinaryFile): # get bounds of filter bins for akey in tally.filters.keys(): - idx = tally.filters.keys().index(akey) - filtmax[n_filters - idx] = tally.filters[akey].length + idx = tally.filters.keys().index(akey) + filtmax[n_filters - idx] = tally.filters[akey].length # compute bin info i = 0 while i < n_filters: - # compute indices for filter combination - filters[:,n_filters - i - 1] = np.floor((np.arange(n_bins) % np.prod(filtmax[0:i+2]))/(np.prod(filtmax[0:i+1]))) + 1 + # compute indices for filter combination + filters[:,n_filters - i - 1] = np.floor((np.arange(n_bins) % + np.prod(filtmax[0:i+2]))/(np.prod(filtmax[0:i+1]))) + 1 - # append in dictionary bin with filter - data.update({tally.filters.keys()[n_filters - i - 1]:filters[:,n_filters - i - 1]}) + # append in dictionary bin with filter + data.update({tally.filters.keys()[n_filters - i - 1]: + filters[:,n_filters - i - 1]}) - # check for mesh - if tally.filters.keys()[n_filters - i - 1] == 'mesh': - mesh_idx = tally.filters['mesh'].bins - dims = self.meshes[tally.filters['mesh'].bins[0] - 1].dimension - dims.reverse() - dims = np.asarray(dims) - if score_str == 'current': - dims += 1 - meshmax[1:4] = dims - mesh_bins = np.zeros((n_bins,3)) - mesh_bins[:,2] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:2]))/(np.prod(meshmax[0:1]))) + 1 - mesh_bins[:,1] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:3]))/(np.prod(meshmax[0:2]))) + 1 - mesh_bins[:,0] = np.floor(((filters[:,n_filters - i - 1] - 1) % np.prod(meshmax[0:4]))/(np.prod(meshmax[0:3]))) + 1 - data.update({'mesh':zip(mesh_bins[:,0],mesh_bins[:,1],mesh_bins[:,2])}) - i += 1 + # check for mesh + if tally.filters.keys()[n_filters - i - 1] == 'mesh': + mesh_idx = tally.filters['mesh'].bins + dims = self.meshes[tally.filters['mesh'].bins[0] - 1].dimension + dims.reverse() + dims = np.asarray(dims) + if score_str == 'current': + dims += 1 + meshmax[1:4] = dims + mesh_bins = np.zeros((n_bins,3)) + mesh_bins[:,2] = np.floor(((filters[:,n_filters - i - 1] - 1) % + np.prod(meshmax[0:2]))/(np.prod(meshmax[0:1]))) + 1 + mesh_bins[:,1] = np.floor(((filters[:,n_filters - i - 1] - 1) % + np.prod(meshmax[0:3]))/(np.prod(meshmax[0:2]))) + 1 + mesh_bins[:,0] = np.floor(((filters[:,n_filters - i - 1] - 1) % + np.prod(meshmax[0:4]))/(np.prod(meshmax[0:3]))) + 1 + data.update({'mesh':zip(mesh_bins[:,0],mesh_bins[:,1], + mesh_bins[:,2])}) + i += 1 # add in maximum bin filters and order b = tally.filters.keys() b.reverse() filtmax = list(filtmax[1:]) try: - idx = b.index('mesh') - filtmax[idx] = np.max(mesh_bins[:,2]) - filtmax.insert(idx,np.max(mesh_bins[:,1])) - filtmax.insert(idx,np.max(mesh_bins[:,0])) + idx = b.index('mesh') + filtmax[idx] = np.max(mesh_bins[:,2]) + filtmax.insert(idx,np.max(mesh_bins[:,1])) + filtmax.insert(idx,np.max(mesh_bins[:,0])) - except ValueError: pass + except ValueError: + pass data.update({'bin_order':b,'bin_max':filtmax}) return data From 4850a63a100f4f2c90a8dd5083a1a2b8b40d70ed Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 2 Apr 2013 14:25:41 -0400 Subject: [PATCH 08/12] added spacing between args --- src/utils/statepoint.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index c4e952aece..9ab8db445e 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -440,7 +440,7 @@ class StatePoint(object): # has been called already. return t.results[filter_index, score_index] - def extract_results(self,tally_id,score_str): + def extract_results(self, tally_id, score_str): # get tally try: From 4316904cb45a4a09b091bbb09264380243186f1c Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 2 Apr 2013 17:07:35 -0400 Subject: [PATCH 09/12] fixed indent error for statepoint revision --- src/utils/statepoint.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 9ab8db445e..fbac9ef246 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -184,10 +184,10 @@ class StatePoint(object): self.entropy = self._get_double( self.current_batch*self.gen_per_batch, path='entropy') if self.revision >= 8: - self.k_col_abs = self._get_double(path='k_col_abs')[0] - self.k_col_tra = self._get_double(path='k_col_tra')[0] - self.k_abs_tra = self._get_double(path='k_abs_tra')[0] - self.k_combined = self._get_double(2, path='k_combined') + self.k_col_abs = self._get_double(path='k_col_abs')[0] + self.k_col_tra = self._get_double(path='k_col_tra')[0] + self.k_abs_tra = self._get_double(path='k_abs_tra')[0] + self.k_combined = self._get_double(2, path='k_combined') # Read number of meshes n_meshes = self._get_int(path='tallies/n_meshes')[0] From 8737ce22dcbcf6587bd59cbac3652b5aad546dc4 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Wed, 3 Apr 2013 10:56:20 -0400 Subject: [PATCH 10/12] added list conversion when using HDF5 for mesh reading --- src/utils/statepoint.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index fbac9ef246..6732607ea8 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -505,8 +505,7 @@ class StatePoint(object): # check for mesh if tally.filters.keys()[n_filters - i - 1] == 'mesh': - mesh_idx = tally.filters['mesh'].bins - dims = self.meshes[tally.filters['mesh'].bins[0] - 1].dimension + dims = list(self.meshes[tally.filters['mesh'].bins[0] - 1].dimension) dims.reverse() dims = np.asarray(dims) if score_str == 'current': From 8a784331c1fcf81541d1ae515536255d12707cdc Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Mon, 8 Apr 2013 08:29:53 -0400 Subject: [PATCH 11/12] added docstring for extract_results method --- src/utils/statepoint.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index 6732607ea8..f8dd76ada6 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -441,6 +441,18 @@ class StatePoint(object): return t.results[filter_index, score_index] def extract_results(self, tally_id, score_str): + """Returns a tally results dictionary given a tally_id and score string. + + Parameters + ---------- + tally_id : int + Index for the tally in StatePoint.tallies list + + score_str : string + Corresponds to the string entered for a score in tallies.xml. + For a flux score extraction it would be 'score' + + """ # get tally try: From fec7c049dfaa3c99bd657c55dd12e57ca34f1bcd Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Mon, 8 Apr 2013 08:31:43 -0400 Subject: [PATCH 12/12] replaced while loop with for loop construction --- src/utils/statepoint.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py index f8dd76ada6..e91d474254 100644 --- a/src/utils/statepoint.py +++ b/src/utils/statepoint.py @@ -504,8 +504,7 @@ class StatePoint(object): filtmax[n_filters - idx] = tally.filters[akey].length # compute bin info - i = 0 - while i < n_filters: + for i in range(n_filters): # compute indices for filter combination filters[:,n_filters - i - 1] = np.floor((np.arange(n_bins) %