From 5fb8959c2ff41edf2772544e215af24a98644c9f Mon Sep 17 00:00:00 2001 From: Ralf Gommers Date: Sat, 1 Feb 2025 19:58:05 +0100 Subject: [PATCH] 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. --- pywt/_dwt.py | 24 --- pywt/_extensions/_dwt.pyx | 293 +++++++++++++++---------------- pywt/_extensions/_pywt.pxd | 18 +- pywt/_extensions/_swt.pyx | 161 ++++++++--------- pywt/_extensions/c/common.h | 9 +- pywt/_extensions/c/convolution.c | 22 ++- pywt/_extensions/c/convolution.h | 24 ++- pywt/_extensions/c/wt.c | 22 ++- pywt/_extensions/c/wt.h | 24 ++- pywt/_extensions/c_wt.pxd | 139 +++++++-------- pywt/_extensions/config.pxi.in | 3 - pywt/_extensions/meson.build | 12 +- pywt/_pytesttester.py | 5 - 13 files changed, 344 insertions(+), 412 deletions(-) delete mode 100644 pywt/_extensions/config.pxi.in diff --git a/pywt/_dwt.py b/pywt/_dwt.py index 1af28d945..3d021827a 100644 --- a/pywt/_dwt.py +++ b/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') diff --git a/pywt/_extensions/_dwt.pyx b/pywt/_extensions/_dwt.pyx index 0fba4291d..e661f3b61 100644 --- a/pywt/_extensions/_dwt.pyx +++ b/pywt/_extensions/_dwt.pyx @@ -7,8 +7,6 @@ from ._pywt cimport _check_dtype cimport numpy as np import numpy as np -include "config.pxi" - np.import_array() cpdef dwt_max_level(size_t data_len, size_t filter_len): @@ -57,27 +55,26 @@ cpdef dwt_single(cdata_t[::1] data, Wavelet wavelet, MODE mode): if ( retval_a < 0 or retval_d < 0): raise RuntimeError("C dwt failed.") - IF HAVE_C99_CPLX: - if cdata_t is np.complex128_t: - cA = np.zeros(output_len, np.complex128) - cD = np.zeros(output_len, np.complex128) - with nogil: - retval_a = c_wt.double_complex_dec_a(&data[0], data_size, wavelet.w, - cA.data, output_len, mode) - retval_d = c_wt.double_complex_dec_d(&data[0], data_size, wavelet.w, - cD.data, output_len, mode) - if ( retval_a < 0 or retval_d < 0): - raise RuntimeError("C dwt failed.") - elif cdata_t is np.complex64_t: - cA = np.zeros(output_len, np.complex64) - cD = np.zeros(output_len, np.complex64) - with nogil: - retval_a = c_wt.float_complex_dec_a(&data[0], data_size, wavelet.w, - cA.data, output_len, mode) - retval_d = c_wt.float_complex_dec_d(&data[0], data_size, wavelet.w, - cD.data, output_len, mode) - if ( retval_a < 0 or retval_d < 0): - raise RuntimeError("C dwt failed.") + if cdata_t is np.complex128_t: + cA = np.zeros(output_len, np.complex128) + cD = np.zeros(output_len, np.complex128) + with nogil: + retval_a = c_wt.double_complex_dec_a(&data[0], data_size, wavelet.w, + cA.data, output_len, mode) + retval_d = c_wt.double_complex_dec_d(&data[0], data_size, wavelet.w, + cD.data, output_len, mode) + if ( retval_a < 0 or retval_d < 0): + raise RuntimeError("C dwt failed.") + elif cdata_t is np.complex64_t: + cA = np.zeros(output_len, np.complex64) + cD = np.zeros(output_len, np.complex64) + with nogil: + retval_a = c_wt.float_complex_dec_a(&data[0], data_size, wavelet.w, + cA.data, output_len, mode) + retval_d = c_wt.float_complex_dec_d(&data[0], data_size, wavelet.w, + cD.data, output_len, mode) + if ( retval_a < 0 or retval_d < 0): + raise RuntimeError("C dwt failed.") return (cA, cD) @@ -141,37 +138,37 @@ cpdef dwt_axis(np.ndarray data, Wavelet wavelet, MODE mode, unsigned int axis=0) 0, common.DWT_TRANSFORM) if retval: raise RuntimeError("C wavelet transform failed") - IF HAVE_C99_CPLX: - if data.dtype == np.complex64: - with nogil: - retval = c_wt.float_complex_downcoef_axis( data.data, data_info, - cA.data, output_info, - wavelet.w, axis, common.COEF_APPROX, mode, - 0, common.DWT_TRANSFORM) - if retval: - raise RuntimeError("C wavelet transform failed") - with nogil: - retval = c_wt.float_complex_downcoef_axis( data.data, data_info, - cD.data, output_info, - wavelet.w, axis, common.COEF_DETAIL, mode, - 0, common.DWT_TRANSFORM) - if retval: - raise RuntimeError("C wavelet transform failed") - elif data.dtype == np.complex128: - with nogil: - retval = c_wt.double_complex_downcoef_axis( data.data, data_info, - cA.data, output_info, - wavelet.w, axis, common.COEF_APPROX, mode, - 0, common.DWT_TRANSFORM) - if retval: - raise RuntimeError("C wavelet transform failed") - with nogil: - retval = c_wt.double_complex_downcoef_axis( data.data, data_info, - cD.data, output_info, - wavelet.w, axis, common.COEF_DETAIL, mode, + + if data.dtype == np.complex64: + with nogil: + retval = c_wt.float_complex_downcoef_axis( data.data, data_info, + cA.data, output_info, + wavelet.w, axis, common.COEF_APPROX, mode, + 0, common.DWT_TRANSFORM) + if retval: + raise RuntimeError("C wavelet transform failed") + with nogil: + retval = c_wt.float_complex_downcoef_axis( data.data, data_info, + cD.data, output_info, + wavelet.w, axis, common.COEF_DETAIL, mode, + 0, common.DWT_TRANSFORM) + if retval: + raise RuntimeError("C wavelet transform failed") + elif data.dtype == np.complex128: + with nogil: + retval = c_wt.double_complex_downcoef_axis( data.data, data_info, + cA.data, output_info, + wavelet.w, axis, common.COEF_APPROX, mode, 0, common.DWT_TRANSFORM) - if retval: - raise RuntimeError("C wavelet transform failed") + if retval: + raise RuntimeError("C wavelet transform failed") + with nogil: + retval = c_wt.double_complex_downcoef_axis( data.data, data_info, + cD.data, output_info, + wavelet.w, axis, common.COEF_DETAIL, mode, + 0, common.DWT_TRANSFORM) + if retval: + raise RuntimeError("C wavelet transform failed") if retval == -5: raise TypeError("Array must be floating point, not {}" @@ -220,25 +217,25 @@ cpdef idwt_single(np.ndarray cA, np.ndarray cD, Wavelet wavelet, MODE mode): wavelet.w, mode) if retval < 0: raise RuntimeError("C idwt failed.") - IF HAVE_C99_CPLX: - if cA.dtype == np.complex128: - rec = np.zeros(rec_len, dtype=np.complex128) - with nogil: - retval = c_wt.double_complex_idwt(cA.data, input_len, - cD.data, input_len, - rec.data, rec_len, - wavelet.w, mode) - if retval < 0: - raise RuntimeError("C idwt failed.") - elif cA.dtype == np.complex64: - rec = np.zeros(rec_len, dtype=np.complex64) - with nogil: - retval = c_wt.float_complex_idwt(cA.data, input_len, - cD.data, input_len, - rec.data, rec_len, - wavelet.w, mode) - if retval < 0: - raise RuntimeError("C idwt failed.") + + if cA.dtype == np.complex128: + rec = np.zeros(rec_len, dtype=np.complex128) + with nogil: + retval = c_wt.double_complex_idwt(cA.data, input_len, + cD.data, input_len, + rec.data, rec_len, + wavelet.w, mode) + if retval < 0: + raise RuntimeError("C idwt failed.") + elif cA.dtype == np.complex64: + rec = np.zeros(rec_len, dtype=np.complex64) + with nogil: + retval = c_wt.float_complex_idwt(cA.data, input_len, + cD.data, input_len, + rec.data, rec_len, + wavelet.w, mode) + if retval < 0: + raise RuntimeError("C idwt failed.") return rec @@ -311,23 +308,23 @@ cpdef idwt_axis(np.ndarray coefs_a, np.ndarray coefs_d, wavelet.w, axis, mode) if retval: raise RuntimeError("C inverse wavelet transform failed") - IF HAVE_C99_CPLX: - if output.dtype == np.complex128: - with nogil: - retval = c_wt.double_complex_idwt_axis( data_a, a_info_p, - data_d, d_info_p, - output.data, output_info, - wavelet.w, axis, mode) - if retval: - raise RuntimeError("C inverse wavelet transform failed") - elif output.dtype == np.complex64: - with nogil: - retval = c_wt.float_complex_idwt_axis( data_a, a_info_p, - data_d, d_info_p, - output.data, output_info, - wavelet.w, axis, mode) - if retval: - raise RuntimeError("C inverse wavelet transform failed") + + elif output.dtype == np.complex128: + with nogil: + retval = c_wt.double_complex_idwt_axis( data_a, a_info_p, + data_d, d_info_p, + output.data, output_info, + wavelet.w, axis, mode) + if retval: + raise RuntimeError("C inverse wavelet transform failed") + elif output.dtype == np.complex64: + with nogil: + retval = c_wt.float_complex_idwt_axis( data_a, a_info_p, + data_d, d_info_p, + output.data, output_info, + wavelet.w, axis, mode) + if retval: + raise RuntimeError("C inverse wavelet transform failed") if retval == -5: raise TypeError("Array must be floating point, not {}" @@ -386,35 +383,36 @@ cpdef upcoef(bint do_rec_a, cdata_t[::1] coeffs, Wavelet wavelet, int level, &rec[0], rec_len) if retval < 0: raise RuntimeError("C rec_d failed.") - IF HAVE_C99_CPLX: - if cdata_t is np.complex128_t: - rec = np.zeros(rec_len, dtype=np.complex128) - if do_rec_a or i > 0: - with nogil: - retval = c_wt.double_complex_rec_a(&coeffs[0], coeffs_size, wavelet.w, - &rec[0], rec_len) - if retval < 0: - raise RuntimeError("C rec_a failed.") - else: - with nogil: - retval = c_wt.double_complex_rec_d(&coeffs[0], coeffs_size, wavelet.w, - &rec[0], rec_len) - if retval < 0: - raise RuntimeError("C rec_d failed.") - elif cdata_t is np.complex64_t: - rec = np.zeros(rec_len, dtype=np.complex64) - if do_rec_a or i > 0: - with nogil: - retval = c_wt.float_complex_rec_a(&coeffs[0], coeffs_size, wavelet.w, - &rec[0], rec_len) - if retval < 0: - raise RuntimeError("C rec_a failed.") - else: - with nogil: - retval = c_wt.float_complex_rec_d(&coeffs[0], coeffs_size, wavelet.w, - &rec[0], rec_len) - if retval < 0: - raise RuntimeError("C rec_d failed.") + + elif cdata_t is np.complex128_t: + rec = np.zeros(rec_len, dtype=np.complex128) + if do_rec_a or i > 0: + with nogil: + retval = c_wt.double_complex_rec_a(&coeffs[0], coeffs_size, wavelet.w, + &rec[0], rec_len) + if retval < 0: + raise RuntimeError("C rec_a failed.") + else: + with nogil: + retval = c_wt.double_complex_rec_d(&coeffs[0], coeffs_size, wavelet.w, + &rec[0], rec_len) + if retval < 0: + raise RuntimeError("C rec_d failed.") + elif cdata_t is np.complex64_t: + rec = np.zeros(rec_len, dtype=np.complex64) + if do_rec_a or i > 0: + with nogil: + retval = c_wt.float_complex_rec_a(&coeffs[0], coeffs_size, wavelet.w, + &rec[0], rec_len) + if retval < 0: + raise RuntimeError("C rec_a failed.") + else: + with nogil: + retval = c_wt.float_complex_rec_d(&coeffs[0], coeffs_size, wavelet.w, + &rec[0], rec_len) + if retval < 0: + raise RuntimeError("C rec_d failed.") + # TODO: this algorithm needs some explaining coeffs = rec @@ -476,35 +474,36 @@ cpdef downcoef(bint do_dec_a, cdata_t[::1] data, Wavelet wavelet, MODE mode, int &coeffs[0], output_len, mode) if retval < 0: raise RuntimeError("C dec_d failed.") - IF HAVE_C99_CPLX: - if cdata_t is np.complex128_t: - coeffs = np.zeros(output_len, dtype=np.complex128) - if do_dec_a or (i < level - 1): - with nogil: - retval = c_wt.double_complex_dec_a(&data[0], data_size, wavelet.w, - &coeffs[0], output_len, mode) - if retval < 0: - raise RuntimeError("C dec_a failed.") - else: - with nogil: - retval = c_wt.double_complex_dec_d(&data[0], data_size, wavelet.w, - &coeffs[0], output_len, mode) - if retval < 0: - raise RuntimeError("C dec_d failed.") - elif cdata_t is np.complex64_t: - coeffs = np.zeros(output_len, dtype=np.complex64) - if do_dec_a or (i < level - 1): - with nogil: - retval = c_wt.float_complex_dec_a(&data[0], data_size, wavelet.w, - &coeffs[0], output_len, mode) - if retval < 0: - raise RuntimeError("C dec_a failed.") - else: - with nogil: - retval = c_wt.float_complex_dec_d(&data[0], data_size, wavelet.w, - &coeffs[0], output_len, mode) - if retval < 0: - raise RuntimeError("C dec_d failed.") + + elif cdata_t is np.complex128_t: + coeffs = np.zeros(output_len, dtype=np.complex128) + if do_dec_a or (i < level - 1): + with nogil: + retval = c_wt.double_complex_dec_a(&data[0], data_size, wavelet.w, + &coeffs[0], output_len, mode) + if retval < 0: + raise RuntimeError("C dec_a failed.") + else: + with nogil: + retval = c_wt.double_complex_dec_d(&data[0], data_size, wavelet.w, + &coeffs[0], output_len, mode) + if retval < 0: + raise RuntimeError("C dec_d failed.") + elif cdata_t is np.complex64_t: + coeffs = np.zeros(output_len, dtype=np.complex64) + if do_dec_a or (i < level - 1): + with nogil: + retval = c_wt.float_complex_dec_a(&data[0], data_size, wavelet.w, + &coeffs[0], output_len, mode) + if retval < 0: + raise RuntimeError("C dec_a failed.") + else: + with nogil: + retval = c_wt.float_complex_dec_d(&data[0], data_size, wavelet.w, + &coeffs[0], output_len, mode) + if retval < 0: + raise RuntimeError("C dec_d failed.") + data = coeffs return coeffs diff --git a/pywt/_extensions/_pywt.pxd b/pywt/_extensions/_pywt.pxd index 8485495ef..7ad9f7f07 100644 --- a/pywt/_extensions/_pywt.pxd +++ b/pywt/_extensions/_pywt.pxd @@ -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 diff --git a/pywt/_extensions/_swt.pyx b/pywt/_extensions/_swt.pyx index 26e9e3349..590e5d243 100644 --- a/pywt/_extensions/_swt.pyx +++ b/pywt/_extensions/_swt.pyx @@ -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( + data.data, data_info, + 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( + data.data, data_info, + 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( + data.data, data_info, + 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( + data.data, data_info, + 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( - data.data, data_info, - 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( - data.data, data_info, - 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( - data.data, data_info, - 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( - data.data, data_info, - 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)) diff --git a/pywt/_extensions/c/common.h b/pywt/_extensions/c/common.h index 7fa05f4c1..a5f48b507 100644 --- a/pywt/_extensions/c/common.h +++ b/pywt/_extensions/c/common.h @@ -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 ##### */ diff --git a/pywt/_extensions/c/convolution.c b/pywt/_extensions/c/convolution.c index d64e275db..6f85660e8 100644 --- a/pywt/_extensions/c/convolution.c +++ b/pywt/_extensions/c/convolution.c @@ -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 */ diff --git a/pywt/_extensions/c/convolution.h b/pywt/_extensions/c/convolution.h index 03e55cb21..4402a5691 100644 --- a/pywt/_extensions/c/convolution.h +++ b/pywt/_extensions/c/convolution.h @@ -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 */ diff --git a/pywt/_extensions/c/wt.c b/pywt/_extensions/c/wt.c index 3abb123a2..cee04bf4b 100644 --- a/pywt/_extensions/c/wt.c +++ b/pywt/_extensions/c/wt.c @@ -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 */ diff --git a/pywt/_extensions/c/wt.h b/pywt/_extensions/c/wt.h index a2e2ea8c7..daad7ca04 100644 --- a/pywt/_extensions/c/wt.h +++ b/pywt/_extensions/c/wt.h @@ -24,19 +24,17 @@ #undef REAL_TYPE #undef TYPE -#ifdef HAVE_C99_COMPLEX - #define TYPE float_complex - #define REAL_TYPE float - #include "wt.template.h" - #undef REAL_TYPE - #undef TYPE - - #define TYPE double_complex - #define REAL_TYPE double - #include "wt.template.h" - #undef REAL_TYPE - #undef TYPE -#endif +#define TYPE float_complex +#define REAL_TYPE float +#include "wt.template.h" +#undef REAL_TYPE +#undef TYPE + +#define TYPE double_complex +#define REAL_TYPE double +#include "wt.template.h" +#undef REAL_TYPE +#undef TYPE #endif /* REAL_TYPE */ #endif /* TYPE */ diff --git a/pywt/_extensions/c_wt.pxd b/pywt/_extensions/c_wt.pxd index b0d33854b..ecbf3d60b 100644 --- a/pywt/_extensions/c_wt.pxd +++ b/pywt/_extensions/c_wt.pxd @@ -7,8 +7,6 @@ from .common cimport (MODE, pywt_index_t, ArrayInfo, Coefficient, DiscreteTransformType) from .wavelet cimport DiscreteWavelet, ContinuousWavelet -include "config.pxi" - cdef extern from "c/wt.h": # Cython does not know the 'restrict' keyword cdef int double_downcoef_axis(const double * const input, const ArrayInfo input_info, @@ -86,83 +84,82 @@ cdef extern from "c/wt.h": cdef int float_swt_d(const float * const input, pywt_index_t input_len, const DiscreteWavelet * const wavelet, float * const output, pywt_index_t output_len, int level) nogil - IF HAVE_C99_CPLX: - # complex variants only available if the compiler supports C99 complex - cdef int double_complex_downcoef_axis(const double complex * const input, const ArrayInfo input_info, - double complex * const output, const ArrayInfo output_info, - const DiscreteWavelet * const wavelet, const size_t axis, - const Coefficient detail, const MODE dwt_mode, - const size_t swt_level, - const DiscreteTransformType transform) nogil - cdef int double_complex_idwt_axis(const double complex * const coefs_a, const ArrayInfo * const a_info, - const double complex * const coefs_d, const ArrayInfo * const d_info, + # complex variants only available if the compiler supports C99 complex + cdef int double_complex_downcoef_axis(const double complex * const input, const ArrayInfo input_info, double complex * const output, const ArrayInfo output_info, const DiscreteWavelet * const wavelet, const size_t axis, - const MODE mode) nogil - cdef int double_complex_dec_a(const double complex * const input, const size_t input_len, - const DiscreteWavelet * const wavelet, - double complex * const output, const size_t output_len, - const MODE mode) nogil - cdef int double_complex_dec_d(const double complex * const input, const size_t input_len, - const DiscreteWavelet * const wavelet, - double complex * const output, const size_t output_len, + const Coefficient detail, const MODE dwt_mode, + const size_t swt_level, + const DiscreteTransformType transform) nogil + cdef int double_complex_idwt_axis(const double complex * const coefs_a, const ArrayInfo * const a_info, + const double complex * const coefs_d, const ArrayInfo * const d_info, + double complex * const output, const ArrayInfo output_info, + const DiscreteWavelet * const wavelet, const size_t axis, const MODE mode) nogil + cdef int double_complex_dec_a(const double complex * const input, const size_t input_len, + const DiscreteWavelet * const wavelet, + double complex * const output, const size_t output_len, + const MODE mode) nogil + cdef int double_complex_dec_d(const double complex * const input, const size_t input_len, + const DiscreteWavelet * const wavelet, + double complex * const output, const size_t output_len, + const MODE mode) nogil + + cdef int double_complex_rec_a(const double complex * const coeffs_a, const size_t coeffs_len, + const DiscreteWavelet * const wavelet, + double complex * const output, const size_t output_len) nogil + cdef int double_complex_rec_d(const double complex * const coeffs_d, const size_t coeffs_len, + const DiscreteWavelet * const wavelet, + double complex * const output, const size_t output_len) nogil + + cdef int double_complex_idwt(const double complex * const coeffs_a, const size_t coeffs_a_len, + const double complex * const coeffs_d, const size_t coeffs_d_len, + double complex * const output, const size_t output_len, + const DiscreteWavelet * const wavelet, const MODE mode) nogil + + cdef int double_complex_swt_a(const double complex * const input, size_t input_len, const DiscreteWavelet * const wavelet, + double complex * const output, size_t output_len, int level) nogil + cdef int double_complex_swt_d(const double complex * const input, size_t input_len, const DiscreteWavelet * const wavelet, + double complex * const output, size_t output_len, int level) nogil - cdef int double_complex_rec_a(const double complex * const coeffs_a, const size_t coeffs_len, - const DiscreteWavelet * const wavelet, - double complex * const output, const size_t output_len) nogil - cdef int double_complex_rec_d(const double complex * const coeffs_d, const size_t coeffs_len, - const DiscreteWavelet * const wavelet, - double complex * const output, const size_t output_len) nogil - - cdef int double_complex_idwt(const double complex * const coeffs_a, const size_t coeffs_a_len, - const double complex * const coeffs_d, const size_t coeffs_d_len, - double complex * const output, const size_t output_len, - const DiscreteWavelet * const wavelet, const MODE mode) nogil - - cdef int double_complex_swt_a(const double complex * const input, size_t input_len, const DiscreteWavelet * const wavelet, - double complex * const output, size_t output_len, int level) nogil - cdef int double_complex_swt_d(const double complex * const input, size_t input_len, const DiscreteWavelet * const wavelet, - double complex * const output, size_t output_len, int level) nogil - - - - cdef int float_complex_downcoef_axis(const float complex * const input, const ArrayInfo input_info, - float complex * const output, const ArrayInfo output_info, - const DiscreteWavelet * const wavelet, const size_t axis, - const Coefficient detail, const MODE dwt_mode, - const size_t swt_level, - const DiscreteTransformType transform) nogil - cdef int float_complex_idwt_axis(const float complex * const coefs_a, const ArrayInfo * const a_info, - const float complex * const coefs_d, const ArrayInfo * const d_info, + + + cdef int float_complex_downcoef_axis(const float complex * const input, const ArrayInfo input_info, float complex * const output, const ArrayInfo output_info, const DiscreteWavelet * const wavelet, const size_t axis, - const MODE mode) nogil - cdef int float_complex_dec_a(const float complex * const input, const size_t input_len, - const DiscreteWavelet * const wavelet, - float complex * const output, const size_t output_len, - const MODE mode) nogil - cdef int float_complex_dec_d(const float complex * const input, const size_t input_len, - const DiscreteWavelet * const wavelet, - float complex * const output, const size_t output_len, + const Coefficient detail, const MODE dwt_mode, + const size_t swt_level, + const DiscreteTransformType transform) nogil + cdef int float_complex_idwt_axis(const float complex * const coefs_a, const ArrayInfo * const a_info, + const float complex * const coefs_d, const ArrayInfo * const d_info, + float complex * const output, const ArrayInfo output_info, + const DiscreteWavelet * const wavelet, const size_t axis, const MODE mode) nogil + cdef int float_complex_dec_a(const float complex * const input, const size_t input_len, + const DiscreteWavelet * const wavelet, + float complex * const output, const size_t output_len, + const MODE mode) nogil + cdef int float_complex_dec_d(const float complex * const input, const size_t input_len, + const DiscreteWavelet * const wavelet, + float complex * const output, const size_t output_len, + const MODE mode) nogil + + cdef int float_complex_rec_a(const float complex * const coeffs_a, const size_t coeffs_len, + const DiscreteWavelet * const wavelet, + float complex * const output, const size_t output_len) nogil + cdef int float_complex_rec_d(const float complex * const coeffs_d, const size_t coeffs_len, + const DiscreteWavelet * const wavelet, + float complex * const output, const size_t output_len) nogil + + cdef int float_complex_idwt(const float complex * const coeffs_a, const size_t coeffs_a_len, + const float complex * const coeffs_d, const size_t coeffs_d_len, + float complex * const output, const size_t output_len, + const DiscreteWavelet * const wavelet, const MODE mode) nogil - cdef int float_complex_rec_a(const float complex * const coeffs_a, const size_t coeffs_len, - const DiscreteWavelet * const wavelet, - float complex * const output, const size_t output_len) nogil - cdef int float_complex_rec_d(const float complex * const coeffs_d, const size_t coeffs_len, - const DiscreteWavelet * const wavelet, - float complex * const output, const size_t output_len) nogil - - cdef int float_complex_idwt(const float complex * const coeffs_a, const size_t coeffs_a_len, - const float complex * const coeffs_d, const size_t coeffs_d_len, - float complex * const output, const size_t output_len, - const DiscreteWavelet * const wavelet, const MODE mode) nogil - - cdef int float_complex_swt_a(const float complex * const input, size_t input_len, const DiscreteWavelet* const wavelet, - float complex * const output, size_t output_len, int level) nogil - cdef int float_complex_swt_d(const float complex * const input, size_t input_len, const DiscreteWavelet* const wavelet, - float complex * const output, size_t output_len, int level) nogil + cdef int float_complex_swt_a(const float complex * const input, size_t input_len, const DiscreteWavelet* const wavelet, + float complex * const output, size_t output_len, int level) nogil + cdef int float_complex_swt_d(const float complex * const input, size_t input_len, const DiscreteWavelet* const wavelet, + float complex * const output, size_t output_len, int level) nogil cdef extern from "c/cwt.h": # Cython does not know the 'restrict' keyword diff --git a/pywt/_extensions/config.pxi.in b/pywt/_extensions/config.pxi.in deleted file mode 100644 index f2d7e1794..000000000 --- a/pywt/_extensions/config.pxi.in +++ /dev/null @@ -1,3 +0,0 @@ -# Autogenerated file containing Cython compile-time defines - -DEF HAVE_C99_CPLX = @have_c99_complex_cy@ diff --git a/pywt/_extensions/meson.build b/pywt/_extensions/meson.build index 8ebf71ea2..30f14368d 100644 --- a/pywt/_extensions/meson.build +++ b/pywt/_extensions/meson.build @@ -42,13 +42,6 @@ print(incdir) np_dep = declare_dependency(include_directories: inc_np, compile_args: numpy_nodepr_api) endif -config_pxi = configure_file( - input: 'config.pxi.in', - output: 'config.pxi', - configuration: cdata, - install: false -) - sources = [ 'c/common.c', 'c/convolution.c', @@ -70,10 +63,7 @@ _cython_tree = [ ] # avoid compiler warnings: tell Cython to use C99 complex types -c_args = ['-DPY_EXTENSION'] -if have_c99_complex - c_args += ['-DHAVE_C99_COMPLEX', '-DCYTHON_CCOMPLEX=1'] -endif +c_args = ['-DPY_EXTENSION', '-DCYTHON_CCOMPLEX=1'] libc_wt = static_library('c_wt', sources, diff --git a/pywt/_pytesttester.py b/pywt/_pytesttester.py index 6344fa8e3..298ea3a09 100644 --- a/pywt/_pytesttester.py +++ b/pywt/_pytesttester.py @@ -33,12 +33,7 @@ def _show_pywt_info(): import pywt - from pywt._c99_config import _have_c99_complex print(f"PyWavelets version {pywt.__version__}") - if _have_c99_complex: - print("Compiled with C99 complex support.") - else: - print("Compiled without C99 complex support.") class PytestTester: