From 81753eaf4004a62d276e18ba613189276d575658 Mon Sep 17 00:00:00 2001
From: Manuel Guenther <manuel.guenther@idiap.ch>
Date: Mon, 24 Mar 2014 16:36:23 +0100
Subject: [PATCH] Made documentation more beautiful using the novel
 xbob.extension/documentation.h

---
 buildout.cfg                    |   1 +
 doc/conf.py                     |   9 +-
 doc/index.rst                   |   5 +
 doc/py_api.rst                  |  12 +
 setup.py                        |  10 +-
 xbob/math/__init__.py           |   6 +-
 xbob/math/histogram.cpp         | 154 +++----
 xbob/math/lp_interior_point.cpp | 501 ++++++++++-----------
 xbob/math/main.cpp              | 774 +++++++++++++++-----------------
 9 files changed, 723 insertions(+), 749 deletions(-)

diff --git a/buildout.cfg b/buildout.cfg
index 3ae5d14..8afd587 100644
--- a/buildout.cfg
+++ b/buildout.cfg
@@ -21,6 +21,7 @@ prefixes = /idiap/group/torch5spro/nightlies/last/bob/linux-x86_64-release
 [sources]
 xbob.extension = git https://github.com/bioidiap/xbob.extension branch=prototype
 xbob.blitz = git https://github.com/bioidiap/xbob.blitz
+numpydoc = git git@github.com:numpy/numpydoc.git
 
 [scripts]
 recipe = xbob.buildout:scripts
diff --git a/doc/conf.py b/doc/conf.py
index 22d011f..e735999 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -31,8 +31,14 @@ extensions = [
   'sphinx.ext.autosummary',
   'sphinx.ext.doctest',
   'sphinx.ext.intersphinx',
+  'numpydoc',
   ]
 
+# See: https://github.com/phn/pytpm/issues/3
+numpydoc_show_class_members = False
+# See... figured out myself :-(
+numpydoc_class_members_toctree = False
+
 # The viewcode extension appeared only on Sphinx >= 1.0.0
 import sphinx
 if sphinx.__version__ >= "1.0":
@@ -143,7 +149,7 @@ html_favicon = ''
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
-html_static_path = ['_static']
+# html_static_path = ['_static']
 
 # If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
 # using the given strftime format.
@@ -250,6 +256,7 @@ autodoc_member_order = 'bysource'
 autodoc_default_flags = [
   'members',
   'undoc-members',
+  'private-members', 
   'inherited-members',
   'show-inheritance',
   ]
diff --git a/doc/index.rst b/doc/index.rst
index 4fcf242..d989a58 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -28,4 +28,9 @@ Indices and tables
 * :ref:`modindex`
 * :ref:`search`
 
+To-Do list
+----------
+
+.. todolist::
+
 .. include:: links.rst
diff --git a/doc/py_api.rst b/doc/py_api.rst
index 627617b..9e2879b 100644
--- a/doc/py_api.rst
+++ b/doc/py_api.rst
@@ -8,5 +8,17 @@
 
 This section includes information for using the pure Python API of ``xbob.math``.
 
+Summary
+.......
+
+.. autosummary::
+  xbob.math.LPInteriorPoint
+  xbob.math.LPInteriorPointShortstep
+  xbob.math.LPInteriorPointLongstep
+
+
+Details
+.......
+
 .. automodule:: xbob.math
    :members:
diff --git a/setup.py b/setup.py
index 79ed6cf..b5772d8 100644
--- a/setup.py
+++ b/setup.py
@@ -4,8 +4,9 @@
 # Mon 16 Apr 08:18:08 2012 CEST
 
 from setuptools import setup, find_packages, dist
-dist.Distribution(dict(setup_requires=['xbob.blitz']))
+dist.Distribution(dict(setup_requires=['xbob.blitz','xbob.extension']))
 from xbob.blitz.extension import Extension
+import xbob.extension
 
 import os
 package_dir = os.path.dirname(os.path.realpath(__file__))
@@ -33,6 +34,8 @@ setup(
     install_requires=[
       'setuptools',
       'xbob.blitz',
+      'xbob.extension',
+      'numpydoc',
     ],
 
     namespace_packages=[
@@ -61,8 +64,9 @@ setup(
         packages = packages,
         version = version,
         include_dirs = include_dirs,
-        ),
-      ],
+#        define_macros = [('XBOB_SHORT_DOCSTRINGS',1)],
+      ),
+    ],
 
     entry_points={
       },
diff --git a/xbob/math/__init__.py b/xbob/math/__init__.py
index 1975a66..ea61aa6 100644
--- a/xbob/math/__init__.py
+++ b/xbob/math/__init__.py
@@ -1,6 +1,10 @@
 from ._library import __version__
 from ._library import *
 
+# To get the full API documentation automatically
+__all__ = dir()
+
+
 def get_config():
   """Returns a string containing the configuration information.
   """
@@ -21,4 +25,4 @@ def get_config():
   return retval.strip()
 
 # gets sphinx autodoc done right - don't remove it
-__all__ = [_ for _ in dir() if not _.startswith('_')]
+__all__ = [_ for _ in dir() if not _.startswith('_')]
\ No newline at end of file
diff --git a/xbob/math/histogram.cpp b/xbob/math/histogram.cpp
index 6126e1e..0fc0ceb 100644
--- a/xbob/math/histogram.cpp
+++ b/xbob/math/histogram.cpp
@@ -1,7 +1,7 @@
 /**
  * @author Manuel Guenther <Manuel.Guenther@idiap.ch>
  * @author Andre Anjos <andre.anjos@idiap.ch>
- * @date Tue  3 Dec 14:23:42 2013 CET 
+ * @date Tue  3 Dec 14:23:42 2013 CET
  *
  * @brief Binds fast versions of some histogram measures
  *
@@ -89,45 +89,45 @@ static PyObject* py_histogram_intersection_1
 }
 
 template <typename T1> PyObject* py_histogram_intersection_2_inner(
-    PyBlitzArrayObject* index1, PyBlitzArrayObject* value1, 
+    PyBlitzArrayObject* index1, PyBlitzArrayObject* value1,
     PyBlitzArrayObject* index2, PyBlitzArrayObject* value2) {
 
   switch(value1->type_num) {
 
     case NPY_UINT8:
       return PyBlitzArrayCxx_FromCScalar(bob::math::histogram_intersection(
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1), 
-            *PyBlitzArrayCxx_AsBlitz<uint8_t,1>(value1), 
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2), 
-            *PyBlitzArrayCxx_AsBlitz<uint8_t,1>(value2))); 
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1),
+            *PyBlitzArrayCxx_AsBlitz<uint8_t,1>(value1),
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2),
+            *PyBlitzArrayCxx_AsBlitz<uint8_t,1>(value2)));
 
     case NPY_UINT16:
       return PyBlitzArrayCxx_FromCScalar(bob::math::histogram_intersection(
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1), 
-            *PyBlitzArrayCxx_AsBlitz<uint16_t,1>(value1), 
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2), 
-            *PyBlitzArrayCxx_AsBlitz<uint16_t,1>(value2))); 
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1),
+            *PyBlitzArrayCxx_AsBlitz<uint16_t,1>(value1),
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2),
+            *PyBlitzArrayCxx_AsBlitz<uint16_t,1>(value2)));
 
     case NPY_INT32:
       return PyBlitzArrayCxx_FromCScalar(bob::math::histogram_intersection(
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1), 
-            *PyBlitzArrayCxx_AsBlitz<int32_t,1>(value1), 
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2), 
-            *PyBlitzArrayCxx_AsBlitz<int32_t,1>(value2))); 
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1),
+            *PyBlitzArrayCxx_AsBlitz<int32_t,1>(value1),
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2),
+            *PyBlitzArrayCxx_AsBlitz<int32_t,1>(value2)));
 
     case NPY_INT64:
       return PyBlitzArrayCxx_FromCScalar(bob::math::histogram_intersection(
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1), 
-            *PyBlitzArrayCxx_AsBlitz<int64_t,1>(value1), 
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2), 
-            *PyBlitzArrayCxx_AsBlitz<int64_t,1>(value2))); 
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1),
+            *PyBlitzArrayCxx_AsBlitz<int64_t,1>(value1),
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2),
+            *PyBlitzArrayCxx_AsBlitz<int64_t,1>(value2)));
 
     case NPY_FLOAT64:
       return PyBlitzArrayCxx_FromCScalar(bob::math::histogram_intersection(
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1), 
-            *PyBlitzArrayCxx_AsBlitz<double,1>(value1), 
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2), 
-            *PyBlitzArrayCxx_AsBlitz<double,1>(value2))); 
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1),
+            *PyBlitzArrayCxx_AsBlitz<double,1>(value1),
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2),
+            *PyBlitzArrayCxx_AsBlitz<double,1>(value2)));
 
     default:
       break;
@@ -151,10 +151,10 @@ static PyObject* py_histogram_intersection_2(PyObject*, PyObject* args, PyObject
   PyBlitzArrayObject* value2 = 0;
 
   if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&O&O&O&",
-        kwlist, 
-        &PyBlitzArray_Converter, &index1, 
+        kwlist,
+        &PyBlitzArray_Converter, &index1,
         &PyBlitzArray_Converter, &value1,
-        &PyBlitzArray_Converter, &index2, 
+        &PyBlitzArray_Converter, &index2,
         &PyBlitzArray_Converter, &value2
         )) return 0;
 
@@ -176,7 +176,7 @@ static PyObject* py_histogram_intersection_2(PyObject*, PyObject* args, PyObject
   }
 
   // input arrays must be 1d
-  if (index1->ndim != 1 || index2->ndim != 1 || 
+  if (index1->ndim != 1 || index2->ndim != 1 ||
       value1->ndim != 1 || value2->ndim != 1) {
     PyErr_SetString(PyExc_TypeError, "all input arrays must be 1D");
     return 0;
@@ -237,7 +237,7 @@ static PyObject* py_histogram_intersection_2(PyObject*, PyObject* args, PyObject
  * Note: Dispatcher function.
  */
 PyObject* py_histogram_intersection (PyObject*, PyObject* args, PyObject* kwargs) {
-  
+
   Py_ssize_t nargs = args?PyTuple_Size(args):0 + kwargs?PyDict_Size(kwargs):0;
 
   PyObject* retval = 0;
@@ -325,45 +325,45 @@ static PyObject* py_chi_square_1
 }
 
 template <typename T1> PyObject* py_chi_square_2_inner(
-    PyBlitzArrayObject* index1, PyBlitzArrayObject* value1, 
+    PyBlitzArrayObject* index1, PyBlitzArrayObject* value1,
     PyBlitzArrayObject* index2, PyBlitzArrayObject* value2) {
 
   switch(value1->type_num) {
 
     case NPY_UINT8:
       return PyBlitzArrayCxx_FromCScalar(bob::math::chi_square(
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1), 
-            *PyBlitzArrayCxx_AsBlitz<uint8_t,1>(value1), 
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2), 
-            *PyBlitzArrayCxx_AsBlitz<uint8_t,1>(value2))); 
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1),
+            *PyBlitzArrayCxx_AsBlitz<uint8_t,1>(value1),
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2),
+            *PyBlitzArrayCxx_AsBlitz<uint8_t,1>(value2)));
 
     case NPY_UINT16:
       return PyBlitzArrayCxx_FromCScalar(bob::math::chi_square(
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1), 
-            *PyBlitzArrayCxx_AsBlitz<uint16_t,1>(value1), 
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2), 
-            *PyBlitzArrayCxx_AsBlitz<uint16_t,1>(value2))); 
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1),
+            *PyBlitzArrayCxx_AsBlitz<uint16_t,1>(value1),
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2),
+            *PyBlitzArrayCxx_AsBlitz<uint16_t,1>(value2)));
 
     case NPY_INT32:
       return PyBlitzArrayCxx_FromCScalar(bob::math::chi_square(
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1), 
-            *PyBlitzArrayCxx_AsBlitz<int32_t,1>(value1), 
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2), 
-            *PyBlitzArrayCxx_AsBlitz<int32_t,1>(value2))); 
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1),
+            *PyBlitzArrayCxx_AsBlitz<int32_t,1>(value1),
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2),
+            *PyBlitzArrayCxx_AsBlitz<int32_t,1>(value2)));
 
     case NPY_INT64:
       return PyBlitzArrayCxx_FromCScalar(bob::math::chi_square(
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1), 
-            *PyBlitzArrayCxx_AsBlitz<int64_t,1>(value1), 
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2), 
-            *PyBlitzArrayCxx_AsBlitz<int64_t,1>(value2))); 
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1),
+            *PyBlitzArrayCxx_AsBlitz<int64_t,1>(value1),
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2),
+            *PyBlitzArrayCxx_AsBlitz<int64_t,1>(value2)));
 
     case NPY_FLOAT64:
       return PyBlitzArrayCxx_FromCScalar(bob::math::chi_square(
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1), 
-            *PyBlitzArrayCxx_AsBlitz<double,1>(value1), 
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2), 
-            *PyBlitzArrayCxx_AsBlitz<double,1>(value2))); 
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1),
+            *PyBlitzArrayCxx_AsBlitz<double,1>(value1),
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2),
+            *PyBlitzArrayCxx_AsBlitz<double,1>(value2)));
 
     default:
       break;
@@ -387,10 +387,10 @@ static PyObject* py_chi_square_2(PyObject*, PyObject* args, PyObject* kwds) {
   PyBlitzArrayObject* value2 = 0;
 
   if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&O&O&O&",
-        kwlist, 
-        &PyBlitzArray_Converter, &index1, 
+        kwlist,
+        &PyBlitzArray_Converter, &index1,
         &PyBlitzArray_Converter, &value1,
-        &PyBlitzArray_Converter, &index2, 
+        &PyBlitzArray_Converter, &index2,
         &PyBlitzArray_Converter, &value2
         )) return 0;
 
