[ae91fce] | 1 | #!/usr/bin/env python |
---|
| 2 | r""" |
---|
| 3 | Multiple scattering calculator |
---|
| 4 | |
---|
| 5 | Calculate multiple scattering using 2D FFT convolution. |
---|
| 6 | |
---|
| 7 | Usage: |
---|
| 8 | |
---|
[a09d55d] | 9 | -p, --probability: the scattering probability |
---|
[ae91fce] | 10 | -q, --qmax: that max q that you care about |
---|
| 11 | -w, --window: the extension window (q is calculated for qmax*window) |
---|
| 12 | -n, --nq: the number of mesh points (dq = qmax*window/nq) |
---|
| 13 | -r, --random: generate a random parameter set |
---|
| 14 | -2, --2d: perform the calculation for an oriented pattern |
---|
| 15 | model_name |
---|
| 16 | model_par=value ... |
---|
| 17 | |
---|
| 18 | Assume the probability of scattering is $p$. After each scattering event, |
---|
| 19 | $1-p$ neutrons will leave the system and go to the detector, and the remaining |
---|
| 20 | $p$ will scatter again. |
---|
| 21 | |
---|
[a09d55d] | 22 | Let the scattering probability for $n$ scattering event at $q$ be $f_n(q)$, |
---|
[ae91fce] | 23 | where |
---|
[a09d55d] | 24 | .. math:: f_1(q) = \frac{I_1(q)}{\int I_1(q) dq} |
---|
[ae91fce] | 25 | for $I_1(q)$, the single scattering from the system. After two scattering |
---|
[a09d55d] | 26 | events, the scattering probability will be the convolution of the first |
---|
| 27 | scattering and itself, or $f_2(q) = (f_1*f_1)(q)$. After $n$ events it will be |
---|
| 28 | $f_n(q) = (f_1 * \cdots * f_1)(q)$. The total scattering is calculated |
---|
| 29 | as the weighted sum of $f_k$, with weights following the Poisson distribution |
---|
| 30 | .. math:: P(k; \lambda) = \frac{\lambda^k e^{-\lambda}}{k!} |
---|
| 31 | for $\lambda$ determined by the total thickness divided by the mean |
---|
| 32 | free path between scattering, giving |
---|
| 33 | .. math:: |
---|
| 34 | I(q) = \sum_{k=0}^\infty P(k; \lambda) f_k(q) |
---|
| 35 | The $k=0$ term is ignored since it involves no scattering. |
---|
| 36 | We cut the series when cumulative probability is less than cutoff $C=99\%$, |
---|
| 37 | which is $\max n$ such that |
---|
[ae91fce] | 38 | .. math:: |
---|
[a09d55d] | 39 | \frac{\sum_{k=1}^n \frac{P(k; \lambda)}{1 - P(0; \lambda)} < C |
---|
[ae91fce] | 40 | |
---|
| 41 | Using the convolution theorem, where |
---|
| 42 | $F = \mathcal{F}(f)$ is the Fourier transform, |
---|
[a09d55d] | 43 | .. math:: f * g = \mathcal{F}^{-1}\{\mathcal{F}\{f\} \cdot \mathcal{F}\{g\}\} |
---|
[ae91fce] | 44 | so |
---|
[a09d55d] | 45 | .. math:: f * \ldots * f = \mathcal{F}^{-1}\{ F^n \} |
---|
[ae91fce] | 46 | Since the Fourier transform is a linear operator, we can move the polynomial |
---|
| 47 | expression for the convolution into the transform, giving |
---|
| 48 | .. math:: |
---|
[a09d55d] | 49 | I(q) = \mathcal{F}^{-1}\left\{ \sum_{k=1}^{n} P(k; \lambda) F^k \right\} |
---|
| 50 | In the dilute limit $L \rightarrow 0$ only the $k=1$ term is active, |
---|
| 51 | and so |
---|
| 52 | .. math:: |
---|
| 53 | P(1; \lambda) = \lambda e{-\lambda} = \int I_1(q) dq |
---|
| 54 | therefore we compute |
---|
| 55 | .. math:: |
---|
| 56 | I(q) = \int I_1(q) dq \mathcal{F}^{-1}\left\{ |
---|
| 57 | \sum_{l=1}^{n} \frac{P(k; \lambda)}{P(1; \lambda))} F^k \right\} |
---|
[ae91fce] | 58 | |
---|
[a09d55d] | 59 | For speed we may use the fast fourier transform with a power of two. |
---|
| 60 | The resulting $I(q)$ will be linearly spaced and likely heavily oversampled. |
---|
| 61 | The usual pinhole or slit resolution calculation can performed from these |
---|
| 62 | calculated values. |
---|
[bc248f8] | 63 | |
---|
| 64 | By default single precision OpenCL is used for the calculation. Set the |
---|
| 65 | environment variable *SAS_OPENCL=none* to use double precision numpy FFT |
---|
| 66 | instead. The OpenCL versions is about 10x faster on an elderly Mac with |
---|
| 67 | Intel HD 4000 graphics. The single precision numerical artifacts don't |
---|
| 68 | seem to seriously impact overall accuracy, though they look pretty bad. |
---|
[ae91fce] | 69 | """ |
---|
| 70 | |
---|
[a09d55d] | 71 | from __future__ import print_function, division |
---|
[ae91fce] | 72 | |
---|
| 73 | import argparse |
---|
| 74 | import time |
---|
[bc248f8] | 75 | import os.path |
---|
[ae91fce] | 76 | |
---|
| 77 | import numpy as np |
---|
[bc248f8] | 78 | from numpy import pi |
---|
[a09d55d] | 79 | from scipy.special import gamma |
---|
[ae91fce] | 80 | |
---|
| 81 | from sasmodels import core |
---|
| 82 | from sasmodels import compare |
---|
| 83 | from sasmodels import resolution2d |
---|
| 84 | from sasmodels.resolution import Resolution, bin_edges |
---|
| 85 | from sasmodels.data import empty_data1D, empty_data2D, plot_data |
---|
| 86 | from sasmodels.direct_model import call_kernel |
---|
[bc248f8] | 87 | import sasmodels.kernelcl |
---|
| 88 | |
---|
| 89 | # TODO: select fast and accurate fft library |
---|
| 90 | # clFFT: https://github.com/clMathLibraries/clFFT (AMD's OpenCL) |
---|
| 91 | # - gpyfft: https://github.com/geggo/gpyfft (wraps clFFT) |
---|
| 92 | # - arrayfire: https://github.com/arrayfire (wraps clFFT and much more; +cuda) |
---|
| 93 | # pyFFT: https://github.com/fjarri-attic/pyfft (based on Apple's OpenCL; +cuda) |
---|
| 94 | # - Reikna: https://github.com/fjarri/reikna (evolved from pyfft) |
---|
| 95 | # genFFT: https://software.intel.com/en-us/articles/genFFT (Intel's OpenCL) |
---|
| 96 | # VexCL: https://vexcl.readthedocs.io/en/latest/ (c++ library) |
---|
| 97 | # TODO: switch to fftw when opencl is not available |
---|
| 98 | |
---|
| 99 | try: |
---|
| 100 | import pyfft.cl |
---|
| 101 | import pyopencl.array as cl_array |
---|
| 102 | HAVE_OPENCL = sasmodels.kernelcl.use_opencl() |
---|
| 103 | except ImportError: |
---|
| 104 | HAVE_OPENCL = False |
---|
| 105 | PRECISION = np.dtype('f' if HAVE_OPENCL else 'd') # 'f' or 'd' |
---|
| 106 | USE_FAST = True # OpenCL faster, less accurate math |
---|
| 107 | |
---|
| 108 | class NumpyCalculator: |
---|
| 109 | def __init__(self, dims=None, dtype=PRECISION): |
---|
| 110 | self.dtype = dtype |
---|
| 111 | self.complex_dtype = np.dtype('F') if dtype == np.dtype('f') else np.dtype('D') |
---|
| 112 | pass |
---|
| 113 | |
---|
| 114 | def fft(self, Iq): |
---|
[a09d55d] | 115 | #t0 = time.time() |
---|
[bc248f8] | 116 | Iq = np.asarray(Iq, self.dtype) |
---|
| 117 | result = np.fft.fft2(Iq) |
---|
| 118 | #print("fft time", time.time()-t0) |
---|
| 119 | return result |
---|
[a09d55d] | 120 | |
---|
[bc248f8] | 121 | def ifft(self, fourier_frame): |
---|
| 122 | #t0 = time.time() |
---|
| 123 | fourier_frame = np.asarray(fourier_frame, self.complex_dtype) |
---|
| 124 | result = np.fft.ifft2(fourier_frame) |
---|
| 125 | #print("ifft time", time.time()-t0) |
---|
| 126 | return result |
---|
| 127 | |
---|
| 128 | def multiple_scattering(self, Iq, p, coverage=0.99): |
---|
| 129 | r""" |
---|
| 130 | Compute multiple scattering for I(q) given scattering probability p. |
---|
| 131 | |
---|
| 132 | Given a probability p of scattering with the thickness, the expected |
---|
| 133 | number of scattering events, $\lambda$ is $-\log(1 - p)$, giving a |
---|
| 134 | Poisson weighted sum of single, double, triple, etc. scattering patterns. |
---|
| 135 | The number of patterns used is based on coverage (default 99%). |
---|
| 136 | """ |
---|
| 137 | #t0 = time.time() |
---|
| 138 | coeffs = scattering_coeffs(p, coverage) |
---|
| 139 | poly = np.asarray(coeffs[::-1], dtype=self.dtype) |
---|
| 140 | scale = np.sum(Iq) |
---|
| 141 | frame = _forward_shift(Iq/scale, dtype=self.dtype) |
---|
| 142 | F = np.fft.fft2(frame) |
---|
| 143 | F_convolved = F * np.polyval(poly, F) |
---|
| 144 | frame = np.fft.ifft2(F_convolved) |
---|
| 145 | result = scale * _inverse_shift(frame.real, dtype=self.dtype) |
---|
| 146 | #print("multiscat time", time.time()-t0) |
---|
| 147 | return result |
---|
| 148 | |
---|
| 149 | # polyval1(c, x) computes (...((c0 x + c1) x + c2) x ... + cn) x |
---|
| 150 | # where c is an array of length *degree* and x is an array of |
---|
| 151 | # complex values (type double2) of length *n*. 2-D arrays can of |
---|
| 152 | # course be treated as 1-D arrays of length *nx* X *ny*. |
---|
| 153 | # When compiling with sasmodels.kernelcl.compile_model the double precision |
---|
| 154 | # types are converted to single precision as needed. See the code in |
---|
| 155 | # sasmodels.generate._convert_type for details. |
---|
| 156 | POLYVAL1_KERNEL = """ |
---|
| 157 | kernel void polyval1( |
---|
| 158 | const int degree, |
---|
| 159 | global const double *coeff, |
---|
| 160 | const int n, |
---|
| 161 | global double2 *array) |
---|
| 162 | { |
---|
| 163 | int index = get_global_id(0); |
---|
| 164 | if (index < n) { |
---|
| 165 | const double2 x = array[index]; |
---|
| 166 | double2 total = coeff[0]; |
---|
| 167 | for (int k=1; k < degree; k++) { |
---|
| 168 | total = fma(total, x, coeff[k]); |
---|
| 169 | } |
---|
| 170 | array[index] = total * x; |
---|
| 171 | } |
---|
| 172 | } |
---|
| 173 | """ |
---|
[ae91fce] | 174 | |
---|
[bc248f8] | 175 | class OpenclCalculator(NumpyCalculator): |
---|
| 176 | polyval1f = None |
---|
| 177 | polyval1d = None |
---|
| 178 | def __init__(self, dims, dtype=PRECISION): |
---|
| 179 | env = sasmodels.kernelcl.environment() |
---|
| 180 | context = env.get_context(dtype) |
---|
| 181 | if dtype == np.dtype('f'): |
---|
| 182 | if self.polyval1f is None: |
---|
| 183 | program = sasmodels.kernelcl.compile_model( |
---|
| 184 | context, POLYVAL1_KERNEL, dtype, fast=USE_FAST) |
---|
| 185 | # Assume context is always the same for a given dtype |
---|
| 186 | OpenclCalculator.polyval1f = program.polyval1 |
---|
| 187 | self.dtype = dtype |
---|
| 188 | self.complex_dtype = np.dtype('F') |
---|
| 189 | self.polyval1 = self.polyval1f |
---|
[ae91fce] | 190 | else: |
---|
[bc248f8] | 191 | if self.polyval1d is None: |
---|
| 192 | program = sasmodels.kernelcl.compile_model( |
---|
| 193 | context, POLYVAL1_KERNEL, dtype, fast=False) |
---|
| 194 | # Assume context is always the same for a given dtype |
---|
| 195 | OpenclCalculator.polyval1d = program.polyval1 |
---|
| 196 | self.dtype = dtype |
---|
| 197 | self.complex_type = np.dtype('D') |
---|
| 198 | self.polyval1 = self.polyval1d |
---|
| 199 | self.queue = env.get_queue(dtype) |
---|
| 200 | self.plan = pyfft.cl.Plan(dims, queue=self.queue) |
---|
| 201 | |
---|
| 202 | def fft(self, Iq): |
---|
| 203 | # forward transform |
---|
| 204 | #t0 = time.time() |
---|
| 205 | data = np.asarray(Iq, self.complex_dtype) |
---|
| 206 | gpu_data = cl_array.to_device(self.queue, data) |
---|
| 207 | self.plan.execute(gpu_data.data) |
---|
| 208 | result = gpu_data.get() |
---|
| 209 | #print("fft time", time.time()-t0) |
---|
| 210 | return result |
---|
| 211 | |
---|
| 212 | def ifft(self, fourier_frame): |
---|
| 213 | # inverse transform |
---|
| 214 | #t0 = time.time() |
---|
| 215 | data = np.asarray(fourier_frame, self.complex_dtype) |
---|
| 216 | gpu_data = cl_array.to_device(self.queue, data) |
---|
| 217 | self.plan.execute(gpu_data.data, inverse=True) |
---|
| 218 | result = gpu_data.get() |
---|
| 219 | #print("ifft time", time.time()-t0) |
---|
| 220 | return result |
---|
| 221 | |
---|
| 222 | def multiple_scattering(self, Iq, p, coverage=0.99): |
---|
| 223 | #t0 = time.time() |
---|
| 224 | coeffs = scattering_coeffs(p, coverage) |
---|
| 225 | scale = np.sum(Iq) |
---|
| 226 | poly = np.asarray(coeffs[::-1], self.dtype) |
---|
| 227 | frame = _forward_shift(Iq/scale, dtype=self.complex_dtype) |
---|
| 228 | gpu_data = cl_array.to_device(self.queue, frame) |
---|
| 229 | gpu_poly = cl_array.to_device(self.queue, poly) |
---|
| 230 | self.plan.execute(gpu_data.data) |
---|
| 231 | degree, n = poly.shape[0], frame.shape[0]*frame.shape[1] |
---|
| 232 | self.polyval1( |
---|
| 233 | self.queue, [n], None, |
---|
| 234 | np.int32(degree), gpu_poly.data, np.int32(n), gpu_data.data) |
---|
| 235 | self.plan.execute(gpu_data.data, inverse=True) |
---|
| 236 | frame = gpu_data.get() |
---|
| 237 | #result = scale * _inverse_shift(frame.real, dtype=self.dtype) |
---|
| 238 | result = scale * _inverse_shift(frame.real, dtype=self.dtype) |
---|
| 239 | #print("multiscat time", time.time()-t0) |
---|
| 240 | return result |
---|
| 241 | |
---|
| 242 | Calculator = OpenclCalculator if HAVE_OPENCL else NumpyCalculator |
---|
| 243 | |
---|
| 244 | def scattering_powers(Iq, n, dtype='f', transform=None): |
---|
| 245 | r""" |
---|
| 246 | Calculate the scattering powers up to n. |
---|
| 247 | |
---|
| 248 | This includes 1 even though it should just be Iq itself |
---|
| 249 | |
---|
| 250 | The frames are unweighted; to weight scale by $\lambda^k e^{-\lambda}/k!$. |
---|
[54a5024] | 251 | """ |
---|
[bc248f8] | 252 | if transform is None: |
---|
| 253 | nx, ny = Iq.shape |
---|
| 254 | transform = Calculator(dims=(nx*2, ny*2), dtype=dtype) |
---|
[a09d55d] | 255 | scale = np.sum(Iq) |
---|
[bc248f8] | 256 | frame = _forward_shift(Iq/scale, dtype=dtype) |
---|
| 257 | F = transform.fft(frame) |
---|
| 258 | powers = [scale * _inverse_shift(transform.ifft(F**(k+1)).real, dtype=dtype) |
---|
| 259 | for k in range(n)] |
---|
| 260 | return powers |
---|
[a09d55d] | 261 | |
---|
[bc248f8] | 262 | def scattering_coeffs(p, coverage=0.99): |
---|
[a09d55d] | 263 | L = -np.log(1-p) |
---|
| 264 | num_scatter = truncated_poisson_invcdf(coverage, L) |
---|
[bc248f8] | 265 | coeffs = [L**k/gamma(k+2) for k in range(num_scatter)] |
---|
| 266 | return coeffs |
---|
[a09d55d] | 267 | |
---|
| 268 | def truncated_poisson_invcdf(p, L): |
---|
[bc248f8] | 269 | r""" |
---|
[a09d55d] | 270 | Return smallest k such that cdf(k; L) > p from the truncated Poisson |
---|
| 271 | probability excluding k=0 |
---|
| 272 | """ |
---|
| 273 | # pmf(k; L) = L**k * exp(-L) / (k! * (1-exp(-L)) |
---|
| 274 | cdf = 0 |
---|
| 275 | pmf = -np.exp(-L) / np.expm1(-L) |
---|
| 276 | k = 0 |
---|
| 277 | while cdf < p: |
---|
| 278 | k += 1 |
---|
| 279 | pmf *= L/k |
---|
| 280 | cdf += pmf |
---|
| 281 | return k |
---|
| 282 | |
---|
[bc248f8] | 283 | def _forward_shift(Iq, dtype=PRECISION): |
---|
[a09d55d] | 284 | # Prepare padded array and forward transform |
---|
| 285 | nq = Iq.shape[0] |
---|
[ae91fce] | 286 | half_nq = nq//2 |
---|
[bc248f8] | 287 | frame = np.zeros((2*nq, 2*nq), dtype=dtype) |
---|
[a09d55d] | 288 | frame[:half_nq, :half_nq] = Iq[half_nq:, half_nq:] |
---|
| 289 | frame[-half_nq:, :half_nq] = Iq[:half_nq, half_nq:] |
---|
| 290 | frame[:half_nq, -half_nq:] = Iq[half_nq:, :half_nq] |
---|
| 291 | frame[-half_nq:, -half_nq:] = Iq[:half_nq, :half_nq] |
---|
[bc248f8] | 292 | return frame |
---|
[a09d55d] | 293 | |
---|
[bc248f8] | 294 | def _inverse_shift(frame, dtype=PRECISION): |
---|
[a09d55d] | 295 | # Invert the transform and recover the transformed data |
---|
[bc248f8] | 296 | nq = frame.shape[0]//2 |
---|
[ae91fce] | 297 | half_nq = nq//2 |
---|
[bc248f8] | 298 | Iq = np.empty((nq, nq), dtype=dtype) |
---|
[a09d55d] | 299 | Iq[half_nq:, half_nq:] = frame[:half_nq, :half_nq] |
---|
| 300 | Iq[:half_nq, half_nq:] = frame[-half_nq:, :half_nq] |
---|
| 301 | Iq[half_nq:, :half_nq] = frame[:half_nq, -half_nq:] |
---|
| 302 | Iq[:half_nq, :half_nq] = frame[-half_nq:, -half_nq:] |
---|
| 303 | return Iq |
---|
[ae91fce] | 304 | |
---|
[bc248f8] | 305 | |
---|
| 306 | class MultipleScattering(Resolution): |
---|
| 307 | def __init__(self, qmin=None, qmax=None, nq=None, window=2, |
---|
| 308 | probability=None, coverage=0.99, |
---|
| 309 | is2d=False, resolution=None, |
---|
| 310 | dtype=PRECISION): |
---|
| 311 | r""" |
---|
| 312 | Compute multiple scattering using Fourier convolution. |
---|
| 313 | |
---|
| 314 | The fourier steps are determined by *qmax*, the maximum $q$ value |
---|
| 315 | desired, *nq* the number of $q$ steps and *window*, the amount |
---|
| 316 | of padding around the circular convolution. The $q$ spacing |
---|
| 317 | will be $\Delta q = 2 q_\mathrm{max} w / n_q$. If *nq* is not |
---|
| 318 | given it will use $n_q = 2^k$ such that $\Delta q < q_\mathrm{min}$. |
---|
| 319 | |
---|
| 320 | *probability* is related to the expected number of scattering |
---|
| 321 | events in the sample $\lambda$ as $p = 1 = e^{-\lambda}$. As a |
---|
| 322 | hack to allow probability to be a fitted parameter, the "value" |
---|
| 323 | can be a function that takes no parameters and returns the current |
---|
| 324 | value of the probability. *coverage* determines how many scattering |
---|
| 325 | steps to consider. The default is 0.99, which sets $n$ such that |
---|
| 326 | $1 \ldots n$ covers 99% of the Poisson probability mass function. |
---|
| 327 | |
---|
| 328 | *is2d* is True then 2D scattering is used, otherwise it accepts |
---|
| 329 | and returns 1D scattering. |
---|
| 330 | |
---|
| 331 | *resolution* is the resolution function to apply after multiple |
---|
| 332 | scattering. If present, then the resolution $q$ vectors will provide |
---|
| 333 | default values for *qmin*, *qmax* and *nq*. |
---|
| 334 | """ |
---|
| 335 | # Infer qmin, qmax from instrument resolution calculator, if present |
---|
| 336 | if resolution is not None: |
---|
| 337 | is2d = hasattr(resolution, 'qx_data') |
---|
| 338 | if is2d: |
---|
| 339 | # 2D data |
---|
| 340 | if qmax is None: |
---|
| 341 | qx_calc, qy_calc = resolution.q_calc |
---|
| 342 | qmax = np.sqrt(np.max(qx_calc**2 + qy_calc**2)) |
---|
| 343 | if qmin is None and nq is None: |
---|
| 344 | qx, qy = resolution.data.x_bins, resolution.data.y_bins |
---|
| 345 | if qx and qy: |
---|
| 346 | dx = (np.max(qx) - np.min(qx)) / len(qx) |
---|
| 347 | dy = (np.max(qy) - np.min(qy)) / len(qy) |
---|
| 348 | else: |
---|
| 349 | qx, qy = resolution.data.qx_data, resolution.data.qy_data |
---|
| 350 | steps = np.sqrt(len(qx)) |
---|
| 351 | dx = (np.max(qx) - np.min(qx)) / steps |
---|
| 352 | dy = (np.max(qy) - np.min(qy)) / steps |
---|
| 353 | qmin = min(dx, dy) |
---|
| 354 | else: |
---|
| 355 | # 1D data |
---|
| 356 | if qmax is None: |
---|
| 357 | qmax = np.max(resolution.q_calc) |
---|
| 358 | if qmin is None and nq is None: |
---|
| 359 | qmin = np.min(np.abs(resolution.q_calc)) |
---|
| 360 | |
---|
| 361 | # estimate nq from qmin, qmax if not given explicitly |
---|
| 362 | q_range = qmax * window |
---|
| 363 | if nq is None: |
---|
| 364 | nq = 2**np.ceil(np.log2(q_range/qmin)) |
---|
| 365 | nq = int(nq) |
---|
| 366 | # Compute available qmin based on nq |
---|
| 367 | qmin = 2*q_range / nq |
---|
| 368 | #print(nq) |
---|
| 369 | |
---|
| 370 | # remember input parameters |
---|
| 371 | self.qmax = qmax |
---|
| 372 | self.qmin = qmin |
---|
| 373 | self.nq = nq |
---|
| 374 | self.probability = probability |
---|
| 375 | self.coverage = coverage |
---|
| 376 | self.is2d = is2d |
---|
| 377 | self.window = window |
---|
| 378 | self.resolution = resolution |
---|
| 379 | |
---|
| 380 | # Determine the q values to calculate |
---|
| 381 | q = np.linspace(-q_range, q_range, nq) |
---|
| 382 | qx, qy = np.meshgrid(q, q) |
---|
| 383 | if is2d: |
---|
| 384 | q_calc = (qx.flatten(), qy.flatten()) |
---|
| 385 | else: |
---|
| 386 | # For 1-D patterns, compute q from the center to the corners and |
---|
| 387 | # interpolate from there into the individual pixels. Given that |
---|
| 388 | # nq represents the points in [-qmax*windows, qmax*window], |
---|
| 389 | # then using 2*sqrt(2)*nq/2 will oversample the points in q by |
---|
| 390 | # a factor of two relative to the pixels. |
---|
| 391 | q_range_to_corner = np.sqrt(2.) * q_range |
---|
| 392 | nq_to_corner = 10*int(np.ceil(np.sqrt(2.) * nq)) |
---|
| 393 | q_to_corner = np.linspace(0, q_range_to_corner, nq_to_corner+1)[1:] |
---|
| 394 | q_calc = (q_to_corner,) |
---|
| 395 | # Remember the q radii of the calculated points |
---|
| 396 | self._radius = np.sqrt(qx**2 + qy**2) |
---|
| 397 | #self._q = q_to_corner |
---|
| 398 | self._q_steps = q |
---|
| 399 | self.q_calc = q_calc |
---|
| 400 | |
---|
| 401 | # TODO: use cleaner data representation than that from sasview |
---|
| 402 | # Resolution function forwards underlying q data (for plotting, etc?) |
---|
| 403 | if is2d: |
---|
| 404 | if resolution is not None: |
---|
| 405 | # forward resolution function info to multiscattering |
---|
| 406 | self.qx_data = resolution.qx_data |
---|
| 407 | self.qy_data = resolution.qy_data |
---|
| 408 | else: |
---|
| 409 | # no underlying resolution function, but make it look like there is |
---|
| 410 | self.qx_data, self.qy_data = q_calc |
---|
| 411 | else: |
---|
| 412 | # 1-D radial profile is determined by the q values we need to |
---|
| 413 | # compute, either for the calculated q values for the resolution |
---|
| 414 | # function (if any) or for the raw q values desired |
---|
| 415 | self._q = np.linspace(qmin, qmax, nq//(2*window)) |
---|
| 416 | self._edges = bin_edges(self._q) |
---|
| 417 | self._norm, _ = np.histogram(self._radius, bins=self._edges) |
---|
| 418 | if resolution is not None: |
---|
| 419 | self.q = resolution.q |
---|
| 420 | else: |
---|
| 421 | # no underlying resolution function, but make it look like there is |
---|
| 422 | self.q = self._q |
---|
| 423 | |
---|
| 424 | # Prepare the multiple scattering calculator (either numpy or OpenCL) |
---|
| 425 | self.transform = Calculator((2*nq, 2*nq), dtype=dtype) |
---|
| 426 | |
---|
| 427 | def apply(self, theory): |
---|
| 428 | if self.is2d: |
---|
| 429 | Iq_calc = theory |
---|
| 430 | else: |
---|
| 431 | Iq_calc = np.interp(self._radius, self.q_calc[0], theory) |
---|
| 432 | Iq_calc = Iq_calc.reshape(self.nq, self.nq) |
---|
| 433 | |
---|
| 434 | probability = self.probability() if callable(self.probability) else self.probability |
---|
| 435 | coverage = self.coverage |
---|
| 436 | #t0 = time.time() |
---|
| 437 | Iqxy = self.transform.multiple_scattering(Iq_calc, probability, coverage) |
---|
| 438 | #print("multiple scattering calc time", time.time()-t0) |
---|
| 439 | |
---|
| 440 | # remember the intermediate result in case we want to see it later |
---|
| 441 | self.Iqxy = Iqxy |
---|
| 442 | |
---|
| 443 | if self.is2d: |
---|
| 444 | if self.resolution is not None: |
---|
| 445 | Iqxy = self.resolution.apply(Iqxy) |
---|
| 446 | return Iqxy |
---|
| 447 | else: |
---|
| 448 | # remember the intermediate result in case we want to see it later |
---|
| 449 | Iq = self.radial_profile(Iqxy) |
---|
| 450 | self.Iq = Iq |
---|
| 451 | if self.resolution is not None: |
---|
| 452 | q = self._q |
---|
| 453 | Iq_res = np.interp(np.abs(self.resolution.q_calc), q, self.Iq) |
---|
| 454 | _ = """ |
---|
| 455 | k = 6 |
---|
| 456 | print("q theory", q[:k]) |
---|
| 457 | print("Iq theory", theory[:k]) |
---|
| 458 | print("interp NaN?", np.any(np.isnan(Iq_calc))) |
---|
| 459 | print("convolved NaN?", np.any(np.isnan(Iqxy))) |
---|
| 460 | print("Iq intgrated", self.Iq[:k]) |
---|
| 461 | print("radius", self._radius[self.nq/2,self.nq/2:self.nq/2+k]) |
---|
| 462 | print("q edges", self._edges[:k+1]) |
---|
| 463 | print("Iq norm", self._norm[:k]) |
---|
| 464 | print("res q", self.resolution.q_calc[:k]) |
---|
| 465 | print("Iq res", Iq_res[:k]) |
---|
| 466 | #print(Iq) |
---|
| 467 | #print(Iq_res) |
---|
| 468 | """ |
---|
| 469 | Iq = self.resolution.apply(Iq_res) |
---|
| 470 | return Iq |
---|
| 471 | |
---|
| 472 | def radial_profile(self, Iqxy): |
---|
| 473 | # circular average, no anti-aliasing |
---|
| 474 | Iq = np.histogram(self._radius, bins=self._edges, weights=Iqxy)[0]/self._norm |
---|
| 475 | return Iq |
---|
| 476 | |
---|
| 477 | |
---|
| 478 | def annular_average(qxy, Iqxy, qbins): |
---|
| 479 | """ |
---|
| 480 | Compute annular average of points at |
---|
| 481 | """ |
---|
| 482 | qxy, Iqxy = qxy.flatten(), Iqxy.flatten() |
---|
| 483 | index = np.argsort(qxy) |
---|
| 484 | qxy, Iqxy = qxy[index], Iqxy[index] |
---|
| 485 | print(qxy.shape, Iqxy.shape, index.shape, qbins.shape) |
---|
| 486 | #values = rebin(np.vstack((0., qxy)), Iqxy, qbins) |
---|
| 487 | integral = np.cumsum(Iqxy) |
---|
| 488 | Io = np.diff(np.interp(qbins, qxy, integral, left=0.)) |
---|
| 489 | # normalize by area of annulus |
---|
| 490 | # TODO: doesn't properly account for box. |
---|
| 491 | # Need to chop off chords that are greater than the box edges |
---|
| 492 | # (left, right, top and bottom), then add back in the corners |
---|
| 493 | # chopped off by both. https://en.wikipedia.org/wiki/Circular_segment |
---|
| 494 | norms = np.diff(pi*qbins**2) |
---|
| 495 | return Io/norms |
---|
| 496 | |
---|
| 497 | def rebin(x, I, xo): |
---|
| 498 | """ |
---|
| 499 | Rebin from edges *x*, bins I into edges *xo*. |
---|
| 500 | |
---|
| 501 | *x* and *xo* should be monotonically increasing. |
---|
| 502 | |
---|
| 503 | If *x* has duplicate values, then all corresponding values at I(x) will |
---|
| 504 | be effectively summed into the same bin. If *xo* has duplicate values, |
---|
| 505 | the first bin will contain the entire contents and subsequent bins will |
---|
| 506 | contain zeros. |
---|
| 507 | """ |
---|
| 508 | integral = np.cumsum(I) |
---|
| 509 | Io = np.diff(np.interp(xo, x[1:], integral, left=0.)) |
---|
| 510 | return Io |
---|
| 511 | |
---|
| 512 | |
---|
[ae91fce] | 513 | def parse_pars(model, opts): |
---|
| 514 | # type: (ModelInfo, argparse.Namespace) -> Dict[str, float] |
---|
| 515 | |
---|
| 516 | seed = np.random.randint(1000000) if opts.random and opts.seed < 0 else opts.seed |
---|
| 517 | compare_opts = { |
---|
| 518 | 'info': (model.info, model.info), |
---|
| 519 | 'use_demo': False, |
---|
| 520 | 'seed': seed, |
---|
| 521 | 'mono': True, |
---|
| 522 | 'magnetic': False, |
---|
| 523 | 'values': opts.pars, |
---|
[bc248f8] | 524 | #'show_pars': True, |
---|
| 525 | 'show_pars': False, |
---|
[ae91fce] | 526 | 'is2d': opts.is2d, |
---|
| 527 | } |
---|
| 528 | pars, pars2 = compare.parse_pars(compare_opts) |
---|
| 529 | return pars |
---|
| 530 | |
---|
| 531 | |
---|
| 532 | def main(): |
---|
| 533 | parser = argparse.ArgumentParser( |
---|
| 534 | description="Compute multiple scattering", |
---|
| 535 | formatter_class=argparse.ArgumentDefaultsHelpFormatter, |
---|
| 536 | ) |
---|
[a09d55d] | 537 | parser.add_argument('-p', '--probability', type=float, default=0.1, help="scattering probability") |
---|
[ae91fce] | 538 | parser.add_argument('-n', '--nq', type=int, default=1024, help='number of mesh points') |
---|
| 539 | parser.add_argument('-q', '--qmax', type=float, default=0.5, help='max q') |
---|
| 540 | parser.add_argument('-w', '--window', type=float, default=2.0, help='q calc = q max * window') |
---|
| 541 | parser.add_argument('-2', '--2d', dest='is2d', action='store_true', help='oriented sample') |
---|
| 542 | parser.add_argument('-s', '--seed', default=-1, help='random pars with given seed') |
---|
| 543 | parser.add_argument('-r', '--random', action='store_true', help='random pars with random seed') |
---|
[a09d55d] | 544 | parser.add_argument('-o', '--outfile', type=str, default="", help='random pars with random seed') |
---|
[049e02d] | 545 | parser.add_argument('model', type=str, help='sas model name such as cylinder') |
---|
| 546 | parser.add_argument('pars', type=str, nargs='*', help='model parameters such as radius=30') |
---|
[ae91fce] | 547 | opts = parser.parse_args() |
---|
| 548 | assert opts.nq%2 == 0, "require even # points" |
---|
| 549 | |
---|
| 550 | model = core.load_model(opts.model) |
---|
| 551 | pars = parse_pars(model, opts) |
---|
[bc248f8] | 552 | res = MultipleScattering(qmax=opts.qmax, nq=opts.nq, window=opts.window, |
---|
| 553 | probability=opts.probability, is2d=opts.is2d) |
---|
[ae91fce] | 554 | kernel = model.make_kernel(res.q_calc) |
---|
| 555 | #print(pars) |
---|
| 556 | bg = pars.get('background', 0.0) |
---|
| 557 | pars['background'] = 0.0 |
---|
[bc248f8] | 558 | theory = call_kernel(kernel, pars) |
---|
| 559 | Iq = res.apply(theory) + bg |
---|
| 560 | plot_and_save_powers(res, theory, Iq, outfile=opts.outfile, background=bg) |
---|
| 561 | |
---|
| 562 | def plot_and_save_powers(res, theory, result, plot=True, outfile="", background=0.): |
---|
[ae91fce] | 563 | import pylab |
---|
[bc248f8] | 564 | probability, coverage = res.probability, res.coverage |
---|
| 565 | weights = scattering_coeffs(probability, coverage) |
---|
| 566 | |
---|
| 567 | # cribbed from MultipleScattering.apply |
---|
| 568 | if res.is2d: |
---|
| 569 | Iq_calc = theory |
---|
| 570 | else: |
---|
| 571 | Iq_calc = np.interp(res._radius, res.q_calc[0], theory) |
---|
| 572 | Iq_calc = Iq_calc.reshape(res.nq, res.nq) |
---|
| 573 | |
---|
| 574 | # Compute the scattering powers for 1, 2, ... n scattering events |
---|
| 575 | powers = scattering_powers(Iq_calc, len(weights)) |
---|
| 576 | |
---|
| 577 | #plotxy(Iqxy); import pylab; pylab.figure() |
---|
| 578 | if res.is2d: |
---|
| 579 | if outfile: |
---|
| 580 | data = np.vstack([Ipower.flatten() for Ipower in powers]).T |
---|
| 581 | np.savetxt(outfile + "_powers.txt", data) |
---|
| 582 | data = np.vstack(Iq_calc).T |
---|
| 583 | np.savetxt(outfile + ".txt", data) |
---|
| 584 | if plot: |
---|
| 585 | plotxy((res._q_steps, res._q_steps), Iq_calc) |
---|
| 586 | pylab.title("single scattering F") |
---|
| 587 | for k, v in enumerate(powers[1:]): |
---|
| 588 | pylab.figure() |
---|
| 589 | plotxy((res._q_steps, res._q_steps), v+background) |
---|
| 590 | pylab.title("multiple scattering F^%d" % (k+2)) |
---|
| 591 | pylab.figure() |
---|
| 592 | plotxy((res._q_steps, res._q_steps), res.Iqxy+background) |
---|
| 593 | pylab.title("total scattering for p=%g" % probability) |
---|
[ae91fce] | 594 | else: |
---|
[bc248f8] | 595 | q = res._q |
---|
| 596 | Iq_powers = [res.radial_profile(Iqxy) for Iqxy in powers] |
---|
| 597 | if outfile: |
---|
| 598 | data = np.vstack([q, theory, res.Iq]).T |
---|
| 599 | np.savetxt(outfile + ".txt", data) |
---|
| 600 | # circular average, no anti-aliasing for individual powers |
---|
| 601 | data = np.vstack([q] + Iq_powers).T |
---|
| 602 | np.savetxt(outfile + "_powers.txt", data) |
---|
| 603 | if plot: |
---|
| 604 | # Plot 2D pattern for total scattering |
---|
| 605 | plotxy((res._q_steps, res._q_steps), res.Iqxy+background) |
---|
| 606 | pylab.title("total scattering for p=%g" % probability) |
---|
| 607 | pylab.figure() |
---|
| 608 | |
---|
| 609 | # Plot 1D pattern for partial scattering |
---|
| 610 | pylab.loglog(q, res.Iq+background, label="total for p=%g"%probability) |
---|
| 611 | #new_annulus = annular_average(res._radius, res.Iqxy, res._edges) |
---|
| 612 | #pylab.loglog(q, new_annulus+background, label="new total for p=%g"%probability) |
---|
| 613 | for n, (w, Ipower) in enumerate(zip(weights, Iq_powers)): |
---|
| 614 | pylab.loglog(q, w*Ipower+background, label="scattering^%d"%(n+1)) |
---|
| 615 | pylab.legend() |
---|
| 616 | pylab.title('total scattering for p=%g' % probability) |
---|
[ae91fce] | 617 | pylab.show() |
---|
| 618 | |
---|
[bc248f8] | 619 | def plotxy(q, Iq): |
---|
[ae91fce] | 620 | import pylab |
---|
[bc248f8] | 621 | # q is a tuple of (q,) or (qx, qy) |
---|
| 622 | if len(q) == 1: |
---|
| 623 | pylab.loglog(q[0], Iq) |
---|
[ae91fce] | 624 | else: |
---|
[bc248f8] | 625 | data = Iq.copy() |
---|
| 626 | data[Iq <= 0] = np.min(Iq[Iq > 0])/2 |
---|
[ae91fce] | 627 | pylab.imshow(np.log10(data)) |
---|
| 628 | |
---|
| 629 | if __name__ == "__main__": |
---|
| 630 | main() |
---|