Changeset d9633b1 in sasmodels


Ignore:
Timestamp:
Mar 13, 2015 2:09:47 PM (9 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
db8756e
Parents:
af1d68c
Message:

Add tests for pinhole resolution calculator

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/resolution.py

    r49d1d42f rd9633b1  
     1""" 
     2Define the resolution functions for the data. 
     3 
     4This defines classes for 1D and 2D resolution calculations. 
     5""" 
    16from scipy.special import erf 
    27from numpy import sqrt 
     
    49 
    510SLIT_SMEAR_POINTS = 500 
    6  
    7 class MatrixResolution: 
    8     def apply(self, Iq): 
    9         return np.dot(Iq, self.resolution_matrix) 
    10  
    11 class Pinhole1D(MatrixResolution): 
    12     def __init__(self, q, q_width): 
     11MINIMUM_RESOLUTION = 1e-8 
     12# TODO: Q_EXTEND_STEPS is much bigger than necessary 
     13Q_EXTEND_STEPS = 30   # number of extra q points above and below 
     14 
     15 
     16class Resolution1D(object): 
     17    """ 
     18    Abstract base class defining a 1D resolution function. 
     19 
     20    *q* is the set of q values at which the data is measured. 
     21 
     22    *q_calc* is the set of q values at which the theory needs to be evaluated. 
     23    This may extend and interpolate the q values. 
     24 
     25    *apply* is the method to call with I(q_calc) to compute the resolution 
     26    smeared theory I(q). 
     27    """ 
     28    q = None 
     29    q_calc = None 
     30    def apply(self, Iq_calc): 
     31        """ 
     32        Smear *Iq_calc* by the resolution function, returning *Iq*. 
     33        """ 
     34        raise NotImplementedError("Subclass does not define the apply function") 
     35 
     36class Perfect1D(Resolution1D): 
     37    """ 
     38    Resolution function to use when there is no actual resolution smearing 
     39    to be applied.  It has the same interface as the other resolution 
     40    functions, but returns the identity function. 
     41    """ 
     42    def __init__(self, q): 
     43        self.q_calc = self.q = q 
     44 
     45    def apply(self, Iq_calc): 
     46        return Iq_calc 
     47 
     48class Pinhole1D(Resolution1D): 
     49    r""" 
     50    Pinhole aperture with q-dependent gaussian resolution. 
     51 
     52    *q* points at which the data is measured. 
     53 
     54    *dq* gaussian 1-sigma resolution at each data point. 
     55 
     56    *q_calc* is the list of points to calculate, or None if this should 
     57    be autogenerated from the *q,dq*. 
     58    """ 
     59    def __init__(self, q, q_width, q_calc=None): 
     60        #TODO: maybe add min_step=np.inf 
     61        #*min_step* is the minimum point spacing to use when computing the 
     62        #underlying model.  It should be on the order of 
     63        #$\tfrac{1}{10}\tfrac{2\pi}{d_\text{max}}$ to make sure that fringes 
     64        #are computed with sufficient density to avoid aliasing effects. 
     65 
     66        # Protect against calls with q_width=0.  The extend_q function will 
     67        # not extend the q if q_width is 0, but q_width must be non-zero when 
     68        # constructing the weight matrix to avoid division by zero errors. 
     69        # In practice this should never be needed, since resolution should 
     70        # default to Perfect1D if the pinhole geometry is not defined. 
    1371        self.q, self.q_width = q, q_width 
    14         self.q_calc = pinhole_extend_q(q,q_width) 
    15         self.resolution_matrix = \ 
    16             pinhole_resolution(self.q_calc, self.q, self.q_width) 
    17  
    18 class Slit1D(MatrixResolution): 
     72        self.q_calc = pinhole_extend_q(q, q_width) if q_calc is None else q_calc 
     73        self.weight_matrix = pinhole_resolution(self.q_calc, 
     74                self.q, np.maximum(q_width, MINIMUM_RESOLUTION)) 
     75 
     76    def apply(self, Iq_calc): 
     77        return apply_resolution_matrix(self.weight_matrix, Iq_calc) 
     78 
     79class Slit1D(Resolution1D): 
     80    """ 
     81    Slit aperture with a complicated resolution function. 
     82 
     83    *q* points at which the data is measured. 
     84 
     85    *qx_width* slit width 
     86 
     87    *qy_height* slit height 
     88    """ 
    1989    def __init__(self, q, qx_width, qy_width): 
    2090        if np.isscalar(qx_width): 
     
    2292        if np.isscalar(qy_width): 
    2393            qy_width = qy_width*np.ones(len(q)) 
    24         self.q, self.qx_width, self.qy_width = q, qx_width, qy_width 
     94        self.q, self.qx_width, self.qy_width = [ 
     95            v.flatten() for v in (q, qx_width, qy_width)] 
    2596        self.q_calc = slit_extend_q(q, qx_width, qy_width) 
    26         self.resolution_matrix = \ 
     97        self.weight_matrix = \ 
    2798            slit_resolution(self.q_calc, self.q, self.qx_width, self.qy_width) 
     99 
     100    def apply(self, Iq_calc): 
     101        return apply_resolution_matrix(self.weight_matrix, Iq_calc) 
     102 
     103 
     104def apply_resolution_matrix(weight_matrix, Iq_calc): 
     105    """ 
     106    Apply the resolution weight matrix to the computed theory function. 
     107    """ 
     108    #print "apply shapes",Iq_calc.shape, self.weight_matrix.shape 
     109    Iq = np.dot(Iq_calc[None,:], weight_matrix) 
     110    #print "result shape",Iq.shape 
     111    return Iq.flatten() 
    28112 
    29113def pinhole_resolution(q_calc, q, q_width): 
     
    38122    """ 
    39123    edges = bin_edges(q_calc) 
    40     edges[edges<0.] = 0. # clip edges below zero 
     124    edges[edges<0.0] = 0.0 # clip edges below zero 
     125    index = (q_width>0.0) # treat perfect resolution differently 
    41126    G = erf( (edges[:,None] - q[None,:]) / (sqrt(2.0)*q_width)[None,:] ) 
    42127    weights = G[1:] - G[:-1] 
    43     weights /= np.sum(weights, axis=1) 
     128    # w = np.sum(weights, axis=0); print "w",w.shape, w 
     129    weights /= np.sum(weights, axis=0)[None,:] 
    44130    return weights 
    45131 
    46132 
    47133def pinhole_extend_q(q, q_width): 
    48     return q 
     134    """ 
     135    Given *q* and *q_width*, find a set of sampling points *q_calc* so 
     136    that each point I(q) has sufficient support from the underlying 
     137    function. 
     138    """ 
     139    # If using min_step, you could do something like: 
     140    #     q_calc = np.arange(min, max, min_step) 
     141    #     index = np.searchsorted(q_calc, q) 
     142    #     q_calc[index] = q 
     143    # A refinement would be to assign q to the nearest neighbour.  A further 
     144    # refinement would be to guard multiple q points between q_calc points, 
     145    # either by removing duplicates or by only moving the values nearest the 
     146    # edges.  For really sparse samples, you will want to remove all remaining 
     147    # points that are not within three standard deviations of a measured q. 
     148    q_min, q_max = np.min(q), np.max(q) 
     149    extended_min, extended_max = np.min(q - 3*q_width), np.max(q + 3*q_width) 
     150    nbins_low, nbins_high = Q_EXTEND_STEPS, Q_EXTEND_STEPS 
     151    if extended_min < q_min: 
     152        q_low = np.linspace(extended_min, q_min, nbins_low+1)[:-1] 
     153        q_low = q_low[q_low>0.0] 
     154    else: 
     155        q_low = [] 
     156    if extended_max > q_max: 
     157        q_high = np.linspace(q_max, extended_max, nbins_high+1)[:-1] 
     158    else: 
     159        q_high = [] 
     160    return np.concatenate([q_low, q, q_high]) 
    49161 
    50162 
    51163def slit_resolution(q_calc, q, qx_width, qy_width): 
    52164    edges = bin_edges(q_calc) # Note: requires q > 0 
    53     edges[edges<0.] = 0.0 # clip edges below zero 
     165    edges[edges<0.0] = 0.0 # clip edges below zero 
    54166    qy_min, qy_max = 0.0, edges[-1] 
    55167 
     
    80192        qy_sq = qy**2 
    81193        weights += (sqrt(E_sq[1:]-qy_sq) - sqrt(qy_sq - E_sq[:-1]))*in_x 
    82     weights /= np.sum(weights, axis=1) 
     194    w = np.sum(weights, axis=1); print "w",w.shape, w 
     195    weights /= np.sum(weights, axis=1)[:,None] 
    83196    return weights 
    84197 
     
    89202 
    90203def bin_edges(x): 
     204    """ 
     205    Determine bin edges from bin centers, assuming that edges are centered 
     206    between the bins. 
     207 
     208    Note: this uses the arithmetic mean, which may not be appropriate for 
     209    log-scaled data. 
     210    """ 
    91211    if len(x) < 2 or (np.diff(x)<0).any(): 
    92212        raise ValueError("Expected bins to be an increasing set") 
     
    98218    return edges 
    99219 
     220############################################################################ 
     221# unit tests 
     222############################################################################ 
     223import unittest 
     224 
     225class ResolutionTest(unittest.TestCase): 
     226 
     227    def setUp(self): 
     228        self.x = 0.001*np.arange(1,11) 
     229        self.y = self.Iq(self.x) 
     230 
     231    def Iq(self, q): 
     232        "Linear function for resolution unit test" 
     233        return 12.0 - 1000.0*q 
     234 
     235    def test_perfect(self): 
     236        """ 
     237        Perfect resolution and no smearing. 
     238        """ 
     239        resolution = Perfect1D(self.x) 
     240        Iq_calc = self.Iq(resolution.q_calc) 
     241        output = resolution.apply(Iq_calc) 
     242        np.testing.assert_equal(output, self.y) 
     243 
     244    def test_slit_zero(self): 
     245        """ 
     246        Slit smearing with perfect resolution. 
     247        """ 
     248        resolution = Slit1D(self.x, 0., 0.) 
     249        Iq_calc = self.Iq(resolution.q_calc) 
     250        output = resolution.apply(Iq_calc) 
     251        np.testing.assert_equal(output, self.y) 
     252 
     253    @unittest.skip("slit smearing is still broken") 
     254    def test_slit(self): 
     255        """ 
     256        Slit smearing with height 0.005 
     257        """ 
     258        resolution = Slit1D(self.x, 0., 0.005) 
     259        Iq_calc = self.Iq(resolution.q_calc) 
     260        output = resolution.apply(Iq_calc) 
     261        # The following commented line was the correct output for even bins [see smearer.cpp for details] 
     262        #answer = [ 9.666,  9.056,  8.329,  7.494,  6.642,  5.721,  4.774,  3.824,  2.871, 2.   ] 
     263        answer = [ 9.0618,  8.6401,  8.1186,  7.1391,  6.1528,  5.5555,     4.5584,  3.5606,  2.5623, 2.    ] 
     264        np.testing.assert_allclose(output, answer, atol=1e-4) 
     265 
     266    def test_pinhole_zero(self): 
     267        """ 
     268        Pinhole smearing with perfect resolution 
     269        """ 
     270        resolution = Pinhole1D(self.x, 0.0*self.x) 
     271        Iq_calc = self.Iq(resolution.q_calc) 
     272        output = resolution.apply(Iq_calc) 
     273        np.testing.assert_equal(output, self.y) 
     274 
     275    def test_pinhole(self): 
     276        """ 
     277        Pinhole smearing with dQ = 0.001 [Note: not dQ/Q = 0.001] 
     278        """ 
     279        resolution = Pinhole1D(self.x, 0.001*np.ones_like(self.x), 
     280                               q_calc=self.x) 
     281        Iq_calc = 12.0-1000.0*resolution.q_calc 
     282        output = resolution.apply(Iq_calc) 
     283        answer = [ 10.44785079,   9.84991299,   8.98101708, 
     284                  7.99906585,   6.99998311,   6.00001689, 
     285                  5.00093415,   4.01898292,   3.15008701,   2.55214921] 
     286        np.testing.assert_allclose(output, answer, atol=1e-8) 
     287 
     288# Q, dQ, I(Q) for Igor default sphere model. 
     289# combines CMSphere5.txt and CMSphere5smearsphere.txt from sasview/tests 
     290# TODO: move test data into its own file? 
     291SPHERE_RESOLUTION_TEST_DATA = """\ 
     2920.001278 0.0002847 2538.41176383 
     2930.001562 0.0002905 2536.91820405 
     2940.001846 0.0002956 2535.13182479 
     2950.002130 0.0003017 2533.06217813 
     2960.002414 0.0003087 2530.70378586 
     2970.002698 0.0003165 2528.05024192 
     2980.002982 0.0003249 2525.10408349 
     2990.003266 0.0003340 2521.86667499 
     3000.003550 0.0003437 2518.33907750 
     3010.003834 0.0003539 2514.52246995 
     3020.004118 0.0003646 2510.41798319 
     3030.004402 0.0003757 2506.02690988 
     3040.004686 0.0003872 2501.35067884 
     3050.004970 0.0003990 2496.38678318 
     3060.005253 0.0004112 2491.16237596 
     3070.005537 0.0004237 2485.63911673 
     3080.005821 0.0004365 2479.83657083 
     3090.006105 0.0004495 2473.75676948 
     3100.006389 0.0004628 2467.40145990 
     3110.006673 0.0004762 2460.77293372 
     3120.006957 0.0004899 2453.86724627 
     3130.007241 0.0005037 2446.69623838 
     3140.007525 0.0005177 2439.25775219 
     3150.007809 0.0005318 2431.55421398 
     3160.008093 0.0005461 2423.58785521 
     3170.008377 0.0005605 2415.36158137 
     3180.008661 0.0005750 2406.87009473 
     3190.008945 0.0005896 2398.12841186 
     3200.009229 0.0006044 2389.13360806 
     3210.009513 0.0006192 2379.88958042 
     3220.009797 0.0006341 2370.39776774 
     3230.010080 0.0006491 2360.69528793 
     3240.010360 0.0006641 2350.85169027 
     3250.010650 0.0006793 2340.42023633 
     3260.010930 0.0006945 2330.11206013 
     3270.011220 0.0007097 2319.20109972 
     3280.011500 0.0007251 2308.43503981 
     3290.011780 0.0007404 2297.44820179 
     3300.012070 0.0007558 2285.83853677 
     3310.012350 0.0007713 2274.41290746 
     3320.012640 0.0007868 2262.36219581 
     3330.012920 0.0008024 2250.51169731 
     3340.013200 0.0008180 2238.45596231 
     3350.013490 0.0008336 2225.76495666 
     3360.013770 0.0008493 2213.29618391 
     3370.014060 0.0008650 2200.19110751 
     3380.014340 0.0008807 2187.34050325 
     3390.014620 0.0008965 2174.30529864 
     3400.014910 0.0009123 2160.61632548 
     3410.015190 0.0009281 2147.21038112 
     3420.015470 0.0009440 2133.62023580 
     3430.015760 0.0009598 2119.37907426 
     3440.016040 0.0009757 2105.45234903 
     3450.016330 0.0009916 2090.86319102 
     3460.016610 0.0010080 2076.60576032 
     3470.016890 0.0010240 2062.19214565 
     3480.017180 0.0010390 2047.10550219 
     3490.017460 0.0010550 2032.38715621 
     3500.017740 0.0010710 2017.52560123 
     3510.018030 0.0010880 2001.99124318 
     3520.018310 0.0011040 1986.84662060 
     3530.018600 0.0011200 1971.03389745 
     3540.018880 0.0011360 1955.61395119 
     3550.019160 0.0011520 1940.08291563 
     3560.019450 0.0011680 1923.87672225 
     3570.019730 0.0011840 1908.10656374 
     3580.020020 0.0012000 1891.66297192 
     3590.020300 0.0012160 1875.66789021 
     3600.020580 0.0012320 1859.56357196 
     3610.020870 0.0012490 1842.79468290 
     3620.021150 0.0012650 1826.50064489 
     3630.021430 0.0012810 1810.11533702 
     3640.021720 0.0012970 1793.06840882 
     3650.022000 0.0013130 1776.51153580 
     3660.022280 0.0013290 1759.87201249 
     3670.022570 0.0013460 1742.57354412 
     3680.022850 0.0013620 1725.79397319 
     3690.023140 0.0013780 1708.35831550 
     3700.023420 0.0013940 1691.45256069 
     3710.023700 0.0014110 1674.48561783 
     3720.023990 0.0014270 1656.86525366 
     3730.024270 0.0014430 1639.79847285 
     3740.024550 0.0014590 1622.68887088 
     3750.024840 0.0014760 1604.96421100 
     3760.025120 0.0014920 1587.85768129 
     3770.025410 0.0015080 1569.99297335 
     3780.025690 0.0015240 1552.84580279 
     3790.025970 0.0015410 1535.54074115 
     3800.026260 0.0015570 1517.75249337 
     3810.026540 0.0015730 1500.40115023 
     3820.026820 0.0015900 1483.03632237 
     3830.027110 0.0016060 1465.05942429 
     3840.027390 0.0016220 1447.67682181 
     3850.027670 0.0016390 1430.46495191 
     3860.027960 0.0016550 1412.49232282 
     3870.028240 0.0016710 1395.13182318 
     3880.028520 0.0016880 1377.93439837 
     3890.028810 0.0017040 1359.99528971 
     3900.029090 0.0017200 1342.67274512 
     3910.029370 0.0017370 1325.55375609 
     392""" 
     393 
     394class Pinhole1DTest(unittest.TestCase): 
     395 
     396    def setUp(self): 
     397        # sample q, dq, I(q) calculated by NIST Igor SANS package 
     398        self.data = np.loadtxt(SPHERE_RESOLUTION_TEST_DATA.split('\n')).T 
     399        self.pars = dict(scale=1.0, background=0.01, radius=60.0, 
     400                         solvent_sld=6.3, sld=1) 
     401 
     402    def Iq(self, q, dq): 
     403        from sasmodels import core 
     404        from sasmodels.models import sphere 
     405 
     406        model = core.load_model(sphere) 
     407        resolution = Pinhole1D(q, dq) 
     408        kernel = core.make_kernel(model, [resolution.q_calc]) 
     409        Iq_calc = core.call_kernel(kernel, self.pars) 
     410        result = resolution.apply(Iq_calc) 
     411        return result 
     412 
     413    def test_sphere(self): 
     414        """ 
     415        Compare pinhole resolution smearing with NIST Igor SANS 
     416        """ 
     417        q, dq, answer = self.data 
     418        output = self.Iq(q, dq) 
     419        np.testing.assert_allclose(output, answer, rtol=0.006) 
     420 
     421    #TODO: is all sas data sampled densely enough to support resolution calcs? 
     422    @unittest.skip("suppress sparse data test; not supported by current code") 
     423    def test_sphere_sparse(self): 
     424        """ 
     425        Compare pinhole resolution smearing with NIST Igor SANS on sparse data 
     426        """ 
     427        q, dq, answer = self.data[:, ::20]  # Take every nth point 
     428        output = self.Iq(q, dq) 
     429        np.testing.assert_allclose(output, answer, rtol=0.006) 
     430 
    100431 
    101432############################################################################ 
Note: See TracChangeset for help on using the changeset viewer.