Example 1 (using Chaospy (3.0.5))

$ x :nbsphinx-math:`in [-2, 2] `$

$ y = x - 0.5x ^ 3 + 0.1x ^ 4 $

[1]:
# true model
def model(x):
    return x - 0.5 * x ** 3 + 0.1 * x ** 4
[2]:
import numpy as np
# set the domain of interest
pdom = np . empty ((2 , 1) )
pdom[0 , 0] = -2.0
pdom[1 , 0] = 2.0
[3]:
# load the UQ interface
from surrogate.chaospy_model import UQChaospy as UQ

# setup the UQ analysis ( proj : projection method )
uq = UQ(pdom, order=4, method='proj')

# create the training data
xtrain = uq.generate_quadrature()
ytrain = model(xtrain)

# train the surrogate model
uq.fit(xtrain, ytrain )

# predict the model response on validation points
xval = np.random.uniform(pdom [0 , 0], pdom[1 , 0], (5 , 1))
yval = uq.predict(xval)
---------------------------------------------------------------------------
InvalidParameterError                     Traceback (most recent call last)
Cell In[3], line 8
      5 uq = UQ(pdom, order=4, method='proj')
      7 # create the training data
----> 8 xtrain = uq.generate_quadrature()
      9 ytrain = model(xtrain)
     11 # train the surrogate model

File ~/pyOPALTools/surrogate/chaospy_model.py:184, in UQChaospy.generate_quadrature(self, rule)
    178 self._hasNodes = True
    180 self._nodes, self._weights = cp.generate_quadrature(self._order,
    181                                                     self._distribution,
    182                                                     rule=rule)
--> 184 xx = self._unscale(self._nodes[0, :].T, self._pdom[:, 0].T)
    185 for i in range(1, self._npar):
    186     xx = np.hstack((xx, self._unscale(self._nodes[i, :].T, self._pdom[:, i])))

File ~/pyOPALTools/surrogate/chaospy_model.py:297, in UQChaospy._unscale(self, x, pdom)
    296 def _unscale(self, x, pdom):
--> 297     LU_transform = MinMaxScaler(feature_range=(pdom)).fit([[-1], [1]])
    298     return LU_transform.transform(x.reshape(1, -1)).T

File /usr/local/lib/python3.8/dist-packages/sklearn/preprocessing/_data.py:435, in MinMaxScaler.fit(self, X, y)
    433 # Reset internal state before fitting
    434 self._reset()
--> 435 return self.partial_fit(X, y)

File /usr/local/lib/python3.8/dist-packages/sklearn/base.py:1145, in _fit_context.<locals>.decorator.<locals>.wrapper(estimator, *args, **kwargs)
   1140 partial_fit_and_fitted = (
   1141     fit_method.__name__ == "partial_fit" and _is_fitted(estimator)
   1142 )
   1144 if not global_skip_validation and not partial_fit_and_fitted:
-> 1145     estimator._validate_params()
   1147 with config_context(
   1148     skip_parameter_validation=(
   1149         prefer_skip_nested_validation or global_skip_validation
   1150     )
   1151 ):
   1152     return fit_method(estimator, *args, **kwargs)

File /usr/local/lib/python3.8/dist-packages/sklearn/base.py:638, in BaseEstimator._validate_params(self)
    630 def _validate_params(self):
    631     """Validate types and values of constructor parameters
    632
    633     The expected type and values must be defined in the `_parameter_constraints`
   (...)
    636     accepted constraints.
    637     """
--> 638     validate_parameter_constraints(
    639         self._parameter_constraints,
    640         self.get_params(deep=False),
    641         caller_name=self.__class__.__name__,
    642     )

File /usr/local/lib/python3.8/dist-packages/sklearn/utils/_param_validation.py:95, in validate_parameter_constraints(parameter_constraints, params, caller_name)
     89 else:
     90     constraints_str = (
     91         f"{', '.join([str(c) for c in constraints[:-1]])} or"
     92         f" {constraints[-1]}"
     93     )
---> 95 raise InvalidParameterError(
     96     f"The {param_name!r} parameter of {caller_name} must be"
     97     f" {constraints_str}. Got {param_val!r} instead."
     98 )

