==========================
Wrapping an existing model
==========================
SALib performs sensitivity analysis for any model that can be expressed in the form of :math:`f(X) = Y`,
where :math:`X` is a matrix of inputs (often referred to as the model's factors).
The analysis methods are independent of the model and can be applied non-intrusively such
that it does not matter what :math:`f` is.
Typical model implementations take the form of :math:`f(a, b, c, ...) = Y`. In other words, each model
factor is supplied as a separate argument to the function. In such cases it is necessary to
write a wrapper to allow use with SALib. This is illustrated here with a simple linear function:
.. code:: python
def linear(a, b, x):
"""Return y = a + b + x"""
return a + b + x
As SALib expects a (numpy) matrix of factors, we simply "wrap" the function above like so:
.. code:: python
def wrapped_linear(X, func=linear):
"""g(X) = Y, where X := [a b x] and g(X) := f(X)"""
# We transpose to obtain each column (the model factors) as separate variables
a, b, x = X.T
# Then call the original model
return func(a, b, x)
.. note:: **Wrapped function is an argument**
Note here that the model being "wrapped" is also passed in as an argument.
This will be revisited further down below.
.. tip:: **Interfacing with external models/programs**
Here we showcase interacting with models written in Python.
If the model is an external program, this is where interfacing code
would be written.
A pragmatic approach could be to use `subprocess `_
to start the external program, then read in the results.
Constants, which SALib should not consider, can be expressed by defining default keyword arguments
(for flexibility) or otherwise defined within the wrapper function itself.
.. code:: python
def wrapped_linear_w_constant(X, a=10, func=linear):
"""f(X, a) = Y, where X := [b x] and a = 10"""
# We transpose to obtain each column as separate variables
b, x = X.T
# Then call the original model
return func(a, b, x)
Note that the first argument to any function provided to SALib is assumed to be
a numpy array of shape :math:`N*D`, where :math:`D` is the number of model
factors (dimensions) and :math:`N` is the number of their combinations. The
argument name is, by convention, denoted as :code:`X`. This is to maximize
compatibility with all methods provided in SALib as they expect the first
argument to hold the model factor values. Using :py:func:`functools.partial`
from the `functools` package to create wrappers can be useful.
In this example, the model (:code:`linear()`) can be used with both scalar
inputs or `numpy` arrays. In cases where `a`, `b` or `x` are a vector of
inputs, `numpy` will automatically vectorize the calculation. There are many
cases where the model is not (or cannot be easily) expressed in a vectorizable
form. In such cases, simply apply a :code:`for` loop as in the example below.
.. code:: python
import numpy as np
from SALib import ProblemSpec
def linear(a: float, b: float, x: float) -> float:
return a + b * x
def wrapped_linear(X: np.ndarray, func=linear) -> np.ndarray:
N, D = X.shape
results = np.empty(N)
for i in range(N):
a, b, x = X[i, :]
results[i] = func(a, b, x)
return results
sp = ProblemSpec({
'names': ['a', 'b', 'x'],
'bounds': [
[-1, 0],
[-1, 0],
[-1, 1],
],
})
(
sp.sample_sobol(2**6)
.evaluate(wrapped_linear)
.analyze_sobol()
)
sp.to_df()
# [ ST ST_conf
# a 0.173636 0.072142
# b 0.167933 0.059599
# x 0.654566 0.208328,
# S1 S1_conf
# a 0.182788 0.111548
# b 0.179003 0.145714
# x 0.664727 0.241977,
# S2 S2_conf
# (a, b) -0.022070 0.185510
# (a, x) -0.010781 0.186743
# (b, x) -0.014616 0.279925]
Use of the core SALib functions equivalent to the previous example are shown
below:
.. code:: python
problem = {
'names': ['a', 'b', 'x'],
'bounds': [
[-1, 0],
[-1, 0],
[-1, 1],
],
'num_vars': 3
}
X = saltelli.sample(problem, 64)
Y = np.empty(params.shape[0])
for i in range(params.shape[0]):
Y[i] = wrapped_linear(params[i, :])
res = sobol.analyze(problem, Y)
res.to_df()
# [ ST ST_conf
# a 0.165854 0.054096
# b 0.165854 0.053200
# x 0.665366 0.192756,
# S1 S1_conf
# a 0.167805 0.121550
# b 0.167805 0.125178
# x 0.665366 0.230872,
# S2 S2_conf
# (a, b) -2.775558e-17 0.180493
# (a, x) -3.902439e-03 0.202343
# (b, x) -3.902439e-03 0.232957]
Parallel evaluation and analysis
--------------------------------
Here we expand on some technical details that enable parallel evaluation and
analysis. We noted earlier that the model being "wrapped" is also passed in
as an argument. This is to facilitate parallel evaluation, as the arguments
to the wrapper are passed on to workers. The approach works be using Python's
`mutable default argument `_
behavior.
A further consideration is that imported modules/packages are not made
available to workers in cases where functions are defined in the same file
SALib is used in. Running the previous example with :code:`.evaluate(wrapped_linear, nprocs=2)`
will fail with :code:`NameError: name 'np' is not defined`.
The quick fix is to re-import the required packages within the model function
itself:
.. code:: python
def wrapped_linear(X: np.ndarray, func=linear) -> np.ndarray:
import numpy as np # re-import necessary packages
N, D = X.shape
results = np.empty(N)
for i in range(N):
a, b, x = X[i, :]
results[i] = func(a, b, x)
return results
This can, however, get unwieldy for complicated models. The recommended best
practice is to separate implementation (i.e., model definitions) from its use.
Simply moving the model functions into a separate file is enough for this
example, such that the project structure is something like:
::
project_directory
|-- model_definition.py
└── analysis.py
.. tip:: **Project structure**
The project structure shown above is for example purposes only.
It is highly recommended that a standardized directory structure,
such as https://github.com/drivendata/cookiecutter-data-science,
be adopted to improve usability and reproducibility.
Here, :code:`model_definitions.py` holds the model definitions:
.. code:: python
import numpy as np
def linear(a: float, b: float, x: float) -> float:
return a + b * x
def wrapped_linear(X: np.ndarray, func=linear) -> np.ndarray:
N, D = X.shape
results = np.empty(N)
for i in range(N):
a, b, x = X[i, :]
results[i] = func(a, b, x)
return results
and :code:`analysis.py` contains use of SALib:
.. code:: python
from model_definition import wrapped_linear
sp = ProblemSpec({
'names': ['a', 'b', 'x'],
'bounds': [
[-1, 0],
[-1, 0],
[-1, 1],
],
})
(
sp.sample_sobol(2**6)
.evaluate(wrapped_linear, nprocs=2)
.analyze_sobol(nprocs=2)
)
.. note:: **Multi-processing**
Some interactive Python consoles, including earlier versions of
IPython, may appear to hang on Windows when utilizing parallel
evaluation and analysis. In such cases, the recommended workaround
is to wrap use of SALib with a :code:`__main__` check to ensure
it is only run in the top-level environment.
.. code:: python
if __name__ == "__main__":
(
sp.sample_sobol(2**6)
.evaluate(wrapped_linear, nprocs=2)
.analyze_sobol(nprocs=2)
)