source: sasmodels/sasmodels/multiscat.py @ 49d1f8b8

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

move multiscat into sasmodels namespace to make it easier to run the example

  • Property mode set to 100644
File size: 24.9 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
64By default single precision OpenCL is used for the calculation.  Set the
65environment variable *SAS_OPENCL=none* to use double precision numpy FFT
66instead.  The OpenCL versions is about 10x faster on an elderly Mac with
67Intel HD 4000 graphics.  The single precision numerical artifacts don't
68seem to seriously impact overall accuracy, though they look pretty bad.
69"""
70
71from __future__ import print_function, division
72
73import argparse
74import time
75import os.path
76
77import numpy as np
78from numpy import pi
79from scipy.special import gamma
80
81from sasmodels import core
82from sasmodels import compare
83from sasmodels import resolution2d
84from sasmodels.resolution import Resolution, bin_edges
85from sasmodels.data import empty_data1D, empty_data2D, plot_data
86from sasmodels.direct_model import call_kernel
87import 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
99try:
100    import pyfft.cl
101    import pyopencl.array as cl_array
102    HAVE_OPENCL = sasmodels.kernelcl.use_opencl()
103except ImportError:
104    HAVE_OPENCL = False
105PRECISION = np.dtype('f' if HAVE_OPENCL else 'd')  # 'f' or 'd'
106USE_FAST = True  # OpenCL faster, less accurate math
107
108class 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.
156POLYVAL1_KERNEL = """
157kernel 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
175class 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
242Calculator = OpenclCalculator if HAVE_OPENCL else NumpyCalculator
243
244def 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
262def 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
268def 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
283def _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
294def _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
306class 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
478def 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
497def 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
513def 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
532def 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
562def 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
619def 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
629if __name__ == "__main__":
630    main()
Note: See TracBrowser for help on using the repository browser.