Sensitivity analysis with SALib

We have got the single parts now for the sensitivity analysis. We are now using the global sensitivity analysis methods of the Python package SALib, available on:

https://github.com/jdherman/SALib

As a start, we will test the sensitivity of the model at each drillhole position separately. As parameters, we will use the parameters of the fault events: dip, dip direction, and slip.


In [35]:
from IPython.core.display import HTML
css_file = 'pynoddy.css'
HTML(open(css_file, "r").read())


Out[35]:

In [1]:
%matplotlib inline

In [12]:
import sys, os
import matplotlib.pyplot as plt
import numpy as np
# adjust some settings for matplotlib
from matplotlib import rcParams
# print rcParams
rcParams['font.size'] = 15
# determine path of repository to set paths corretly below
os.chdir(r'/Users/flow/git/pynoddy/docs/notebooks/')
repo_path = os.path.realpath('../..')
import pynoddy.history
import pynoddy.experiment
rcParams.update({'font.size': 20})

Model set-up

We use the two-fault model from previous examples and assign parameter ranges with a dictionary:


In [3]:
reload(pynoddy.history)
reload(pynoddy.experiment)
# Start again with the original model
his_filename = "two_faults_sensi.his"
sa = pynoddy.experiment.SensitivityAnalysis(history = his_filename)

# Initialise list
param_stats = []

# Add one entry as dictionary with relevant properties:

# for event 2:
param_stats.append({'event' : 2, 'parameter' : 'Dip', 'min' : 55., 'max' : 65., 
                    'type' : 'normal', 'stdev' : 10., 'mean' : 60., 'initial' : 60.})
param_stats.append({'event' : 2, 'parameter' : 'Dip Direction', 'min' : 85., 'max' : 95., 
                    'type' : 'normal', 'stdev' : 10., 'mean' : 90., 'initial' : 90.})
param_stats.append({'event' : 2, 'parameter' : 'Slip', 'min' : 900., 'max' : 1100., 
                    'type' : 'normal', 'stdev' : 500, 'mean' : 1000., 'initial' : 1000.})
# for event 3:
param_stats.append({'event' : 3, 'parameter' : 'Dip', 'min' : 55., 'max' : 65., 
                    'type' : 'normal', 'stdev' : 10., 'mean' : 60., 'initial' : 60.})
param_stats.append({'event' : 3, 'parameter' : 'Dip Direction', 'min' : 265., 'max' : 275., 
                    'type' : 'normal', 'stdev' : 10., 'mean' : 270., 'initial' : 270.})
param_stats.append({'event' : 3, 'parameter' : 'Slip', 'min' : 900., 'max' : 1100., 
                    'type' : 'normal', 'stdev' : 500, 'mean' : 1000., 'initial' : 1000.})

sa.set_parameter_statistics(param_stats)


 STRATIGRAPHY
 FAULT
 FAULT

Define sampling lines

As before, we need to define points in the model (or lines) which we want to evaluate the sensitivity for:


In [4]:
# sa.add_sampling_line(5000, 3500, label = 'centre')
sa.add_sampling_line(2500, 3500, label = 'left')
# sa.add_sampling_line(7500, 3500, label = 'right')
# sa.add_sampling_line(4000, 3500, label = 'compare')

And, again, we "freeze" the base state for later comparison and distance caluclations:


In [5]:
sa.freeze()

Setting-up the parameter set

For use with SALib, we have to define a parameter set as a text file (maybe there is a different way directly in Python - something to figure out for the future). The sensitivity object has a method to do that automatically:


In [6]:
param_file = "params_file_tmp.txt"
sa.create_params_file(filename = param_file)

We now invoke the methods of the SALib library to generate parameter data sets that are required for the type of sensitivity analysis that we want to perform:


In [7]:
# import SALib method
from SALib.sample import saltelli

In [8]:
param_values = saltelli.sample(10, param_file, calc_second_order = True)

The object 'param_values' is a list of samples for the parameters that are defined in the model, in the order of appearance in param_file, e.g.:


In [9]:
param_values[0]


Out[9]:
array([   57.19726562,    85.96679688,  1003.7109375 ,    61.76757812,
         267.80273438,  1081.4453125 ])

Calculating distances for all parameter sets

We now need to create a model realisation for each of these parameter sets and calculate the distance between the realisation and the base model at the position of the defined sampling lines. As we are not (always) interested in keeping the results of all realisations, those steps are combined and only the calculated distance is retained (per default):


In [10]:
distances = sa.determine_distances(param_values = param_values)

In [45]:
# plot(sa.get_model_lines(model_type = 'base'))
plt.plot(sa.get_model_lines(model_type = 'current'))


Out[45]:
[<matplotlib.lines.Line2D at 0x111b8a2d0>]

In [11]:
# Just to check if we actualy did get different models:
plt.plot(distances, '.-k')
plt.title("Model distances")
plt.xlabel("Sensitivity step")
plt.ylabel("Distance")


Out[11]:
<matplotlib.text.Text at 0x104471710>

Sensitivity analysis

We can now analyse the sensitivity of the modelled stratigraphy along the defined vertical lines ("drillholes") with respect to the model parameters:


In [47]:
# save results
results_file = 'dist_tmp.txt'
np.savetxt(results_file, distances, delimiter=' ')

In [48]:
from SALib.analyze import sobol

In [49]:
Si = sobol.analyze(param_file, results_file, 
                   column = 0, 
                   conf_level = 0.95,
                   calc_second_order = True, 
                   print_to_console=False)

In [50]:
# create composite matrix for sensitivities
n_params = 6
comp_matrix = np.ndarray(shape = (n_params,n_params))
for j in range(n_params):
    for i in range(n_params):
        if i == j:
            comp_matrix[i,j] = Si['S1'][i]
        else:
            comp_matrix[i,j] = Si['S2'][i,j]
            comp_matrix[j,i] = Si['S2'][i,j]
            
# print comp_matrix

# define labels for figure: phi = dip, d = dip direction, s = slip, subscript = fault event
label_names = ["","$\phi_1$", "$d_1$", "$s_1$", "$\phi_2$", "$d_2$", "$s_2$"]

# Create figure
fig = plt.figure()
ax = fig.add_subplot(111)
im = ax.imshow(comp_matrix, interpolation='nearest', cmap='RdBu_r', 
               vmax = np.max(np.abs(comp_matrix)),
                vmin = -np.max(np.abs(comp_matrix)),

)

ax.yaxis.set_ticks_position("both")
ax.xaxis.set_ticks_position("top")

ax.set_xticklabels(label_names)
ax.set_yticklabels(label_names)

# ax.set_title("Sensitivities")

ax.set_xlabel("Parameter Sensitivities")

fig.colorbar(im)

plt.tight_layout()
# plt.savefig("two_fault_sensi.png")


Reading parameter ranges from file

So, now that we have all the required ingredients for the sensitivity analysis, we can make life a bit easier with more automation. First, instead of defining parameters in a dictionary as above, we can actually read them in from a csv file (e.g. saved from Excel as Windows-csv file).

In order to be read in correctly, the header should contain the labels:

  • 'event' : event id
  • 'parameter' : Noddy parameter ('Dip', 'Dip Direction', etc.)
  • 'min' : minimum value
  • 'max' : maximum value
  • 'initial' : initial value

In addition, it is possible to define PDF type and parameters. For now, the following settings are supported:

  • 'type' = 'normal'
  • 'stdev' : standard deviation
  • 'mean' : mean value (default: 'initial' value)

We can read in the parameters simply with:


In [51]:
reload(pynoddy.history)
reload(pynoddy.experiment)
# Start again with the original model
his_filename = "two_faults_sensi.his"
sa = pynoddy.experiment.SensitivityAnalysis(history = his_filename)
sa.load_parameter_file("params_fault_model.csv")


 STRATIGRAPHY
 FAULT
 FAULT

The only further aspect we need to define are the sampling lines:


In [52]:
# sa.add_sampling_line(5000, 3500, label = 'centre')
sa.add_sampling_line(2500, 3500, label = 'left')
# sa.add_sampling_line(7500, 3500, label = 'right')
# sa.add_sampling_line(4000, 3500, label = 'compare')

And then we know everything to perform the sensitivity analysis. The relevant steps are combined in one method:


In [53]:
sa.perform_analsis(10)

In [54]:
sa.plot_distances()



In [55]:
sa.plot_sensitivity_matrix()



In [56]:
# for event 2:
param_stats.append({'event' : 2, 'parameter' : 'Dip', 'min' : 55., 'max' : 65., 
                    'type' : 'normal', 'stdev' : 10., 'mean' : 60., 'initial' : 60.})
param_stats.append({'event' : 2, 'parameter' : 'Dip Direction', 'min' : 85., 'max' : 95., 
                    'type' : 'normal', 'stdev' : 10., 'mean' : 90., 'initial' : 90.})
param_stats.append({'event' : 2, 'parameter' : 'Slip', 'min' : 900., 'max' : 1100., 
                    'type' : 'normal', 'stdev' : 500, 'mean' : 1000., 'initial' : 1000.})
