|
Hi, I looked more closely : getFirstOrderIndices throws an is caught at the same time :( Martinez does some normalization, an a stddev component is zero at the 800th model evaluation. I can work around this by adding a random noise: return self.zref + h + np.array(ot.Normal([0.0]*400, [0.01]*400, ot.IdentityMatrix(400)).getRealization()) I'll if something cleaner is possible.
j
De : roy <[email protected]>
Envoyé : mercredi 31 mai 2017 14:37 À : Julien Schueller | Phimeca Cc : [email protected] Objet : Re: [ot-users] sobol with multi-output Hi Julien,
Thanks for your prompt reply. Here is the MNWE.
Sorry for the delay I had to «crack» my code.
##########################################################
import numpy as np
import openturns as ot
import otwrapy as otw
import logging
class Wrapper(ot.OpenTURNSPythonFunction):
"""Wrap predictor with OpenTURNS."""
def __init__(self, f, p_len, output_len, block=False):
"""Initialize the wrapper."""
super(Wrapper, self).__init__(p_len, output_len)
self.surrogate = f
self._exec = self.func
def func(self, coords):
"""Evaluate the surrogate at a given point."""
return self.surrogate(coords)
class Channel_Flow(object):
r"""Channel Flow class.
.. math:: \frac{dh}{ds}=\mathcal{F}(h)=I\frac{1-(h/h_n)^{-10/3}}{1-(h/h_c)^{-3}}
with :math:`h_c=\left(\frac{q^2}{g}\right)^{1/3}, h_n=\left(\frac{q^2}{IK_s^2}\right)^{3/10}`.
"""
logger = logging.getLogger(__name__)
def __init__(self, dx=100., length=40000., width=500.):
"""Initialize the geometrical configuration.
:param float dx: discretization
:param float length: canal length
:param float width: canal width
"""
self.w = width
self.I = 5e-4
self.g = 9.8
self.dx = dx
self.length = length
self.x = np.arange(self.dx, self.length + 1, self.dx)
self.d_out = len(self.x)
self.d_in = 2
self.dl = int(self.length // self.dx)
self.hinit = 10.
self.zref = - self.x * self.I
# Sensitivity
self.s_first = np.array([0.92925829, 0.05243018])
self.s_second = np.array([[0., 0.01405351], [0.01405351, 0.]])
self.s_total = np.array([0.93746788, 0.05887997])
self.logger.info("Using function Channel Flow with: dx={}, length={}, "
"width={}".format(dx, length, width))
def __call__(self, x):
"""Call function.
:param list x: inputs [Ks, Q]
:return: Water height along the channel
:rtype: np.array 1D
"""
ks, q = x
hc = np.power((q ** 2) / (self.g * self.w ** 2), 1. / 3.)
hn = np.power((q ** 2) / (self.I * self.w ** 2 * ks ** 2), 3. / 10.)
h = self.hinit * np.ones(self.dl)
for i in range(2, self.dl + 1):
h[self.dl - i] = h[self.dl - i + 1] - self.dx * self.I\
* ((1 - np.power(h[self.dl - i + 1] / hn, -10. / 3.))
/ (1 - np.power(h[self.dl - i + 1] / hc, -3.)))
return self.zref + h
f = Channel_Flow()
p_len = 2
output_len = 400
n_cpus = 2
points_sample = 2000
distribution = ot.ComposedDistribution([ot.Uniform(15., 60.),
ot.Normal(4035., 400.)],
ot.IndependentCopula(2))
wrapper = Wrapper(f, p_len, output_len)
model = otw.Parallelizer(wrapper, backend='pathos', n_cpus=n_cpus)
input_design = ot.SobolIndicesAlgorithmImplementation.Generate(
distribution, points_sample, True)
output_design = model(input_design)
ot.ResourceMap.SetAsBool('MartinezSensitivityAlgorithm-UseAsmpytoticInterval', True)
###### Martinez, Saltelli ######
# sobol = ot.SaltelliSensitivityAlgorithm(input_design, output_design, points_sample)
sobol = ot.MartinezSensitivityAlgorithm(input_design, output_design, points_sample)
###### Martinez, Saltelli ######
indices = [[], []]
for i in range(output_len):
try:
indices[0].append(np.array(sobol.getFirstOrderIndices(i)))
except TypeError:
indices[0].append(np.zeros(p_len))
try:
indices[1].append(np.array(sobol.getTotalOrderIndices(i)))
except TypeError:
indices[1].append(np.zeros(p_len))
print("First order: {}".format(indices[0]))
print("Total: {}".format(indices[1]))
##########################################################
The try block is here to overcome the following. It is due to some limit conditions. The value being fix, the calcul can be done for this index.
[34m[1mWRN - The estimated total order Sobol index (0) is lesser than first order index . You may increase the sampling size. HERE we have: S_0=0.996425, ST_0=0.95792[0m
Traceback (most recent call last):
File "/Users/roy/Desktop/test.py", line 103, in <module>
indices[0].append(np.array(sobol.getFirstOrderIndices(i)))
File "/Users/roy/Applications/miniconda3/envs/batman3/lib/python3.6/site-packages/openturns/uncertainty.py", line 1442, in getFirstOrderIndices
return _uncertainty.SobolIndicesAlgorithmImplementation_getFirstOrderIndices(self, marginalIndex)
TypeError: InvalidArgumentException : Error: cannot divide by 0.
If you run this, the output is a list of zeros but if you run it with Saltelli, I get the expected output.
Pamphile ROY
Chercheur doctorant en Quantification d’Incertitudes CERFACS - Toulouse (31) - France +33 (0) 5 61 19 31 57 +33 (0) 7 86 43 24 22
|
_______________________________________________ OpenTURNS users mailing list [email protected] http://openturns.org/mailman/listinfo/users
