Defining a Custom Function

This is a short tutorial on using the Parameter Space Dimension Reduction library on the (in)famous borehole test function. Although this function is provided as a demo (with analytic gradients), here we provide a complete example setting up a function in this library.

Standard imports

To begin, we import numpy and psdr.

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import psdr

Problem Setup

To define a problem for use with this library, there are two basic steps. First we need to define the domain on which we will be working. As engineering applications frequently have poorly scaled variables with different units, defining a domain allows us to correct this effect when evaluating the importance of a particular set of parameters. Second, we need to define the function, posed on this domain, that we are interested in studying. Depending on the technique we use, we need to either have access to the function itself or to samples of the function. Here we define the function and then sample it as appropreate for each technique.

Defining the domain

The borehole function has eight variables:

Parameter and Domain Description
\(r_w \in [0.05, 0.15]\) radius of borehole (m)
\(r \in [100, 50\times 10^3]\) radius of influence (m)
\(T_u \in [63700, 115600]\) transmissivity of upper aquifer (\(m^2\)/yr)
\(H_u \in [990, 1110]\) potentiometric head of upper aquifer (m)
$T_:nbsphinx-math:ell `:nbsphinx-math:in [63.1, 116] `$ transmissivity of lower aquifer (\(m^2\)/yr)
\(H_\ell \in [700, 820]\) potentiometric head of lower aquifer (m)
\(L \in [1120, 1680]\) length of borehole (m)
\(K_w \in [9855, 12045]\) hydraulic conductivity of borehole (m/yr)

To begin, we first setup a domain for each variable.

rw_domain = psdr.BoxDomain(0.05, 0.15)
r_domain = psdr.BoxDomain(100, 50e3)
Tu_domain = psdr.BoxDomain(63700, 115600)
Hu_domain = psdr.BoxDomain(990, 1110)
Tell_domain = psdr.BoxDomain(63.1, 116)
Hell_domain = psdr.BoxDomain(700, 820)
L_domain = psdr.BoxDomain(1120, 1680)
Kw_domain = psdr.BoxDomain(9855, 12045)

Then, we combine these domains for each variable into a domain for the entire borehole function:

borehole_domain = psdr.TensorProductDomain([
    rw_domain, r_domain, Tu_domain, Hu_domain,
    Tell_domain, Hell_domain, L_domain, Kw_domain])

This class instance provides many functions to interact with the domain. For example, we can sample a random point in this domain.

array([1.26086551e-01, 2.13484764e+04, 9.81072917e+04, 1.06116661e+03,
       8.43097799e+01, 8.12394991e+02, 1.48176131e+03, 1.00089628e+04])

Defining the Function

The borehole function is:

\[f(r_w, r, T_u, H_u, T_\ell, H_\ell, L, K_w) := \frac{ 2\pi T_u (H_u - H_\ell)}{\log(r/r_w) ( 1 + \frac{2L T_u}{\log(r/r_w) r_w^2 K_w} + \frac{T_u}{T_\ell} )}.\]

So we simply define this function that takes an array-like input and returns the output.

def borehole_func(x):
    r_w = x[0]
    r = x[1]
    T_u = x[2]
    H_u = x[3]
    T_l = x[4]
    H_l = x[5]
    L = x[6]
    K_w = x[7]
    return 2*np.pi*T_u*(H_u - H_l)/(np.log(r/r_w)*(1 + 2*L*T_u/(np.log(r/r_w)*r_w**2*K_w) + T_u/T_l))

Wrapping the Function

For convience, we can wrap a function in the Function environment. This does several important things. First, it automatically converts function on the application domain to a function on the normalized domain. Further, as illustrated here, it can automatically estimate gradients using a finite difference approximation. However, if gradients are availble, these can also be provided.

borehole = psdr.Function(borehole_func, borehole_domain, fd_grad = True)

Active Subspace

X = borehole.domain.sample(200)
grads = borehole.grad(X)
active = psdr.ActiveSubspace()
fX = borehole(X)
active.shadow_plot(X, fX, dim = 1)
<matplotlib.axes._subplots.AxesSubplot at 0x7fd1e1bdc048>
active.shadow_plot(X, fX, dim = 2)
<matplotlib.axes._subplots.AxesSubplot at 0x7fd1df78b3c8>

Constructing a Ridge Approximation

X = borehole.domain.sample(200)
fX = borehole(X)
pra = psdr.PolynomialRidgeApproximation(degree = 5, subspace_dimension = 1, norm = 2, bound = 'upper'), fX)
pra.shadow_plot(X, fX);

Of course, we can also construct ridge approximations with higher dimensional subspaces.

pra2 = psdr.PolynomialRidgeApproximation(degree = 5, subspace_dimension = 2), fX)
array([[ 9.14198112e-01,  4.00956972e-01],
       [ 3.40707457e-04, -3.28473218e-02],
       [ 2.65737966e-04,  4.08062840e-02],
       [ 2.24130490e-01, -5.34285853e-01],
       [ 5.61918566e-03, -5.59284153e-02],
       [-2.23417157e-01,  5.70980932e-01],
       [-2.26785340e-01,  4.06213378e-01],
       [ 1.12377742e-01, -2.38477492e-01]])
pra2.shadow_plot(X, fX)
<matplotlib.axes._subplots.AxesSubplot at 0x7fd1df987eb8>
[ ]: