Wrapping an existing model#

SALib performs sensitivity analysis for any model that can be expressed in the form of \(f(X) = Y\), where \(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 \(f\) is.

Typical model implementations take the form of \(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:

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:

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.

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 \(N*D\), where \(D\) is the number of model factors (dimensions) and \(N\) is the number of their combinations. The argument name is, by convention, denoted as 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 functools.partial() from the functools package to create wrappers can be useful.

In this example, the model (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 for loop as in the example below.

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:

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 .evaluate(wrapped_linear, nprocs=2) will fail with NameError: name 'np' is not defined.

The quick fix is to re-import the required packages within the model function itself:

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 drivendata/cookiecutter-data-science, be adopted to improve usability and reproducibility.

Here, model_definitions.py holds the model definitions:

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 analysis.py contains use of SALib:

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 __main__ check to ensure it is only run in the top-level environment.

if __name__ == "__main__":
    (
        sp.sample_sobol(2**6)
        .evaluate(wrapped_linear, nprocs=2)
        .analyze_sobol(nprocs=2)
    )