Sensivity analysis with python using SALib

In [62]:

    

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


    

In [5]:

    

# STATE THE PROBLEM DICTIONNARY
# what will be varying (=inputs) ? in what bounds ?

problem = {
    'num_vars': 3,
    'names': ['x1', 'x2', 'x3'],
    'bounds': [[-3.14159265359, 3.14159265359],
               [-3.14159265359, 3.14159265359],
               [-3.14159265359, 3.14159265359]]
}


    

In [6]:

    

# say we want to draw 150 samples

num_samples = 150


    

In [7]:

    

# we draw a Latin Hypercube Sampling (LHS) that is fitted for an RBD FAST analysis
# (other sampling metods available in the library though)

from SALib.sample.latin import sample

all_samples = sample(problem, num_samples)


    

In [12]:

    

# this is where you define your own model, procedure, experiment...
from SALib.test_functions import Ishigami
import numpy as np

def run_model(x1, x2, x3):
    """
    COPY HERE YOUR OWN CODE
    
    the function takes as input 1 sample
    returns 1 or more output
    """
    
    # Delete from HERE =======================================
    
    # As an example, we'll look at the famous Ishigami function
    # (A and B are specific parameters for the Ishigami function)
    A = 7
    B = 0.1
    y = (np.sin(x1) +
         A * np.sin(x2)**2 + 
         B * x3 ** 4 * np.sin(x1))

    # ========= TO HERE and replace with your own piece of code
    
    return y


    

In [19]:

    

# run your model, procedure, experiment for each sample of your sampling

# unpack all_samples into 3 vectors x1, x2, x3
x1, x2, x3 = all_samples.T

# run the model, all samples at the same time
ishigami_results = run_model(x1, x2, x3)


    

In [17]:

    

from SALib.analyze import rbd_fast


    

In [20]:

    

# Let us look at the analyze method

rbd_fast.analyze(problem=problem, Y=ishigami_results, X=all_samples)


    

    
Out[20]:





{'S1': [0.3248982740323958, 0.5273048950859327, -0.03802652771291007],
 'names': ['x1', 'x2', 'x3']}

In [21]:

    

# storing the first order indices of the analyze method
si1 = rbd_fast.analyze(problem=problem, Y=ishigami_results, X=all_samples)['S1']


    

In [23]:

    

# make nice plots with the indices (looks good on your presentations)
# do not use the plotting tools of SALib, they are made for the method of Morris ...

fig, ax = plt.subplots()
fig.set_size_inches(18, 5)

ax.tick_params(labelsize=18)

# ===== X-AXIS =====
ax.set_xticks(np.arange(problem['num_vars']))
ax.set_xticklabels(problem['names'])
ax.set_xlim(xmin=-0.5, xmax=problem['num_vars'] - 0.5)

# ===== Y-AXIS =====
ax.set_ylabel('RBD-FAST\nSensitivity indices\n', fontsize=25)

# ===== BARS REPRESENTING THE SENSITIVITY INDICES =====
ax.bar(np.arange(problem['num_vars']), si1,
       color='DarkSeaGreen');

#in striped grey : not significant indices 
ax.fill_between(x=[-0.5, 5.5], y1=-0.1, y2=0.1, color='grey', alpha=0.2, hatch='//', edgecolor='white');


    

In [29]:

    

# take a closer look to undestand interactions (looks even better on your presentations)
# this part can be done without analyzing with SALib, just with the ouput from the samples

fig, ax = plt.subplots()
fig.set_size_inches(8, 7)

ax.tick_params(labelsize=15)
ax.set_title('Output of the model studied (Ishigami)\n'
             'according to the value of the 2 most influent parameters',
             fontsize=20)

# ===== SCATTER =====
size = np.ones(num_samples) * 75

sc = ax.scatter(x1, x2,
                c=ishigami_results,
                s=size,
                vmin=ishigami_results.min(),
                vmax=ishigami_results.max(),
                cmap='seismic',
                edgecolor=None)


# ===== X-AXIS =====
ax.set_xlim(xmin=problem['bounds'][0][0], xmax=problem['bounds'][0][1])
ax.set_xlabel('Parameter 1', fontsize=20)


