Changeset a09d55d in sasmodels
 Timestamp:
 Feb 8, 2018 5:08:55 PM (6 years ago)
 Branches:
 master, core_shell_microgels, magnetic_model, ticket1257vesicleproduct, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
 Children:
 54a5024
 Parents:
 aadec17
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

explore/multiscat.py
r049e02d ra09d55d 7 7 Usage: 8 8 9  f, fraction: the scattering probability9 p, probability: the scattering probability 10 10 q, qmax: that max q that you care about 11 11 w, window: the extension window (q is calculated for qmax*window) … … 20 20 $p$ will scatter again. 21 21 22 Let the scattering probability for a single scattering event at $q$ be $f(q)$,22 Let the scattering probability for $n$ scattering event at $q$ be $f_n(q)$, 23 23 where 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} 25 25 for $I_1(q)$, the single scattering from the system. After two scattering 26 events, the scattering will be the convolution of the first scattering27 and itself, or $(f*f)(q)$. After $n$ events it will be28 $ (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 $(=1p)$ plus the number that scattered only once 32 $(=(1p)f)$ plus the number that scattered only twice $(=(1p)(f*f))$, etc., 33 so 34 .. math:: 35 I(q) = (1p)\sum_{k=0}^{\infty} f^{*k}(q) 36 Since convolution is linear, the contributions from the individual terms 37 w ill fall off as $p^k$, so we can cut off the series38 at $n = \lceil \ln C / \ln p \rceil$ for some small fraction $C$ (default is 39 $C = 0.001$). 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 40 41 41 Using the convolution theorem, where 42 42 $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\}\} 47 44 so 48 .. math:: 49 50 f * \ldots * f = \mathcal{F}^{1}\{ F^n \} 51 45 .. math:: f * \ldots * f = \mathcal{F}^{1}\{ F^n \} 52 46 Since the Fourier transform is a linear operator, we can move the polynomial 53 47 expression for the convolution into the transform, giving 54 48 .. math:: 55 56 I(q) = \mathcal{F}^{1}\left\{ (1p) \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\} 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. 65 63 """ 66 64 67 from __future__ import print_function 65 from __future__ import print_function, division 68 66 69 67 import argparse … … 71 69 72 70 import numpy as np 71 from scipy.special import gamma 73 72 74 73 from sasmodels import core … … 81 80 82 81 class 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): 84 83 self.qmax = qmax 85 84 self.nq = nq 86 85 self.power = power 87 self. fraction = fraction86 self.probability = probability 88 87 self.is2d = is2d 89 88 self.window = window … … 107 106 self.q_range = q_range 108 107 109 def apply(self, theory ):110 t0 = time.time()108 def apply(self, theory, background=0.0, outfile="", plot=False): 109 #t0 = time.time() 111 110 if self.is2d: 112 111 Iq_calc = theory … … 117 116 #plotxy(Iq_calc); import pylab; pylab.figure() 118 117 if self.power > 0: 119 Iqxy = autoconv_n(Iq_calc, self.power) 118 Iqxy = scattering_power(Iq_calc, self.power) 119 powers = [] 120 120 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 123 128 #plotxy(Iqxy); import pylab; pylab.figure() 124 129 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: 126 137 import pylab 127 138 plotxy(Iq_calc) … … 129 140 pylab.figure() 130 141 131 return Iqxy 142 return Iqxy + background 132 143 else: 133 144 # circular average, no antialiasing 134 145 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 antialiasing for individual powers 155 data = np.vstack([q_corners] + Iq_powers).T 156 np.savetxt(outfile + "_powers.txt", data) 157 if plot: 137 158 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)) 144 166 pylab.legend() 145 167 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 171 def 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 177 def 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(1p) 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**(k1)/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 204 def 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! * (1exp(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 219 def _forward_fft(Iq): 220 # Prepare padded array and forward transform 221 nq = Iq.shape[0] 155 222 half_nq = nq//2 156 223 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 = (1a)f + (1a)af^2 + (1a)a^2f^3 + ... 168 # = (1a)f[1 + af + (af)^2 + (af)^3 + ...] 169 multiple_scattering = ( 170 (1frac)*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 231 def _inverse_fft(fourier_frame): 232 # Invert the transform and recover the transformed data 233 nq = fourier_frame.shape[0]//2 198 234 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_scatter1): 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 259 242 260 243 def parse_pars(model, opts): … … 281 264 formatter_class=argparse.ArgumentDefaultsHelpFormatter, 282 265 ) 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") 285 268 parser.add_argument('n', 'nq', type=int, default=1024, help='number of mesh points') 286 269 parser.add_argument('q', 'qmax', type=float, default=0.5, help='max q') … … 289 272 parser.add_argument('s', 'seed', default=1, help='random pars with given seed') 290 273 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') 291 275 parser.add_argument('model', type=str, help='sas model name such as cylinder') 292 276 parser.add_argument('pars', type=str, nargs='*', help='model parameters such as radius=30') … … 296 280 model = core.load_model(opts.model) 297 281 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, 299 283 window=opts.window, power=opts.power) 300 284 kernel = model.make_kernel(res.q_calc) … … 303 287 pars['background'] = 0.0 304 288 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) 310 290 plotxy(Iq) 311 291 import pylab … … 313 293 pylab.title('scattering power %d'%opts.power) 314 294 else: 315 pylab.title('multiple scattering with fraction %g'%opts. fraction)295 pylab.title('multiple scattering with fraction %g'%opts.probability) 316 296 pylab.show() 317 297
Note: See TracChangeset
for help on using the changeset viewer.