1 | #!/usr/bin/env python |
---|
2 | r""" |
---|
3 | Multiple scattering calculator |
---|
4 | |
---|
5 | Calculate multiple scattering using 2D FFT convolution. |
---|
6 | |
---|
7 | Usage: |
---|
8 | |
---|
9 | -p, --probability: the scattering probability |
---|
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 | |
---|
22 | Let the scattering probability for $n$ scattering event at $q$ be $f_n(q)$, |
---|
23 | where |
---|
24 | .. math:: f_1(q) = \frac{I_1(q)}{\int I_1(q) dq} |
---|
25 | for $I_1(q)$, the single scattering from the system. After two scattering |
---|
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 |
---|
38 | .. math:: |
---|
39 | \frac{\sum_{k=1}^n \frac{P(k; \lambda)}{1 - P(0; \lambda)} < C |
---|
40 | |
---|
41 | Using the convolution theorem, where |
---|
42 | $F = \mathcal{F}(f)$ is the Fourier transform, |
---|
43 | .. math:: f * g = \mathcal{F}^{-1}\{\mathcal{F}\{f\} \cdot \mathcal{F}\{g\}\} |
---|
44 | so |
---|
45 | .. math:: f * \ldots * f = \mathcal{F}^{-1}\{ F^n \} |
---|
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:: |
---|
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\} |
---|
58 | |
---|
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. |
---|
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. |
---|
69 | """ |
---|
70 | |
---|
71 | from __future__ import print_function, division |
---|
72 | |
---|
73 | import argparse |
---|
74 | import time |
---|
75 | import os.path |
---|
76 | |
---|
77 | import numpy as np |
---|
78 | from numpy import pi |
---|
79 | from scipy.special import gamma |
---|
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 |
---|
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): |
---|
115 | #t0 = time.time() |
---|
116 | Iq = np.asarray(Iq, self.dtype) |
---|
117 | result = np.fft.fft2(Iq) |
---|
118 | #print("fft time", time.time()-t0) |
---|
119 | return result |
---|
120 | |
---|
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("numpy 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 | """ |
---|
174 | |
---|
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 |
---|
190 | else: |
---|
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("OpenCL 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!$. |
---|
251 | """ |
---|
252 | if transform is None: |
---|
253 | nx, ny = Iq.shape |
---|
254 | transform = Calculator(dims=(nx*2, ny*2), dtype=dtype) |
---|
255 | scale = np.sum(Iq) |
---|
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 |
---|
261 | |
---|
262 | def scattering_coeffs(p, coverage=0.99): |
---|
263 | L = -np.log(1-p) |
---|
264 | num_scatter = truncated_poisson_invcdf(coverage, L) |
---|
265 | coeffs = [L**k/gamma(k+2) for k in range(num_scatter)] |
---|
266 | return coeffs |
---|
267 | |
---|
268 | def truncated_poisson_invcdf(p, L): |
---|
269 | r""" |
---|
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 | |
---|
283 | def _forward_shift(Iq, dtype=PRECISION): |
---|
284 | # Prepare padded array and forward transform |
---|
285 | nq = Iq.shape[0] |
---|
286 | half_nq = nq//2 |
---|
287 | frame = np.zeros((2*nq, 2*nq), dtype=dtype) |
---|
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] |
---|
292 | return frame |
---|
293 | |
---|
294 | def _inverse_shift(frame, dtype=PRECISION): |
---|
295 | # Invert the transform and recover the transformed data |
---|
296 | nq = frame.shape[0]//2 |
---|
297 | half_nq = nq//2 |
---|
298 | Iq = np.empty((nq, nq), dtype=dtype) |
---|
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 |
---|
304 | |
---|
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 | |
---|
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, |
---|
524 | #'show_pars': True, |
---|
525 | 'show_pars': False, |
---|
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 | ) |
---|
537 | parser.add_argument('-p', '--probability', type=float, default=0.1, help="scattering probability") |
---|
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') |
---|
544 | parser.add_argument('-o', '--outfile', type=str, default="", help='random pars with random seed') |
---|
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') |
---|
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) |
---|
552 | res = MultipleScattering(qmax=opts.qmax, nq=opts.nq, window=opts.window, |
---|
553 | probability=opts.probability, is2d=opts.is2d) |
---|
554 | kernel = model.make_kernel(res.q_calc) |
---|
555 | #print(pars) |
---|
556 | bg = pars.get('background', 0.0) |
---|
557 | pars['background'] = 0.0 |
---|
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.): |
---|
563 | import pylab |
---|
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) |
---|
594 | else: |
---|
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) |
---|
617 | pylab.show() |
---|
618 | |
---|
619 | def plotxy(q, Iq): |
---|
620 | import pylab |
---|
621 | # q is a tuple of (q,) or (qx, qy) |
---|
622 | if len(q) == 1: |
---|
623 | pylab.loglog(q[0], Iq) |
---|
624 | else: |
---|
625 | data = Iq.copy() |
---|
626 | data[Iq <= 0] = np.min(Iq[Iq > 0])/2 |
---|
627 | pylab.imshow(np.log10(data)) |
---|
628 | |
---|
629 | if __name__ == "__main__": |
---|
630 | main() |
---|