summaryrefslogtreecommitdiffstats
path: root/src/boost/libs/python/example/numpy/wrap.cpp
blob: 33f71c85ade8149b7c178cbb26a5becb650ced54 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
// Copyright Jim Bosch 2011-2012.
// Distributed under the Boost Software License, Version 1.0.
//    (See accompanying file LICENSE_1_0.txt or copy at
//          http://www.boost.org/LICENSE_1_0.txt)

/**
 *  A simple example showing how to wrap a couple of C++ functions that
 *  operate on 2-d arrays into Python functions that take NumPy arrays
 *  as arguments.
 *        
 *  If you find have a lot of such functions to wrap, you may want to
 *  create a C++ array type (or use one of the many existing C++ array
 *  libraries) that maps well to NumPy arrays and create Boost.Python
 *  converters.  There's more work up front than the approach here,
 *  but much less boilerplate per function.  See the "Gaussian" example
 *  included with Boost.NumPy for an example of custom converters, or
 *  take a look at the "ndarray" project on GitHub for a more complete,
 *  high-level solution.
 *
 *  Note that we're using embedded Python here only to make a convenient
 *  self-contained example; you could just as easily put the wrappers
 *  in a regular C++-compiled module and imported them in regular
 *  Python.  Again, see the Gaussian demo for an example.
 */

#include <boost/python/numpy.hpp>
#include <boost/scoped_array.hpp>
#include <iostream>

namespace p = boost::python;
namespace np = boost::python::numpy;

// This is roughly the most efficient way to write a C/C++ function that operates
// on a 2-d NumPy array - operate directly on the array by incrementing a pointer
// with the strides.
void fill1(double * array, int rows, int cols, int row_stride, int col_stride) {
    double * row_iter = array;
    double n = 0.0; // just a counter we'll fill the array with.
    for (int i = 0; i < rows; ++i, row_iter += row_stride) {
        double * col_iter = row_iter;
        for (int j = 0; j < cols; ++j, col_iter += col_stride) {
            *col_iter = ++n;
        }
    }
}

// Here's a simple wrapper function for fill1.  It requires that the passed
// NumPy array be exactly what we're looking for - no conversion from nested
// sequences or arrays with other data types, because we want to modify it
// in-place.
void wrap_fill1(np::ndarray const & array) {
    if (array.get_dtype() != np::dtype::get_builtin<double>()) {
        PyErr_SetString(PyExc_TypeError, "Incorrect array data type");
        p::throw_error_already_set();
    }
    if (array.get_nd() != 2) {
        PyErr_SetString(PyExc_TypeError, "Incorrect number of dimensions");
        p::throw_error_already_set();
    }
    fill1(reinterpret_cast<double*>(array.get_data()),
          array.shape(0), array.shape(1),
          array.strides(0) / sizeof(double), array.strides(1) / sizeof(double));
}

// Another fill function that takes a double**.  This is less efficient, because
// it's not the native NumPy data layout, but it's common enough in C/C++ that
// it's worth its own example.  This time we don't pass the strides, and instead
// in wrap_fill2 we'll require the C_CONTIGUOUS bitflag, which guarantees that
// the column stride is 1 and the row stride is the number of columns.  That
// restricts the arrays that can be passed to fill2 (it won't work on most
// subarray views or transposes, for instance).
void fill2(double ** array, int rows, int cols) {
    double n = 0.0; // just a counter we'll fill the array with.
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            array[i][j] = ++n;
        }
    }    
}
// Here's the wrapper for fill2; it's a little more complicated because we need
// to check the flags and create the array of pointers.
void wrap_fill2(np::ndarray const & array) {
    if (array.get_dtype() != np::dtype::get_builtin<double>()) {
        PyErr_SetString(PyExc_TypeError, "Incorrect array data type");
        p::throw_error_already_set();
    }
    if (array.get_nd() != 2) {
        PyErr_SetString(PyExc_TypeError, "Incorrect number of dimensions");
        p::throw_error_already_set();
    }
    if (!(array.get_flags() & np::ndarray::C_CONTIGUOUS)) {
        PyErr_SetString(PyExc_TypeError, "Array must be row-major contiguous");
        p::throw_error_already_set();
    }
    double * iter = reinterpret_cast<double*>(array.get_data());
    int rows = array.shape(0);
    int cols = array.shape(1);
    boost::scoped_array<double*> ptrs(new double*[rows]);
    for (int i = 0; i < rows; ++i, iter += cols) {
        ptrs[i] = iter;
    }
    fill2(ptrs.get(), array.shape(0), array.shape(1));
}

BOOST_PYTHON_MODULE(example) {
    np::initialize();  // have to put this in any module that uses Boost.NumPy
    p::def("fill1", wrap_fill1);
    p::def("fill2", wrap_fill2);
}

int main(int argc, char **argv)
{
    // This line makes our module available to the embedded Python intepreter.
# if PY_VERSION_HEX >= 0x03000000
    PyImport_AppendInittab("example", &PyInit_example);
# else
    PyImport_AppendInittab("example", &initexample);
# endif
    // Initialize the Python runtime.
    Py_Initialize();

    PyRun_SimpleString(
        "import example\n"
        "import numpy\n"
        "z1 = numpy.zeros((5,6), dtype=float)\n"
        "z2 = numpy.zeros((4,3), dtype=float)\n"
        "example.fill1(z1)\n"
        "example.fill2(z2)\n"
        "print z1\n"
        "print z2\n"
    );
    Py_Finalize();
}