Added support for complex numbers.

eigency/conversions.pxd

@@ -30,3 +30,13 @@ cdef api np.ndarray[char, ndim=2] ndarray_char_F(char *data, long rows, long col
 cdef api np.ndarray[char, ndim=2] ndarray_copy_char_C(const char *data, long rows, long cols, long outer_stride, long inner_stride)
 cdef api np.ndarray[char, ndim=2] ndarray_copy_char_F(const char *data, long rows, long cols, long outer_stride, long inner_stride)
 
+cdef api np.ndarray[np.complex128_t, ndim=2] ndarray_complex_double_C(np.complex128_t *data, long rows, long cols, long outer_stride, long inner_stride)
+cdef api np.ndarray[np.complex128_t, ndim=2] ndarray_complex_double_F(np.complex128_t *data, long rows, long cols, long outer_stride, long inner_stride)
+cdef api np.ndarray[np.complex128_t, ndim=2] ndarray_copy_complex_double_C(const np.complex128_t *data, long rows, long cols, long outer_stride, long inner_stride)
+cdef api np.ndarray[np.complex128_t, ndim=2] ndarray_copy_complex_double_F(const np.complex128_t *data, long rows, long cols, long outer_stride, long inner_stride)
+
+cdef api np.ndarray[np.complex64_t, ndim=2] ndarray_complex_float_C(np.complex64_t *data, long rows, long cols, long outer_stride, long inner_stride)
+cdef api np.ndarray[np.complex64_t, ndim=2] ndarray_complex_float_F(np.complex64_t *data, long rows, long cols, long outer_stride, long inner_stride)
+cdef api np.ndarray[np.complex64_t, ndim=2] ndarray_copy_complex_float_C(const np.complex64_t *data, long rows, long cols, long outer_stride, long inner_stride)
+cdef api np.ndarray[np.complex64_t, ndim=2] ndarray_copy_complex_float_F(const np.complex64_t *data, long rows, long cols, long outer_stride, long inner_stride)
+

eigency/conversions.pyx

@@ -140,3 +140,48 @@ cdef np.ndarray[char, ndim=2] ndarray_copy_char_F(const char *data, long rows, l
     return np.copy(as_strided(np.asarray(mem_view, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
 
 
+@cython.boundscheck(False)
+cdef np.ndarray[np.complex128_t, ndim=2] ndarray_complex_double_C(np.complex128_t *data, long rows, long cols, long row_stride, long col_stride):
+    cdef np.complex128_t[:,:] mem_view = <np.complex128_t[:rows,:cols]>data
+    cdef int itemsize = np.dtype('complex128').itemsize
+    return as_strided(np.asarray(mem_view, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+@cython.boundscheck(False)
+cdef np.ndarray[np.complex128_t, ndim=2] ndarray_complex_double_F(np.complex128_t *data, long rows, long cols, long row_stride, long col_stride):
+    cdef np.complex128_t[::1,:] mem_view = <np.complex128_t[:rows:1,:cols]>data
+    cdef int itemsize = np.dtype('complex128').itemsize
+    return as_strided(np.asarray(mem_view, order="F"), strides=[row_stride*itemsize, col_stride*itemsize])
+
+@cython.boundscheck(False)
+cdef np.ndarray[np.complex128_t, ndim=2] ndarray_copy_complex_double_C(const np.complex128_t *data, long rows, long cols, long row_stride, long col_stride):
+    cdef np.complex128_t[:,:] mem_view = <np.complex128_t[:rows,:cols]>data
+    cdef int itemsize = np.dtype('complex128').itemsize
+    return np.copy(as_strided(np.asarray(mem_view, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+@cython.boundscheck(False)
+cdef np.ndarray[np.complex128_t, ndim=2] ndarray_copy_complex_double_F(const np.complex128_t *data, long rows, long cols, long row_stride, long col_stride):
+    cdef np.complex128_t[::1,:] mem_view = <np.complex128_t[:rows:1,:cols]>data
+    cdef int itemsize = np.dtype('complex128').itemsize
+    return np.copy(as_strided(np.asarray(mem_view, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
+
+
+@cython.boundscheck(False)
+cdef np.ndarray[np.complex64_t, ndim=2] ndarray_complex_float_C(np.complex64_t *data, long rows, long cols, long row_stride, long col_stride):
+    cdef np.complex64_t[:,:] mem_view = <np.complex64_t[:rows,:cols]>data
+    cdef int itemsize = np.dtype('complex64').itemsize
+    return as_strided(np.asarray(mem_view, order="C"), strides=[row_stride*itemsize, col_stride*itemsize])
+@cython.boundscheck(False)
+cdef np.ndarray[np.complex64_t, ndim=2] ndarray_complex_float_F(np.complex64_t *data, long rows, long cols, long row_stride, long col_stride):
+    cdef np.complex64_t[::1,:] mem_view = <np.complex64_t[:rows:1,:cols]>data
+    cdef int itemsize = np.dtype('complex64').itemsize
+    return as_strided(np.asarray(mem_view, order="F"), strides=[row_stride*itemsize, col_stride*itemsize])
+
+@cython.boundscheck(False)
+cdef np.ndarray[np.complex64_t, ndim=2] ndarray_copy_complex_float_C(const np.complex64_t *data, long rows, long cols, long row_stride, long col_stride):
+    cdef np.complex64_t[:,:] mem_view = <np.complex64_t[:rows,:cols]>data
+    cdef int itemsize = np.dtype('complex64').itemsize
+    return np.copy(as_strided(np.asarray(mem_view, order="C"), strides=[row_stride*itemsize, col_stride*itemsize]))
+@cython.boundscheck(False)
+cdef np.ndarray[np.complex64_t, ndim=2] ndarray_copy_complex_float_F(const np.complex64_t *data, long rows, long cols, long row_stride, long col_stride):
+    cdef np.complex64_t[::1,:] mem_view = <np.complex64_t[:rows:1,:cols]>data
+    cdef int itemsize = np.dtype('complex64').itemsize
+    return np.copy(as_strided(np.asarray(mem_view, order="F"), strides=[row_stride*itemsize, col_stride*itemsize]))
+

eigency/eigency_cpp.h

@@ -2,6 +2,11 @@
 
 #include <iostream>
 #include <stdexcept>
+#include <complex>
+
+typedef ::std::complex< double > __pyx_t_double_complex;
+typedef ::std::complex< float > __pyx_t_float_complex;
+
 #include "conversions_api.h"
 
 #ifndef EIGENCY_CPP
@@ -117,6 +122,37 @@ inline PyArrayObject *_ndarray_copy<char>(const char *data, long rows, long cols
         return ndarray_copy_char_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
 }
 
+template<>
+inline PyArrayObject *_ndarray_view<std::complex<double> >(std::complex<double> *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
+    if (is_row_major)
+        return ndarray_complex_double_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
+    else
+        return ndarray_complex_double_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
+}
+template<>
+inline PyArrayObject *_ndarray_copy<std::complex<double> >(const std::complex<double> *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
+    if (is_row_major)
+        return ndarray_copy_complex_double_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
+    else
+        return ndarray_copy_complex_double_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
+}
+
+template<>
+inline PyArrayObject *_ndarray_view<std::complex<float> >(std::complex<float> *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
+    if (is_row_major)
+        return ndarray_complex_float_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
+    else
+        return ndarray_complex_float_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
+}
+template<>
+inline PyArrayObject *_ndarray_copy<std::complex<float> >(const std::complex<float> *data, long rows, long cols, bool is_row_major, long outer_stride, long inner_stride) {
+    if (is_row_major)
+        return ndarray_copy_complex_float_C(data, rows, cols, outer_stride>0?outer_stride:cols, inner_stride>0?inner_stride:1);
+    else
+        return ndarray_copy_complex_float_F(data, rows, cols, inner_stride>0?inner_stride:1, outer_stride>0?outer_stride:rows);
+}
+
+
 template <typename Derived>
 inline PyArrayObject *ndarray(Eigen::PlainObjectBase<Derived> &m) {
     import_eigency__conversions();

setup.py

@@ -14,7 +14,8 @@
 
 extensions = [
     Extension("eigency.conversions", ["eigency/conversions.pyx"],
-              include_dirs = [np.get_include(), __eigen_dir__]
+              include_dirs = [np.get_include(), __eigen_dir__],
+              language="c++"
     ),
     Extension("eigency.core", ["eigency/core.pyx"],
               include_dirs = [np.get_include(), __eigen_dir__],
@@ -30,13 +31,13 @@
 
 dist = setup(
     name = __package_name__,
-    version = "1.67",
+    version = "1.68",
     description = "Cython interface between the numpy arrays and the Matrix/Array classes of the Eigen C++ library",
     long_description=long_description,
     author = "Wouter Boomsma",
     author_email = "wb@bio.ku.dk",
     url = "https://github.com/wouterboomsma/eigency",
-    download_url = "https://github.com/wouterboomsma/eigency/tarball/1.67",
+    download_url = "https://github.com/wouterboomsma/eigency/tarball/1.68",
     ext_modules = cythonize(extensions),
     packages = [__package_name__]
 )

tests/eigency_tests/eigency_tests.pyx

@@ -21,6 +21,8 @@ cdef extern from "eigency_tests/eigency_tests_cpp.h":
 
      cdef void _function_w_mat_arg "function_w_mat_arg"(Map[MatrixXd] &)
 
+     cdef void _function_w_complex_mat_arg "function_w_complex_mat_arg"(Map[MatrixXcd] &)
+
      cdef void _function_w_fullspec_arg "function_w_fullspec_arg" (FlattenedMap[Array, double, Dynamic, _1] &)
 
      cdef VectorXd _function_w_vec_retval "function_w_vec_retval" ()
@@ -71,6 +73,10 @@ def function_w_vec_arg_no_map2(np.ndarray array):
 def function_w_mat_arg(np.ndarray array):
     return _function_w_mat_arg(Map[MatrixXd](array))
 
+# Function with complex matrix argument. 
+def function_w_complex_mat_arg(np.ndarray array):
+    return _function_w_complex_mat_arg(Map[MatrixXcd](array))
+
 # Function using a full Map specification, rather than the convenience typedefs
 # Note that since cython does not support nested fused types, the Map has been
 # flattened to include all arguments at once

tests/eigency_tests/eigency_tests_cpp.cpp

@@ -25,6 +25,10 @@ void function_w_mat_arg(Eigen::Map<Eigen::MatrixXd> &mat) {
     mat(0,0) = 0.;
 }
 
+void function_w_complex_mat_arg(Eigen::Map<Eigen::MatrixXcd> &mat) {
+    mat(0,0) = 0.+0.j;
+}
+
 void function_w_fullspec_arg(Eigen::Map<Eigen::Array<double, Eigen::Dynamic, 1> > &vec) {
     vec[0] = 0.;
 }

tests/eigency_tests/eigency_tests_cpp.h

@@ -13,6 +13,8 @@ void function_w_vec_arg_no_map2(const Eigen::VectorXd &vec);
 
 void function_w_mat_arg(Eigen::Map<Eigen::MatrixXd> &mat);
 
+void function_w_complex_mat_arg(Eigen::Map<Eigen::MatrixXcd> &mat);
+
 void function_w_fullspec_arg(Eigen::Map<Eigen::Array<double, Eigen::Dynamic, 1> > &vec);
 
 Eigen::VectorXd function_w_vec_retval();

tests/run_tests.py

@@ -37,6 +37,12 @@ def test_function_w_mat_arg(self):
         # Shared memory test: Verify that first entry was set to 0 by C++ code.
         self.assertAlmostEqual(x[0], 0.)        
 
+    def test_function_w_complex_mat_arg(self):
+        x = np.array([1.+1.j, 2.+1.j, 3.+1.j, 4.+1.j])
+        eigency_tests.function_w_complex_mat_arg(x.reshape([2,2]))
+        # Shared memory test: Verify that first entry was set to 0.+0.j by C++ code.
+        self.assertAlmostEqual(x[0], 0.+0.j)        
+
     def test_funcion_w_fullspec_arg(self):
         x = np.array([1., 2., 3., 4.])
         eigency_tests.function_w_fullspec_arg(x)