From ae168e5aec911e279defb15ee07f126fd88fb5b7 Mon Sep 17 00:00:00 2001 From: Caitlyn Wolf Date: Sat, 29 Oct 2022 18:54:13 +0200 Subject: [PATCH 01/10] fixed switched nomenclature for length and width --- sasmodels/resolution.py | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/sasmodels/resolution.py b/sasmodels/resolution.py index 27d7b255..0b60e46f 100644 --- a/sasmodels/resolution.py +++ b/sasmodels/resolution.py @@ -151,7 +151,7 @@ def __init__(self, q, q_length, q_width=0., q_calc=None): # Build weight matrix from calculated q values self.weight_matrix = \ - slit_resolution(self.q_calc, self.q, q_length, q_width) + slit_resolution(self.q_calc, self.q, q_width, q_length) self.q_calc = abs(self.q_calc) def apply(self, theory): @@ -332,7 +332,6 @@ def slit_resolution(q_calc, q, width, length, n_length=30): """ - # np.set_printoptions(precision=6, linewidth=10000) # The current algorithm is a midpoint rectangle rule. @@ -349,29 +348,29 @@ def slit_resolution(q_calc, q, width, length, n_length=30): # up qi in q_calc, then weighting the result by the relative # distance to the neighbouring points. weights[i, :] = (q_calc == qi) - elif l == 0: - weights[i, :] = _q_perp_weights(q_edges, qi, w) elif w == 0: - in_x = 1.0 * ((q_calc >= qi-l) & (q_calc <= qi+l)) - abs_x = 1.0*(q_calc < abs(qi - l)) if qi < l else 0. - #print(qi - l, qi + l) + weights[i, :] = _q_perp_weights(q_edges, qi, l) + elif l == 0: + in_x = 1.0 * ((q_calc >= qi-w) & (q_calc <= qi+w)) + abs_x = 1.0*(q_calc < abs(qi - w)) if qi < w else 0. + #print(qi - w, qi + w) #print(in_x + abs_x) - weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*l) + weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*w) else: for k in range(-n_length, n_length+1): - weights[i, :] += _q_perp_weights(q_edges, qi+k*l/n_length, w) + weights[i, :] += _q_perp_weights(q_edges, qi+k*w/n_length, l) weights[i, :] /= 2*n_length + 1 return weights.T -def _q_perp_weights(q_edges, qi, w): +def _q_perp_weights(q_edges, qi, l): # Convert bin edges from q to u - u_limit = np.sqrt(qi**2 + w**2) + u_limit = np.sqrt(qi**2 + l**2) u_edges = q_edges**2 - qi**2 u_edges[q_edges < abs(qi)] = 0. u_edges[q_edges > u_limit] = u_limit**2 - qi**2 - weights = np.diff(np.sqrt(u_edges))/w + weights = np.diff(np.sqrt(u_edges))/l #print("i, qi",i,qi,qi+width) #print(q_calc) #print(weights) From fd03bdaca12b3384d3a02ca8d0635fb4a8da66ef Mon Sep 17 00:00:00 2001 From: Caitlyn Wolf Date: Mon, 31 Oct 2022 21:20:32 +0100 Subject: [PATCH 02/10] fixed single letter variable names of 'w' and 'l' --- sasmodels/resolution.py | 62 ++++++++++++++++++++--------------------- 1 file changed, 31 insertions(+), 31 deletions(-) diff --git a/sasmodels/resolution.py b/sasmodels/resolution.py index 0b60e46f..16ad9e33 100644 --- a/sasmodels/resolution.py +++ b/sasmodels/resolution.py @@ -340,37 +340,37 @@ def slit_resolution(q_calc, q, width, length, n_length=30): weights = np.zeros((len(q), len(q_calc)), 'd') #print(q_calc) - for i, (qi, w, l) in enumerate(zip(q, width, length)): - if w == 0. and l == 0.: + for i, (qi, wi, li) in enumerate(zip(q, width, length)): + if wi == 0. and li == 0.: # Perfect resolution, so return the theory value directly. # Note: assumes that q is a subset of q_calc. If qi need not be # in q_calc, then we can do a weighted interpolation by looking # up qi in q_calc, then weighting the result by the relative # distance to the neighbouring points. weights[i, :] = (q_calc == qi) - elif w == 0: - weights[i, :] = _q_perp_weights(q_edges, qi, l) - elif l == 0: - in_x = 1.0 * ((q_calc >= qi-w) & (q_calc <= qi+w)) - abs_x = 1.0*(q_calc < abs(qi - w)) if qi < w else 0. + elif wi == 0: + weights[i, :] = _q_perp_weights(q_edges, qi, li) + elif li == 0: + in_x = 1.0 * ((q_calc >= qi-wi) & (q_calc <= qi+wi)) + abs_x = 1.0*(q_calc < abs(qi - wi)) if qi < wi else 0. #print(qi - w, qi + w) #print(in_x + abs_x) - weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*w) + weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*wi) else: for k in range(-n_length, n_length+1): - weights[i, :] += _q_perp_weights(q_edges, qi+k*w/n_length, l) + weights[i, :] += _q_perp_weights(q_edges, qi+k*wi/n_length, li) weights[i, :] /= 2*n_length + 1 return weights.T -def _q_perp_weights(q_edges, qi, l): +def _q_perp_weights(q_edges, qi, length): # Convert bin edges from q to u - u_limit = np.sqrt(qi**2 + l**2) + u_limit = np.sqrt(qi**2 + length**2) u_edges = q_edges**2 - qi**2 u_edges[q_edges < abs(qi)] = 0. u_edges[q_edges > u_limit] = u_limit**2 - qi**2 - weights = np.diff(np.sqrt(u_edges))/l + weights = np.diff(np.sqrt(u_edges))/length #print("i, qi",i,qi,qi+width) #print(q_calc) #print(weights) @@ -566,8 +566,8 @@ def romberg_slit_1d(q, width, length, form, pars): width = [width]*len(q) if np.isscalar(length): length = [length]*len(q) - _int_w = lambda w, qi: eval_form(sqrt(qi**2 + w**2), form, pars) - _int_l = lambda l, qi: eval_form(abs(qi+l), form, pars) + _int_w = lambda wi, qi: eval_form(sqrt(qi**2 + wi**2), form, pars) + _int_l = lambda li, qi: eval_form(abs(qi+li), form, pars) # If both width and length are defined, then it is too slow to use dblquad. # Instead use trapz on a fixed grid, interpolated into the I(Q) for # the extended Q range. @@ -575,23 +575,23 @@ def romberg_slit_1d(q, width, length, form, pars): q_calc = slit_extend_q(q, np.asarray(width), np.asarray(length)) Iq = eval_form(q_calc, form, pars) result = np.empty(len(q)) - for i, (qi, w, l) in enumerate(zip(q, width, length)): - if l == 0.: - total = romberg(_int_w, 0, w, args=(qi,), + for i, (qi, wi, li) in enumerate(zip(q, width, length)): + if li == 0.: + total = romberg(_int_w, 0, wi, args=(qi,), divmax=100, vec_func=True, tol=0, rtol=1e-8) - result[i] = total/w - elif w == 0.: - total = romberg(_int_l, -l, l, args=(qi,), + result[i] = total/wi + elif wi == 0.: + total = romberg(_int_l, -li, li, args=(qi,), divmax=100, vec_func=True, tol=0, rtol=1e-8) - result[i] = total/(2*l) + result[i] = total/(2*li) else: - w_grid = np.linspace(0, w, 21)[None, :] - l_grid = np.linspace(-l, l, 23)[:, None] + w_grid = np.linspace(0, wi, 21)[None, :] + l_grid = np.linspace(-li, li, 23)[:, None] u_sub = sqrt((qi+l_grid)**2 + w_grid**2) f_at_u = np.interp(u_sub, q_calc, Iq) #print(np.trapz(Iu, w_grid, axis=1)) total = np.trapz(np.trapz(f_at_u, w_grid, axis=1), l_grid[:, 0]) - result[i] = total / (2*l*w) + result[i] = total / (2*li*wi) # from scipy.integrate import dblquad # r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w, # args=(qi,)) @@ -1160,13 +1160,13 @@ def demo_slit_1d(): Show example of slit smearing. """ q = np.logspace(-4, np.log10(0.2), 100) - w = l = 0. - #w = 0.000000277790 - w = 0.0277790 - #l = 0.00277790 - #l = 0.0277790 - resolution = Slit1D(q, w, l) - _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w, l)) + width = length = 0. + #width = 0.000000277790 + width = 0.0277790 + #length = 0.00277790 + #length = 0.0277790 + resolution = Slit1D(q, width, length) + _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(width, length)) def demo(): """ From 81ffb2e92dd64b706b1cb9b44ced8dbf9b4bdd7b Mon Sep 17 00:00:00 2001 From: Caitlyn Wolf Date: Tue, 27 Jun 2023 13:14:35 -0400 Subject: [PATCH 03/10] Fixed inverse q-spacing in the geometric extrapolation function --- sasmodels/resolution.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sasmodels/resolution.py b/sasmodels/resolution.py index 27d7b255..24dc3bf8 100644 --- a/sasmodels/resolution.py +++ b/sasmodels/resolution.py @@ -500,7 +500,7 @@ def geometric_extrapolation(q, q_min, q_max, points_per_decade=None): if points_per_decade is None: log_delta_q = (len(q) - 1) / (log(q[-1]) - log(q[0])) else: - log_delta_q = log(10.) / points_per_decade + log_delta_q = points_per_decade / log(10.) if q_min < q[0]: if q_min < 0: q_min = q[0]*MINIMUM_ABSOLUTE_Q From f07dc5df87e9d7aabf9a36befcc2921acddf06ff Mon Sep 17 00:00:00 2001 From: Caitlyn Wolf Date: Tue, 27 Jun 2023 13:38:51 -0400 Subject: [PATCH 04/10] fixing typo from resolving merge conflict to ensure inverse term of log_delta_q is calculated when points_per_decade is specified --- sasmodels/resolution.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sasmodels/resolution.py b/sasmodels/resolution.py index 05750718..2827ea3a 100644 --- a/sasmodels/resolution.py +++ b/sasmodels/resolution.py @@ -514,7 +514,7 @@ def geometric_extrapolation(q, q_min, q_max, points_per_decade=None): else: log_delta_q = DEFAULT_POINTS_PER_DECADE / log(10.) else: - points_per_decade / log_delta_q = log(10.) + log_delta_q = points_per_decade / log(10.) if q_min < data_min: if q_min < 0: q_min = data_min*MINIMUM_ABSOLUTE_Q From 563ad5b6d53e981cfca864d67a3cb0b58c7dd5fa Mon Sep 17 00:00:00 2001 From: Caitlyn Wolf Date: Tue, 16 Jan 2024 17:31:00 -0500 Subject: [PATCH 05/10] switched order of length and width in slit_resolution method for consistency with Slit1D --- sasmodels/resolution.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sasmodels/resolution.py b/sasmodels/resolution.py index 5c3ec860..5c693de5 100644 --- a/sasmodels/resolution.py +++ b/sasmodels/resolution.py @@ -157,7 +157,7 @@ def __init__(self, q, q_length=None, q_width=None, q_calc=None): # Build weight matrix from calculated q values self.weight_matrix = ( - slit_resolution(self.q_calc, self.q, q_width, q_length) + slit_resolution(self.q_calc, self.q, q_length, q_width) ) self.q_calc = abs(self.q_calc) @@ -216,7 +216,7 @@ def pinhole_resolution(q_calc, q, q_width, nsigma=PINHOLE_N_SIGMA): return weights -def slit_resolution(q_calc, q, width, length, n_length=30): +def slit_resolution(q_calc, q, length, width, n_length=30): r""" Build a weight matrix to compute *I_s(q)* from *I(q_calc)*, given $q_\perp$ = *width* (in the high-resolution axis) and $q_\parallel$ From d0dfa73ff5463831624e4529e23bdc7700c17c36 Mon Sep 17 00:00:00 2001 From: Caitlyn Wolf Date: Tue, 16 Jan 2024 22:38:45 -0500 Subject: [PATCH 06/10] updated romberg_slit_1d to use length, width ordering of arguments, this also led to updating ordering of length and width arguments for slit_extend_q --- sasmodels/resolution.py | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/sasmodels/resolution.py b/sasmodels/resolution.py index 5c693de5..ed82e0b1 100644 --- a/sasmodels/resolution.py +++ b/sasmodels/resolution.py @@ -146,7 +146,7 @@ def __init__(self, q, q_length=None, q_width=None, q_calc=None): self.q = q.flatten() self.q_calc = ( - slit_extend_q(q, q_width, q_length) + slit_extend_q(q, q_length, q_width) if q_calc is None else np.sort(q_calc) ) @@ -398,13 +398,13 @@ def pinhole_extend_q(q, q_width, nsigma=PINHOLE_N_SIGMA): return linear_extrapolation(q, q_min, q_max) -def slit_extend_q(q, width, length): +def slit_extend_q(q, length, width): """ Given *q*, *width* and *length*, find a set of sampling points *q_calc* so that each point I(q) has sufficient support from the underlying function. """ - q_min, q_max = np.min(q-length), np.max(np.sqrt((q+length)**2 + width**2)) + q_min, q_max = np.min(q-width), np.max(np.sqrt((q+width)**2 + length**2)) return geometric_extrapolation(q, q_min, q_max) @@ -560,7 +560,7 @@ def gaussian(q, q0, dq, nsigma=2.5): return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)/(1-two_tail_density) -def romberg_slit_1d(q, width, length, form, pars): +def romberg_slit_1d(q, length, width, form, pars): """ Romberg integration for slit resolution. @@ -580,31 +580,31 @@ def romberg_slit_1d(q, width, length, form, pars): width = [width]*len(q) if np.isscalar(length): length = [length]*len(q) - _int_w = lambda wi, qi: eval_form(sqrt(qi**2 + wi**2), form, pars) - _int_l = lambda li, qi: eval_form(abs(qi+li), form, pars) + _int_l = lambda li, qi: eval_form(sqrt(qi**2 + li**2), form, pars) + _int_w = lambda wi, qi: eval_form(abs(qi+wi), form, pars) # If both width and length are defined, then it is too slow to use dblquad. # Instead use trapz on a fixed grid, interpolated into the I(Q) for # the extended Q range. #_int_wh = lambda w, h, qi: eval_form(sqrt((qi+h)**2 + w**2), form, pars) - q_calc = slit_extend_q(q, np.asarray(width), np.asarray(length)) + q_calc = slit_extend_q(q, np.asarray(length), np.asarray(width)) Iq = eval_form(q_calc, form, pars) result = np.empty(len(q)) for i, (qi, wi, li) in enumerate(zip(q, width, length)): - if li == 0.: - total = romberg(_int_w, 0, wi, args=(qi,), + if wi == 0.: + total = romberg(_int_l, 0, li, args=(qi,), divmax=100, vec_func=True, tol=0, rtol=1e-8) - result[i] = total/wi - elif wi == 0.: - total = romberg(_int_l, -li, li, args=(qi,), + result[i] = total/li + elif li == 0.: + total = romberg(_int_w, -wi, wi, args=(qi,), divmax=100, vec_func=True, tol=0, rtol=1e-8) - result[i] = total/(2*li) + result[i] = total/(2*wi) else: - w_grid = np.linspace(0, wi, 21)[None, :] - l_grid = np.linspace(-li, li, 23)[:, None] - u_sub = sqrt((qi+l_grid)**2 + w_grid**2) + l_grid = np.linspace(0, li, 21)[None, :] + w_grid = np.linspace(-wi, wi, 23)[:, None] + u_sub = sqrt((qi+w_grid)**2 + l_grid**2) f_at_u = np.interp(u_sub, q_calc, Iq) #print(np.trapz(Iu, w_grid, axis=1)) - total = np.trapz(np.trapz(f_at_u, w_grid, axis=1), l_grid[:, 0]) + total = np.trapz(np.trapz(f_at_u, l_grid, axis=1), w_grid[:, 0]) result[i] = total / (2*li*wi) # from scipy.integrate import dblquad # r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w, From 0e0b63d4bfe4628d9565e361a3070b6ac2bf5215 Mon Sep 17 00:00:00 2001 From: Caitlyn Wolf Date: Tue, 16 Jan 2024 22:51:50 -0500 Subject: [PATCH 07/10] fixed typo in the doc string for resolution.slit_resolution --- sasmodels/resolution.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sasmodels/resolution.py b/sasmodels/resolution.py index ed82e0b1..0110bbd0 100644 --- a/sasmodels/resolution.py +++ b/sasmodels/resolution.py @@ -282,7 +282,7 @@ def slit_resolution(q_calc, q, length, width, n_length=30): where $I(u_j)$ is the value at the mid-point, and $\Delta u_j$ is the difference between consecutive edges which have been first converted to $u$. Only $u_j \in [0, \Delta q_\perp]$ are used, which corresponds - to $q_j \in \left[q, \sqrt{q^2 + \Delta q_\perp}\right]$, so + to $q_j \in \left[q, \sqrt{q^2 + \Delta q_\perp^2}\right]$, so .. math:: From aa2b7b2cb0b3aae5e04002052e496a1756d78db3 Mon Sep 17 00:00:00 2001 From: Caitlyn Wolf Date: Wed, 17 Jan 2024 14:44:07 -0500 Subject: [PATCH 08/10] fixed width, length ordering in slit smearing fucntions to length, width --- sasmodels/resolution.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/sasmodels/resolution.py b/sasmodels/resolution.py index 0110bbd0..62e6a019 100644 --- a/sasmodels/resolution.py +++ b/sasmodels/resolution.py @@ -1146,7 +1146,7 @@ def _eval_demo_1d(resolution, title): if isinstance(resolution, Slit1D): width, length = resolution.q_width, resolution.q_length - Iq_romb = romberg_slit_1d(resolution.q, width, length, model, pars) + Iq_romb = romberg_slit_1d(resolution.q, length, width, model, pars) else: dq = resolution.q_width Iq_romb = romberg_pinhole_1d(resolution.q, dq, model, pars) @@ -1176,11 +1176,11 @@ def demo_slit_1d(): q = np.logspace(-4, np.log10(0.2), 100) width = length = 0. #width = 0.000000277790 - width = 0.0277790 + length = 0.0277790 #length = 0.00277790 #length = 0.0277790 - resolution = Slit1D(q, width, length) - _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(width, length)) + resolution = Slit1D(q, length, width) + _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(length, width)) def demo(): """ From e57e7763087ad6ac4d59812b94d9feb6974d26b3 Mon Sep 17 00:00:00 2001 From: Caitlyn Wolf Date: Wed, 17 Jan 2024 15:42:10 -0500 Subject: [PATCH 09/10] testing opencl ubuntu tests on github ci --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index ff0441e5..a57c8212 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -8,7 +8,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - os: [macos-latest, ubuntu-latest, windows-latest] + os: [macos-latest, ubuntu-22.04, windows-latest] python-version: ["3.7", "3.8", "3.9", "3.10"] steps: From 98053f4fe246668e5247487785c0c7aeb4a1773b Mon Sep 17 00:00:00 2001 From: Caitlyn Wolf Date: Thu, 18 Jan 2024 11:09:48 -0500 Subject: [PATCH 10/10] switched the log_delta_q parameter in the first place so the variable name is accurate --- sasmodels/resolution.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/sasmodels/resolution.py b/sasmodels/resolution.py index 2827ea3a..619296af 100644 --- a/sasmodels/resolution.py +++ b/sasmodels/resolution.py @@ -510,20 +510,20 @@ def geometric_extrapolation(q, q_min, q_max, points_per_decade=None): data_min, data_max = q[0], q[-1] if points_per_decade is None: if data_max > data_min: - log_delta_q = (len(q) - 1) / (log(data_max) - log(data_min)) + log_delta_q = (log(data_max) - log(data_min)) / (len(q) - 1) else: - log_delta_q = DEFAULT_POINTS_PER_DECADE / log(10.) + log_delta_q = log(10.) / DEFAULT_POINTS_PER_DECADE else: - log_delta_q = points_per_decade / log(10.) + log_delta_q = log(10.) / points_per_decade if q_min < data_min: if q_min < 0: q_min = data_min*MINIMUM_ABSOLUTE_Q - n_low = int(np.ceil(log_delta_q * (log(q[0])-log(q_min)))) + n_low = int(np.ceil((log(q[0])-log(q_min)) / log_delta_q)) q_low = np.logspace(log10(q_min), log10(q[0]), n_low+1)[:-1] else: q_low = [] if q_max > data_max: - n_high = int(np.ceil(log_delta_q * (log(q_max)-log(data_max)))) + n_high = int(np.ceil((log(q_max)-log(data_max)) / log_delta_q)) q_high = np.logspace(log10(data_max), log10(q_max), n_high+1)[1:] else: q_high = []