[d9633b1] | 1 | """ |
---|
| 2 | Define the resolution functions for the data. |
---|
| 3 | |
---|
| 4 | This defines classes for 1D and 2D resolution calculations. |
---|
| 5 | """ |
---|
[3fdb4b6] | 6 | from scipy.special import erf |
---|
| 7 | from numpy import sqrt |
---|
| 8 | import numpy as np |
---|
| 9 | |
---|
[63b32bb] | 10 | SLIT_SMEAR_POINTS = 500 |
---|
[d9633b1] | 11 | MINIMUM_RESOLUTION = 1e-8 |
---|
| 12 | # TODO: Q_EXTEND_STEPS is much bigger than necessary |
---|
| 13 | Q_EXTEND_STEPS = 30 # number of extra q points above and below |
---|
| 14 | |
---|
| 15 | |
---|
| 16 | class 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 | |
---|
| 36 | class 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 | |
---|
| 48 | class Pinhole1D(Resolution1D): |
---|
| 49 | r""" |
---|
| 50 | Pinhole aperture with q-dependent gaussian resolution. |
---|
[63b32bb] | 51 | |
---|
[d9633b1] | 52 | *q* points at which the data is measured. |
---|
[49d1d42f] | 53 | |
---|
[d9633b1] | 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. |
---|
[49d1d42f] | 71 | self.q, self.q_width = q, q_width |
---|
[d9633b1] | 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 | |
---|
| 79 | class Slit1D(Resolution1D): |
---|
| 80 | """ |
---|
| 81 | Slit aperture with a complicated resolution function. |
---|
| 82 | |
---|
| 83 | *q* points at which the data is measured. |
---|
[49d1d42f] | 84 | |
---|
[d9633b1] | 85 | *qx_width* slit width |
---|
| 86 | |
---|
| 87 | *qy_height* slit height |
---|
| 88 | """ |
---|
[49d1d42f] | 89 | def __init__(self, q, qx_width, qy_width): |
---|
| 90 | if np.isscalar(qx_width): |
---|
| 91 | qx_width = qx_width*np.ones(len(q)) |
---|
| 92 | if np.isscalar(qy_width): |
---|
| 93 | qy_width = qy_width*np.ones(len(q)) |
---|
[d9633b1] | 94 | self.q, self.qx_width, self.qy_width = [ |
---|
| 95 | v.flatten() for v in (q, qx_width, qy_width)] |
---|
[49d1d42f] | 96 | self.q_calc = slit_extend_q(q, qx_width, qy_width) |
---|
[d9633b1] | 97 | self.weight_matrix = \ |
---|
[49d1d42f] | 98 | slit_resolution(self.q_calc, self.q, self.qx_width, self.qy_width) |
---|
| 99 | |
---|
[d9633b1] | 100 | def apply(self, Iq_calc): |
---|
| 101 | return apply_resolution_matrix(self.weight_matrix, Iq_calc) |
---|
| 102 | |
---|
| 103 | |
---|
| 104 | def 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() |
---|
| 112 | |
---|
[63b32bb] | 113 | def pinhole_resolution(q_calc, q, q_width): |
---|
[3fdb4b6] | 114 | """ |
---|
| 115 | Compute the convolution matrix *W* for pinhole resolution 1-D data. |
---|
| 116 | |
---|
| 117 | Each row *W[i]* determines the normalized weight that the corresponding |
---|
| 118 | points *q_calc* contribute to the resolution smeared point *q[i]*. Given |
---|
| 119 | *W*, the resolution smearing can be computed using *dot(W,q)*. |
---|
| 120 | |
---|
| 121 | *q_calc* must be increasing. |
---|
| 122 | """ |
---|
| 123 | edges = bin_edges(q_calc) |
---|
[d9633b1] | 124 | edges[edges<0.0] = 0.0 # clip edges below zero |
---|
| 125 | index = (q_width>0.0) # treat perfect resolution differently |
---|
[63b32bb] | 126 | G = erf( (edges[:,None] - q[None,:]) / (sqrt(2.0)*q_width)[None,:] ) |
---|
| 127 | weights = G[1:] - G[:-1] |
---|
[d9633b1] | 128 | # w = np.sum(weights, axis=0); print "w",w.shape, w |
---|
| 129 | weights /= np.sum(weights, axis=0)[None,:] |
---|
[3fdb4b6] | 130 | return weights |
---|
| 131 | |
---|
[49d1d42f] | 132 | |
---|
| 133 | def pinhole_extend_q(q, q_width): |
---|
[d9633b1] | 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]) |
---|
[49d1d42f] | 161 | |
---|
| 162 | |
---|
[3fdb4b6] | 163 | def slit_resolution(q_calc, q, qx_width, qy_width): |
---|
| 164 | edges = bin_edges(q_calc) # Note: requires q > 0 |
---|
[d9633b1] | 165 | edges[edges<0.0] = 0.0 # clip edges below zero |
---|
[63b32bb] | 166 | qy_min, qy_max = 0.0, edges[-1] |
---|
[3fdb4b6] | 167 | |
---|
[49d1d42f] | 168 | # Make q_calc into a row vector, and q, qx_width, qy_width into columns |
---|
| 169 | # Make weights match [ q_calc X q ] |
---|
[3fdb4b6] | 170 | weights = np.zeros((len(q),len(q_calc)),'d') |
---|
[49d1d42f] | 171 | q_calc = q_calc[None,:] |
---|
| 172 | q, qx_width, qy_width, edges = [ |
---|
| 173 | v[:,None] for v in (q, qx_width, qy_width, edges)] |
---|
| 174 | |
---|
[63b32bb] | 175 | # Loop for width (height is analytical). |
---|
[3fdb4b6] | 176 | # Condition: height >>> width, otherwise, below is not accurate enough. |
---|
[63b32bb] | 177 | # Smear weight numerical iteration for width>0 when height>0. |
---|
[3fdb4b6] | 178 | # When width = 0, the numerical iteration will be skipped. |
---|
| 179 | # The resolution calculation for the height is done by direct integration, |
---|
[63b32bb] | 180 | # assuming the I(q'=sqrt(q_j^2-(q+shift_w)^2)) is constant within |
---|
| 181 | # a q' bin, [q_high, q_low]. |
---|
| 182 | # In general, this weight numerical iteration for width>0 might be a rough |
---|
| 183 | # approximation, but it must be good enough when height >>> width. |
---|
[49d1d42f] | 184 | E_sq = edges**2 |
---|
[63b32bb] | 185 | y_points = SLIT_SMEAR_POINTS if np.any(qy_width>0) else 1 |
---|
| 186 | qy_step = 0 if y_points == 1 else qy_width/(y_points-1) |
---|
| 187 | for k in range(-y_points+1,y_points): |
---|
| 188 | qy = np.clip(q + qy_step*k, qy_min, qy_max) |
---|
[3fdb4b6] | 189 | qx_low = qy |
---|
| 190 | qx_high = sqrt(qx_low**2 + qx_width**2) |
---|
[49d1d42f] | 191 | in_x = (q_calc>=qx_low)*(q_calc<=qx_high) |
---|
| 192 | qy_sq = qy**2 |
---|
| 193 | weights += (sqrt(E_sq[1:]-qy_sq) - sqrt(qy_sq - E_sq[:-1]))*in_x |
---|
[d9633b1] | 194 | w = np.sum(weights, axis=1); print "w",w.shape, w |
---|
| 195 | weights /= np.sum(weights, axis=1)[:,None] |
---|
[63b32bb] | 196 | return weights |
---|
[3fdb4b6] | 197 | |
---|
[49d1d42f] | 198 | |
---|
| 199 | def slit_extend_q(q, qx_width, qy_width): |
---|
| 200 | return q |
---|
| 201 | |
---|
| 202 | |
---|
[3fdb4b6] | 203 | def bin_edges(x): |
---|
[d9633b1] | 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 | """ |
---|
[3fdb4b6] | 211 | if len(x) < 2 or (np.diff(x)<0).any(): |
---|
| 212 | raise ValueError("Expected bins to be an increasing set") |
---|
| 213 | edges = np.hstack([ |
---|
| 214 | x[0] - 0.5*(x[1] - x[0]), # first point minus half first interval |
---|
| 215 | 0.5*(x[1:] + x[:-1]), # mid points of all central intervals |
---|
| 216 | x[-1] + 0.5*(x[-1] - x[-2]), # last point plus half last interval |
---|
| 217 | ]) |
---|
| 218 | return edges |
---|
[49d1d42f] | 219 | |
---|
[d9633b1] | 220 | ############################################################################ |
---|
| 221 | # unit tests |
---|
| 222 | ############################################################################ |
---|
| 223 | import unittest |
---|
| 224 | |
---|
| 225 | class 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? |
---|
| 291 | SPHERE_RESOLUTION_TEST_DATA = """\ |
---|
| 292 | 0.001278 0.0002847 2538.41176383 |
---|
| 293 | 0.001562 0.0002905 2536.91820405 |
---|
| 294 | 0.001846 0.0002956 2535.13182479 |
---|
| 295 | 0.002130 0.0003017 2533.06217813 |
---|
| 296 | 0.002414 0.0003087 2530.70378586 |
---|
| 297 | 0.002698 0.0003165 2528.05024192 |
---|
| 298 | 0.002982 0.0003249 2525.10408349 |
---|
| 299 | 0.003266 0.0003340 2521.86667499 |
---|
| 300 | 0.003550 0.0003437 2518.33907750 |
---|
| 301 | 0.003834 0.0003539 2514.52246995 |
---|
| 302 | 0.004118 0.0003646 2510.41798319 |
---|
| 303 | 0.004402 0.0003757 2506.02690988 |
---|
| 304 | 0.004686 0.0003872 2501.35067884 |
---|
| 305 | 0.004970 0.0003990 2496.38678318 |
---|
| 306 | 0.005253 0.0004112 2491.16237596 |
---|
| 307 | 0.005537 0.0004237 2485.63911673 |
---|
| 308 | 0.005821 0.0004365 2479.83657083 |
---|
| 309 | 0.006105 0.0004495 2473.75676948 |
---|
| 310 | 0.006389 0.0004628 2467.40145990 |
---|
| 311 | 0.006673 0.0004762 2460.77293372 |
---|
| 312 | 0.006957 0.0004899 2453.86724627 |
---|
| 313 | 0.007241 0.0005037 2446.69623838 |
---|
| 314 | 0.007525 0.0005177 2439.25775219 |
---|
| 315 | 0.007809 0.0005318 2431.55421398 |
---|
| 316 | 0.008093 0.0005461 2423.58785521 |
---|
| 317 | 0.008377 0.0005605 2415.36158137 |
---|
| 318 | 0.008661 0.0005750 2406.87009473 |
---|
| 319 | 0.008945 0.0005896 2398.12841186 |
---|
| 320 | 0.009229 0.0006044 2389.13360806 |
---|
| 321 | 0.009513 0.0006192 2379.88958042 |
---|
| 322 | 0.009797 0.0006341 2370.39776774 |
---|
| 323 | 0.010080 0.0006491 2360.69528793 |
---|
| 324 | 0.010360 0.0006641 2350.85169027 |
---|
| 325 | 0.010650 0.0006793 2340.42023633 |
---|
| 326 | 0.010930 0.0006945 2330.11206013 |
---|
| 327 | 0.011220 0.0007097 2319.20109972 |
---|
| 328 | 0.011500 0.0007251 2308.43503981 |
---|
| 329 | 0.011780 0.0007404 2297.44820179 |
---|
| 330 | 0.012070 0.0007558 2285.83853677 |
---|
| 331 | 0.012350 0.0007713 2274.41290746 |
---|
| 332 | 0.012640 0.0007868 2262.36219581 |
---|
| 333 | 0.012920 0.0008024 2250.51169731 |
---|
| 334 | 0.013200 0.0008180 2238.45596231 |
---|
| 335 | 0.013490 0.0008336 2225.76495666 |
---|
| 336 | 0.013770 0.0008493 2213.29618391 |
---|
| 337 | 0.014060 0.0008650 2200.19110751 |
---|
| 338 | 0.014340 0.0008807 2187.34050325 |
---|
| 339 | 0.014620 0.0008965 2174.30529864 |
---|
| 340 | 0.014910 0.0009123 2160.61632548 |
---|
| 341 | 0.015190 0.0009281 2147.21038112 |
---|
| 342 | 0.015470 0.0009440 2133.62023580 |
---|
| 343 | 0.015760 0.0009598 2119.37907426 |
---|
| 344 | 0.016040 0.0009757 2105.45234903 |
---|
| 345 | 0.016330 0.0009916 2090.86319102 |
---|
| 346 | 0.016610 0.0010080 2076.60576032 |
---|
| 347 | 0.016890 0.0010240 2062.19214565 |
---|
| 348 | 0.017180 0.0010390 2047.10550219 |
---|
| 349 | 0.017460 0.0010550 2032.38715621 |
---|
| 350 | 0.017740 0.0010710 2017.52560123 |
---|
| 351 | 0.018030 0.0010880 2001.99124318 |
---|
| 352 | 0.018310 0.0011040 1986.84662060 |
---|
| 353 | 0.018600 0.0011200 1971.03389745 |
---|
| 354 | 0.018880 0.0011360 1955.61395119 |
---|
| 355 | 0.019160 0.0011520 1940.08291563 |
---|
| 356 | 0.019450 0.0011680 1923.87672225 |
---|
| 357 | 0.019730 0.0011840 1908.10656374 |
---|
| 358 | 0.020020 0.0012000 1891.66297192 |
---|
| 359 | 0.020300 0.0012160 1875.66789021 |
---|
| 360 | 0.020580 0.0012320 1859.56357196 |
---|
| 361 | 0.020870 0.0012490 1842.79468290 |
---|
| 362 | 0.021150 0.0012650 1826.50064489 |
---|
| 363 | 0.021430 0.0012810 1810.11533702 |
---|
| 364 | 0.021720 0.0012970 1793.06840882 |
---|
| 365 | 0.022000 0.0013130 1776.51153580 |
---|
| 366 | 0.022280 0.0013290 1759.87201249 |
---|
| 367 | 0.022570 0.0013460 1742.57354412 |
---|
| 368 | 0.022850 0.0013620 1725.79397319 |
---|
| 369 | 0.023140 0.0013780 1708.35831550 |
---|
| 370 | 0.023420 0.0013940 1691.45256069 |
---|
| 371 | 0.023700 0.0014110 1674.48561783 |
---|
| 372 | 0.023990 0.0014270 1656.86525366 |
---|
| 373 | 0.024270 0.0014430 1639.79847285 |
---|
| 374 | 0.024550 0.0014590 1622.68887088 |
---|
| 375 | 0.024840 0.0014760 1604.96421100 |
---|
| 376 | 0.025120 0.0014920 1587.85768129 |
---|
| 377 | 0.025410 0.0015080 1569.99297335 |
---|
| 378 | 0.025690 0.0015240 1552.84580279 |
---|
| 379 | 0.025970 0.0015410 1535.54074115 |
---|
| 380 | 0.026260 0.0015570 1517.75249337 |
---|
| 381 | 0.026540 0.0015730 1500.40115023 |
---|
| 382 | 0.026820 0.0015900 1483.03632237 |
---|
| 383 | 0.027110 0.0016060 1465.05942429 |
---|
| 384 | 0.027390 0.0016220 1447.67682181 |
---|
| 385 | 0.027670 0.0016390 1430.46495191 |
---|
| 386 | 0.027960 0.0016550 1412.49232282 |
---|
| 387 | 0.028240 0.0016710 1395.13182318 |
---|
| 388 | 0.028520 0.0016880 1377.93439837 |
---|
| 389 | 0.028810 0.0017040 1359.99528971 |
---|
| 390 | 0.029090 0.0017200 1342.67274512 |
---|
| 391 | 0.029370 0.0017370 1325.55375609 |
---|
| 392 | """ |
---|
| 393 | |
---|
| 394 | class 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 | |
---|
[49d1d42f] | 431 | |
---|
[db8756e] | 432 | def main(): |
---|
| 433 | """ |
---|
| 434 | Run tests given is sys.argv. |
---|
| 435 | |
---|
| 436 | Returns 0 if success or 1 if any tests fail. |
---|
| 437 | """ |
---|
| 438 | import sys |
---|
| 439 | import xmlrunner |
---|
| 440 | |
---|
| 441 | suite = unittest.TestSuite() |
---|
| 442 | suite.addTest(unittest.defaultTestLoader.loadTestsFromModule(sys.modules[__name__])) |
---|
| 443 | |
---|
| 444 | runner = xmlrunner.XMLTestRunner(output='logs') |
---|
| 445 | result = runner.run(suite) |
---|
| 446 | return 1 if result.failures or result.errors else 0 |
---|
| 447 | |
---|
| 448 | |
---|
[49d1d42f] | 449 | ############################################################################ |
---|
| 450 | # usage demo |
---|
| 451 | ############################################################################ |
---|
| 452 | |
---|
| 453 | def _eval_demo_1d(resolution, title): |
---|
| 454 | from sasmodels import core |
---|
| 455 | from sasmodels.models import cylinder |
---|
| 456 | ## or alternatively: |
---|
| 457 | # cylinder = core.load_model_definition('cylinder') |
---|
| 458 | model = core.load_model(cylinder) |
---|
| 459 | |
---|
| 460 | kernel = core.make_kernel(model, [resolution.q_calc]) |
---|
| 461 | Iq_calc = core.call_kernel(kernel, {'length':210, 'radius':500}) |
---|
| 462 | Iq = resolution.apply(Iq_calc) |
---|
| 463 | |
---|
| 464 | import matplotlib.pyplot as plt |
---|
| 465 | plt.loglog(resolution.q_calc, Iq_calc, label='unsmeared') |
---|
| 466 | plt.loglog(resolution.q, Iq, label='smeared', hold=True) |
---|
| 467 | plt.legend() |
---|
| 468 | plt.title(title) |
---|
| 469 | plt.xlabel("Q (1/Ang)") |
---|
| 470 | plt.ylabel("I(Q) (1/cm)") |
---|
| 471 | |
---|
| 472 | def demo_pinhole_1d(): |
---|
| 473 | q = np.logspace(-3,-1,400) |
---|
| 474 | dq = 0.1*q |
---|
| 475 | resolution = Pinhole1D(q, dq) |
---|
| 476 | _eval_demo_1d(resolution, title="10% dQ/Q Pinhole Resolution") |
---|
| 477 | |
---|
| 478 | def demo_slit_1d(): |
---|
| 479 | q = np.logspace(-3,-1,400) |
---|
| 480 | qx_width = 0.005 |
---|
| 481 | qy_width = 0.0 |
---|
| 482 | resolution = Slit1D(q, qx_width, qy_width) |
---|
| 483 | _eval_demo_1d(resolution, title="0.005 Qx Slit Resolution") |
---|
| 484 | |
---|
| 485 | def demo(): |
---|
| 486 | import matplotlib.pyplot as plt |
---|
| 487 | plt.subplot(121) |
---|
| 488 | demo_pinhole_1d() |
---|
| 489 | plt.subplot(122) |
---|
| 490 | demo_slit_1d() |
---|
| 491 | plt.show() |
---|
| 492 | |
---|
| 493 | |
---|
| 494 | if __name__ == "__main__": |
---|
[db8756e] | 495 | #demo() |
---|
| 496 | main() |
---|
[49d1d42f] | 497 | |
---|
| 498 | |
---|