# source:sasmodels/explore/multiscat.py@e309e23

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since e309e23 was e309e23, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

multiscat: improved plots

• Property mode set to 100644
File size: 11.6 KB
Line
1#!/usr/bin/env python
2r"""
3Multiple scattering calculator
4
5Calculate multiple scattering using 2D FFT convolution.
6
7Usage:
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
18Assume 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
22Let the scattering probability for $n$ scattering event at $q$ be $f_n(q)$,
23where
24.. math:: f_1(q) = \frac{I_1(q)}{\int I_1(q) dq}
25for $I_1(q)$, the single scattering from the system. After two scattering
26events, the scattering probability will be the convolution of the first
27scattering 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
29as the weighted sum of $f_k$, with weights following the Poisson distribution
30.. math:: P(k; \lambda) = \frac{\lambda^k e^{-\lambda}}{k!}
31for $\lambda$ determined by the total thickness divided by the mean
32free path between scattering, giving
33.. math::
34    I(q) = \sum_{k=0}^\infty P(k; \lambda) f_k(q)
35The $k=0$ term is ignored since it involves no scattering.
36We cut the series when cumulative probability is less than cutoff $C=99\%$,
37which is $\max n$ such that
38.. math::
39    \frac{\sum_{k=1}^n \frac{P(k; \lambda)}{1 - P(0; \lambda)} < C
40
41Using 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\}\}
44so
45.. math:: f * \ldots * f = \mathcal{F}^{-1}\{ F^n \}
46Since the Fourier transform is a linear operator, we can move the polynomial
47expression 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\}
50In the dilute limit $L \rightarrow 0$ only the $k=1$ term is active,
51and so
52.. math::
53    P(1; \lambda) = \lambda e{-\lambda} = \int I_1(q) dq
54therefore 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
59For speed we may use the fast fourier transform with a power of two.
60The resulting $I(q)$ will be linearly spaced and likely heavily oversampled.
61The usual pinhole or slit resolution calculation can performed from these
62calculated values.
63"""
64
65from __future__ import print_function, division
66
67import argparse
68import time
69
70import numpy as np
71from scipy.special import gamma
72
73from sasmodels import core
74from sasmodels import compare
75from sasmodels import resolution2d
76from sasmodels.resolution import Resolution, bin_edges
77from sasmodels.data import empty_data1D, empty_data2D, plot_data
78from sasmodels.direct_model import call_kernel
79
80
81class MultipleScattering(Resolution):
82    def __init__(self, qmax, nq, probability, is2d, window=2, power=0):
83        self.qmax = qmax
84        self.nq = nq
85        self.power = power
86        self.probability = probability
87        self.is2d = is2d
88        self.window = window
89
90        q_range = qmax * window
91        q = np.linspace(-q_range, q_range, nq)
92        qx, qy = np.meshgrid(q, q)
93
94        if is2d:
95            q_calc = [qx.flatten(), qy.flatten()]
96        else:
97            q_range_corners = np.sqrt(2.) * q_range
98            nq_corners = int(np.sqrt(2.) * nq/2)
99            q_corners = np.linspace(0, q_range_corners, nq_corners+1)[1:]
100            q_calc = [q_corners]
101            self._qxy = np.sqrt(qx**2 + qy**2)
102            self._edges = bin_edges(q_corners)
103            self._norm = np.histogram(self._qxy, bins=self._edges)[0]
104
105        self.q_calc = q_calc
106        self.q_range = q_range
107
108    def apply(self, theory, background=0.0, outfile="", plot=False):
109        #t0 = time.time()
110        if self.is2d:
111            Iq_calc = theory
112        else:
113            q_corners = self.q_calc[0]
114            Iq_calc = np.interp(self._qxy, q_corners, theory)
115        Iq_calc = Iq_calc.reshape(self.nq, self.nq)
116        #plotxy(Iq_calc); import pylab; pylab.figure()
117        if self.power > 0:
118            Iqxy = scattering_power(Iq_calc, self.power)
119            powers = []
120        else:
121            Iqxy, powers = multiple_scattering(
122                Iq_calc, self.probability,
123                coverage=0.99,
124                return_powers=(outfile!="") or plot,
125                )
126        #print("multiple scattering calc time", time.time()-t0)
127
128        #plotxy(Iqxy); import pylab; pylab.figure()
129        if self.is2d:
130            if outfile and powers:
131                data = np.vstack([Ipower.flatten() for Ipower in powers]).T
132                np.savetxt(outfile + "_powers.txt", data)
133            if outfile:
134                data = np.vstack(Iq_calc).T
135                np.savetxt(outfile + ".txt", data)
136            if plot:
137                import pylab
138                plotxy(Iq_calc)
139                pylab.title("single scattering")
140                pylab.figure()
141
142            return Iqxy + background
143        else:
144            # circular average, no anti-aliasing
145            Iq = np.histogram(self._qxy, bins=self._edges, weights=Iqxy)[0]/self._norm
146            if powers:
147                Iq_powers = [np.histogram(self._qxy, bins=self._edges, weights=Ipower)[0]/self._norm
148                             for Ipower in powers]#  or powers[:5] for the first five
149
150            if outfile:
151                data = np.vstack([q_corners, theory, Iq]).T
152                np.savetxt(outfile + ".txt", data)
153            if outfile and powers:
154                # circular average, no anti-aliasing for individual powers
155                data = np.vstack([q_corners] + Iq_powers).T
156                np.savetxt(outfile + "_powers.txt", data)
157            if plot:
158                import pylab
159                plotxy(Iqxy)
160                pylab.title("multiple scattering")
161                pylab.figure()
162            if plot and powers:
163                import pylab
164                L = -np.log(1-self.probability)
165                pylab.loglog(q_corners, Iq, label="total %g"%self.probability)
166                for n, Ipower in enumerate(Iq_powers):
167                    k = n+1
168                    w = L**(k-1)/gamma(k+1)
169                    pylab.loglog(q_corners, w*Ipower, label="scattering**%d"%k)
170                pylab.legend()
171                pylab.figure()
172
173            return q_corners, Iq + background
174
175def scattering_power(Iq, n):
176    """
177    Calculate the nth scattering power as a distribution.  To get the
178    weighted contribution, scale by $\lambda^k e^{-\lambda}/k!$.
179    """
180    scale = np.sum(Iq)
181    F = _forward_fft(Iq/scale)
182    result = _inverse_fft(F**n)
183    return result
184
185def multiple_scattering(Iq, p, coverage=0.99, return_powers=False):
186    """
187    Compute multiple scattering for I(q) given scattering probability p.
188
189    Given a probability p of scattering with the thickness, the expected
190    number of scattering events, $\lambda$ is $-\log(1 - p)$, giving a
191    Poisson weighted sum of single, double, triple, etc. scattering patterns.
192    The number of patterns used is based on coverage (default 99%).
193    """
194    L = -np.log(1-p)
195    num_scatter = truncated_poisson_invcdf(coverage, L)
196
197    # Compute multiple scattering via convolution.
198    scale = np.sum(Iq)
199    F = _forward_fft(Iq/scale)
200    coeffs = [L**(k-1)/gamma(k+1) for k in range(1, num_scatter+1)]
201    multiple_scattering = F * np.polyval(coeffs[::-1], F)
202    result = scale * _inverse_fft(multiple_scattering)
203
204    if return_powers:
205        powers = [scale * _inverse_fft(F**k) for k in range(1, num_scatter+1)]
206    else:
207        powers = []
208
209    return result, powers
210
211def truncated_poisson_invcdf(p, L):
212    """
213    Return smallest k such that cdf(k; L) > p from the truncated Poisson
214    probability excluding k=0
215    """
216    # pmf(k; L) = L**k * exp(-L) / (k! * (1-exp(-L))
217    cdf = 0
218    pmf = -np.exp(-L) / np.expm1(-L)
219    k = 0
220    while cdf < p:
221        k += 1
222        pmf *= L/k
223        cdf += pmf
224    return k
225
226def _forward_fft(Iq):
227    # Prepare padded array and forward transform
228    nq = Iq.shape[0]
229    half_nq = nq//2
230    frame = np.zeros((2*nq, 2*nq))
231    frame[:half_nq, :half_nq] = Iq[half_nq:, half_nq:]
232    frame[-half_nq:, :half_nq] = Iq[:half_nq, half_nq:]
233    frame[:half_nq, -half_nq:] = Iq[half_nq:, :half_nq]
234    frame[-half_nq:, -half_nq:] = Iq[:half_nq, :half_nq]
235    fourier_frame = np.fft.fft2(frame)
236    return fourier_frame
237
238def _inverse_fft(fourier_frame):
239    # Invert the transform and recover the transformed data
240    nq = fourier_frame.shape[0]//2
241    half_nq = nq//2
242    frame = np.fft.ifft2(fourier_frame).real
243    Iq = np.empty((nq, nq))
244    Iq[half_nq:, half_nq:] = frame[:half_nq, :half_nq]
245    Iq[:half_nq, half_nq:] = frame[-half_nq:, :half_nq]
246    Iq[half_nq:, :half_nq] = frame[:half_nq, -half_nq:]
247    Iq[:half_nq, :half_nq] = frame[-half_nq:, -half_nq:]
248    return Iq
249
250def parse_pars(model, opts):
251    # type: (ModelInfo, argparse.Namespace) -> Dict[str, float]
252
253    seed = np.random.randint(1000000) if opts.random and opts.seed < 0 else opts.seed
254    compare_opts = {
255        'info': (model.info, model.info),
256        'use_demo': False,
257        'seed': seed,
258        'mono': True,
259        'magnetic': False,
260        'values': opts.pars,
261        'show_pars': True,
262        'is2d': opts.is2d,
263    }
264    pars, pars2 = compare.parse_pars(compare_opts)
265    return pars
266
267
268def main():
269    parser = argparse.ArgumentParser(
270        description="Compute multiple scattering",
271        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
272        )
273    parser.add_argument('-k', '--power', type=int, default=0, help="show pattern for nth scattering")
274    parser.add_argument('-p', '--probability', type=float, default=0.1, help="scattering probability")
275    parser.add_argument('-n', '--nq', type=int, default=1024, help='number of mesh points')
276    parser.add_argument('-q', '--qmax', type=float, default=0.5, help='max q')
277    parser.add_argument('-w', '--window', type=float, default=2.0, help='q calc = q max * window')
278    parser.add_argument('-2', '--2d', dest='is2d', action='store_true', help='oriented sample')
279    parser.add_argument('-s', '--seed', default=-1, help='random pars with given seed')
280    parser.add_argument('-r', '--random', action='store_true', help='random pars with random seed')
281    parser.add_argument('-o', '--outfile', type=str, default="", help='random pars with random seed')
282    parser.add_argument('model', type=str, help='sas model name such as cylinder')
284    opts = parser.parse_args()
285    assert opts.nq%2 == 0, "require even # points"
286
288    pars = parse_pars(model, opts)
289    res = MultipleScattering(opts.qmax, opts.nq, opts.probability, opts.is2d,
290                             window=opts.window, power=opts.power)
291    kernel = model.make_kernel(res.q_calc)
292    #print(pars)
293    bg = pars.get('background', 0.0)
294    pars['background'] = 0.0
295    Iq_calc = call_kernel(kernel, pars)
296    Iq = res.apply(Iq_calc, background=bg, outfile=opts.outfile, plot=True)
297    plotxy(Iq)
298    import pylab
299    if opts.power > 0:
300        pylab.title('scattering power %d'%opts.power)
301    else:
302        pylab.title('multiple scattering with fraction %g'%opts.probability)
303    pylab.show()
304
305def plotxy(Iq):
306    import pylab
307    if isinstance(Iq, tuple):
308        q, Iq = Iq
309        pylab.loglog(q, Iq)
310    else:
311        data = Iq+0.
312        data[Iq <= 0] = np.min(Iq[Iq>0])/2
313        pylab.imshow(np.log10(data))
314    #import pylab; pylab.show()
315
316if __name__ == "__main__":
317    main()
Note: See TracBrowser for help on using the repository browser.