.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_reliability_sensitivity/reliability/plot_nais.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_reliability_sensitivity_reliability_plot_nais.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_reliability_sensitivity_reliability_plot_nais.py:


Non parametric Adaptive Importance Sampling (NAIS)
==================================================

.. GENERATED FROM PYTHON SOURCE LINES 6-33

The objective is to evaluate a probability from the Non parametric Adaptive Importance Sampling (NAIS) technique.

We consider the four-branch function :math:`g : \mathbb{R}^2 \rightarrow \mathbb{R}` defined by:

.. math::
  \begin{align*}
  g(\vect{X}) = \min \begin{pmatrix}5+0.1(x_1-x_2)^2-\frac{(x_1+x_2)}{\sqrt{2}}\\
  5+0.1(x_1-x_2)^2+\frac{(x_1+x_2)}{\sqrt{2}}\\
  (x_1-x_2)+ \frac{9}{\sqrt{2}}\\
  (x_2-x_1)+ \frac{9}{\sqrt{2}}
  \end{pmatrix}
  \end{align*}

and the input random vector :math:`\vect{X} = (X_1, X_2)` which follows the standard 2-dimensional Normal distribution:

.. math::
  \begin{align*}
  \vect{X} \sim  \mathcal{N}(\mu = [0, 0], \sigma = [1,1], corr = \mat{I}_2)
  \end{align*}

We want to evaluate the probability:

.. math::
  \begin{align*}
  p = \mathbb{P} ( g(\vect{X}) \leq 0 )
  \end{align*}


.. GENERATED FROM PYTHON SOURCE LINES 36-37

First, import the python modules:

.. GENERATED FROM PYTHON SOURCE LINES 39-43

.. code-block:: Python

    import openturns as ot
    from openturns.viewer import View
    import math








.. GENERATED FROM PYTHON SOURCE LINES 44-46

Create the probabilistic model :math:`Y = g(\vect{X})`
------------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 48-49

Create the input random vector :math:`\vect{X}`:

.. GENERATED FROM PYTHON SOURCE LINES 51-53

.. code-block:: Python

    X = ot.RandomVector(ot.Normal(2))








.. GENERATED FROM PYTHON SOURCE LINES 54-55

Create the function :math:`g` from a :class:`~openturns.PythonFunction`:

.. GENERATED FROM PYTHON SOURCE LINES 57-73

.. code-block:: Python



    def fourBranch(x):
        x1 = x[0]
        x2 = x[1]

        g1 = 5 + 0.1 * (x1 - x2) ** 2 - (x1 + x2) / math.sqrt(2)
        g2 = 5 + 0.1 * (x1 - x2) ** 2 + (x1 + x2) / math.sqrt(2)
        g3 = (x1 - x2) + 9 / math.sqrt(2)
        g4 = (x2 - x1) + 9 / math.sqrt(2)

        return [min((g1, g2, g3, g4))]


    g = ot.PythonFunction(2, 1, fourBranch)








.. GENERATED FROM PYTHON SOURCE LINES 74-75

Draw the function :math:`g` to help to understand the shape of the limit state function:

.. GENERATED FROM PYTHON SOURCE LINES 77-83

.. code-block:: Python

    graph = ot.Graph("Four Branch function", "x1", "x2", True, "upper right")
    drawfunction = g.draw([-8] * 2, [8] * 2, [100] * 2)
    graph.add(drawfunction)
    view = View(graph)





.. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_nais_001.png
   :alt: Four Branch function
   :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_nais_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 84-85

In order to be able to get the NAIS samples used in the algorithm, it is necessary to transform the :class:`~openturns.PythonFunction` into a :class:`~openturns.MemoizeFunction`:

.. GENERATED FROM PYTHON SOURCE LINES 87-89

.. code-block:: Python

    g = ot.MemoizeFunction(g)








.. GENERATED FROM PYTHON SOURCE LINES 90-91

