Skip to content

Commit

Permalink
Automate running sensitivity analysis using UQTk
Browse files Browse the repository at this point in the history
Handle renaming of compsets for SP, CROP and FATES simulations
  • Loading branch information
dmricciuto committed May 19, 2021
1 parent 7db62b1 commit 8bfaa1c
Show file tree
Hide file tree
Showing 5 changed files with 257 additions and 34 deletions.
31 changes: 31 additions & 0 deletions UQTk_scripts/prepare_env.x
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
#!/bin/bash

export UQTK_INS=/home/dmricciuto/models/Uncertainty-Quantification/Param-UQ-Workflow/UQTk_v3.0.2-install

if [[ $HOSTNAME == *"titan"* ]] ; then

printf "Seems to be on Titan: loading python_anaconda/2.3.0, gcc, cmake/2.8.11.2 \n"
module load python_anaconda/2.3.0
module load gcc
module load cmake/2.8.11.2

elif [[ $HOSTNAME == "eos"* ]] ; then


printf "Seems to be on Eos: loading python_anaconda/2.3.0, gcc, cmake/2.8.11.2 \n"
module load python_anaconda/2.3.0
module load gcc
module load cmake/2.8.11.2


elif [[ $HOSTNAME == "arsenal" ]]; then

echo "Setting up arsenal"


elif [[ $HOSTNAME == "lofty" ]]; then

echo "Setting up lofty"


fi
152 changes: 152 additions & 0 deletions UQTk_scripts/run_sensitivity.x
Original file line number Diff line number Diff line change
@@ -0,0 +1,152 @@
#!/bin/bash -e

###############################################################################################################
# Set necessary paths and variables
###############################################################################################################
UQSCRIPTS_PATH=`dirname "$0"`

if [ ! -f prepare_env.x ]; then
echo "=== Sourcing ${UQSCRIPTS_PATH}/prepare_env.x "
echo "=== please edit for your computer-specific paths."
source ${UQSCRIPTS_PATH}/prepare_env.x
else
echo "Sourcing ./prepare_env.x : please edit for your computer-specific paths."
source ./prepare_env.x
fi


## Set other relevant paths and variable
#export UQTKv2_PATH=${PWD}/UQTk_v2.2/src_cpp
#export UQTK_INS=${UQTKv2_PATH}
export UQTK_BIN=${UQTK_INS}/bin
export UQPC=${UQTK_INS}/examples/uqpc

###############################################################################################################
# Hardwired parameters to play with
###############################################################################################################

# Directory where ensemble is run
# Expects files param_range.txt, ptrain.dat, ytrain.dat, pval.dat, yval.dat in this directory
PRERUN_DIR=${UQSCRIPTS_PATH}/data #alm_5d_example or alm_5d_example_nz, or create your own folder with prerun data

# Surrogate construction method
# Options are lsq (least-squares) or bcs (Bayesian compressive sensing)
METHOD=bcs

# Select the order of the surrogate
ORDER=4

# Number of cores
NCORES=1

###############################################################################################################
# Sanity checks
###############################################################################################################

# Sanity check on method
if [[ "$METHOD" != "lsq" && "$METHOD" != "bcs" ]]; then
echo "METHOD=$METHOD is not recognized (choose lsq or bcs)"
exit
fi

###############################################################################################################
# Copy from previous runs
###############################################################################################################