InvalidParameterError: The 'feature_range' parameter of MinMaxScaler must be an instance of 'tuple'. Got array([-2.,  2.]) instead.
[4]:
import matplotlib.pyplot as plt

plt.figure(figsize=(12, 6), dpi=150)

# evaluate true model
x = np.linspace(pdom[0, 0], pdom[1, 0], 100)
plt.plot(x, model(x), label='true model')

# plot training points
plt.plot(xtrain, ytrain, 'ro', label='training points')

# plot validation points
plt.plot(xval, yval, 'bo', label='predicted points')

plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid(linestyle='dashed')
plt.legend()
plt.show()
plt.close()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[4], line 10
      7 plt.plot(x, model(x), label='true model')
      9 # plot training points
---> 10 plt.plot(xtrain, ytrain, 'ro', label='training points')
     12 # plot validation points
     13 plt.plot(xval, yval, 'bo', label='predicted points')

NameError: name 'xtrain' is not defined
_images/uq_example_4_1.svg

Example: Ishigami function (using Chaospy (3.0.5))

$ m(x_1, x_2, x_3) = \sin`(x_1) + a :nbsphinx-math:cdot`:nbsphinx-math:sin`^2(x_2) + b :nbsphinx-math:cdot x_3^4:nbsphinx-math:sin`(x_1) $

with \(x_i\in [-\pi, \pi]\)

[5]:
def ishigami(x1, x2, x3, a, b):
    import numpy as np
    return np.sin(x1) + a * np.sin(x2) ** 2 + b * x3 ** 4 * np.sin(x1)
[6]:
import numpy as np
l = -np.pi
u =  np.pi

order = 9

a = 7
b = 0.05

nsamples = 440

np.random.seed(42)
x1 = np.random.uniform(l, u, nsamples)

np.random.seed(47)
x2 = np.random.uniform(l, u, nsamples)

np.random.seed(53)
x3 = np.random.uniform(l, u, nsamples)

xtrain = np.vstack((x1, x2, x3)).T

y_true = ishigami(xtrain[:, 0], xtrain[:, 1], xtrain[:, 2], a, b)
[7]:
# Analytical solution
D1  = b * np.pi ** 4 / 5.0 + b ** 2 * np.pi ** 8 / 50.0 + 0.5
D2  = a ** 2 / 8.0
D13 = (1.0 / 18 - 1.0 / 50.0) * b ** 2 * np.pi ** 8
D   = D1 + D2 + D13

print ( D1 / D )
print ( D2 / D )
print ( 0 )
print ( (D1 + 0 + D13) / D )
print ( (D2 + 0 + 0) / D )
print ( (0  + 0 + D13) / D )
0.21851856442701334
0.686894643648693
0
0.3131053563513071
0.686894643648693
0.09458679192429374
[8]:
from chaospy_model import UQChaospy
import numpy as np

pdom = np.empty((2, 3))

for i in range(3):
    pdom[0, i] = l
    pdom[1, i] = u

ch = UQChaospy(pdom, order)

ch.fit(xtrain, y_true)

ytrain = ch.predict(xtrain)
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[8], line 1
----> 1 from chaospy_model import UQChaospy
      2 import numpy as np
      4 pdom = np.empty((2, 3))

ModuleNotFoundError: No module named 'chaospy_model'
[9]:
sens_m = ch.main_sensitivity()
sens_t = ch.total_sensitivity()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[9], line 1
----> 1 sens_m = ch.main_sensitivity()
      2 sens_t = ch.total_sensitivity()

NameError: name 'ch' is not defined
[10]:
print ( sens_m )
print ( sens_t )
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[10], line 1
----> 1 print ( sens_m )
      2 print ( sens_t )

NameError: name 'sens_m' is not defined
[11]:
from bootstrap.bootstrap import bootstrap_sobol

s_m, s_t = bootstrap_sobol(xtrain, y_true, ch, n_boot=20, seed=42)

print ( s_m )
print ( s_t )
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[11], line 1
----> 1 from bootstrap.bootstrap import bootstrap_sobol
      3 s_m, s_t = bootstrap_sobol(xtrain, y_true, ch, n_boot=20, seed=42)
      5 print ( s_m )

ModuleNotFoundError: No module named 'bootstrap'