Create the output random vector :math:`Y = g(\vect{X})`:

.. GENERATED FROM PYTHON SOURCE LINES 93-95

.. code-block:: Python

    Y = ot.CompositeRandomVector(g, X)








.. GENERATED FROM PYTHON SOURCE LINES 96-98

Create the event :math:`\{ Y = g(\vect{X}) \leq 0 \}`
-----------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 100-103

.. code-block:: Python

    threshold = 0.0
    myEvent = ot.ThresholdEvent(Y, ot.Less(), threshold)








.. GENERATED FROM PYTHON SOURCE LINES 104-106

Evaluate the probability with the NAIS technique
------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 108-111

.. code-block:: Python

    quantileLevel = 0.1
    algo = ot.NAIS(myEvent, quantileLevel)








.. GENERATED FROM PYTHON SOURCE LINES 112-113

In order to get all the inputs and outputs that realize the event, you have to mention it now:

.. GENERATED FROM PYTHON SOURCE LINES 115-117

.. code-block:: Python

    algo.setKeepSample(True)








.. GENERATED FROM PYTHON SOURCE LINES 118-119

Now you can run the algorithm.

.. GENERATED FROM PYTHON SOURCE LINES 121-127

.. code-block:: Python

    algo.run()
    result = algo.getResult()
    proba = result.getProbabilityEstimate()
    print("Proba NAIS = ", proba)
    print("Current coefficient of variation = ", result.getCoefficientOfVariation())





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Proba NAIS =  6.783040767375951e-06
    Current coefficient of variation =  0.09906313104733465




.. GENERATED FROM PYTHON SOURCE LINES 128-129

The length of the confidence interval of level :math:`95\%` is:

.. GENERATED FROM PYTHON SOURCE LINES 131-134

.. code-block:: Python

    length95 = result.getConfidenceLength()
    print("Confidence length (0.95) = ", result.getConfidenceLength())





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Confidence length (0.95) =  2.6339926841138084e-06




.. GENERATED FROM PYTHON SOURCE LINES 135-136

which enables to build the confidence interval:

.. GENERATED FROM PYTHON SOURCE LINES 138-146

.. code-block:: Python

    print(
        "Confidence interval (0.95) = [",
        proba - length95 / 2,
        ", ",
        proba + length95 / 2,
        "]",
    )





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Confidence interval (0.95) = [ 5.466044425319046e-06 ,  8.100037109432855e-06 ]




.. GENERATED FROM PYTHON SOURCE LINES 147-148

You can also get the successive thresholds used by the algorithm:

.. GENERATED FROM PYTHON SOURCE LINES 150-153

.. code-block:: Python

    levels = algo.getThresholdPerStep()
    print("Levels of g = ", levels)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Levels of g =  [3.40173,2.03082,0.448448,0]




.. GENERATED FROM PYTHON SOURCE LINES 154-156

Draw the NAIS samples used by the algorithm
---------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 158-159

You can get the number :math:`N_s` of steps with:

.. GENERATED FROM PYTHON SOURCE LINES 159-162

.. code-block:: Python

    Ns = algo.getStepsNumber()
    print("Number of steps= ", Ns)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Number of steps=  4




.. GENERATED FROM PYTHON SOURCE LINES 163-164

Get all the inputs where :math:`g` was evaluated at each step

.. GENERATED FROM PYTHON SOURCE LINES 164-168

.. code-block:: Python

    list_subSamples = list()
    for step in range(Ns):
        list_subSamples.append(algo.getInputSample(step))








.. GENERATED FROM PYTHON SOURCE LINES 169-170

The following graph draws each NAIS sample and the frontier :math:`g(x_1, x_2) = l_i` where :math:`l_i` is the threshold at the step :math:`i`:

.. GENERATED FROM PYTHON SOURCE LINES 172-180

