diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 2d95cc91..b14ca609 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -41,6 +41,12 @@ jobs: # Linux arm64 - runs-on: ubuntu-22.04-arm python-version: "3.14" + # Linux arm64, MSVC-style struct complex types forced, to exercise + # the operator-free complex code path on a C99-complex compiler + # (see PYWT_TEST_STRUCT_COMPLEX in pywt/_extensions/c/common.h) + - runs-on: ubuntu-22.04-arm + python-version: "3.14" + OPTIONS_NAME: "struct-complex" # Linux amd64 - runs-on: ubuntu-latest python-version: "3.12" @@ -100,6 +106,8 @@ jobs: elif [ "${REFGUIDE_CHECK}" == "1" ]; then pip install sphinx numpydoc scipy-doctest pip install . -v + elif [ "${OPTIONS_NAME}" == "struct-complex" ]; then + pip install . -v -Csetup-args=-Dc_args=-DPYWT_TEST_STRUCT_COMPLEX else pip install . -v fi diff --git a/.gitignore b/.gitignore index cd7ca61c..c64d84c2 100644 --- a/.gitignore +++ b/.gitignore @@ -36,8 +36,6 @@ doc/source/regression/*.ipynb # Project working files # Expanded Cython pywt/_extensions/*.[ch] -pywt/_extensions/_c99_config.py -pywt/_extensions/config.pxi cythonize.dat pywt/version.py build.log diff --git a/meson.build b/meson.build index 5da23bac..4bbfdafe 100644 --- a/meson.build +++ b/meson.build @@ -11,7 +11,7 @@ project( ], ) -cython = find_program('cython') +cython = find_program('cython', version: '>=3.2.5') py = import('python').find_installation(pure: false) py_dep = py.dependency() diff --git a/pyproject.toml b/pyproject.toml index 6aca3680..5ab200fb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,7 +8,7 @@ build-backend = "mesonpy" requires = [ "meson-python>=0.18.0", - "Cython>=3.1.3", + "Cython>=3.2.5", # numpy requirement for wheel builds for distribution on PyPI - building # against 2.x yields wheels that are also compatible with numpy 1.x at diff --git a/pywt/_c99_config.py.in b/pywt/_c99_config.py.in deleted file mode 100644 index ad2dc41c..00000000 --- a/pywt/_c99_config.py.in +++ /dev/null @@ -1,3 +0,0 @@ -# Autogenerated file containing compile-time definitions - -_have_c99_complex = @have_c99_complex@ diff --git a/pywt/_dwt.py b/pywt/_dwt.py index cf01b61a..151b8787 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 d4177615..9a8c0167 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): @@ -56,28 +54,26 @@ cpdef dwt_single(const cdata_t[::1] data, Wavelet wavelet, MODE mode): cD.data, output_len, 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.") + elif 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 +137,36 @@ 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, + elif 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 +215,24 @@ 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.") + elif 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 +305,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 +380,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 +471,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 8485495e..7ad9f7f0 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 26e9e334..7b8db479 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 7fa05f4c..25117f64 100644 --- a/pywt/_extensions/c/common.h +++ b/pywt/_extensions/c/common.h @@ -8,12 +8,6 @@ #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 - /* ##### Typedefs ##### */ #ifdef PY_EXTENSION @@ -39,6 +33,191 @@ #include #include +/* ##### Complex types ##### + * + * Portable complex types and accessors, following NumPy's approach (see + * numpy/_core/include/numpy/npy_common.h and npy_math.h): MSVC has no C99 + * `_Complex` and provides the struct-based `_Fcomplex`/`_Dcomplex` instead, + * with no arithmetic operators. All complex arithmetic in the C sources + * therefore goes through the inline helpers below, which only rely on the + * layout (two consecutive reals) that all representations share. + * + * For templating, the typedef names must not contain spaces. + * + * PYWT_TEST_STRUCT_COMPLEX is a testing hook: it swaps in plain structs with + * the same layout and no operators, so that MSVC-like semantics can be + * compiled and tested with any compiler (used by the struct-complex CI job). + * The structs match the ones NumPy's headers use when compiled as C++, where + * C99 `_Complex` is also unavailable. + */ + +#ifdef PYWT_TEST_STRUCT_COMPLEX + +typedef struct +{ + float _Val[2]; +} float_complex; + +typedef struct +{ + double _Val[2]; +} double_complex; + +#else + +#include + +#if defined(_MSC_VER) && !defined(__INTEL_COMPILER) && !defined(__INTEL_LLVM_COMPILER) +typedef _Fcomplex float_complex; +typedef _Dcomplex double_complex; +#else +typedef float _Complex float_complex; +typedef double _Complex double_complex; +#endif + +#endif + +static inline float pywt_crealf(const float_complex z) +{ +#ifdef PYWT_TEST_STRUCT_COMPLEX + return z._Val[0]; +#else + return crealf(z); +#endif +} + +static inline float pywt_cimagf(const float_complex z) +{ +#ifdef PYWT_TEST_STRUCT_COMPLEX + return z._Val[1]; +#else + return cimagf(z); +#endif +} + +static inline void pywt_csetrealf(float_complex *z, const float r) +{ + ((float *) z)[0] = r; +} + +static inline void pywt_csetimagf(float_complex *z, const float i) +{ + ((float *) z)[1] = i; +} + +static inline double pywt_creal(const double_complex z) +{ +#ifdef PYWT_TEST_STRUCT_COMPLEX + return z._Val[0]; +#else + return creal(z); +#endif +} + +static inline double pywt_cimag(const double_complex z) +{ +#ifdef PYWT_TEST_STRUCT_COMPLEX + return z._Val[1]; +#else + return cimag(z); +#endif +} + +static inline void pywt_csetreal(double_complex *z, const double r) +{ + ((double *) z)[0] = r; +} + +static inline void pywt_csetimag(double_complex *z, const double i) +{ + ((double *) z)[1] = i; +} + +/* Arithmetic helpers used by the type-templated convolution kernels; the + * CAT(TYPE, _op) spellings must exist for every TYPE the templates are + * instantiated with. The complex variants only need a real second factor in + * _mul because the wavelet filters are always real. */ + +static inline float float_zero(void) { return 0.0f; } +static inline float float_add(const float a, const float b) { return a + b; } +static inline float float_sub(const float a, const float b) { return a - b; } +static inline float float_mul(const float a, const float b) { return a * b; } + +static inline double double_zero(void) { return 0.0; } +static inline double double_add(const double a, const double b) { return a + b; } +static inline double double_sub(const double a, const double b) { return a - b; } +static inline double double_mul(const double a, const double b) { return a * b; } + +static inline float_complex float_complex_zero(void) +{ + float_complex z; + pywt_csetrealf(&z, 0.0f); + pywt_csetimagf(&z, 0.0f); + return z; +} + +static inline float_complex float_complex_add(const float_complex a, + const float_complex b) +{ + float_complex z; + pywt_csetrealf(&z, pywt_crealf(a) + pywt_crealf(b)); + pywt_csetimagf(&z, pywt_cimagf(a) + pywt_cimagf(b)); + return z; +} + +static inline float_complex float_complex_sub(const float_complex a, + const float_complex b) +{ + float_complex z; + pywt_csetrealf(&z, pywt_crealf(a) - pywt_crealf(b)); + pywt_csetimagf(&z, pywt_cimagf(a) - pywt_cimagf(b)); + return z; +} + +static inline float_complex float_complex_mul(const float a, + const float_complex b) +{ + float_complex z; + pywt_csetrealf(&z, a * pywt_crealf(b)); + pywt_csetimagf(&z, a * pywt_cimagf(b)); + return z; +} + +static inline double_complex double_complex_zero(void) +{ + double_complex z; + pywt_csetreal(&z, 0.0); + pywt_csetimag(&z, 0.0); + return z; +} + +static inline double_complex double_complex_add(const double_complex a, + const double_complex b) +{ + double_complex z; + pywt_csetreal(&z, pywt_creal(a) + pywt_creal(b)); + pywt_csetimag(&z, pywt_cimag(a) + pywt_cimag(b)); + return z; +} + +static inline double_complex double_complex_sub(const double_complex a, + const double_complex b) +{ + double_complex z; + pywt_csetreal(&z, pywt_creal(a) - pywt_creal(b)); + pywt_csetimag(&z, pywt_cimag(a) - pywt_cimag(b)); + return z; +} + +static inline double_complex double_complex_mul(const double a, + const double_complex b) +{ + double_complex z; + pywt_csetreal(&z, a * pywt_creal(b)); + pywt_csetimag(&z, a * pywt_cimag(b)); + return z; +} + /* standard c memory management */ #define wtmalloc(size) malloc(size) #define wtfree(ptr) free(ptr) diff --git a/pywt/_extensions/c/convolution.c b/pywt/_extensions/c/convolution.c index d64e275d..6f85660e 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 03e55cb2..4402a569 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/convolution.template.c b/pywt/_extensions/c/convolution.template.c index 4153200e..b6003a63 100644 --- a/pywt/_extensions/c/convolution.template.c +++ b/pywt/_extensions/c/convolution.template.c @@ -23,6 +23,18 @@ #define restrict __restrict__ #endif +/* Type-generic arithmetic, implemented in common.h. MSVC's struct-based + * complex types have no arithmetic operators, so every operation on TYPE + * values must go through these. ADD_PRODUCT(s, f, x) means `s += f * x`, + * SUB_PRODUCT(s, f, x) means `s -= f * x`. The multiplier `f` is always + * REAL_TYPE; the filters stay real for complex data. */ +#define ZERO CAT(TYPE, _zero)() +#define ADD(a, b) CAT(TYPE, _add)((a), (b)) +#define SUB(a, b) CAT(TYPE, _sub)((a), (b)) +#define MUL(f, x) CAT(TYPE, _mul)((f), (x)) +#define ADD_PRODUCT(s, f, x) ((s) = ADD((s), MUL((f), (x)))) +#define SUB_PRODUCT(s, f, x) ((s) = SUB((s), MUL((f), (x)))) + /* This file contains several functions for computing the convolution of a * signal with a filter. The general scheme is: * output[o] = sum(filter[j] * input[i-j] for j = [0..F) and i = [0..N)) @@ -50,74 +62,74 @@ int CAT(TYPE, _downsampling_convolution_periodization)(const TYPE * const restri const size_t padding = (step - (N % step)) % step; for (; i < F && i < N; i += step, ++o) { - TYPE sum = 0; + TYPE sum = ZERO; size_t j; size_t k_start = 0; for (j = 0; j <= i; j += fstep) - sum += filter[j] * input[i-j]; + ADD_PRODUCT(sum, filter[j], input[i-j]); if (fstep > 1) k_start = j - (i + 1); while (j < F){ size_t k; for (k = k_start; k < padding && j < F; k += fstep, j += fstep) - sum += filter[j] * input[N-1]; + ADD_PRODUCT(sum, filter[j], input[N-1]); for (k = k_start; k < N && j < F; k += fstep, j += fstep) - sum += filter[j] * input[N-1-k]; + ADD_PRODUCT(sum, filter[j], input[N-1-k]); } output[o] = sum; } for(; i < N; i+=step, ++o){ - TYPE sum = 0; + TYPE sum = ZERO; size_t j; for(j = 0; j < F; j += fstep) - sum += input[i-j]*filter[j]; + ADD_PRODUCT(sum, filter[j], input[i-j]); output[o] = sum; } for (; i < F && i < N + F/2; i += step, ++o) { - TYPE sum = 0; + TYPE sum = ZERO; size_t j = 0; size_t k_start = 0; while (i-j >= N){ size_t k; // for simplicity, not using fstep here for (k = 0; k < padding && i-j >= N; ++k, ++j) - sum += filter[i-N-j] * input[N-1]; + ADD_PRODUCT(sum, filter[i-N-j], input[N-1]); for (k = 0; k < N && i-j >= N; ++k, ++j) - sum += filter[i-N-j] * input[k]; + ADD_PRODUCT(sum, filter[i-N-j], input[k]); } if (fstep > 1) j += (fstep - j % fstep) % fstep; // move to next non-zero entry for (; j <= i; j += fstep) - sum += filter[j] * input[i-j]; + ADD_PRODUCT(sum, filter[j], input[i-j]); if (fstep > 1) k_start = j - (i + 1); while (j < F){ size_t k; for (k = k_start; k < padding && j < F; k += fstep, j += fstep) - sum += filter[j] * input[N-1]; + ADD_PRODUCT(sum, filter[j], input[N-1]); for (k = k_start; k < N && j < F; k += fstep, j += fstep) - sum += filter[j] * input[N-1-k]; + ADD_PRODUCT(sum, filter[j], input[N-1-k]); } output[o] = sum; } for(; i < N + F/2; i += step, ++o){ - TYPE sum = 0; + TYPE sum = ZERO; size_t j = 0; while (i-j >= N){ // for simplicity, not using fstep here size_t k; for (k = 0; k < padding && i-j >= N; ++k, ++j) - sum += filter[i-N-j] * input[N-1]; + ADD_PRODUCT(sum, filter[i-N-j], input[N-1]); for (k = 0; k < N && i-j >= N; ++k, ++j) - sum += filter[i-N-j] * input[k]; + ADD_PRODUCT(sum, filter[i-N-j], input[k]); } if (fstep > 1) j += (fstep - j % fstep) % fstep; // move to next non-zero entry for (; j < F; j += fstep) - sum += filter[j] * input[i-j]; + ADD_PRODUCT(sum, filter[j], input[i-j]); output[o] = sum; } return 0; @@ -144,19 +156,19 @@ int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, cons // left boundary overhang for(; i < F && i < N; i+=step, ++o){ - TYPE sum = 0; + TYPE sum = ZERO; size_t j; for(j = 0; j <= i; ++j) - sum += filter[j]*input[i-j]; + ADD_PRODUCT(sum, filter[j], input[i-j]); switch(mode) { case MODE_SYMMETRIC: while (j < F){ size_t k; for(k = 0; k < N && j < F; ++j, ++k) - sum += filter[j]*input[k]; + ADD_PRODUCT(sum, filter[j], input[k]); for(k = 0; k < N && j < F; ++k, ++j) - sum += filter[j]*input[N-1-k]; + ADD_PRODUCT(sum, filter[j], input[N-1-k]); } break; case MODE_ANTISYMMETRIC: @@ -164,34 +176,34 @@ int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, cons while (j < F){ size_t k; for(k = 0; k < N && j < F; ++j, ++k) - sum -= filter[j]*input[k]; + SUB_PRODUCT(sum, filter[j], input[k]); for(k = 0; k < N && j < F; ++k, ++j) - sum += filter[j]*input[N-1-k]; + ADD_PRODUCT(sum, filter[j], input[N-1-k]); } break; case MODE_REFLECT: while (j < F){ size_t k; for(k = 1; k < N && j < F; ++j, ++k) - sum += filter[j]*input[k]; + ADD_PRODUCT(sum, filter[j], input[k]); for(k = 1; k < N && j < F; ++k, ++j) - sum += filter[j]*input[N-1-k]; + ADD_PRODUCT(sum, filter[j], input[N-1-k]); } break; case MODE_ANTIREFLECT:{ // whole-sample anti-symmetric size_t k; TYPE le = input[0]; // current left edge value - TYPE tmp = 0; + TYPE tmp = ZERO; while (j < F) { for(k = 1; k < N && j < F; ++j, ++k){ - tmp = le - (input[k] - input[0]); - sum += filter[j]*tmp; + tmp = SUB(le, SUB(input[k], input[0])); + ADD_PRODUCT(sum, filter[j], tmp); } le = tmp; for(k = 1; k < N && j < F; ++j, ++k){ - tmp = le + (input[N-1-k] - input[N-1]); - sum += filter[j]*tmp; + tmp = ADD(le, SUB(input[N-1-k], input[N-1])); + ADD_PRODUCT(sum, filter[j], tmp); } le = tmp; } @@ -199,19 +211,19 @@ int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, cons } case MODE_CONSTANT_EDGE: for(; j < F; ++j) - sum += filter[j]*input[0]; + ADD_PRODUCT(sum, filter[j], input[0]); break; case MODE_SMOOTH:{ size_t k; for(k = 1; j < F; ++j, ++k) - sum += filter[j]*(input[0] + k * (input[0] - input[1])); + ADD_PRODUCT(sum, filter[j], ADD(input[0], MUL((REAL_TYPE)k, SUB(input[0], input[1])))); break; } case MODE_PERIODIC: while (j < F){ size_t k; for(k = 0; k < N && j < F; ++k, ++j) - sum += filter[j]*input[N-1-k]; + ADD_PRODUCT(sum, filter[j], input[N-1-k]); } break; case MODE_ZEROPAD: @@ -223,16 +235,16 @@ int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, cons // center (if input equal or wider than filter: N >= F) for(; i < N; i+=step, ++o){ - TYPE sum = 0; + TYPE sum = ZERO; size_t j; for(j = 0; j < F; ++j) - sum += input[i-j]*filter[j]; + ADD_PRODUCT(sum, filter[j], input[i-j]); output[o] = sum; } // center (if filter is wider than input: F > N) for(; i < F; i+=step, ++o){ - TYPE sum = 0; + TYPE sum = ZERO; size_t j = 0; switch(mode) { @@ -244,9 +256,9 @@ int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, cons while (i - j >= N){ size_t k; for(k = 0; k < N && i-j >= N; ++j, ++k) - sum += filter[i-N-j]*input[N-1-k]; + ADD_PRODUCT(sum, filter[i-N-j], input[N-1-k]); for(k = 0; k < N && i-j >= N; ++j, ++k) - sum += filter[i-N-j]*input[k]; + ADD_PRODUCT(sum, filter[i-N-j], input[k]); } break; case MODE_ANTISYMMETRIC: @@ -254,34 +266,34 @@ int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, cons while (i - j >= N){ size_t k; for(k = 0; k < N && i-j >= N; ++j, ++k) - sum -= filter[i-N-j]*input[N-1-k]; + SUB_PRODUCT(sum, filter[i-N-j], input[N-1-k]); for(k = 0; k < N && i-j >= N; ++j, ++k) - sum += filter[i-N-j]*input[k]; + ADD_PRODUCT(sum, filter[i-N-j], input[k]); } break; case MODE_REFLECT: while (i - j >= N){ size_t k; for(k = 1; k < N && i-j >= N; ++j, ++k) - sum += filter[i-N-j]*input[N-1-k]; + ADD_PRODUCT(sum, filter[i-N-j], input[N-1-k]); for(k = 1; k < N && i-j >= N; ++j, ++k) - sum += filter[i-N-j]*input[k]; + ADD_PRODUCT(sum, filter[i-N-j], input[k]); } break; case MODE_ANTIREFLECT:{ // whole-sample anti-symmetric size_t k; TYPE re = input[N-1]; // current right edge value - TYPE tmp = 0; + TYPE tmp = ZERO; while (i - j >= N) { for(k = 1; k < N && i-j >= N; ++j, ++k){ - tmp = re - (input[N-1-k] - input[N-1]); - sum += filter[i-N-j]*tmp; + tmp = SUB(re, SUB(input[N-1-k], input[N-1])); + ADD_PRODUCT(sum, filter[i-N-j], tmp); } re = tmp; for(k = 1; k < N && i-j >= N; ++j, ++k){ - tmp = re + (input[k] - input[0]); - sum += filter[i-N-j]*tmp; + tmp = ADD(re, SUB(input[k], input[0])); + ADD_PRODUCT(sum, filter[i-N-j], tmp); } re = tmp; } @@ -289,19 +301,19 @@ int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, cons } case MODE_CONSTANT_EDGE: for(; i-j >= N; ++j) - sum += filter[j]*input[N-1]; + ADD_PRODUCT(sum, filter[j], input[N-1]); break; case MODE_SMOOTH:{ size_t k; for(k = i - N + 1; i-j >= N; ++j, --k) - sum += filter[j]*(input[N-1] + k * (input[N-1] - input[N-2])); + ADD_PRODUCT(sum, filter[j], ADD(input[N-1], MUL((REAL_TYPE)k, SUB(input[N-1], input[N-2])))); break; } case MODE_PERIODIC: while (i-j >= N){ size_t k; for (k = 0; k < N && i-j >= N; ++j, ++k) - sum += filter[i-N-j]*input[k]; + ADD_PRODUCT(sum, filter[i-N-j], input[k]); } break; case MODE_ZEROPAD: @@ -311,16 +323,16 @@ int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, cons } for(; j <= i; ++j) - sum += filter[j]*input[i-j]; + ADD_PRODUCT(sum, filter[j], input[i-j]); switch(mode) { case MODE_SYMMETRIC: while (j < F){ size_t k; for(k = 0; k < N && j < F; ++j, ++k) - sum += filter[j]*input[k]; + ADD_PRODUCT(sum, filter[j], input[k]); for(k = 0; k < N && j < F; ++k, ++j) - sum += filter[j]*input[N-1-k]; + ADD_PRODUCT(sum, filter[j], input[N-1-k]); } break; case MODE_ANTISYMMETRIC: @@ -328,34 +340,34 @@ int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, cons while (j < F){ size_t k; for(k = 0; k < N && j < F; ++j, ++k) - sum -= filter[j]*input[k]; + SUB_PRODUCT(sum, filter[j], input[k]); for(k = 0; k < N && j < F; ++k, ++j) - sum += filter[j]*input[N-1-k]; + ADD_PRODUCT(sum, filter[j], input[N-1-k]); } break; case MODE_REFLECT: while (j < F){ size_t k; for(k = 1; k < N && j < F; ++j, ++k) - sum += filter[j]*input[k]; + ADD_PRODUCT(sum, filter[j], input[k]); for(k = 1; k < N && j < F; ++k, ++j) - sum += filter[j]*input[N-1-k]; + ADD_PRODUCT(sum, filter[j], input[N-1-k]); } break; case MODE_ANTIREFLECT:{ // whole-sample anti-symmetric size_t k; TYPE le = input[0]; // current left edge value - TYPE tmp = 0; + TYPE tmp = ZERO; while (j < F) { for(k = 1; k < N && j < F; ++j, ++k){ - tmp = le - (input[k] - input[0]); - sum += filter[j]*tmp; + tmp = SUB(le, SUB(input[k], input[0])); + ADD_PRODUCT(sum, filter[j], tmp); } le = tmp; for(k = 1; k < N && j < F; ++j, ++k){ - tmp = le + (input[N-1-k] - input[N-1]); - sum += filter[j]*tmp; + tmp = ADD(le, SUB(input[N-1-k], input[N-1])); + ADD_PRODUCT(sum, filter[j], tmp); } le = tmp; } @@ -363,19 +375,19 @@ int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, cons } case MODE_CONSTANT_EDGE: for(; j < F; ++j) - sum += filter[j]*input[0]; + ADD_PRODUCT(sum, filter[j], input[0]); break; case MODE_SMOOTH:{ size_t k; for(k = 1; j < F; ++j, ++k) - sum += filter[j]*(input[0] + k * (input[0] - input[1])); + ADD_PRODUCT(sum, filter[j], ADD(input[0], MUL((REAL_TYPE)k, SUB(input[0], input[1])))); break; } case MODE_PERIODIC: while (j < F){ size_t k; for(k = 0; k < N && j < F; ++k, ++j) - sum += filter[j]*input[N-1-k]; + ADD_PRODUCT(sum, filter[j], input[N-1-k]); } break; case MODE_ZEROPAD: @@ -387,7 +399,7 @@ int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, cons // right boundary overhang for(; i < N+F-1; i += step, ++o){ - TYPE sum = 0; + TYPE sum = ZERO; size_t j = 0; switch(mode) { case MODE_SYMMETRIC: @@ -395,9 +407,9 @@ int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, cons while (i - j >= N){ size_t k; for(k = 0; k < N && i-j >= N; ++j, ++k) - sum += filter[i-N-j]*input[N-1-k]; + ADD_PRODUCT(sum, filter[i-N-j], input[N-1-k]); for(k = 0; k < N && i-j >= N; ++j, ++k) - sum += filter[i-N-j]*input[k]; + ADD_PRODUCT(sum, filter[i-N-j], input[k]); } break; case MODE_ANTISYMMETRIC: @@ -405,36 +417,36 @@ int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, cons while (i - j >= N){ size_t k; for(k = 0; k < N && i-j >= N; ++j, ++k) - sum -= filter[i-N-j]*input[N-1-k]; + SUB_PRODUCT(sum, filter[i-N-j], input[N-1-k]); for(k = 0; k < N && i-j >= N; ++j, ++k) - sum += filter[i-N-j]*input[k]; + ADD_PRODUCT(sum, filter[i-N-j], input[k]); } break; case MODE_REFLECT: while (i - j >= N){ size_t k; for(k = 1; k < N && i-j >= N; ++j, ++k) - sum += filter[i-N-j]*input[N-1-k]; + ADD_PRODUCT(sum, filter[i-N-j], input[N-1-k]); for(k = 1; k < N && i-j >= N; ++j, ++k) - sum += filter[i-N-j]*input[k]; + ADD_PRODUCT(sum, filter[i-N-j], input[k]); } break; case MODE_ANTIREFLECT:{ // whole-sample anti-symmetric size_t k; TYPE re = input[N-1]; //current right edge value - TYPE tmp = 0; + TYPE tmp = ZERO; while (i - j >= N) { //first reflection for(k = 1; k < N && i-j >= N; ++j, ++k){ - tmp = re - (input[N-1-k] - input[N-1]); - sum += filter[i-N-j]*tmp; + tmp = SUB(re, SUB(input[N-1-k], input[N-1])); + ADD_PRODUCT(sum, filter[i-N-j], tmp); } re = tmp; //second reflection for(k = 1; k < N && i-j >= N; ++j, ++k){ - tmp = re + (input[k] - input[0]); - sum += filter[i-N-j]*tmp; + tmp = ADD(re, SUB(input[k], input[0])); + ADD_PRODUCT(sum, filter[i-N-j], tmp); } re = tmp; } @@ -442,19 +454,19 @@ int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, cons } case MODE_CONSTANT_EDGE: for(; i-j >= N; ++j) - sum += filter[j]*input[N-1]; + ADD_PRODUCT(sum, filter[j], input[N-1]); break; case MODE_SMOOTH:{ size_t k; for(k = i - N + 1; i-j >= N; ++j, --k) - sum += filter[j]*(input[N-1] + k * (input[N-1] - input[N-2])); + ADD_PRODUCT(sum, filter[j], ADD(input[N-1], MUL((REAL_TYPE)k, SUB(input[N-1], input[N-2])))); break; } case MODE_PERIODIC: while (i-j >= N){ size_t k; for (k = 0; k < N && i-j >= N; ++j, ++k) - sum += filter[i-N-j]*input[k]; + ADD_PRODUCT(sum, filter[i-N-j], input[k]); } break; case MODE_ZEROPAD: @@ -463,7 +475,7 @@ int CAT(TYPE, _downsampling_convolution)(const TYPE * const restrict input, cons break; } for(; j < F; ++j) - sum += filter[j]*input[i-j]; + ADD_PRODUCT(sum, filter[j], input[i-j]); output[o] = sum; } return 0; @@ -492,32 +504,32 @@ int CAT(TYPE, _upsampling_convolution_full)(const TYPE * const restrict input, c for(; i < N && i < F/2; ++i, o += 2){ size_t j; for(j = 0; j <= i; ++j){ - output[o] += filter[j*2] * input[i-j]; - output[o+1] += filter[j*2+1] * input[i-j]; + ADD_PRODUCT(output[o], filter[j*2], input[i-j]); + ADD_PRODUCT(output[o+1], filter[j*2+1], input[i-j]); } } for(; i < N; ++i, o += 2){ size_t j; for(j = 0; j < F/2; ++j){ - output[o] += filter[j*2] * input[i-j]; - output[o+1] += filter[j*2+1] * input[i-j]; + ADD_PRODUCT(output[o], filter[j*2], input[i-j]); + ADD_PRODUCT(output[o+1], filter[j*2+1], input[i-j]); } } for(; i < F/2; ++i, o += 2){ size_t j; for(j = i-(N-1); j <= i; ++j){ - output[o] += filter[j*2] * input[i-j]; - output[o+1] += filter[j*2+1] * input[i-j]; + ADD_PRODUCT(output[o], filter[j*2], input[i-j]); + ADD_PRODUCT(output[o+1], filter[j*2+1], input[i-j]); } } for(; i < N+F/2; ++i, o += 2){ size_t j; for(j = i-(N-1); j < F/2; ++j){ - output[o] += filter[j*2] * input[i-j]; - output[o+1] += filter[j*2+1] * input[i-j]; + ADD_PRODUCT(output[o], filter[j*2], input[i-j]); + ADD_PRODUCT(output[o+1], filter[j*2+1], input[i-j]); } } return 0; @@ -545,19 +557,19 @@ static int CAT(TYPE, _upsampling_convolution_valid_sf_periodization)(const TYPE while(j <= start-1){ size_t k; for (k = 0; k < N && j <= start-1; ++k, ++j){ - output[2*N-1] += filter[2*(start-1-j)] * input[k]; - output[0] += filter[2*(start-1-j)+1] * input[k]; + ADD_PRODUCT(output[2*N-1], filter[2*(start-1-j)], input[k]); + ADD_PRODUCT(output[0], filter[2*(start-1-j)+1], input[k]); } } for (; j <= N+start-1 && j < F/2; ++j){ - output[2*N-1] += filter[2*j] * input[N+start-1-j]; - output[0] += filter[2*j+1] * input[N+start-1-j]; + ADD_PRODUCT(output[2*N-1], filter[2*j], input[N+start-1-j]); + ADD_PRODUCT(output[0], filter[2*j+1], input[N+start-1-j]); } while (j < F / 2){ size_t k; for (k = 0; k < N && j < F/2; ++k, ++j){ - output[2*N-1] += filter[2*j] * input[N-1-k]; - output[0] += filter[2*j+1] * input[N-1-k]; + ADD_PRODUCT(output[2*N-1], filter[2*j], input[N-1-k]); + ADD_PRODUCT(output[0], filter[2*j+1], input[N-1-k]); } } @@ -567,14 +579,14 @@ static int CAT(TYPE, _upsampling_convolution_valid_sf_periodization)(const TYPE for (; i < F/2 && i < N; ++i, o += 2){ size_t j = 0; for(; j <= i; ++j){ - output[o] += filter[2*j] * input[i-j]; - output[o+1] += filter[2*j+1] * input[i-j]; + ADD_PRODUCT(output[o], filter[2*j], input[i-j]); + ADD_PRODUCT(output[o+1], filter[2*j+1], input[i-j]); } while (j < F/2){ size_t k; for(k = 0; k < N && j < F/2; ++k, ++j){ - output[o] += filter[2*j] * input[N-1-k]; - output[o+1] += filter[2*j+1] * input[N-1-k]; + ADD_PRODUCT(output[o], filter[2*j], input[N-1-k]); + ADD_PRODUCT(output[o+1], filter[2*j+1], input[N-1-k]); } } } @@ -582,8 +594,8 @@ static int CAT(TYPE, _upsampling_convolution_valid_sf_periodization)(const TYPE for (; i < N; ++i, o += 2){ size_t j; for(j = 0; j < F/2; ++j){ - output[o] += filter[2*j] * input[i-j]; - output[o+1] += filter[2*j+1] * input[i-j]; + ADD_PRODUCT(output[o], filter[2*j], input[i-j]); + ADD_PRODUCT(output[o+1], filter[2*j+1], input[i-j]); } } @@ -592,19 +604,19 @@ static int CAT(TYPE, _upsampling_convolution_valid_sf_periodization)(const TYPE while(i-j >= N){ size_t k; for (k = 0; k < N && i-j >= N; ++k, ++j){ - output[o] += filter[2*(i-N-j)] * input[k]; - output[o+1] += filter[2*(i-N-j)+1] * input[k]; + ADD_PRODUCT(output[o], filter[2*(i-N-j)], input[k]); + ADD_PRODUCT(output[o+1], filter[2*(i-N-j)+1], input[k]); } } for (; j <= i && j < F/2; ++j){ - output[o] += filter[2*j] * input[i-j]; - output[o+1] += filter[2*j+1] * input[i-j]; + ADD_PRODUCT(output[o], filter[2*j], input[i-j]); + ADD_PRODUCT(output[o+1], filter[2*j+1], input[i-j]); } while (j < F / 2){ size_t k; for (k = 0; k < N && j < F/2; ++k, ++j){ - output[o] += filter[2*j] * input[N-1-k]; - output[o+1] += filter[2*j+1] * input[N-1-k]; + ADD_PRODUCT(output[o], filter[2*j], input[N-1-k]); + ADD_PRODUCT(output[o+1], filter[2*j+1], input[N-1-k]); } } } @@ -614,13 +626,13 @@ static int CAT(TYPE, _upsampling_convolution_valid_sf_periodization)(const TYPE while(i-j >= N){ size_t k; for (k = 0; k < N && i-j >= N; ++k, ++j){ - output[o] += filter[2*(i-N-j)] * input[k]; - output[o+1] += filter[2*(i-N-j)+1] * input[k]; + ADD_PRODUCT(output[o], filter[2*(i-N-j)], input[k]); + ADD_PRODUCT(output[o+1], filter[2*(i-N-j)+1], input[k]); } } for (; j <= i && j < F/2; ++j){ - output[o] += filter[2*j] * input[i-j]; - output[o+1] += filter[2*j+1] * input[i-j]; + ADD_PRODUCT(output[o], filter[2*j], input[i-j]); + ADD_PRODUCT(output[o+1], filter[2*j+1], input[i-j]); } } @@ -653,15 +665,15 @@ int CAT(TYPE, _upsampling_convolution_valid_sf)(const TYPE * const restrict inpu { size_t o, i; for(o = 0, i = F/2 - 1; i < N; ++i, o += 2){ - TYPE sum_even = 0; - TYPE sum_odd = 0; + TYPE sum_even = ZERO; + TYPE sum_odd = ZERO; size_t j; for(j = 0; j < F/2; ++j){ - sum_even += filter[j*2] * input[i-j]; - sum_odd += filter[j*2+1] * input[i-j]; + ADD_PRODUCT(sum_even, filter[j*2], input[i-j]); + ADD_PRODUCT(sum_odd, filter[j*2+1], input[i-j]); } - output[o] += sum_even; - output[o+1] += sum_odd; + output[o] = ADD(output[o], sum_even); + output[o+1] = ADD(output[o+1], sum_odd); } } return 0; @@ -676,6 +688,12 @@ int CAT(TYPE, _upsampled_filter_convolution)(const TYPE * const restrict input, return -1; } +#undef ZERO +#undef ADD +#undef SUB +#undef MUL +#undef ADD_PRODUCT +#undef SUB_PRODUCT #undef restrict #endif /* REAL_TYPE */ #endif /* TYPE */ diff --git a/pywt/_extensions/c/wt.c b/pywt/_extensions/c/wt.c index 3abb123a..cee04bf4 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 a2e2ea8c..daad7ca0 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 b0d33854..ee86dab2 100644 --- a/pywt/_extensions/c_wt.pxd +++ b/pywt/_extensions/c_wt.pxd @@ -7,7 +7,15 @@ from .common cimport (MODE, pywt_index_t, ArrayInfo, Coefficient, DiscreteTransformType) from .wavelet cimport DiscreteWavelet, ContinuousWavelet -include "config.pxi" +cdef extern from "c/common.h": + # Opaque complex types; _Fcomplex/_Dcomplex on MSVC, C99 complex + # elsewhere (see common.h). They match the layout of the corresponding + # NumPy and Cython complex types, but are distinct types to Cython, so + # pointers must be cast explicitly at call sites. + ctypedef struct float_complex: + pass + ctypedef struct double_complex: + pass cdef extern from "c/wt.h": # Cython does not know the 'restrict' keyword @@ -86,83 +94,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, - double complex * const output, const ArrayInfo output_info, + # Complex variants; the filters stay real, only the data is 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, - float complex * const output, const ArrayInfo output_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 f2d7e179..00000000 --- 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 8ebf71ea..8f9808e6 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', @@ -57,7 +50,12 @@ sources = [ 'c/wt.c', ] -# pyx -> c transpile output depends on copied __init__.py and pxd files +# No generated Cython sources are involved anymore, but these copies are +# still needed for correct rebuilds: the `depends: _cython_tree` on the +# generator is what re-cythonizes the .pyx files when a .pxd changes. +# Cython's own depfile output (`cython -M`, used by Meson's native Cython +# support) omits module-namespace cimports like `from . cimport c_wt` +# (seen with Cython 3.2.9), so switching to it would cause stale builds. _cython_tree = [ fs.copyfile('__init__.py'), fs.copyfile('common.pxd'), @@ -69,11 +67,7 @@ _cython_tree = [ fs.copyfile('wavelets_list.pxi'), ] -# 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 libc_wt = static_library('c_wt', sources, @@ -101,7 +95,7 @@ pyx_files = [ ['_pywt', fs.copyfile('_pywt.pyx')], ] -cy_deps = declare_dependency(sources: [__init__py, _cython_tree])#, config_pxi]) +cy_deps = declare_dependency(sources: [__init__py, _cython_tree]) foreach pyx_file: pyx_files py.extension_module(pyx_file[0], [cython_gen.process(pyx_file[1])], diff --git a/pywt/_multidim.py b/pywt/_multidim.py index 5a30b80c..361c2cf2 100644 --- a/pywt/_multidim.py +++ b/pywt/_multidim.py @@ -12,7 +12,6 @@ import numpy as np -from ._c99_config import _have_c99_complex from ._extensions._dwt import dwt_axis, idwt_axis from ._utils import AxisError, _modes_per_axis, _wavelets_per_axis @@ -162,10 +161,6 @@ def dwtn(data, wavelet, mode='symmetric', axes=None): """ data = np.asarray(data) - if not _have_c99_complex and np.iscomplexobj(data): - real = dwtn(data.real, wavelet, mode, axes) - imag = dwtn(data.imag, wavelet, mode, axes) - return {k: real[k] + 1j * imag[k] for k in real} if data.dtype == np.dtype('object'): raise TypeError("Input must be a numeric array-like") @@ -256,12 +251,6 @@ def idwtn(coeffs, wavelet, mode='symmetric', axes=None): # Raise error for invalid key combinations coeffs = _fix_coeffs(coeffs) - if (not _have_c99_complex and - any(np.iscomplexobj(v) for v in coeffs.values())): - real_coeffs = {k: v.real for k, v in coeffs.items()} - imag_coeffs = {k: v.imag for k, v in coeffs.items()} - return (idwtn(real_coeffs, wavelet, mode, axes) + - 1j * idwtn(imag_coeffs, wavelet, mode, axes)) # key length matches the number of axes transformed ndim_transform = max(len(key) for key in coeffs) diff --git a/pywt/_pytesttester.py b/pywt/_pytesttester.py index 6344fa8e..298ea3a0 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: diff --git a/pywt/_swt.py b/pywt/_swt.py index 8b279d4e..5588b8f7 100644 --- a/pywt/_swt.py +++ b/pywt/_swt.py @@ -3,7 +3,6 @@ import numpy as np -from ._c99_config import _have_c99_complex from ._extensions._dwt import idwt_single from ._extensions._pywt import Modes, Wavelet, _check_dtype from ._extensions._swt import swt as _swt @@ -114,20 +113,6 @@ def swt(data, wavelet, level=None, start_level=0, axis=-1, Cambridge University Press, 2000. """ - if not _have_c99_complex and np.iscomplexobj(data): - data = np.asarray(data) - kwargs = {"wavelet": wavelet, "level": level, "start_level": start_level, - "trim_approx": trim_approx, "axis": axis, "norm": norm} - coeffs_real = swt(data.real, **kwargs) - coeffs_imag = swt(data.imag, **kwargs) - if not trim_approx: - coeffs_cplx = [] - for (cA_r, cD_r), (cA_i, cD_i) in zip(coeffs_real, coeffs_imag): - coeffs_cplx.append((cA_r + 1j*cA_i, cD_r + 1j*cD_i)) - else: - coeffs_cplx = [cr + 1j*ci - for (cr, ci) in zip(coeffs_real, coeffs_imag)] - return coeffs_cplx # accept array_like input; make a copy to ensure a contiguous array dt = _check_dtype(data) @@ -200,16 +185,6 @@ def iswt(coeffs, wavelet, norm=False, axis=-1): return iswtn(coeffs_nd, wavelet, axes=(axis,), norm=norm) elif axis != 0 and axis != -1: raise AxisError("Axis greater than data dimensions") - if not _have_c99_complex and np.iscomplexobj(cA): - if trim_approx: - coeffs_real = [c.real for c in coeffs] - coeffs_imag = [c.imag for c in coeffs] - else: - coeffs_real = [(ca.real, cd.real) for ca, cd in coeffs] - coeffs_imag = [(ca.imag, cd.imag) for ca, cd in coeffs] - kwargs = {"wavelet": wavelet, "norm": norm} - y = iswt(coeffs_real, **kwargs) - return y + 1j * iswt(coeffs_imag, **kwargs) if trim_approx: coeffs = coeffs[1:] @@ -440,20 +415,6 @@ def iswt2(coeffs, wavelet, norm=False, axes=(-2, -1)): coeffs_nd = [{'aa': a, 'da': h, 'ad': v, 'dd': d} for a, (h, v, d) in coeffs] return iswtn(coeffs_nd, wavelet, axes=axes, norm=norm) - if not _have_c99_complex and np.iscomplexobj(cA): - if trim_approx: - coeffs_real = [cA.real] - coeffs_real += [(h.real, v.real, d.real) for h, v, d in coeffs[1:]] - coeffs_imag = [cA.imag] - coeffs_imag += [(h.imag, v.imag, d.imag) for h, v, d in coeffs[1:]] - else: - coeffs_real = [(a.real, (h.real, v.real, d.real)) - for a, (h, v, d) in coeffs] - coeffs_imag = [(a.imag, (h.imag, v.imag, d.imag)) - for a, (h, v, d) in coeffs] - kwargs = {"wavelet": wavelet, "norm": norm} - y = iswt2(coeffs_real, **kwargs) - return y + 1j * iswt2(coeffs_imag, **kwargs) if trim_approx: coeffs = coeffs[1:] @@ -613,21 +574,6 @@ def swtn(data, wavelet, level, start_level=0, axes=None, trim_approx=False, """ data = np.asarray(data) - if not _have_c99_complex and np.iscomplexobj(data): - kwargs = {"wavelet": wavelet, "level": level, "start_level": start_level, - "trim_approx": trim_approx, "axes": axes, "norm": norm} - real = swtn(data.real, **kwargs) - imag = swtn(data.imag, **kwargs) - if trim_approx: - cplx = [real[0] + 1j * imag[0]] - offset = 1 - else: - cplx = [] - offset = 0 - for rdict, idict in zip(real[offset:], imag[offset:]): - cplx.append( - {k: rdict[k] + 1j * idict[k] for k in rdict}) - return cplx if data.dtype == np.dtype('object'): raise TypeError("Input must be a numeric array-like") @@ -719,19 +665,6 @@ def iswtn(coeffs, wavelet, axes=None, norm=False): trim_approx = not isinstance(coeffs[0], dict) cA = coeffs[0] if trim_approx else coeffs[0]['a'*ndim_transform] - if not _have_c99_complex and np.iscomplexobj(cA): - if trim_approx: - coeffs_real = [coeffs[0].real] - coeffs_imag = [coeffs[0].imag] - coeffs = coeffs[1:] - else: - coeffs_real = [] - coeffs_imag = [] - coeffs_real += [{k: v.real for k, v in c.items()} for c in coeffs] - coeffs_imag += [{k: v.imag for k, v in c.items()} for c in coeffs] - kwargs = {"wavelet": wavelet, "axes": axes, "norm": norm} - y = iswtn(coeffs_real, **kwargs) - return y + 1j * iswtn(coeffs_imag, **kwargs) if trim_approx: coeffs = coeffs[1:] diff --git a/pywt/meson.build b/pywt/meson.build index 8c66156f..57e9ecd5 100644 --- a/pywt/meson.build +++ b/pywt/meson.build @@ -35,26 +35,32 @@ generate_version = custom_target( cc = meson.get_compiler('c') -# Check for C99 complex types support. MSVC shouldn't have it, other compilers -# most likely will. -cdata = configuration_data() -have_c99_complex = cc.has_type('complex float', prefix: '#include ') -if have_c99_complex - cdata.set('have_c99_complex', 'True') - cdata.set('have_c99_complex_cy', '1') -else - cdata.set('have_c99_complex', 'False') - cdata.set('have_c99_complex_cy', '0') +# Check for complex support (same check as numpy/_core/meson.build; MSVC has +# no C99 complex and provides struct-based types instead, see +# pywt/_extensions/c/common.h) +if not cc.has_header('complex.h') + error('"complex.h" header not found') endif -_c99_config = configure_file( - input: '_c99_config.py.in', - output: '_c99_config.py', - configuration: cdata, - install: true, - install_dir: py.get_install_dir() / 'pywt', - install_tag: 'python-runtime', -) +if cc.get_id() == 'msvc' or cc.get_id() == 'clang-cl' + complex_types_to_check = [ + ['PYWT_SIZEOF_COMPLEX_FLOAT', '_Fcomplex'], + ['PYWT_SIZEOF_COMPLEX_DOUBLE', '_Dcomplex'], + ] +else + complex_types_to_check = [ + ['PYWT_SIZEOF_COMPLEX_FLOAT', 'complex float'], + ['PYWT_SIZEOF_COMPLEX_DOUBLE', 'complex double'], + ] +endif +cdata = configuration_data() +foreach symbol_type: complex_types_to_check + if not cc.has_type(symbol_type[1], prefix: '#include ') + t = symbol_type[1] + error(f'"complex.h" header does not include complex type @t@') + endif + cdata.set(symbol_type[0], cc.sizeof(symbol_type[1], prefix: '#include ')) +endforeach install_subdir('data', install_dir: py.get_install_dir() / 'pywt') install_subdir('tests', install_dir: py.get_install_dir() / 'pywt')