Skip to content

Commit

Permalink
Updating plot script to handle error bar colors with series colors an…
Browse files Browse the repository at this point in the history 
…d add hairline outlines to markers for enhanced visibility.
  • Loading branch information
@sskutnik
sskutnik committed Jul 17, 2016
1 parent 9f4d5e1 commit 2bebfaf
Showing 1 changed file with 160 additions and 49 deletions.
209 changes: 160 additions & 49 deletions plot_UNF_benchmarks.py
View file
Original file line number Diff line number Diff line change
@@ -1,11 +1,13 @@
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
#import numpy as np
import numpy as np
import seaborn as sns
import re
import itertools

sns.set(style="ticks",font_scale=1.5)
sns.set_context(rc={'lines.markeredgewidth': 1.0, 'lines.markeredgecolor': 'k', 'lines.linecolor': 'k'})

IMAGE_DEST='./images/'
TEX_DEST='./latex/'
Expand Down Expand Up @@ -91,7 +93,9 @@ def plotIsos(df_iso, labels, seriesInfo, pltTitle, filename, sigmaLabel="Sigma",
plt.savefig(IMAGE_DEST + filename + '.' + OUTPUT_FORMAT)
plt.close()

def plotIsosFP(df_iso, assemName, seriesInfo, filename, burnups={},plotAspect=1.05):
def plotIsosFP(df_iso, assemName, seriesInfo, filename, burnups={},plotAspect=1.05,colWrap=GLOB_COL_WRAP, fontScale=1.5):

sns.set(style="ticks",font_scale=fontScale)

isoNumbers = []
idx = 0
Expand All @@ -100,24 +104,34 @@ def plotIsosFP(df_iso, assemName, seriesInfo, filename, burnups={},plotAspect=1.
idx = idx + 1
isoDict = dict(isoNumbers)

isoMarkers = ['o', 'o', 's', '^', 'v']
isoMarkers = ['o', 'o', 's', '^', 'v', '*', 'H']
edgeColors = ['k'] * len(df_iso['Evaluation'].unique())
edgeWidths = ['1'] * len(df_iso['Evaluation'].unique())
df_iso["IsoNum"] = df_iso["Isotope"].map(isoDict)

g = sns.factorplot(data=df_iso,kind='point',col='Sample',x='Isotope',y='C/E',hue='Evaluation',legend=False,scale=1.0,size=4,
col_wrap=GLOB_COL_WRAP,linestyles='None',aspect=plotAspect, markers=isoMarkers)
g.map(plt.errorbar,"IsoNum","C/E",'Sigma',data=df_iso,ls='',ms=0,zorder=0,elinewidth=2,color="gray")
palette = itertools.cycle(sns.color_palette())
ebDict = { }

for eval in df_iso['Evaluation'].unique():
ebDict[eval] = palette.next()

g = sns.factorplot(data=df_iso,kind='point',col='Sample',x='Isotope',y='C/E',hue='Evaluation',legend=False,scale=1.75,size=4,
col_wrap=colWrap,linestyles='None',aspect=plotAspect,markers=isoMarkers)

for ax in g.axes:
for mk in ax.collections:
mk.set_edgecolor('black')
mk.set_linewidth(0.075)

g.map_dataframe(facetErrorBars,"IsoNum","C/E","Sigma",ebDict,data=df_iso, hue='Evaluation',zorder=0,elinewidth=1.5,markeredgecolor='k',ls='')

g.set_xlabels('Isotope').set_titles("{col_name}")
#plt.legend(loc='lower left')
#f = plt.figure()
#f.subplots_adjust(right=0.8)
#plt.legend(loc='center left', bbox_to_anchor=(1.0, 0.5),title="Evaluation")
plt.legend(bbox_to_anchor=(1.0, 0.90), loc=2, borderaxespad=0., title="Evaluation")
plt.legend(bbox_to_anchor=(1.1, 0.90), loc=2, borderaxespad=0., title="Evaluation")
sns.despine(top=True)

titles = []
for sample in df_iso["Sample"].unique():
tmpTitle = "Sample {:s}".format(sample)
#print(sample, burnups[sample])
if(sample in burnups):
tmpTitle = tmpTitle + "\nBurnup: {:.2f} GWd/MTU".format(burnups[sample])
titles.append(tmpTitle)
Expand All @@ -134,6 +148,65 @@ def plotIsosFP(df_iso, assemName, seriesInfo, filename, burnups={},plotAspect=1.
plt.savefig(IMAGE_DEST + filename + '.' + OUTPUT_FORMAT)
plt.close()

def facetErrorBars(x, y, sigma, cDict, color=None, label=None, **kwargs):
# Add error bars to a FacetGrid / FactorPlot and match the bar color to the hue
data = kwargs.pop("data")
hue = kwargs.pop("hue")

ebColors = []
for eval in data[hue].unique():
ebColor = cDict[eval]
xVals = data[x].loc[data[hue] == eval].values
yVals = data[y].loc[data[hue] == eval].values
yErrs = data[sigma].loc[data[hue] == eval].values
plt.errorbar(xVals, yVals, xerr=None, yerr=yErrs, ecolor=ebColor, color=None,**kwargs)

def plotIsoBuFP(df_iso, labels, seriesInfo, isoName, sigmaLabel="Sigma", plotAspect=1.05, plotFit=False):

plt.clf()
isoMarkers = ['o', 'o', 's', '^', 'v', '*', 'H']

for idx in range(0,len(labels)):
# print("Label = ", labels[idx])
# Need to down-select either sample or evaluation as a series - plot C/E values for this...
df_plot = df_iso.loc[(df_iso["Evaluation"] == labels[idx])]
eBar = plt.errorbar(df_plot["Burnup"].values, df_plot["C/E"].values,df_plot[sigmaLabel].values,marker=isoMarkers[idx],label=labels[idx],ls='',markeredgecolor='k', markeredgewidth=0.15, ms=9)
if(plotFit):
eColor = eBar[0].get_color()
fitParams = np.polyfit(df_plot["Burnup"].values, df_plot["C/E"].values, 1, w=1.0/df_plot["Sigma"].values, full=False)


buVals = df_plot["Burnup"].values
ybar = np.mean(df_plot["C/E"].values)
ssreg = np.sum(((np.poly1d(fitParams)(buVals)-ybar)/df_plot["Sigma"])**2)
sstot = np.sum(((df_plot["C/E"].values-ybar)/df_plot["Sigma"])**2)
#print(ssreg,sstot)
r2 = ssreg/sstot
print("Evaluation: ", labels[idx], "Slope: {:.4f}, {:.4f}, R^2={:.4f}".format(fitParams[0], fitParams[1], r2))

buVals.sort()
plt.plot(buVals, np.poly1d(fitParams)(buVals), color=eColor,lw=0.75,ls='--')

sns.despine(top=True)
plt.xlabel("Burnup (GWd/MTU)")
plt.ylabel("C/E")
plt.title(isoName + " C/E variation with burnup")
plt.xlim(np.floor(np.min(df_iso["Burnup"].values)-2.0), np.ceil(np.max(df_iso["Burnup"].values))+2.0)

lgd = plt.legend(bbox_to_anchor=(1.05, 0.5), loc='center left', borderaxespad=0., title="Evaluation")

filename = isoName + "_CE_BU"

#
if SHOW_PLOTS:
plt.tight_layout()
plt.show()
else:
#plt.savefig(IMAGE_DEST + filename + '.' + OUTPUT_FORMAT)
plt.savefig(IMAGE_DEST + filename + '.' + OUTPUT_FORMAT,bbox_extra_artists=(lgd,), bbox_inches='tight')
plt.close()


def generateLatexTables(df_iso, df_measured, filename):
#texTable = "\\begin{table}[htb] \n\\centering \n"
for sample in df_iso["Sample"].unique():
Expand Down Expand Up @@ -184,121 +257,159 @@ def isoToLatex(isoName):
latexIso = RX_ISO3.sub(r'${}^{\2}\mathrm{\1}$',isoName)
#print(latexIso)
return latexIso


#================================================
# Measurement data
#===============================================

labels = ["ENDF VII.0","ENDF VII.1","ENDF VII.0 + mod-Wright","ENDF 7.0 + mod-Mughabghab", "ENDF VII.1 + mod-Mughabghab"]
labels = ["ENDF VII.0","ENDF VII.1","ENDF VII.0 + mod-Wright","ENDF VII.0 + mod-Mughabghab", "ENDF VII.1 + mod-Mughabghab", "JEFF-3.2 + mod-Mughabghab"]
# Default series meta-info
defaultSeriesInfo = [{ 'markersize':8,'marker':'o','label':'ENDF VII.0 (nominal)' },
{ 'markersize':8,'marker':'o','label':'ENDF/VII.1 (nominal)' },
{ 'markersize':8,'marker':'^','label':'VII.0 + mod. Wright' },
{ 'markersize':8,'marker':'v','label':'VII.0 + mod. Mughabghab' },
{ 'markersize':8,'marker':'>','label':'VII.1 + mod. Mughabghab' },
{ 'markersize':8,'marker':'h','label':'JEFF-3.2' },
{ 'markersize':8,'marker':'d','label':'JEFF-3.2 + mod. Mugabghab' }
]

df_measurements = pd.read_csv('Measurements.csv')
cleanupIsoNames(df_measurements)
df_measurements.sort_values(["Assembly","Sample","Isotope"])
#print(df_measurements[["Sample","Isotope","Element","Inventory"]])


#================================================
# TMI-1 NJ05YU
#================================================

tmiTitle = "TMI-1 NJ05YU"
euFilename = "TMI1_Eu_CE"
gdFilename = "TMI1_Gd_CE"
df_TMI1 = pd.read_csv('TMI1_inventories.csv')
cleanupIsoNames(df_TMI1)

#df_TMI1.sort_values(["Assembly","Sample","Isotope","Evaluation"])
TMI_measured = df_measurements[df_measurements["Assembly"] == "TMI-1 NJ05YU"]
calculateCE(df_TMI1, TMI_measured, "mg / gIU")

#pd.set_option('display.max_rows', 500)
#print(df_TMI1[["Sample","Evaluation","Isotope","Element","C/E"]])

TMI_burnups = {'AG536-C2D1': 53.5, 'AG536-C2D2':52.7, 'AG616-A1':45.9, 'AG616-B1':55.0, 'AG616-B2':52.4}
df_TMI1["Burnup"] = df_TMI1["Sample"].map(TMI_burnups)

isEu = df_TMI1.loc[(df_TMI1["Element"] == "Eu") & (df_TMI1["Isotope"] != "Eu-152")]
isGd = df_TMI1.loc[(df_TMI1["Element"] == "Gd") & (df_TMI1["Isotope"] != "Gd-152") & (df_TMI1["Isotope"] != "Gd-157")]
isEu_TMI = df_TMI1.loc[(df_TMI1["Element"] == "Eu") & (df_TMI1["Isotope"] != "Eu-152") & (df_TMI1["Evaluation"] != "JEFF-3.2")]
isGd_TMI = df_TMI1.loc[(df_TMI1["Element"] == "Gd") & (df_TMI1["Isotope"] != "Gd-152") & (df_TMI1["Isotope"] != "Gd-157") & (df_TMI1["Evaluation"] != "JEFF-3.2")]

tmiTitle = "TMI-1 NJ05YU"
euFilename = "TMI1_Eu_CE"
gdFilename = "TMI1_Gd_CE"
#pd.set_option('display.max_rows', 500)
#print(df_TMI1[["Sample","Evaluation","Isotope","Element","C/E"]])

TMI_burnups = {'AG536-C2D1': 53.5, 'AG536-C2D2':52.7, 'AG616-A1':45.9, 'AG616-B1':55.0, 'AG616-B2':52.4}
plotIsosFP(isEu, tmiTitle, defaultSeriesInfo, euFilename, TMI_burnups, plotAspect=1.25)
plotIsosFP(isGd, tmiTitle, defaultSeriesInfo, gdFilename, TMI_burnups, plotAspect=1.25)
plotIsosFP(isEu_TMI, tmiTitle, defaultSeriesInfo, euFilename, TMI_burnups, plotAspect=1.25,colWrap=3)
plotIsosFP(isGd_TMI, tmiTitle, defaultSeriesInfo, gdFilename, TMI_burnups, plotAspect=1.25,colWrap=3)
generateLatexTables(df_TMI1, TMI_measured, "TMI1_NJ05YU")

#================================================
# Calvert Cliffs
#================================================

ccTitle = "Calvert Cliffs"
euFilename = "CC_Eu_CE"
gdFilename = "CC_Gd_CE"
df_CC = pd.read_csv('CC_inventories.csv')
cleanupIsoNames(df_CC)

CC_measured = df_measurements[df_measurements["Assembly"] == "Calvert Cliffs"]
calculateCE(df_CC, CC_measured, "g / g Uinit")

isEu = df_CC.loc[(df_CC["Element"] == "Eu") & (df_CC["Isotope"] != "Eu-152") & (df_CC["Isotope"] != "Eu-151")]
isGd = df_CC.loc[(df_CC["Element"] == "Gd") & (df_CC["Isotope"] != "Gd-160")]
ccTitle = "Calvert Cliffs"
euFilename = "CC_Eu_CE"
gdFilename = "CC_Gd_CE"
CC_burnups = {"87-63": 44.34, "87-72":37.12, "87-81":27.35}
df_CC["Burnup"] = df_CC["Sample"].map(CC_burnups)

isEu_CC = df_CC.loc[(df_CC["Element"] == "Eu") & (df_CC["Isotope"] != "Eu-152") & (df_CC["Isotope"] != "Eu-151") & (df_CC["Evaluation"] != "JEFF-3.2")]
isGd_CC = df_CC.loc[(df_CC["Element"] == "Gd") & (df_CC["Isotope"] != "Gd-160") & (df_CC["Evaluation"] != "JEFF-3.2")]

plotIsosFP(isEu, ccTitle, defaultSeriesInfo, euFilename, CC_burnups, plotAspect=1.25)
plotIsosFP(isGd, ccTitle, defaultSeriesInfo, gdFilename, CC_burnups, plotAspect=1.25)
plotIsosFP(isEu_CC, ccTitle, defaultSeriesInfo, euFilename, CC_burnups, plotAspect=1.25)
plotIsosFP(isGd_CC, ccTitle, defaultSeriesInfo, gdFilename, CC_burnups, plotAspect=1.25)
generateLatexTables(df_CC, CC_measured, "CalvertCliffs")

#================================================
# REBUS
#================================================

rebusTitle = "REBUS"
euFilename = "REBUS_Eu_CE"
df_REBUS = pd.read_csv('REBUS_inventories.csv')
cleanupIsoNames(df_REBUS)

REBUS_measured = df_measurements[df_measurements["Assembly"] == "REBUS"]
calculateCE(df_REBUS, REBUS_measured, "g / gUI")

isEu = df_REBUS.loc[(df_REBUS["Element"] == "Eu") & (df_REBUS["Isotope"] != "Eu-152") & (df_REBUS["Isotope"] != "Eu-151")]
rebusTitle = "REBUS"
euFilename = "REBUS_Eu_CE"
REBUS_burnups = {"GKN-II":59.656}
plotIsosFP(isEu, rebusTitle, defaultSeriesInfo, euFilename, REBUS_burnups)
df_REBUS["Burnup"] = df_REBUS["Sample"].map(REBUS_burnups)

isEu_REBUS = df_REBUS.loc[(df_REBUS["Element"] == "Eu") & (df_REBUS["Isotope"] != "Eu-152") & (df_REBUS["Isotope"] != "Eu-151") & (df_REBUS["Evaluation"] != "JEFF-3.2")]

plotIsosFP(isEu_REBUS, rebusTitle, defaultSeriesInfo, euFilename, REBUS_burnups, fontScale=1.25)
generateLatexTables(df_REBUS, REBUS_measured, "REBUS")

#================================================
# ARIANE
#================================================

ARIANETitle = "ARIANE"
euFilename = "ARIANE_Eu_CE"
df_ARIANE = pd.read_csv('ARIANE_inventories.csv')
cleanupIsoNames(df_ARIANE)

ARIANE_measured = df_measurements[df_measurements["Assembly"] == "ARIANE"]
calculateCE(df_ARIANE, ARIANE_measured, "g / gUI")
#(df_ARIANE["Isotope"] != "Eu-152") & (df_ARIANE["Isotope"] != "Eu-151")
isEu = df_ARIANE.loc[(df_ARIANE["Element"] == "Eu") & (df_ARIANE["Isotope"] != "Eu-152") & (df_ARIANE["Isotope"] != "Eu-151")]
ARIANETitle = "ARIANE"
euFilename = "ARIANE_Eu_CE"

ARIANE_burnups = {"GU1":53.331}
plotIsosFP(isEu, ARIANETitle, defaultSeriesInfo, euFilename, ARIANE_burnups)
df_ARIANE["Burnup"] = df_ARIANE["Sample"].map(ARIANE_burnups)

isEu_ARIANE = df_ARIANE.loc[(df_ARIANE["Element"] == "Eu") & (df_ARIANE["Isotope"] != "Eu-152") & (df_ARIANE["Isotope"] != "Eu-151") & (df_ARIANE["Evaluation"] != "JEFF-3.2")]

plotIsosFP(isEu_ARIANE, ARIANETitle, defaultSeriesInfo, euFilename, ARIANE_burnups, fontScale=1.25)
generateLatexTables(df_ARIANE, ARIANE_measured, "ARIANE")


#================================================
# CANDU 28-element / Pickering-A 19558C
#================================================

CANDUTitle = "Pickering-A 19558C (CANDU 28-element)"
euFilename = "CANDU_Eu_CE"
df_CANDU = pd.read_csv('CANDU_inventories.csv')
cleanupIsoNames(df_CANDU)

CANDU_measured = df_measurements[df_measurements["Assembly"] == "CANDU28"]
calculateCE(df_CANDU, CANDU_measured, "Calculated")

#(df_CANDU["Isotope"] != "Eu-152") & (df_CANDU["Isotope"] != "Eu-151")
isEu = df_CANDU.loc[(df_CANDU["Element"] == "Eu") & (df_CANDU["Isotope"] != "Eu-152") & (df_CANDU["Isotope"] != "Eu-151")]
CANDUTitle = "Pickering-A 19558C (CANDU 28-element)"
euFilename = "CANDU_Eu_CE"
CANDU_burnups = {"Pickering-A 19558C":9.208}
plotIsosFP(isEu, CANDUTitle, defaultSeriesInfo, euFilename, CANDU_burnups)
df_CANDU["Burnup"] = df_CANDU["Sample"].map(CANDU_burnups)

isEu_CANDU = df_CANDU.loc[(df_CANDU["Element"] == "Eu") & (df_CANDU["Isotope"] != "Eu-152") & (df_CANDU["Isotope"] != "Eu-151") & (df_CANDU["Evaluation"] != "JEFF-3.2")]

plotIsosFP(isEu_CANDU, CANDUTitle, defaultSeriesInfo, euFilename, CANDU_burnups, fontScale=1.25)
generateLatexTables(df_CANDU, CANDU_measured, "CANDU")

#plotIsos(isEu, labels, defaultSeriesInfo, CANDUTitle, euFilename, CANDU_burnups)
#plotIsos(isEu, labels, defaultSeriesInfo, CANDUTitle, euFilename, CANDU_burnups)

#================================================
# Burnup trend analysis
#================================================

pd.set_option('display.max_rows', 500)

df_Master = pd.concat([df_TMI1, df_CC, df_REBUS, df_ARIANE], ignore_index=True)
df_Master_Eu153 = df_Master.loc[(df_Master["Isotope"] == "Eu-153") & ((df_Master["Evaluation"] == "ENDF VII.0") | (df_Master["Evaluation"] == "ENDF VII.1"))]
df_Master_Eu154 = df_Master.loc[(df_Master["Isotope"] == "Eu-154") & (df_Master["Evaluation"] != "JEFF-3.2")]
df_Master_Eu155 = df_Master.loc[(df_Master["Isotope"] == "Eu-155") & (df_Master["Evaluation"] != "JEFF-3.2")]
df_Master_Gd154 = df_Master.loc[(df_Master["Isotope"] == "Gd-154") & (df_Master["Evaluation"] != "JEFF-3.2")]
df_Master_Gd155 = df_Master.loc[(df_Master["Isotope"] == "Gd-155") & (df_Master["Evaluation"] != "JEFF-3.2")]


labels_Eu153 = labels[0:2]
#["ENDF VII.0","ENDF VII.1","ENDF VII.0 + mod-Wright","ENDF VII.0 + mod-Mughabghab", "ENDF VII.1 + mod-Mughabghab", "JEFF-3.2 + mod-Mughabghab"]

plotIsoBuFP(df_Master_Eu153, labels_Eu153, defaultSeriesInfo, "Eu-153",plotAspect=2.5)
plotIsoBuFP(df_Master_Eu154, labels, defaultSeriesInfo, "Eu-154",plotAspect=2.5, plotFit=True)
plotIsoBuFP(df_Master_Eu155, labels, defaultSeriesInfo, "Eu-155",plotAspect=2.5)
plotIsoBuFP(df_Master_Gd154, labels, defaultSeriesInfo, "Gd-154",plotAspect=2.5)
plotIsoBuFP(df_Master_Gd155, labels, defaultSeriesInfo, "Gd-155",plotAspect=2.5)
#print(df_EuMaster[["Sample","Isotope","Evaluation","Burnup","C/E"]])

0 comments on commit 2bebfaf

Please sign in to comment.