Changeset 80ba1a2 in sasview for src/sas/models


Ignore:
Timestamp:
Nov 12, 2015 6:11:10 PM (9 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
d430ee8, 662d8d87
Parents:
bc873053
Message:

Fix slit resolution for q_parallel and mixed q_perp and q_parallel. Fixes #474.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/sas/models/resolution.py

    rbc873053 r80ba1a2  
    1010 
    1111MINIMUM_RESOLUTION = 1e-8 
     12 
     13 
     14# When extrapolating to -q, what is the minimum positive q relative to q_min 
     15# that we wish to calculate? 
     16MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION = 0.01 
    1217 
    1318class Resolution(object): 
     
    148153 
    149154 
    150 def slit_resolution(q_calc, q, width, height): 
     155def slit_resolution(q_calc, q, width, height, n_height=30): 
    151156    r""" 
    152157    Build a weight matrix to compute *I_s(q)* from *I(q_calc)*, given 
    153     $q_v$ = *width* and $q_h$ = *height*. 
    154  
    155     *width* and *height* are scalars since current instruments use the 
    156     same slit settings for all measurement points. 
     158    $q_\perp$ = *width* and $q_\parallel$ = *height*.  *n_height* is 
     159    is the number of steps to use in the integration over $q_\parallel$ 
     160    when both $q_\perp$ and $q_\parallel$ are non-zero. 
     161 
     162    Each $q$ can have an independent width and height value even though 
     163    current instruments use the same slit setting for all measured points. 
    157164 
    158165    If slit height is large relative to width, use: 
     
    160167    .. math:: 
    161168 
    162         I_s(q_o) = \frac{1}{\Delta q_v} 
    163             \int_0^{\Delta q_v} I(\sqrt{q_o^2 + u^2} du 
     169        I_s(q_i) = \frac{1}{\Delta q_\perp} 
     170            \int_0^{\Delta q_\perp} I(\sqrt{q_i^2 + q_\perp^2} dq_\perp 
    164171 
    165172    If slit width is large relative to height, use: 
     
    167174    .. math:: 
    168175 
    169         I_s(q_o) = \frac{1}{2 \Delta q_v} 
    170             \int_{-\Delta q_v}^{\Delta q_v} I(u) du 
     176        I_s(q_i) = \frac{1}{2 \Delta q_\parallel} 
     177            \int_{-\Delta q_\parallel}^{\Delta q_\parallel} 
     178                I(|q_i + q_\parallel|) dq_\parallel 
     179 
     180    For a mixture of slit width and height use: 
     181 
     182    .. math:: 
     183 
     184        I_s(q_i) = \frac{1}{2 \Delta q_\parallel \Delta q_\perp} 
     185            \int_{-\Delta q_\parallel)^{\Delta q_parallel} 
     186            \int_0^[\Delta q_\perp} 
     187                I(\sqrt{(q_i + q_\parallel)^2 + q_\perp^2}) 
     188                dq_\perp dq_\parallel 
     189 
     190 
     191    Algorithm 
     192    --------- 
     193 
     194    We are using the mid-point integration rule to assign weights to each 
     195    element of a weight matrix $W$ so that 
     196 
     197    .. math:: 
     198 
     199        I_s(q) = W I(q_\text{calc}) 
     200 
     201    If *q_calc* is at the mid-point, we can infer the bin edges from the 
     202    pairwise averages of *q_calc*, adding the missing edges before 
     203    *q_calc[0]* and after *q_calc[-1]*. 
     204 
     205    For $q_\parallel = 0$, the smeared value can be computed numerically 
     206    using the $u$ substitution 
     207 
     208    .. math:: 
     209 
     210        u_j = \sqrt{q_j^2 - q^2} 
     211 
     212    This gives 
     213 
     214    .. math:: 
     215 
     216        I_s(q) \approx \sum_j I(u_j) \Delta u_j 
     217 
     218    where $I(u_j)$ is the value at the mid-point, and $\Delta u_j$ is the 
     219    difference between consecutive edges which have been first converted 
     220    to $u$.  Only $u_j \in [0, \Delta q_\perp]$ are used, which corresponds 
     221    to $q_j \in [q, \sqrt{q^2 + \Delta q_\perp}]$, so 
     222 
     223    .. math:: 
     224 
     225        W_{ij} = \frac{1}{\Delta q_\perp} \Delta u_j 
     226               = \frac{1}{\Delta q_\perp} 
     227                    \sqrt{q_{j+1}^2 - q_i^2} - \sqrt{q_j^2 - q_i^2} 
     228            \text{if} q_j \in [q_i, \sqrt{q_i^2 + q_\perp^2}] 
     229 
     230    where $I_s(q_i)$ is the theory function being computed and $q_j$ are the 
     231    mid-points between the calculated values in *q_calc*.  We tweak the 
     232    edges of the initial and final intervals so that they lie on integration 
     233    limits. 
     234 
     235    (To be precise, the transformed midpoint $u(q_j)$ is not necessarily the 
     236    midpoint of the edges $u((q_{j-1}+q_j)/2)$ and $u((q_j + q_{j+1})/2)$, 
     237    but it is at least in the interval, so the approximation is going to be 
     238    a little better than the left or right Riemann sum, and should be 
     239    good enough for our purposes.) 
     240 
     241    For $q_\perp = 0$, the $u$ substitution is simpler: 
     242 
     243    .. math:: 
     244 
     245        u_j = |q_j - q| 
     246 
     247    so 
     248 
     249    .. math:: 
     250 
     251        W_ij = \frac{1}{2 \Delta q_\parallel} \Delta u_j 
     252            = \frac{1}{2 \Delta q_\parallel} (q_{j+1} - q_j) 
     253            \text{if} q_j \in [q-\Delta q_\parallel, q+\Delta q_\parallel] 
     254 
     255    However, we need to support cases were $u_j < 0$, which means using 
     256    $2 (q_{j+1} - q_j)$ when $q_j \in [0, q_\parallel-q_i]$.  This is not 
     257    an issue for $q_i > q_\parallel$. 
     258 
     259    For bot $q_\perp > 0$ and $q_\parallel > 0$ we perform a 2 dimensional 
     260    integration with 
     261 
     262    .. math:: 
     263 
     264        u_jk = \sqrt{q_j^2 - (q + (k\Delta q_\parallel/L))^2} 
     265            \text{for} k = -L \ldots L 
     266 
     267    for $L$ = *n_height*.  This gives 
     268 
     269    .. math:: 
     270 
     271        W_{ij} = \frac{1}{2 \Delta q_\perp q_\parallel} 
     272            \sum_{k=-L}^L \Delta u_jk (\frac{\Delta q_\parallel}{2 L + 1} 
     273 
     274 
    171275    """ 
    172276    #np.set_printoptions(precision=6, linewidth=10000) 
     
    177281    weights = np.zeros((len(q), len(q_calc)), 'd') 
    178282 
     283    #print q_calc 
    179284    for i, (qi, w, h) in enumerate(zip(q, width, height)): 
    180285        if w == 0. and h == 0.: 
     
    185290            # distance to the neighbouring points. 
    186291            weights[i, :] = (q_calc == qi) 
    187         elif h == 0 or w > 100.0 * h: 
    188             # Convert bin edges from q to u 
    189             u_limit = np.sqrt(qi**2 + w**2) 
    190             u_edges = q_edges**2 - qi**2 
    191             u_edges[q_edges < qi] = 0. 
    192             u_edges[q_edges > u_limit] = u_limit**2 - qi**2 
    193             weights[i, :] = np.diff(np.sqrt(u_edges))/w 
    194             #print "i, qi",i,qi,qi+width 
    195             #print q_calc 
    196             #print weights[i,:] 
     292        elif h == 0: 
     293            weights[i, :] = _q_perp_weights(q_edges, qi, w) 
    197294        elif w == 0: 
    198             in_x = (q_calc >= qi-h) * (q_calc <= qi+h) 
    199             weights[i,:] = in_x * np.diff(q_edges) / (2*h) 
     295            in_x = 1.0 * ((q_calc >= qi-h) & (q_calc <= qi+h)) 
     296            abs_x = 1.0*(q_calc < abs(qi - h)) if qi < h else 0. 
     297            #print qi - h, qi + h 
     298            #print in_x + abs_x 
     299            weights[i,:] = (in_x + abs_x) * np.diff(q_edges) / (2*h) 
     300        else: 
     301            L = n_height 
     302            for k in range(-L, L+1): 
     303                weights[i,:] += _q_perp_weights(q_edges, qi+k*h/L, w) 
     304            weights[i,:] /= 2*L + 1 
    200305 
    201306    return weights.T 
     307 
     308 
     309def _q_perp_weights(q_edges, qi, w): 
     310    # Convert bin edges from q to u 
     311    u_limit = np.sqrt(qi**2 + w**2) 
     312    u_edges = q_edges**2 - qi**2 
     313    u_edges[q_edges < abs(qi)] = 0. 
     314    u_edges[q_edges > u_limit] = u_limit**2 - qi**2 
     315    weights = np.diff(np.sqrt(u_edges))/w 
     316    #print "i, qi",i,qi,qi+width 
     317    #print q_calc 
     318    #print weights 
     319    return weights 
    202320 
    203321 
     
    218336    function. 
    219337    """ 
    220     q_min, q_max = np.min(q), np.max(np.sqrt(q**2 + width**2)) 
     338    q_min, q_max = np.min(q-height), np.max(np.sqrt((q+height)**2 + width**2)) 
     339 
    221340    return geometric_extrapolation(q, q_min, q_max) 
    222341 
     
    270389    q = np.sort(q) 
    271390    if q_min < q[0]: 
    272         if q_min <= 0: q_min = q[0]/10 
     391        if q_min <= 0: q_min = q_min*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 
    273392        n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1]>q[0] else 15 
    274393        q_low = np.linspace(q_min, q[0], n_low+1)[:-1] 
     
    320439        log_delta_q = log(10.) / points_per_decade 
    321440    if q_min < q[0]: 
    322         if q_min < 0: q_min = q[0]/10 
     441        if q_min < 0: q_min = q[0]*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 
    323442        n_low = log_delta_q * (log(q[0])-log(q_min)) 
    324443        q_low  = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1] 
     
    373492    _int_w = lambda w, qi: eval_form(sqrt(qi**2 + w**2), form, pars) 
    374493    _int_h = lambda h, qi: eval_form(abs(qi+h), form, pars) 
    375     _int_wh = lambda w, h, qi: eval_form(sqrt((qi+h)**2 + w**2), form, pars) 
     494    # If both width and height are defined, then it is too slow to use dblquad. 
     495    # Instead use trapz on a fixed grid, interpolated into the I(Q) for 
     496    # the extended Q range. 
     497    #_int_wh = lambda w, h, qi: eval_form(sqrt((qi+h)**2 + w**2), form, pars) 
     498    q_calc = slit_extend_q(q, np.asarray(width), np.asarray(height)) 
     499    Iq = eval_form(q_calc, form, pars) 
    376500    result = np.empty(len(q)) 
    377501    for i, (qi, w, h) in enumerate(zip(q, width, height)): 
     
    385509            result[i] = r/(2*h) 
    386510        else: 
    387             r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: width, 
     511            w_grid = np.linspace(0, w, 21)[None,:] 
     512            h_grid = np.linspace(-h, h, 23)[:,None] 
     513            u = sqrt((qi+h_grid)**2 + w_grid**2) 
     514            Iu = np.interp(u, q_calc, Iq) 
     515            #print np.trapz(Iu, w_grid, axis=1) 
     516            Is = np.trapz(np.trapz(Iu, w_grid, axis=1), h_grid[:,0]) 
     517            result[i] = Is / (2*h*w) 
     518            """ 
     519            r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w, 
    388520                             args=(qi,)) 
    389521            result[i] = r/(w*2*h) 
     522            """ 
    390523 
    391524    # r should be [float, ...], but it is [array([float]), array([float]),...] 
     
    9201053 
    9211054def demo_slit_1d(): 
    922     q = np.logspace(-4, np.log10(0.2), 400) 
     1055    q = np.logspace(-4, np.log10(0.2), 100) 
    9231056    w = h = 0. 
    924     #w = 0.005 
     1057    #w = 0.000000277790 
    9251058    w = 0.0277790 
    9261059    #h = 0.00277790 
     1060    #h = 0.0277790 
    9271061    resolution = Slit1D(q, w, h) 
    9281062    _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w,h)) 
     
    9401074 
    9411075if __name__ == "__main__": 
    942     #demo() 
    943     main() 
     1076    demo() 
     1077    #main() 
Note: See TracChangeset for help on using the changeset viewer.