# ===== Y-AXIS =====
ax.set_ylim(ymin=problem['bounds'][1][0], ymax=problem['bounds'][1][1])
ax.set_ylabel('Parameter 2', fontsize=20)

# ===== COLORBAR =====
ticks = np.arange(ishigami_results.min(), ishigami_results.max(), 5)
cb = plt.colorbar(sc, ticks=ticks);
cb.ax.set_yticklabels([str(round(i,1)) for i in ticks])

fig.tight_layout();


    

In [53]:

    

def conv_study(n, Y, X):
    # take n samples among the num_samples, without replacement
    subset = np.random.choice(num_samples, size=n, replace=False)
    return rbd_fast.analyze(problem=problem,
                            Y=Y[subset],
                            X=X[subset])['S1']

all_indices = np.array([conv_study(n=n, Y=ishigami_result, X=all_samples)
                        for n in np.arange(50, num_samples + 1, 5)])


    

In [54]:

    

# convergence check
fig, ax = plt.subplots()
fig.set_size_inches(15,8)

ax.set_title('Convergence of the sensitivity indices', fontsize=20)

ax.tick_params(labelsize=16)

ax.set_xlim(xmin=0, xmax=(num_samples - 50)//10)
ax.set_xticks(np.arange(0,(num_samples - 50)//10 +1,5))
ax.set_xticklabels([str(i) for i in range(50,num_samples+1,50)])

ax.set_ylim(ymin=-0.15, ymax=1)


for p,param in enumerate(problem['names']):
    ax.plot(all_indices[:,p], linewidth=3, label=param)

ax.fill_between(x=[0,(num_samples - 50)//10], y1=-0.15, y2=0.1, color='grey', alpha=0.2, hatch='//', edgecolor='white',
              label='REMINDER : not significant')
    
ax.legend(fontsize=20, ncol=4);


    

In [59]:

    

def bootstrap(problem, Y, X):
    """
    Calculate confidence intervals of rbd-fast indices
    1000 draws
    returns 95% confidence intervals of the 1000 indices
    
    problem : dictionnary as SALib uses it
    X : SA input(s)
    Y : SA output(s)
    """
    all_indices = []
    
    for i in range(1000):
        X_new = np.zeros(X.shape)
        Y_new = np.zeros(Y.shape).flatten()
        # draw with replacement
        tirage_indices = np.random.randint(0, high=Y.shape[0], size = Y.shape[0])
        
        for j, index in enumerate(tirage_indices):
            X_new[j,:] = X[index, :]
            Y_new[j] = Y[index]
            
        all_indices.append(rbd_fast.analyze(problem=problem, Y=Y_new, X=X_new)['S1'])
    
    means = np.array([i.mean() for i in np.array(all_indices).T])
    stds = np.array([i.std() for i in np.array(all_indices).T])
    return np.array([means - 2 * stds, means + 2 * stds])


    

In [60]:

    

# Get bootstrap confidence intervals for each index

bootstrap_conf = bootstrap(problem=problem, Y=ishigami_results, X=all_samples)


    

In [61]:

    

# make nice plots with the indices (looks good on your presentations)
# do not use the plotting tools of SALib, they are made for the method of Morris ...

fig, ax = plt.subplots()
fig.set_size_inches(8,4)

ax.tick_params(labelsize=18)

# ===== X-AXIS =====
ax.set_xticks(np.arange(problem['num_vars']))
ax.set_xticklabels(problem['names'])
ax.set_xlim(xmin=-0.5, xmax=problem['num_vars'] - 0.5)

# ===== Y-AXIS =====
ax.set_ylabel('RBD-FAST\nSensitivity indices\n', fontsize=25)

# ===== BARS REPRESENTING THE SENSITIVITY INDICES =====
ax.bar(np.arange(problem['num_vars']), si1,
         color='DarkSeaGreen');

# ===== LINES REPRESENTING THE BOOTSTRAP "CONFIDENCE INTERVALS" =====
for j in range(problem['num_vars']):
    ax.plot([j, j], [bootstrap_conf[0, j], bootstrap_conf[1, j]],
              'k')

#ax.fill_between(x=[-0.5, 5.5], y1=-0.1, y2=0.1, color='grey', alpha=0.2, hatch='//', edgecolor='white');