Hi, Thanks for the hint. I will stick with Saltelli then.
Pamphile De: "Julien Schueller | Phimeca" <[email protected]> À: "roy" <[email protected]> Cc: "users" <[email protected]> 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 <[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 Le 31 mai 2017 à 14:07, Julien Schueller | Phimeca < [email protected] > a écrit : Hello, Do you have a minimal non-working example ? j De : [email protected] < [email protected] > de la part de roy < [email protected] > Envoyé : mercredi 31 mai 2017 13:48 À : [email protected] 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
_______________________________________________ OpenTURNS users mailing list [email protected] http://openturns.org/mailman/listinfo/users
