WIP: ENH: enable complex dtype support unconditionally

This also gets rid of usage of the deprecated IF macro in Cython code,
which resolves a bunch of build warnings.

Author: rgommers

pywt/_dwt.py

@@ -2,7 +2,6 @@
 
 import numpy as np
 
-from ._c99_config import _have_c99_complex
 from ._extensions._dwt import downcoef as _downcoef
 from ._extensions._dwt import dwt_axis, dwt_single, idwt_axis, idwt_single
 from ._extensions._dwt import dwt_coeff_len as _dwt_coeff_len
@@ -161,12 +160,6 @@ def dwt(data, wavelet, mode='symmetric', axis=-1):
     array([-0.70710678, -0.70710678, -0.70710678])
 
     """
-    if not _have_c99_complex and np.iscomplexobj(data):
-        data = np.asarray(data)
-        cA_r, cD_r = dwt(data.real, wavelet, mode, axis)
-        cA_i, cD_i = dwt(data.imag, wavelet, mode, axis)
-        return (cA_r + 1j*cA_i, cD_r + 1j*cD_i)
-
     # accept array_like input; make a copy to ensure a contiguous array
     dt = _check_dtype(data)
     data = np.asarray(data, dtype=dt, order='C')
@@ -241,17 +234,6 @@ def idwt(cA, cD, wavelet, mode='symmetric', axis=-1):
         raise ValueError("At least one coefficient parameter must be "
                          "specified.")
 
-    # for complex inputs: compute real and imaginary separately then combine
-    if not _have_c99_complex and (np.iscomplexobj(cA) or np.iscomplexobj(cD)):
-        if cA is None:
-            cD = np.asarray(cD)
-            cA = np.zeros_like(cD)
-        elif cD is None:
-            cA = np.asarray(cA)
-            cD = np.zeros_like(cA)
-        return (idwt(cA.real, cD.real, wavelet, mode, axis) +
-                1j*idwt(cA.imag, cD.imag, wavelet, mode, axis))
-
     if cA is not None:
         dt = _check_dtype(cA)
         cA = np.asarray(cA, dtype=dt, order='C')
@@ -328,9 +310,6 @@ def downcoef(part, data, wavelet, mode='symmetric', level=1):
     upcoef
 
     """
-    if not _have_c99_complex and np.iscomplexobj(data):
-        return (downcoef(part, data.real, wavelet, mode, level) +
-                1j*downcoef(part, data.imag, wavelet, mode, level))
     # accept array_like input; make a copy to ensure a contiguous array
     dt = _check_dtype(data)
     data = np.asarray(data, dtype=dt, order='C')
