num_util.cpp
Go to the documentation of this file.
1 // Copyright 2006 Phil Austin (http://www.eos.ubc.ca/personal/paustin)
2 // Distributed under the Boost Software License, Version 1.0. (See
3 // accompanying file LICENSE_1_0.txt or copy at
4 // http://www.boost.org/LICENSE_1_0.txt)
5 
6 #define PY_ARRAY_UNIQUE_SYMBOL HippoPyArrayHandle
7 #define NO_IMPORT_ARRAY
8 #include "num_util.h"
9 
10 namespace { const char* rcsid = "$Id: num_util.cpp,v 1.12 2007/05/31 21:17:55 pfkeb Exp $"; }
11 
12 using namespace boost::python;
13 
14 namespace num_util{
15 
16  //specializations for use by makeNum
17 
18 
19  template <>
20  PyArray_TYPES getEnum<unsigned char>(void)
21  {
22  return PyArray_UBYTE;
23  }
24 
25 
26  template <>
27  PyArray_TYPES getEnum<signed char>(void)
28  {
29 #ifdef HAVE_NUMPY
30  return PyArray_BYTE;
31 #else
32  return PyArray_SBYTE;
33 #endif
34  }
35 
36  template <>
37  PyArray_TYPES getEnum<short>(void)
38  {
39  return PyArray_SHORT;
40  }
41 
42  template <>
43  PyArray_TYPES getEnum<unsigned short>(void)
44  {
45  return PyArray_USHORT;
46  }
47 
48 
49  template <>
50  PyArray_TYPES getEnum<unsigned int>(void)
51  {
52  return PyArray_UINT;
53  }
54 
55  template <>
56  PyArray_TYPES getEnum<int>(void)
57  {
58  return PyArray_INT;
59  }
60 
61  template <>
62  PyArray_TYPES getEnum<long>(void)
63  {
64  return PyArray_LONG;
65  }
66 
67 #ifdef HAVE_NUMPY
68  template <>
69  PyArray_TYPES getEnum<unsigned long>(void)
70  {
71  return PyArray_ULONG;
72  }
73 
74  template <>
75  PyArray_TYPES getEnum<long long>(void)
76  {
77  return PyArray_LONGLONG;
78  }
79 
80  template <>
81  PyArray_TYPES getEnum<unsigned long long>(void)
82  {
83  return PyArray_ULONGLONG;
84  }
85 #endif
86 
87  template <>
88  PyArray_TYPES getEnum<float>(void)
89  {
90  return PyArray_FLOAT;
91  }
92 
93  template <>
94  PyArray_TYPES getEnum<double>(void)
95  {
96  return PyArray_DOUBLE;
97  }
98 
99 #ifdef HAVE_NUMPY
100  template <>
101  PyArray_TYPES getEnum<long double>(void)
102  {
103  return PyArray_LONGDOUBLE;
104  }
105 #endif
106  template <>
107  PyArray_TYPES getEnum<std::complex<float> >(void)
108  {
109  return PyArray_CFLOAT;
110  }
111 
112 
113  template <>
114  PyArray_TYPES getEnum<std::complex<double> >(void)
115  {
116  return PyArray_CDOUBLE;
117  }
118 
119 #ifdef HAVE_NUMPY
120  template <>
121  PyArray_TYPES getEnum<std::complex<long double> >(void)
122  {
123  return PyArray_CLONGDOUBLE;
124  }
125 #endif
126 
127  template <> boost::python::numeric::array makeNum<double> (double * data, std::vector<int> dims){
128  intp total = std::accumulate(dims.begin(),dims.end(),1,std::multiplies<intp>());
129  boost::python::object obj(boost::python::handle<>(PyArray_FromDims(dims.size(),&dims[0], PyArray_DOUBLE)));
130 #ifdef HAVE_NUMPY
131  void *arr_data = PyArray_DATA((PyArrayObject*) obj.ptr());
132  memcpy(arr_data, data, PyArray_ITEMSIZE((PyArrayObject*) obj.ptr()) * total);
133 #else
134  char *arr_data = ((PyArrayObject*) obj.ptr())->data;
135  memcpy(arr_data, data, sizeof(double) * total); // copies the input data to
136  // PyArrayObject->data
137 #endif
138  return boost::python::extract<boost::python::numeric::array>(obj);
139  }
140 
141 
142 typedef KindStringMap::value_type KindStringMapEntry;
144  {
145  KindStringMapEntry(PyArray_UBYTE, "PyArray_UBYTE"),
146 #ifdef HAVE_NUMPY
147  KindStringMapEntry(PyArray_BYTE, "PyArray_BYTE"),
148 #else
149  KindStringMapEntry(PyArray_SBYTE, "PyArray_SBYTE"),
150 #endif
151  KindStringMapEntry(PyArray_SHORT, "PyArray_SHORT"),
152  KindStringMapEntry(PyArray_INT, "PyArray_INT"),
153  KindStringMapEntry(PyArray_LONG, "PyArray_LONG"),
154  KindStringMapEntry(PyArray_FLOAT, "PyArray_FLOAT"),
155  KindStringMapEntry(PyArray_DOUBLE, "PyArray_DOUBLE"),
156  KindStringMapEntry(PyArray_CFLOAT, "PyArray_CFLOAT"),
157  KindStringMapEntry(PyArray_CDOUBLE,"PyArray_CDOUBLE"),
158  KindStringMapEntry(PyArray_OBJECT, "PyArray_OBJECT"),
159  KindStringMapEntry(PyArray_NTYPES, "PyArray_NTYPES"),
160  KindStringMapEntry(PyArray_NOTYPE ,"PyArray_NOTYPE")
161  };
162 
163 typedef KindCharMap::value_type KindCharMapEntry;
165  {
166  KindCharMapEntry(PyArray_UBYTE, 'B'),
167 #ifdef HAVE_NUMPY
168  KindCharMapEntry(PyArray_BYTE, 'b'),
169 #else
170  KindCharMapEntry(PyArray_SBYTE, 'b'),
171 #endif
172  KindCharMapEntry(PyArray_SHORT, 'h'),
173  KindCharMapEntry(PyArray_INT, 'i'),
174  KindCharMapEntry(PyArray_LONG, 'l'),
175  KindCharMapEntry(PyArray_FLOAT, 'f'),
176  KindCharMapEntry(PyArray_DOUBLE, 'd'),
177  KindCharMapEntry(PyArray_CFLOAT, 'F'),
178  KindCharMapEntry(PyArray_CDOUBLE,'D'),
179  KindCharMapEntry(PyArray_OBJECT, 'O')
180  };
181 
182 typedef KindTypeMap::value_type KindTypeMapEntry;
184  {
185  KindTypeMapEntry('B',PyArray_UBYTE),
186 #ifdef HAVE_NUMPY
187  KindTypeMapEntry('b',PyArray_BYTE),
188 #else
189  KindTypeMapEntry('b',PyArray_SBYTE),
190 #endif
191  KindTypeMapEntry('h',PyArray_SHORT),
192  KindTypeMapEntry('i',PyArray_INT),
193  KindTypeMapEntry('l',PyArray_LONG),
194  KindTypeMapEntry('f',PyArray_FLOAT),
195  KindTypeMapEntry('d',PyArray_DOUBLE),
196  KindTypeMapEntry('F',PyArray_CFLOAT),
197  KindTypeMapEntry('D',PyArray_CDOUBLE),
198  KindTypeMapEntry('O',PyArray_OBJECT)
199  };
200 
201 
205 
206 
207 using namespace boost::python;
208 
211 
214 
217 
218 //Create a numarray referencing Python sequence object
219 numeric::array makeNum(object x){
220  if (!PySequence_Check(x.ptr())){
221  PyErr_SetString(PyExc_ValueError, "expected a sequence");
222  throw_error_already_set();
223  }
224  object obj(handle<>
225  (PyArray_ContiguousFromObject(x.ptr(),PyArray_NOTYPE,0,0)));
227  return extract<numeric::array>(obj);
228 }
229 
230 //Create a one-dimensional Numeric array of length n and Numeric type t
231 numeric::array makeNum(int n, PyArray_TYPES t=PyArray_DOUBLE){
232  object obj(handle<>(PyArray_FromDims(1, &n, t)));
233  return extract<numeric::array>(obj);
234 }
235 
236 //Create a Numeric array with dimensions dimens and Numeric type t
237 numeric::array makeNum(std::vector<int> dimens,
238  PyArray_TYPES t=PyArray_DOUBLE){
239  object obj(handle<>(PyArray_FromDims(dimens.size(), &dimens[0], t)));
240  return extract<numeric::array>(obj);
241 }
242 
243 numeric::array makeNum(const numeric::array& arr){
244  //Returns a reference of arr by calling numeric::array copy constructor.
245  //The copy constructor increases arr's reference count.
246  return numeric::array(arr);
247 }
248 
249 PyArray_TYPES type(numeric::array arr){
250 #ifdef HAVE_NUMPY
251  return PyArray_TYPES(PyArray_TYPE(arr.ptr()));
252 #else
253  return char2type(arr.typecode());
254 #endif
255 }
256 
257 void check_type(boost::python::numeric::array arr,
258  PyArray_TYPES expected_type){
259  PyArray_TYPES actual_type = type(arr);
260  if (actual_type != expected_type) {
261  std::ostringstream stream;
262  stream << "expected Numeric type " << kindstrings[expected_type]
263  << ", found Numeric type " << kindstrings[actual_type] << std::ends;
264  PyErr_SetString(PyExc_TypeError, stream.str().c_str());
265  throw_error_already_set();
266  }
267  return;
268 }
269 
270 //Return the number of dimensions
271 int rank(numeric::array arr){
272  //std::cout << "inside rank" << std::endl;
273  if(!PyArray_Check(arr.ptr())){
274  PyErr_SetString(PyExc_ValueError, "expected a PyArrayObject");
275  throw_error_already_set();
276  }
277 #ifdef HAVE_NUMPY
278  return PyArray_NDIM(arr.ptr());
279 #else
280  return ((PyArrayObject*) arr.ptr())->nd;
281 #endif
282 }
283 
284 void check_rank(boost::python::numeric::array arr, int expected_rank){
285  int actual_rank = rank(arr);
286  if (actual_rank != expected_rank) {
287  std::ostringstream stream;
288  stream << "expected rank " << expected_rank
289  << ", found rank " << actual_rank << std::ends;
290  PyErr_SetString(PyExc_RuntimeError, stream.str().c_str());
291  throw_error_already_set();
292  }
293  return;
294 }
295 
296 intp size(numeric::array arr)
297 {
298  if(!PyArray_Check(arr.ptr())){
299  PyErr_SetString(PyExc_ValueError, "expected a PyArrayObject");
300  throw_error_already_set();
301  }
302  return PyArray_Size(arr.ptr());
303 }
304 
305 void check_size(boost::python::numeric::array arr, intp expected_size){
306  intp actual_size = size(arr);
307  if (actual_size != expected_size) {
308  std::ostringstream stream;
309  stream << "expected size " << expected_size
310  << ", found size " << actual_size << std::ends;
311  PyErr_SetString(PyExc_RuntimeError, stream.str().c_str());
312  throw_error_already_set();
313  }
314  return;
315 }
316 
317 std::vector<intptr_t> shape(numeric::array arr){
318  std::vector<intptr_t> out_dims;
319  if(!PyArray_Check(arr.ptr())){
320  PyErr_SetString(PyExc_ValueError, "expected a PyArrayObject");
321  throw_error_already_set();
322  }
323 #ifdef HAVE_NUMPY
324  intp* dims_ptr = PyArray_DIMS(arr.ptr());
325 #else
326  int* dims_ptr = ((PyArrayObject*) arr.ptr())->dimensions;
327 #endif
328  int the_rank = rank(arr);
329  for (int i = 0; i < the_rank; i++){
330  out_dims.push_back(*(dims_ptr + i));
331  }
332  return out_dims;
333 }
334 
335 intp get_dim(boost::python::numeric::array arr, int dimnum){
336  assert(dimnum >= 0);
337  int the_rank=rank(arr);
338  if(the_rank < dimnum){
339  std::ostringstream stream;
340  stream << "Error: asked for length of dimension ";
341  stream << dimnum << " but rank of array is " << the_rank << std::ends;
342  PyErr_SetString(PyExc_RuntimeError, stream.str().c_str());
343  throw_error_already_set();
344  }
345  std::vector<intptr_t> actual_dims = shape(arr);
346  return actual_dims[dimnum];
347 }
348 
349 void check_shape(boost::python::numeric::array arr, std::vector<intptr_t> expected_dims){
350  std::vector<intptr_t> actual_dims = shape(arr);
351  if (actual_dims != expected_dims) {
352  std::ostringstream stream;
353  stream << "expected dimensions " << vector_str(expected_dims)
354  << ", found dimensions " << vector_str(actual_dims) << std::ends;
355  PyErr_SetString(PyExc_RuntimeError, stream.str().c_str());
356  throw_error_already_set();
357  }
358  return;
359 }
360 
361 void check_dim(boost::python::numeric::array arr, int dimnum, intp dimsize){
362  std::vector<intptr_t> actual_dims = shape(arr);
363  if(actual_dims[dimnum] != dimsize){
364  std::ostringstream stream;
365  stream << "Error: expected dimension number ";
366  stream << dimnum << " to be length " << dimsize;
367  stream << ", but found length " << actual_dims[dimnum] << std::ends;
368  PyErr_SetString(PyExc_RuntimeError, stream.str().c_str());
369  throw_error_already_set();
370  }
371  return;
372 }
373 
374 bool iscontiguous(numeric::array arr)
375 {
376  // return arr.iscontiguous();
377  return PyArray_ISCONTIGUOUS(arr.ptr());
378 }
379 
380 void check_contiguous(numeric::array arr)
381 {
382  if (!iscontiguous(arr)) {
383  PyErr_SetString(PyExc_RuntimeError, "expected a contiguous array");
384  throw_error_already_set();
385  }
386  return;
387 }
388 
389 void* data(numeric::array arr){
390  if(!PyArray_Check(arr.ptr())){
391  PyErr_SetString(PyExc_ValueError, "expected a PyArrayObject");
392  throw_error_already_set();
393  }
394 #ifdef HAVE_NUMPY
395  return PyArray_DATA(arr.ptr());
396 #else
397  return ((PyArrayObject*) arr.ptr())->data;
398 #endif
399 }
400 
401 //Copy data into the array
402 void copy_data(boost::python::numeric::array arr, char* new_data){
403  char* arr_data = (char*) data(arr);
404  intp nbytes = PyArray_NBYTES(arr.ptr());
405  for (intp i = 0; i < nbytes; i++) {
406  arr_data[i] = new_data[i];
407  }
408  return;
409 }
410 
411 //Return a clone of this array
412 numeric::array clone(numeric::array arr){
413 #ifdef HAVE_NUMPY
414  object obj(handle<>(PyArray_NewCopy((PyArrayObject*)arr.ptr(),PyArray_CORDER)));
415 #else
416  object obj(handle<>(PyArray_Copy((PyArrayObject*)arr.ptr())));
417 #endif
418  return makeNum(obj);
419 }
420 
421 
422 //Return a clone of this array with a new type
423 numeric::array astype(boost::python::numeric::array arr, PyArray_TYPES t){
424  return (numeric::array) arr.astype(type2char(t));
425 }
426 
427 std::vector<intp> strides(numeric::array arr){
428  std::vector<intp> out_strides;
429  if(!PyArray_Check(arr.ptr())){
430  PyErr_SetString(PyExc_ValueError, "expected a PyArrayObject");
431  throw_error_already_set();
432  }
433 #ifdef HAVE_NUMPY
434  intp* strides_ptr = PyArray_STRIDES(arr.ptr());
435 #else
436  int* strides_ptr = ((PyArrayObject*) arr.ptr())->strides;
437 #endif
438  intp the_rank = rank(arr);
439  for (intp i = 0; i < the_rank; i++){
440  out_strides.push_back(*(strides_ptr + i));
441  }
442  return out_strides;
443 }
444 
445 int refcount(numeric::array arr){
446  return REFCOUNT(arr.ptr());
447 }
448 
449 void check_PyArrayElementType(object newo){
450 #ifdef HAVE_NUMPY
451  PyArray_TYPES theType=PyArray_TYPES(PyArray_TYPE(newo.ptr()));
452 #else
453  PyArrayObject* PyArrayPtr = (PyArrayObject*) newo.ptr();
454  PyArray_TYPES theType=(PyArray_TYPES) PyArrayPtr->descr->type_num;
455 #endif
456  if(theType == PyArray_OBJECT){
457  std::ostringstream stream;
458  stream << "array elments have been cast to PyArray_OBJECT, "
459  << "numhandle can only accept arrays with numerical elements"
460  << std::ends;
461  PyErr_SetString(PyExc_TypeError, stream.str().c_str());
462  throw_error_already_set();
463  }
464  return;
465 }
466 
467 std::string type2string(PyArray_TYPES t_type){
468  return kindstrings[t_type];
469 }
470 
471 char type2char(PyArray_TYPES t_type){
472  return kindchars[t_type];
473 }
474 
475 PyArray_TYPES char2type(char e_type){
476  return kindtypes[e_type];
477 }
478 
479 template <class T>
480 inline std::string vector_str(const std::vector<T>& vec)
481 {
482  std::ostringstream stream;
483  stream << "(" << vec[0];
484 
485  for(std::size_t i = 1; i < vec.size(); i++){
486  stream << ", " << vec[i];
487  }
488  stream << ")";
489  return stream.str();
490 }
491 
492 inline void check_size_match(std::vector<intp> dims, intp n)
493 {
494  intp total = std::accumulate(dims.begin(),dims.end(),1,std::multiplies<intp>());
495  if (total != n) {
496  std::ostringstream stream;
497  stream << "expected array size " << n
498  << ", dimensions give array size " << total << std::ends;
499  PyErr_SetString(PyExc_TypeError, stream.str().c_str());
500  throw_error_already_set();
501  }
502  return;
503 }
504 
505 } //namespace num_util
506 

Generated for HippoDraw Class Library by doxygen