[Pyublas] Issues with 3d arrays

Anton Loukianov

2012-04-23 16:56:19 UTC

Hello,

I am trying to write a Euclidean distance function that takes two
numpy arrays. In trying to vectorize it from 2 to 3 dimensional
arrays, I ran into a bug that I can't seem to figure out. I would
appreciate any help.

----- C++ code -----
#include <pyublas/numpy.hpp>
#include <stdexcept>
#include <cmath>

typedef pyublas::numpy_vector<double> py_vec;

inline double square(double a) { return a*a; }

inline double dist(double ax, double ay, double bx, double by)
{
return std::sqrt(square(ax - bx) + square(ay - by));
}

py_vec distance2d(py_vec a, py_vec b)
{
if(a.size() != b.size()) throw std::invalid_argument("mismatched
argument lengths");
if(a.ndim() < 2 || b.ndim() < 2) throw std::invalid_argument("need
at least 2d arrays");
if(a.ndim() > 3 || b.ndim() > 3) throw
std::invalid_argument("can't handle greater than 3d arrays");

py_vec z(a.size()/2);

const npy_intp* dims = a.dims();

if(a.ndim() == 2) {
for (int i = 0; i < dims[0]; i++) {
double ax = a.sub(i,0);
double ay = a.sub(i,1);
double bx = b.sub(i,0);
double by = b.sub(i,1);

z[i] = dist(ax, ay, bx, by);
}
} else {
for (int i = 0; i < dims[0]; i++) {
for (int j = 0; j < dims[1]; j++) {
double ax = a.sub(i,j,0);
double ay = a.sub(i,j,1);
double bx = b.sub(i,j,0);
double by = b.sub(i,j,1);

z[i] = dist(ax, ay, bx, by);
}
}
}
return z;
}
BOOST_PYTHON_MODULE(_osc_helper)
{
boost::python::def("distance2d", distance2d);
}
----- End C++ code -----

The idea is that it should behave like this from python:

----- Python code -----
2d array case (this seems to work)

import pyublas; import _osc_helper; import numpy as np;
a = tile([0, 1], (5, 1)).astype(float)
_osc_helper.distance2d(a, zeros((5, 2), dtype=float))

array([ 1., 1., 1., 1., 1.])

import pyublas; import _osc_helper; import numpy as np;
a = tile([0, 1], (5, 2, 1)).astype(float)
_osc_helper.distance2d(a, zeros((5, 2, 2), dtype=float))

array([ 1.00000000e+000, 1.00000000e+000, 1.00000000e+000,
1.00000000e+000, 1.00000000e+000, 0.00000000e+000,
0.00000000e+000, 0.00000000e+000, -5.70857814e+168,
3.02147214e-306])
----- End python code -----

Am I doing something horribly wrong?

Thank you in advance,
Anton Loukianov