From 3df8b5563d5b06ec37840d19e8ba0661c15c4e33 Mon Sep 17 00:00:00 2001 From: Saransh Singh Date: Thu, 28 May 2026 14:00:39 -0700 Subject: [PATCH 1/4] better implementation of the numba accelerated exp1 and wofz scipy.special functions Also remove unused constants and include more tests for wofz and exp1 functions. We also must enable numba JIT coverage instead of disabling numba JIT now. --- .github/workflows/test.yml | 2 +- hexrd/core/constants.py | 43 --- hexrd/core/fitting/peakfunctions.py | 91 +---- hexrd/core/fitting/special.py | 444 +++++++++++++++++++++++ hexrd/core/fitting/utils.py | 96 +---- hexrd/powder/wppf/peakfunctions.py | 12 +- tests/core/fitting/test_fitting_utils.py | 9 - tests/core/fitting/test_peakfunctions.py | 111 +++--- 8 files changed, 510 insertions(+), 298 deletions(-) create mode 100644 hexrd/core/fitting/special.py diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 438f8fe7f..011359a9c 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -72,7 +72,7 @@ jobs: HEXRD_EXAMPLE_REPO_PATH: ${{ github.workspace }}/examples NUMBA_DISABLE_JIT: 0 run: | - NUMBA_DISABLE_JIT=1 coverage run --source hexrd -m pytest tests/ + NUMBA_JIT_COVERAGE=1 coverage run --source hexrd -m pytest tests/ coverage combine coverage xml -i if: ${{ matrix.config.os == 'ubuntu-latest'}} diff --git a/hexrd/core/constants.py b/hexrd/core/constants.py index 126690451..a7eac458a 100644 --- a/hexrd/core/constants.py +++ b/hexrd/core/constants.py @@ -138,49 +138,6 @@ ] ).astype(np.float64) -c_coeff_exp1exp = np.array( - [ - 0.999999584, - -0.249992399, - 0.055514994, - -0.010315766, - 0.001535370, - -0.000142164, - ] -).astype(np.complex128) - -cnum_exp1exp = np.array( - [ - 1.0, - 99.0, - 3952.0, - 82544.0, - 979524.0, - 6712860.0, - 25815840.0, - 51369120.0, - 44339040.0, - 10628640.0, - 0.0, - ] -).astype(np.complex128) - -cden_exp1exp = np.array( - [ - 1.0, - 100.0, - 4050.0, - 86400.0, - 1058400.0, - 7620480.0, - 31752000.0, - 72576000.0, - 81648000.0, - 36288000.0, - 3628800.0, - ] -).astype(np.complex128) - """ >> @AUTHOR: Saransh Singh, Lawrence Livermore National Lab, diff --git a/hexrd/core/fitting/peakfunctions.py b/hexrd/core/fitting/peakfunctions.py index a70789507..ca2f131a0 100644 --- a/hexrd/core/fitting/peakfunctions.py +++ b/hexrd/core/fitting/peakfunctions.py @@ -29,12 +29,9 @@ from numba import njit import copy from hexrd.core import constants -from hexrd.core.constants import ( - c_erf, - cnum_exp1exp, - cden_exp1exp, - c_coeff_exp1exp, -) +from hexrd.core.fitting.special import erfc, exp1_complex_numba, wofz + +exp1exp = exp1_complex_numba gauss_width_fact = constants.sigma_to_fwhm lorentz_width_fact = 2.0 @@ -49,88 +46,6 @@ 'pink_beam_dcs': 8, } -""" -cutom function to compute the complementary error function -based on rational approximation of the convergent Taylor -series. coefficients found in -Formula 7.1.26 -Handbook of Mathematical Functions, -Abramowitz and Stegun -Error is < 1.5e-7 for all x -""" - - -@njit(cache=True, nogil=True) -def erfc(x): - # save the sign of x - sign = np.sign(x) - x = np.abs(x) - - # constants - a1, a2, a3, a4, a5, p = c_erf - - # A&S formula 7.1.26 - t = 1.0 / (1.0 + p * x) - y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * np.exp(-x * x) - erf = sign * y # erf(-x) = -erf(x) - return 1.0 - erf - - -""" -cutom function to compute the exponential integral -based on Padé approximation of exponential integral -function. coefficients found in pg. 231 Abramowitz -and Stegun, eq. 5.1.53 -""" - - -@njit(cache=True, nogil=True) -def exp1exp_under1(x): - f = np.zeros(x.shape).astype(np.complex128) - for i in range(6): - xx = x ** (i + 1) - f += c_coeff_exp1exp[i] * xx - - return (f - np.log(x) - np.euler_gamma) * np.exp(x) - - -""" -cutom function to compute the exponential integral -based on Padé approximation of exponential integral -function. coefficients found in pg. 415 Y. Luke, The -special functions and their approximations, vol 2 -(1969) Elsevier -""" - - -@njit(cache=True, nogil=True) -def exp1exp_over1(x): - num = np.zeros(x.shape).astype(np.complex128) - den = np.zeros(x.shape).astype(np.complex128) - - for i in range(11): - p = 10 - i - if p != 0: - xx = x**p - num += cnum_exp1exp[i] * xx - den += cden_exp1exp[i] * xx - else: - num += cnum_exp1exp[i] - den += cden_exp1exp[i] - - return (num / den) * (1.0 / x) - - -@njit(cache=True, nogil=True) -def exp1exp(x): - mask = np.sign(x.real) * np.abs(x) > 1.0 - - f = np.zeros(x.shape).astype(np.complex128) - f[mask] = exp1exp_over1(x[mask]) - f[~mask] = exp1exp_under1(x[~mask]) - - return f - # ============================================================================= # 1-D Gaussian Functions diff --git a/hexrd/core/fitting/special.py b/hexrd/core/fitting/special.py new file mode 100644 index 000000000..21d28f03c --- /dev/null +++ b/hexrd/core/fitting/special.py @@ -0,0 +1,444 @@ +import math +import numpy as np +from numba import njit, vectorize, complex128, float64 +from hexrd.core.constants import c_erf + +_EULER_GAMMA = 0.577215664901532860606512090082402431 +_EPS = 2.220446049250313e-16 +_TINY = 1.0e-300 +_MAX_ITER = 1000 + +_SMALL_SERIES_CUTOFF = 4.0 +_LEFT_SERIES_CUTOFF = 100.0 +_MAX_LOG = 709.782712893384 # approximately log(max float64) + + +# Fixed N=11 modified trapezoidal-rule approximation. +# h = sqrt(pi / (N + 1)), H = pi / h +_NTRAP = 11 +_H_STEP = math.sqrt(math.pi / (_NTRAP + 1.0)) +_H_CAP = math.pi / _H_STEP + +_MAX_EXP = 709.782712893384 # approx log(float64 max) +_MIN_EXP = -745.1332191019411 # exp(x) underflows to 0 below this + + +@njit(cache=True, nopython=True) +def erfc(x): + # save the sign of x + sign = np.sign(x) + x = np.abs(x) + + # constants + a1, a2, a3, a4, a5, p = c_erf + + # A&S formula 7.1.26 + t = 1.0 / (1.0 + p * x) + y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * np.exp(-x * x) + erf = sign * y # erf(-x) = -erf(x) + return 1.0 - erf + + +@njit(cache=True, nopython=True) +def _complex_log_principal_numba(z): + """ + Principal complex log, written manually so signed zero in imag + is handled through atan2. + """ + return complex( + math.log(math.hypot(z.real, z.imag)), + math.atan2(z.imag, z.real), + ) + + +@njit(cache=True, nopython=True) +def _exp_neg_complex_numba(z): + """ + Compute exp(-z) with simple overflow/underflow guards. + """ + x = z.real + y = z.imag + + if x > 745.0: + return complex(0.0, 0.0) + + if x < -_MAX_LOG: + # The true result is generally outside float64 range. + return complex(math.inf, math.nan) + + scale = math.exp(-x) + return complex(scale * math.cos(y), -scale * math.sin(y)) + + +@njit(cache=True, nopython=True) +def _exp1_series_numba(z): + """ + Logarithmic power series: + + E1(z) = -gamma - log(z) - sum_{k=1}^inf (-z)^k / (k k!) + + Valid on the principal branch, away from z = 0. + """ + term = -z + s = term + + for k in range(2, _MAX_ITER + 1): + term *= -z / k + add = term / k + s += add + + if abs(add) <= _EPS * max(1.0, abs(s)): + break + + return complex(-_EULER_GAMMA, 0.0) - _complex_log_principal_numba(z) - s + + +@njit(cache=True, nopython=True) +def _exp1_contfrac_numba(z): + """ + Continued fraction via modified Lentz iteration. + Good away from the origin and away from the difficult part + of the negative-real branch cut. + """ + b = z + 1.0 + c = complex(1.0 / _TINY, 0.0) + d = 1.0 / b + h = d + + for i in range(1, _MAX_ITER + 1): + a = -float(i * i) + b += 2.0 + + denom = a * d + b + if abs(denom) < _TINY: + denom = complex(_TINY, 0.0) + d = 1.0 / denom + + c = b + a / c + if abs(c) < _TINY: + c = complex(_TINY, 0.0) + + delta = c * d + h *= delta + + if abs(delta - 1.0) <= _EPS: + break + + return h * _exp_neg_complex_numba(z) + + +@njit(cache=True, nopython=True) +def exp1_complex_scalar_numba(z): + """ + Numba-compatible complex exponential integral E1(z). + + This is the complex analogue of scipy.special.exp1 for complex128 input. + For exact negative-real complex inputs, the side of the branch cut is + selected by the sign of the imaginary zero: + complex(-x, +0.0) -> -i*pi side + complex(-x, -0.0) -> +i*pi side + """ + x = z.real + y = z.imag + + if math.isnan(x) or math.isnan(y): + return complex(math.nan, math.nan) + + if x == 0.0 and y == 0.0: + return complex(math.inf, 0.0) + + if math.isinf(x) or math.isinf(y): + if x == math.inf and not math.isinf(y): + return complex(0.0, 0.0) + return complex(math.nan, math.nan) + + az = abs(z) + + near_negative_real_axis = ( + x < 0.0 and az <= _LEFT_SERIES_CUTOFF and abs(y) <= max(2.0, 0.25 * abs(x)) + ) + + if az <= _SMALL_SERIES_CUTOFF or near_negative_real_axis: + return _exp1_series_numba(z) + + out = _exp1_contfrac_numba(z) + + # If the continued fraction path is used exactly on the negative-real + # branch cut, insert the correct signed-zero branch jump. + if x < 0.0 and y == 0.0: + if math.copysign(1.0, y) < 0.0: + out += complex(0.0, math.pi) + else: + out += complex(0.0, -math.pi) + + return out + + +@vectorize([complex128(complex128)], nopython=True) +def exp1_complex_numba(z): + """ + Vectorized complex ufunc version. + """ + return exp1_complex_scalar_numba(z) + + +@njit(cache=True, nopython=True) +def exp1_real_scalar_numba(x): + """ + Real-valued SciPy-like wrapper. + + For real negative x, scipy.special.exp1(x) returns nan. + For complex negative x + 0j, use exp1_complex_scalar_numba instead. + """ + if math.isnan(x): + return math.nan + + if x < 0.0: + return math.nan + + if x == 0.0: + return math.inf + + if math.isinf(x): + return 0.0 + + return exp1_complex_scalar_numba(complex(x, 0.0)).real + + +@vectorize([float64(float64)], nopython=True) +def exp1_real_numba(x): + """ + Vectorized real ufunc version. + """ + return exp1_real_scalar_numba(x) + + +# def exp1exp(x): +# if x.dtype == float: +# return exp1_real_numba(x) +# elif x.dtype == complex: +# return exp1_complex_numba(x) + + +@njit(cache=True, nopython=True) +def _signed_inf_from_trig(v): + if v > 0.0: + return math.inf + if v < 0.0: + return -math.inf + return 0.0 + + +@njit(cache=True, nopython=True) +def _cexp_parts(a, b): + """ + Compute exp(a + i b) with simple overflow/underflow guards. + """ + if a > _MAX_EXP: + c = math.cos(b) + s = math.sin(b) + return complex(_signed_inf_from_trig(c), _signed_inf_from_trig(s)) + + if a < _MIN_EXP: + return complex(0.0, 0.0) + + ea = math.exp(a) + return complex(ea * math.cos(b), ea * math.sin(b)) + + +@njit(cache=True, nopython=True) +def _exp_neg_z2(z): + """ + Compute exp(-z*z), where z = x + i y. + -z^2 = (y^2 - x^2) - 2 i x y + """ + x = z.real + y = z.imag + return _cexp_parts(y * y - x * x, -2.0 * x * y) + + +@njit(cache=True, nopython=True) +def _trap_correction(z, plus): + """ + Correction term: + + 2 exp(-z^2) / (1 + exp(-2 i H z)) if plus=True + 2 exp(-z^2) / (1 - exp(-2 i H z)) if plus=False + + Uses a reciprocal form when exp(-2 i H z) would be very large. + """ + x = z.real + y = z.imag + + log_abs_b = 2.0 * _H_CAP * y + + if log_abs_b > 40.0: + # B = exp(-2 i H z), inv_b = 1/B = exp(2 i H z) + inv_b = _cexp_parts( + -2.0 * _H_CAP * y, + 2.0 * _H_CAP * x, + ) + + # 2 exp(-z^2) / B = 2 exp(-z^2 + 2 i H z) + a_over_b = 2.0 * _cexp_parts( + y * y - x * x - 2.0 * _H_CAP * y, + -2.0 * x * y + 2.0 * _H_CAP * x, + ) + + if plus: + return a_over_b / (1.0 + inv_b) + else: + return a_over_b / (inv_b - 1.0) + + else: + a = 2.0 * _exp_neg_z2(z) + b = _cexp_parts( + log_abs_b, + -2.0 * _H_CAP * x, + ) + + if plus: + return a / (1.0 + b) + else: + return a / (1.0 - b) + + +@njit(cache=True, nopython=True) +def _wofz_midpoint(z): + """ + Midpoint-rule approximation used in part of the first quadrant. + """ + z2 = z * z + az = complex(0.0, 2.0 / _H_CAP) * z + + h0 = 0.5 * _H_STEP + t2 = h0 * h0 + s = math.exp(-t2) / (z2 - t2) + + for k in range(1, _NTRAP + 1): + t = _H_STEP * (k + 0.5) + t2 = t * t + s += math.exp(-t2) / (z2 - t2) + + return az * s + + +@njit(cache=True, nogil=True) +def _wofz_modified_midpoint(z): + """ + Modified midpoint-rule approximation. + """ + return _wofz_midpoint(z) + _trap_correction(z, True) + + +@njit(cache=True, nopython=True) +def _wofz_modified_trapezoid(z): + """ + Modified trapezoidal-rule approximation. + """ + z2 = z * z + az = complex(0.0, 2.0 / _H_CAP) * z + + t = _H_STEP * _NTRAP + t2 = t * t + s = math.exp(-t2) / (z2 - t2) + + for k in range(1, _NTRAP): + t = _H_STEP * k + t2 = t * t + s += math.exp(-t2) / (z2 - t2) + + return complex(0.0, 1.0 / _H_CAP) / z + az * s + _trap_correction(z, False) + + +@njit(cache=True, nopython=True) +def _wofz_first_quadrant(z): + """ + Evaluate w(z) for Re(z) >= 0, Im(z) >= 0. + """ + x = z.real + y = z.imag + + rzh = x / _H_STEP + buff = abs(rzh - math.floor(rzh) - 0.5) + + if y >= max(x, _H_CAP): + return _wofz_midpoint(z) + + elif y < x and buff <= 0.25: + return _wofz_modified_trapezoid(z) + + else: + return _wofz_modified_midpoint(z) + + +@njit(cache=True, nopython=True) +def wofz_complex_scalar_numba(z): + """ + Numba-compatible scalar complex Faddeeva function. + + Approximate replacement for: + + scipy.special.wofz(z) + + for complex128-like finite inputs. + """ + x = z.real + y = z.imag + + if math.isnan(x) or math.isnan(y): + return complex(math.nan, math.nan) + + # Basic SciPy-like limiting behavior for common infinite cases. + # For ordinary finite inputs, the main algorithm below is used. + if math.isinf(x) or math.isinf(y): + if y == -math.inf and x == 0.0: + return complex(math.inf, 0.0) + if y == -math.inf: + return complex(math.nan, math.nan) + return complex(0.0, 0.0) + + # Map arbitrary z to the first quadrant using symmetries. + xneg = x < 0.0 + yneg = y < 0.0 + not_both = xneg != yneg + + zw = z + + if xneg: + zw = -zw + + if not_both: + zw = complex(zw.real, -zw.imag) + + w = _wofz_first_quadrant(zw) + + if not_both: + w = complex(w.real, -w.imag) + + if yneg: + # w(z) = 2 exp(-z^2) - w(-z) + w = 2.0 * _exp_neg_z2(z) - w + + return w + + +@njit(cache=True, nopython=True) +def wofz_real_scalar_numba(x): + """ + Convenience wrapper for real x, returning complex w(x). + """ + return wofz_complex_scalar_numba(complex(x, 0.0)) + + +@vectorize( + [complex128(complex128), complex128(float64)], + nopython=True, +) +def wofz(z): + """ + Vectorized ufunc version. + + Works on complex128 arrays and float64 arrays. + """ + return wofz_complex_scalar_numba(z) diff --git a/hexrd/core/fitting/utils.py b/hexrd/core/fitting/utils.py index 5d0ce8b9f..ed93fe1fd 100644 --- a/hexrd/core/fitting/utils.py +++ b/hexrd/core/fitting/utils.py @@ -2,15 +2,12 @@ import numpy as np from numba import njit +from hexrd.core.fitting.special import erfc, exp1_complex_numba -from hexrd.core.constants import ( - c_erf, - cnum_exp1exp, - cden_exp1exp, - c_coeff_exp1exp, -) from hexrd.core.matrixutil import uniqueVectors +# alias for exp1 function +exp1exp = exp1_complex_numba # ============================================================================= # LMFIT Parameter munging utilities @@ -125,93 +122,6 @@ def _set_peak_center_bounds(params, window_range, min_sep=0.01): prev_peak = curr_peak -# ============================================================================= -# DCS-related utilities -# ============================================================================= - -""" -cutom function to compute the complementary error function -based on rational approximation of the convergent Taylor -series. coefficients found in -Formula 7.1.26 -Handbook of Mathematical Functions, -Abramowitz and Stegun -Error is < 1.5e-7 for all x -""" - - -@njit(cache=True, nogil=True) -def erfc(x): - # save the sign of x - sign = np.sign(x) - x = np.abs(x) - - # constants - a1, a2, a3, a4, a5, p = c_erf - - # A&S formula 7.1.26 - t = 1.0 / (1.0 + p * x) - y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * np.exp(-x * x) - erf = sign * y # erf(-x) = -erf(x) - return 1.0 - erf - - -""" -cutom function to compute the exponential integral -based on Padé approximation of exponential integral -function. coefficients found in pg. 231 Abramowitz -and Stegun, eq. 5.1.53 -""" - - -@njit(cache=True, nogil=True) -def exp1exp_under1(x): - f = np.zeros(x.shape).astype(np.complex128) - for i in range(6): - xx = x ** (i + 1) - f += c_coeff_exp1exp[i] * xx - - return (f - np.log(x) - np.euler_gamma) * np.exp(x) - - -""" -cutom function to compute the exponential integral -based on Padé approximation of exponential integral -function. coefficients found in pg. 415 Y. Luke, The -special functions and their approximations, vol 2 -(1969) Elsevier -""" - - -@njit(cache=True, nogil=True) -def exp1exp_over1(x): - num = np.zeros(x.shape).astype(np.complex128) - den = np.zeros(x.shape).astype(np.complex128) - - for i in range(11): - p = 10 - i - if p != 0: - xx = x**p - num += cnum_exp1exp[i] * xx - den += cden_exp1exp[i] * xx - else: - num += cnum_exp1exp[i] - den += cden_exp1exp[i] - - return (num / den) * (1.0 / x) - - -@njit(cache=True, nogil=True) -def exp1exp(x): - mask = np.sign(x.real) * np.abs(x) > 1.0 - - f = np.zeros(x.shape).astype(np.complex128) - f[mask] = exp1exp_over1(x[mask]) - f[~mask] = exp1exp_under1(x[~mask]) - - return f - - @njit(cache=True, nogil=True) def _calc_alpha(alpha, x0): a0, a1 = alpha diff --git a/hexrd/powder/wppf/peakfunctions.py b/hexrd/powder/wppf/peakfunctions.py index fcb283ad3..8c14faa34 100644 --- a/hexrd/powder/wppf/peakfunctions.py +++ b/hexrd/powder/wppf/peakfunctions.py @@ -29,8 +29,9 @@ import copy from hexrd.core import constants from numba import vectorize, float64, njit, prange -from hexrd.core.fitting.peakfunctions import erfc, exp1exp +from hexrd.core.fitting.special import erfc, exp1_complex_numba +exp1exp = exp1_complex_numba # from scipy.special import erfc, exp1 # addr = get_cython_function_address("scipy.special.cython_special", "exp1") @@ -521,7 +522,7 @@ def _gaussian_pink_beam(alpha, beta, fwhm_g, tth, tth_list): t1 = erfc(y) t2 = erfc(z) g = np.zeros(tth_list.shape) - zmask = np.abs(del_tth) > 5.0 + zmask = np.abs(del_tth) > 100.0 g[~zmask] = (0.5 * (alpha * beta) / (alpha + beta)) * np.exp(u[~zmask]) * t1[ ~zmask ] + np.exp(v[~zmask]) * t2[~zmask] @@ -550,7 +551,12 @@ def _lorentzian_pink_beam(alpha, beta, fwhm_l, tth, tth_list): # f1 = exp1(p) # f2 = exp1(q) - y = -(alpha * beta) / (np.pi * (alpha + beta)) * (f1 + f2).imag + y = np.zeros(tth_list.shape) + zmask = np.abs(del_tth) > 100.0 + + y[~zmask] = ( + -(alpha * beta) / (np.pi * (alpha + beta)) * (f1[~zmask] + f2[~zmask]).imag + ) mask = np.isnan(y) y[mask] = 0.0 diff --git a/tests/core/fitting/test_fitting_utils.py b/tests/core/fitting/test_fitting_utils.py index d168b1540..ea28c6634 100644 --- a/tests/core/fitting/test_fitting_utils.py +++ b/tests/core/fitting/test_fitting_utils.py @@ -83,15 +83,6 @@ def test_peak_center_bounds(mock_uniq, params): def test_math_funcs(): - assert np.isclose(utils.erfc(np.array([0.0, 100.0]))[0], 1.0) - assert np.isclose(utils.erfc(np.array([-100.0]))[0], 2.0) - - assert not np.isnan(utils.exp1exp_under1(np.array([0.1 + 0j]))[0]) - assert not np.isnan(utils.exp1exp_over1(np.array([2.0 + 0j]))[0]) - assert not np.any( - np.isnan(utils.exp1exp(np.array([0.1, 2.0], dtype=np.complex128))) - ) - assert utils._calc_alpha(np.array([1.0, 2.0]), 0.0) == 1.0 assert utils._calc_beta(np.array([3.0, 4.0]), 0.0) == 3.0 diff --git a/tests/core/fitting/test_peakfunctions.py b/tests/core/fitting/test_peakfunctions.py index c44bd7961..39e41e44f 100644 --- a/tests/core/fitting/test_peakfunctions.py +++ b/tests/core/fitting/test_peakfunctions.py @@ -1,8 +1,9 @@ # tests/test_peakfunctions_simplified.py import numpy as np import pytest -from scipy.special import erfc as scipy_erfc +from scipy.special import erfc as scipy_erfc, wofz as scipy_wofz, exp1 as scipy_exp1 import hexrd.core.fitting.peakfunctions as pf +import hexrd.core.fitting.special as sp # ----------------- small utilities ----------------- @@ -21,18 +22,41 @@ def finite_diff_derivative(func, p, x, idx, rel_step=1e-6): # ----------------- special functions ----------------- -def test_erfc_and_exp1exp_smoke(): - vals = np.array([-5.0, -1.0, 0.0, 0.5, 1.0, 3.2], dtype=float) - got = pf.erfc(vals) +def test_erfc(): + # test real values + vals = np.random.uniform(-5, 5, 5) + got = sp.erfc(vals) assert got.shape == vals.shape assert np.all((got >= 0.0) & (got <= 2.0)) assert np.allclose(got, scipy_erfc(vals), rtol=1e-6, atol=1e-6) - x_pos = np.array([0.1, 0.5, 1.0, 2.0], dtype=float) - for fn in (pf.exp1exp_under1, pf.exp1exp_over1, pf.exp1exp): - arr = fn(x_pos) - assert arr.shape == x_pos.shape - assert np.all(np.isfinite(arr)) + +def test_wofz(): + # test complex values + vals = np.random.uniform(-5, 5, 5) + 1j * np.random.uniform(-5, 5, 5) + got = sp.wofz(vals) + assert got.shape == vals.shape + assert np.allclose(got, scipy_wofz(vals), rtol=1e-6, atol=1e-6) + + +def test_exp1(): + # test real values + vals = np.random.uniform(-5, 5, 5) + got = sp.exp1_real_numba(vals) + assert got.shape == vals.shape + # filter out nans. these are infs in scipy + mask = ~np.isnan(got) + assert np.allclose(got[mask], scipy_exp1(vals)[mask], rtol=1e-6, atol=1e-6) + + # test complex values + vals = np.random.uniform(-5, 5, 5) + 1j * np.random.uniform(-5, 5, 5) + got = sp.exp1_complex_numba(vals) + scipy_got = scipy_exp1(vals) + assert got.shape == vals.shape + + # filter out nans. these are infs in scipy + mask = ~np.isnan(got) + assert np.allclose(got[mask], scipy_got[mask], rtol=1e-6, atol=1e-6) # ----------------- 1D Gaussian ----------------- @@ -80,12 +104,8 @@ def test_gaussian1d_full_and_derivatives(): np.testing.assert_allclose(dmat[0:3, :], analytic_no_bg) for i in range(3): - fd = finite_diff_derivative( - pf._gaussian1d_no_bg, p[:3], x, i, rel_step=1e-6 - ) - np.testing.assert_allclose( - analytic_no_bg[i, :], fd, rtol=1e-5, atol=1e-8 - ) + fd = finite_diff_derivative(pf._gaussian1d_no_bg, p[:3], x, i, rel_step=1e-6) + np.testing.assert_allclose(analytic_no_bg[i, :], fd, rtol=1e-5, atol=1e-8) # ----------------- Lorentzian ----------------- @@ -113,9 +133,7 @@ def test_lorentzian_unit_and_full_and_deriv(): np.testing.assert_allclose(no_bg_out, A * expected_unit, rtol=1e-12) full_out = pf.lorentzian1d(np.array([A, x0, FWHM, c0, c1]), x) - np.testing.assert_allclose( - full_out, A * expected_unit + c0 + c1 * x, rtol=1e-12 - ) + np.testing.assert_allclose(full_out, A * expected_unit + c0 + c1 * x, rtol=1e-12) assert int(np.argmax(no_bg_out)) == int(np.argmin(np.abs(x - x0))) p = np.array([2.5, -0.3, 0.9], dtype=float) @@ -191,25 +209,19 @@ def test_pvoigt_and_split_variants(): FWHM_plus = 1.6 n_minus = 0.2 n_plus = 0.7 - p_no_bg = np.array( - [A, x0, FWHM_minus, FWHM_plus, n_minus, n_plus], dtype=float - ) + p_no_bg = np.array([A, x0, FWHM_minus, FWHM_plus, n_minus, n_plus], dtype=float) x_small = np.linspace(-2.0, 2.0, 401) out_no_bg = pf._split_pvoigt1d_no_bg(p_no_bg, x_small) expected = np.zeros_like(x_small) xr = x_small >= x0 xl = x_small < x0 - expected[xr] = A * pf._unit_pvoigt1d( - np.array([x0, FWHM_plus, n_plus]), x_small[xr] - ) + expected[xr] = A * pf._unit_pvoigt1d(np.array([x0, FWHM_plus, n_plus]), x_small[xr]) expected[xl] = A * pf._unit_pvoigt1d( np.array([x0, FWHM_minus, n_minus]), x_small[xl] ) np.testing.assert_allclose(out_no_bg, expected) np.testing.assert_allclose( - pf.split_pvoigt1d( - np.hstack([p_no_bg, np.array([0.05, -0.01])]), x_small - ), + pf.split_pvoigt1d(np.hstack([p_no_bg, np.array([0.05, -0.01])]), x_small), expected + 0.05 - 0.01 * x_small, ) @@ -217,8 +229,7 @@ def test_pvoigt_and_split_variants(): x_sym = np.linspace(-2.0, 2.0, 201) np.testing.assert_allclose( pf._split_pvoigt1d_no_bg(p_sym, x_sym), - p_sym[0] - * pf._unit_pvoigt1d(np.array([p_sym[1], p_sym[2], p_sym[4]]), x_sym), + p_sym[0] * pf._unit_pvoigt1d(np.array([p_sym[1], p_sym[2], p_sym[4]]), x_sym), ) p_zeroA = np.array( @@ -245,9 +256,7 @@ def test_calc_alpha_beta_and_mixing_factor_behavior_and_wrappers(): assert np.allclose(got_beta, expected_beta, atol=1e-12) for fwhm_g, fwhm_l in [(0.1, 0.2), (0.5, 0.5), (1.0, 2.0), (2.5, 0.3)]: - eta, fwhm = pf._mixing_factor_pv( - np.float64(fwhm_g), np.float64(fwhm_l) - ) + eta, fwhm = pf._mixing_factor_pv(np.float64(fwhm_g), np.float64(fwhm_l)) assert np.isfinite(eta) and np.isfinite(fwhm) assert 0.0 <= eta <= 1.0 poly = ( @@ -283,9 +292,7 @@ def test_calc_alpha_beta_and_mixing_factor_behavior_and_wrappers(): L = pf._lorentzian_pink_beam(p_l, x) expected_full = eta * L + (1.0 - eta) * G + b0 + b1 * x out_wrapper = pf.pink_beam_dcs(p_full, x) - np.testing.assert_allclose( - out_wrapper, expected_full, rtol=1e-8, atol=1e-12 - ) + np.testing.assert_allclose(out_wrapper, expected_full, rtol=1e-8, atol=1e-12) x_lm = np.linspace(-0.8, 0.8, 101) out_lm = pf.pink_beam_dcs_lmfit( @@ -293,21 +300,15 @@ def test_calc_alpha_beta_and_mixing_factor_behavior_and_wrappers(): ) alpha_lm = pf._calc_alpha(np.array([alpha0, alpha1], dtype=float), x0) beta_lm = pf._calc_beta(np.array([beta0, beta1], dtype=float), x0) - p_g_lm = np.hstack( - ([A, x0], np.array([alpha_lm, beta_lm, fwhm_g], dtype=float)) - ) - p_l_lm = np.hstack( - ([A, x0], np.array([alpha_lm, beta_lm, fwhm_l], dtype=float)) - ) + p_g_lm = np.hstack(([A, x0], np.array([alpha_lm, beta_lm, fwhm_g], dtype=float))) + p_l_lm = np.hstack(([A, x0], np.array([alpha_lm, beta_lm, fwhm_l], dtype=float))) eta_lm, _ = pf._mixing_factor_pv(fwhm_g, fwhm_l) expected_lm = eta_lm * pf._lorentzian_pink_beam(p_l_lm, x_lm) + ( 1.0 - eta_lm ) * pf._gaussian_pink_beam(p_g_lm, x_lm) np.testing.assert_allclose(out_lm, expected_lm, rtol=1e-8, atol=1e-12) - p = np.array( - [10.0, 0.5, 0.01, 0.02, 0.03, 0.01, 0.15, 0.20, 1.0, 0.1], dtype=float - ) + p = np.array([10.0, 0.5, 0.01, 0.02, 0.03, 0.01, 0.15, 0.20, 1.0, 0.1], dtype=float) x_small = np.linspace(-1.0, 2.0, 7) y_no_bg = pf._pink_beam_dcs_no_bg(p[:-2], x_small) y_full = pf.pink_beam_dcs(p, x_small) @@ -346,9 +347,7 @@ def test_tanh_and_2d_coord_transform_and_gauss2d_rot_and_3d(): assert np.allclose(out2, out2[0, 0]) p_rot = np.array([1.0, 0.0, 0.0, 1.0, 1.0, np.pi / 3]) - out_rot = pf._gaussian2d_rot_no_bg( - p_rot, np.zeros((3, 3)), np.zeros((3, 3)) - ) + out_rot = pf._gaussian2d_rot_no_bg(p_rot, np.zeros((3, 3)), np.zeros((3, 3))) assert out_rot.shape == (3, 3) and np.isfinite(out_rot).all() p_full = np.array([1.0, 0.0, 0.0, 1.0, 1.0, 2.0, 0.5, -0.5]) @@ -359,20 +358,14 @@ def test_tanh_and_2d_coord_transform_and_gauss2d_rot_and_3d(): p_gr = np.array([1.0, 0.0, 0.0, 1.0, 1.0, np.pi / 4, 1.0, 0.3, 0.2]) out_gr = pf.gaussian2d_rot(p_gr, np.zeros((3, 3)), np.zeros((3, 3))) - base_r = pf._gaussian2d_rot_no_bg( - p_gr[:6], np.zeros((3, 3)), np.zeros((3, 3)) - ) + base_r = pf._gaussian2d_rot_no_bg(p_gr[:6], np.zeros((3, 3)), np.zeros((3, 3))) assert np.allclose(out_gr, base_r + 1.0 + 0.3 * 0 - 0.2 * 0) xq = np.array([[-1.0, 1.0], [-1.0, 1.0]]) yq = np.array([[-1.0, -1.0], [1.0, 1.0]]) p2d = np.array([1.0, 0.0, 0.0, 1.0, 2.0, 1.5, 2.5, 0.1, 0.2, 0.3, 0.4]) out2d = pf._split_pvoigt2d_no_bg(p2d, xq, yq) - assert ( - out2d.shape == (2, 2) - and np.isfinite(out2d).all() - and (out2d >= 0).all() - ) + assert out2d.shape == (2, 2) and np.isfinite(out2d).all() and (out2d >= 0).all() out2drot = pf._split_pvoigt2d_rot_no_bg( np.hstack((p2d, np.pi / 6)), np.zeros((3, 3)), np.zeros((3, 3)) @@ -380,9 +373,7 @@ def test_tanh_and_2d_coord_transform_and_gauss2d_rot_and_3d(): assert out2drot.shape == (3, 3) and np.isfinite(out2drot).all() p_split_bg = np.hstack((p2d, np.pi / 7, 1.0, 0.2, -0.1)) - out_sp_bg = pf.split_pvoigt2d_rot( - p_split_bg, np.zeros((4, 4)), np.zeros((4, 4)) - ) + out_sp_bg = pf.split_pvoigt2d_rot(p_split_bg, np.zeros((4, 4)), np.zeros((4, 4))) base_sp = pf._split_pvoigt2d_rot_no_bg( p_split_bg[:12], np.zeros((4, 4)), np.zeros((4, 4)) ) @@ -417,9 +408,7 @@ def _build_peak_row(pktype): if pktype == "split_pvoigt": return np.array([0.9, 0.0, 0.7, 1.1, 0.25, 0.75], dtype=float) if pktype == "pink_beam_dcs": - return np.array( - [1.3, 0.1, 0.01, 0.002, 0.02, -0.0003, 0.2, 0.4], dtype=float - ) + return np.array([1.3, 0.1, 0.01, 0.002, 0.02, -0.0003, 0.2, 0.4], dtype=float) raise ValueError From ae914536211d029c991ba78a77c760d0fa6a8c88 Mon Sep 17 00:00:00 2001 From: Patrick Avery Date: Mon, 8 Jun 2026 09:55:50 +0800 Subject: [PATCH 2/4] Upload coverage for Python3.11 only The numba flag that uploads coverage info for compiled functions only works for Python3.11. It should get fixed upstream sometime. Until then, there's no issue with just using Python3.11 for code coverage measurements. Signed-off-by: Patrick Avery --- .github/workflows/test.yml | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 011359a9c..e5703bbb0 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -10,6 +10,11 @@ jobs: pytest: name: ${{ matrix.config.name }} runs-on: ${{ matrix.config.os }} + # Single Python version that runs the coverage job. NUMBA_JIT_COVERAGE is + # incompatible with coverage.py's sys.monitoring tracer on Python 3.12+, + # so coverage only runs here; other versions run plain pytest. + env: + COVERAGE_PYTHON: '3.11' strategy: fail-fast: false matrix: @@ -64,7 +69,8 @@ jobs: HEXRD_EXAMPLE_REPO_PATH: ${{ github.workspace }}/examples run: | pytest -s tests/ - if: ${{ matrix.config.os != 'ubuntu-latest'}} + # Run plain tests everywhere except the single job that runs coverage. + if: ${{ matrix.config.os != 'ubuntu-latest' || matrix.python-version != env.COVERAGE_PYTHON }} working-directory: hexrd - name: Run tests with codecov @@ -75,7 +81,7 @@ jobs: NUMBA_JIT_COVERAGE=1 coverage run --source hexrd -m pytest tests/ coverage combine coverage xml -i - if: ${{ matrix.config.os == 'ubuntu-latest'}} + if: ${{ matrix.config.os == 'ubuntu-latest' && matrix.python-version == env.COVERAGE_PYTHON }} working-directory: hexrd - name: Upload coverage to Codecov @@ -84,4 +90,4 @@ jobs: token: ${{ secrets.CODECOV_TOKEN }} file: coverage.xml working-directory: hexrd - if: ${{ matrix.config.os == 'ubuntu-latest'}} + if: ${{ matrix.config.os == 'ubuntu-latest' && matrix.python-version == env.COVERAGE_PYTHON }} From 6d2d2b46b64fa9a03beaa0880369218602c555d3 Mon Sep 17 00:00:00 2001 From: Patrick Avery Date: Mon, 8 Jun 2026 10:09:11 +0800 Subject: [PATCH 3/4] Add type annotations and seed tests for rng We need to seed the tests for reproducible behavior. In theory, the tests should always pass, but seeding removes any flakiness. We also removed an unused import. Signed-off-by: Patrick Avery --- hexrd/core/fitting/peakfunctions.py | 2 +- hexrd/core/fitting/special.py | 72 ++++++++++++------------ tests/core/fitting/test_peakfunctions.py | 3 + 3 files changed, 40 insertions(+), 37 deletions(-) diff --git a/hexrd/core/fitting/peakfunctions.py b/hexrd/core/fitting/peakfunctions.py index ca2f131a0..b3d97daeb 100644 --- a/hexrd/core/fitting/peakfunctions.py +++ b/hexrd/core/fitting/peakfunctions.py @@ -29,7 +29,7 @@ from numba import njit import copy from hexrd.core import constants -from hexrd.core.fitting.special import erfc, exp1_complex_numba, wofz +from hexrd.core.fitting.special import erfc, exp1_complex_numba exp1exp = exp1_complex_numba diff --git a/hexrd/core/fitting/special.py b/hexrd/core/fitting/special.py index 21d28f03c..18ac00a98 100644 --- a/hexrd/core/fitting/special.py +++ b/hexrd/core/fitting/special.py @@ -23,8 +23,8 @@ _MIN_EXP = -745.1332191019411 # exp(x) underflows to 0 below this -@njit(cache=True, nopython=True) -def erfc(x): +@njit(cache=True, nogil=True) +def erfc(x: np.ndarray) -> np.ndarray: # save the sign of x sign = np.sign(x) x = np.abs(x) @@ -39,8 +39,8 @@ def erfc(x): return 1.0 - erf -@njit(cache=True, nopython=True) -def _complex_log_principal_numba(z): +@njit(cache=True, nogil=True) +def _complex_log_principal_numba(z: complex) -> complex: """ Principal complex log, written manually so signed zero in imag is handled through atan2. @@ -51,8 +51,8 @@ def _complex_log_principal_numba(z): ) -@njit(cache=True, nopython=True) -def _exp_neg_complex_numba(z): +@njit(cache=True, nogil=True) +def _exp_neg_complex_numba(z: complex) -> complex: """ Compute exp(-z) with simple overflow/underflow guards. """ @@ -70,8 +70,8 @@ def _exp_neg_complex_numba(z): return complex(scale * math.cos(y), -scale * math.sin(y)) -@njit(cache=True, nopython=True) -def _exp1_series_numba(z): +@njit(cache=True, nogil=True) +def _exp1_series_numba(z: complex) -> complex: """ Logarithmic power series: @@ -93,8 +93,8 @@ def _exp1_series_numba(z): return complex(-_EULER_GAMMA, 0.0) - _complex_log_principal_numba(z) - s -@njit(cache=True, nopython=True) -def _exp1_contfrac_numba(z): +@njit(cache=True, nogil=True) +def _exp1_contfrac_numba(z: complex) -> complex: """ Continued fraction via modified Lentz iteration. Good away from the origin and away from the difficult part @@ -127,8 +127,8 @@ def _exp1_contfrac_numba(z): return h * _exp_neg_complex_numba(z) -@njit(cache=True, nopython=True) -def exp1_complex_scalar_numba(z): +@njit(cache=True, nogil=True) +def exp1_complex_scalar_numba(z: complex) -> complex: """ Numba-compatible complex exponential integral E1(z). @@ -175,15 +175,15 @@ def exp1_complex_scalar_numba(z): @vectorize([complex128(complex128)], nopython=True) -def exp1_complex_numba(z): +def exp1_complex_numba(z: complex) -> complex: """ Vectorized complex ufunc version. """ return exp1_complex_scalar_numba(z) -@njit(cache=True, nopython=True) -def exp1_real_scalar_numba(x): +@njit(cache=True, nogil=True) +def exp1_real_scalar_numba(x: float) -> float: """ Real-valued SciPy-like wrapper. @@ -206,7 +206,7 @@ def exp1_real_scalar_numba(x): @vectorize([float64(float64)], nopython=True) -def exp1_real_numba(x): +def exp1_real_numba(x: float) -> float: """ Vectorized real ufunc version. """ @@ -220,8 +220,8 @@ def exp1_real_numba(x): # return exp1_complex_numba(x) -@njit(cache=True, nopython=True) -def _signed_inf_from_trig(v): +@njit(cache=True, nogil=True) +def _signed_inf_from_trig(v: float) -> float: if v > 0.0: return math.inf if v < 0.0: @@ -229,8 +229,8 @@ def _signed_inf_from_trig(v): return 0.0 -@njit(cache=True, nopython=True) -def _cexp_parts(a, b): +@njit(cache=True, nogil=True) +def _cexp_parts(a: float, b: float) -> complex: """ Compute exp(a + i b) with simple overflow/underflow guards. """ @@ -246,8 +246,8 @@ def _cexp_parts(a, b): return complex(ea * math.cos(b), ea * math.sin(b)) -@njit(cache=True, nopython=True) -def _exp_neg_z2(z): +@njit(cache=True, nogil=True) +def _exp_neg_z2(z: complex) -> complex: """ Compute exp(-z*z), where z = x + i y. -z^2 = (y^2 - x^2) - 2 i x y @@ -257,8 +257,8 @@ def _exp_neg_z2(z): return _cexp_parts(y * y - x * x, -2.0 * x * y) -@njit(cache=True, nopython=True) -def _trap_correction(z, plus): +@njit(cache=True, nogil=True) +def _trap_correction(z: complex, plus: bool) -> complex: """ Correction term: @@ -303,8 +303,8 @@ def _trap_correction(z, plus): return a / (1.0 - b) -@njit(cache=True, nopython=True) -def _wofz_midpoint(z): +@njit(cache=True, nogil=True) +def _wofz_midpoint(z: complex) -> complex: """ Midpoint-rule approximation used in part of the first quadrant. """ @@ -324,15 +324,15 @@ def _wofz_midpoint(z): @njit(cache=True, nogil=True) -def _wofz_modified_midpoint(z): +def _wofz_modified_midpoint(z: complex) -> complex: """ Modified midpoint-rule approximation. """ return _wofz_midpoint(z) + _trap_correction(z, True) -@njit(cache=True, nopython=True) -def _wofz_modified_trapezoid(z): +@njit(cache=True, nogil=True) +def _wofz_modified_trapezoid(z: complex) -> complex: """ Modified trapezoidal-rule approximation. """ @@ -351,8 +351,8 @@ def _wofz_modified_trapezoid(z): return complex(0.0, 1.0 / _H_CAP) / z + az * s + _trap_correction(z, False) -@njit(cache=True, nopython=True) -def _wofz_first_quadrant(z): +@njit(cache=True, nogil=True) +def _wofz_first_quadrant(z: complex) -> complex: """ Evaluate w(z) for Re(z) >= 0, Im(z) >= 0. """ @@ -372,8 +372,8 @@ def _wofz_first_quadrant(z): return _wofz_modified_midpoint(z) -@njit(cache=True, nopython=True) -def wofz_complex_scalar_numba(z): +@njit(cache=True, nogil=True) +def wofz_complex_scalar_numba(z: complex) -> complex: """ Numba-compatible scalar complex Faddeeva function. @@ -423,8 +423,8 @@ def wofz_complex_scalar_numba(z): return w -@njit(cache=True, nopython=True) -def wofz_real_scalar_numba(x): +@njit(cache=True, nogil=True) +def wofz_real_scalar_numba(x: float) -> complex: """ Convenience wrapper for real x, returning complex w(x). """ @@ -435,7 +435,7 @@ def wofz_real_scalar_numba(x): [complex128(complex128), complex128(float64)], nopython=True, ) -def wofz(z): +def wofz(z: complex) -> complex: """ Vectorized ufunc version. diff --git a/tests/core/fitting/test_peakfunctions.py b/tests/core/fitting/test_peakfunctions.py index 39e41e44f..c6f55b8f8 100644 --- a/tests/core/fitting/test_peakfunctions.py +++ b/tests/core/fitting/test_peakfunctions.py @@ -24,6 +24,7 @@ def finite_diff_derivative(func, p, x, idx, rel_step=1e-6): def test_erfc(): # test real values + np.random.seed(0) vals = np.random.uniform(-5, 5, 5) got = sp.erfc(vals) assert got.shape == vals.shape @@ -33,6 +34,7 @@ def test_erfc(): def test_wofz(): # test complex values + np.random.seed(0) vals = np.random.uniform(-5, 5, 5) + 1j * np.random.uniform(-5, 5, 5) got = sp.wofz(vals) assert got.shape == vals.shape @@ -41,6 +43,7 @@ def test_wofz(): def test_exp1(): # test real values + np.random.seed(0) vals = np.random.uniform(-5, 5, 5) got = sp.exp1_real_numba(vals) assert got.shape == vals.shape From 929af0c91d85bb1cc7f58de8f69ee14d3115afca Mon Sep 17 00:00:00 2001 From: Patrick Avery Date: Mon, 8 Jun 2026 15:35:06 +0800 Subject: [PATCH 4/4] Add missing exp in exp1exp function This is completely necessary for exp1exp to work correctly. This commit also removes the unnecessary threshold changes. Signed-off-by: Patrick Avery --- hexrd/core/fitting/peakfunctions.py | 4 +--- hexrd/core/fitting/special.py | 11 +++++++++++ hexrd/core/fitting/utils.py | 5 +---- hexrd/powder/wppf/peakfunctions.py | 12 +++--------- 4 files changed, 16 insertions(+), 16 deletions(-) diff --git a/hexrd/core/fitting/peakfunctions.py b/hexrd/core/fitting/peakfunctions.py index b3d97daeb..32095a328 100644 --- a/hexrd/core/fitting/peakfunctions.py +++ b/hexrd/core/fitting/peakfunctions.py @@ -29,9 +29,7 @@ from numba import njit import copy from hexrd.core import constants -from hexrd.core.fitting.special import erfc, exp1_complex_numba - -exp1exp = exp1_complex_numba +from hexrd.core.fitting.special import erfc, exp1exp gauss_width_fact = constants.sigma_to_fwhm lorentz_width_fact = 2.0 diff --git a/hexrd/core/fitting/special.py b/hexrd/core/fitting/special.py index 18ac00a98..ffb25c0d3 100644 --- a/hexrd/core/fitting/special.py +++ b/hexrd/core/fitting/special.py @@ -182,6 +182,17 @@ def exp1_complex_numba(z: complex) -> complex: return exp1_complex_scalar_numba(z) +@njit(cache=True, nogil=True) +def exp1exp(z: np.ndarray) -> np.ndarray: + """ + E1(z) * exp(z), the combination used by the back-to-back-exponential + pink-beam Lorentzian (Von Dreele et al., J. Appl. Cryst. (2021) 54, 3-6). + The exp(z) factor cancels the exponential growth of E1, keeping the + result well-behaved across the peak. + """ + return np.exp(z) * exp1_complex_numba(z) + + @njit(cache=True, nogil=True) def exp1_real_scalar_numba(x: float) -> float: """ diff --git a/hexrd/core/fitting/utils.py b/hexrd/core/fitting/utils.py index ed93fe1fd..ccbaa5026 100644 --- a/hexrd/core/fitting/utils.py +++ b/hexrd/core/fitting/utils.py @@ -2,13 +2,10 @@ import numpy as np from numba import njit -from hexrd.core.fitting.special import erfc, exp1_complex_numba +from hexrd.core.fitting.special import erfc, exp1exp from hexrd.core.matrixutil import uniqueVectors -# alias for exp1 function -exp1exp = exp1_complex_numba - # ============================================================================= # LMFIT Parameter munging utilities # ============================================================================= diff --git a/hexrd/powder/wppf/peakfunctions.py b/hexrd/powder/wppf/peakfunctions.py index 8c14faa34..309b51037 100644 --- a/hexrd/powder/wppf/peakfunctions.py +++ b/hexrd/powder/wppf/peakfunctions.py @@ -29,9 +29,8 @@ import copy from hexrd.core import constants from numba import vectorize, float64, njit, prange -from hexrd.core.fitting.special import erfc, exp1_complex_numba +from hexrd.core.fitting.special import erfc, exp1exp -exp1exp = exp1_complex_numba # from scipy.special import erfc, exp1 # addr = get_cython_function_address("scipy.special.cython_special", "exp1") @@ -522,7 +521,7 @@ def _gaussian_pink_beam(alpha, beta, fwhm_g, tth, tth_list): t1 = erfc(y) t2 = erfc(z) g = np.zeros(tth_list.shape) - zmask = np.abs(del_tth) > 100.0 + zmask = np.abs(del_tth) > 5.0 g[~zmask] = (0.5 * (alpha * beta) / (alpha + beta)) * np.exp(u[~zmask]) * t1[ ~zmask ] + np.exp(v[~zmask]) * t2[~zmask] @@ -551,12 +550,7 @@ def _lorentzian_pink_beam(alpha, beta, fwhm_l, tth, tth_list): # f1 = exp1(p) # f2 = exp1(q) - y = np.zeros(tth_list.shape) - zmask = np.abs(del_tth) > 100.0 - - y[~zmask] = ( - -(alpha * beta) / (np.pi * (alpha + beta)) * (f1[~zmask] + f2[~zmask]).imag - ) + y = -(alpha * beta) / (np.pi * (alpha + beta)) * (f1 + f2).imag mask = np.isnan(y) y[mask] = 0.0