@@ -412,7 +412,7 @@ static PyObject* py_chi_square_2(PyObject*, PyObject* args, PyObject* kwds) {
   }
 
   // input arrays must be 1d
-  if (index1->ndim != 1 || index2->ndim != 1 || 
+  if (index1->ndim != 1 || index2->ndim != 1 ||
       value1->ndim != 1 || value2->ndim != 1) {
     PyErr_SetString(PyExc_TypeError, "all input arrays must be 1D");
     return 0;
@@ -549,45 +549,45 @@ static PyObject* py_kullback_leibler_1
 }
 
 template <typename T1> PyObject* py_kullback_leibler_2_inner(
-    PyBlitzArrayObject* index1, PyBlitzArrayObject* value1, 
+    PyBlitzArrayObject* index1, PyBlitzArrayObject* value1,
     PyBlitzArrayObject* index2, PyBlitzArrayObject* value2) {
 
   switch(value1->type_num) {
 
     case NPY_UINT8:
       return PyBlitzArrayCxx_FromCScalar(bob::math::kullback_leibler(
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1), 
-            *PyBlitzArrayCxx_AsBlitz<uint8_t,1>(value1), 
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2), 
-            *PyBlitzArrayCxx_AsBlitz<uint8_t,1>(value2))); 
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1),
+            *PyBlitzArrayCxx_AsBlitz<uint8_t,1>(value1),
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2),
+            *PyBlitzArrayCxx_AsBlitz<uint8_t,1>(value2)));
 
     case NPY_UINT16:
       return PyBlitzArrayCxx_FromCScalar(bob::math::kullback_leibler(
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1), 
-            *PyBlitzArrayCxx_AsBlitz<uint16_t,1>(value1), 
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2), 
-            *PyBlitzArrayCxx_AsBlitz<uint16_t,1>(value2))); 
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1),
+            *PyBlitzArrayCxx_AsBlitz<uint16_t,1>(value1),
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2),
+            *PyBlitzArrayCxx_AsBlitz<uint16_t,1>(value2)));
 
     case NPY_INT32:
       return PyBlitzArrayCxx_FromCScalar(bob::math::kullback_leibler(
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1), 
-            *PyBlitzArrayCxx_AsBlitz<int32_t,1>(value1), 
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2), 
-            *PyBlitzArrayCxx_AsBlitz<int32_t,1>(value2))); 
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1),
+            *PyBlitzArrayCxx_AsBlitz<int32_t,1>(value1),
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2),
+            *PyBlitzArrayCxx_AsBlitz<int32_t,1>(value2)));
 
     case NPY_INT64:
       return PyBlitzArrayCxx_FromCScalar(bob::math::kullback_leibler(
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1), 
-            *PyBlitzArrayCxx_AsBlitz<int64_t,1>(value1), 
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2), 
-            *PyBlitzArrayCxx_AsBlitz<int64_t,1>(value2))); 
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1),
+            *PyBlitzArrayCxx_AsBlitz<int64_t,1>(value1),
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2),
+            *PyBlitzArrayCxx_AsBlitz<int64_t,1>(value2)));
 
     case NPY_FLOAT64:
       return PyBlitzArrayCxx_FromCScalar(bob::math::kullback_leibler(
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1), 
-            *PyBlitzArrayCxx_AsBlitz<double,1>(value1), 
-            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2), 
-            *PyBlitzArrayCxx_AsBlitz<double,1>(value2))); 
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index1),
+            *PyBlitzArrayCxx_AsBlitz<double,1>(value1),
+            *PyBlitzArrayCxx_AsBlitz<T1,1>(index2),
+            *PyBlitzArrayCxx_AsBlitz<double,1>(value2)));
 
     default:
       break;
@@ -611,10 +611,10 @@ static PyObject* py_kullback_leibler_2(PyObject*, PyObject* args, PyObject* kwds
   PyBlitzArrayObject* value2 = 0;
 
   if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&O&O&O&",
-        kwlist, 
-        &PyBlitzArray_Converter, &index1, 
+        kwlist,
+        &PyBlitzArray_Converter, &index1,
         &PyBlitzArray_Converter, &value1,
-        &PyBlitzArray_Converter, &index2, 
+        &PyBlitzArray_Converter, &index2,
         &PyBlitzArray_Converter, &value2
         )) return 0;
 
@@ -636,7 +636,7 @@ static PyObject* py_kullback_leibler_2(PyObject*, PyObject* args, PyObject* kwds
   }
 
   // input arrays must be 1d
