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.

[1]:

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.

[2]:

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:

[3]:

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.

[4]:

borehole_domain.sample(1)

[4]:

array([6.66855202e-02, 1.10622146e+04, 9.07114369e+04, 1.10078797e+03,
       1.04110193e+02, 7.41295866e+02, 1.53865350e+03, 1.02057674e+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.

[5]:

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.

[6]:

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

Active Subspace¶

[7]:

X = borehole.domain.sample(200)
grads = borehole.grad(X)

[8]:

active = psdr.ActiveSubspace()
active.fit(grads)

[9]:

fX = borehole(X)
active.shadow_plot(X, fX, dim = 1)

[9]:

<AxesSubplot:xlabel='active coordinate $\\mathbf{u}^\\top \\mathbf{x}$', ylabel='$f(\\mathbf{x})$'>

[10]:

active.shadow_plot(X, fX, dim = 2)

[10]:

<AxesSubplot:xlabel='active coordinate 1 $\\mathbf{u}_1^\\top \\mathbf{x}$', ylabel='active coordinate 2 $\\mathbf{u}_2^\\top \\mathbf{x}$'>

Constructing a Ridge Approximation¶

[11]:

X = borehole.domain.sample(200)
fX = borehole(X)

[12]:

pra = psdr.PolynomialRidgeApproximation(degree = 5, subspace_dimension = 1, norm = 2, bound = 'upper')
pra.fit(X, fX)

/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "
/home/docs/checkouts/readthedocs.org/user_builds/psdr/envs/latest/lib/python3.7/site-packages/cvxpy/problems/problem.py:1061: UserWarning: Solution may be inaccurate. Try another solver, adjusting the solver settings, or solve with verbose=True for more information.
  "Solution may be inaccurate. Try another solver, "

[13]:

pra.shadow_plot(X, fX);

[13]:

<AxesSubplot:xlabel='active coordinate $\\mathbf{u}^\\top \\mathbf{x}$', ylabel='$f(\\mathbf{x})$'>

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

[14]:

pra2 = psdr.PolynomialRidgeApproximation(degree = 5, subspace_dimension = 2)
pra2.fit(X, fX)

[15]:

pra2.U

[15]:

array([[ 0.9184318 ,  0.39405712],
       [ 0.00254579, -0.0480603 ],
       [ 0.00099252, -0.05184537],
       [ 0.22154035, -0.51814463],
       [ 0.00712601, -0.07182107],
       [-0.21655779,  0.50038448],
       [-0.22177025,  0.49907056],
       [ 0.10613828, -0.25813383]])

[16]:

pra2.shadow_plot(X, fX)

[16]:

<AxesSubplot:xlabel='active coordinate 1 $\\mathbf{u}_1^\\top \\mathbf{x}$', ylabel='active coordinate 2 $\\mathbf{u}_2^\\top \\mathbf{x}$'>

[ ]: