diff --git a/setup.py b/setup.py
index c14fbbcdcfc710cf75080152bf43ac7d0f83d52f..9ad58995a366e3d9f8a37ec702c7bb27d4ea9513 100644
--- a/setup.py
+++ b/setup.py
@@ -44,6 +44,7 @@ setup(
         [
           "xbob/math/histogram.cpp",
           "xbob/math/linsolve.cpp",
+          "xbob/math/pavx.cpp",
           "xbob/math/main.cpp",
           ],
         packages = packages,
diff --git a/xbob/math/histogram.cpp b/xbob/math/histogram.cpp
index aa257cf8dbd88e6236e5828c4f2d2c57e04f2ca7..21e22303b0a745469b175b1d7364924e097854e3 100644
--- a/xbob/math/histogram.cpp
+++ b/xbob/math/histogram.cpp
@@ -43,32 +43,42 @@ static PyObject* py_histogram_intersection_1
 
   PyObject* retval = 0;
 
-  switch(h1->type_num) {
+  try {
 
-    case NPY_UINT8:
-      retval = PyBlitzArrayCxx_FromCScalar(bob::math::histogram_intersection(*PyBlitzArrayCxx_AsBlitz<uint8_t,1>(h1), *PyBlitzArrayCxx_AsBlitz<uint8_t,1>(h2)));
-      break;
+    switch(h1->type_num) {
 
-    case NPY_UINT16:
-      retval = PyBlitzArrayCxx_FromCScalar(bob::math::histogram_intersection(*PyBlitzArrayCxx_AsBlitz<uint16_t,1>(h1), *PyBlitzArrayCxx_AsBlitz<uint16_t,1>(h2)));
-      break;
+      case NPY_UINT8:
+        retval = PyBlitzArrayCxx_FromCScalar(bob::math::histogram_intersection(*PyBlitzArrayCxx_AsBlitz<uint8_t,1>(h1), *PyBlitzArrayCxx_AsBlitz<uint8_t,1>(h2)));
+        break;
 
-    case NPY_INT32:
-      retval = PyBlitzArrayCxx_FromCScalar(bob::math::histogram_intersection(*PyBlitzArrayCxx_AsBlitz<int32_t,1>(h1), *PyBlitzArrayCxx_AsBlitz<int32_t,1>(h2)));
-      break;
+      case NPY_UINT16:
+        retval = PyBlitzArrayCxx_FromCScalar(bob::math::histogram_intersection(*PyBlitzArrayCxx_AsBlitz<uint16_t,1>(h1), *PyBlitzArrayCxx_AsBlitz<uint16_t,1>(h2)));
+        break;
 
-    case NPY_INT64:
-      retval = PyBlitzArrayCxx_FromCScalar(bob::math::histogram_intersection(*PyBlitzArrayCxx_AsBlitz<int64_t,1>(h1), *PyBlitzArrayCxx_AsBlitz<int64_t,1>(h2)));
-      break;
+      case NPY_INT32:
+        retval = PyBlitzArrayCxx_FromCScalar(bob::math::histogram_intersection(*PyBlitzArrayCxx_AsBlitz<int32_t,1>(h1), *PyBlitzArrayCxx_AsBlitz<int32_t,1>(h2)));
+        break;
 
-    case NPY_FLOAT64:
-      retval = PyBlitzArrayCxx_FromCScalar(bob::math::histogram_intersection(*PyBlitzArrayCxx_AsBlitz<double,1>(h1), *PyBlitzArrayCxx_AsBlitz<double,1>(h2)));
-      break;
+      case NPY_INT64:
+        retval = PyBlitzArrayCxx_FromCScalar(bob::math::histogram_intersection(*PyBlitzArrayCxx_AsBlitz<int64_t,1>(h1), *PyBlitzArrayCxx_AsBlitz<int64_t,1>(h2)));
+        break;
 
-    default:
-      PyErr_Format(PyExc_TypeError, "Histogram intersection currently not implemented for type '%s'", PyBlitzArray_TypenumAsString(h1->type_num));
+      case NPY_FLOAT64:
+        retval = PyBlitzArrayCxx_FromCScalar(bob::math::histogram_intersection(*PyBlitzArrayCxx_AsBlitz<double,1>(h1), *PyBlitzArrayCxx_AsBlitz<double,1>(h2)));
+        break;
+
+      default:
+        PyErr_Format(PyExc_TypeError, "Histogram intersection currently not implemented for type '%s'", PyBlitzArray_TypenumAsString(h1->type_num));
+
+    }
 
   }
+  catch (std::exception& e) {
+    PyErr_SetString(PyExc_RuntimeError, e.what());
+  }
+  catch (...) {
+    PyErr_SetString(PyExc_RuntimeError, "histogram intersection failed: unknown exception caught");
+  }
 
   Py_DECREF(h1);
   Py_DECREF(h2);
@@ -178,36 +188,46 @@ static PyObject* py_histogram_intersection_2(PyObject*, PyObject* args, PyObject
 
   PyObject* retval = 0;
 
-  switch(index1->type_num) {
+  try {
 
-    case NPY_UINT8:
-      retval = py_histogram_intersection_2_inner<uint8_t>(index1, value1,
-          index2, value2);
-      break;
+    switch(index1->type_num) {
 
-    case NPY_UINT16:
-      retval = py_histogram_intersection_2_inner<uint16_t>(index1, value1,
-          index2, value2);
-      break;
+      case NPY_UINT8:
+        retval = py_histogram_intersection_2_inner<uint8_t>(index1, value1,
+            index2, value2);
+        break;
 
-    case NPY_INT32:
-      retval = py_histogram_intersection_2_inner<int32_t>(index1, value1,
-          index2, value2);
-      break;
+      case NPY_UINT16:
+        retval = py_histogram_intersection_2_inner<uint16_t>(index1, value1,
+            index2, value2);
+        break;
 
-    case NPY_INT64:
-      retval = py_histogram_intersection_2_inner<int64_t>(index1, value1,
-          index2, value2);
-      break;
+      case NPY_INT32:
+        retval = py_histogram_intersection_2_inner<int32_t>(index1, value1,
+            index2, value2);
+        break;
 
-    case NPY_FLOAT64:
-      retval = py_histogram_intersection_2_inner<double>(index1, value1,
-          index2, value2);
-      break;
+      case NPY_INT64:
+        retval = py_histogram_intersection_2_inner<int64_t>(index1, value1,
+            index2, value2);
+        break;
 
-    default:
-      PyErr_Format(PyExc_TypeError, "Histogram intersection currently not implemented for index type '%s'", PyBlitzArray_TypenumAsString(index1->type_num));
+      case NPY_FLOAT64:
+        retval = py_histogram_intersection_2_inner<double>(index1, value1,
+            index2, value2);
+        break;
 
+      default:
+        PyErr_Format(PyExc_TypeError, "Histogram intersection currently not implemented for index type '%s'", PyBlitzArray_TypenumAsString(index1->type_num));
+
+    }
+
+  }
+  catch (std::exception& e) {
+    PyErr_SetString(PyExc_RuntimeError, e.what());
+  }
+  catch (...) {
+    PyErr_SetString(PyExc_RuntimeError, "histogram intersection failed: unknown exception caught");
   }
 
   Py_DECREF(index1);
diff --git a/xbob/math/linsolve.h b/xbob/math/linsolve.h
index 09b98ff1f17829617987888f8c129e9194ccad5f..515159fe640020e38bb6362dc814af52158d26bd 100644
--- a/xbob/math/linsolve.h
+++ b/xbob/math/linsolve.h
@@ -2,7 +2,7 @@
  * @author Andre Anjos <andre.anjos@idiap.ch>
  * @date Wed  4 Dec 15:26:54 2013 
  *
- * @brief 
+ * @brief Declaration of components relevant for main.cpp
  */
 
 #include <Python.h>
diff --git a/xbob/math/main.cpp b/xbob/math/main.cpp
index b888098ac2827791ea15af64b16694d572f01f52..a4bb4bbaacc4c97d596764d27a0c99a1e40bde32 100644
--- a/xbob/math/main.cpp
+++ b/xbob/math/main.cpp
@@ -15,6 +15,7 @@
 
 #include "histogram.h"
 #include "linsolve.h"
+#include "pavx.h"
 
 PyDoc_STRVAR(s_histogram_intersection_str, "histogram_intersection");
 PyDoc_STRVAR(s_histogram_intersection_doc,
@@ -85,9 +86,11 @@ linsolve_(A, x, b) -> None\n\
 Solves the linear system :py:math:`Ax=b` and returns the result in ``x``.\n\
 This method uses LAPACK's ``dgesv`` generic solver.\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 are well-behaved (contiguous, c-style, memory aligned).\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\
@@ -121,9 +124,11 @@ Solves the linear system :py: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\
-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 are well-behaved (contiguous, c-style, memory aligned).\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\
@@ -157,9 +162,11 @@ Solves the linear system :py: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\
-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 are well-behaved (contiguous, c-style, memory aligned).\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\
@@ -168,6 +175,70 @@ 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\
+");
+
 static PyMethodDef module_methods[] = {
     {
       s_histogram_intersection_str,
@@ -223,6 +294,30 @@ static PyMethodDef module_methods[] = {
       METH_VARARGS|METH_KEYWORDS,
       s_linsolve_cg_sympos_nocheck_doc
     },
+    {
+      s_pavx_str,
+      (PyCFunction)py_pavx,
+      METH_VARARGS|METH_KEYWORDS,
+      s_pavx_doc
+    },
+    {
+      s_pavx_nocheck_str,
+      (PyCFunction)py_pavx_nocheck,
+      METH_VARARGS|METH_KEYWORDS,
+      s_pavx_nocheck_doc
+    },
+    {
+      s_pavx_width_str,
+      (PyCFunction)py_pavx_width,
+      METH_VARARGS|METH_KEYWORDS,
+      s_pavx_width_doc
+    },
+    {
+      s_pavx_width_height_str,
+      (PyCFunction)py_pavx_width_height,
+      METH_VARARGS|METH_KEYWORDS,
+      s_pavx_width_height_doc
+    },
     {0}  /* Sentinel */
 };
 
diff --git a/xbob/math/pavx.cpp b/xbob/math/pavx.cpp
index 491791552ba50fd33e92fe11d0ccdf1b334df3b7..067ff048046a08571246fa8b93a1fc713f1eb8aa 100644
--- a/xbob/math/pavx.cpp
+++ b/xbob/math/pavx.cpp
@@ -4,86 +4,232 @@
  * @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
  *
  * @brief Binds the Pool-Adjacent-Violators Algorithm
- *
- * Copyright (C) 2011-2013 Idiap Research Institute, Martigny, Switzerland
- * 
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, version 3 of the License.
- * 
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- * 
- * You should have received a copy of the GNU General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include "pavx.h"
+#include <xbob.blitz/cppapi.h>
 #include "bob/math/pavx.h"
-
-#include "bob/python/ndarray.h"
 #include "bob/core/cast.h"
 
-using namespace boost::python;
+PyObject* py_pavx (PyObject*, PyObject* args, PyObject* kwds) {
 
-static const char* C_PAVX_DOC = "Applies the Pool-Adjacent-Violators Algorithm. The input and output arrays should have the same size. This is a simplified C++ port of the isotonic regression code made available at http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html.";
-static const char* C_PAVX__DOC = "Applies the Pool-Adjacent-Violators Algorithm. The input and output arrays should have the same size. Arguments are not checked! This is a simplified C++ port of the isotonic regression code made available at http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html.";
-static const char* P_PAVX_DOC = "Applies the Pool-Adjacent-Violators Algorithm. The output array is allocated and returned. This is a simplified C++ port of the isotonic regression code made available at http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html.";
-static const char* P_PAVXWIDTH_DOC = "Applies the Pool-Adjacent-Violators Algorithm. The input and output arrays should have the same size. The width array is returned. This is a simplified C++ port of the isotonic regression code made available at http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html.";
-static const char* P_PAVXWIDTHHEIGHT_DOC = "Applies the Pool-Adjacent-Violators Algorithm. The input and output arrays should have the same size. The width and height arrays are returned. This is a simplified C++ port of the isotonic regression code made available at http://www.imsv.unibe.ch/content/staff/personalhomepages/duembgen/software/isotonicregression/index_eng.html.";
+  /* Parses input arguments in a single shot */
+  static const char* const_kwlist[] = { "input", "output", 0 /* Sentinel */ };
+  static char** kwlist = const_cast<char**>(const_kwlist);
 
-static void c_pavx_(bob::python::const_ndarray y, bob::python::ndarray ghat) 
-{
-  blitz::Array<double,1> ghat_ = ghat.bz<double,1>();
-  bob::math::pavx_(y.bz<double,1>(), ghat_);
-}
+  PyBlitzArrayObject* input = 0;
+  PyBlitzArrayObject* output = 0;
 
-static void c_pavx(bob::python::const_ndarray y, bob::python::ndarray ghat) 
-{
-  blitz::Array<double,1> ghat_ = ghat.bz<double,1>();
-  bob::math::pavx(y.bz<double,1>(), ghat_);
-}
+  if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&|O&",
+        kwlist, 
+        &PyBlitzArray_Converter, &input, 
+        &PyBlitzArray_OutputConverter, &output
+        )) return 0;
 
-static object p_pavx(bob::python::const_ndarray y) 
-{
-  const bob::core::array::typeinfo& info = y.type();
-  blitz::Array<double,1> y_ = y.bz<double,1>();
-  bob::python::ndarray ghat(bob::core::array::t_float64, info.shape[0]);
-  blitz::Array<double,1> ghat_ = ghat.bz<double,1>();
-  bob::math::pavx(y.bz<double,1>(), ghat_);
-  return ghat.self();
-}
+  // can only handle 1D arrays
+  if (input->ndim != 1 || (output && output->ndim != 1)) {
+    PyErr_SetString(PyExc_TypeError, "input and output arrays should be one-dimensional");
+    Py_DECREF(input);
+    Py_XDECREF(output);
+    return 0;
+  }
+
+  // can only handle float arrays
+  if (input->type_num != NPY_FLOAT64 || (output && output->type_num != NPY_FLOAT64)) {
+    PyErr_SetString(PyExc_TypeError, "input and output arrays data types should be float (i.e. `numpy.float64' equivalents)");
+    Py_DECREF(input);
+    Py_XDECREF(output);
+    return 0;
+  }
+
+  // if output was not provided by user
+  bool returns_output = false;
+  if (!output) {
+    output = (PyBlitzArrayObject*)PyBlitzArray_SimpleNew(NPY_FLOAT64, input->ndim, input->shape);
+    if (!output) {
+      Py_DECREF(input);
+      return 0;
+    }
+    returns_output = true;
+  }
+
+  try {
+    bob::math::pavx(*PyBlitzArrayCxx_AsBlitz<double,1>(input),
+        *PyBlitzArrayCxx_AsBlitz<double,1>(output));
+  }
+  catch (std::exception& e) {
+    PyErr_SetString(PyExc_RuntimeError, e.what());
+  }
+  catch (...) {
+    PyErr_SetString(PyExc_RuntimeError, "pavx failed: unknown exception caught");
+  }
+
+  Py_DECREF(input);
+  if (returns_output) return (PyObject*)output;
+  Py_DECREF(output);
+  Py_RETURN_NONE;
 
-static object p_pavxWidth(bob::python::const_ndarray y, bob::python::ndarray ghat) 
-{
-  blitz::Array<double,1> ghat_ = ghat.bz<double,1>();
-  blitz::Array<size_t,1> w_ = bob::math::pavxWidth(y.bz<double,1>(), ghat_);
-  bob::python::ndarray width(bob::core::array::t_uint64, w_.extent(0));
-  blitz::Array<uint64_t,1> width_ = width.bz<uint64_t,1>();
-  width_ = bob::core::array::cast<uint64_t>(w_); 
-  return width.self();
 }
 
-static object p_pavxWidthHeight(bob::python::const_ndarray y, bob::python::ndarray ghat) 
-{
-  blitz::Array<double,1> ghat_ = ghat.bz<double,1>();
-  std::pair<blitz::Array<size_t,1>,blitz::Array<double,1> > pair = bob::math::pavxWidthHeight(y.bz<double,1>(), ghat_);
-  bob::python::ndarray width(bob::core::array::t_uint64, pair.first.extent(0));
-  blitz::Array<uint64_t,1> width_ = width.bz<uint64_t,1>();
-  width_ = bob::core::array::cast<uint64_t>(pair.first);
-  bob::python::ndarray height(bob::core::array::t_float64, pair.second.extent(0));
-  blitz::Array<double,1> height_ = height.bz<double,1>();
-  height_ = pair.second;
-  return make_tuple(width, height);
+PyObject* py_pavx_nocheck (PyObject*, PyObject* args, PyObject* kwds) {
+
+  /* Parses input arguments in a single shot */
+  static const char* const_kwlist[] = { "input", "output", 0 /* Sentinel */ };
+  static char** kwlist = const_cast<char**>(const_kwlist);
+
+  PyBlitzArrayObject* input = 0;
+  PyBlitzArrayObject* output = 0;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&O&",
+        kwlist, 
+        &PyBlitzArray_Converter, &input, 
+        &PyBlitzArray_OutputConverter, &output
+        )) return 0;
+
+  // can only handle 1D arrays
+  if (input->ndim != 1 || output->ndim != 1) {
+    PyErr_SetString(PyExc_TypeError, "input and output arrays should be one-dimensional");
+    Py_DECREF(input);
+    Py_DECREF(output);
+    return 0;
+  }
+
+  // can only handle float arrays
+  if (input->type_num != NPY_FLOAT64 || output->type_num != NPY_FLOAT64) {
+    PyErr_SetString(PyExc_TypeError, "input and output arrays data types should be float (i.e. `numpy.float64' equivalents)");
+    Py_DECREF(input);
+    Py_DECREF(output);
+    return 0;
+  }
+
+  try {
+    bob::math::pavx_(*PyBlitzArrayCxx_AsBlitz<double,1>(input),
+        *PyBlitzArrayCxx_AsBlitz<double,1>(output));
+  }
+  catch (std::exception& e) {
+    PyErr_SetString(PyExc_RuntimeError, e.what());
+  }
+  catch (...) {
+    PyErr_SetString(PyExc_RuntimeError, "pavx failed: unknown exception caught");
+  }
+
+  Py_DECREF(input);
+  Py_DECREF(output);
+  Py_RETURN_NONE;
 }
 
-void bind_math_pavx()
-{
-  def("pavx_", &c_pavx, (arg("input"), arg("output")), C_PAVX__DOC);
-  def("pavx", &c_pavx_, (arg("input"), arg("output")), C_PAVX_DOC);
-  def("pavx", &p_pavx, (arg("input")), P_PAVX_DOC);
-  def("pavxWidth", &p_pavxWidth, (arg("input"), arg("output")), P_PAVXWIDTH_DOC);
-  def("pavxWidthHeight", &p_pavxWidthHeight, (arg("input"), arg("output")), P_PAVXWIDTHHEIGHT_DOC);
+PyObject* py_pavx_width (PyObject*, PyObject* args, PyObject* kwds) {
+
+  /* Parses input arguments in a single shot */
+  static const char* const_kwlist[] = { "input", "output", 0 /* Sentinel */ };
+  static char** kwlist = const_cast<char**>(const_kwlist);
+
+  PyBlitzArrayObject* input = 0;
+  PyBlitzArrayObject* output = 0;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&O&",
+        kwlist,
+        &PyBlitzArray_Converter, &input,
+        &PyBlitzArray_OutputConverter, &output
+        )) return 0;
+
+  // can only handle 1D arrays
+  if (input->ndim != 1 || output->ndim != 1) {
+    PyErr_SetString(PyExc_TypeError, "input and output arrays should be one-dimensional");
+    Py_DECREF(input);
+    Py_DECREF(output);
+    return 0;
+  }
+
+  // can only handle float arrays
+  if (input->type_num != NPY_FLOAT64 || output->type_num != NPY_FLOAT64) {
+    PyErr_SetString(PyExc_TypeError, "input and output arrays data types should be float (i.e. `numpy.float64' equivalents)");
+    Py_DECREF(input);
+    Py_DECREF(output);
+    return 0;
+  }
+
+  PyObject* retval = 0;
+
+  try {
+    blitz::Array<size_t,1> width =
+      bob::math::pavxWidth(*PyBlitzArrayCxx_AsBlitz<double,1>(input),
+          *PyBlitzArrayCxx_AsBlitz<double,1>(output));
+    blitz::Array<uint64_t,1> uwidth = bob::core::array::cast<uint64_t>(width);
+    retval = PyBlitzArrayCxx_NewFromArray(uwidth);
+  }
+  catch (std::exception& e) {
+    PyErr_SetString(PyExc_RuntimeError, e.what());
+  }
+  catch (...) {
+    PyErr_SetString(PyExc_RuntimeError, "pavx failed: unknown exception caught");
+  }
+
+  Py_DECREF(input);
+  Py_DECREF(output);
+  return retval;
 }
 
+PyObject* py_pavx_width_height (PyObject*, PyObject* args, PyObject* kwds) {
+
+  /* Parses input arguments in a single shot */
+  static const char* const_kwlist[] = { "input", "output", 0 /* Sentinel */ };
+  static char** kwlist = const_cast<char**>(const_kwlist);
+
+  PyBlitzArrayObject* input = 0;
+  PyBlitzArrayObject* output = 0;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&O&",
+        kwlist, 
+        &PyBlitzArray_Converter, &input, 
+        &PyBlitzArray_OutputConverter, &output
+        )) return 0;
+
+  // can only handle 1D arrays
+  if (input->ndim != 1 || output->ndim != 1) {
+    PyErr_SetString(PyExc_TypeError, "input and output arrays should be one-dimensional");
+    Py_DECREF(input);
+    Py_DECREF(output);
+    return 0;
+  }
+
+  // can only handle float arrays
+  if (input->type_num != NPY_FLOAT64 || output->type_num != NPY_FLOAT64) {
+    PyErr_SetString(PyExc_TypeError, "input and output arrays data types should be float (i.e. `numpy.float64' equivalents)");
+    Py_DECREF(input);
+    Py_DECREF(output);
+    return 0;
+  }
+
+  PyObject* width = 0;
+  PyObject* height = 0;
+
+  try {
+    std::pair<blitz::Array<size_t,1>,blitz::Array<double,1>> width_height =
+      bob::math::pavxWidthHeight(*PyBlitzArrayCxx_AsBlitz<double,1>(input),
+          *PyBlitzArrayCxx_AsBlitz<double,1>(output));
+    blitz::Array<uint64_t,1> uwidth = bob::core::array::cast<uint64_t>(width_height.first);
+    width = PyBlitzArrayCxx_NewFromArray(uwidth);
+    if (width) {
+      height = PyBlitzArrayCxx_NewFromArray(width_height.second);
+      if (!height) Py_CLEAR(width);
+    }
+  }
+  catch (std::exception& e) {
+    PyErr_SetString(PyExc_RuntimeError, e.what());
+  }
+  catch (...) {
+    PyErr_SetString(PyExc_RuntimeError, "pavx failed: unknown exception caught");
+  }
+
+  Py_DECREF(input);
+  Py_DECREF(output);
+
+  if (!height) return 0;
+
+  // creates the return pair and returns
+  PyObject* retval = PyTuple_New(2);
+  PyTuple_SET_ITEM(retval, 0, width);
+  PyTuple_SET_ITEM(retval, 1, height);
+  return retval;
+}
diff --git a/xbob/math/pavx.h b/xbob/math/pavx.h
new file mode 100644
index 0000000000000000000000000000000000000000..c9e1b9082f1f2cd62bdc5805b2d55adb33723353
--- /dev/null
+++ b/xbob/math/pavx.h
@@ -0,0 +1,13 @@
+/**
+ * @author Andre Anjos <andre.anjos@idiap.ch>
+ * @date Wed  4 Dec 17:46:18 2013 
+ *
+ * @brief Declaration of components relevant for main.cpp
+ */
+
+#include <Python.h>
+
+PyObject* py_pavx(PyObject*, PyObject* args, PyObject* kwds);
+PyObject* py_pavx_nocheck(PyObject*, PyObject* args, PyObject* kwds);
+PyObject* py_pavx_width(PyObject*, PyObject* args, PyObject* kwds);
+PyObject* py_pavx_width_height(PyObject*, PyObject* args, PyObject* kwds);
diff --git a/xbob/math/stats.cpp b/xbob/math/stats.cpp
index fa0c97abc77c9a72297017e6f9d8ee69993d8ee1..32ce880baa6926fb6cb82209cde4e3aa46575806 100644
--- a/xbob/math/stats.cpp
+++ b/xbob/math/stats.cpp
@@ -6,18 +6,6 @@
  * @brief Python bindings to statistical methods
  *
  * Copyright (C) 2011-2013 Idiap Research Institute, Martigny, Switzerland
- * 
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, version 3 of the License.
- * 
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- * 
- * You should have received a copy of the GNU General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
 #include <boost/python.hpp>
diff --git a/xbob/math/test/test_pavx.py b/xbob/math/test/test_pavx.py
index 905ccdb4c60ee46b4484a3f84ade005f607c62a0..ffe0c4b3dfd990620a61acad009438a5445e085d 100644
--- a/xbob/math/test/test_pavx.py
+++ b/xbob/math/test/test_pavx.py
@@ -12,17 +12,16 @@ import os, sys
 from xbob.math import pavx, pavx_, pavxWidth, pavxWidthHeight
 import numpy
 
-def pavx_internal_test(y, ghat_ref, w_ref, h_ref):
+def pavx_check(y, ghat_ref, w_ref, h_ref):
   """Make a full test for a given sample"""
 
-  ghat = numpy.ndarray(dtype=numpy.float64, shape=y.shape)
-  pavx(y, ghat)
+  ghat = pavx(y)
   assert numpy.all(numpy.abs(ghat - ghat_ref) < 1e-4)
   pavx_(y, ghat)
   assert numpy.all(numpy.abs(ghat - ghat_ref) < 1e-4)
   w=pavxWidth(y, ghat)
-  assert numpy.all(numpy.abs(ghat - ghat_ref) < 1e-4)
   assert numpy.all(numpy.abs(w - w_ref) < 1e-4)
+  assert numpy.all(numpy.abs(ghat - ghat_ref) < 1e-4)
   ret=pavxWidthHeight(y, ghat)
   assert numpy.all(numpy.abs(ghat - ghat_ref) < 1e-4)
   assert numpy.all(numpy.abs(ret[0] - w_ref) < 1e-4)
@@ -45,7 +44,7 @@ def test_pavx_sample1():
     118.1671, 138.3151, 141.9755, 145.7352, 157.9881,
     168.6932, 175.2756])
 
-  pavx_internal_test(y, ghat_ref, w_ref, h_ref)
+  pavx_check(y, ghat_ref, w_ref, h_ref)
 
 def test_pavx_sample2():
 
@@ -63,4 +62,4 @@ def test_pavx_sample2():
     98.5769, 102.3841, 132.1742, 140.2783, 150.7383,
     154.7762,  180.8819])
 
-  pavx_internal_test(y, ghat_ref, w_ref, h_ref)
+  pavx_check(y, ghat_ref, w_ref, h_ref)