# for event 3:
param_stats.append({'event' : 3, 'parameter' : 'Dip', 'min' : 55., 'max' : 65., 
                    'type' : 'normal', 'stdev' : 10., 'mean' : 60., 'initial' : 60.})
param_stats.append({'event' : 3, 'parameter' : 'Dip Direction', 'min' : 265., 'max' : 275., 
                    'type' : 'normal', 'stdev' : 10., 'mean' : 270., 'initial' : 270.})
param_stats.append({'event' : 3, 'parameter' : 'Slip', 'min' : 900., 'max' : 1100., 
                    'type' : 'normal', 'stdev' : 500, 'mean' : 1000., 'initial' : 1000.})

In [57]:
sa.param_stats


Out[57]:
[{'event': 2.0,
  'initial': 60.0,
  'max': 65.0,
  'mean': 60.0,
  'min': 55.0,
  'parameter': 'Dip',
  'stdev': 10.0,
  'type': 'normal'},
 {'event': 2.0,
  'initial': 90.0,
  'max': 95.0,
  'mean': 90.0,
  'min': 85.0,
  'parameter': 'Dip Direction',
  'stdev': 10.0,
  'type': 'normal'},
 {'event': 2.0,
  'initial': 1000.0,
  'max': 1100.0,
  'mean': 1000.0,
  'min': 900.0,
  'parameter': 'Slip',
  'stdev': 500.0,
  'type': 'normal'},
 {'event': 3.0,
  'initial': 60.0,
  'max': 65.0,
  'mean': 60.0,
  'min': 55.0,
  'parameter': 'Dip',
  'stdev': 10.0,
  'type': 'normal'},
 {'event': 3.0,
  'initial': 270.0,
  'max': 275.0,
  'mean': 270.0,
  'min': 265.0,
  'parameter': 'Dip Direction',
  'stdev': 10.0,
  'type': 'normal'},
 {'event': 3.0,
  'initial': 1000.0,
  'max': 1100.0,
  'mean': 1000.0,
  'min': 900.0,
  'parameter': 'Slip',
  'stdev': 500.0,
  'type': 'normal'}]

In [58]:
sa.plot_section(model_type = "base")



In [13]:
plt.plot(sa.get_drillhole_data(4000, 3500))
plt.plot(sa.get_model_lines())


Out[13]:
[<matplotlib.lines.Line2D at 0x10452bdd0>]

In [60]:
reload(pynoddy.history)
reload(pynoddy.experiment)
sa2 = pynoddy.experiment.Experiment(history = "two_faults_sensi.his")


 STRATIGRAPHY
 FAULT
 FAULT

In [61]:
sa2.write_history("test5.his")

In [62]:
nm = pynoddy.history.NoddyHistory(history = "two_faults_sensi.his")
# nm.determine_events()
nm.write_history("test6.his")


 STRATIGRAPHY
 FAULT
 FAULT

In [63]:
param_values[0]


Out[63]:
array([   57.19726562,    85.96679688,  1003.7109375 ,    61.76757812,
         267.80273438,  1081.4453125 ])

In [64]:
reload(pynoddy.history)
reload(pynoddy.experiment)
# Start again with the original model
his_filename = "two_faults_sensi.his"
sa = pynoddy.experiment.SensitivityAnalysis(history = his_filename)

sa.freeze()
# sa.change_event_params({3 : {'Slip' : 500.}})
sa.change_event_params({3 : {'Dip' : 15.}})


 STRATIGRAPHY
 FAULT
 FAULT
{3: {'Dip': 15.0}}

In [65]:
fig = plt.figure(figsize = (12,6))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
sa.plot_section(ax = ax1, colorbar = False, title = "")
sa.plot_section(ax = ax2, model_type = "base", colorbar = False, title = "")



In [66]:
sa.change_event_params({3 : {'Slip' : 100.}})
sa.plot_section()


{3: {'Slip': 100.0}}

In [14]:
# sa.add_sampling_line(5000, 3500, label = 'centre')
# sa.add_sampling_line(2500, 3500, label = 'left')
sa.add_sampling_line(7500, 3500, label = 'right')
# sa.add_sampling_line(4000, 3500, label = 'compare')
plt.plot(sa.get_model_lines(), 'k')
plt.plot(sa.get_model_lines(model_type = "base"), 'b')


Out[14]:
[<matplotlib.lines.Line2D at 0x10414b910>]

In [68]:
pwd


Out[68]:
u'/Users/flow/git/pynoddy/docs/notebooks'

In [68]: