-
Notifications
You must be signed in to change notification settings - Fork 1
/
plotSEN.py
108 lines (90 loc) · 3.53 KB
/
plotSEN.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
import pandas as pd
from matplotlib.backends.backend_pdf import PdfPages
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
import os
mpl.rcParams['font.sans-serif'] = 'Univers 57 Condensed'
mpl.rcParams['font.serif'] = 'Times'
mpl.rcParams['pdf.compression'] = 0
mpl.rcParams['pdf.fonttype'] = 42
#--figure text sizes
mpl.rcParams['legend.fontsize'] = 12
mpl.rcParams['axes.labelsize'] = 14
mpl.rcParams['xtick.labelsize'] = 12
mpl.rcParams['ytick.labelsize'] = 12
def SEN_reader(insenfile):
# function to read in a SEN file (insenfile) and return a two-level dictionary
# that includes first level keys as obsgroups and second level keys as parameters with composite sensitivities
# read in the entire SEN file first
parsfound = False
parnames = list()
obsgroupnames = list()
obsinds = dict()
insen = open(insenfile, 'r').readlines()
infilelen = len(insen)
allsen = dict()
# find where the parameter names are
for i, line in enumerate(insen):
if 'parameter name' in line.lower():
startpars = i + 1
parsfound = True
if parsfound:
tmp = line.strip()
if len(tmp) == 0:
endpars = i
break
# read in the parameter names
for i in np.arange(startpars, endpars):
parnames.append(insen[i].strip().split()[0])
NPAR = len(parnames)
# now find the sensitivities at the end of optimization
for i, line in enumerate(insen):
if 'COMPLETION OF OPTIMISATION PROCESS' in line.upper():
completion_line = i + 1
break
# now go through the file and read in the parameter
for i in np.arange(completion_line, len(insen)):
if 'composite sensitivities' in insen[i].lower() and 'prior' not in insen[i].lower():
tmpobsgrpname = insen[i].strip().split()[-2].replace('"', '')
if 'regul' not in tmpobsgrpname.lower():
obsgroupnames.append(tmpobsgrpname)
obsinds[obsgroupnames[-1]] = [i+3, i+3+NPAR]
for cobgrp in obsgroupnames:
strow = obsinds[cobgrp][0] + 1
endrow = obsinds[cobgrp][1] + 1
pars = []
sens = []
for i in np.arange(strow,endrow):
tmp = insen[i].strip().split()
pars.append(tmp[0])
sens.append(tmp[-1])
allsen[cobgrp] = dict(zip(pars,sens))
return allsen
###################################
######## MAIN STARTS HERE #########
###################################
#workingpath = '/Users/mnfienen/Documents/MODELING_CENTER/NACP/NACP15_ss_SENSITIVITY/'
workingpath = os.getcwd()
insenfile = os.path.join(workingpath, 'nacp15_ss_reg.sen')
sendict = SEN_reader(insenfile)
numpars_plot = 15
allsenplot = PdfPages(os.path.join(workingpath,'nacp15_ss_SENSITIVITIY.pdf'))
groupkeyfile = os.path.join(workingpath, 'GroupKeyNACP.txt')
groupkeydata = np.genfromtxt(groupkeyfile, dtype=None, names=True, delimiter=',')
groupkeys = dict(zip(groupkeydata['group'],groupkeydata['name']))
SENdf = pd.DataFrame(sendict)
for c in SENdf.columns:
print 'plotting sensitivity for group --> {0}'.format(c)
s = SENdf.ix[:,c].astype(float)
s.sort(ascending = False)
s_plot = s[0:numpars_plot]
s_plot /= max(s_plot)
fig = s_plot.plot(kind='bar', rot = 45)
fig.set_title('{0}: {1}'.format(c, groupkeys[c]))
plt.xlabel('Parameter')
plt.ylabel('Normalized Composite Sensitivity')
plt.tight_layout()
allsenplot.savefig()
plt.close()
allsenplot.close()