Changeset a09d55d in sasmodels for explore


Ignore:
Timestamp:
Feb 8, 2018 5:08:55 PM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
54a5024
Parents:
aadec17
Message:

multiscat: use correct weightings for the multiscattering contributions

File:
1 edited

Legend:

Unmodified
Added
Removed

  • explore/multiscat.py

    r049e02d ra09d55d  
    77Usage: 
    88 
    9     -f, --fraction: the scattering probability 
     9    -p, --probability: the scattering probability 
    1010    -q, --qmax: that max q that you care about 
    1111    -w, --window: the extension window (q is calculated for qmax*window) 
    … …  
    2020$p$ will scatter again. 
    2121 
    22 Let the scattering probability for a single scattering event at $q$ be $f(q)$, 
     22Let the scattering probability for $n$ scattering event at $q$ be $f_n(q)$, 
    2323where 
    24 .. math:: f(q) = p \frac{I_1(q)}{\int I_1(q) dq} 
     24.. math:: f_1(q) = \frac{I_1(q)}{\int I_1(q) dq} 
    2525for $I_1(q)$, the single scattering from the system. After two scattering 
    26 events, the scattering will be the convolution of the first scattering 
    27 and itself, or $(f*f)(q)$.  After $n$ events it will be 
    28 $(f* \cdots * f)(q)$. 
    29  
    30 The total scattering after multiple events will then be the number that didn't 
    31 scatter the first time $(=1-p)$ plus the number that scattered only once 
    32 $(=(1-p)f)$ plus the number that scattered only twice $(=(1-p)(f*f))$, etc., 
    33 so 
    34 .. math:: 
    35     I(q) = (1-p)\sum_{k=0}^{\infty} f^{*k}(q) 
    36 Since convolution is linear, the contributions from the individual terms 
    37 will fall off as $p^k$, so we can cut off the series 
    38 at $n = \lceil \ln C / \ln p \rceil$ for some small fraction $C$ (default is 
    39 $C = 0.001$). 
     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 
    4040 
    4141Using the convolution theorem, where 
    4242$F = \mathcal{F}(f)$ is the Fourier transform, 
    43 .. math:: 
    44  
    45     f * g = \mathcal{F}^{-1}\{\mathcal{F}\{f\} \cdot \mathcal{F}\{g\}\} 
    46  
     43.. math:: f * g = \mathcal{F}^{-1}\{\mathcal{F}\{f\} \cdot \mathcal{F}\{g\}\} 
    4744so 
    48 .. math:: 
    49  
    50     f * \ldots * f = \mathcal{F}^{-1}\{ F^n \} 
    51  
     45.. math:: f * \ldots * f = \mathcal{F}^{-1}\{ F^n \} 
    5246Since the Fourier transform is a linear operator, we can move the polynomial 
    5347expression for the convolution into the transform, giving 
    5448.. math:: 
    55  
    56     I(q) = \mathcal{F}^{-1}\left\{ (1-p) \sum_{k=0}^{n} F^k \right\} 
    57  
    58 We drop the transmission term, $k=0$, and rescale the result by the 
    59 total scattering $\int I_1(q) dq$. 
    60  
    61 For speed we use the fast fourier transform for a power of two.  The resulting 
    62 $I(q)$ will be linearly spaced and likely heavily oversampled.  The usual 
    63 pinhole or slit resolution calculation can performed from these calculated 
    64 values. 
     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. 
    6563""" 
    6664 
    67 from __future__ import print_function 
     65from __future__ import print_function, division 
    6866 
    6967import argparse 
    … …  
    7169 
    7270import numpy as np 
     71from scipy.special import gamma 
    7372 
    7473from sasmodels import core 
    … …  
    8180 
    8281class MultipleScattering(Resolution): 
    83     def __init__(self, qmax, nq, fraction, is2d, window=2, power=0): 
     82    def __init__(self, qmax, nq, probability, is2d, window=2, power=0): 
    8483        self.qmax = qmax 
    8584        self.nq = nq 
    8685        self.power = power 
    87         self.fraction = fraction 
     86        self.probability = probability 
    8887        self.is2d = is2d 
    8988        self.window = window 
    … …  
    107106        self.q_range = q_range 
    108107 
    109     def apply(self, theory): 
    110         t0 = time.time() 
     108    def apply(self, theory, background=0.0, outfile="", plot=False): 
     109        #t0 = time.time() 
    111110        if self.is2d: 
    112111            Iq_calc = theory 
    … …  
    117116        #plotxy(Iq_calc); import pylab; pylab.figure() 
    118117        if self.power > 0: 
    119             Iqxy = autoconv_n(Iq_calc, self.power) 
     118            Iqxy = scattering_power(Iq_calc, self.power) 
     119            powers = [] 
    120120        else: 
    121             Iqxy = multiple_scattering(Iq_calc, self.fraction, cutoff=0.001) 
    122         print("multiple scattering calc time", time.time()-t0) 
     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 
    123128        #plotxy(Iqxy); import pylab; pylab.figure() 
    124129        if self.is2d: 
    125             if 1: 
     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: 
    126137                import pylab 
    127138                plotxy(Iq_calc) 
    … …  
    129140                pylab.figure() 
    130141 
    131             return Iqxy 
     142            return Iqxy + background 
    132143        else: 
    133144            # circular average, no anti-aliasing 
    134145            Iq = np.histogram(self._qxy, bins=self._edges, weights=Iqxy)[0]/self._norm 
    135  
    136             if 1: 
     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: 
    137158                import pylab 
    138                 pylab.loglog(q_corners, theory, label="single scattering") 
    139                 if self.power > 0: 
    140                     label = "scattering power %d"%self.power 
    141                 else: 
    142                     label = "scattering fraction %d"%self.fraction 
    143                 pylab.loglog(q_corners, Iq, label=label) 
     159                plotxy(Iqxy) 
     160                pylab.title("multiple scattering") 
     161                pylab.figure() 
     162            if plot and powers: 
     163                import pylab 
     164                for k, Ipower in enumerate(Iq_powers): 
     165                    pylab.loglog(q_corners, Ipower, label="scattering=%d"%(k+1)) 
    144166                pylab.legend() 
    145167                pylab.figure() 
    146                 return Iqxy 
    147             return q_corners, Iq 
    148  
    149 def multiple_scattering(Iq_calc, frac, cutoff=0.001): 
    150     #plotxy(Iq_calc) 
    151     num_scatter = int(np.ceil(np.log(cutoff)/np.log(frac))) 
    152  
    153     # Prepare padded array for transform 
    154     nq = Iq_calc.shape[0] 
     168 
     169            return q_corners, Iq + background 
     170 
     171def scattering_power(Iq, n): 
     172    scale = np.sum(Iq) 
     173    F = _forward_fft(Iq/scale) 
     174    result = scale * _inverse_fft(F**n) 
     175    return result 
     176 
     177def multiple_scattering(Iq, p, coverage=0.99, return_powers=False): 
     178    """ 
     179    Compute multiple scattering for I(q) given scattering probability p. 
     180 
     181    Given a probability p of scattering with the thickness, the expected 
     182    number of scattering events, $\lambda$ is $-\log(1 - p)$, giving a 
     183    Poisson weighted sum of single, double, triple, ... scattering patterns. 
     184    The number of patterns used is based on coverage (default 99%).  If 
     185    return_poi 
     186    """ 
     187    L = -np.log(1-p) 
     188    num_scatter = truncated_poisson_invcdf(coverage, L) 
     189 
     190    # Compute multiple scattering via convolution. 
     191    scale = np.sum(Iq) 
     192    F = _forward_fft(Iq/scale) 
     193    coeffs = [L**(k-1)/gamma(k+1) for k in range(1, num_scatter+1)] 
     194    multiple_scattering = F * np.polyval(coeffs[::-1], F) 
     195    result = scale * _inverse_fft(multiple_scattering) 
     196 
     197    if return_powers: 
     198        powers = [scale * _inverse_fft(F**k) for k in range(1, num_scatter+1)] 
     199    else: 
     200        powers = [] 
     201 
     202    return result, powers 
     203 
     204def truncated_poisson_invcdf(p, L): 
     205    """ 
     206    Return smallest k such that cdf(k; L) > p from the truncated Poisson 
     207    probability excluding k=0 
     208    """ 
     209    # pmf(k; L) = L**k * exp(-L) / (k! * (1-exp(-L)) 
     210    cdf = 0 
     211    pmf = -np.exp(-L) / np.expm1(-L) 
     212    k = 0 
     213    while cdf < p: 
     214        k += 1 
     215        pmf *= L/k 
     216        cdf += pmf 
     217    return k 
     218 
     219def _forward_fft(Iq): 
     220    # Prepare padded array and forward transform 
     221    nq = Iq.shape[0] 
    155222    half_nq = nq//2 
    156223    frame = np.zeros((2*nq, 2*nq)) 
    157     frame[:half_nq, :half_nq] = Iq_calc[half_nq:, half_nq:] 
    158     frame[-half_nq:, :half_nq] = Iq_calc[:half_nq, half_nq:] 
    159     frame[:half_nq, -half_nq:] = Iq_calc[half_nq:, :half_nq] 
    160     frame[-half_nq:, -half_nq:] = Iq_calc[:half_nq, :half_nq] 
    161     #plotxy(frame) 
    162  
    163     # Compute multiple scattering via convolution. 
    164     scale = np.sum(Iq_calc) 
    165     fourier_frame = np.fft.fft2(frame/scale) 
    166     #plotxy(abs(frame)) 
    167     # total = (1-a)f + (1-a)af^2 + (1-a)a^2f^3 + ... 
    168     #       = (1-a)f[1 + af + (af)^2 + (af)^3 + ...] 
    169     multiple_scattering = ( 
    170         (1-frac)*fourier_frame 
    171         *np.polyval(np.ones(num_scatter), frac*fourier_frame)) 
    172     conv_frame = scale*np.fft.ifft2(multiple_scattering).real 
    173  
    174     # Recover the transformed data 
    175     #plotxy(conv_frame) 
    176     Iq_conv = np.empty((nq, nq)) 
    177     Iq_conv[half_nq:, half_nq:] = conv_frame[:half_nq, :half_nq] 
    178     Iq_conv[:half_nq, half_nq:] = conv_frame[-half_nq:, :half_nq] 
    179     Iq_conv[half_nq:, :half_nq] = conv_frame[:half_nq, -half_nq:] 
    180     Iq_conv[:half_nq, :half_nq] = conv_frame[-half_nq:, -half_nq:] 
    181     #plotxy(Iq_conv) 
    182     return Iq_conv 
    183  
    184 def multiple_scattering_cl(Iq_calc, frac, cutoff=0.001): 
    185     raise NotImplementedError("no support for opencl calculations at this time") 
    186  
    187     import pyopencl as cl 
    188     import pyopencl.array as cla 
    189     from gpyfft.fft import FFT 
    190     context = cl.create_some_context() 
    191     queue = cl.CommandQueue(context) 
    192  
    193     #plotxy(Iq_calc) 
    194     num_scatter = int(np.ceil(np.log(cutoff)/np.log(frac))) 
    195  
    196     # Prepare padded array for transform 
    197     nq = Iq_calc.shape[0] 
     224    frame[:half_nq, :half_nq] = Iq[half_nq:, half_nq:] 
     225    frame[-half_nq:, :half_nq] = Iq[:half_nq, half_nq:] 
     226    frame[:half_nq, -half_nq:] = Iq[half_nq:, :half_nq] 
     227    frame[-half_nq:, -half_nq:] = Iq[:half_nq, :half_nq] 
     228    fourier_frame = np.fft.fft2(frame) 
     229    return fourier_frame 
     230 
     231def _inverse_fft(fourier_frame): 
     232    # Invert the transform and recover the transformed data 
     233    nq = fourier_frame.shape[0]//2 
    198234    half_nq = nq//2 
    199     frame = np.zeros((2*nq, 2*nq), dtype='float32') 
    200     frame[:half_nq, :half_nq] = Iq_calc[half_nq:, half_nq:] 
    201     frame[-half_nq:, :half_nq] = Iq_calc[:half_nq, half_nq:] 
    202     frame[:half_nq, -half_nq:] = Iq_calc[half_nq:, :half_nq] 
    203     frame[-half_nq:, -half_nq:] = Iq_calc[:half_nq, :half_nq] 
    204     #plotxy(frame) 
    205  
    206     # Compute multiple scattering via convolution (OpenCL operations) 
    207     frame_gpu = cla.to_device(queue, frame) 
    208     fourier_frame_gpu = cla.zeros(frame.shape, dtype='complex64') 
    209     scale = frame_gpu.sum() 
    210     frame_gpu /= scale 
    211     transform = FFT(context, queue, frame_gpu, fourier_frame_gpu, axes=(0,1)) 
    212     event, = transform.enqueue() 
    213     event.wait() 
    214     fourier_frame_gpu *= frac 
    215     multiple_scattering_gpu = fourier_frame_gpu.copy() 
    216     for _ in range(num_scatter-1): 
    217         multiple_scattering_gpu += 1 
    218         multiple_scattering_gpu *= fourier_frame_gpu 
    219     multiple_scattering_gpu *= (1 - frac)/frac 
    220     transform = FFT(context, queue, multiple_scattering_gpu, frame_gpu, axes=(0,1)) 
    221     event, = transform.enqueue(forward=False) 
    222     event.wait() 
    223     conv_frame = cla.from_device(queue, frame_gpu) 
    224  
    225     # Recover the transformed data 
    226     #plotxy(conv_frame) 
    227     Iq_conv = np.empty((nq, nq)) 
    228     Iq_conv[half_nq:, half_nq:] = conv_frame[:half_nq, :half_nq] 
    229     Iq_conv[:half_nq, half_nq:] = conv_frame[-half_nq:, :half_nq] 
    230     Iq_conv[half_nq:, :half_nq] = conv_frame[:half_nq, -half_nq:] 
    231     Iq_conv[:half_nq, :half_nq] = conv_frame[-half_nq:, -half_nq:] 
    232     #plotxy(Iq_conv) 
    233     return Iq_conv 
    234  
    235 def autoconv_n(Iq_calc, power): 
    236     # Compute multiple scattering via convolution. 
    237     #plotxy(Iq_calc) 
    238     scale = np.sum(Iq_calc) 
    239     nq = Iq_calc.shape[0] 
    240     frame = np.zeros((2*nq, 2*nq)) 
    241     half_nq = nq//2 
    242     frame[:half_nq, :half_nq] = Iq_calc[half_nq:, half_nq:] 
    243     frame[-half_nq:, :half_nq] = Iq_calc[:half_nq, half_nq:] 
    244     frame[:half_nq, -half_nq:] = Iq_calc[half_nq:, :half_nq] 
    245     frame[-half_nq:, -half_nq:] = Iq_calc[:half_nq, :half_nq] 
    246     #plotxy(frame) 
    247     fourier_frame = np.fft.fft2(frame/scale) 
    248     #plotxy(abs(frame)) 
    249     fourier_frame = fourier_frame**power 
    250     conv_frame = scale*np.fft.ifft2(fourier_frame).real 
    251     #plotxy(conv_frame) 
    252     Iq_conv = np.empty((nq, nq)) 
    253     Iq_conv[half_nq:, half_nq:] = conv_frame[:half_nq, :half_nq] 
    254     Iq_conv[:half_nq, half_nq:] = conv_frame[-half_nq:, :half_nq] 
    255     Iq_conv[half_nq:, :half_nq] = conv_frame[:half_nq, -half_nq:] 
    256     Iq_conv[:half_nq, :half_nq] = conv_frame[-half_nq:, -half_nq:] 
    257     #plotxy(Iq_conv) 
    258     return Iq_conv 
     235    frame = np.fft.ifft2(fourier_frame).real 
     236    Iq = np.empty((nq, nq)) 
     237    Iq[half_nq:, half_nq:] = frame[:half_nq, :half_nq] 
     238    Iq[:half_nq, half_nq:] = frame[-half_nq:, :half_nq] 
     239    Iq[half_nq:, :half_nq] = frame[:half_nq, -half_nq:] 
     240    Iq[:half_nq, :half_nq] = frame[-half_nq:, -half_nq:] 
     241    return Iq 
    259242 
    260243def parse_pars(model, opts): 
    … …  
    281264        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
    282265        ) 
    283     parser.add_argument('-p', '--power', type=int, default=0, help="show pattern for nth scattering") 
    284     parser.add_argument('-f', '--fraction', type=float, default=0.1, help="scattering fraction") 
     266    parser.add_argument('-k', '--power', type=int, default=0, help="show pattern for nth scattering") 
     267    parser.add_argument('-p', '--probability', type=float, default=0.1, help="scattering probability") 
    285268    parser.add_argument('-n', '--nq', type=int, default=1024, help='number of mesh points') 
    286269    parser.add_argument('-q', '--qmax', type=float, default=0.5, help='max q') 
    … …  
    289272    parser.add_argument('-s', '--seed', default=-1, help='random pars with given seed') 
    290273    parser.add_argument('-r', '--random', action='store_true', help='random pars with random seed') 
     274    parser.add_argument('-o', '--outfile', type=str, default="", help='random pars with random seed') 
    291275    parser.add_argument('model', type=str, help='sas model name such as cylinder') 
    292276    parser.add_argument('pars', type=str, nargs='*', help='model parameters such as radius=30') 
    … …  
    296280    model = core.load_model(opts.model) 
    297281    pars = parse_pars(model, opts) 
    298     res = MultipleScattering(opts.qmax, opts.nq, opts.fraction, opts.is2d, 
     282    res = MultipleScattering(opts.qmax, opts.nq, opts.probability, opts.is2d, 
    299283                             window=opts.window, power=opts.power) 
    300284    kernel = model.make_kernel(res.q_calc) 
    … …  
    303287    pars['background'] = 0.0 
    304288    Iq_calc = call_kernel(kernel, pars) 
    305     t0 = time.time() 
    306     for i in range(10): 
    307         Iq_calc = call_kernel(kernel, pars) 
    308     print("single scattering calc time", (time.time()-t0)/10) 
    309     Iq = res.apply(Iq_calc) + bg 
     289    Iq = res.apply(Iq_calc, background=bg, outfile=opts.outfile, plot=True) 
    310290    plotxy(Iq) 
    311291    import pylab 
    … …  
    313293        pylab.title('scattering power %d'%opts.power) 
    314294    else: 
    315         pylab.title('multiple scattering with fraction %g'%opts.fraction) 
     295        pylab.title('multiple scattering with fraction %g'%opts.probability) 
    316296    pylab.show() 
    317297 
Note: See TracChangeset for help on using the changeset viewer.