@@ -387,9 +366,6 @@ def upcoef(part, coeffs, wavelet, level=1, take=0):
     array([ 1.,  2.,  3.,  4.,  5.,  6.])
 
     """
-    if not _have_c99_complex and np.iscomplexobj(coeffs):
-        return (upcoef(part, coeffs.real, wavelet, level, take) +
-                1j*upcoef(part, coeffs.imag, wavelet, level, take))
     # accept array_like input; make a copy to ensure a contiguous array
     dt = _check_dtype(coeffs)
     coeffs = np.asarray(coeffs, dtype=dt, order='C')

pywt/_extensions/_dwt.pyx

@@ -3,25 +3,17 @@ cimport numpy as np
 
 np.import_array()
 
-include "config.pxi"
-
 ctypedef Py_ssize_t pywt_index_t
 
 ctypedef fused data_t:
     np.float32_t
     np.float64_t
 
-cdef int have_c99_complex
-IF HAVE_C99_CPLX:
-    ctypedef fused cdata_t:
-        np.float32_t
-        np.float64_t
-        np.complex64_t
-        np.complex128_t
-    have_c99_complex = 1
-ELSE:
-    ctypedef data_t cdata_t
-    have_c99_complex = 0
+ctypedef fused cdata_t:
+    np.float32_t
+    np.float64_t
+    np.complex64_t
+    np.complex128_t
 
 cdef public class Wavelet [type WaveletType, object WaveletObject]:
     cdef wavelet.DiscreteWavelet* w

pywt/_extensions/_pywt.pxd

@@ -10,8 +10,6 @@ cimport numpy as np
 from .common cimport pywt_index_t
 from ._pywt cimport c_wavelet_from_object, cdata_t, Wavelet, _check_dtype
 
-include "config.pxi"
-
 np.import_array()
 
 
@@ -99,21 +97,20 @@ def swt(cdata_t[::1] data, Wavelet wavelet, size_t level, size_t start_level,
                                 &cD[0], output_len, i)
             if retval < 0:
                 raise RuntimeError("C swt failed.")
-        IF HAVE_C99_CPLX:
-            if cdata_t is np.complex128_t:
-                cD = np.zeros(output_len, dtype=np.complex128)
-                with nogil:
-                    retval = c_wt.double_complex_swt_d(&data[0], data_size, wavelet.w,
-                                     &cD[0], output_len, i)
-                if retval < 0:
-                    raise RuntimeError("C swt failed.")
-            elif cdata_t is np.complex64_t:
-                cD = np.zeros(output_len, dtype=np.complex64)
-                with nogil:
-                    retval = c_wt.float_complex_swt_d(&data[0], data_size, wavelet.w,
-                                    &cD[0], output_len, i)
-                if retval < 0:
-                    raise RuntimeError("C swt failed.")
+        elif cdata_t is np.complex128_t:
+            cD = np.zeros(output_len, dtype=np.complex128)
+            with nogil:
+                retval = c_wt.double_complex_swt_d(&data[0], data_size, wavelet.w,
+                                 &cD[0], output_len, i)
+            if retval < 0:
+                raise RuntimeError("C swt failed.")
+        elif cdata_t is np.complex64_t:
+            cD = np.zeros(output_len, dtype=np.complex64)
+            with nogil:
+                retval = c_wt.float_complex_swt_d(&data[0], data_size, wavelet.w,
+                                &cD[0], output_len, i)
+            if retval < 0:
+                raise RuntimeError("C swt failed.")
 
         # alloc memory, decompose A
         if cdata_t is np.float64_t:
@@ -130,21 +127,20 @@ def swt(cdata_t[::1] data, Wavelet wavelet, size_t level, size_t start_level,
                                 &cA[0], output_len, i)
             if retval < 0:
                 raise RuntimeError("C swt failed.")
-        IF HAVE_C99_CPLX:
-            if cdata_t is np.complex128_t:
-                cA = np.zeros(output_len, dtype=np.complex128)
-                with nogil:
-                    retval = c_wt.double_complex_swt_a(&data[0], data_size, wavelet.w,
-                                     &cA[0], output_len, i)
-                if retval < 0:
-                    raise RuntimeError("C swt failed.")
-            elif cdata_t is np.complex64_t:
-                cA = np.zeros(output_len, dtype=np.complex64)
-                with nogil:
-                    retval = c_wt.float_complex_swt_a(&data[0], data_size, wavelet.w,
-                                    &cA[0], output_len, i)
-                if retval < 0:
-                    raise RuntimeError("C swt failed.")
+        elif cdata_t is np.complex128_t:
+            cA = np.zeros(output_len, dtype=np.complex128)
+            with nogil:
+                retval = c_wt.double_complex_swt_a(&data[0], data_size, wavelet.w,
+                                 &cA[0], output_len, i)
+            if retval < 0:
+                raise RuntimeError("C swt failed.")
+        elif cdata_t is np.complex64_t:
+            cA = np.zeros(output_len, dtype=np.complex64)
+            with nogil:
+                retval = c_wt.float_complex_swt_a(&data[0], data_size, wavelet.w,
+                                &cA[0], output_len, i)
+            if retval < 0:
+                raise RuntimeError("C swt failed.")
 
         data = cA
         if not trim_approx:
@@ -253,58 +249,57 @@ cpdef swt_axis(np.ndarray data, Wavelet wavelet, size_t level,
             if retval:
                 raise RuntimeError(
                     "C wavelet transform failed with error code %d" % retval)
+        elif data.dtype == np.complex128:
+            cA = np.zeros(output_shape, dtype=np.complex128)
+            with nogil:
+                retval = c_wt.double_complex_downcoef_axis(
+                    <double complex *> data.data, data_info,
+                    <double complex *> cA.data, output_info,
+                    wavelet.w, axis,
+                    common.COEF_APPROX, common.MODE_PERIODIZATION,
+                    i, common.SWT_TRANSFORM)
+            if retval:
+                raise RuntimeError(
+                    "C wavelet transform failed with error code %d" %
+                    retval)
+            cD = np.zeros(output_shape, dtype=np.complex128)
+            with nogil:
+                retval = c_wt.double_complex_downcoef_axis(
+                    <double complex *> data.data, data_info,
+                    <double complex *> cD.data, output_info,
+                    wavelet.w, axis,
+                    common.COEF_DETAIL, common.MODE_PERIODIZATION,
+                    i, common.SWT_TRANSFORM)
+            if retval:
+                raise RuntimeError(
+                    "C wavelet transform failed with error code %d" %
+                    retval)
+        elif data.dtype == np.complex64:
+            cA = np.zeros(output_shape, dtype=np.complex64)
+            with nogil:
+                retval = c_wt.float_complex_downcoef_axis(
+                    <float complex *> data.data, data_info,
+                    <float complex *> cA.data, output_info,
+                    wavelet.w, axis,
+                    common.COEF_APPROX, common.MODE_PERIODIZATION,
+                    i, common.SWT_TRANSFORM)
+            if retval:
+                raise RuntimeError(
+                    "C wavelet transform failed with error code %d" %
+                    retval)
+            cD = np.zeros(output_shape, dtype=np.complex64)
+            with nogil:
+                retval = c_wt.float_complex_downcoef_axis(
+                    <float complex *> data.data, data_info,
+                    <float complex *> cD.data, output_info,
+                    wavelet.w, axis,
+                    common.COEF_DETAIL, common.MODE_PERIODIZATION,
+                    i, common.SWT_TRANSFORM)
+            if retval:
+                raise RuntimeError(
+                    "C wavelet transform failed with error code %d" %
+                    retval)
 
-        IF HAVE_C99_CPLX:
-            if data.dtype == np.complex128:
-                cA = np.zeros(output_shape, dtype=np.complex128)
-                with nogil:
-                    retval = c_wt.double_complex_downcoef_axis(
-                        <double complex *> data.data, data_info,
-                        <double complex *> cA.data, output_info,
-                        wavelet.w, axis,
-                        common.COEF_APPROX, common.MODE_PERIODIZATION,
-                        i, common.SWT_TRANSFORM)
-                if retval:
-                    raise RuntimeError(
-                        "C wavelet transform failed with error code %d" %
-                        retval)
-                cD = np.zeros(output_shape, dtype=np.complex128)
-                with nogil:
-                    retval = c_wt.double_complex_downcoef_axis(
-                        <double complex *> data.data, data_info,
-                        <double complex *> cD.data, output_info,
-                        wavelet.w, axis,
-                        common.COEF_DETAIL, common.MODE_PERIODIZATION,
-                        i, common.SWT_TRANSFORM)
-                if retval:
-                    raise RuntimeError(
-                        "C wavelet transform failed with error code %d" %
-                        retval)
-            elif data.dtype == np.complex64:
-                cA = np.zeros(output_shape, dtype=np.complex64)
-                with nogil:
-                    retval = c_wt.float_complex_downcoef_axis(
-                        <float complex *> data.data, data_info,
-                        <float complex *> cA.data, output_info,
-                        wavelet.w, axis,
-                        common.COEF_APPROX, common.MODE_PERIODIZATION,
-                        i, common.SWT_TRANSFORM)
-                if retval:
-                    raise RuntimeError(
-                        "C wavelet transform failed with error code %d" %
-                        retval)
-                cD = np.zeros(output_shape, dtype=np.complex64)
-                with nogil:
-                    retval = c_wt.float_complex_downcoef_axis(
-                        <float complex *> data.data, data_info,
-                        <float complex *> cD.data, output_info,
-                        wavelet.w, axis,
-                        common.COEF_DETAIL, common.MODE_PERIODIZATION,
-                        i, common.SWT_TRANSFORM)
-                if retval:
-                    raise RuntimeError(
-                        "C wavelet transform failed with error code %d" %
-                        retval)
         if retval == -5:
             raise TypeError("Array must be floating point, not {}"
                             .format(data.dtype))

pywt/_extensions/_swt.pyx

@@ -8,11 +8,10 @@
 
 #pragma once
 
-#ifdef HAVE_C99_COMPLEX
-    /* For templating, we need typedefs without spaces for complex types. */
-    typedef float _Complex float_complex;
-    typedef double _Complex double_complex;
-#endif
+/* For templating, we need typedefs without spaces for complex types. */
+/* FIXME: needs more portable complex types here */
+typedef float _Complex float_complex;
+typedef double _Complex double_complex;
 
 /* ##### Typedefs ##### */
 

pywt/_extensions/c/common.h

@@ -20,19 +20,17 @@
 #undef REAL_TYPE
 #undef TYPE
 
-#ifdef HAVE_C99_COMPLEX
-    #define TYPE float_complex
-    #define REAL_TYPE float
-    #include "convolution.template.c"
-    #undef REAL_TYPE
-    #undef TYPE
+#define TYPE float_complex
+#define REAL_TYPE float
+#include "convolution.template.c"
+#undef REAL_TYPE
+#undef TYPE
 
-    #define TYPE double_complex
-    #define REAL_TYPE double
-    #include "convolution.template.c"
-    #undef REAL_TYPE
-    #undef TYPE
-#endif
+#define TYPE double_complex
+#define REAL_TYPE double
+#include "convolution.template.c"
+#undef REAL_TYPE
+#undef TYPE
 
 #endif /* REAL_TYPE */
 #endif /* TYPE */

pywt/_extensions/c/convolution.c

@@ -22,19 +22,17 @@
 #undef REAL_TYPE
 #undef TYPE
 
-#ifdef HAVE_C99_COMPLEX
-    #define TYPE float_complex
-    #define REAL_TYPE float
-    #include "convolution.template.h"
-    #undef REAL_TYPE
-    #undef TYPE
-
-    #define TYPE double_complex
-    #define REAL_TYPE double
-    #include "convolution.template.h"
-    #undef REAL_TYPE
-    #undef TYPE
-#endif
+#define TYPE float_complex
+#define REAL_TYPE float
+#include "convolution.template.h"
+#undef REAL_TYPE
+#undef TYPE
+
+#define TYPE double_complex
+#define REAL_TYPE double
+#include "convolution.template.h"
+#undef REAL_TYPE
+#undef TYPE
 
 #endif /* REAL_TYPE */
 #endif /* TYPE */

pywt/_extensions/c/convolution.h

@@ -20,19 +20,17 @@
 #undef REAL_TYPE
 #undef TYPE
 
-#ifdef HAVE_C99_COMPLEX
-    #define TYPE float_complex
-    #define REAL_TYPE float
-    #include "wt.template.c"
-    #undef REAL_TYPE
-    #undef TYPE
+#define TYPE float_complex
+#define REAL_TYPE float
+#include "wt.template.c"
+#undef REAL_TYPE
+#undef TYPE
 
-    #define TYPE double_complex
-    #define REAL_TYPE double
-    #include "wt.template.c"
-    #undef REAL_TYPE
-    #undef TYPE
-#endif
+#define TYPE double_complex
+#define REAL_TYPE double
+#include "wt.template.c"
+#undef REAL_TYPE
+#undef TYPE
 
 #endif /* REAL_TYPE */
 #endif /* TYPE */