.. code-block:: Python

    graph = ot.Graph()
    graph.setAxes(True)
    graph.setGrid(True)
    graph.setTitle("NAIS sampling: samples")
    graph.setXTitle(r"$x_1$")
    graph.setYTitle(r"$x_2$")
    graph.setLegendPosition("lower left")








.. GENERATED FROM PYTHON SOURCE LINES 181-182

Add all the NAIS samples:

.. GENERATED FROM PYTHON SOURCE LINES 184-190

.. code-block:: Python

    for i in range(Ns):
        cloud = ot.Cloud(list_subSamples[i])
        # cloud.setPointStyle("dot")
        graph.add(cloud)
    col = graph.getColors()








.. GENERATED FROM PYTHON SOURCE LINES 191-192

Add the frontiers :math:`g(x_1, x_2) = l_i` where :math:`l_i` is the threshold at the step :math:`i`:

.. GENERATED FROM PYTHON SOURCE LINES 194-204

.. code-block:: Python

    gIsoLines = g.draw([-5] * 2, [5] * 2, [128] * 2)
    dr = gIsoLines.getDrawable(0)
    for i, lv in enumerate(levels):
        dr.setLevels([lv])
        dr.setLineStyle("solid")
        dr.setLegend(r"$g(X) = $" + str(round(lv, 2)))
        dr.setLineWidth(3)
        dr.setColor(col[i])
        graph.add(dr)








.. GENERATED FROM PYTHON SOURCE LINES 205-207

.. code-block:: Python

    _ = View(graph)




.. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_nais_002.png
   :alt: NAIS sampling: samples
   :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_nais_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 208-212

Draw the frontiers only
-----------------------

The following graph enables to understand the progression of the algorithm:

.. GENERATED FROM PYTHON SOURCE LINES 214-233

.. code-block:: Python

    graph = ot.Graph()
    graph.setAxes(True)
    graph.setGrid(True)
    dr = gIsoLines.getDrawable(0)
    for i, lv in enumerate(levels):
        dr.setLevels([lv])
        dr.setLineStyle("solid")
        dr.setLegend(r"$g(X) = $" + str(round(lv, 2)))
        dr.setLineWidth(3)
        graph.add(dr)

    graph.setColors(col)
    graph.setLegendPosition("lower left")
    graph.setTitle("NAIS sampling: thresholds")
    graph.setXTitle(r"$x_1$")
    graph.setYTitle(r"$x_2$")

    _ = View(graph)




.. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_nais_003.png
   :alt: NAIS sampling: thresholds
   :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_nais_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 234-237

Get all the input and output points that realized the event
-----------------------------------------------------------
The following lines are possible only if you have mentioned that you wanted to keep samples with the method *algo.setKeepSample(True)*

.. GENERATED FROM PYTHON SOURCE LINES 239-245

.. code-block:: Python

    select = ot.NAIS.EVENT1  # points that realize the event
    step = Ns - 1  # get the working sample from last iteration
    inputEventSample = algo.getInputSample(step, select)
    outputEventSample = algo.getOutputSample(step, select)
    print("Number of event realizations = ", inputEventSample.getSize())





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Number of event realizations =  441




.. GENERATED FROM PYTHON SOURCE LINES 246-247

Draw them! They are all in the event space.

.. GENERATED FROM PYTHON SOURCE LINES 249-263

.. code-block:: Python

    graph = ot.Graph()
    graph.setAxes(True)
    graph.setGrid(True)
    cloud = ot.Cloud(inputEventSample)
    cloud.setPointStyle("bullet")
    graph.add(cloud)
    gIsoLines = g.draw([-5] * 2, [5] * 2, [1000] * 2)
    dr = gIsoLines.getDrawable(0)
    dr.setLevels([0.0])
    dr.setColor("red")
    graph.add(dr)
    _ = View(graph)

    View.ShowAll()



.. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_nais_004.png
   :alt: plot nais
   :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_nais_004.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 2.842 seconds)


.. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_nais.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_nais.ipynb <plot_nais.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_nais.py <plot_nais.py>`