echo "Getting the simulation results from ${PRERUN_DIR}"
for file in ${PRERUN_DIR}/*; do
rm -f `basename $file`
cp $file .
done

# Scale
${UQPC}/scale.x ptrain.dat from param_range.txt qtrain.dat
${UQPC}/scale.x pval.dat from param_range.txt qval.dat

NTRAIN=`awk 'END{print NR}' ptrain.dat`
NVAL=`awk 'END{print NR}' pval.dat`
SAMPLING=rand
# Number of outputs
NOUT=`awk 'NR==1{print NF}' ytrain.dat`
echo "Number of output QoIs : $NOUT "


###############################################################################################################
# Compute the surrogate
###############################################################################################################

rm -rf results.pk

if [ "$NCORES" -eq "1" ]; then
echo "Running in a serial mode"
${UQPC}/uq_pc.py -r offline_post -p param_range.txt -m $METHOD -s $SAMPLING -n $NTRAIN -v $NVAL -t $ORDER

else
echo "Running in parallel on $NCORES cores"
echo "Creating tasks file"
echo -n"" > tasks
rm -rf task_*
for ((outid=1;outid<=$NOUT;outid++)); do
mkdir task_$outid
for f in *train.dat *val.dat param_range.txt; do
ln -sf ../$f task_$outid/$f
done
let outidd=$NOUT+1-$outid #$outid
echo "${UQPC}/uq_pc.py out.log -a $outidd -r offline_post -p param_range.txt -m $METHOD -s $SAMPLING -n $NTRAIN -v $NVAL -t $ORDER" >> tasks
done
echo "Running tasks in task_* folders"
${UQTK_INS}/PyUQTk/multirun/multirun.py tasks $NCORES > uqmulti.log

${UQPC}/join_results.py
fi

if [ -f results.pk ]; then
echo "Results are packed into results.pk"
else
echo "Something is wrong. The file results.pk not found."
fi

###############################################################################################################
# Plot results
###############################################################################################################

# Plot data-versus-model for surrogate accuracy assessment
${UQPC}/plot.py dm training validation

# Plot runId versus data and model for surrogate accuracy assessment
${UQPC}/plot.py idm training
${UQPC}/plot.py idm validation

# Plot sensitivities (multi-output bars)
${UQPC}/plot.py sens main
${UQPC}/plot.py sens total

# Plot total sensitivities (matrix plots for all outputs and most relevant inputs)
${UQPC}/plot.py sensmat total

# Plot main and joint sensitivities for all outputs (circular plots)
${UQPC}/plot.py senscirc

# Plot high-d representation of multiindex
${UQPC}/plot.py mindex

# Compute and plot output PDFs
if [[ "$METHOD" != "bcs" ]]; then # requires full multiindex to work, therefore won't compute it with bcs
${UQPC}/plot.py pdf
fi

###############################################################################################################
# Cleanup
###############################################################################################################

# Save the input-outputs/sensitivities and logs just in case
tar -czf logs.tar $(ls -d *.log 2>/dev/null)
tar -czf inouts_sens.tar *train.dat *val.dat *sens*dat mindex_output*dat pcf_output*dat

# Clean the rest
rm -rf tmp cfs *.dat *.log pcf *.txt





43 changes: 31 additions & 12 deletions manage_ensemble.py
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,9 @@
help = 'CNP mode - initialize P pools')
parser.add_option("--site", dest="site", default='parm_list', \
help = 'Site name')
parser.add_option('--run_sa', dest="run_sa", default=True, action="store_true", \
help = 'Run sensitivity analysis using UQTk')

(options, args) = parser.parse_args()

options.n = int(options.n)
Expand Down Expand Up @@ -96,17 +99,20 @@ def postproc(myvars, myyear_start, myyear_end, myday_start, myday_end, myavg, \
myindex = mypft[index]
hol_add = 17
if (os.path.exists(fname)):
mydata = nffun.getvar(fname,v)
if (len(mydata) < 365):
mydata = np.zeros([365,34], np.float)+np.NaN
mydata = nffun.getvar(fname,v)
if (len(mydata) < 10):
npy = 1
elif (len(mydata) >= 365): #does not currently allow hourly
npy = 365
else:
mydata = np.zeros([365,34], np.float)+np.NaN
print(fname)
mydata = np.zeros([npy,34], np.float)+np.NaN
#get output and average over days/years
n_days = myday_end[index]-myday_start[index]+1
ndays_total = ndays_total + n_days
#get number of timesteps per output file

if (('20TR' in case or (not '1850' in case)) and (not 'ED' in case)): #Transient assume daily ouput
if (npy == 365):
for d in range(myday_start[index]-1,myday_end[index]):
if ('US-SPR' in case):
output.append(0.25*(mydata[d][myindex+hol_add]*myfactor[index] \
Expand All @@ -115,7 +121,7 @@ def postproc(myvars, myyear_start, myyear_end, myday_start, myday_end, myavg, \
else:
output.append(mydata[d][myindex]*myfactor[index] \
+myoffset[index])
else: #Assume annual output (ignore days)
elif (npy == 1): #Assume annual output (ignore days)
for d in range(myday_start[index]-1,myday_end[index]): #28-38 was myindex
if ('SCPF' in v):
output.append(sum(mydata[0,28:38])/10.0*myfactor[index]+myoffset[index])
Expand All @@ -125,6 +131,8 @@ def postproc(myvars, myyear_start, myyear_end, myday_start, myday_end, myavg, \
data[thiscol] = sum(output[(i*myavg[index]):((i+1)*myavg[index])])/myavg[index]
thiscol=thiscol+1
index=index+1

#get the parameters
if (options.microbe):
pfname = rundir+'microbepar_in'
pnum=0
Expand Down Expand Up @@ -176,17 +184,24 @@ def postproc(myvars, myyear_start, myyear_end, myday_start, myday_end, myavg, \
elif (int(ppfts[pnum]) > 0):
parms[pnum] = mydata[int(ppfts[pnum])]
elif (int(ppfts[pnum]) == 0):
parms[pnum] = mydata[int(ppfts[pnum])]
try:
parms[pnum] = mydata[int(ppfts[pnum])]
except:
parms[pnum] = mydata
else:
parms[pnum] = mydata[...]
try:
parms[pnum] = mydata[0]
except:
parms[pnum] = mydata
else: #Regular parameter file
mydata = nffun.getvar(pfname,p)
if (int(ppfts[pnum]) > 0):
parms[pnum] = mydata[int(ppfts[pnum])]
elif(int(ppfts[pnum]) == 0):
parms[pnum] = mydata[0]
else:
parms[pnum] = mydata
elif(int(ppfts[pnum]) <= 0):
try:
parms[pnum] = mydata[0]
except:
parms[pnum] = mydata
pnum=pnum+1

return ierr
Expand Down Expand Up @@ -521,6 +536,10 @@ def calc_costfucntion(constraints, thisjob, runroot, case):
myoutput.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_sa):
os.system('cp UQTk_scripts/*.x '+UQ_output+'/')
os.chdir(UQ_output)
os.system('./run_sensitivity.x')

MPI.Finalize()

Expand Down
8 changes: 5 additions & 3 deletions runcase.py
Original file line number Diff line number Diff line change
Expand Up @@ -316,10 +316,12 @@
#check whether model named clm or elm
if (os.path.exists(options.csmdir+'/components/elm')):
mylsm='ELM'
model_name='elm'
if (options.res == 'CLM_USRDAT'):
options.res = 'ELM_USRDAT'
else:
mylsm='CLM'
model_name='clm2'

#machine info: cores per node
ppn=1
Expand Down Expand Up @@ -1724,17 +1726,17 @@
cmd = mpicmd+' -np '+str(np_total)+' python manage_ensemble.py ' \
+'--case '+casename+' --runroot '+runroot+' --n_ensemble '+str(nsamples)+' --ens_file '+ \
options.ensemble_file+' --exeroot '+exeroot+' --parm_list '+options.parm_list+' --cnp '+cnp + \
' --site '+options.site
' --site '+options.site+' --model_name '+model_name
elif (('titan' in options.machine or 'eos' in options.machine) and int(options.ninst) == 1):
cmd = 'aprun -n '+str(np_total)+' python manage_ensemble.py ' \
+'--case '+casename+' --runroot '+runroot+' --n_ensemble '+str(nsamples)+' --ens_file '+ \
options.ensemble_file+' --exeroot '+exeroot+' --parm_list '+options.parm_list+' --cnp '+cnp + \
' --site '+options.site
' --site '+options.site+' --model_name '+model_name
elif ('anvil' in options.machine or 'compy' in options.machine or 'cori' in options.machine):
cmd = 'srun -n '+str(np_total)+' python manage_ensemble.py ' \
+'--case '+casename+' --runroot '+runroot+' --n_ensemble '+str(nsamples)+' --ens_file '+ \
options.ensemble_file+' --exeroot '+exeroot+' --parm_list '+options.parm_list+' --cnp '+cnp + \
' --site '+options.site
' --site '+options.site+' --model_name '+model_name
if (options.constraints != ''):
cmd = cmd + ' --constraints '+options.constraints
if (options.postproc_file != ''):
Expand Down
Loading

0 comments on commit 8bfaa1c

Please sign in to comment.