Skip to content

Commit 019d53d

Browse files
dblalockdblalock
committed
succesfully wrap eigen matrices + vectors
1 parent f8a4d88 commit 019d53d

File tree

329 files changed

+91252
-33
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

329 files changed

+91252
-33
lines changed

‎eigen.i‎

Lines changed: 296 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,296 @@
1+
/*
2+
* The Biomechanical ToolKit
3+
* Copyright (c) 2009-2013, Arnaud Barré
4+
* All rights reserved.
5+
*
6+
* Redistribution and use in source and binary forms, with or without
7+
* modification, are permitted provided that the following conditions
8+
* are met:
9+
*
10+
* * Redistributions of source code must retain the above
11+
* copyright notice, this list of conditions and the following
12+
* disclaimer.
13+
* * Redistributions in binary form must reproduce the above
14+
* copyright notice, this list of conditions and the following
15+
* disclaimer in the documentation and/or other materials
16+
* provided with the distribution.
17+
* * Neither the name(s) of the copyright holders nor the names
18+
* of its contributors may be used to endorse or promote products
19+
* derived from this software without specific prior written
20+
* permission.
21+
*
22+
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
23+
* "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
24+
* LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
25+
* FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
26+
* COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
27+
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
28+
* BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
29+
* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
30+
* CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31+
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
32+
* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
33+
* POSSIBILITY OF SUCH DAMAGE.
34+
*/
35+
36+
%{
37+
#define SWIG_FILE_WITH_INIT
38+
#include "Eigen/Core"
39+
%}
40+
41+
%include "numpy.i"
42+
43+
%init
44+
%{
45+
import_array();
46+
%}
47+
48+
%fragment("Eigen_Fragments", "header", fragment="NumPy_Fragments")
49+
%{
50+
template <typename T> int NumPyType() {return -1;};
51+
52+
template <class Derived>
53+
void ConvertFromNumpyToEigenMatrix(Eigen::MatrixBase<Derived>* out, PyObject* in)
54+
{
55+
int rows = 0;
56+
int cols = 0;
57+
// Check object type
58+
if (!is_array(in))
59+
{
60+
PyErr_SetString(PyExc_ValueError, "The given input is not known as a NumPy array or matrix.");
61+
return;
62+
}
63+
// Check data type
64+
else if (array_type(in) != NumPyType<typename Derived::Scalar>())
65+
{
66+
PyErr_SetString(PyExc_ValueError, "Type mismatch between NumPy and Eigen objects.");
67+
return;
68+
}
69+
// Check dimensions
70+
else if (array_numdims(in) > 2)
71+
{
72+
PyErr_SetString(PyExc_ValueError, "Eigen only support 1D or 2D array.");
73+
return;
74+
}
75+
else if (array_numdims(in) == 1)
76+
{
77+
rows = array_size(in,0);
78+
cols = 1;
79+
if ((Derived::RowsAtCompileTime != Eigen::Dynamic) && (Derived::RowsAtCompileTime != rows))
80+
{
81+
PyErr_SetString(PyExc_ValueError, "Row dimension mismatch between NumPy and Eigen objects (1D).");
82+
return;
83+
}
84+
else if ((Derived::ColsAtCompileTime != Eigen::Dynamic) && (Derived::ColsAtCompileTime != 1))
85+
{
86+
PyErr_SetString(PyExc_ValueError, "Column dimension mismatch between NumPy and Eigen objects (1D).");
87+
return;
88+
}
89+
}
90+
else if (array_numdims(in) == 2)
91+
{
92+
rows = array_size(in,0);
93+
cols = array_size(in,1);
94+
if ((Derived::RowsAtCompileTime != Eigen::Dynamic) && (Derived::RowsAtCompileTime != array_size(in,0)))
95+
{
96+
PyErr_SetString(PyExc_ValueError, "Row dimension mismatch between NumPy and Eigen objects (2D).");
97+
return;
98+
}
99+
else if ((Derived::ColsAtCompileTime != Eigen::Dynamic) && (Derived::ColsAtCompileTime != array_size(in,1)))
100+
{
101+
PyErr_SetString(PyExc_ValueError, "Column dimension mismatch between NumPy and Eigen objects (2D).");
102+
return;
103+
}
104+
}
105+
// Extract data
106+
int isNewObject = 0;
107+
PyArrayObject* temp = obj_to_array_contiguous_allow_conversion(in, array_type(in), &isNewObject);
108+
if (temp == NULL)
109+
{
110+
PyErr_SetString(PyExc_ValueError, "Impossible to convert the input into a Python array object.");
111+
return;
112+
}
113+
out->derived().setZero(rows, cols);
114+
typename Derived::Scalar* data = static_cast<typename Derived::Scalar*>(array_data(temp));
115+
for (int i = 0; i != rows; ++i)
116+
for (int j = 0; j != cols; ++j)
117+
out->coeffRef(i,j) = data[i*cols+j];
118+
};
119+
120+
// Copies values from Eigen type into an existing NumPy type
121+
template <class Derived>
122+
void CopyFromEigenToNumPyMatrix(PyObject* out, Eigen::MatrixBase<Derived>* in)
123+
{
124+
int rows = 0;
125+
int cols = 0;
126+
// Check object type
127+
if (!is_array(out))
128+
{
129+
PyErr_SetString(PyExc_ValueError, "The given input is not known as a NumPy array or matrix.");
130+
return;
131+
}
132+
// Check data type
133+
else if (array_type(out) != NumPyType<typename Derived::Scalar>())
134+
{
135+
PyErr_SetString(PyExc_ValueError, "Type mismatch between NumPy and Eigen objects.");
136+
return;
137+
}
138+
// Check dimensions
139+
else if (array_numdims(out) > 2)
140+
{
141+
PyErr_SetString(PyExc_ValueError, "Eigen only support 1D or 2D array.");
142+
return;
143+
}
144+
else if (array_numdims(out) == 1)
145+
{
146+
rows = array_size(out,0);
147+
cols = 1;
148+
if ((Derived::RowsAtCompileTime != Eigen::Dynamic) && (Derived::RowsAtCompileTime != rows))
149+
{
150+
PyErr_SetString(PyExc_ValueError, "Row dimension mismatch between NumPy and Eigen objects (1D).");
151+
return;
152+
}
153+
else if ((Derived::ColsAtCompileTime != Eigen::Dynamic) && (Derived::ColsAtCompileTime != 1))
154+
{
155+
PyErr_SetString(PyExc_ValueError, "Column dimension mismatch between NumPy and Eigen objects (1D).");
156+
return;
157+
}
158+
}
159+
else if (array_numdims(out) == 2)
160+
{
161+
rows = array_size(out,0);
162+
cols = array_size(out,1);
163+
}
164+
165+
if (in->cols() != cols || in->rows() != rows) {
166+
/// TODO: be forgiving and simply create or resize the array
167+
PyErr_SetString(PyExc_ValueError, "Dimension mismatch between NumPy and Eigen object (return argument).");
168+
return;
169+
}
170+
171+
// Extract data
172+
int isNewObject = 0;
173+
PyArrayObject* temp = obj_to_array_contiguous_allow_conversion(out, array_type(out), &isNewObject);
174+
if (temp == NULL)
175+
{
176+
PyErr_SetString(PyExc_ValueError, "Impossible to convert the input into a Python array object.");
177+
return;
178+
}
179+
180+
typename Derived::Scalar* data = static_cast<typename Derived::Scalar*>(array_data(out));
181+
182+
for (int i = 0; i != in->rows(); ++i) {
183+
for (int j = 0; j != in->cols(); ++j) {
184+
data[i*in->cols()+j] = in->coeff(i,j);
185+
}
186+
}
187+
};
188+
189+
template <class Derived>
190+
void ConvertFromEigenToNumPyMatrix(PyObject** out, Eigen::MatrixBase<Derived>* in)
191+
{
192+
// vector (1D)
193+
if (in->cols() == 1) {
194+
npy_intp dims[1] = {in->rows()};
195+
*out = PyArray_SimpleNew(1, dims, NumPyType<typename Derived::Scalar>());
196+
typename Derived::Scalar* data = static_cast<typename Derived::Scalar*>(array_data(*out));
197+
for (int i = 0; i != dims[0]; ++i)
198+
data[i] = in->coeff(i, 1);
199+
return;
200+
}
201+
// matrix (2D)
202+
npy_intp dims[2] = {in->rows(), in->cols()};
203+
*out = PyArray_SimpleNew(2, dims, NumPyType<typename Derived::Scalar>());
204+
typename Derived::Scalar* data = static_cast<typename Derived::Scalar*>(array_data(*out));
205+
for (int i = 0; i != dims[0]; ++i)
206+
for (int j = 0; j != dims[1]; ++j)
207+
data[i*dims[1]+j] = in->coeff(i,j);
208+
};
209+
210+
// these funcs define the mapping between c types and numpy types;
211+
// add more as needed
212+
template<> int NumPyType<double>() {return NPY_DOUBLE;};
213+
template<> int NumPyType<float>() {return NPY_FLOAT;};
214+
template<> int NumPyType<int>() {return NPY_INT;};
215+
template<> int NumPyType<long>() {return NPY_LONG;};
216+
%}
217+
218+
// ----------------------------------------------------------------------------
219+
// Macro to create the typemap for Eigen classes
220+
// ----------------------------------------------------------------------------
221+
%define %eigen_typemaps(CLASS)
222+
223+
// Argout: const & (Disabled and prevents calling of the non-const typemap)
224+
%typemap(argout, fragment="Eigen_Fragments") const CLASS & ""
225+
226+
// Argout: & (for returning values to in-out arguments)
227+
%typemap(argout, fragment="Eigen_Fragments") CLASS &
228+
{
229+
// Argout: &
230+
CopyFromEigenToNumPyMatrix<CLASS>($input, $1);
231+
}
232+
233+
// In: (nothing: no constness)
234+
%typemap(in, fragment="Eigen_Fragments") CLASS (CLASS temp)
235+
{
236+
ConvertFromNumpyToEigenMatrix<CLASS>(&temp, $input);
237+
$1 = temp;
238+
}
239+
// In: const&
240+
%typemap(in, fragment="Eigen_Fragments") CLASS const& (CLASS temp)
241+
{
242+
// In: const&
243+
ConvertFromNumpyToEigenMatrix<CLASS>(&temp, $input);
244+
$1 = &temp;
245+
}
246+
// In: & (not yet implemented)
247+
%typemap(in, fragment="Eigen_Fragments") CLASS & (CLASS temp)
248+
{
249+
// In: non-const&
250+
ConvertFromNumpyToEigenMatrix<CLASS>(&temp, $input);
251+
252+
$1 = &temp;
253+
}
254+
// In: const* (not yet implemented)
255+
%typemap(in, fragment="Eigen_Fragments") CLASS const*
256+
{
257+
PyErr_SetString(PyExc_ValueError, "The input typemap for const pointer is not yet implemented. Please report this problem to the developer.");
258+
}
259+
// In: * (not yet implemented)
260+
%typemap(in, fragment="Eigen_Fragments") CLASS *
261+
{
262+
PyErr_SetString(PyExc_ValueError, "The input typemap for non-const pointer is not yet implemented. Please report this problem to the developer.");
263+
}
264+
265+
// Out: (nothing: no constness)
266+
%typemap(out, fragment="Eigen_Fragments") CLASS
267+
{
268+
ConvertFromEigenToNumPyMatrix<CLASS>(&$result, &$1);
269+
}
270+
// Out: const
271+
%typemap(out, fragment="Eigen_Fragments") CLASS const
272+
{
273+
ConvertFromEigenToNumPyMatrix<CLASS>(&$result, &$1);
274+
}
275+
// Out: const&
276+
%typemap(out, fragment="Eigen_Fragments") CLASS const&
277+
{
278+
ConvertFromEigenToNumPyMatrix<CLASS>(&$result, $1);
279+
}
280+
// Out: & (not yet implemented)
281+
%typemap(out, fragment="Eigen_Fragments") CLASS &
282+
{
283+
PyErr_SetString(PyExc_ValueError, "The output typemap for non-const reference is not yet implemented. Please report this problem to the developer.");
284+
}
285+
// Out: const* (not yet implemented)
286+
%typemap(out, fragment="Eigen_Fragments") CLASS const*
287+
{
288+
PyErr_SetString(PyExc_ValueError, "The output typemap for const pointer is not yet implemented. Please report this problem to the developer.");
289+
}
290+
// Out: * (not yet implemented)
291+
%typemap(out, fragment="Eigen_Fragments") CLASS *
292+
{
293+
PyErr_SetString(PyExc_ValueError, "The output typemap for non-const pointer is not yet implemented. Please report this problem to the developer.");
294+
}
295+
296+
%enddef

