Skip to content
Snippets Groups Projects
Commit 77636441 authored by André Anjos's avatar André Anjos :speech_balloon:
Browse files

PAVX tests now passing

parent d9dca7fc
No related branches found
No related tags found
No related merge requests found
......@@ -44,6 +44,7 @@ setup(
[
"xbob/math/histogram.cpp",
"xbob/math/linsolve.cpp",
"xbob/math/pavx.cpp",
"xbob/math/main.cpp",
],
packages = packages,
......
......@@ -43,6 +43,8 @@ static PyObject* py_histogram_intersection_1
PyObject* retval = 0;
try {
switch(h1->type_num) {
case NPY_UINT8:
......@@ -70,6 +72,14 @@ static PyObject* py_histogram_intersection_1
}
}
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);
return retval;
......@@ -178,6 +188,8 @@ static PyObject* py_histogram_intersection_2(PyObject*, PyObject* args, PyObject
PyObject* retval = 0;
try {
switch(index1->type_num) {
case NPY_UINT8:
......@@ -210,6 +222,14 @@ static PyObject* py_histogram_intersection_2(PyObject*, PyObject* args, PyObject
}
}
catch (std::exception& e) {
PyErr_SetString(PyExc_RuntimeError, e.what());
}
catch (...) {
PyErr_SetString(PyExc_RuntimeError, "histogram intersection failed: unknown exception caught");
}
Py_DECREF(index1);
Py_DECREF(value1);
Py_DECREF(index2);
......
......@@ -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>
......
......@@ -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,
......@@ -84,10 +85,12 @@ linsolve_(A, x, b) -> None\n\
\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\
.. 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 are well-behaved (contiguous, c-style, memory aligned).\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\
......@@ -120,10 +123,12 @@ linsolve_sympos_(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 ``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 are well-behaved (contiguous, c-style, memory aligned).\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\
......@@ -156,10 +161,12 @@ linsolve_cg_sympos_(A, x, b) -> None\n\
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\
.. 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 are well-behaved (contiguous, c-style, memory aligned).\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 */
};
......
......@@ -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) {
/* 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;
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.";
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 && 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;
}
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_);
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");
}
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_);
Py_DECREF(input);
if (returns_output) return (PyObject*)output;
Py_DECREF(output);
Py_RETURN_NONE;
}
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();
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;
}
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();
// 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;
}
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);
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");
}
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);
Py_DECREF(input);
Py_DECREF(output);
Py_RETURN_NONE;
}
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;
}
/**
* @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);
......@@ -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>
......
......@@ -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)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment