Changes in sasmodels/resolution.py [2d81cfe:0b9c6df] in sasmodels


Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/resolution.py

    r2d81cfe r0b9c6df  
    2020MINIMUM_RESOLUTION = 1e-8 
    2121MINIMUM_ABSOLUTE_Q = 0.02  # relative to the minimum q in the data 
     22PINHOLE_N_SIGMA = 2.5 # From: Barker & Pedersen 1995 JAC 
    2223 
    2324class Resolution(object): 
    … …  
    6566    *q_calc* is the list of points to calculate, or None if this should 
    6667    be estimated from the *q* and *q_width*. 
    67     """ 
    68     def __init__(self, q, q_width, q_calc=None, nsigma=3): 
     68 
     69    *nsigma* is the width of the resolution function.  Should be 2.5. 
     70    See :func:`pinhole_resolution` for details. 
     71    """ 
     72    def __init__(self, q, q_width, q_calc=None, nsigma=PINHOLE_N_SIGMA): 
    6973        #*min_step* is the minimum point spacing to use when computing the 
    7074        #underlying model.  It should be on the order of 
    … …  
    8286 
    8387        # Protect against models which are not defined for very low q.  Limit 
    84         # the smallest q value evaluated (in absolute) to 0.02*min 
     88        # the smallest q value evaluated to 0.02*min.  Note that negative q 
     89        # values are trimmed even for broad resolution.  Although not possible 
     90        # from the geometry, they may appear since we are using a truncated 
     91        # gaussian to represent resolution rather than a skew distribution. 
    8592        cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q) 
    86         self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff] 
     93        self.q_calc = self.q_calc[self.q_calc >= cutoff] 
    8794 
    8895        # Build weight matrix from calculated q values 
    8996        self.weight_matrix = pinhole_resolution( 
    90             self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION)) 
    91         self.q_calc = abs(self.q_calc) 
     97            self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION), 
     98            nsigma=nsigma) 
    9299 
    93100    def apply(self, theory): 
    … …  
    101108    *q* points at which the data is measured. 
    102109 
    103     *dqx* slit width in qx 
    104  
    105     *dqy* slit height in qy 
     110    *qx_width* slit width in qx 
     111 
     112    *qy_width* slit height in qy 
    106113 
    107114    *q_calc* is the list of points to calculate, or None if this should 
    … …  
    154161 
    155162 
    156 def pinhole_resolution(q_calc, q, q_width): 
    157     """ 
     163def pinhole_resolution(q_calc, q, q_width, nsigma=PINHOLE_N_SIGMA): 
     164    r""" 
    158165    Compute the convolution matrix *W* for pinhole resolution 1-D data. 
    159166 
    … …  
    162169    *W*, the resolution smearing can be computed using *dot(W,q)*. 
    163170 
     171    Note that resolution is limited to $\pm 2.5 \sigma$.[1]  The true resolution 
     172    function is a broadened triangle, and does not extend over the entire 
     173    range $(-\infty, +\infty)$.  It is important to impose this limitation 
     174    since some models fall so steeply that the weighted value in gaussian 
     175    tails would otherwise dominate the integral. 
     176 
    164177    *q_calc* must be increasing.  *q_width* must be greater than zero. 
     178 
     179    [1] Barker, J. G., and J. S. Pedersen. 1995. Instrumental Smearing Effects 
     180    in Radially Symmetric Small-Angle Neutron Scattering by Numerical and 
     181    Analytical Methods. Journal of Applied Crystallography 28 (2): 105--14. 
     182    https://doi.org/10.1107/S0021889894010095. 
    165183    """ 
    166184    # The current algorithm is a midpoint rectangle rule.  In the test case, 
    … …  
    170188    cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :]) 
    171189    weights = cdf[1:] - cdf[:-1] 
     190    # Limit q range to +/- 2.5 sigma 
     191    qhigh = q + nsigma*q_width 
     192    #qlow = q - nsigma*q_width  # linear limits 
     193    qlow = q*q/qhigh  # log limits 
     194    weights[q_calc[:, None] < qlow[None, :]] = 0. 
     195    weights[q_calc[:, None] > qhigh[None, :]] = 0. 
    172196    weights /= np.sum(weights, axis=0)[None, :] 
    173197    return weights 
    … …  
    494518 
    495519 
    496 def gaussian(q, q0, dq): 
    497     """ 
    498     Return the Gaussian resolution function. 
     520def gaussian(q, q0, dq, nsigma=2.5): 
     521    """ 
     522    Return the truncated Gaussian resolution function. 
    499523 
    500524    *q0* is the center, *dq* is the width and *q* are the points to evaluate. 
    501525    """ 
    502     return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq) 
     526    # Calculate the density of the tails; the resulting gaussian needs to be 
     527    # scaled by this amount in ordere to integrate to 1.0 
     528    two_tail_density = 2 * (1 + erf(-nsigma/sqrt(2)))/2 
     529    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq)/(1-two_tail_density) 
    503530 
    504531 
    … …  
    558585 
    559586 
    560 def romberg_pinhole_1d(q, q_width, form, pars, nsigma=5): 
     587def romberg_pinhole_1d(q, q_width, form, pars, nsigma=2.5): 
    561588    """ 
    562589    Romberg integration for pinhole resolution. 
    … …  
    678705        np.testing.assert_equal(output, self.y) 
    679706 
     707    # TODO: turn pinhole/slit demos into tests 
     708 
     709    @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    680710    def test_pinhole(self): 
    681711        """ 
    … …  
    686716        theory = 12.0-1000.0*resolution.q_calc 
    687717        output = resolution.apply(theory) 
     718        # Note: answer came from output of previous run.  Non-integer 
     719        # values at ends come from the fact that q_calc does not 
     720        # extend beyond q, and so the weights don't balance. 
    688721        answer = [ 
    689             10.44785079, 9.84991299, 8.98101708, 
    690             7.99906585, 6.99998311, 6.00001689, 
    691             5.00093415, 4.01898292, 3.15008701, 2.55214921, 
     722            10.47037734, 9.86925860, 
     723            9., 8., 7., 6., 5., 4., 
     724            3.13074140, 2.52962266, 
    692725            ] 
    693726        np.testing.assert_allclose(output, answer, atol=1e-8) 
    … …  
    732765        self._compare(q, output, answer, 1e-6) 
    733766 
     767    @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    734768    def test_pinhole(self): 
    735769        """ 
    … …  
    746780        self._compare(q, output, answer, 3e-4) 
    747781 
     782    @unittest.skip("suppress comparison with old version; pinhole calc changed") 
    748783    def test_pinhole_romberg(self): 
    749784        """ 
    … …  
    761796        #                     2*np.pi/pars['radius']/200) 
    762797        #tol = 0.001 
    763         ## The default 3 sigma and no extra points gets 1% 
     798        ## The default 2.5 sigma and no extra points gets 1% 
    764799        q_calc = None  # type: np.ndarray 
    765800        tol = 0.01 
    … …  
    10801115 
    10811116    if isinstance(resolution, Slit1D): 
    1082         width, height = resolution.dqx, resolution.dqy 
     1117        width, height = resolution.qx_width, resolution.qy_width 
    10831118        Iq_romb = romberg_slit_1d(resolution.q, width, height, model, pars) 
    10841119    else: 
Note: See TracChangeset for help on using the changeset viewer.