‎example.i‎

Lines changed: 39 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -3,49 +3,70 @@
33
%{
44
#define SWIG_FILE_WITH_INIT
55
#include <vector>
6-
#include <stdio.h> // TODO remove
76
#include "src/include/public_interface.hpp"
7+
#include "src/include/public_interface_eigen.hpp"
88
%}
99

1010
%include "numpy.i"
1111
%init %{
1212
import_array();
1313
%}
14+
%include <eigen.i>
1415

15-
// namespace std {
16-
// %template(LenVector) vector<length_t>;
17-
// %template(IntVector) vector<int>;
18-
// %template(DoubleVector) vector<double>;
19-
// %template(NeighborVector) vector<Neighbor>;
20-
//}
16+
// ------------------------------------------------
17+
// vector typemaps
18+
// ------------------------------------------------
19+
20+
%define %np_vector_typemaps(DTYPE, NPY_DPTYE)
2121

2222
namespace std {
23-
%typemap(out, fragment="NumPy_Fragments") vector<int> {
24-
// prepare resulting array
25-
// npy_intp dims[] = {$1->rows(), $1->cols()};
23+
// hmm...apparently telling SWIG to try to optimize this breaks it
24+
// %typemap(out, fragment="NumPy_Fragments", optimal="1") vector<DTYPE> {
25+
%typemap(out, fragment="NumPy_Fragments") vector<DTYPE> {
26+
// create python array of appropriate shape
2627
npy_intp sz = static_cast<npy_intp>($1.size());
2728
npy_intp dims[] = {sz};
28-
PyObject * out_array = PyArray_SimpleNew(1, dims, NPY_INT);
29+
PyObject* out_array = PyArray_SimpleNew(1, dims, NPY_DPTYE);
2930

3031
if (! out_array) {
3132
PyErr_SetString(PyExc_ValueError,
32-
"vector<int>: unable to create the output array.");
33+
"vector wrap: unable to create the output array.");
3334
return NULL;
3435
}
3536

36-
printf("calling std::vector<int> typemap\n");
37-
3837
// copy data from vect into numpy array
39-
int* out_data = (int*) array_data(out_array);
38+
DTYPE* out_data = (DTYPE*) array_data(out_array);
4039
for (size_t i = 0; i < sz; i++) {
41-
// out_data[i] = (*$1)[i];
42-
out_data[i] = static_cast<int>($1[i]);
40+
out_data[i] = static_cast<DTYPE>($1[i]);
4341
}
4442

4543
$result = out_array;
4644
}
4745
}
4846