-  if (index1->ndim != 1 || index2->ndim != 1 || 
+  if (index1->ndim != 1 || index2->ndim != 1 ||
       value1->ndim != 1 || value2->ndim != 1) {
     PyErr_SetString(PyExc_TypeError, "all input arrays must be 1D");
     return 0;
diff --git a/xbob/math/lp_interior_point.cpp b/xbob/math/lp_interior_point.cpp
index f50df2b..a99ac4a 100644
--- a/xbob/math/lp_interior_point.cpp
+++ b/xbob/math/lp_interior_point.cpp
@@ -15,32 +15,22 @@
 #include <bob/math/LPInteriorPoint.h>
 #include <structmember.h>
 
+#include <xbob.extension/documentation.h>
+
 /************************************************
  * Implementation of LPInteriorPoint base class *
  ************************************************/
 
-PyDoc_STRVAR(s_lpinteriorpoint_str, XBOB_EXT_MODULE_PREFIX ".LPInteriorPoint");
-
-PyDoc_STRVAR(s_lpinteriorpoint_doc,
-"Base class to solve a linear program using interior point methods.\n\
-For more details about the algorithms,please refer to the following\n\
-book: *\"Primal-Dual Interior-Point Methods\", Stephen J. Wright,\n\
-ISBN: 978-0898713824, Chapter 5, \"Path-Following Algorithms\"*.\n\
-\n\
-.. warning::\n\
-\n\
-   You cannot instantiate an object of this type directly, you must\n\
-   use it through one of the inherited types.\n\
-\n\
-The primal linear program (LP) is defined as follows:\n\
-\n\
-   min transpose(c)*x, s.t. A*x=b, x>=0\n\
-\n\
-The dual formulation is:\n\
-\n\
-   min transpose(b)*lambda, s.t. transpose(A)*lambda+mu=c\n\
-\n\
-");
+static auto s_lpinteriorpoint = xbob::extension::ClassDoc(
+  XBOB_EXT_MODULE_PREFIX ".LPInteriorPoint",
+  "Base class to solve a linear program using interior point methods.",
+  "For more details about the algorithms,please refer to the following book: *'Primal-Dual Interior-Point Methods', Stephen J. Wright, ISBN: 978-0898713824, Chapter 5, 'Path-Following Algorithms'*. "
+  ".. warning:: You cannot instantiate an object of this type directly, you must use it through one of the inherited types.\n"
+  "The primal linear program (LP) is defined as follows:\n"
+  ".. math:: \\min c^T*x \\text{, s.t. } A*x=b, x>=0\n"
+  "The dual formulation is:\n"
+  ".. math:: \\min b^T*\\lambda \\text{, s.t. } A^T*\\lambda+\\mu=c"
+);
 
 /* Type definition for PyBobMathLpInteriorPointObject */
 typedef struct {
@@ -59,11 +49,13 @@ static int PyBobMathLpInteriorPoint_init(PyBobMathLpInteriorPointObject* self, P
 
 }
 
-PyDoc_STRVAR(s_M_str, "m");
-PyDoc_STRVAR(s_M_doc,
-"The first dimension of the problem/A matrix"
+static auto s_M = xbob::extension::VariableDoc(
+  "m",
+  "int",
+  "The first dimension of the problem/A matrix"
 );
 
+
 static PyObject* PyBobMathLpInteriorPoint_getM (PyBobMathLpInteriorPointObject* self,
     void* /*closure*/) {
   return Py_BuildValue("n", self->base->getDimM());
@@ -91,9 +83,10 @@ static int PyBobMathLpInteriorPoint_setM (PyBobMathLpInteriorPointObject* self,
 
 }
 
-PyDoc_STRVAR(s_N_str, "n");
-PyDoc_STRVAR(s_N_doc,
-"The second dimension of the problem/A matrix"
+static auto s_N = xbob::extension::VariableDoc(
+  "n",
+  "int",
+  "The second dimension of the problem/A matrix"
 );
 
 static PyObject* PyBobMathLpInteriorPoint_getN (PyBobMathLpInteriorPointObject* self,
@@ -123,9 +116,10 @@ static int PyBobMathLpInteriorPoint_setN (PyBobMathLpInteriorPointObject* self,
 
 }
 
-PyDoc_STRVAR(s_epsilon_str, "epsilon");
-PyDoc_STRVAR(s_epsilon_doc,
-"The precision to determine whether an equality constraint is fulfilled or not"
+static auto s_epsilon = xbob::extension::VariableDoc(
+  "epsilon",
+  "float",
+  "The precision to determine whether an equality constraint is fulfilled or not"
 );
 
 static PyObject* PyBobMathLpInteriorPoint_getEpsilon (PyBobMathLpInteriorPointObject* self, void* /*closure*/) {
@@ -154,9 +148,10 @@ static int PyBobMathLpInteriorPoint_setEpsilon (PyBobMathLpInteriorPointObject*
 
 }
 
-PyDoc_STRVAR(s_lambda_str, "lambda_");
-PyDoc_STRVAR(s_lambda_doc,
-"The value of the lambda dual variable (read-only)"
+static auto s_lambda = xbob::extension::VariableDoc(
+  "lambda_",
+  "float",
+  "The value of the :math:`\\lambda` dual variable (read-only)"
 );
 
 static PyObject* PyBobMathLpInteriorPoint_lambda (PyBobMathLpInteriorPointObject* self) {
@@ -171,9 +166,10 @@ static PyObject* PyBobMathLpInteriorPoint_lambda (PyBobMathLpInteriorPointObject
   return PyBlitzArray_NUMPY_WRAP(retval);
 }
 
-PyDoc_STRVAR(s_mu_str, "mu");
-PyDoc_STRVAR(s_mu_doc,
-"The value of the mu dual variable (read-only)"
+static auto s_mu = xbob::extension::VariableDoc(
+  "mu",
+  "float",
+  "The value of the :math:`\\mu` dual variable (read-only)"
 );
 
 static PyObject* PyBobMathLpInteriorPoint_mu (PyBobMathLpInteriorPointObject* self) {
@@ -190,49 +186,51 @@ static PyObject* PyBobMathLpInteriorPoint_mu (PyBobMathLpInteriorPointObject* se
 
 static PyGetSetDef PyBobMathLpInteriorPoint_getseters[] = {
     {
-      s_M_str,
+      s_M.name(),
       (getter)PyBobMathLpInteriorPoint_getM,
       (setter)PyBobMathLpInteriorPoint_setM,
-      s_M_doc,
+      s_M.doc(),
       0
     },
     {
-      s_N_str,
+      s_N.name(),
       (getter)PyBobMathLpInteriorPoint_getN,
       (setter)PyBobMathLpInteriorPoint_setN,
-      s_N_doc,
+      s_N.doc(),
       0
     },
     {
-      s_epsilon_str,
+      s_epsilon.name(),
       (getter)PyBobMathLpInteriorPoint_getEpsilon,
       (setter)PyBobMathLpInteriorPoint_setEpsilon,
-      s_epsilon_doc,
+      s_epsilon.doc(),
       0
     },
     {
-      s_lambda_str,
+      s_lambda.name(),
       (getter)PyBobMathLpInteriorPoint_lambda,
       0,
-      s_lambda_doc,
+      s_lambda.doc(),
       0
     },
     {
-      s_mu_str,
+      s_mu.name(),
       (getter)PyBobMathLpInteriorPoint_mu,
       0,
-      s_mu_doc,
+      s_mu.doc(),
       0
     },
     {0}  /* Sentinel */
 };
 
-PyDoc_STRVAR(s_reset_str, "reset");
-PyDoc_STRVAR(s_reset_doc,
-"o.reset(M, N) -> None\n\
-\n\
-Resets the size of the problem (M and N correspond to the dimensions of the\n\
-A matrix");
+static auto s_reset = xbob::extension::FunctionDoc(
+    "reset",
+    "Resets the size of the problem (M and N correspond to the dimensions of the A matrix)"
+  )
+  .add_prototype("M, N")
+  .add_parameter("M", "int", "The new first dimension of the problem/A matrix")
+  .add_parameter("N", "int", "The new second dimension of the problem/A matrix")
+;
 
 static PyObject* PyBobMathLpInteriorPoint_reset
 (PyBobMathLpInteriorPointObject* self, PyObject *args, PyObject* kwds) {
@@ -262,12 +260,15 @@ static PyObject* PyBobMathLpInteriorPoint_reset
 
 }
 
-PyDoc_STRVAR(s_solve_str, "solve");
-PyDoc_STRVAR(s_solve_doc,
-"o.solve(A, b, c, x0, [lambda, mu]) -> x\n\
-\n\
-Solves an LP problem\n\
-");
+static auto s_solve = xbob::extension::FunctionDoc(
+    "solve",
+    "Solves an LP problem"
+  )
+  .add_prototype("A, b, c, x0, lambda, mu", "x")
+  .add_parameter("lambda", "?, optional", ".. todo:: Document parameter labmda")
+  .add_parameter("mu", "?, optional", ".. todo:: Document parameter mu")
+;
+
 
 static PyObject* PyBobMathLpInteriorPoint_solve
 (PyBobMathLpInteriorPointObject* self, PyObject *args, PyObject* kwds) {
@@ -380,14 +381,13 @@ static PyObject* PyBobMathLpInteriorPoint_solve
 
 }
 
-PyDoc_STRVAR(s_is_feasible_str, "is_feasible");
-PyDoc_STRVAR(s_is_feasible_doc,
-"o.is_feasible(A, b, c, x, lambda, mu) -> bool\n\
-\n\
-Checks if a primal-dual point (x,lambda,mu) belongs to the set of\n\
-feasible point (i.e. fulfill the constraints)\n\
-\n\
-");
+static auto s_is_feasible = xbob::extension::FunctionDoc(
+    "is_feasible",
+    "Checks if a primal-dual point (x, lambda, mu) belongs to the set of feasible points (i.e. fulfills the constraints)."
+  )
+  .add_prototype("A, b, c, x, lambda, mu", "test")
+  .add_return("test", "bool", "``True`` if (x, labmda, mu) belongs to the set of feasible points, otherwise ``False``")
+;
 
 static PyObject* PyBobMathLpInteriorPoint_is_feasible
 (PyBobMathLpInteriorPointObject* self, PyObject *args, PyObject* kwds) {
@@ -476,14 +476,15 @@ static PyObject* PyBobMathLpInteriorPoint_is_feasible
 
 }
 
-PyDoc_STRVAR(s_is_in_v_str, "is_in_v");
-PyDoc_STRVAR(s_is_in_v_doc,
-"o.is_in_v(x, mu, theta) -> bool\n\
-\n\
-Checks if a primal-dual point (x,lambda,mu) belongs to the V2\n\
-neighborhood of the central path.\n\
-\n\
-");
+static auto s_is_in_v = xbob::extension::FunctionDoc(
+    "is_in_v",
+    "Checks if a primal-dual point (x, lambda, mu) belongs to the V2 neighborhood of the central path.",
+    ".. todo:: This documentation seems wrong since lambda is not in the list of parameters."
+  )
+  .add_prototype("x, mu, theta", "test")
+  .add_return("test", "bool", "``True`` if (x, labmda, mu) belongs to the V2 neighborhood of the central path, otherwise ``False``")
+;
+
 
 static PyObject* PyBobMathLpInteriorPoint_is_in_v
 (PyBobMathLpInteriorPointObject* self, PyObject *args, PyObject* kwds) {
@@ -539,14 +540,13 @@ static PyObject* PyBobMathLpInteriorPoint_is_in_v
 
 }
 
-PyDoc_STRVAR(s_is_in_v_s_str, "is_in_v_s");
-PyDoc_STRVAR(s_is_in_v_s_doc,
-"o.is_in_v_s(A, b, c, x, lambda, mu) -> bool\n\
-\n\
-Checks if a primal-dual point (x,lambda,mu) belongs to the V\n\
-neighborhood of the central path and the set of feasible points.\n\
-\n\
-");
+static auto s_is_in_v_s = xbob::extension::FunctionDoc(
+    "is_in_v_s",
+    "Checks if a primal-dual point (x,lambda,mu) belongs to the V neighborhood of the central path and the set of feasible points."
+  )
+  .add_prototype("A, b, c, x, lambda, mu", "test")
+  .add_return("test", "bool", "``True`` if (x, labmda, mu) belongs to the V neighborhood of the central path and the set of feasible points, otherwise ``False``")
+;
 
 static PyObject* PyBobMathLpInteriorPoint_is_in_v_s
 (PyBobMathLpInteriorPointObject* self, PyObject *args, PyObject* kwds) {
@@ -638,14 +638,12 @@ static PyObject* PyBobMathLpInteriorPoint_is_in_v_s
 
 }
 
-PyDoc_STRVAR(s_initialize_dual_lambda_mu_str, "initialize_dual_lambda_mu");
-PyDoc_STRVAR(s_initialize_dual_lambda_mu_doc,
-"o.initialize_dual_lambda_mu(A, c) -> None\n\
-\n\
-Initializes the dual variables ``lambda`` and ``mu`` by minimizing the\n\
-logarithmic barrier function.\n\
-\n\
-");
+static auto s_initialize_dual_lambda_mu = xbob::extension::FunctionDoc(
+    "initialize_dual_lambda_mu",
+    "Initializes the dual variables ``lambda`` and ``mu`` by minimizing the logarithmic barrier function."
+  )
+  .add_prototype("A, c")
+;
 
 static PyObject* PyBobMathLpInteriorPoint_initialize_dual_lambda_mu
 (PyBobMathLpInteriorPointObject* self, PyObject *args, PyObject* kwds) {
@@ -697,40 +695,40 @@ static PyObject* PyBobMathLpInteriorPoint_initialize_dual_lambda_mu
 
 static PyMethodDef PyBobMathLpInteriorPoint_methods[] = {
     {
-      s_reset_str,
+      s_reset.name(),
       (PyCFunction)PyBobMathLpInteriorPoint_reset,
       METH_VARARGS|METH_KEYWORDS,
-      s_reset_doc
+      s_reset.doc()
     },
     {
-      s_solve_str,
+      s_solve.name(),
       (PyCFunction)PyBobMathLpInteriorPoint_solve,
       METH_VARARGS|METH_KEYWORDS,
-      s_solve_doc
+      s_solve.doc()
     },
     {
-      s_is_feasible_str,
+      s_is_feasible.name(),
       (PyCFunction)PyBobMathLpInteriorPoint_is_feasible,
       METH_VARARGS|METH_KEYWORDS,
-      s_is_feasible_doc
+      s_is_feasible.doc()
     },
     {
-      s_is_in_v_str,
+      s_is_in_v.name(),
       (PyCFunction)PyBobMathLpInteriorPoint_is_in_v,
       METH_VARARGS|METH_KEYWORDS,
-      s_is_in_v_doc
+      s_is_in_v.doc()
     },
     {
-      s_is_in_v_s_str,
+      s_is_in_v_s.name(),
       (PyCFunction)PyBobMathLpInteriorPoint_is_in_v_s,
       METH_VARARGS|METH_KEYWORDS,
-      s_is_in_v_s_doc
+      s_is_in_v_s.doc()
     },
     {
-      s_initialize_dual_lambda_mu_str,
+      s_initialize_dual_lambda_mu.name(),
       (PyCFunction)PyBobMathLpInteriorPoint_initialize_dual_lambda_mu,
       METH_VARARGS|METH_KEYWORDS,
-      s_initialize_dual_lambda_mu_doc
+      s_initialize_dual_lambda_mu.doc()
     },
     {0}  /* Sentinel */
 };
@@ -767,7 +765,7 @@ static PyObject* PyBobMathLpInteriorPoint_RichCompare (PyBobMathLpInteriorPointO
 
 PyTypeObject PyBobMathLpInteriorPoint_Type = {
     PyVarObject_HEAD_INIT(0, 0)
-    s_lpinteriorpoint_str,                             /*tp_name*/
+    s_lpinteriorpoint.name(),                          /*tp_name*/
     sizeof(PyBobMathLpInteriorPointObject),            /*tp_basicsize*/
     0,                                                 /*tp_itemsize*/
     0,                                                 /*tp_dealloc*/
@@ -786,7 +784,7 @@ PyTypeObject PyBobMathLpInteriorPoint_Type = {
     0,                                                 /*tp_setattro*/
     0,                                                 /*tp_as_buffer*/
     Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE,          /*tp_flags*/
-    s_lpinteriorpoint_doc,                             /* tp_doc */
+    s_lpinteriorpoint.doc(),                           /* tp_doc */
     0,		                                             /* tp_traverse */
     0,		                                             /* tp_clear */
     (richcmpfunc)PyBobMathLpInteriorPoint_RichCompare, /* tp_richcompare */
@@ -810,37 +808,29 @@ PyTypeObject PyBobMathLpInteriorPoint_Type = {
  * Implementation of LPInteriorPointShortstep class *
  ****************************************************/
 
-PyDoc_STRVAR(s_lpinteriorpointshortstep_str, XBOB_EXT_MODULE_PREFIX ".LPInteriorPointShortstep");
-
-PyDoc_STRVAR(s_lpinteriorpointshortstep_doc,
-"LPInteriorPointShortstep(M, N, theta, epsilon) -> new LPInteriorPointShortstep\n\
-LPInteriorPointShortstep(solver) -> new LPInteriorPointShortstep\n\
-\n\
-A Linear Program solver based on a short step interior point method.\n\
-\n\
-See :py:class:`LPInteriorPoint` for more details on the base class.\n\
-\n\
-Objects of this class can be initialized in two different ways: a\n\
-detailed constructor with the parameters described below or a copy\n\
-constructor, that deep-copies the input object and creates a new\n\
-object (**not** a new reference to the same object).\n\
-\n\
-Constructor parameters:\n\
-\n\
-M\n\
-  (int) first dimension of the A matrix\n\
-\n\
-N\n\
-  (int) second dimension of the A matrix\n\
-\n\
-theta\n\
-  (float) theta The value defining the size of the V2 neighborhood\n\
-\n\
-epsilon\n\
-  (float) The precision to determine whether an equality constraint\n\
-  is fulfilled or not.\n\
-\n\
-");
+const auto s_lpinteriorpointshortstep = xbob::extension::ClassDoc(
+    XBOB_EXT_MODULE_PREFIX ".LPInteriorPointShortstep",
+    "A Linear Program solver based on a short step interior point method.\n"
+    "See :py:class:`LPInteriorPoint` for more details on the base class."
+  )
+  .add_constructor(xbob::extension::FunctionDoc(
+      "LPInteriorPointShortstep",
+      "Objects of this class can be initialized in two different ways: "
+      "a detailed constructor with the parameters described below or "
+      "a copy constructor that deep-copies the input object and creates a new object (**not** a new reference to the same object)."
+    )
+    .add_prototype("M, N, theta, epsilon", "")
+    .add_prototype("solver", "")
+    .add_parameter("M", "int", "first dimension of the A matrix")
+    .add_parameter("N", "int", "second dimension of the A matrix")
+    .add_parameter("theta", "float", "The value defining the size of the V2 neighborhood")
+    .add_parameter("epsilon", "float", "The precision to determine whether an equality constraint is fulfilled or not.")
+    .add_parameter("solver", "LPInteriorPointShortstep", "The solver to make a deep copy of")
+  )
+  .highlight(s_solve)
+  .highlight(s_mu)
+  .highlight(s_lambda)
+;
 
 /* Type definition for PyBobMathLpInteriorPointObject */
 typedef struct {
@@ -866,7 +856,7 @@ static int PyBobMathLpInteriorPointShortstep_init1(PyBobMathLpInteriorPointShort
   if (!PyArg_ParseTupleAndKeywords(args, kwds, "O", kwlist, &solver)) return -1;
 
   if (!PyBobMathLpInteriorPointShortstep_Check(solver)) {
-    PyErr_Format(PyExc_TypeError, "copy-constructor for %s requires an object of the same type, not %s", s_lpinteriorpointshortstep_str, solver->ob_type->tp_name);
+    PyErr_Format(PyExc_TypeError, "copy-constructor for %s requires an object of the same type, not %s", s_lpinteriorpointshortstep.name(), solver->ob_type->tp_name);
     return -1;
   }
 
@@ -879,7 +869,7 @@ static int PyBobMathLpInteriorPointShortstep_init1(PyBobMathLpInteriorPointShort
     PyErr_SetString(PyExc_RuntimeError, ex.what());
   }
   catch (...) {
-    PyErr_Format(PyExc_RuntimeError, "cannot deep-copy object of type %s: unknown exception caught", s_lpinteriorpointshortstep_str);
+    PyErr_Format(PyExc_RuntimeError, "cannot deep-copy object of type %s: unknown exception caught", s_lpinteriorpointshortstep.name());
   }
 
   self->parent.base = self->base;
@@ -912,7 +902,7 @@ static int PyBobMathLpInteriorPointShortstep_init4(PyBobMathLpInteriorPointShort
     return -1;
   }
   catch (...) {
-    PyErr_Format(PyExc_RuntimeError, "cannot construct object of type %s: unknown exception caught", s_lpinteriorpointshortstep_str);
+    PyErr_Format(PyExc_RuntimeError, "cannot construct object of type %s: unknown exception caught", s_lpinteriorpointshortstep.name());
     return -1;
   }
 
@@ -938,7 +928,7 @@ static int PyBobMathLpInteriorPointShortstep_init(PyBobMathLpInteriorPointShorts
 
     default:
 
-      PyErr_Format(PyExc_RuntimeError, "number of arguments mismatch - %s requires 1 or 4 arguments, but you provided %" PY_FORMAT_SIZE_T "d (see help)", s_lpinteriorpointshortstep_str, nargs);
+      PyErr_Format(PyExc_RuntimeError, "number of arguments mismatch - %s requires 1 or 4 arguments, but you provided %" PY_FORMAT_SIZE_T "d (see help)", s_lpinteriorpointshortstep.name(), nargs);
 
   }
 
@@ -955,9 +945,10 @@ static void PyBobMathLpInteriorPointShortstep_delete (PyBobMathLpInteriorPointSh
 
 }
 
-PyDoc_STRVAR(s_theta_str, "theta");
-PyDoc_STRVAR(s_theta_doc,
-"The value theta used to define a V2 neighborhood"
+static auto s_theta = xbob::extension::VariableDoc(
+  "theta",
+  "float",
+  "The value theta used to define a V2 neighborhood"
 );
 
 static PyObject* PyBobMathLpInteriorPointShortstep_getTheta (PyBobMathLpInteriorPointShortstepObject* self, void* /*closure*/) {
@@ -978,7 +969,7 @@ static int PyBobMathLpInteriorPointShortstep_setTheta (PyBobMathLpInteriorPointS
     return -1;
   }
   catch (...) {
-    PyErr_Format(PyExc_RuntimeError, "cannot reset `theta' of %s: unknown exception caught", s_lpinteriorpointshortstep_str);
+    PyErr_Format(PyExc_RuntimeError, "cannot reset `theta' of %s: unknown exception caught", s_lpinteriorpointshortstep.name());
     return -1;
   }
 
@@ -988,10 +979,10 @@ static int PyBobMathLpInteriorPointShortstep_setTheta (PyBobMathLpInteriorPointS
 
 static PyGetSetDef PyBobMathLpInteriorPointShortstep_getseters[] = {
     {
-      s_theta_str,
+      s_theta.name(),
       (getter)PyBobMathLpInteriorPointShortstep_getTheta,
       (setter)PyBobMathLpInteriorPointShortstep_setTheta,
-      s_theta_doc,
+      s_theta.doc(),
       0
     },
     {0}  /* Sentinel */
@@ -1002,7 +993,7 @@ static PyObject* PyBobMathLpInteriorPointShortstep_RichCompare
 
   if (!PyBobMathLpInteriorPointShortstep_Check(other)) {
     PyErr_Format(PyExc_TypeError, "cannot compare `%s' with `%s'",
-        s_lpinteriorpointshortstep_str, other->ob_type->tp_name);
+        s_lpinteriorpointshortstep.name(), other->ob_type->tp_name);
     return 0;
   }
 
@@ -1024,9 +1015,10 @@ static PyObject* PyBobMathLpInteriorPointShortstep_RichCompare
 
 }
 
+
 PyTypeObject PyBobMathLpInteriorPointShortstep_Type = {
     PyVarObject_HEAD_INIT(0, 0)
-    s_lpinteriorpointshortstep_str,                             /*tp_name*/
+    s_lpinteriorpointshortstep.name(),                          /*tp_name*/
     sizeof(PyBobMathLpInteriorPointShortstepObject),            /*tp_basicsize*/
     0,                                                          /*tp_itemsize*/
     (destructor)PyBobMathLpInteriorPointShortstep_delete,       /*tp_dealloc*/
@@ -1045,7 +1037,7 @@ PyTypeObject PyBobMathLpInteriorPointShortstep_Type = {
     0,                                                          /*tp_setattro*/
     0,                                                          /*tp_as_buffer*/
     Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE,                   /*tp_flags*/
-    s_lpinteriorpointshortstep_doc,                             /* tp_doc */
+    s_lpinteriorpointshortstep.doc(),                           /* tp_doc */
     0,		                                                      /* tp_traverse */
     0,		                                                      /* tp_clear */
     (richcmpfunc)PyBobMathLpInteriorPointShortstep_RichCompare, /* tp_richcompare */
@@ -1069,41 +1061,27 @@ PyTypeObject PyBobMathLpInteriorPointShortstep_Type = {
  * Implementation of LPInteriorPointPredictorCorrector class *
  *************************************************************/
 
-PyDoc_STRVAR(s_lpinteriorpointpredictorcorrector_str, XBOB_EXT_MODULE_PREFIX ".LPInteriorPointPredictorCorrector");
-
-PyDoc_STRVAR(s_lpinteriorpointpredictorcorrector_doc,
-"LPInteriorPointPredictorCorrector(M, N, theta_pred, theta_corr, epsilon) -> new LPInteriorPointPredictorCorrector\n\
-LPInteriorPointPredictorCorrector(solver) -> new LPInteriorPointPredictorCorrector\n\
-\n\
-A Linear Program solver based on a predictor-corrector interior point\n\
-method.\n\
-\n\
-See :py:class:`LPInteriorPoint` for more details on the base class.\n\
-\n\
-Objects of this class can be initialized in two different ways: a\n\
-detailed constructor with the parameters described below or a copy\n\
-constructor, that deep-copies the input object and creates a new\n\
-object (**not** a new reference to the same object).\n\
-\n\
-Constructor parameters:\n\
-\n\
-M\n\
-  (int) first dimension of the A matrix\n\
-\n\
-N\n\
-  (int) second dimension of the A matrix\n\
-\n\
-theta_pred\n\
-  (float) the value theta_pred used to define a V2 neighborhood\n\
-\n\
-theta_corr\n\
-  (float) the value theta_corr used to define a V2 neighborhood\n\
-\n\
-epsilon\n\
-  (float) The precision to determine whether an equality constraint\n\
-  is fulfilled or not.\n\
-\n\
-");
+static auto s_lpinteriorpointpredictorcorrector = xbob::extension::ClassDoc(
+    XBOB_EXT_MODULE_PREFIX ".LPInteriorPointPredictorCorrector",
+    "A Linear Program solver based on a predictor-corrector interior point method.",
+    "See :py:class:`LPInteriorPoint` for more details on the base class."
+  )
+  .add_constructor(xbob::extension::FunctionDoc(
+      "LPInteriorPointPredictorCorrector",
+      "Objects of this class can be initialized in two different ways: "
+      "a detailed constructor with the parameters described below or "
+      "a copy constructor, that deep-copies the input object and creates a new object (**not** a new reference to the same object)."
+    )
+    .add_prototype("M, N, theta_pred, theta_corr, epsilon", "")
+    .add_prototype("solver", "")
+    .add_parameter("M", "int", "first dimension of the A matrix")
+    .add_parameter("N", "int", "second dimension of the A matrix")
+    .add_parameter("theta_pred", "float", "the value theta_pred used to define a V2 neighborhood")
+    .add_parameter("theta_corr", "float", "the value theta_corr used to define a V2 neighborhood")
+    .add_parameter("epsilon", "float", "the precision to determine whether an equality constraint is fulfilled or not")
+    .add_parameter("solver", "LPInteriorPointPredictorCorrector", "the solver to make a deep copy of")
+  )
+;
 
 /* Type definition for PyBobMathLpInteriorPointObject */
 typedef struct {
@@ -1129,7 +1107,7 @@ static int PyBobMathLpInteriorPointPredictorCorrector_init1(PyBobMathLpInteriorP
   if (!PyArg_ParseTupleAndKeywords(args, kwds, "O", kwlist, &solver)) return -1;
 
   if (!PyBobMathLpInteriorPointPredictorCorrector_Check(solver)) {
-    PyErr_Format(PyExc_TypeError, "copy-constructor for %s requires an object of the same type, not %s", s_lpinteriorpointpredictorcorrector_str, solver->ob_type->tp_name);
+    PyErr_Format(PyExc_TypeError, "copy-constructor for %s requires an object of the same type, not %s", s_lpinteriorpointpredictorcorrector.name(), solver->ob_type->tp_name);
     return -1;
   }
 
@@ -1142,7 +1120,7 @@ static int PyBobMathLpInteriorPointPredictorCorrector_init1(PyBobMathLpInteriorP
     PyErr_SetString(PyExc_RuntimeError, ex.what());
   }
   catch (...) {
-    PyErr_Format(PyExc_RuntimeError, "cannot deep-copy object of type %s: unknown exception caught", s_lpinteriorpointpredictorcorrector_str);
+    PyErr_Format(PyExc_RuntimeError, "cannot deep-copy object of type %s: unknown exception caught", s_lpinteriorpointpredictorcorrector.name());
   }
 
   self->parent.base = self->base;
@@ -1176,7 +1154,7 @@ static int PyBobMathLpInteriorPointPredictorCorrector_init5(PyBobMathLpInteriorP
     return -1;
   }
   catch (...) {
-    PyErr_Format(PyExc_RuntimeError, "cannot construct object of type %s: unknown exception caught", s_lpinteriorpointpredictorcorrector_str);
+    PyErr_Format(PyExc_RuntimeError, "cannot construct object of type %s: unknown exception caught", s_lpinteriorpointpredictorcorrector.name());
     return -1;
   }
 
@@ -1202,7 +1180,7 @@ static int PyBobMathLpInteriorPointPredictorCorrector_init(PyBobMathLpInteriorPo
 
     default:
 
-      PyErr_Format(PyExc_RuntimeError, "number of arguments mismatch - %s requires 1 or 5 arguments, but you provided %" PY_FORMAT_SIZE_T "d (see help)", s_lpinteriorpointpredictorcorrector_str, nargs);
+      PyErr_Format(PyExc_RuntimeError, "number of arguments mismatch - %s requires 1 or 5 arguments, but you provided %" PY_FORMAT_SIZE_T "d (see help)", s_lpinteriorpointpredictorcorrector.name(), nargs);
 
   }
 
@@ -1219,9 +1197,10 @@ static void PyBobMathLpInteriorPointPredictorCorrector_delete (PyBobMathLpInteri
 
 }
 
-PyDoc_STRVAR(s_theta_pred_str, "theta_pred");
-PyDoc_STRVAR(s_theta_pred_doc,
-"The value theta_pred used to define a V2 neighborhood"
+static auto s_theta_pred = xbob::extension::VariableDoc(
+  "theta_pred",
+  "float",
+  "The value theta_pred used to define a V2 neighborhood"
 );
 
 static PyObject* PyBobMathLpInteriorPointPredictorCorrector_getThetaPred (PyBobMathLpInteriorPointPredictorCorrectorObject* self, void* /*closure*/) {
@@ -1241,7 +1220,7 @@ static int PyBobMathLpInteriorPointPredictorCorrector_setThetaPred (PyBobMathLpI
     return -1;
   }
   catch (...) {
-    PyErr_Format(PyExc_RuntimeError, "cannot reset `theta_pred' of %s: unknown exception caught", s_lpinteriorpointpredictorcorrector_str);
+    PyErr_Format(PyExc_RuntimeError, "cannot reset `theta_pred' of %s: unknown exception caught", s_lpinteriorpointpredictorcorrector.name());
     return -1;
   }
 
@@ -1249,9 +1228,10 @@ static int PyBobMathLpInteriorPointPredictorCorrector_setThetaPred (PyBobMathLpI
 
 }
 
-PyDoc_STRVAR(s_theta_corr_str, "theta_corr");
-PyDoc_STRVAR(s_theta_corr_doc,
-"The value theta_corr used to define a V2 neighborhood"
+static auto s_theta_corr = xbob::extension::VariableDoc(
+  "theta_corr",
+  "float",
+  "The value theta_corr used to define a V2 neighborhood"
 );
 
 static PyObject* PyBobMathLpInteriorPointPredictorCorrector_getThetaCorr (PyBobMathLpInteriorPointPredictorCorrectorObject* self, void* /*closure*/) {
@@ -1271,7 +1251,7 @@ static int PyBobMathLpInteriorPointPredictorCorrector_setThetaCorr (PyBobMathLpI
     return -1;
   }
   catch (...) {
-    PyErr_Format(PyExc_RuntimeError, "cannot reset `theta_corr' of %s: unknown exception caught", s_lpinteriorpointpredictorcorrector_str);
+    PyErr_Format(PyExc_RuntimeError, "cannot reset `theta_corr' of %s: unknown exception caught", s_lpinteriorpointpredictorcorrector.name());
     return -1;
   }
 
@@ -1281,17 +1261,17 @@ static int PyBobMathLpInteriorPointPredictorCorrector_setThetaCorr (PyBobMathLpI
 
 static PyGetSetDef PyBobMathLpInteriorPointPredictorCorrector_getseters[] = {
     {
-      s_theta_pred_str,
+      s_theta_pred.name(),
       (getter)PyBobMathLpInteriorPointPredictorCorrector_getThetaPred,
       (setter)PyBobMathLpInteriorPointPredictorCorrector_setThetaPred,
-      s_theta_pred_doc,
+      s_theta_pred.doc(),
       0
     },
     {
-      s_theta_corr_str,
+      s_theta_corr.name(),
       (getter)PyBobMathLpInteriorPointPredictorCorrector_getThetaCorr,
       (setter)PyBobMathLpInteriorPointPredictorCorrector_setThetaCorr,
-      s_theta_corr_doc,
+      s_theta_corr.doc(),
       0
     },
     {0}  /* Sentinel */
@@ -1302,7 +1282,7 @@ static PyObject* PyBobMathLpInteriorPointPredictorCorrector_RichCompare
 
   if (!PyBobMathLpInteriorPointPredictorCorrector_Check(other)) {
     PyErr_Format(PyExc_TypeError, "cannot compare `%s' with `%s'",
-        s_lpinteriorpointpredictorcorrector_str, other->ob_type->tp_name);
+        s_lpinteriorpointpredictorcorrector.name(), other->ob_type->tp_name);
     return 0;
   }
 
@@ -1326,7 +1306,7 @@ static PyObject* PyBobMathLpInteriorPointPredictorCorrector_RichCompare
 
 PyTypeObject PyBobMathLpInteriorPointPredictorCorrector_Type = {
     PyVarObject_HEAD_INIT(0, 0)
-    s_lpinteriorpointpredictorcorrector_str,                    /*tp_name*/
+    s_lpinteriorpointpredictorcorrector.name(),                 /*tp_name*/
     sizeof(PyBobMathLpInteriorPointPredictorCorrectorObject),   /*tp_basicsize*/
     0,                                                          /*tp_itemsize*/
     (destructor)PyBobMathLpInteriorPointPredictorCorrector_delete, /*tp_dealloc*/
@@ -1345,7 +1325,7 @@ PyTypeObject PyBobMathLpInteriorPointPredictorCorrector_Type = {
     0,                                                          /*tp_setattro*/
     0,                                                          /*tp_as_buffer*/
     Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE,                   /*tp_flags*/
-    s_lpinteriorpointpredictorcorrector_doc,                    /* tp_doc */
+    s_lpinteriorpointpredictorcorrector.doc(),                  /* tp_doc */
     0,		                                                      /* tp_traverse */
     0,		                                                      /* tp_clear */
     (richcmpfunc)PyBobMathLpInteriorPointPredictorCorrector_RichCompare, /* tp_richcompare */
@@ -1369,40 +1349,27 @@ PyTypeObject PyBobMathLpInteriorPointPredictorCorrector_Type = {
  * Implementation of LPInteriorPointLongstep class *
  ****************************************************/
 
-PyDoc_STRVAR(s_lpinteriorpointlongstep_str, XBOB_EXT_MODULE_PREFIX ".LPInteriorPointLongstep");
-
-PyDoc_STRVAR(s_lpinteriorpointlongstep_doc,
-"LPInteriorPointLongstep(M, N, gamma, sigma, epsilon) -> new LPInteriorPointLongstep\n\
-LPInteriorPointLongstep(solver) -> new LPInteriorPointLongstep\n\
-\n\
-A Linear Program solver based on a long step interior point method.\n\
-\n\
-See :py:class:`LPInteriorPoint` for more details on the base class.\n\
-\n\
-Objects of this class can be initialized in two different ways: a\n\
-detailed constructor with the parameters described below or a copy\n\
-constructor, that deep-copies the input object and creates a new\n\
-object (**not** a new reference to the same object).\n\
-\n\
-Constructor parameters:\n\
-\n\
-M\n\
-  (int) first dimension of the A matrix\n\
-\n\
-N\n\
-  (int) second dimension of the A matrix\n\
-\n\
-gamma\n\
-  (float) The value gamma used to define a V-inf neighborhood\n\
-\n\
-sigma\n\
-  (float) The value sigma used to define a V-inf neighborhood\n\
-\n\
-epsilon\n\
-  (float) The precision to determine whether an equality constraint\n\
-  is fulfilled or not.\n\
-\n\
-");
+static auto s_lpinteriorpointlongstep = xbob::extension::ClassDoc(
+    XBOB_EXT_MODULE_PREFIX ".LPInteriorPointLongstep",
+    "A Linear Program solver based on a long step interior point method.",
+    "See :py:class:`LPInteriorPoint` for more details on the base class."
+  )
+  .add_constructor(xbob::extension::FunctionDoc(
+      "LPInteriorPointLongstep",
+      "Objects of this class can be initialized in two different ways: "
+      "a detailed constructor with the parameters described below or "
+      "a copy constructor, that deep-copies the input object and creates a new object (**not** a new reference to the same object)"
+    )
+    .add_prototype("M, N, gamma, sigma, epsilon", "")
+    .add_prototype("solver", "")
+    .add_parameter("M", "int", "first dimension of the A matrix")
+    .add_parameter("N", "int", "second dimension of the A matrix")
+    .add_parameter("gamma", "float", "the value gamma used to define a V-inf neighborhood")
+    .add_parameter("sigma", "float", "the value sigma used to define a V-inf neighborhood")
+    .add_parameter("epsilon", "float", "the precision to determine whether an equality constraint is fulfilled or not")
+    .add_parameter("solver", "LPInteriorPointLongstep", "the solver to make a deep copy of")
+  )
+;
 
 typedef struct {
   PyBobMathLpInteriorPointObject parent;
@@ -1427,7 +1394,7 @@ static int PyBobMathLpInteriorPointLongstep_init1(PyBobMathLpInteriorPointLongst
   if (!PyArg_ParseTupleAndKeywords(args, kwds, "O", kwlist, &solver)) return -1;
 
   if (!PyBobMathLpInteriorPointLongstep_Check(solver)) {
-    PyErr_Format(PyExc_TypeError, "copy-constructor for %s requires an object of the same type, not %s", s_lpinteriorpointlongstep_str, solver->ob_type->tp_name);
+    PyErr_Format(PyExc_TypeError, "copy-constructor for %s requires an object of the same type, not %s", s_lpinteriorpointlongstep.name(), solver->ob_type->tp_name);
     return -1;
   }
 
@@ -1440,7 +1407,7 @@ static int PyBobMathLpInteriorPointLongstep_init1(PyBobMathLpInteriorPointLongst
     PyErr_SetString(PyExc_RuntimeError, ex.what());
   }
   catch (...) {
-    PyErr_Format(PyExc_RuntimeError, "cannot deep-copy object of type %s: unknown exception caught", s_lpinteriorpointlongstep_str);
+    PyErr_Format(PyExc_RuntimeError, "cannot deep-copy object of type %s: unknown exception caught", s_lpinteriorpointlongstep.name());
   }
 
   self->parent.base = self->base;
@@ -1474,7 +1441,7 @@ static int PyBobMathLpInteriorPointLongstep_init5(PyBobMathLpInteriorPointLongst
     return -1;
   }
   catch (...) {
-    PyErr_Format(PyExc_RuntimeError, "cannot construct object of type %s: unknown exception caught", s_lpinteriorpointlongstep_str);
+    PyErr_Format(PyExc_RuntimeError, "cannot construct object of type %s: unknown exception caught", s_lpinteriorpointlongstep.name());
     return -1;
   }
   self->parent.base = self->base;
@@ -1499,7 +1466,7 @@ static int PyBobMathLpInteriorPointLongstep_init(PyBobMathLpInteriorPointLongste
 
     default:
 
-      PyErr_Format(PyExc_RuntimeError, "number of arguments mismatch - %s requires 1 or 5 arguments, but you provided %" PY_FORMAT_SIZE_T "d (see help)", s_lpinteriorpointlongstep_str, nargs);
+      PyErr_Format(PyExc_RuntimeError, "number of arguments mismatch - %s requires 1 or 5 arguments, but you provided %" PY_FORMAT_SIZE_T "d (see help)", s_lpinteriorpointlongstep.name(), nargs);
 
   }
 
@@ -1516,9 +1483,10 @@ static void PyBobMathLpInteriorPointLongstep_delete (PyBobMathLpInteriorPointLon
 
 }
 
-PyDoc_STRVAR(s_gamma_str, "gamma");
-PyDoc_STRVAR(s_gamma_doc,
-"The value gamma used to define a V-Inf neighborhood"
+static auto s_gamma = xbob::extension::VariableDoc(
+  "gamma",
+  "float",
+  "The value gamma used to define a V-Inf neighborhood"
 );
 
 static PyObject* PyBobMathLpInteriorPointLongstep_getGamma (PyBobMathLpInteriorPointLongstepObject* self, void* /*closure*/) {
@@ -1539,7 +1507,7 @@ static int PyBobMathLpInteriorPointLongstep_setGamma (PyBobMathLpInteriorPointLo
     return -1;
   }
   catch (...) {
-    PyErr_Format(PyExc_RuntimeError, "cannot reset `gamma' of %s: unknown exception caught", s_lpinteriorpointlongstep_str);
+    PyErr_Format(PyExc_RuntimeError, "cannot reset `gamma' of %s: unknown exception caught", s_lpinteriorpointlongstep.name());
     return -1;
   }
 
@@ -1547,9 +1515,10 @@ static int PyBobMathLpInteriorPointLongstep_setGamma (PyBobMathLpInteriorPointLo
 
 }
 
-PyDoc_STRVAR(s_sigma_str, "sigma");
-PyDoc_STRVAR(s_sigma_doc,
-"The value sigma used to define a V-Inf neighborhood"
+static auto s_sigma = xbob::extension::VariableDoc(
+  "sigma",
+  "float",
+  "The value sigma used to define a V-Inf neighborhood"
 );
 
 static PyObject* PyBobMathLpInteriorPointLongstep_getSigma (PyBobMathLpInteriorPointLongstepObject* self, void* /*closure*/) {
@@ -1570,7 +1539,7 @@ static int PyBobMathLpInteriorPointLongstep_setSigma (PyBobMathLpInteriorPointLo
     return -1;
   }
   catch (...) {
-    PyErr_Format(PyExc_RuntimeError, "cannot reset `sigma' of %s: unknown exception caught", s_lpinteriorpointlongstep_str);
+    PyErr_Format(PyExc_RuntimeError, "cannot reset `sigma' of %s: unknown exception caught", s_lpinteriorpointlongstep.name());
     return -1;
   }
 
@@ -1580,17 +1549,17 @@ static int PyBobMathLpInteriorPointLongstep_setSigma (PyBobMathLpInteriorPointLo
 
 static PyGetSetDef PyBobMathLpInteriorPointLongstep_getseters[] = {
     {
-      s_gamma_str,
+      s_gamma.name(),
       (getter)PyBobMathLpInteriorPointLongstep_getGamma,
       (setter)PyBobMathLpInteriorPointLongstep_setGamma,
-      s_gamma_doc,
+      s_gamma.doc(),
       0
     },
     {
-      s_sigma_str,
+      s_sigma.name(),
       (getter)PyBobMathLpInteriorPointLongstep_getSigma,
       (setter)PyBobMathLpInteriorPointLongstep_setSigma,
-      s_sigma_doc,
+      s_sigma.doc(),
       0
     },
     {0}  /* Sentinel */
@@ -1601,7 +1570,7 @@ static PyObject* PyBobMathLpInteriorPointLongstep_RichCompare
 
   if (!PyBobMathLpInteriorPointLongstep_Check(other)) {
     PyErr_Format(PyExc_TypeError, "cannot compare `%s' with `%s'",
-        s_lpinteriorpointlongstep_str, other->ob_type->tp_name);
+        s_lpinteriorpointlongstep.name(), other->ob_type->tp_name);
     return 0;
   }
 
@@ -1623,14 +1592,14 @@ static PyObject* PyBobMathLpInteriorPointLongstep_RichCompare
 
 }
 
-PyDoc_STRVAR(s_is_in_vinf_str, "is_in_v");
-PyDoc_STRVAR(s_is_in_vinf_doc,
-"o.is_in_v(x, mu, gamma) -> bool\n\
-\n\
-Checks if a primal-dual point (x,lambda,mu) belongs to the V-Inf\n\
-neighborhood of the central path.\n\
-\n\
-");
+static auto s_is_in_vinf = xbob::extension::FunctionDoc(
+    "is_in_v",
+    "Checks if a primal-dual point (x, lambda, mu) belongs to the V-Inf neighborhood of the central path.",
+    ".. todo:: This documenation looks wrong since lambda is not part of the parameters"
+  )
+  .add_prototype("x, mu, gamma", "test")
+  .add_return("test", "bool", "``True`` if (x, lambda, mu) belong to the  V-Inf neighborhood of the central path, otherwise ``False``")
+;
 
 static PyObject* PyBobMathLpInteriorPoint_is_in_vinf
 (PyBobMathLpInteriorPointObject* self, PyObject *args, PyObject* kwds) {
@@ -1688,17 +1657,17 @@ static PyObject* PyBobMathLpInteriorPoint_is_in_vinf
 
 static PyMethodDef PyBobMathLpInteriorPointLongstep_methods[] = {
     {
-      s_is_in_vinf_str,
+      s_is_in_vinf.name(),
       (PyCFunction)PyBobMathLpInteriorPoint_is_in_vinf,
       METH_VARARGS|METH_KEYWORDS,
-      s_is_in_vinf_doc
+      s_is_in_vinf.doc()
     },
     {0}  /* Sentinel */
 };
 
 PyTypeObject PyBobMathLpInteriorPointLongstep_Type = {
     PyVarObject_HEAD_INIT(0, 0)
-    s_lpinteriorpointlongstep_str,                              /*tp_name*/
+    s_lpinteriorpointlongstep.name(),                           /*tp_name*/
     sizeof(PyBobMathLpInteriorPointLongstepObject),             /*tp_basicsize*/
     0,                                                          /*tp_itemsize*/
     (destructor)PyBobMathLpInteriorPointLongstep_delete,        /*tp_dealloc*/
@@ -1717,7 +1686,7 @@ PyTypeObject PyBobMathLpInteriorPointLongstep_Type = {
     0,                                                          /*tp_setattro*/
     0,                                                          /*tp_as_buffer*/
     Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE,                   /*tp_flags*/
-    s_lpinteriorpointlongstep_doc,                              /* tp_doc */
+    s_lpinteriorpointlongstep.doc(),                            /* tp_doc */
     0,		                                                      /* tp_traverse */
     0,		                                                      /* tp_clear */
     (richcmpfunc)PyBobMathLpInteriorPointLongstep_RichCompare,  /* tp_richcompare */
diff --git a/xbob/math/main.cpp b/xbob/math/main.cpp
index 55731fe..b8432fd 100644
--- a/xbob/math/main.cpp
+++ b/xbob/math/main.cpp
@@ -11,6 +11,8 @@
 #include <xbob.blitz/capi.h>
 #include <xbob.blitz/cleanup.h>
 
+#include <xbob.extension/documentation.h>
+
 #include "histogram.h"
 #include "linsolve.h"
 #include "pavx.h"
@@ -18,484 +20,454 @@
 #include "scatter.h"
 #include "lp_interior_point.h"
 
-PyDoc_STRVAR(s_histogram_intersection_str, "histogram_intersection");
-PyDoc_STRVAR(s_histogram_intersection_doc,
-"histogram_intersection(h1, h2) -> scalar\n\
-histogram_intersection(index_1, value_1, index_2, value_2) -> scalar\n\
-\n\
-Computes the histogram intersection between the given histograms, which\n\
-might be of singular dimension only. The histogram intersection defines\n\
-a similarity measure, so higher values are better.\n\
-\n\
-You can use this method in two different formats. The first interface\n\
-accepts non-sparse histograms. The second interface accepts sparse\n\
-histograms represented by index and values.\n\
-"
-);
-
-PyDoc_STRVAR(s_chi_square_str, "chi_square");
-PyDoc_STRVAR(s_chi_square_doc,
-"chi_square(h1, h2) -> scalar\n\
-chi_square(index_1, value_1, index_2, value_2) -> scalar\n\
-\n\
-Computes the chi square distance between the given histograms, which\n\
-might be of singular dimension only. The chi square function is a \n\
-distance measure, so lower values are better.\n\
-\n\
-You can use this method in two different formats. The first interface\n\
-accepts non-sparse histograms. The second interface accepts sparse\n\
-histograms represented by index and values.\n\
-"
-);
-
-PyDoc_STRVAR(s_kullback_leibler_str, "kullback_leibler");
-PyDoc_STRVAR(s_kullback_leibler_doc,
-"kullback_leibler(h1, h2) -> scalar\n\
-kullback_leibler(index_1, value_1, index_2, value_2) -> scalar\n\
-\n\
-Computes the Kullback-Leibler histogram divergence between the given\n\
-histograms, which might be of singular dimension only. The\n\
-Kullback-Leibler divergence is a distance measure, so lower values\n\
-are better.\n\
-\n\
-You can use this method in two different formats. The first interface\n\
-accepts non-sparse histograms. The second interface accepts sparse\n\
-histograms represented by index and values.\n\
-"
-);
-
-PyDoc_STRVAR(s_linsolve_str, "linsolve");
-PyDoc_STRVAR(s_linsolve_doc,
-"linsolve(A, b) -> array\n\
-linsolve(A, x, b) -> None\n\
-\n\
-Solves the linear system :math:`Ax=b` and returns the result in ``x``.\n\
-This method uses LAPACK's ``dgesv`` generic solver.\n\
-\n\
-You can use this method in two different formats. The first interface\n\
-accepts the matrices ``A`` and ``b`` returning ``x``. The second one\n\
-accepts a pre-allocated ``x`` matrix and sets it with the linear system\n\
-solution.\n\
-"
-);
-
-PyDoc_STRVAR(s_linsolve_nocheck_str, "linsolve_");
-PyDoc_STRVAR(s_linsolve_nocheck_doc,
-"linsolve_(A, b) -> array\n\
-linsolve_(A, x, b) -> None\n\
-\n\
-Solves the linear system :math:`Ax=b` and returns the result in ``x``.\n\
-This method uses LAPACK's ``dgesv`` generic solver.\n\
-\n\
-.. warning::\n\
-\n\
-   THIS VARIANT DOES NOT PERFORM ANY CHECKS ON THE INPUT MATRICES AND IS,\n\
-   FASTER THEN THE VARIANT NOT ENDING IN ``_``. Use it when you are sure\n\
-   your input matrices sizes match.\n\
-\n\
-You can use this method in two different formats. The first interface\n\
-accepts the matrices ``A`` and ``b`` returning ``x``. The second one\n\
-accepts a pre-allocated ``x`` matrix and sets it with the linear system\n\
-solution.\n\
-"
-);
-
-PyDoc_STRVAR(s_linsolve_sympos_str, "linsolve_sympos");
-PyDoc_STRVAR(s_linsolve_sympos_doc,
-"linsolve_sympos(A, b) -> array\n\
-linsolve_sympos(A, x, b) -> None\n\
-\n\
-Solves the linear system :math:`Ax=b` and returns the result in ``x``.\n\
-This method uses LAPACK's ``dposv`` solver, assuming ``A`` is a symmetric.\n\
-positive definite matrix.\n\
-\n\
-You can use this method in two different formats. The first interface\n\
-accepts the matrices ``A`` and ``b`` returning ``x``. The second one\n\
-accepts a pre-allocated ``x`` matrix and sets it with the linear system\n\
-solution.\n\
-"
-);
-
-PyDoc_STRVAR(s_linsolve_sympos_nocheck_str, "linsolve_sympos_");
-PyDoc_STRVAR(s_linsolve_sympos_nocheck_doc,
-"linsolve_sympos_(A, b) -> array\n\
-linsolve_sympos_(A, x, b) -> None\n\
-\n\
-Solves the linear system :math:`Ax=b` and returns the result in ``x``.\n\
-This method uses LAPACK's ``dposv`` solver, assuming ``A`` is a symmetric.\n\
-positive definite matrix.\n\
-\n\
-.. warning::\n\
-\n\
-   THIS VARIANT DOES NOT PERFORM ANY CHECKS ON THE INPUT MATRICES AND IS,\n\
-   FASTER THEN THE VARIANT NOT ENDING IN ``_``. Use it when you are sure\n\
-   your input matrices sizes match.\n\
-\n\
-You can use this method in two different formats. The first interface\n\
-accepts the matrices ``A`` and ``b`` returning ``x``. The second one\n\
-accepts a pre-allocated ``x`` matrix and sets it with the linear system\n\
-solution.\n\
-"
-);
-
-PyDoc_STRVAR(s_linsolve_cg_sympos_str, "linsolve_cg_sympos");
-PyDoc_STRVAR(s_linsolve_cg_sympos_doc,
-"linsolve_cg_sympos(A, b) -> array\n\
-linsolve_cg_sympos(A, x, b) -> None\n\
-\n\
-Solves the linear system :math:`Ax=b` and returns the result in ``x``.\n\
-This method solves the linear system via conjugate gradients and assumes\n\
-``A`` is a symmetric positive definite matrix.\n\
-\n\
-You can use this method in two different formats. The first interface\n\
-accepts the matrices ``A`` and ``b`` returning ``x``. The second one\n\
-accepts a pre-allocated ``x`` matrix and sets it with the linear system\n\
-solution.\n\
-"
-);
-
-PyDoc_STRVAR(s_linsolve_cg_sympos_nocheck_str, "linsolve_cg_sympos_");
-PyDoc_STRVAR(s_linsolve_cg_sympos_nocheck_doc,
-"linsolve_cg_sympos_(A, b) -> array\n\
-linsolve_cg_sympos_(A, x, b) -> None\n\
-\n\
-Solves the linear system :math:`Ax=b` and returns the result in ``x``.\n\
-This method solves the linear system via conjugate gradients and assumes\n\
-``A`` is a symmetric positive definite matrix.\n\
-\n\
-.. warning::\n\
-\n\
-   THIS VARIANT DOES NOT PERFORM ANY CHECKS ON THE INPUT MATRICES AND IS,\n\
-   FASTER THEN THE VARIANT NOT ENDING IN ``_``. Use it when you are sure\n\
-   your input matrices sizes match.\n\
-\n\
-You can use this method in two different formats. The first interface\n\
-accepts the matrices ``A`` and ``b`` returning ``x``. The second one\n\
-accepts a pre-allocated ``x`` matrix and sets it with the linear system\n\
-solution.\n\
-"
-);
-
-PyDoc_STRVAR(s_pavx_str, "pavx");
-PyDoc_STRVAR(s_pavx_doc,
-"pavx(input, output) -> None\n\
-pavx(input) -> array\n\
-\n\
-Applies the Pool-Adjacent-Violators Algorithm to ``input``. The ``input``\n\
-and ``output`` arrays should have the same size. This is a simplified\n\
-C++ port of the isotonic regression code made available at the `University\n\
-of Bern website <http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html>`_.\n\
-\n\
-You can use this method in two different formats. The first interface\n\
-accepts the 1D float arrays ``input`` and ``output``. The second one\n\
-accepts the input array ``input`` and allocates a new ``output`` array\n\
-which is returned. In such a case, the ``output`` is a 1D float array\n\
-with the same length as ``input``.\n\
-");
-
-PyDoc_STRVAR(s_pavx_nocheck_str, "pavx_");
-PyDoc_STRVAR(s_pavx_nocheck_doc,
-"pavx(input, output) -> None\n\
-\n\
-Applies the Pool-Adjacent-Violators Algorithm to ``input`` and places the\n\
-result on ``output``. The ``input`` and ``output`` arrays should be 1D\n\
-float arrays with the same length.\n\
-\n\
-This is a simplified C++ port of the isotonic regression code\n\
-made available at the `University of Bern website <http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html>`_.\n\
-\n\
-.. warning::\n\
-\n\
-   THIS VARIANT DOES NOT PERFORM ANY CHECKS ON THE INPUT MATRICES AND IS,\n\
-   FASTER THEN THE VARIANT NOT ENDING IN ``_``. Use it when you are sure\n\
-   your input and output vector sizes match.\n\
-\n\
-");
-
-PyDoc_STRVAR(s_pavx_width_str, "pavxWidth");
-PyDoc_STRVAR(s_pavx_width_doc,
-"pavxWidth(input, output) -> array\n\
-\n\
-Applies the Pool-Adjacent-Violators Algorithm to ``input`` and places the\n\
-result on ``output``. The ``input`` and ``output`` arrays should be 1D\n\
-float arrays with the same length.\n\
-\n\
-The width array (64-bit unsigned integer 1D) is returned and has the\n\
-same size as ``input`` and ``output``.\n\
-");
-
-PyDoc_STRVAR(s_pavx_width_height_str, "pavxWidthHeight");
-PyDoc_STRVAR(s_pavx_width_height_doc,
-"pavxWidthHeight(input, output) -> (array, array)\n\
-\n\
-Applies the Pool-Adjacent-Violators Algorithm to ``input`` and sets the\n\
-result on ``output``. The ``input`` and ``output`` arrays should be 1D\n\
-float arrays of the same length.\n\
-\n\
-This is a simplified C++ port of the isotonic regression code\n\
-made available at the `University of Bern website <http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html>`_.\n\
-\n\
-The width and height arrays are returned. The width array is a 64-bit\n\
-**unsigned integer** 1D array, while the height array (second component\n\
-of the returned tuple) is a 64-bit **float** 1D array of the same size.\n\
-");
-
-PyDoc_STRVAR(s_norminv_str, "norminv");
-PyDoc_STRVAR(s_norminv_doc,
-"norminv(p, mu, sigma) -> scalar\n\
-\n\
-Computes the inverse normal cumulative distribution for a probability\n\
-``p``, given a distribution with mean ``mu`` and standard deviation\n\
-``sigma``. The value ``p`` must lie in the range [0,1].\n\
-\n\
-Reference: `<http://home.online.no/~pjacklam/notes/invnorm/>`_\n\
-");
-
-PyDoc_STRVAR(s_normsinv_str, "normsinv");
-PyDoc_STRVAR(s_normsinv_doc,
-"normsinv(p) -> scalar\n\
-\n\
-Computes the inverse normal cumulative distribution for a probability\n\
-``p``, given a distribution with mean 0.0 and standard deviation 1.0.\n\
-It is equivalent as calling :py:func:`norminv(p, 0, 1)`. The value\n\
-``p`` must lie in the range [0,1].\n\
-\n\
-Reference: `<http://home.online.no/~pjacklam/notes/invnorm/>`_\n\
-");
-
-PyDoc_STRVAR(s_scatter_str, "scatter");
-PyDoc_STRVAR(s_scatter_doc,
-"scatter(a, [s, [m]]) -> tuple\n\
-\n\
-Computes the scatter matrix of a 2D array *considering data is organized\n\
-row-wise* (each sample is a row, each feature is a column). The\n\
-resulting array ``s`` is squared with extents equal to the\n\
-number of columns in ``a``. The resulting array ``m`` is a 1D array\n\
-with the row means of ``a``. This method supports only 32 or 64-bit\n\
-float arrays as input.\n\
-\n\
-This function supports many calling modes, but you should provide, at\n\
-least, the input data matrix ``a``. All non-provided arguments will be\n\
-allocated internally and returned.\n\
-");
-
-PyDoc_STRVAR(s_scatter_nocheck_str, "scatter_");
-PyDoc_STRVAR(s_scatter_nocheck_doc,
-"scatter_(a, s, m) -> None\n\
-\n\
-Computes the scatter matrix of a 2D array *considering data is organized\n\
-row-wise* (each sample is a row, each feature is a column). The\n\
-resulting array ``s`` is squared with extents equal to the\n\
-number of columns in ``a``. The resulting array ``m`` is a 1D array\n\
-with the row means of ``a``. This method supports only 32 or 64-bit\n\
-float arrays as input.\n\
-\n\
-.. warning::\n\
-\n\
-   THIS VARIANT DOES NOT PERFORM ANY CHECKS ON THE INPUT MATRICES AND IS,\n\
-   FASTER THEN THE VARIANT NOT ENDING IN ``_``. Use it when you are sure\n\
-   your input matrices sizes match.\n\
-");
-
-PyDoc_STRVAR(s_scatters_str, "scatters");
-PyDoc_STRVAR(s_scatters_doc,
-"scatters(data, [sw, [sb, [m]]]) -> tuple\n\
-\n\
-Computes the within-class (``sw``) and between-class (``sb``) scatter\n\
-matrices of a set of 2D arrays considering data is organized row-wise\n\
-(each sample is a row, each feature is a column).\n\
-\n\
-This function supports many calling modes, but you should provide, at\n\
-least, the input data matrices ``data``, which should be sequence of\n\
-2D 32-bit or 64-bit float values in which every row of each matrix\n\
-represents one observation for a given class. **Every matrix in\n\
-``data`` should have exactly the same number of columns.**\n\
-\n\
-The returned values ``sw`` and ``sb`` are square matrices with the\n\
-same number of rows and columns as the number of columns in ``data``.\n\
-The returned value ``m`` (last call variant) is a 1D array with the\n\
-same length as number of columns in each ``data`` matrix and represents\n\
-the ensemble mean with no prior (i.e., biased towards classes with more\n\
-samples.\n\
-\n\
-Strategy implemented:\n\
-\n\
-1. Evaluate the overall mean (``m``), class means (:math:`m_k`) and the\n\
-   total class counts (:math:`N`).\n\
-2. Evaluate ``sw`` and ``sb`` using normal loops.\n\
-\n\
-Note that ``sw`` and ``sb``, in this implementation, will be normalized\n\
-by N-1 (number of samples) and K (number of classes). This procedure\n\
-makes the eigen values scaled by (N-1)/K, effectively increasing their\n\
-values. The main motivation for this normalization are numerical\n\
-precision concerns with the increasing number of samples causing a\n\
-rather large Sw matrix. A normalization strategy mitigates this\n\
-problem. The eigen vectors will see no effect on this normalization as\n\
-they are normalized in the euclidean sense (:math:`||a|| = 1`) so that\n\
-does not change those.\n\
-");
-
-PyDoc_STRVAR(s_scatters_nocheck_str, "scatters_");
-PyDoc_STRVAR(s_scatters_nocheck_doc,
-"scatters_(data, sw, sb) -> None\n\
-scatters_(data, sw, sb, m) -> None\n\
-\n\
-Computes the within-class (``sw``) and between-class (``sb``) scatter\n\
-matrices of a set of 2D arrays considering data is organized row-wise\n\
-(each sample is a row, each feature is a column).\n\
-\n\
-.. warning::\n\
-\n\
-   THIS VARIANT DOES NOT PERFORM ANY CHECKS ON THE INPUT MATRICES AND IS,\n\
-   FASTER THEN THE VARIANT NOT ENDING IN ``_``. Use it when you are sure\n\
-   your input matrices sizes match.\n\
-\n\
-This function supports many calling modes, but you should provide, at\n\
-least, the input data matrices ``data``, which should be sequence of\n\
-2D 32-bit or 64-bit float values in which every row of each matrix\n\
-represents one observation for a given class. **Every matrix in\n\
-``data`` should have exactly the same number of columns.**\n\
-\n\
-The returned values ``sw`` and ``sb`` are square matrices with the\n\
-same number of rows and columns as the number of columns in ``data``.\n\
-The returned value ``m`` (last call variant) is a 1D array with the\n\
-same length as number of columns in each ``data`` matrix and represents\n\
-the ensemble mean with no prior (i.e., biased towards classes with more\n\
-samples. **In this variant, you should pre-allocate all output matrices\n\
-so that scatters (and the overall mean) are stored on your provided\n\
-arrays**.\n\
-\n\
-Strategy implemented:\n\
-\n\
-1. Evaluate the overall mean (``m``), class means (:math:`m_k`) and the\n\
-   total class counts (:math:`N`).\n\
-2. Evaluate ``sw`` and ``sb`` using normal loops.\n\
-\n\
-Note that ``sw`` and ``sb``, in this implementation, will be normalized\n\
-by N-1 (number of samples) and K (number of classes). This procedure\n\
-makes the eigen values scaled by (N-1)/K, effectively increasing their\n\
-values. The main motivation for this normalization are numerical\n\
-precision concerns with the increasing number of samples causing a\n\
-rather large Sw matrix. A normalization strategy mitigates this\n\
-problem. The eigen vectors will see no effect on this normalization as\n\
-they are normalized in the euclidean sense (:math:`||a|| = 1`) so that\n\
-does not change those.\n\
-");
+static xbob::extension::FunctionDoc s_histogram_intersection = xbob::extension::FunctionDoc(
+    "histogram_intersection",
+    "Computes the histogram intersection between the given histograms, which might be of singular dimension only.",
+    "The histogram intersection is computed as follows:\n"
+    ".. math:: sim(h_1,h_2) = \\sum_i \\min \\{h_{1i}, h_{2i}\\}\n"
+    "The histogram intersection defines a similarity measure, so higher values are better. "
+    "You can use this method in two different formats. "
+    "The first interface accepts non-sparse histograms. "
+    "The second interface accepts sparse histograms represented by indexes and values.\n"
+    ".. note:: Histograms are given as two matrices, one with the indexes and one with the data. All data points that for which no index exists are considered to be zero.\n"
+    ".. note:: In general, histogram intersection with sparse histograms needs more time to be computed."
+  )
+  .add_prototype("h1, h2", "sim")
+  .add_prototype("index_1, value_1, index_2, value_2", "sim")
+  .add_parameter("h1, h2", "array_like (1D)", "Histograms to compute the histogram intersection for")
+  .add_parameter("index_1, index_2", "array_like (int, 1D)", "Indices of the sparse histograms value_1 and value_2")
+  .add_parameter("value_1, value_2", "array_like (1D)", "Sparse histograms to compute the histogram intersection for")
+  .add_return("sim", "float", "The histogram intersection value for the given histograms.")
+;
+
+static xbob::extension::FunctionDoc s_chi_square = xbob::extension::FunctionDoc(
+    "chi_square",
+    "Computes the chi square distance between the given histograms, which might be of singular dimension only.",
+    "The chi square distance is computed as follows:\n"
+    ".. math:: dist(h_1,h_2) = \\sum_i \\frac{(h_{1i} - h_{2i})^2}{h_{1i} + h_{2i}}\n"
+    "Chi square defines a distance metric, so lower values are better. "
+    "You can use this method in two different formats. "
+    "The first interface accepts non-sparse histograms. "
+    "The second interface accepts sparse histograms represented by indexes and values.\n"
+    ".. note:: Histograms are given as two matrices, one with the indexes and one with the data. All data points that for which no index exists are considered to be zero.\n"
+    ".. note:: In general, histogram intersection with sparse histograms needs more time to be computed."
+  )
+  .add_prototype("h1, h2", "dist")
+  .add_prototype("index_1, value_1, index_2, value_2", "dist")
+  .add_parameter("h1, h2", "array_like (1D)", "Histograms to compute the chi square distance for")
+  .add_parameter("index_1, index_2", "array_like (int, 1D)", "Indices of the sparse histograms value_1 and value_2")
+  .add_parameter("value_1, value_2", "array_like (1D)", "Sparse histograms to compute the chi square distance for")
+  .add_return("dist", "float", "The chi square distance value for the given histograms.")
+;
+
+static xbob::extension::FunctionDoc s_kullback_leibler = xbob::extension::FunctionDoc(
+    "kullback_leibler",
+    "Computes the Kullback-Leibler histogram divergence between the given histograms, which might be of singular dimension only.",
+    "The chi square distance is inspired by `link <http://www.informatik.uni-freiburg.de/~tipaldi/FLIRTLib/HistogramDistances_8hpp_source.html>`_ and computed as follows:\n"
+    ".. math:: dist(h_1,h_2) = \\sum_i (h_{1i} - h_{2i}) * \\log (h_{1i} / h_{2i})\n"
+    "The Kullback-Leibler divergence defines a distance metric, so lower values are better. "
+    "You can use this method in two different formats. "
+    "The first interface accepts non-sparse histograms. "
+    "The second interface accepts sparse histograms represented by indexes and values.\n"
+    ".. note:: Histograms are given as two matrices, one with the indexes and one with the data. All data points that for which no index exists are considered to be zero.\n"
+    ".. note:: In general, histogram intersection with sparse histograms needs more time to be computed."
+  )
+  .add_prototype("h1, h2", "dist")
+  .add_prototype("index_1, value_1, index_2, value_2", "dist")
+  .add_parameter("h1, h2", "array_like (1D)", "Histograms to compute the Kullback-Leibler divergence for")
+  .add_parameter("index_1, index_2", "array_like (int, 1D)", "Indices of the sparse histograms value_1 and value_2")
+  .add_parameter("value_1, value_2", "array_like (1D)", "Sparse histograms to compute the Kullback-Leibler divergence for")
+  .add_return("dist", "float", "The Kullback-Leibler divergence value for the given histograms.")
+;
+
+static xbob::extension::FunctionDoc s_linsolve = xbob::extension::FunctionDoc(
+  "linsolve",
+  "Solves the linear system :math:`Ax=b` and returns the result in :math:`x`.",
+  "This method uses LAPACK's ``dgesv`` generic solver. "
+  "You can use this method in two different formats. "
+  "The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
+  "The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
+  )
+  .add_prototype("A, b", "x")
+  .add_prototype("A, b, x")
+  .add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
+  .add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
+  .add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
+  .add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
+;
+
+static xbob::extension::FunctionDoc s_linsolve_nocheck = xbob::extension::FunctionDoc(
+  "linsolve_",
+  "Solves the linear system :math:`Ax=b` and returns the result in :math:`x`.",
+  ".. warning:: This variant does not perform any checks on the input matrices and is faster then :func:`linsolve`. "
+  "Use it when you are sure your input matrices sizes match.\n"
+  "This method uses LAPACK's ``dgesv`` generic solver. "
+  "You can use this method in two different formats. "
+  "The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
+  "The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution.\n"
+  )
+  .add_prototype("A, b", "x")
+  .add_prototype("A, b, x")
+  .add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
+  .add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
+  .add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
+  .add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
+;
+
+static xbob::extension::FunctionDoc s_linsolve_sympos = xbob::extension::FunctionDoc(
+  "linsolve_sympos",
+  "Solves the linear system :math:`Ax=b` and returns the result in :math:`x` for symmetric :math:`A` matrix.",
+  "This method uses LAPACK's ``dposv`` solver, assuming :math:`A` is a symmetric positive definite matrix. "
+  "You can use this method in two different formats. "
+  "The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
+  "The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
+  )
+  .add_prototype("A, b", "x")
+  .add_prototype("A, b, x")
+  .add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
+  .add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
+  .add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
+  .add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
+;
+
+static xbob::extension::FunctionDoc s_linsolve_sympos_nocheck = xbob::extension::FunctionDoc(
+  "linsolve_sympos_",
+  ".. warning:: This variant does not perform any checks on the input matrices and is faster then :func:`linsolve_sympos`. "
+  "Use it when you are sure your input matrices sizes match.\n"
+  "Solves the linear system :math:`Ax=b` and returns the result in :math:`x` for symmetric :math:`A` matrix.",
+  "This method uses LAPACK's ``dposv`` solver, assuming :math:`A` is a symmetric positive definite matrix. "
+  "You can use this method in two different formats. "
+  "The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
+  "The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
+  )
+  .add_prototype("A, b", "x")
+  .add_prototype("A, b, x")
+  .add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
+  .add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
+  .add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
+  .add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
+;
+
+static xbob::extension::FunctionDoc s_linsolve_cg_sympos = xbob::extension::FunctionDoc(
+  "linsolve_cg_sympos",
+  "Solves the linear system :math:`Ax=b` using conjugate gradients and returns the result in :math:`x` for symmetric :math:`A` matrix.",
+  "This method uses the conjugate gradient solver, assuming :math:`A` is a symmetric positive definite matrix. "
+  "You can use this method in two different formats. "
+  "The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
+  "The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
+  )
+  .add_prototype("A, b", "x")
+  .add_prototype("A, b, x")
+  .add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
+  .add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
+  .add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
+  .add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
+;
+
+static xbob::extension::FunctionDoc s_linsolve_cg_sympos_nocheck = xbob::extension::FunctionDoc(
+  "linsolve_cg_sympos_",
+  "Solves the linear system :math:`Ax=b` using conjugate gradients and returns the result in :math:`x` for symmetric :math:`A` matrix.",
+  ".. warning:: This variant does not perform any checks on the input matrices and is faster then :func:`linsolve_cg_sympos`. "
+  "Use it when you are sure your input matrices sizes match.\n"
+  "This method uses the conjugate gradient solver, assuming :math:`A` is a symmetric positive definite matrix. "
+  "You can use this method in two different formats. "
+  "The first interface accepts the matrices :math:`A` and :math:`b` returning :math:`x`. "
+  "The second one accepts a pre-allocated :math:`x` vector and sets it with the linear system solution."
+  )
+  .add_prototype("A, b", "x")
+  .add_prototype("A, b, x")
+  .add_parameter("A", "array_like (2D)", "The matrix :math:`A` of the linear system")
+  .add_parameter("b", "array_like (1D)", "The vector :math:`b` of the linear system")
+  .add_parameter("x", "array_like (1D)", "The result vector :math:`x`, as parameter")
+  .add_return("x", "array_like (1D)", "The result vector :math:`x`, as return value")
+;
+
+static xbob::extension::FunctionDoc s_pavx = xbob::extension::FunctionDoc(
+  "pavx",
+  "Applies the Pool-Adjacent-Violators Algorithm",
+  "Applies the Pool-Adjacent-Violators Algorithm to ``input``. "
+  "This is a simplified C++ port of the isotonic regression code made available at the `University of Bern website <http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html>`_.\n"
+  "You can use this method in two different formats. "
+  "The first interface accepts the ``input`` and ``output``. "
+  "The second one accepts the input array ``input`` and allocates a new ``output`` array, which is returned. "
+  )
+  .add_prototype("input, output")
+  .add_prototype("input", "output")
+  .add_parameter("input", "array_like (float, 1D)", "The input matrix for the PAV algorithm.")
+  .add_parameter("output", "array_like (float, 1D)", "The output matrix, must be of the same size as ``input``")
+  .add_return("output", "array_like (float, 1D)", "The output matrix; will be created in the same size as ``input``")
+;
+
+static xbob::extension::FunctionDoc s_pavx_nocheck = xbob::extension::FunctionDoc(
+  "pavx_",
+  "Applies the Pool-Adjacent-Violators Algorithm",
+  ".. warning:: This variant does not perform any checks on the input matrices and is faster then :func:`pavx`. "
+  "Use it when you are sure your input matrices sizes match.\n"
+  "Applies the Pool-Adjacent-Violators Algorithm to ``input``. "
+  "This is a simplified C++ port of the isotonic regression code made available at the `University of Bern website <http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html>`_.\n"
+  "You can use this method in two different formats. "
+  "The first interface accepts the ``input`` and ``output``. "
+  "The second one accepts the input array ``input`` and allocates a new ``output`` array, which is returned. "
+  )
+  .add_prototype("input, output")
+  .add_prototype("input", "output")
+  .add_parameter("input", "array_like (float, 1D)", "The input matrix for the PAV algorithm.")
+  .add_parameter("output", "array_like (float, 1D)", "The output matrix, must be of the same size as ``input``")
+  .add_return("output", "array_like (float, 1D)", "The output matrix; will be created in the same size as ``input``")
+;
+
+static xbob::extension::FunctionDoc s_pavx_width = xbob::extension::FunctionDoc(
+  "pavxWidth",
+  "Applies the Pool-Adjacent-Violators Algorithm and returns the width.",
+  "Applies the Pool-Adjacent-Violators Algorithm to ``input``. "
+  "This is a simplified C++ port of the isotonic regression code made available at the `University of Bern website <http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html>`_."
+  )
+  .add_prototype("input, output", "width")
+  .add_parameter("input", "array_like (float, 1D)", "The input matrix for the PAV algorithm.")
+  .add_parameter("output", "array_like (float, 1D)", "The output matrix, must be of the same size as ``input``")
+  .add_return("width", "array_like (uint64, 1D)", "The width matrix will be created in the same size as ``input``\n.. todo:: Explain, what width means in this case")
+;
+
+static xbob::extension::FunctionDoc s_pavx_width_height = xbob::extension::FunctionDoc(
+  "pavxWidthHeight",
+  "Applies the Pool-Adjacent-Violators Algorithm and returns the width and the height.",
+  "Applies the Pool-Adjacent-Violators Algorithm to ``input``. "
+  "This is a simplified C++ port of the isotonic regression code made available at the `University of Bern website <http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html>`_."
+  )
+  .add_prototype("input, output", "width, height")
+  .add_parameter("input", "array_like (float, 1D)", "The input matrix for the PAV algorithm.")
+  .add_parameter("output", "array_like (float, 1D)", "The output matrix, must be of the same size as ``input``")
+  .add_return("width", "array_like (uint64, 1D)", "The width matrix will be created in the same size as ``input``\n.. todo:: Explain, what width means in this case")
+  .add_return("height", "array_like (float, 1D)", "The height matrix will be created in the same size as ``input``\n.. todo:: Explain, what height means in this case")
+;
+
+static xbob::extension::FunctionDoc s_norminv = xbob::extension::FunctionDoc(
+  "norminv",
+  "Computes the inverse normal cumulative distribution",
+  "Computes the inverse normal cumulative distribution for a probability :math:`p`, given a distribution with mean :math:`\\mu` and standard deviation :math:`\\sigma`. "
+  "Reference: http://home.online.no/~pjacklam/notes/invnorm/"
+  )
+  .add_prototype("p, mu, sigma", "inv")
+  .add_parameter("p", "float", "The value to get the inverse distribution of, must lie in the range :math:`[0,1]`")
+  .add_parameter("mu", "float", "The mean :math:`\\mu` of the normal distribution")
+  .add_parameter("sigma", "float", "The standard deviation :math:`\\sigma` of the normal distribution")
+  .add_return("inv", "float", "The inverse of the normal distribution")
+;
+
+static xbob::extension::FunctionDoc s_normsinv = xbob::extension::FunctionDoc(
+  "normsinv",
+  "Computes the inverse normal cumulative distribution",
+  "Computes the inverse normal cumulative distribution for a probability :math:`p`, given a distribution with mean :math:`\\mu=0` and standard deviation :math:`\\sigma=1`. "
+  "It is equivalent as calling ``norminv(p, 0, 1)`` (see :func:`norminv`). "
+  "Reference: http://home.online.no/~pjacklam/notes/invnorm/"
+  )
+  .add_prototype("p", "inv")
+  .add_parameter("p", "float", "The value to get the inverse distribution of, must lie in the range :math:`[0,1]`")
+  .add_return("inv", "float", "The inverse of the normal distribution")
+;
+
+static xbob::extension::FunctionDoc s_scatter = xbob::extension::FunctionDoc(
+  "scatter",
+  "Computes scatter matrix of a 2D array.",
+  "Computes the scatter matrix of a 2D array *considering data is organized row-wise* (each sample is a row, each feature is a column). "
+  "The resulting array ``s`` is squared with extents equal to the number of columns in ``a``. "
+  "The resulting array ``m`` is a 1D array with the row means of ``a``. "
+  "This function supports many calling modes, but you should provide, at least, the input data matrix ``a``. "
+  "All non-provided arguments will be allocated internally and returned."
+  )
+  .add_prototype("a", "s, m")
+  .add_prototype("a, s", "m")
+  .add_prototype("a, m", "s")
+  .add_prototype("a, s, m")
+  .add_parameter("a", "array_like (float, 2D)", "The sample matrix, *considering data is organized row-wise* (each sample is a row, each feature is a column)")
+  .add_parameter("s", "array_like (float, 2D)", "The scatter matrix, squared with extents equal to the number of columns in ``a``")
+  .add_parameter("m", "array_like (float,1D)", "The mean matrix, with with the row means of ``a``")
+  .add_return("s", "array_like (float, 2D)", "The scatter matrix, squared with extents equal to the number of columns in ``a``")
+  .add_return("m", "array_like (float, 1D)", "The mean matrix, with with the row means of ``a``")
+;
+
+static xbob::extension::FunctionDoc s_scatter_nocheck = xbob::extension::FunctionDoc(
+  "scatter_",
+  "Computes scatter matrix of a 2D array.",
+  ".. warning:: This variant does not perform any checks on the input matrices and is faster then :func:`scatter`. "
+  "Use it when you are sure your input matrices sizes match.\n"
+  "Computes the scatter matrix of a 2D array *considering data is organized row-wise* (each sample is a row, each feature is a column). "
+  "The resulting array ``s`` is squared with extents equal to the number of columns in ``a``. "
+  "The resulting array ``m`` is a 1D array with the row means of ``a``. "
+  "This function supports many calling modes, but you should provide, at least, the input data matrix ``a``. "
+  "All non-provided arguments will be allocated internally and returned."
+  )
+  .add_prototype("a, s, m")
+  .add_parameter("a", "array_like (float, 2D)", "The sample matrix, *considering data is organized row-wise* (each sample is a row, each feature is a column)")
+  .add_parameter("s", "array_like (float, 2D)", "The scatter matrix, squared with extents equal to the number of columns in ``a``")
+  .add_parameter("m", "array_like (float,1D)", "The mean matrix, with with the row means of ``a``")
+;
+
+
+static xbob::extension::FunctionDoc s_scatters = xbob::extension::FunctionDoc(
+  "scatters",
+  "Computes :math:`S_w` and :math:`S_b` scatter matrices of a set of 2D arrays.",
+  "Computes the within-class :math:`S_w` and between-class :math:`S_b` scatter matrices of a set of 2D arrays considering data is organized row-wise (each sample is a row, each feature is a column), and each matrix contains data of one class. "
+  "Computes the scatter matrix of a 2D array *considering data is organized row-wise* (each sample is a row, each feature is a column). "
+  "The implemented strategy is:\n"
+  "1. Evaluate the overall mean (``m``), class means (:math:`m_k`) and the  total class counts (:math:`N`).\n"
+  "2. Evaluate ``sw`` and ``sb`` using normal loops.\n"
+  "Note that in this implementation, ``sw`` and ``sb`` will be normalized by N-1 (number of samples) and K (number of classes). "
+  "This procedure makes the eigen values scaled by (N-1)/K, effectively increasing their values. "
+  "The main motivation for this normalization are numerical precision concerns with the increasing number of samples causing a rather large :math:`S_w` matrix. "
+  "A normalization strategy mitigates this problem. "
+  "The eigen vectors will see no effect on this normalization as they are normalized in the euclidean sense (:math:`||a|| = 1`) so that does not change those.\n"
+  "This function supports many calling modes, but you should provide, at least, the input ``data``. "
+  "All non-provided arguments will be allocated internally and returned."
+  )
+  .add_prototype("data", "sw, sb, m")
+  .add_prototype("data, sw, sb", "m")
+  .add_prototype("data, sw, sb, m")
+  .add_parameter("data", "[array_like (float, 2D)]", "The list of sample matrices. "
+      "In each sample matrix the data is organized row-wise (each sample is a row, each feature is a column). "
+      "Each matrix stores the data of a particular class. "
+      "**Every matrix in ``data`` must have exactly the same number of columns.**")
+  .add_parameter("sw", "array_like (float, 2D)", "The within-class scatter matrix :math:`S_w`, squared with extents equal to the number of columns in ``data``")
+  .add_parameter("sb", "array_like (float, 2D)", "The between-class scatter matrix :math:`S_b`, squared with extents equal to the number of columns in ``data``")
+  .add_parameter("m", "array_like (float,1D)", "The mean matrix, representing the ensemble mean with no prior (i.e., biased towards classes with more samples)")
+  .add_return("sw", "array_like (float, 2D)", "The within-class scatter matrix :math:`S_w`")
+  .add_return("sb", "array_like (float, 2D)", "The between-class scatter matrix :math:`S_b`")
+  .add_return("m", "array_like (float, 1D)", "The mean matrix, representing the ensemble mean with no prior (i.e., biased towards classes with more samples)")
+;
+
+static xbob::extension::FunctionDoc s_scatters_nocheck = xbob::extension::FunctionDoc(
+  "scatters_",
+  "Computes :math:`S_w` and :math:`S_b` scatter matrices of a set of 2D arrays.",
+  ".. warning:: This variant does not perform any checks on the input matrices and is faster then :func:`scatters`. "
+  "Use it when you are sure your input matrices sizes match.\n"
+  "For a detailed description of the function, please see :func:`scatters`."
+  )
+  .add_prototype("data, sw, sb, m")
+  .add_prototype("data, sw, sb")
+  .add_parameter("data", "[array_like (float, 2D)]", "The list of sample matrices. "
+      "In each sample matrix the data is organized row-wise (each sample is a row, each feature is a column). "
+      "Each matrix stores the data of a particular class. "
+      "**Every matrix in ``data`` must have exactly the same number of columns.**")
+  .add_parameter("sw", "array_like (float, 2D)", "The within-class scatter matrix :math:`S_w`, squared with extents equal to the number of columns in ``data``")
+  .add_parameter("sb", "array_like (float, 2D)", "The between-class scatter matrix :math:`S_b`, squared with extents equal to the number of columns in ``data``")
+  .add_parameter("m", "array_like (float,1D)", "The mean matrix, representing the ensemble mean with no prior (i.e., biased towards classes with more samples)")
+;
+
 
 static PyMethodDef module_methods[] = {
     {
-      s_histogram_intersection_str,
+      s_histogram_intersection.name(),
       (PyCFunction)py_histogram_intersection,
       METH_VARARGS|METH_KEYWORDS,
-      s_histogram_intersection_doc
+      s_histogram_intersection.doc()
     },
     {
-      s_chi_square_str,
+      s_chi_square.name(),
       (PyCFunction)py_chi_square,
       METH_VARARGS|METH_KEYWORDS,
-      s_chi_square_doc
+      s_chi_square.doc()
     },
     {
-      s_kullback_leibler_str,
+      s_kullback_leibler.name(),
       (PyCFunction)py_kullback_leibler,
       METH_VARARGS|METH_KEYWORDS,
-      s_kullback_leibler_doc
+      s_kullback_leibler.doc()
     },
     {
-      s_linsolve_str,
+      s_linsolve.name(),
       (PyCFunction)py_linsolve,
       METH_VARARGS|METH_KEYWORDS,
-      s_linsolve_doc
+      s_linsolve.doc()
     },
     {
-      s_linsolve_nocheck_str,
+      s_linsolve_nocheck.name(),
       (PyCFunction)py_linsolve_nocheck,
       METH_VARARGS|METH_KEYWORDS,
-      s_linsolve_nocheck_doc
+      s_linsolve_nocheck.doc()
     },
     {
-      s_linsolve_sympos_str,
+      s_linsolve_sympos.name(),
       (PyCFunction)py_linsolve_sympos,
       METH_VARARGS|METH_KEYWORDS,
-      s_linsolve_sympos_doc
+      s_linsolve_sympos.doc()
     },
     {
-      s_linsolve_sympos_nocheck_str,
+      s_linsolve_sympos_nocheck.name(),
       (PyCFunction)py_linsolve_sympos_nocheck,
       METH_VARARGS|METH_KEYWORDS,
-      s_linsolve_sympos_nocheck_doc
+      s_linsolve_sympos_nocheck.doc()
     },
     {
-      s_linsolve_cg_sympos_str,
+      s_linsolve_cg_sympos.name(),
       (PyCFunction)py_linsolve_cg_sympos,
       METH_VARARGS|METH_KEYWORDS,
-      s_linsolve_cg_sympos_doc
+      s_linsolve_cg_sympos.doc()
     },
     {
-      s_linsolve_cg_sympos_nocheck_str,
+      s_linsolve_cg_sympos_nocheck.name(),
       (PyCFunction)py_linsolve_cg_sympos_nocheck,
       METH_VARARGS|METH_KEYWORDS,
-      s_linsolve_cg_sympos_nocheck_doc
+      s_linsolve_cg_sympos_nocheck.doc()
     },
     {
-      s_pavx_str,
+      s_pavx.name(),
       (PyCFunction)py_pavx,
       METH_VARARGS|METH_KEYWORDS,
-      s_pavx_doc
+      s_pavx.doc()
     },
     {
-      s_pavx_nocheck_str,
+      s_pavx_nocheck.name(),
       (PyCFunction)py_pavx_nocheck,
       METH_VARARGS|METH_KEYWORDS,
-      s_pavx_nocheck_doc
+      s_pavx_nocheck.doc()
     },
     {
-      s_pavx_width_str,
+      s_pavx_width.name(),
       (PyCFunction)py_pavx_width,
       METH_VARARGS|METH_KEYWORDS,
-      s_pavx_width_doc
+      s_pavx_width.doc()
     },
     {
-      s_pavx_width_height_str,
+      s_pavx_width_height.name(),
       (PyCFunction)py_pavx_width_height,
       METH_VARARGS|METH_KEYWORDS,
-      s_pavx_width_height_doc
+      s_pavx_width_height.doc()
     },
     {
-      s_norminv_str,
+      s_norminv.name(),
       (PyCFunction)py_norminv,
       METH_VARARGS|METH_KEYWORDS,
-      s_norminv_doc
+      s_norminv.doc()
     },
     {
-      s_normsinv_str,
+      s_normsinv.name(),
       (PyCFunction)py_normsinv,
       METH_VARARGS|METH_KEYWORDS,
-      s_normsinv_doc
+      s_normsinv.doc()
     },
     {
-      s_scatter_str,
+      s_scatter.name(),
       (PyCFunction)py_scatter,
       METH_VARARGS|METH_KEYWORDS,
-      s_scatter_doc
+      s_scatter.doc()
     },
     {
-      s_scatter_nocheck_str,
+      s_scatter_nocheck.name(),
       (PyCFunction)py_scatter_nocheck,
       METH_VARARGS|METH_KEYWORDS,
-      s_scatter_nocheck_doc
+      s_scatter_nocheck.doc()
     },
     {
-      s_scatters_str,
+      s_scatters.name(),
       (PyCFunction)py_scatters,
       METH_VARARGS|METH_KEYWORDS,
-      s_scatters_doc
+      s_scatters.doc()
     },
     {
-      s_scatters_nocheck_str,
+      s_scatters_nocheck.name(),
       (PyCFunction)py_scatters_nocheck,
       METH_VARARGS|METH_KEYWORDS,
-      s_scatters_nocheck_doc
+      s_scatters_nocheck.doc()
     },
     {0}  /* Sentinel */
 };
-- 
GitLab