[9f7fbd9] | 1 | """ |
---|
| 2 | Define the resolution functions for the data. |
---|
| 3 | |
---|
| 4 | This defines classes for 1D and 2D resolution calculations. |
---|
| 5 | """ |
---|
| 6 | from __future__ import division |
---|
| 7 | from scipy.special import erf |
---|
| 8 | from numpy import sqrt, log, log10 |
---|
| 9 | import numpy as np |
---|
| 10 | |
---|
| 11 | MINIMUM_RESOLUTION = 1e-8 |
---|
| 12 | |
---|
| 13 | class Resolution1D(object): |
---|
| 14 | """ |
---|
| 15 | Abstract base class defining a 1D resolution function. |
---|
| 16 | |
---|
| 17 | *q* is the set of q values at which the data is measured. |
---|
| 18 | |
---|
| 19 | *q_calc* is the set of q values at which the theory needs to be evaluated. |
---|
| 20 | This may extend and interpolate the q values. |
---|
| 21 | |
---|
| 22 | *apply* is the method to call with I(q_calc) to compute the resolution |
---|
| 23 | smeared theory I(q). |
---|
| 24 | """ |
---|
| 25 | q = None |
---|
| 26 | q_calc = None |
---|
| 27 | def apply(self, theory): |
---|
| 28 | """ |
---|
| 29 | Smear *theory* by the resolution function, returning *Iq*. |
---|
| 30 | """ |
---|
| 31 | raise NotImplementedError("Subclass does not define the apply function") |
---|
| 32 | |
---|
| 33 | |
---|
| 34 | class Perfect1D(Resolution1D): |
---|
| 35 | """ |
---|
| 36 | Resolution function to use when there is no actual resolution smearing |
---|
| 37 | to be applied. It has the same interface as the other resolution |
---|
| 38 | functions, but returns the identity function. |
---|
| 39 | """ |
---|
| 40 | def __init__(self, q): |
---|
| 41 | self.q_calc = self.q = q |
---|
| 42 | |
---|
| 43 | def apply(self, theory): |
---|
| 44 | return theory |
---|
| 45 | |
---|
| 46 | |
---|
| 47 | class Pinhole1D(Resolution1D): |
---|
| 48 | r""" |
---|
| 49 | Pinhole aperture with q-dependent gaussian resolution. |
---|
| 50 | |
---|
| 51 | *q* points at which the data is measured. |
---|
| 52 | |
---|
| 53 | *q_width* gaussian 1-sigma resolution at each data point. |
---|
| 54 | |
---|
| 55 | *q_calc* is the list of points to calculate, or None if this should |
---|
| 56 | be estimated from the *q* and *q_width*. |
---|
| 57 | """ |
---|
| 58 | def __init__(self, q, q_width, q_calc=None): |
---|
| 59 | #*min_step* is the minimum point spacing to use when computing the |
---|
| 60 | #underlying model. It should be on the order of |
---|
| 61 | #$\tfrac{1}{10}\tfrac{2\pi}{d_\text{max}}$ to make sure that fringes |
---|
| 62 | #are computed with sufficient density to avoid aliasing effects. |
---|
| 63 | |
---|
| 64 | # Protect against calls with q_width=0. The extend_q function will |
---|
| 65 | # not extend the q if q_width is 0, but q_width must be non-zero when |
---|
| 66 | # constructing the weight matrix to avoid division by zero errors. |
---|
| 67 | # In practice this should never be needed, since resolution should |
---|
| 68 | # default to Perfect1D if the pinhole geometry is not defined. |
---|
| 69 | self.q, self.q_width = q, q_width |
---|
| 70 | self.q_calc = pinhole_extend_q(q, q_width) \ |
---|
| 71 | if q_calc is None else np.sort(q_calc) |
---|
| 72 | self.weight_matrix = pinhole_resolution(self.q_calc, |
---|
| 73 | self.q, np.maximum(q_width, MINIMUM_RESOLUTION)) |
---|
| 74 | |
---|
| 75 | def apply(self, theory): |
---|
| 76 | return apply_resolution_matrix(self.weight_matrix, theory) |
---|
| 77 | |
---|
| 78 | |
---|
| 79 | class 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 | |
---|
| 89 | *q_calc* is the list of points to calculate, or None if this should |
---|
| 90 | be estimated from the *q* and *q_width*. |
---|
| 91 | |
---|
| 92 | The *weight_matrix* is computed by :func:`slit1d_resolution` |
---|
| 93 | """ |
---|
| 94 | def __init__(self, q, width, height, q_calc=None): |
---|
| 95 | # TODO: maybe issue warnings rather than raising errors |
---|
| 96 | if not np.isscalar(width): |
---|
| 97 | if np.any(np.diff(width) > 0.0): |
---|
| 98 | raise ValueError("Slit resolution requires fixed width slits") |
---|
| 99 | width = width[0] |
---|
| 100 | if not np.isscalar(height): |
---|
| 101 | if np.any(np.diff(height) > 0.0): |
---|
| 102 | raise ValueError("Slit resolution requires fixed height slits") |
---|
| 103 | height = height[0] |
---|
| 104 | |
---|
| 105 | # Remember what width/height was used even though we won't need them |
---|
| 106 | # after the weight matrix is constructed |
---|
| 107 | self.width, self.height = width, height |
---|
| 108 | |
---|
| 109 | self.q = q.flatten() |
---|
| 110 | self.q_calc = slit_extend_q(q, width, height) \ |
---|
| 111 | if q_calc is None else np.sort(q_calc) |
---|
| 112 | self.weight_matrix = \ |
---|
| 113 | slit_resolution(self.q_calc, self.q, width, height) |
---|
| 114 | |
---|
| 115 | def apply(self, theory): |
---|
| 116 | return apply_resolution_matrix(self.weight_matrix, theory) |
---|
| 117 | |
---|
| 118 | |
---|
| 119 | def apply_resolution_matrix(weight_matrix, theory): |
---|
| 120 | """ |
---|
| 121 | Apply the resolution weight matrix to the computed theory function. |
---|
| 122 | """ |
---|
| 123 | #print "apply shapes", theory.shape, weight_matrix.shape |
---|
| 124 | Iq = np.dot(theory[None,:], weight_matrix) |
---|
| 125 | #print "result shape",Iq.shape |
---|
| 126 | return Iq.flatten() |
---|
| 127 | |
---|
| 128 | |
---|
| 129 | def pinhole_resolution(q_calc, q, q_width): |
---|
| 130 | """ |
---|
| 131 | Compute the convolution matrix *W* for pinhole resolution 1-D data. |
---|
| 132 | |
---|
| 133 | Each row *W[i]* determines the normalized weight that the corresponding |
---|
| 134 | points *q_calc* contribute to the resolution smeared point *q[i]*. Given |
---|
| 135 | *W*, the resolution smearing can be computed using *dot(W,q)*. |
---|
| 136 | |
---|
| 137 | *q_calc* must be increasing. *q_width* must be greater than zero. |
---|
| 138 | """ |
---|
| 139 | # The current algorithm is a midpoint rectangle rule. In the test case, |
---|
| 140 | # neither trapezoid nor Simpson's rule improved the accuracy. |
---|
| 141 | edges = bin_edges(q_calc) |
---|
| 142 | edges[edges<0.0] = 0.0 # clip edges below zero |
---|
| 143 | G = erf( (edges[:,None] - q[None,:]) / (sqrt(2.0)*q_width)[None,:] ) |
---|
| 144 | weights = G[1:] - G[:-1] |
---|
| 145 | weights /= np.sum(weights, axis=0)[None,:] |
---|
| 146 | return weights |
---|
| 147 | |
---|
| 148 | |
---|
| 149 | def slit_resolution(q_calc, q, width, height): |
---|
| 150 | r""" |
---|
| 151 | Build a weight matrix to compute *I_s(q)* from *I(q_calc)*, given |
---|
| 152 | $q_v$ = *width* and $q_h$ = *height*. |
---|
| 153 | |
---|
| 154 | *width* and *height* are scalars since current instruments use the |
---|
| 155 | same slit settings for all measurement points. |
---|
| 156 | |
---|
| 157 | If slit height is large relative to width, use: |
---|
| 158 | |
---|
| 159 | .. math:: |
---|
| 160 | |
---|
| 161 | I_s(q_o) = \frac{1}{\Delta q_v} |
---|
| 162 | \int_0^{\Delta q_v} I(\sqrt{q_o^2 + u^2} du |
---|
| 163 | |
---|
| 164 | If slit width is large relative to height, use: |
---|
| 165 | |
---|
| 166 | .. math:: |
---|
| 167 | |
---|
| 168 | I_s(q_o) = \frac{1}{2 \Delta q_v} |
---|
| 169 | \int_{-\Delta q_v}^{\Delta q_v} I(u) du |
---|
| 170 | """ |
---|
| 171 | if width == 0.0 and height == 0.0: |
---|
| 172 | #print "condition zero" |
---|
| 173 | return 1 |
---|
| 174 | |
---|
| 175 | q_edges = bin_edges(q_calc) # Note: requires q > 0 |
---|
| 176 | q_edges[q_edges<0.0] = 0.0 # clip edges below zero |
---|
| 177 | |
---|
| 178 | #np.set_printoptions(linewidth=10000) |
---|
| 179 | if width <= 100.0 * height or height == 0: |
---|
| 180 | # The current algorithm is a midpoint rectangle rule. In the test case, |
---|
| 181 | # neither trapezoid nor Simpson's rule improved the accuracy. |
---|
| 182 | #print "condition h", q_edges.shape, q.shape, q_calc.shape |
---|
| 183 | weights = np.zeros((len(q), len(q_calc)), 'd') |
---|
| 184 | for i, qi in enumerate(q): |
---|
| 185 | weights[i, :] = np.diff(q_to_u(q_edges, qi)) |
---|
| 186 | weights /= width |
---|
| 187 | weights = weights.T |
---|
| 188 | else: |
---|
| 189 | #print "condition w" |
---|
| 190 | # Make q_calc into a row vector, and q into a column vector |
---|
| 191 | q, q_calc = q[None,:], q_calc[:,None] |
---|
| 192 | in_x = (q_calc >= q-width) * (q_calc <= q+width) |
---|
| 193 | weights = np.diff(q_edges)[:,None] * in_x |
---|
| 194 | |
---|
| 195 | return weights |
---|
| 196 | |
---|
| 197 | |
---|
| 198 | def pinhole_extend_q(q, q_width, nsigma=3): |
---|
| 199 | """ |
---|
| 200 | Given *q* and *q_width*, find a set of sampling points *q_calc* so |
---|
| 201 | that each point I(q) has sufficient support from the underlying |
---|
| 202 | function. |
---|
| 203 | """ |
---|
| 204 | q_min, q_max = np.min(q - nsigma*q_width), np.max(q + nsigma*q_width) |
---|
| 205 | return geometric_extrapolation(q, q_min, q_max) |
---|
| 206 | |
---|
| 207 | |
---|
| 208 | def slit_extend_q(q, width, height): |
---|
| 209 | """ |
---|
| 210 | Given *q*, *width* and *height*, find a set of sampling points *q_calc* so |
---|
| 211 | that each point I(q) has sufficient support from the underlying |
---|
| 212 | function. |
---|
| 213 | """ |
---|
| 214 | height # keep lint happy |
---|
| 215 | q_min, q_max = np.min(q), np.max(np.sqrt(q**2 + width**2)) |
---|
| 216 | return geometric_extrapolation(q, q_min, q_max) |
---|
| 217 | |
---|
| 218 | |
---|
| 219 | def bin_edges(x): |
---|
| 220 | """ |
---|
| 221 | Determine bin edges from bin centers, assuming that edges are centered |
---|
| 222 | between the bins. |
---|
| 223 | |
---|
| 224 | Note: this uses the arithmetic mean, which may not be appropriate for |
---|
| 225 | log-scaled data. |
---|
| 226 | """ |
---|
| 227 | if len(x) < 2 or (np.diff(x)<0).any(): |
---|
| 228 | raise ValueError("Expected bins to be an increasing set") |
---|
| 229 | edges = np.hstack([ |
---|
| 230 | x[0] - 0.5*(x[1] - x[0]), # first point minus half first interval |
---|
| 231 | 0.5*(x[1:] + x[:-1]), # mid points of all central intervals |
---|
| 232 | x[-1] + 0.5*(x[-1] - x[-2]), # last point plus half last interval |
---|
| 233 | ]) |
---|
| 234 | return edges |
---|
| 235 | |
---|
| 236 | |
---|
| 237 | def q_to_u(q, q0): |
---|
| 238 | """ |
---|
| 239 | Convert *q* values to *u* values for the integral computed at *q0*. |
---|
| 240 | """ |
---|
| 241 | # array([value])**2 - value**2 is not always zero |
---|
| 242 | qpsq = q**2 - q0**2 |
---|
| 243 | qpsq[qpsq<0] = 0 |
---|
| 244 | return sqrt(qpsq) |
---|
| 245 | |
---|
| 246 | |
---|
| 247 | def interpolate(q, max_step): |
---|
| 248 | """ |
---|
| 249 | Returns *q_calc* with points spaced at most max_step apart. |
---|
| 250 | """ |
---|
| 251 | step = np.diff(q) |
---|
| 252 | index = step>max_step |
---|
| 253 | if np.any(index): |
---|
| 254 | inserts = [] |
---|
| 255 | for q_i,step_i in zip(q[:-1][index],step[index]): |
---|
| 256 | n = np.ceil(step_i/max_step) |
---|
| 257 | inserts.extend(q_i + np.arange(1,n)*(step_i/n)) |
---|
| 258 | # Extend a couple of fringes beyond the end of the data |
---|
| 259 | inserts.extend(q[-1] + np.arange(1,8)*max_step) |
---|
| 260 | q_calc = np.sort(np.hstack((q,inserts))) |
---|
| 261 | else: |
---|
| 262 | q_calc = q |
---|
| 263 | return q_calc |
---|
| 264 | |
---|
| 265 | |
---|
| 266 | def linear_extrapolation(q, q_min, q_max): |
---|
| 267 | """ |
---|
| 268 | Extrapolate *q* out to [*q_min*, *q_max*] using the step size in *q* as |
---|
| 269 | a guide. Extrapolation below uses steps the same size as the first |
---|
| 270 | interval. Extrapolation above uses steps the same size as the final |
---|
| 271 | interval. |
---|
| 272 | |
---|
| 273 | if *q_min* is zero or less then *q[0]/10* is used instead. |
---|
| 274 | """ |
---|
| 275 | q = np.sort(q) |
---|
| 276 | if q_min < q[0]: |
---|
| 277 | if q_min <= 0: q_min = q[0]/10 |
---|
| 278 | delta = q[1] - q[0] |
---|
| 279 | q_low = np.arange(q_min, q[0], delta) |
---|
| 280 | else: |
---|
| 281 | q_low = [] |
---|
| 282 | if q_max > q[-1]: |
---|
| 283 | delta = q[-1] - q[-2] |
---|
| 284 | q_high = np.arange(q[-1]+delta, q_max, delta) |
---|
| 285 | else: |
---|
| 286 | q_high = [] |
---|
| 287 | return np.concatenate([q_low, q, q_high]) |
---|
| 288 | |
---|
| 289 | |
---|
| 290 | def geometric_extrapolation(q, q_min, q_max): |
---|
| 291 | r""" |
---|
| 292 | Extrapolate *q* to [*q_min*, *q_max*] using geometric steps, with the |
---|
| 293 | average geometric step size in *q* as the step size. |
---|
| 294 | |
---|
| 295 | if *q_min* is zero or less then *q[0]/10* is used instead. |
---|
| 296 | |
---|
| 297 | Starting at $q_1$ and stepping geometrically by $\Delta q$ |
---|
| 298 | to $q_n$ in $n$ points gives a geometric average of: |
---|
| 299 | |
---|
| 300 | .. math:: |
---|
| 301 | |
---|
| 302 | \log \Delta q = (\log q_n - log q_1) / (n - 1) |
---|
| 303 | |
---|
| 304 | From this we can compute the number of steps required to extend $q$ |
---|
| 305 | from $q_n$ to $q_\text{max}$ by $\Delta q$ as: |
---|
| 306 | |
---|
| 307 | .. math:: |
---|
| 308 | |
---|
| 309 | n_\text{extend} = (\log q_\text{max} - \log q_n) / \log \Delta q |
---|
| 310 | |
---|
| 311 | Substituting: |
---|
| 312 | |
---|
| 313 | n_\text{extend} = (n-1) (\log q_\text{max} - \log q_n) |
---|
| 314 | / (\log q_n - log q_1) |
---|
| 315 | """ |
---|
| 316 | q = np.sort(q) |
---|
| 317 | delta_q = (len(q)-1)/(log(q[-1]) - log(q[0])) |
---|
| 318 | if q_min < q[0]: |
---|
| 319 | if q_min < 0: q_min = q[0]/10 |
---|
| 320 | n_low = delta_q * (log(q[0])-log(q_min)) |
---|
| 321 | q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1] |
---|
| 322 | else: |
---|
| 323 | q_low = [] |
---|
| 324 | if q_max > q[-1]: |
---|
| 325 | n_high = delta_q * (log(q_max)-log(q[-1])) |
---|
| 326 | q_high = np.logspace(log10(q[-1]), log10(q_max), np.ceil(n_high)+1)[1:] |
---|
| 327 | else: |
---|
| 328 | q_high = [] |
---|
| 329 | return np.concatenate([q_low, q, q_high]) |
---|
| 330 | |
---|