Add capability to perform GSA using Python SALib

MCMC.py

@@ -52,7 +52,7 @@ def posterior(parms):
 
 #-------------------------------- MCMC ------------------------------------------------------
 
-def MCMC(parms, nevals, type='uniform', nburn=1000, burnsteps=10, default=[]):
+def MCMC(parms, nevals, type='uniform', nburn=1000, burnsteps=10, default_output=[]):
     #Metropolis-Hastings Markov Chain Monte Carlo with adaptive sampling
     post_best = -99999
     post_last = -99999
@@ -66,7 +66,6 @@ def MCMC(parms, nevals, type='uniform', nburn=1000, burnsteps=10, default=[]):
     chain_burn = np.zeros((nparms,nevals))
     output     = np.zeros((model.nobs,nevals))
     mycov      = np.zeros((nparms,nparms))
-
     for p in range(0,nparms):
         #Starting step size = 1% of prior range
         #parm_step[p] = 2.4**2/nparms * (model.pmax[p]-model.pmin[p])
@@ -162,14 +161,14 @@ def MCMC(parms, nevals, type='uniform', nburn=1000, burnsteps=10, default=[]):
                 post_best = post
                 parms_best = parms
                 print(post_best)
-                output_best = model.output
+                output_best = model.output.flatten()
 
         #populate the chain matrix
         for j in range(0,nparms):
             chain[j][i] = parms[j]
         chain[nparms][i] = post_last
         for j in range(0,model.nobs):
-            output[j,i] = model.output[j]
+            output[j,i] = model.output[0,j]
 
         if (i % 1000 == 0):
             print(' -- '+str(i)+' --\n')
@@ -249,19 +248,18 @@ def MCMC(parms, nevals, type='uniform', nburn=1000, burnsteps=10, default=[]):
 
     #make prediction plots
     obs_set = list(set(model.obs_name))
-    print(obs_set)
     for s in range(0,len(obs_set)):
       thisob = [ix for ix, value in enumerate(model.obs_name) if value == obs_set[s]]
       fig = plt.figure()
       ax=fig.add_subplot(111)
-      x = np.cumsum(np.ones([len(thisob)],np.float))
+      x = np.cumsum(np.ones([len(thisob)],float))
       ax.errorbar(x,model.obs[thisob], yerr=model.obs_err[thisob], label='Observations')
       ax.plot(x,output_best[thisob],'r', label = 'Model best')
       ax.plot(x,output_sorted[thisob,int(0.025*(nevals-nburn*burnsteps))], \
                  'k--', label='Model 95% CI')
       ax.plot(x,output_sorted[thisob,int(0.975*(nevals-nburn*burnsteps))],'k--')
       if (options.parm_default != ''):
-        ax.plot(x,default, 'g', label='Default')
+        ax.plot(x,default_output[thisob], 'g', label='Default')
         #plt.xlabel(model.xlabel)
         #plt.ylabel(model.ylabel)
       box = ax.get_position()
@@ -280,10 +278,9 @@ def MCMC(parms, nevals, type='uniform', nburn=1000, burnsteps=10, default=[]):
   parms_default = np.loadtxt(options.parm_default)
   model.run(parms_default)
   default_output = model.output.flatten()
