[ot-users] sobol with multi-output

Pamphile ROY roy at cerfacs.fr
Wed Jun 7 12:30:01 CEST 2017


Hi, 

Thanks for the hint. I will stick with Saltelli then. 

Pamphile 


De: "Julien Schueller | Phimeca" <schueller at phimeca.com> 
À: "roy" <roy at cerfacs.fr> 
Cc: "users" <users at openturns.org> 
Envoyé: Mercredi 7 Juin 2017 09:35:06 
Objet: RE: [ot-users] sobol with multi-output 






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 <roy at cerfacs.fr> 
Envoyé : mercredi 31 mai 2017 14:37 
À : Julien Schueller | Phimeca 
Cc : users at openturns.org 
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 






Le 31 mai 2017 à 14:07, Julien Schueller | Phimeca < schueller at phimeca.com > a écrit : 

Hello, 
Do you have a minimal non-working example ? 
j 


De : users-bounces at openturns.org < users-bounces at openturns.org > de la part de roy < roy at cerfacs.fr > 
Envoyé : mercredi 31 mai 2017 13:48 
À : users at openturns.org 
Objet : [ot-users] sobol with multi-output 
Hi everyone, 

I am using the SobolIndicesAlgorithm class using MartinezSensitivityAlgorithm and SaltelliSensitivityAlgorithm. 

Using the same function : f(x, y) -> R^400 it only works with Saltelli’s class. Using Martinez, I get zeros everywhere. 
I suspect a mis-implementation somewhere. 

Are you aware of this? 

Thanks in advance, 


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 





-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://openturns.org/pipermail/users/attachments/20170607/511db624/attachment.html>


More information about the users mailing list