47+
%enddef
48+
49+
%np_vector_typemaps(int, NPY_INT)
50+
%np_vector_typemaps(long, NPY_LONG)
51+
%np_vector_typemaps(float, NPY_FLOAT)
52+
%np_vector_typemaps(double, NPY_DOUBLE)
53+
// %np_vector_typemaps(SimpleStruct*, NPY_OBJECT) // breaks
54+
55+
// ------------------------------------------------
56+
// eigen typemaps
57+
// ------------------------------------------------
58+
59+
%eigen_typemaps(MatrixXd)
60+
%eigen_typemaps(VectorXd)
61+
// %eigen_typemaps(ArrayXd)
62+
%eigen_typemaps(MatrixXf)
63+
%eigen_typemaps(VectorXf)
64+
%eigen_typemaps(MatrixXi)
65+
%eigen_typemaps(VectorXi)
66+
67+
// ------------------------------------------------
68+
// raw c array typemaps
69+
// ------------------------------------------------
4970

5071
// apply numpy typemaps based on arg types + names
5172
%apply (double* INPLACE_ARRAY1, int DIM1) {(double* inVec, int len)};
@@ -64,6 +85,7 @@ namespace std {
6485
// %apply (int DIM1, double* ARGOUT_ARRAY1) {(int len, double* outVec)};
6586

6687
%include "src/include/public_interface.hpp"
88+
%include "src/include/public_interface_eigen.hpp"
6789

6890
// instantiate templates; note that these must be *after* the appropriate %include
6991
// NOTE: this doesn't seem to work with our custom "apply"s above, so I don't think

0 commit comments

Comments
 (0)