-
   parms = MCMC(model.pdef, int(options.nevals), burnsteps=int(options.burnsteps), \
                           nburn=int(options.nevals)/(2*int(options.burnsteps)), \
-                          default=default_output)
+                          default_output=default_output)
 else:
   #run MCMC
   parms = MCMC(model.pdef, int(options.nevals), burnsteps=int(options.burnsteps), \

manage_ensemble.py

@@ -156,6 +156,7 @@ def postproc(myvars, myyear_start, myyear_end, myday_start, myday_end, myavg, \
         pnum=pnum+1
     else:
       pfname = baserundir+'clm_params_'+str(100000+thisjob)[1:]+'.nc'
+      #pfname_def = baserundir+'clm_params.nc'
       fpfname = baserundir+'fates_params_'+str(100000+thisjob)[1:]+'.nc'
       sfname = baserundir+'surfdata_'+str(100000+thisjob)[1:]+'.nc'
       pnum=0
@@ -376,20 +377,25 @@ def postproc(myvars, myyear_start, myyear_end, myday_start, myday_end, myavg, \
         myoutput = open(UQ_output+'/data/outnames.txt', 'w')
         vlast=''
         for v in myvars:
-          if v != vlast:
-            vcount=0
-            vlast=v
-          if (myvars.count(v) > 1):
-            myoutput.write(v+'_'+str(vcount)+'\n')
-            vcount=vcount+1
-          else:
-            myoutput.write(v+'\n')
+          #if v != vlast:
+          #  vcount=0
+          #  vlast=v
+          #if (myvars.count(v) > 1):
+          #  myoutput.write(v+'_'+str(vcount)+'\n')
+          #  vcount=vcount+1
+          #else:
+          myoutput.write(v+'\n')
           eden_header=eden_header+v+','
         myoutput.close()
+
+        os.system('mkdir -p '+UQ_output+'/GSA')
         myoutput = open(UQ_output+'/data/param_range.txt', 'w')
+        myoutput2 = open(UQ_output+'/GSA/param_range.txt', 'w')
         for p in range(0,len(pmin)):
           myoutput.write(pmin[p]+' '+pmax[p]+'\n')
+          myoutput2.write(pnames[p]+' '+pmin[p]+' '+pmax[p]+'\n')
         myoutput.close()
+        myoutput2.close()
         print(np.hstack((parm_out,data_out)))
         np.savetxt(UQ_output+'/data/foreden.csv', np.hstack((parm_out,data_out)), delimiter=',', header=eden_header[:-1])
         if (options.run_uq):
@@ -406,6 +412,8 @@ def postproc(myvars, myyear_start, myyear_end, myday_start, myday_end, myavg, \
           #os.chdir('../..')
           #Create the surrogate model
           os.system('python surrogate_NN.py --case '+options.casename)
+          #Run the senstivity analysis using SALib
+          os.system('python run_GSA.py --case '+options.casename)
           if (max(myobs_err) > 0):
             #Run the MCMC calibration on surrogate model if data provided
             os.system('python MCMC.py --case '+options.casename+' --parm_list '+options.parm_list)
@@ -457,13 +465,15 @@ def postproc(myvars, myyear_start, myyear_end, myday_start, myday_end, myavg, \
                       lnd_in_new=open('lnd_in','w')
                       pst = str(100+plots[t])[1:]
                       for s in lnd_in_old:
-                        if ('finidat' in s):
+                        if ('finidat =' in s):
                           lnd_in_new.write(" finidat = './"+c+"."+options.model_name+".r.2015-01-01-00000.nc'\n")
                         elif ('metdata_bypass' in s):
                           lnd_in_new.write(s[:-2]+'/plot'+pst+"'\n")
                           if ('CO2' in treatments[t]):
                             lnd_in_new.write(' add_co2 = 500\n')
                             lnd_in_new.write(" startdate_add_co2 = '20160315'\n") 
+                        elif ('landuse_timeseries' in s):
+                          lnd_in_new.write(s.replace('plot07','plot'+pst))
                         else:
                           lnd_in_new.write(s)
                       lnd_in_old.close()

model_surrogate.py

@@ -38,8 +38,8 @@ def __init__(self,case=''):
         obsfile = open(UQdir+'/data/obs.dat')
         i=0
         for s in obsfile:
-          self.obs[i] = np.float(s.split()[0])
-          self.obs_err[i] = np.float(s.split()[1])
+          self.obs[i] = float(s.split()[0])
+          self.obs_err[i] = float(s.split()[1])
           i=i+1
         obsfile.close()
         outnamesfile = open(UQdir+'/data/outnames.txt')
@@ -52,15 +52,24 @@ def __init__(self,case=''):
         self.qoi_good = np.loadtxt(UQdir+'/NN_surrogate/qoi_good.txt').astype(int)
 
     def run(self,parms):
-        parms_nn = np.zeros([1,self.nparms],float)
-        for p in range(0,self.nparms):
-          parms_nn[0,p] = (parms[p]-self.pmin[p])/(self.pmax[p]-self.pmin[p])
-        self.output = np.zeros([self.nobs])
-        output_temp = self.nnmodel.predict(parms_nn).flatten()
+        if ((parms).ndim == 1):
+            nsamples=1
+            theseparms=parms
+        else:
+            nsamples=parms.shape[0]
+        parms_nn = np.zeros([nsamples,self.nparms],float)
+        for n in range(0,nsamples):
+          if nsamples > 1:
+              theseparms = parms[n,:]
+          for p in range(0,self.nparms):
+            parms_nn[n,p] = (theseparms[p]-self.pmin[p])/(self.pmax[p]-self.pmin[p])
+        self.output = np.zeros([nsamples,self.nobs])
+        output_temp = self.nnmodel.predict(parms_nn)
+
         qgood=0
         for q in range(0,self.nobs):
           if (q in self.qoi_good):
-            self.output[q] = output_temp[qgood]*(self.yrange[1,q]-self.yrange[0,q])+self.yrange[0,q]
+            self.output[:,q] = output_temp[:,qgood]*(self.yrange[1,q]-self.yrange[0,q])+self.yrange[0,q]
             qgood=qgood+1
           else:
-            self.output[q] = self.yrange[1,q]
+            self.output[:,q] = self.yrange[1,q]

run_GSA.py

@@ -0,0 +1,80 @@
+import os, math, sys
+import numpy as np
+from optparse import OptionParser
+import model_surrogate as models
+import matplotlib
+import matplotlib.pyplot as plt
+matplotlib.use('Agg')
+
+parser = OptionParser()
+parser.add_option("--case", dest="casename", default="", \
+                  help="Name of case")
+(options, args) = parser.parse_args()
+
+GSAdir = './UQ_output/'+options.casename+'/GSA'
+os.system('mkdir -p '+GSAdir+'/analyses')
+
+os.system('python -m SALib.sample.saltelli -n 8192 -p '+GSAdir+'/param_range.txt -o '+ \
+                  GSAdir+'/Saltelli_samples.txt')
+#Create the model object
+model = models.MyModel(case=options.casename)
+
+#Get the parameter samples for GSA
+samples=np.loadtxt(GSAdir+'/Saltelli_samples.txt')
+
+#Run the NN surrogate model
+model.run(samples)
+
+np.savetxt(GSAdir+'/outputs.txt',model.output)
+
+#Run GSA for each QOI
+sens_main = np.zeros([model.nparms,model.nobs])
+sens_main_unc = np.zeros([model.nparms,model.nobs])
+sens_tot = np.zeros([model.nparms,model.nobs])
+sens_tot_unc = np.zeros([model.nparms,model.nobs])
+
+for n in range(0,model.nobs):
+  print(n)
+  os.system('python -m SALib.analyze.sobol --parallel -p '+GSAdir+'/param_range.txt -Y '+GSAdir+ \
+            '/outputs.txt -c '+str(n)+' > '+GSAdir+'/analyses/analysis_ob'+str(n)+'.txt')
+  myfile = open(GSAdir+'/analyses/analysis_ob'+str(n)+'.txt','r')
+  lnum=0
+  for s in myfile:
+    print(s)
+    if (lnum > 0 and lnum <= model.nparms):
+      sens_tot[lnum-1,n] = float(s.split()[1])
+      sens_tot_unc[lnum-1,n] = float(s.split()[2])
+    elif (lnum > model.nparms+1 and lnum <= model.nparms*2+1):
+      sens_main[lnum-2-model.nparms,n] = float(s.split()[1])
+      sens_main_unc[lnum-2-model.nparms,n] = float(s.split()[2])    
+    lnum=lnum+1
+  myfile.close()
+
+#Plot main sensitivity indices
+fig,ax = plt.subplots()
+x_pos = np.cumsum(np.ones(model.nobs))
+x_labels=model.obs_name.copy()
+for i in range(1,len(x_labels)):
+    if (model.obs_name[i] == model.obs_name[i-1]):
+        x_labels[i]=' '
+ax.bar(x_pos, sens_main[0,:], align='center', alpha=0.5)
+ax.set_xticks(x_pos)
+ax.set_xticklabels(x_labels, rotation=45)
+bottom=sens_main[0,:]
+for p in range(1,model.nparms):
+    ax.bar(x_pos, sens_main[p,:], bottom=bottom)
+    bottom=bottom+sens_main[p,:]
+plt.legend(model.parm_names)
+plt.savefig(GSAdir+'/sens_main.pdf')
+
+#Total sensitivity indices
+fig,ax = plt.subplots()
+ax.bar(x_pos, sens_tot[0,:], align='center', alpha=0.5)
+ax.set_xticklabels(x_labels, rotation=45)
+bottom=sens_tot[0,:]
+for p in range(1,model.nparms):
+    ax.bar(x_pos, sens_tot[p,:], bottom=bottom)
+    bottom=bottom+sens_tot[p,:]
+plt.legend(model.parm_names)
+plt.savefig(GSAdir+'/sens_tot.pdf')
+