P1LagrangeInterpolation can make Python fail
Here is what I get in my console when I try to use the P1LagrangeInterpolation constructor :

_______________________________________________________
A61372@dsp1021827:~$ python
Python 3.10.6 | packaged by conda-forge | (main, Aug 22 2022, 20:35:26) [GCC 10.4.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>> import numpy as np
>> import openturns as ot
>> 
>> tps_max    = 2500
>> mesher     = ot.IntervalMesher([50])
>> interval   = ot.Interval([0.0], [tps_max])
>> outputMesh = mesher.build(interval)
>> 
>> output_dimension = 1
>> 
>> nodes            = ot.Sample(1, np.arange(0,2500,0.8))
>> nombrePasDeTemps = nodes.getDimension()
>> simplices        = [[nodeIndex, nodeIndex + 1] for nodeIndex in range(nombrePasDeTemps - 1)]
>> inputMesh        = ot.Mesh(nodes, simplices)
>> interpolationTf  = ot.P1LagrangeInterpolation(inputMesh, outputMesh, output_dimension)

Erreur de segmentation
A61372@dsp1021827:~$
___________________________________________________________

Is there any workaround for this or another class allowing to interpolate (even roughly) ?

Thank you
diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 63a272d70..cc5fdd485 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -31,8 +31,8 @@ jobs:
         -DSWIG_COMPILE_FLAGS="-O1 -Wno-unused-parameter -Wno-unused-function -Wno-missing-field-initializers" \
         -DCMAKE_UNITY_BUILD=ON -DCMAKE_UNITY_BUILD_BATCH_SIZE=32 .
         make install -j3
-        #make tests -j3
-        OMP_NUM_THREADS=1 ctest -R pyinstall -E GeneralizedParetoFactory_std -j3 --output-on-failure --timeout 100 --schedule-random
+        make tests -j3
+        OMP_NUM_THREADS=1 ctest -E GeneralizedParetoFactory_std -j3 --output-on-failure --timeout 100 --schedule-random
 
   windows:
     runs-on: windows-latest
diff --git a/ChangeLog b/ChangeLog
index 1f751dc9d..7226eb47b 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -35,6 +35,7 @@
  * Deprecated GeneralLinearModelAlgorithm | KrigingAlgorithm ctors using collection of basis
  * Deprecated Drawable.draw|clean & Graph.draw|clean|getRCommand methods for legacy R graphs
  * Deprecated GaussianProcess.GIBBS in favor of GaussianProcess.GALLIGAOGIBBS
+ * Deprecated FunctionalChaosSobolIndices::summary in favor of __str__
 
 === Documentation ===
  * Added example galleries at the end of API doc pages
diff --git a/TODO b/TODO
index b88b9b44c..34fbbbb94 100644
--- a/TODO
+++ b/TODO
@@ -12,3 +12,4 @@ Remove IterativeMoments::getCenteredMoments
 Drop deprecated GeneralLinearModelAlgorithm | KrigingAlgorithm ctors for 1.21
 Drop deprecated Drawable.draw|clean & Graph.draw|clean|getRCommand|makeR*|GetExtensionMap methods
 Drop deprecated GaussianProcess.GIBBS
+Remove FunctionalChaosSobolIndices::summary
diff --git a/lib/src/Base/Algo/EnclosingSimplexAlgorithmImplementation.cxx b/lib/src/Base/Algo/EnclosingSimplexAlgorithmImplementation.cxx
index 653fc5e25..4d59a476c 100644
--- a/lib/src/Base/Algo/EnclosingSimplexAlgorithmImplementation.cxx
+++ b/lib/src/Base/Algo/EnclosingSimplexAlgorithmImplementation.cxx
@@ -105,9 +105,12 @@ void EnclosingSimplexAlgorithmImplementation::setVerticesAndSimplices(const Samp
     }
     for (IndicesCollection::const_iterator cit = simplices_.cbegin_at(i), guard = simplices_.cend_at(i); cit != guard; ++cit)
     {
+      const UnsignedInteger j = *cit;
+      if (j >= vertices_.getSize())
+        throw InvalidArgumentException(HERE) << "the mesh simplices reference more vertices than provided";
       for(UnsignedInteger k = 0; k < dimension; ++k)
       {
-        const Scalar coordinate = vertices_(*cit, k);
+        const Scalar coordinate = vertices_(j, k);
         lower[k] = std::min(lower[k], coordinate);
         upper[k] = std::max(upper[k], coordinate);
       }
diff --git a/lib/src/Base/Stat/ProcessSampleImplementation.cxx b/lib/src/Base/Stat/ProcessSampleImplementation.cxx
index 8c2f6b482..47a2ebe47 100644
--- a/lib/src/Base/Stat/ProcessSampleImplementation.cxx
+++ b/lib/src/Base/Stat/ProcessSampleImplementation.cxx
@@ -523,6 +523,9 @@ ProcessSampleImplementation ProcessSampleImplementation::computeQuantilePerCompo
 /* Extract the sample of values at the given vertex index */
 Sample ProcessSampleImplementation::getSampleAtVertex(const UnsignedInteger index) const
 {
+  const UnsignedInteger verticesNumber = mesh_.getVerticesNumber();
+  if (index >= verticesNumber)
+    throw OutOfBoundException(HERE) << "Index (" << index << ") is not less than vertices number (" << verticesNumber << ")";
   const UnsignedInteger size = getSize();
   const UnsignedInteger dimension = getDimension();
   Sample result(size, dimension);
diff --git a/lib/src/Uncertainty/Algorithm/MetaModel/FunctionalChaos/FunctionalChaosSobolIndices.cxx b/lib/src/Uncertainty/Algorithm/MetaModel/FunctionalChaos/FunctionalChaosSobolIndices.cxx
index 98dcd6e9c..c21d1bcbc 100644
--- a/lib/src/Uncertainty/Algorithm/MetaModel/FunctionalChaos/FunctionalChaosSobolIndices.cxx
+++ b/lib/src/Uncertainty/Algorithm/MetaModel/FunctionalChaos/FunctionalChaosSobolIndices.cxx
@@ -66,6 +66,12 @@ inline bool varianceComparison(const std::pair<UnsignedInteger, Scalar> &a, cons
 }
 
 String FunctionalChaosSobolIndices::summary() const
+{
+  LOGWARN(OSS() << "FunctionalChaosSobolIndices::summary is deprecated, use __str__()");
+  return __str__();
+}
+
+String FunctionalChaosSobolIndices::__str__(const String & /*offset*/) const
 {
   const UnsignedInteger inputDimension = functionalChaosResult_.getDistribution().getDimension();
   const UnsignedInteger outputDimension = functionalChaosResult_.getMetaModel().getOutputDimension();
diff --git a/lib/src/Uncertainty/Algorithm/MetaModel/FunctionalChaos/openturns/FunctionalChaosSobolIndices.hxx b/lib/src/Uncertainty/Algorithm/MetaModel/FunctionalChaos/openturns/FunctionalChaosSobolIndices.hxx
index 3bc73c6f2..867fcda0f 100644
--- a/lib/src/Uncertainty/Algorithm/MetaModel/FunctionalChaos/openturns/FunctionalChaosSobolIndices.hxx
+++ b/lib/src/Uncertainty/Algorithm/MetaModel/FunctionalChaos/openturns/FunctionalChaosSobolIndices.hxx
@@ -48,8 +48,9 @@ public:
 
   /** String converter */
   String __repr__() const override;
+  String __str__(const String & offset = "") const override;
 
-  /** Summary table */
+  /** @deprecated Summary table */
   String summary() const;
 
   /* Here is the interface that all derived class must implement */
diff --git a/lib/src/Uncertainty/Algorithm/MetaModel/LinearModel/LinearModelStepwiseAlgorithm.cxx b/lib/src/Uncertainty/Algorithm/MetaModel/LinearModel/LinearModelStepwiseAlgorithm.cxx
index 8a865dbfb..1f03a75ff 100644
--- a/lib/src/Uncertainty/Algorithm/MetaModel/LinearModel/LinearModelStepwiseAlgorithm.cxx
+++ b/lib/src/Uncertainty/Algorithm/MetaModel/LinearModel/LinearModelStepwiseAlgorithm.cxx
@@ -394,6 +394,13 @@ void LinearModelStepwiseAlgorithm::run()
   const Sample fx(f(inputSample_));
   LOGDEBUG(OSS() << "Total number of columns=" << fx.getDimension());
 
+  // check for null basis term
+  const Point mean = fx.computeMean();
+  const Point stddev = fx.computeStandardDeviation();
+  for (UnsignedInteger i = 0; i < fx.getDimension(); ++ i)
+    if (!(stddev[i] > 0.0) && (mean[i] == 0.0))
+      throw InvalidArgumentException(HERE) << "Null basis term at index " << i;
+
   {
     // Reduce scope of Xt
     const Matrix Xt(fx.getDimension(), size, fx.getImplementation()->getData());
diff --git a/lib/src/Uncertainty/Algorithm/Sensitivity/SobolIndicesAlgorithmImplementation.cxx b/lib/src/Uncertainty/Algorithm/Sensitivity/SobolIndicesAlgorithmImplementation.cxx
index 5246723d5..8280e0e17 100644
--- a/lib/src/Uncertainty/Algorithm/Sensitivity/SobolIndicesAlgorithmImplementation.cxx
+++ b/lib/src/Uncertainty/Algorithm/Sensitivity/SobolIndicesAlgorithmImplementation.cxx
@@ -740,8 +740,6 @@ Graph SobolIndicesAlgorithmImplementation::DrawImportanceFactors(const Point & v
   Point data(dimension);
   /* Normalization */
   for (UnsignedInteger i = 0; i < dimension; ++i) data[i] = values[i] / l1Norm;
-  /* we build the pie */
-  Pie importanceFactorsPie(data);
 
   /* build labels and colors for the pie */
   Description palette(dimension);
@@ -761,9 +759,8 @@ Graph SobolIndicesAlgorithmImplementation::DrawImportanceFactors(const Point & v
     oss << 100.0 * data[i] << "%";
     labels[i] = oss;
   }
-  /* we complete the pie */
-  importanceFactorsPie.setLabels(labels);
-  importanceFactorsPie.buildDefaultPalette();
+  /* we build the pie */
+  Pie importanceFactorsPie(data, labels);
   /* we build the graph with a title */
   Graph importanceFactorsGraph(title);
   /* we embed the pie into the graph */
diff --git a/python/doc/about.rst b/python/doc/about.rst
index b9849ab3e..07a107cc2 100644
--- a/python/doc/about.rst
+++ b/python/doc/about.rst
@@ -43,6 +43,8 @@ People
 - Julien Schueller
 - Kamal Abboud
 - Kieran Delamotte
+- Loic Brevault
+- Mathieu Balesdent
 - Mathieu Couplet
 - Mathieu Lapointe
 - Mathieu Souchaud
diff --git a/python/doc/developer_guide/sphinx_documentation.rst b/python/doc/developer_guide/sphinx_documentation.rst
index a2339ecfe..11dfddfb5 100644
--- a/python/doc/developer_guide/sphinx_documentation.rst
+++ b/python/doc/developer_guide/sphinx_documentation.rst
@@ -29,6 +29,9 @@ matplotlib.sphinxext.plot_directive
     code blocks featured in the pages of this documentation, especially in the
     :ref:`examples section <examples>`.
 
+sphinx_gallery.gen_gallery
+    The `gen_gallery <https://sphinx-gallery.github.io/stable/gen_modules/sphinx_gallery.gen_gallery.html>`_
+    extension from Sphinx-Gallery is used to generate pages from the examples directory.
 
 .. _contribute_docstrings:
 
@@ -125,8 +128,8 @@ SWIG header file::
     OT_Distribution_doc
     ...
 
-Integration to the building suite
----------------------------------
+Integration to the build system
+-------------------------------
 
 The separate docstring SWIG header files are included in the SWIG header files
 of the `openturns repos <https://github.com/openturns>`_, so this does
@@ -146,3 +149,46 @@ SPHINX_FLAGS
 All these targets depend on the rst files located in the sources 
 (:file:`$OT_SOURCE_DIR/python/doc/*.rst`).
 
+Example gallery
+---------------
+
+Example pages are generated from Python scripts in the examples directory
+(:file:`$OT_SOURCE_DIR/python/doc/examples/*.py`).
+
+Each page allows to display notebook-like code or text cells.
+Each cell is delimited by a specific marker and the text cells are written in
+rst format, for example::
+
+    # %%
+    # First cell (text) in *rst* format inside comments
+    # Some formula :math:`\lambda = 4`
+
+    # %%
+    import openturns as ot
+    print("Second cell (code)")
+    dist = ot.Normal()
+
+Output cells are rendered at compilation time, so the script must
+not be too long to run.
+
+Figures are automatically generated from the matplotlib handles,
+for example using the openturns.viewer module::
+
+    # %%
+    from openturns.viewer import View
+    graph = dist.drawPDF()
+    view = View(graph)
+
+Note that sphinx runs the examples in the same Python process so you might want
+to reset the random generator at the beginning of your example for consistent
+results::
+
+    # %%
+    import openturns as ot
+    ot.RandomGenerator.SetSeed(0)
+
+The special variable `sphinx_gallery_thumbnail_number <https://sphinx-gallery.github.io/stable/configuration.html#choosing-thumbnail>`_
+can be used to determine which figure is used as thumbnail in the examples gallery::
+
+    # %%
+    # sphinx_gallery_thumbnail_number = 3
diff --git a/python/doc/examples/calibration/_static/axial-stressed-beam.png b/python/doc/examples/calibration/_static/axial-stressed-beam.png
deleted file mode 100644
index 715392573..000000000
Binary files a/python/doc/examples/calibration/_static/axial-stressed-beam.png and /dev/null differ
diff --git a/python/doc/examples/calibration/_static/pandas.png b/python/doc/examples/calibration/_static/pandas.png
deleted file mode 100644
index cf0e1c28f..000000000
Binary files a/python/doc/examples/calibration/_static/pandas.png and /dev/null differ
diff --git a/python/doc/examples/data_analysis/_static/axial-stressed-beam.png b/python/doc/examples/data_analysis/_static/axial-stressed-beam.png
deleted file mode 100644
index 715392573..000000000
Binary files a/python/doc/examples/data_analysis/_static/axial-stressed-beam.png and /dev/null differ
diff --git a/python/doc/examples/data_analysis/_static/pandas.png b/python/doc/examples/data_analysis/_static/pandas.png
deleted file mode 100644
index cf0e1c28f..000000000
Binary files a/python/doc/examples/data_analysis/_static/pandas.png and /dev/null differ
diff --git a/python/doc/examples/data_analysis/distribution_fitting/plot_estimate_conditional_quantile.py b/python/doc/examples/data_analysis/distribution_fitting/plot_estimate_conditional_quantile.py
index 8b9e33e30..0654c2cc1 100644
--- a/python/doc/examples/data_analysis/distribution_fitting/plot_estimate_conditional_quantile.py
+++ b/python/doc/examples/data_analysis/distribution_fitting/plot_estimate_conditional_quantile.py
@@ -4,12 +4,12 @@ Estimate a conditional quantile
 """
 # sphinx_gallery_thumbnail_number = 8
 # %%
-# From a multivariate data sample we estimate its distribution with kernel smoothing. Here we present a bivariate distribution  X= (X1, X2). We use the `computeConditionalQuantile` method to estimate the 90% quantile :math:`Q_1` of the conditional variable :math:`X1|X2` :
+# From a multivariate data sample, we estimate a distribution with kernel smoothing. Here we present a bivariate distribution  :math:`X= (X_1, X_2)`. We use the `computeConditionalQuantile` method to estimate the 90% quantile :math:`Q_1` of the conditional variable :math:`X_2|X_1` :
 #
 # .. math::
-#    Q_1 : x_2 \mapsto q_{0.9}(X_1|X_2=x_2)
+#    Q_2 : x_1 \mapsto q_{0.9}(X_2|X_1=x_1)
 #
-# We then draw the curve :math:`Q_1 : x_2 \mapsto Q_1(x_2)`. We first start with independent normals then we consider dependent marginals with a Clayton copula.
+# We then draw the curve :math:`Q_2 : x_1 \mapsto Q_2(x_1)`. We first start with independent normals then we consider dependent marginals with a Clayton copula.
 #
 
 # %%
@@ -76,7 +76,7 @@ view = viewer.View(graph)
 
 
 # %%
-# We can compute the conditional quantile of :math:`X_1 | X_2` with the `computeConditionalQuantile` method and draw it after.
+# We can compute the conditional quantile of :math:`X_2 | X_1` with the `computeConditionalQuantile` method and draw it after.
 #
 
 # %%
@@ -99,17 +99,17 @@ yex = [distX.computeConditionalQuantile(0.9, sampleObs[i]) for i in range(N)]
 # %%
 cxy_app = ot.Curve(x, yapp)
 cxy_ex = ot.Curve(x, yex)
-graph = ot.Graph('90% quantile of $X_1 | X_2=x_2$',
-                 '$x_2$', '$Q_1(x_2)$', True, '')
+graph = ot.Graph('90% quantile of $X_2 | X_1=x_1$',
+                 '$x_1$', '$Q_2(x_1)$', True, '')
 graph.add(cxy_app)
 graph.add(cxy_ex)
-graph.setLegends(["$Q_1$ kernel smoothing", "$Q_1$ exact"])
+graph.setLegends(["$Q_2$ kernel smoothing", "$Q_2$ exact"])
 graph.setLegendPosition('bottomright')
 graph.setColors(["red", "blue"])
 view = viewer.View(graph)
 
 # %%
-# In this case the Q_1 quantile is constant because of the independance of the marginals.
+# In this case the :math:`Q_2` quantile is constant because of the independance of the marginals.
 #
 
 # %%
@@ -164,7 +164,7 @@ view = viewer.View(graph)
 
 
 # %%
-# We can compute the conditional quantile of :math:`X_1 | X_2` with the `computeConditionalQuantile` method and draw it after.
+# We can compute the conditional quantile of :math:`X_2 | X_1=x1` with the `computeConditionalQuantile` method and draw it after.
 #
 
 # %%
@@ -187,11 +187,11 @@ yex = [distX.computeConditionalQuantile(0.9, sampleObs[i]) for i in range(N)]
 # %%
 cxy_app = ot.Curve(x, yapp)
 cxy_ex = ot.Curve(x, yex)
-graph = ot.Graph('90% quantile of $X_1 | X_2=x_2$',
-                 '$x_2$', '$Q_1(x_2)$', True, '')
+graph = ot.Graph('90% quantile of $X_2 | X_1=x_1$',
+                 '$x_1$', '$Q_2(x_1)$', True, '')
 graph.add(cxy_app)
 graph.add(cxy_ex)
-graph.setLegends(["$Q_1$ kernel smoothing", "$Q_1$ exact"])
+graph.setLegends(["$Q_2$ kernel smoothing", "$Q_2$ exact"])
 graph.setLegendPosition('bottomright')
 graph.setColors(["red", "blue"])
 view = viewer.View(graph)
diff --git a/python/doc/examples/functional_modeling/_static/axial-stressed-beam.png b/python/doc/examples/functional_modeling/_static/axial-stressed-beam.png
deleted file mode 100644
index 715392573..000000000
Binary files a/python/doc/examples/functional_modeling/_static/axial-stressed-beam.png and /dev/null differ
diff --git a/python/doc/examples/functional_modeling/_static/pandas.png b/python/doc/examples/functional_modeling/_static/pandas.png
deleted file mode 100644
index cf0e1c28f..000000000
Binary files a/python/doc/examples/functional_modeling/_static/pandas.png and /dev/null differ
diff --git a/python/doc/examples/graphs/_static/axial-stressed-beam.png b/python/doc/examples/graphs/_static/axial-stressed-beam.png
deleted file mode 100644
index 715392573..000000000
Binary files a/python/doc/examples/graphs/_static/axial-stressed-beam.png and /dev/null differ
diff --git a/python/doc/examples/graphs/_static/pandas.png b/python/doc/examples/graphs/_static/pandas.png
deleted file mode 100644
index cf0e1c28f..000000000
Binary files a/python/doc/examples/graphs/_static/pandas.png and /dev/null differ
diff --git a/python/doc/examples/meta_modeling/_static/axial-stressed-beam.png b/python/doc/examples/meta_modeling/_static/axial-stressed-beam.png
deleted file mode 100644
index 715392573..000000000
Binary files a/python/doc/examples/meta_modeling/_static/axial-stressed-beam.png and /dev/null differ
diff --git a/python/doc/examples/meta_modeling/_static/pandas.png b/python/doc/examples/meta_modeling/_static/pandas.png
deleted file mode 100644
index cf0e1c28f..000000000
Binary files a/python/doc/examples/meta_modeling/_static/pandas.png and /dev/null differ
diff --git a/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_ishigami.py b/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_ishigami.py
index cb15d8a22..f79925e60 100644
--- a/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_ishigami.py
+++ b/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_ishigami.py
@@ -136,7 +136,7 @@ view = viewer.View(graph)
 
 # %%
 chaosSI = ot.FunctionalChaosSobolIndices(result)
-print(chaosSI.summary())
+print(chaosSI)
 
 # %%
 # We notice the a coefficient with marginal degree equal to 6 has a significant impact on the output variance. Hence, we cannot get a satisfactory polynomial chaos with total degree less that 6.
diff --git a/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_ishigami_grouped_indices.py b/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_ishigami_grouped_indices.py
index 878c62bd6..9fe644217 100644
--- a/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_ishigami_grouped_indices.py
+++ b/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_ishigami_grouped_indices.py
@@ -59,7 +59,7 @@ metamodel = result.getMetaModel()
 
 # %%
 chaosSI = ot.FunctionalChaosSobolIndices(result)
-print(chaosSI.summary())
+print(chaosSI)
 
 # %%
 # We compute the first order indice of the group [0,1].
@@ -80,7 +80,7 @@ chaosSI.getSobolGroupedIndex([0, 1])
 0.279938 + 0.190322 + 0.130033 + 0.12058 + 0.0250262
 
 # %%
-# The difference between the previous sum and the output of `getSobolGroupedIndex` is lower than 0.01, which is the threshold used by the `summary` method.
+# The difference between the previous sum and the output of `getSobolGroupedIndex` is lower than 0.01, which is the threshold used by the `__str__` method.
 
 # %%
 # We compute the total order indice of the group [1,2].
diff --git a/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_functional_chaos.py b/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_functional_chaos.py
index 277e3d782..00a1936cf 100644
--- a/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_functional_chaos.py
+++ b/python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_functional_chaos.py
@@ -113,7 +113,7 @@ view = viewer.View(graph)
 
 # %%
 chaosSI = ot.FunctionalChaosSobolIndices(result)
-print(chaosSI.summary())
+print(chaosSI)
 
 # %%
 # Let us analyse the results of this global sensitivity analysis.
diff --git a/python/doc/examples/numerical_methods/_static/axial-stressed-beam.png b/python/doc/examples/numerical_methods/_static/axial-stressed-beam.png
deleted file mode 100644
index 715392573..000000000
Binary files a/python/doc/examples/numerical_methods/_static/axial-stressed-beam.png and /dev/null differ
diff --git a/python/doc/examples/numerical_methods/_static/pandas.png b/python/doc/examples/numerical_methods/_static/pandas.png
deleted file mode 100644
index cf0e1c28f..000000000
Binary files a/python/doc/examples/numerical_methods/_static/pandas.png and /dev/null differ
diff --git a/python/doc/examples/numerical_methods/optimization/plot_optimization_bonmin.py b/python/doc/examples/numerical_methods/optimization/plot_optimization_bonmin.py
index 9bb038335..5a1b57e45 100644
--- a/python/doc/examples/numerical_methods/optimization/plot_optimization_bonmin.py
+++ b/python/doc/examples/numerical_methods/optimization/plot_optimization_bonmin.py
@@ -5,7 +5,9 @@ Optimization using bonmin
 # %%
 
 # %%
-# In this example we are going to explore mixed-integer non linear problems optimization using OpenTURNS' `bonmin <https://www.coin-or.org/Bonmin/index.html>`_ interface.
+# In this example we are going to explore mixed-integer non linear problems
+# optimization using the `bonmin <https://www.coin-or.org/Bonmin/index.html>`_
+# interface.
 # %%
 import openturns as ot
 ot.Log.Show(ot.Log.NONE)
@@ -16,13 +18,15 @@ for algo in ot.Bonmin.GetAlgorithmNames():
     print(algo)
 
 # %%
-# Details and references on bonmin algorithms are available `here <https://projects.coin-or.org/Bonmin>`_ .
+# Details and references on bonmin algorithms are available
+# `here <https://projects.coin-or.org/Bonmin>`_ .
 
 # %%
 # Setting up and solving a simple problem
 # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 #
-# The following example will demonstrate the use of bonmin "BB" algorithm to solve the following problem:
+# The following example will demonstrate the use of bonmin "BB" algorithm to
+# solve the following problem:
 #
 # .. math::
 #    \min - x_0 - x_1 - x_2
@@ -39,9 +43,12 @@ for algo in ot.Bonmin.GetAlgorithmNames():
 #    x_3 \in \mathbb{N}
 #    \end{array}
 #
-# The theoretical minimum is reached for :math:`x = [1,1,0.5,0]`. At this point, the objective function value is :math:`-2.5`
+# The theoretical minimum is reached for :math:`x = [1,1,0.5,0]`.
+# At this point, the objective function value is :math:`-2.5`
 #
-# N.B.: OpenTURNS requires equality and inequality constraints to be stated as :math:`g(x) = 0` and :math:`h(x) \geq 0`, respectively. Thus the inequalities above will have to be restated to match this requirement:
+# N.B.: OpenTURNS requires equality and inequality constraints to be stated as
+# :math:`g(x) = 0` and :math:`h(x) \geq 0`, respectively. Thus the inequalities
+# above will have to be restated to match this requirement:
 #
 # .. math::
 #    \begin{array}{l}
@@ -64,11 +71,14 @@ bounds = ot.Interval([0, 0, 0, 0], [1, 1e99, 1e99, 5], [
 #    No equality constraint -> nothing to do
 #    Inequality constraints:
 h = ot.SymbolicFunction(['x0', 'x1', 'x2', 'x3'], [
-                        '-(x1-0.5)^2 - (x2-0.5)^2 + 0.25', 'x1 - x0', '-x0 - x2 - x3 + 2'])
+                        '-(x1-0.5)^2 - (x2-0.5)^2 + 0.25',
+                        'x1 - x0', '-x0 - x2 - x3 + 2'])
 
 # Definition of variables types
-variablesType = [ot.OptimizationProblemImplementation.BINARY, ot.OptimizationProblemImplementation.CONTINUOUS,
-                 ot.OptimizationProblemImplementation.CONTINUOUS, ot.OptimizationProblemImplementation.INTEGER]
+CONTINUOUS = ot.OptimizationProblemImplementation.CONTINUOUS
+BINARY = ot.OptimizationProblemImplementation.BINARY
+INTEGER = ot.OptimizationProblemImplementation.INTEGER
+variablesType = [BINARY, CONTINUOUS, CONTINUOUS, INTEGER]
 
 # Setting up Bonmin problem
 problem = ot.OptimizationProblem(objectiveFunction)
@@ -90,6 +100,6 @@ bonminAlgorithm.run()
 
 # Retrieving the results
 result = bonminAlgorithm.getResult()
-print(" -- Optimal point = " + result.getOptimalPoint().__str__())
-print(" -- Optimal value = " + result.getOptimalValue().__str__())
-print(" -- Evaluation number = " + result.getInputSample().getSize().__str__())
+print(" -- Optimal point = " + str(result.getOptimalPoint()))
+print(" -- Optimal value = " + str(result.getOptimalValue()))
+print(" -- Evaluation number = " + str(result.getInputSample().getSize()))
diff --git a/python/doc/examples/probabilistic_modeling/_static/axial-stressed-beam.png b/python/doc/examples/probabilistic_modeling/_static/axial-stressed-beam.png
deleted file mode 100644
index 715392573..000000000
Binary files a/python/doc/examples/probabilistic_modeling/_static/axial-stressed-beam.png and /dev/null differ
diff --git a/python/doc/examples/probabilistic_modeling/_static/pandas.png b/python/doc/examples/probabilistic_modeling/_static/pandas.png
deleted file mode 100644
index cf0e1c28f..000000000
Binary files a/python/doc/examples/probabilistic_modeling/_static/pandas.png and /dev/null differ
diff --git a/python/doc/examples/reliability_sensitivity/_static/axial-stressed-beam.png b/python/doc/examples/reliability_sensitivity/_static/axial-stressed-beam.png
deleted file mode 100644
index 715392573..000000000
Binary files a/python/doc/examples/reliability_sensitivity/_static/axial-stressed-beam.png and /dev/null differ
diff --git a/python/doc/examples/reliability_sensitivity/_static/pandas.png b/python/doc/examples/reliability_sensitivity/_static/pandas.png
deleted file mode 100644
index cf0e1c28f..000000000
Binary files a/python/doc/examples/reliability_sensitivity/_static/pandas.png and /dev/null differ
diff --git a/python/doc/examples/reliability_sensitivity/sensitivity_analysis/plot_functional_chaos_sensitivity.py b/python/doc/examples/reliability_sensitivity/sensitivity_analysis/plot_functional_chaos_sensitivity.py
index e59bd3b22..13894389e 100644
--- a/python/doc/examples/reliability_sensitivity/sensitivity_analysis/plot_functional_chaos_sensitivity.py
+++ b/python/doc/examples/reliability_sensitivity/sensitivity_analysis/plot_functional_chaos_sensitivity.py
@@ -73,7 +73,7 @@ print(result.getRelativeErrors())
 # %%
 # Quick summary of sensitivity analysis
 sensitivityAnalysis = ot.FunctionalChaosSobolIndices(result)
-print(sensitivityAnalysis.summary())
+print(sensitivityAnalysis)
 
 # %%
 # draw Sobol' indices
diff --git a/python/src/Distribution.i b/python/src/Distribution.i
index 6876656d8..e599b3dc5 100644
--- a/python/src/Distribution.i
+++ b/python/src/Distribution.i
@@ -82,6 +82,7 @@ Distribution __rtruediv__(Scalar s) { return self->inverse() * s; }
 
 %pythoncode %{
 from openturns.typ import Interval
+import math
 
 class PythonDistribution(object):
     """
@@ -140,21 +141,20 @@ class SciPyDistribution(PythonDistribution):
     """
     def __init__(self, dist):
         super(SciPyDistribution, self).__init__(1)
-        if dist.__class__.__name__ not in ('rv_frozen', 'rv_continuous_frozen', 'rv_histogram'):
-            raise TypeError('Argument is not a continuous scipy 1D distribution')
+        if dist.__class__.__name__ not in ('rv_frozen', 'rv_continuous_frozen', 'rv_histogram', 'rv_discrete_frozen'):
+            raise TypeError('Argument is not a scipy 1D distribution')
         self._dist = dist
 
         # compute range
         from openturns import ResourceMap
         cdf_epsilon = ResourceMap.GetAsScalar('Distribution-DefaultCDFEpsilon')
-        lb = dist.ppf(0.0)
-        ub = dist.ppf(1.0)
+        lb, ub = dist.support()
         flb = lb != float('-inf')
         fub = ub != float('+inf')
         if not flb:
             lb = dist.ppf(cdf_epsilon)
         if not fub:
-            ub = dist.ppf(1.0 - cdf_epsilon)
+            ub = dist.ppf(0.5 + (0.5 - cdf_epsilon))
         self.__range = Interval([lb], [ub])
         self.__range.setFiniteLowerBound([int(flb)])
         self.__range.setFiniteUpperBound([int(fub)])
@@ -162,6 +162,13 @@ class SciPyDistribution(PythonDistribution):
     def getRange(self):
         return self.__range
 
+    def getSupport(self, interval):
+        # assume discrete distributions are integral
+        intersection = self.getRange().intersect(interval)
+        kMin = math.ceil(intersection.getLowerBound()[0])
+        kMax = math.floor(intersection.getUpperBound()[0])
+        return [[k] for k in range(kMin, kMax + 1)]
+
     def getRealization(self):
         rvs = self._dist.rvs()
         return [rvs]
@@ -171,9 +178,15 @@ class SciPyDistribution(PythonDistribution):
         return rvs.reshape(size, 1)
 
     def computePDF(self, X):
-        # before 1.9, all 1d dist have a pdf method
-        # since 1.9, pdf for continuous, pmf for discrete
-        pdf = self._dist.pdf(X[0])
+        x0 = float(X[0])
+        if self.isContinuous():
+            pdf = self._dist.pdf(x0)
+        else:
+            if x0.is_integer():
+                pdf = self._dist.pmf(x0)
+            else:
+                # newer scipy pmf is a continuous function
+                pdf = 0.0
         return pdf
 
     def computeCDF(self, X):
@@ -202,7 +215,14 @@ class SciPyDistribution(PythonDistribution):
         return [moment]
 
     def isContinuous(self):
-        return True
+        return hasattr(self._dist.dist, 'pdf')
+
+    def isDiscrete(self):
+        return not self.isContinuous()
+
+    def isIntegral(self):
+        # assume discrete distributions are integral
+        return self.isDiscrete()
 
     def computeScalarQuantile(self, p, tail=False):
         q = self._dist.ppf(1.0 - p if tail else p)
diff --git a/python/src/FunctionalChaosSobolIndices_doc.i.in b/python/src/FunctionalChaosSobolIndices_doc.i.in
index 1d154c8d7..2ad25269e 100644
--- a/python/src/FunctionalChaosSobolIndices_doc.i.in
+++ b/python/src/FunctionalChaosSobolIndices_doc.i.in
@@ -155,13 +155,3 @@ Returns
 -------
 s : float
     The grouped Sobol total order index."
-
-// ---------------------------------------------------------------------
-
-%feature("docstring") OT::FunctionalChaosSobolIndices::summary
-"Summary accessor.
-
-Returns
--------
-summary : str
-    A text summary of the sensitivity analysis to be shown in a console."
diff --git a/python/src/InverseWishart_doc.i.in b/python/src/InverseWishart_doc.i.in
index 8f4b33e9d..66eddda4c 100644
--- a/python/src/InverseWishart_doc.i.in
+++ b/python/src/InverseWishart_doc.i.in
@@ -8,14 +8,40 @@ V : 2-d sequence of float
 nu : float, :math:`\nu > p - 1`
     Degrees of freedom.
 
+See Also
+--------
+Wishart
+
+Notes
+-----
+The distribution is of dimension :math:`\frac{p(p+1)}{2}` (without the symmetry terms).
+
 Examples
 --------
 Create a distribution:
 
 >>> import openturns as ot
->>> distribution = ot.InverseWishart(ot.CovarianceMatrix([[1.0]]), 3.0)
+>>> ot.RandomGenerator.SetSeed(0)
+>>> p = 3
+>>> nu = p + 1.0
+>>> distribution = ot.InverseWishart(ot.CovarianceMatrix(p), nu)
+
+Get a realization as covariance matrix:
+
+>>> matrix = distribution.getRealizationAsMatrix()
+>>> print(matrix)
+[[  0.492128 -0.260455 -0.253392 ]
+ [ -0.260455  1.44373  -0.799136 ]
+ [ -0.253392 -0.799136  1.07088  ]]
+
+Get a realization (to be interpreted as the lower triangle terms of the matrix,
+line by line as a flat vector):
+
+>>> flat = distribution.getRealization()
+>>> print(flat)
+[0.241637,0.0606503,0.190742,0.0930118,-0.890494,6.94252]
 
-Draw a sample:
+Draw a sample (also in flat format):
 
 >>> sample = distribution.getSample(5)"
 
diff --git a/python/src/KernelSmoothing_doc.i.in b/python/src/KernelSmoothing_doc.i.in
index 6aa585d6a..c3ba32d5f 100644
--- a/python/src/KernelSmoothing_doc.i.in
+++ b/python/src/KernelSmoothing_doc.i.in
@@ -33,6 +33,7 @@ Examples
 Fit a distribution on data thanks to the kernel smoothing technique:
 
 >>> import openturns as ot
+>>> ot.RandomGenerator.SetSeed(0)
 >>> sample = ot.Gamma(6.0, 1.0).getSample(100)
 >>> ks = ot.KernelSmoothing()
 >>> fittedDist = ks.build(sample)
@@ -46,7 +47,7 @@ Get the bandwidth:
 
 >>> bandwidth = ks.getBandwidth()
 >>> print(bandwidth)
-[1.03374]
+[0.862207]
 
 The bandwidth was guessed by the :meth:`build` method.
 It could also have been provided by the user.
diff --git a/python/src/PythonDistribution.cxx b/python/src/PythonDistribution.cxx
index d1a476e1c..38ea3c4e9 100644
--- a/python/src/PythonDistribution.cxx
+++ b/python/src/PythonDistribution.cxx
@@ -935,6 +935,8 @@ Sample PythonDistribution::getSupport(const Interval & interval) const
 {
   if (PyObject_HasAttrString(pyObj_, const_cast<char *>("getSupport")))
   {
+    if (interval.getDimension() != getDimension())
+      throw InvalidArgumentException(HERE) << "Error: the given interval has a dimension that does not match the distribution dimension.";
     ScopedPyObjectPointer methodName(convert< String, _PyString_ >("getSupport"));
     ScopedPyObjectPointer pyInterval(SWIG_NewPointerObj(new Interval(interval), SWIG_TypeQuery("OT::Interval *"), SWIG_POINTER_OWN));
     ScopedPyObjectPointer callResult(PyObject_CallMethodObjArgs(pyObj_, methodName.get(), pyInterval.get(), NULL));
diff --git a/python/src/Wishart_doc.i.in b/python/src/Wishart_doc.i.in
index de6edea96..cf0ac129d 100644
--- a/python/src/Wishart_doc.i.in
+++ b/python/src/Wishart_doc.i.in
@@ -8,14 +8,40 @@ V : 2-d sequence of float
 nu : float, :math:`\nu > p - 1`
     Degrees of freedom.
 
+See Also
+--------
+InverseWishart
+
+Notes
+-----
+The distribution is of dimension :math:`\frac{p(p+1)}{2}` (without the symmetry terms).
+
 Examples
 --------
 Create a distribution:
 
 >>> import openturns as ot
->>> distribution = ot.Wishart(ot.CovarianceMatrix([[1.0]]), 3.0)
+>>> ot.RandomGenerator.SetSeed(0)
+>>> p = 3
+>>> nu = p + 1.0
+>>> distribution = ot.Wishart(ot.CovarianceMatrix(p), nu)
+
+Get a realization as covariance matrix:
+
+>>> matrix = distribution.getRealizationAsMatrix()
+>>> print(matrix)
+[[ 5.16322 2.73917 3.26582 ]
+ [ 2.73917 2.63329 2.61323 ]
+ [ 3.26582 2.61323 3.65668 ]]
+
+Get a realization (to be interpreted as the lower triangle terms of the matrix,
+line by line as a flat vector):
+
+>>> flat = distribution.getRealization()
+>>> print(flat)
+[5.82928,-5.5291,18.3127,-0.787297,2.42298,0.465375]
 
-Draw a sample:
+Draw a sample (also in flat format):
 
 >>> sample = distribution.getSample(5)"
 
diff --git a/python/src/coupling_tools.py b/python/src/coupling_tools.py
index ad34117bb..258edc3aa 100644
--- a/python/src/coupling_tools.py
+++ b/python/src/coupling_tools.py
@@ -35,7 +35,6 @@ import os
 import shlex
 import subprocess
 import sys
-import warnings
 
 debug = False
 default_encoding = sys.getdefaultencoding()
@@ -46,10 +45,16 @@ def check_param(obj, obj_type):
 
     try:
         obj_type(obj)
-    except:
-        raise AssertionError("error: parameter (" +
-                             str(type(obj)) + ": " + str(obj) +
-                             ") must be of " + str(obj_type) + "!")
+    except Exception:
+        raise AssertionError(
+            "error: parameter ("
+            + str(type(obj))
+            + ": "
+            + str(obj)
+            + ") must be of "
+            + str(obj_type)
+            + "!"
+        )
 
 
 def replace(infile, outfile, tokens, values, formats=None, encoding=default_encoding):
@@ -107,15 +112,15 @@ def replace(infile, outfile, tokens, values, formats=None, encoding=default_enco
         if len(values) != len(formats):
             raise AssertionError("Error: values size is != formats size!")
     else:
-        formats = ['{0}'] * len(values)
+        formats = ["{0}"] * len(values)
 
     inplace = False
     if (outfile is None) or (infile == outfile):
         inplace = True
         outfile = infile + ".temporary_outfile"
 
-    infile_handle = open(infile, 'rb')
-    outfile_handle = open(outfile, 'wb')
+    infile_handle = open(infile, "rb")
+    outfile_handle = open(outfile, "wb")
 
     regex_tokens = []
     found_tokens = []
@@ -130,8 +135,7 @@ def replace(infile, outfile, tokens, values, formats=None, encoding=default_enco
             found = regex_token.search(line)
             while found is not None:
                 found_tokens[i] = True
-                line = line.replace(
-                    found.group(0), formats[i].format(values[i]))
+                line = line.replace(found.group(0), formats[i].format(values[i]))
                 found = regex_token.search(line)
             i += 1
 
@@ -151,14 +155,23 @@ def replace(infile, outfile, tokens, values, formats=None, encoding=default_enco
 
 class OTCalledProcessError(subprocess.CalledProcessError):
     def __str__(self):
-        err_msg = (':\n' + self.stderr[:200].decode()
-                   ) if self.stderr is not None else ''
+        err_msg = (
+            (":\n" + self.stderr[:200].decode()) if self.stderr is not None else ""
+        )
         return super(OTCalledProcessError, self).__str__() + err_msg
 
 
-def execute(cmd, cwd=None, shell=False, executable=None, hide_win=True,
-            check=True, capture_output=False,
-            timeout=None, env=None):
+def execute(
+    cmd,
+    cwd=None,
+    shell=False,
+    executable=None,
+    hide_win=True,
+    check=True,
+    capture_output=False,
+    timeout=None,
+    env=None,
+):
     """
     Launch an external process.
 
@@ -206,24 +219,27 @@ def execute(cmd, cwd=None, shell=False, executable=None, hide_win=True,
     """
 
     # split cmd if not in a shell before passing it to os.execvp()
-    try:
-        import posix
-        posix = True
-    except ImportError:
-        posix = False
+    posix = os.name == "posix"
     process_args = cmd if shell else shlex.split(cmd, posix=posix)
 
     # override startupinfo to hide windows console
     startupinfo = None
-    if hide_win and sys.platform.startswith('win'):
+    if hide_win and sys.platform.startswith("win"):
         startupinfo = subprocess.STARTUPINFO()
         startupinfo.dwFlags |= subprocess.STARTF_USESHOWWINDOW
 
     stdout = subprocess.PIPE if capture_output else None
     stderr = subprocess.PIPE if capture_output else None
-    process = subprocess.Popen(process_args, shell=shell, cwd=cwd,
-                               executable=executable, stdout=stdout, stderr=stderr,
-                               startupinfo=startupinfo, env=env)
+    process = subprocess.Popen(
+        process_args,
+        shell=shell,
+        cwd=cwd,
+        executable=executable,
+        stdout=stdout,
+        stderr=stderr,
+        startupinfo=startupinfo,
+        env=env,
+    )
     stdout_data = None
     stderr_data = None
     try:
@@ -247,7 +263,9 @@ def execute(cmd, cwd=None, shell=False, executable=None, hide_win=True,
     if check and returncode:
         raise OTCalledProcessError(returncode, cmd, stdout_data, stderr_data)
 
-    return subprocess.CompletedProcess(cmd, returncode, stdout=stdout_data, stderr=stderr_data)
+    return subprocess.CompletedProcess(
+        cmd, returncode, stdout=stdout_data, stderr=stderr_data
+    )
 
 
 def get_regex(filename, patterns, encoding=default_encoding):
@@ -286,7 +304,7 @@ def get_regex(filename, patterns, encoding=default_encoding):
     >>> with open('results.out', 'w') as f:
     ...     count = f.write('@E=-9.5E3')
     >>> # parse file with regex
-    >>> ct.get_regex('results.out', patterns=[r'@E=(\R)'])
+    >>> ct.get_regex('results.out', patterns=[r"(\\R)"])
     [-9500.0]
     """
 
@@ -298,14 +316,14 @@ def get_regex(filename, patterns, encoding=default_encoding):
     re_patterns = []
     for pattern in patterns:
         # OT-like shortcuts
-        pattern = pattern.replace(r'\R',
-                                  r'[+-]? *(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?')
-        pattern = pattern.replace(r'\I',
-                                  r'[+-]? *\d+')
+        pattern = pattern.replace(
+            r"\R", r"[+-]? *(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?"
+        )
+        pattern = pattern.replace(r"\I", r"[+-]? *\d+")
 
         re_patterns.append(re.compile(pattern))
 
-    file_handle = open(filename, 'rb')
+    file_handle = open(filename, "rb")
 
     for line in file_handle:
         line = line.decode(encoding)
@@ -320,7 +338,7 @@ def get_regex(filename, patterns, encoding=default_encoding):
 
     for result, pattern in zip(results, patterns):
         if result is None:
-            raise EOFError('error: no pattern (' + pattern + ') found!')
+            raise EOFError("error: no pattern (" + pattern + ") found!")
 
     return results
 
@@ -339,23 +357,19 @@ def get_real_from_line(line):
     value : float or int
         Found value. Raise an exception if nothing found.
     """
-    result = None
-    # print "token=" + str(token) + ', skip_token=' + str(skip_token) + ','\
-    #        'skip_line=' + str(skip_line) + ', skip_col=' + str(skip_col)
 
     # \S*: spaces are allowed at the beginning of the real
-    real_regex = '\s*[+-]? *(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?'
-    # real_regex = '[+-]? *(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?'
+    real_regex = r"\s*[+-]? *(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?"
     re_real = re.compile(real_regex)
 
     # get the value
     match = re_real.match(line)
     if match:
-        result = float(line[match.start():match.end()])
-
-    if result is None:
-        raise EOFError('error: real not found at the beginning of this line: '
-                       '(' + line + ')!')
+        result = float(line[match.start(): match.end()])
+    else:
+        raise EOFError(
+            "error: real not found at the beginning of this line: " "(" + line + ")!"
+        )
 
     return result
 
@@ -373,30 +387,32 @@ def read_line(handle, seek=0, encoding=default_encoding):
     """
     line = handle.readline().decode(encoding)
 
-    if seek < 0 and line != '':
+    if seek < 0 and line != "":
         if handle.tell() > -seek:
             if debug:
-                sys.stderr.write('line before cut ->' + line + '<-\n')
+                sys.stderr.write("line before cut ->" + line + "<-\n")
 
             # cut the end of the line
-            line = line[:len(line) - (handle.tell() + seek)]
+            line = line[: len(line) - (handle.tell() + seek)]
             # after cut, if the line is empty, add a \n to behave like
             # readline() python function
-            if line == '':
-                line = '\n'
+            if line == "":
+                line = "\n"
             # any further read will return ''
             handle.seek(0, os.SEEK_END)
 
             if debug:
-                sys.stderr.write('line after cut ->' + line + '<-\n')
+                sys.stderr.write("line after cut ->" + line + "<-\n")
 
     if debug:
-        sys.stderr.write('read_line: ->' + line + '<-\n')
+        sys.stderr.write("read_line: ->" + line + "<-\n")
 
     return line
 
 
-def get_line_col(filename, skip_line=0, skip_col=0, col_sep=None, seek=0, encoding=default_encoding):
+def get_line_col(
+    filename, skip_line=0, skip_col=0, col_sep=None, seek=0, encoding=default_encoding
+):
     """
     Get a value at specific line/columns coordinates.
 
@@ -442,13 +458,20 @@ def get_line_col(filename, skip_line=0, skip_col=0, col_sep=None, seek=0, encodi
     check_param(skip_col, int)
     check_param(seek, int)
 
-    handle = open(filename, 'rb')
+    handle = open(filename, "rb")
     if seek > 0:
         handle.seek(seek)
 
     if debug:
-        sys.stderr.write('get_line_col(skip_line=' + str(skip_line) +
-                         ', skip_col=' + str(skip_col) + ', seek=' + str(seek) + ')\n')
+        sys.stderr.write(
+            "get_line_col(skip_line="
+            + str(skip_line)
+            + ", skip_col="
+            + str(skip_col)
+            + ", seek="
+            + str(seek)
+            + ")\n"
+        )
 
     # skip line backward
     if skip_line < 0:
@@ -461,7 +484,7 @@ def get_line_col(filename, skip_line=0, skip_col=0, col_sep=None, seek=0, encodi
         # build lines cache
         previous_pos = handle.tell()
         line = read_line(handle, seek, encoding)
-        while line != '':
+        while line != "":
             # append to cache
             lines_cache.append(previous_pos)
             if len(lines_cache) > lines_cache_size:
@@ -471,26 +494,26 @@ def get_line_col(filename, skip_line=0, skip_col=0, col_sep=None, seek=0, encodi
 
             line = read_line(handle, seek, encoding)
             if debug:
-                sys.stderr.write('lines_cache: ' + str(lines_cache) + '\n')
+                sys.stderr.write("lines_cache: " + str(lines_cache) + "\n")
 
         if len(lines_cache) < lines_cache_size:
-            err_msg = 'error: the file has less than ' + str(lines_cache_size) +\
-                      ' lines!'
+            err_msg = (
+                "error: the file has less than " + str(lines_cache_size) + " lines!"
+            )
             handle.close()
             raise EOFError(err_msg)
         else:
             handle.seek(lines_cache[0])
             line_found = read_line(handle, seek, encoding)
             if debug:
-                sys.stderr.write('line found: ->' + line_found + '<-\n')
+                sys.stderr.write("line found: ->" + line_found + "<-\n")
     # skip line forward
     else:
         while skip_line >= 0:
             line = read_line(handle, seek, encoding)
-            if skip_line > 0 and line == '':
+            if skip_line > 0 and line == "":
                 handle.close()
-                raise EOFError(
-                    'error: the file has less lines than "skip_line"!')
+                raise EOFError('error: the file has less lines than "skip_line"!')
             skip_line -= 1
         line_found = line
 
@@ -500,9 +523,8 @@ def get_line_col(filename, skip_line=0, skip_col=0, col_sep=None, seek=0, encodi
     if skip_col != 0:
         try:
             line_found = line_found.split(col_sep)[skip_col]
-        except:
-            raise EOFError('error: value not found on this line: (' +
-                           line_found + ')!')
+        except Exception:
+            raise EOFError("error: value not found on this line: (" + line_found + ")!")
 
     # get the value
     result = get_real_from_line(line_found)
@@ -510,7 +532,15 @@ def get_line_col(filename, skip_line=0, skip_col=0, col_sep=None, seek=0, encodi
     return result
 
 
-def get_value(filename, token=None, skip_token=0, skip_line=0, skip_col=0, col_sep=None, encoding=default_encoding):
+def get_value(
+    filename,
+    token=None,
+    skip_token=0,
+    skip_line=0,
+    skip_col=0,
+    col_sep=None,
+    encoding=default_encoding,
+):
     """
     Get a value from a file using a delimiter and/or offsets.
 
@@ -590,28 +620,37 @@ def get_value(filename, token=None, skip_token=0, skip_line=0, skip_col=0, col_s
     1.3
     """
     if debug:
-        sys.stderr.write("get_value(token=" + str(token) +
-                         ', skip_token=' + str(skip_token) +
-                         ', skip_line=' + str(skip_line) +
-                         ', skip_col=' + str(skip_col) + ')\n')
+        sys.stderr.write(
+            "get_value(token="
+            + str(token)
+            + ", skip_token="
+            + str(skip_token)
+            + ", skip_line="
+            + str(skip_line)
+            + ", skip_col="
+            + str(skip_col)
+            + ")\n"
+        )
 
     # check parameters
     check_param(filename, str)
     if token is not None:
         check_param(token, str)
     elif skip_token != 0:
-        raise AssertionError("error: skip_token parameter needs token "
-                             "parameter to be set!")
+        raise AssertionError(
+            "error: skip_token parameter needs token " "parameter to be set!"
+        )
     check_param(skip_token, int)
     check_param(skip_line, int)
     check_param(skip_col, int)
 
     result = None
     if not token:
-        result = get_line_col(filename, skip_line, skip_col,
-                              col_sep=col_sep, encoding=encoding)
+        result = get_line_col(
+            filename, skip_line, skip_col, col_sep=col_sep, encoding=encoding
+        )
     else:
-        handle = open(filename, 'rb')
+        handle = open(filename, "rb")
 
         re_token = re.compile(token)
 
@@ -626,23 +665,25 @@ def get_value(filename, token=None, skip_token=0, skip_line=0, skip_col=0, col_s
             token_pos_cache = []
 
         line = handle.readline().decode(encoding)
-        while line != '':
+        while line != "":
             for token_match in re_token.finditer(line):
                 if skip_token > 0:
                     # token found but it's not the good one
                     skip_token -= 1
                 elif skip_token == 0:
                     # token found
-                    token_pos = [line_pos + token_match.start(),
-                                 line_pos + token_match.end()]
+                    token_pos = [
+                        line_pos + token_match.start(),
+                        line_pos + token_match.end(),
+                    ]
                     if debug:
-                        sys.stderr.write(
-                            'skip_token == 0, line: ' + line + '\n')
+                        sys.stderr.write("skip_token == 0, line: " + line + "\n")
                     break
                 else:
                     # token wanted in revert order: we first cache them all
-                    token_pos_cache.append([line_pos + token_match.start(),
-                                            line_pos + token_match.end()])
+                    token_pos_cache.append(
+                        [line_pos + token_match.start(), line_pos + token_match.end()]
+                    )
 
             if skip_token >= 0 and token_pos is not None:
                 # token found
@@ -657,14 +698,14 @@ def get_value(filename, token=None, skip_token=0, skip_line=0, skip_col=0, col_s
 
         if token_pos is None:
             handle.close()
-            raise EOFError('error: no token (' + token + ') was found!')
+            raise EOFError("error: no token (" + token + ") was found!")
 
         if skip_line == 0 and skip_col == 0:
             # get the real right after the token
             handle.seek(token_pos[1])
             line = handle.readline().decode(encoding)
             if debug:
-                sys.stderr.write('first token, line_found: ' + line + '\n')
+                sys.stderr.write("first token, line_found: " + line + "\n")
             result = get_real_from_line(line)
             handle.close()  # fixme: no multiple close?
         else:
@@ -678,12 +719,25 @@ def get_value(filename, token=None, skip_token=0, skip_line=0, skip_col=0, col_s
                 seek_pos = token_pos[1]
             handle.close()
             result = get_line_col(
-                filename, skip_line, skip_col, col_sep=col_sep, seek=seek_pos, encoding=encoding)
+                filename,
+                skip_line,
+                skip_col,
+                col_sep=col_sep,
+                seek=seek_pos,
+                encoding=encoding,
+            )
 
     return result
 
 
-def get(filename, tokens=None, skip_tokens=None, skip_lines=None, skip_cols=None, encoding=default_encoding):
+def get(
+    filename,
+    tokens=None,
+    skip_tokens=None,
+    skip_lines=None,
+    skip_cols=None,
+    encoding=default_encoding,
+):
     """
     Get several values from a file using delimiters and/or offsets.
 
@@ -772,11 +826,13 @@ def get(filename, tokens=None, skip_tokens=None, skip_lines=None, skip_cols=None
 
     # possible optimization: get values in one pass
     for n in range(nb_values):
-        results[n] = get_value(filename,
-                               tokens[n],
-                               skip_tokens[n],
-                               skip_lines[n],
-                               skip_cols[n],
-                               encoding=encoding)
+        results[n] = get_value(
+            filename,
+            tokens[n],
+            skip_tokens[n],
+            skip_lines[n],
+            skip_cols[n],
+            encoding=encoding,
+        )
 
     return results
diff --git a/python/src/openturns/PythonWrappingFunctions.hxx b/python/src/openturns/PythonWrappingFunctions.hxx
index 23c564010..d76d342f8 100644
--- a/python/src/openturns/PythonWrappingFunctions.hxx
+++ b/python/src/openturns/PythonWrappingFunctions.hxx
@@ -1847,7 +1847,7 @@ ScopedPyObjectPointer deepCopy(PyObject * pyObj)
   return pyObjDeepCopy;
 }
 
-inline void SignalHandler(int signum)
+inline void SignalHandler(int /*signum*/)
 {
   throw InterruptionException(HERE);
 }
diff --git a/python/test/CMakeLists.txt b/python/test/CMakeLists.txt
index 5c58519f7..ff0989021 100644
--- a/python/test/CMakeLists.txt
+++ b/python/test/CMakeLists.txt
@@ -384,6 +384,7 @@ ot_pyinstallcheck_test (Event_system IGNOREOUT)
 ot_pyinstallcheck_test (WhiteNoise_std)
 ot_pyinstallcheck_test (RandomWalk_std)
 ot_pyinstallcheck_test (DiscreteMarkovChain_std)
+ot_pyinstallcheck_test (AggregatedProcess_std)
 ot_pyinstallcheck_test (ARMAState_std)
 ot_pyinstallcheck_test (ARMACoefficients_std)
 ot_pyinstallcheck_test (ARMA_std)
diff --git a/python/test/t_AggregatedProcess_std.expout b/python/test/t_AggregatedProcess_std.expout
index e4d96692e..8dbc0fe64 100644
--- a/python/test/t_AggregatedProcess_std.expout
+++ b/python/test/t_AggregatedProcess_std.expout
@@ -1,17 +1,29 @@
-myAntecedentProcess =   ARMA process dimension = 1
-X_{0,t} = E_{0,t}
-with E_t ~ Normal(mu = 0, sigma = 1)
-myCompositeProcess =  class=CompositeProcess function=class=FieldFunction name=Unnamed implementation=class=ValueFunction evaluation=class=SymbolicEvaluation name=Unnamed inputVariablesNames=[x] outputVariablesNames=[y0] formulas=[2 * x + 5.0] antecedent=class= ARMA timeGrid=class=RegularGrid name=Unnamed start=0 step=0.1 n=11 coefficients AR=class=ARMACoefficients coefficients MA=class=ARMACoefficients noiseDistribution= class=Normal name=Normal dimension=1 mean=class=Point name=Unnamed dimension=1 values=[0] sigma=class=Point name=Unnamed dimension=1 values=[1] correlationMatrix=class=CorrelationMatrix dimension=1 implementation=class=MatrixImplementation name=Unnamed rows=1 columns=1 values=[1] state= class= ARMAState x= class=Sample name=Unnamed description=[] implementation=class=SampleImplementation name=Unnamed size=0 dimension=1 data=[] epsilon= class=Sample name=Unnamed description=[] implementation=class=SampleImplementation name=Unnamed size=0 dimension=1 data=[]
-One realization= 
-     [          y0       ]
- 0 : [ 0        0.63723  ]
- 1 : [ 0.1      5.70008  ]
- 2 : [ 0.2      4.28999  ]
- 3 : [ 0.3      7.8745   ]
- 4 : [ 0.4      6.62134  ]
- 5 : [ 0.5      6.58631  ]
- 6 : [ 0.6      4.05895  ]
- 7 : [ 0.7      5.52204  ]
- 8 : [ 0.8      0.419876 ]
- 9 : [ 0.9      2.43423  ]
-10 : [ 1        2.37644  ]
+myProcess= AggregatedProcess([WhiteNoise(Normal(mu = 0, sigma = 1))])
+a realization=      [ t         X0        ]
+ 0 : [  0         0.608202 ]
+ 1 : [  0.1      -1.26617  ]
+ 2 : [  0.2      -0.438266 ]
+ 3 : [  0.3       1.20548  ]
+ 4 : [  0.4      -2.18139  ]
+ 5 : [  0.5       0.350042 ]
+ 6 : [  0.6      -0.355007 ]
+ 7 : [  0.7       1.43725  ]
+ 8 : [  0.8       0.810668 ]
+ 9 : [  0.9       0.793156 ]
+10 : [  1        -0.470526 ]
+a marginal process= WhiteNoise(Normal(mu = 0, sigma = 1))
+myProcess= AggregatedProcess([WhiteNoise(Normal(mu = 0, sigma = 1)),ARMA(X_{0,t} = E_{0,t}, E_t ~ Normal(mu = 0, sigma = 1)),GaussianProcess(trend=[x0]->[0.0], covariance=ExponentialModel(scale=[1], amplitude=[1], no spatial correlation))])
+a realization=      [ t          X0         X0         y0         ]
+ 0 : [  0         -1.28289    0.351418  -0.0436123 ]
+ 1 : [  0.1       -1.31178    1.78236    0.190168  ]
+ 2 : [  0.2       -0.0907838  0.0702074  0.299777  ]
+ 3 : [  0.3        0.995793  -0.781366   0.444838  ]
+ 4 : [  0.4       -0.139453  -0.721533   0.195966  ]
+ 5 : [  0.5       -0.560206  -0.241223   0.0142559 ]
+ 6 : [  0.6        0.44549   -1.78796   -0.307618  ]
+ 7 : [  0.7        0.322925   0.40136   -0.16853   ]
+ 8 : [  0.8        0.445785   1.36783    0.685721  ]
+ 9 : [  0.9       -1.03808    1.00434    0.33466   ]
+10 : [  1         -0.856712   0.741548   1.09293   ]
+a marginal process= ARMA(X_{0,t} = E_{0,t}, E_t ~ Normal(mu = 0, sigma = 1))
+another marginal process= AggregatedProcess([WhiteNoise(Normal(mu = 0, sigma = 1)),GaussianProcess(trend=[x0]->[0.0], covariance=ExponentialModel(scale=[1], amplitude=[1], no spatial correlation))])
diff --git a/python/test/t_AggregatedProcess_std.py b/python/test/t_AggregatedProcess_std.py
index f22f80930..a66b1fccd 100755
--- a/python/test/t_AggregatedProcess_std.py
+++ b/python/test/t_AggregatedProcess_std.py
@@ -1,40 +1,35 @@
 #! /usr/bin/env python
 
-from openturns import *
-
-TESTPREAMBLE()
-RandomGenerator.SetSeed(0)
-
-try:
-    # Time grid creation
-    Tmin = 0.0
-    deltaT = 0.1
-    steps = 11
-
-    # Initialization of the time grid timeGrid1
-    timeGrid = RegularGrid(Tmin, deltaT, steps)
-
-    # Default aggregated process creation
-    myProcess = AggregatedProcess()
-
-    # We fix the time grid to the process
-    myProcess.setTimeGrid(timeGrid)
-
-    print("myProcess=", myProcess)
-    print("a realization=", myProcess.getRealization())
-    marginalProcess = myProcess.getMarginal(0)
-    print("a marginal process=", marginalProcess)
-
-    # Parameter aggregated process creation
-    coll = [WhiteNoise(), ARMA(), GaussianProcess()]
-    myProcess = AggregatedProcess(coll)
-    myProcess.setTimeGrid(timeGrid)
-    print("myProcess=", myProcess)
-    print("a realization=", myProcess.getRealization())
-    marginalProcess = myProcess.getMarginal(1)
-    print("a marginal process=", marginalProcess)
-    print("another marginal process=", myProcess.getMarginal([0, 2]))
-
-except:
-    import sys
-    print("t_AggregatedProcess_std.py", sys.exc_info()[0], sys.exc_info()[1])
+import openturns as ot
+
+ot.TESTPREAMBLE()
+
+# Time grid creation
+Tmin = 0.0
+deltaT = 0.1
+steps = 11
+
+# Initialization of the time grid timeGrid1
+timeGrid = ot.RegularGrid(Tmin, deltaT, steps)
+
+# Default aggregated process creation
+myProcess = ot.AggregatedProcess()
+
+# We fix the time grid to the process
+myProcess.setTimeGrid(timeGrid)
+
+print("myProcess=", myProcess)
+print("a realization=", myProcess.getRealization())
+marginalProcess = myProcess.getMarginal(0)
+print("a marginal process=", marginalProcess)
+
+# Parameter aggregated process creation
+coll = [ot.WhiteNoise(), ot.ARMA(), ot.GaussianProcess()]
+myProcess = ot.AggregatedProcess(coll)
+myProcess.setTimeGrid(timeGrid)
+print("myProcess=", myProcess)
+print("a realization=", myProcess.getRealization())
+marginalProcess = myProcess.getMarginal(1)
+print("a marginal process=", marginalProcess)
+print("another marginal process=", myProcess.getMarginal([0, 2]))
+
diff --git a/python/test/t_Distribution_scipy.expout b/python/test/t_Distribution_scipy.expout
index dda583aa2..a9a9770d0 100644
--- a/python/test/t_Distribution_scipy.expout
+++ b/python/test/t_Distribution_scipy.expout
@@ -1,30 +1,37 @@
 distribution= class=PythonDistribution name=SciPyDistribution
-realization= [0.37454]
-sample=     [ v0       ]
-0 : [ 0.950714 ]
-1 : [ 0.731994 ]
-2 : [ 0.598658 ]
-3 : [ 0.156019 ]
-4 : [ 0.155995 ]
-pdf= 1
-cdf= 0.6
-quantile= [0.6]
-quantile (tail)= [0.4]
-scalar quantile=0.6
-scalar quantile (tail)=0.4
-mean= [0.5]
-mean(sampling)= [0.494159]
-std= [0.288675]
-std(sampling)= [0.28763]
+continuous? True
+discrete? False
+integral? False
+realization= [0.49816]
+sample=     [ v0        ]
+0 : [  2.80286  ]
+1 : [  1.92798  ]
+2 : [  1.39463  ]
+3 : [ -0.375925 ]
+4 : [ -0.376022 ]
+pdf@0.1= 0.25
+pdf= 0.25
+cdf= 0.25
+quantile= [0]
+quantile (tail)= [2]
+scalar quantile=0
+scalar quantile (tail)=2
+mean= [1]
+mean(sampling)= [0.976638]
+std= [1.1547]
+std(sampling)= [1.15052]
 skewness= [0]
 skewness(sampling)= [0.0248807]
 kurtosis= [1.8]
 kurtosis(sampling)= [1.81085]
-range= [0, 1]
-parameter= []
-parameter= []
-parameterDesc= []
+range= [-1, 3]
+parameter= [-1,4]
+parameter= [3.5,4]
+parameterDesc= [parameter1,parameter2]
 distribution= class=PythonDistribution name=SciPyDistribution
+continuous? True
+discrete? False
+integral? False
 realization= [-1.65058]
 sample=     [ v0        ]
 0 : [ -0.409878 ]
@@ -32,12 +39,13 @@ sample=     [ v0        ]
 2 : [ -1.20987  ]
 3 : [ -2.3312   ]
 4 : [ -2.33131  ]
-pdf= 0.000502634
-cdf= 0.999936
-quantile= [0.6]
-quantile (tail)= [-8.49004]
-scalar quantile=0.6
-scalar quantile (tail)=-8.49004
+pdf@0.1= 0.0190201
+pdf= 0.034761
+cdf= 0.994614
+quantile= [-7.89492e-16]
+quantile (tail)= [-4.77177]
+scalar quantile=-7.89492e-16
+scalar quantile (tail)=-4.77177
 mean= [-1.53652]
 mean(sampling)= [-1.55154]
 std= [0.870537]
@@ -50,3 +58,88 @@ range= ]-inf (-46.5418), (4.77335) +inf[
 parameter= [2.55,2.25]
 parameter= [3.5,2.25]
 parameterDesc= [parameter1,parameter2]
+distribution= class=PythonDistribution name=SciPyDistribution
+continuous? False
+discrete? True
+integral? True
+realization= [4]
+sample=     [ v0 ]
+0 : [ 8  ]
+1 : [ 6  ]
+2 : [ 5  ]
+3 : [ 3  ]
+4 : [ 3  ]
+pdf@0.1= 0
+pdf= 0.000976562
+cdf= 0.000976562
+quantile= [0]
+quantile (tail)= [9]
+scalar quantile=0
+scalar quantile (tail)=9
+mean= [5]
+mean(sampling)= [4.9684]
+std= [1.58114]
+std(sampling)= [1.57081]
+skewness= [0]
+skewness(sampling)= [-0.00105239]
+kurtosis= [2.8]
+kurtosis(sampling)= [2.8062]
+range= [0, 10]
+support=  0 : [  0 ]
+ 1 : [  1 ]
+ 2 : [  2 ]
+ 3 : [  3 ]
+ 4 : [  4 ]
+ 5 : [  5 ]
+ 6 : [  6 ]
+ 7 : [  7 ]
+ 8 : [  8 ]
+ 9 : [  9 ]
+10 : [ 10 ]
+parameter= [10,0.5]
+parameter= [3.5,0.5]
+parameterDesc= [parameter1,parameter2]
+distribution= class=PythonDistribution name=SciPyDistribution
+continuous? False
+discrete? True
+integral? True
+realization= [0]
+sample=     [ v0 ]
+0 : [ 2  ]
+1 : [ 0  ]
+2 : [ 0  ]
+3 : [ 0  ]
+4 : [ 1  ]
+pdf@0.1= 0
+pdf= 0.548812
+cdf= 0.548812
+quantile= [0]
+quantile (tail)= [0]
+scalar quantile=0
+scalar quantile (tail)=0
+mean= [0.6]
+mean(sampling)= [0.5893]
+std= [0.774597]
+std(sampling)= [0.765169]
+skewness= [1.29099]
+skewness(sampling)= [1.29299]
+kurtosis= [4.66667]
+kurtosis(sampling)= [4.62228]
+range= [0, (13) +inf[
+support=  0 : [  0 ]
+ 1 : [  1 ]
+ 2 : [  2 ]
+ 3 : [  3 ]
+ 4 : [  4 ]
+ 5 : [  5 ]
+ 6 : [  6 ]
+ 7 : [  7 ]
+ 8 : [  8 ]
+ 9 : [  9 ]
+10 : [ 10 ]
+11 : [ 11 ]
+12 : [ 12 ]
+13 : [ 13 ]
+parameter= [0.6]
+parameter= [3.5]
+parameterDesc= [parameter1]
diff --git a/python/test/t_Distribution_scipy.py b/python/test/t_Distribution_scipy.py
index 7679f0bf4..7c33bc89b 100755
--- a/python/test/t_Distribution_scipy.py
+++ b/python/test/t_Distribution_scipy.py
@@ -4,7 +4,8 @@ import openturns as ot
 import scipy.stats as st
 import numpy as np
 
-for scipy_dist in [st.uniform(), st.johnsonsu(2.55, 2.25)]:
+for scipy_dist in [st.uniform(-1.0, 4.0), st.johnsonsu(2.55, 2.25),
+                   st.binom(10, 0.5), st.poisson(0.6)]:
     np.random.seed(42)
 
     # create an openturns distribution
@@ -12,10 +13,14 @@ for scipy_dist in [st.uniform(), st.johnsonsu(2.55, 2.25)]:
     distribution = ot.Distribution(py_dist)
 
     print('distribution=', distribution)
+    print('continuous?', distribution.isContinuous())
+    print('discrete?', distribution.isDiscrete())
+    print('integral?', distribution.isIntegral())
     print('realization=', distribution.getRealization())
     sample = distribution.getSample(10000)
     print('sample=', sample[0:5])
-    point = [0.6]
+    print('pdf@0.1= %.6g' % (distribution.computePDF([0.1])))
+    point = [0.0]
     print('pdf= %.6g' % distribution.computePDF(point))
     cdf = distribution.computeCDF(point)
     print('cdf= %.6g' % cdf)
@@ -33,6 +38,8 @@ for scipy_dist in [st.uniform(), st.johnsonsu(2.55, 2.25)]:
     print('kurtosis=', distribution.getKurtosis())
     print('kurtosis(sampling)=', sample.computeKurtosis())
     print('range=', distribution.getRange())
+    if distribution.isDiscrete():
+        print('support=', distribution.getSupport())
     parameter = distribution.getParameter()
     print('parameter=', distribution.getParameter())
     if len(parameter) > 0:
diff --git a/python/test/t_FunctionalChaosSobolIndices_std.py b/python/test/t_FunctionalChaosSobolIndices_std.py
index 7e42ef54b..8ef99a095 100755
--- a/python/test/t_FunctionalChaosSobolIndices_std.py
+++ b/python/test/t_FunctionalChaosSobolIndices_std.py
@@ -46,4 +46,4 @@ algo.run()
 result = algo.getResult()
 sensitivity = ot.FunctionalChaosSobolIndices(result)
 
-print(sensitivity.summary())
+print(sensitivity)
diff --git a/python/test/t_LinearModelStepwiseAlgorithm_std.py b/python/test/t_LinearModelStepwiseAlgorithm_std.py
index 2a85b670e..9d6db6f8f 100755
--- a/python/test/t_LinearModelStepwiseAlgorithm_std.py
+++ b/python/test/t_LinearModelStepwiseAlgorithm_std.py
@@ -113,3 +113,21 @@ for k in [penalty_AIC, penalty_BIC]:
             print(" Both " + IC)
         print("{0:~^60s}".format(""))
         print(algo_result)
+
+dist = ot.ComposedDistribution([ot.Uniform(0., 10.)]*2)
+f = ot.SymbolicFunction(['x1', 'x2'], ["0.5*x1 + x2+ 7"])
+X = ot.LHSExperiment(dist, 10).generate()
+X = ot.MonteCarloExperiment(dist, 10).generate()
+Y = f(X)
+i_min = [0]
+coll = []
+coll.append(ot.SymbolicFunction(['x1', 'x2'], ["x1"]))
+coll.append(ot.SymbolicFunction(['x1', 'x2'], ["0"]))
+coll.append(ot.SymbolicFunction(['x1', 'x2'], ["x2"]))
+basis = ot.Basis(coll)
+algo = ot.LinearModelStepwiseAlgorithm(X, basis, Y, i_min)
+try:
+    algo.run()
+    assert False, "should not go here"
+except:
+    pass
diff --git a/python/test/t_ProcessSample_std.py b/python/test/t_ProcessSample_std.py
index 6786b0c5f..3af13c117 100755
--- a/python/test/t_ProcessSample_std.py
+++ b/python/test/t_ProcessSample_std.py
@@ -127,3 +127,10 @@ print('quantile at level', q, '=',
 print('empirical CDF at', x, '=', processSample.computeEmpiricalCDF(x).getValues())
 print('complementary empirical CDF at', x, '=',
       processSample.computeEmpiricalCDF(x, True).getValues())
+
+n = mesh.getVerticesNumber()
+try:
+    s = processSample.getSampleAtVertex(n)
+    assert False, "should not go there"
+except:
+    pass
diff --git a/python/test/t_Waarts_noisy_lsf.py b/python/test/t_Waarts_noisy_lsf.py
index 52ce49225..bc37b9b41 100755
--- a/python/test/t_Waarts_noisy_lsf.py
+++ b/python/test/t_Waarts_noisy_lsf.py
@@ -29,7 +29,7 @@ try:
 
     EtatLimite = Function(inputFunction, outputFunction, formulas)
 
-    d * dimim = EtatLimite.getInputDimension()
+    dim = EtatLimite.getInputDimension()
     print(dim)
 
     #
['lib/src/Base/Algo/EnclosingSimplexAlgorithmImplementation.cxx', 'python/src/Wishart_doc.i.in', 'python/src/KernelSmoothing_doc.i.in', 'lib/src/Uncertainty/Algorithm/MetaModel/LinearModel/LinearModelStepwiseAlgorithm.cxx', '.github/workflows/build.yml', 'python/test/t_AggregatedProcess_std.py', 'python/test/t_AggregatedProcess_std.expout', 'python/doc/examples/numerical_methods/optimization/plot_optimization_bonmin.py', 'python/test/t_FunctionalChaosSobolIndices_std.py', 'python/src/InverseWishart_doc.i.in', 'lib/src/Uncertainty/Algorithm/Sensitivity/SobolIndicesAlgorithmImplementation.cxx', 'python/test/t_Waarts_noisy_lsf.py', 'TODO', 'python/doc/developer_guide/sphinx_documentation.rst', 'python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_ishigami_grouped_indices.py', 'python/test/t_ProcessSample_std.py', 'python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_chaos_ishigami.py', 'python/doc/examples/data_analysis/distribution_fitting/plot_estimate_conditional_quantile.py', 'lib/src/Uncertainty/Algorithm/MetaModel/FunctionalChaos/openturns/FunctionalChaosSobolIndices.hxx', 'python/doc/examples/reliability_sensitivity/sensitivity_analysis/plot_functional_chaos_sensitivity.py', 'python/src/PythonDistribution.cxx', 'python/src/openturns/PythonWrappingFunctions.hxx', 'python/test/t_Distribution_scipy.py', 'python/src/Distribution.i', 'python/test/t_Distribution_scipy.expout', 'python/doc/examples/meta_modeling/polynomial_chaos_metamodel/plot_functional_chaos.py', 'lib/src/Uncertainty/Algorithm/MetaModel/FunctionalChaos/FunctionalChaosSobolIndices.cxx', 'ChangeLog', 'python/src/coupling_tools.py', 'python/src/FunctionalChaosSobolIndices_doc.i.in', 'lib/src/Base/Stat/ProcessSampleImplementation.cxx', 'python/test/CMakeLists.txt', 'python/doc/about.rst', 'python/test/t_LinearModelStepwiseAlgorithm_std.py']