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})
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)
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()
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]:
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]:
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]:
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")
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:
In addition, it is possible to define PDF type and parameters. For now, the following settings are supported:
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")
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]:
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]:
In [60]:
reload(pynoddy.history)
reload(pynoddy.experiment)
sa2 = pynoddy.experiment.Experiment(history = "two_faults_sensi.his")
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")
In [63]:
param_values[0]
Out[63]:
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.}})
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()
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]:
In [68]:
pwd
Out[68]:
In [68]: