[ot-users] sobol with multi-output

roy roy at cerfacs.fr
Wed May 31 14:37:21 CEST 2017


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.

WRN - 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
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 <mailto:users-bounces at openturns.org> <users-bounces at openturns.org <mailto:users-bounces at openturns.org>> de la part de roy <roy at cerfacs.fr <mailto:roy at cerfacs.fr>>
> Envoyé : mercredi 31 mai 2017 13:48
> À : users at openturns.org <mailto: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/20170531/bece0729/attachment.html>


More information about the users mailing list