Changeset bc248f8 in sasmodels
- Timestamp:
- Mar 4, 2018 6:43:57 PM (7 years ago)
- Branches:
- master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- c11d09f, 49d1f8b8
- Parents:
- 1ee5ac6
- Files:
-
- 1 added
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
explore/multiscat.py
re309e23 rbc248f8 61 61 The usual pinhole or slit resolution calculation can performed from these 62 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. 63 69 """ 64 70 … … 67 73 import argparse 68 74 import time 75 import os.path 69 76 70 77 import numpy as np 78 from numpy import pi 71 79 from scipy.special import gamma 72 80 … … 77 85 from sasmodels.data import empty_data1D, empty_data2D, plot_data 78 86 from sasmodels.direct_model import call_kernel 79 80 81 class 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()] 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("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 96 190 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): 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 109 204 #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 175 def scattering_power(Iq, n): 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("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!$. 176 251 """ 177 Calculate the nth scattering power as a distribution. To get the178 weighted contribution, scale by $\lambda^k e^{-\lambda}/k!$.179 """252 if transform is None: 253 nx, ny = Iq.shape 254 transform = Calculator(dims=(nx*2, ny*2), dtype=dtype) 180 255 scale = np.sum(Iq) 181 F = _forward_fft(Iq/scale) 182 result = _inverse_fft(F**n) 183 return result 184 185 def 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 """ 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): 194 263 L = -np.log(1-p) 195 264 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 265 coeffs = [L**k/gamma(k+2) for k in range(num_scatter)] 266 return coeffs 210 267 211 268 def truncated_poisson_invcdf(p, L): 212 """269 r""" 213 270 Return smallest k such that cdf(k; L) > p from the truncated Poisson 214 271 probability excluding k=0 … … 224 281 return k 225 282 226 def _forward_ fft(Iq):283 def _forward_shift(Iq, dtype=PRECISION): 227 284 # Prepare padded array and forward transform 228 285 nq = Iq.shape[0] 229 286 half_nq = nq//2 230 frame = np.zeros((2*nq, 2*nq) )287 frame = np.zeros((2*nq, 2*nq), dtype=dtype) 231 288 frame[:half_nq, :half_nq] = Iq[half_nq:, half_nq:] 232 289 frame[-half_nq:, :half_nq] = Iq[:half_nq, half_nq:] 233 290 frame[:half_nq, -half_nq:] = Iq[half_nq:, :half_nq] 234 291 frame[-half_nq:, -half_nq:] = Iq[:half_nq, :half_nq] 235 fourier_frame = np.fft.fft2(frame) 236 return fourier_frame 237 238 def _inverse_fft(fourier_frame): 292 return frame 293 294 def _inverse_shift(frame, dtype=PRECISION): 239 295 # Invert the transform and recover the transformed data 240 nq = f ourier_frame.shape[0]//2296 nq = frame.shape[0]//2 241 297 half_nq = nq//2 242 frame = np.fft.ifft2(fourier_frame).real 243 Iq = np.empty((nq, nq)) 298 Iq = np.empty((nq, nq), dtype=dtype) 244 299 Iq[half_nq:, half_nq:] = frame[:half_nq, :half_nq] 245 300 Iq[:half_nq, half_nq:] = frame[-half_nq:, :half_nq] … … 247 302 Iq[:half_nq, :half_nq] = frame[-half_nq:, -half_nq:] 248 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 249 512 250 513 def parse_pars(model, opts): … … 259 522 'magnetic': False, 260 523 'values': opts.pars, 261 'show_pars': True, 524 #'show_pars': True, 525 'show_pars': False, 262 526 'is2d': opts.is2d, 263 527 } … … 271 535 formatter_class=argparse.ArgumentDefaultsHelpFormatter, 272 536 ) 273 parser.add_argument('-k', '--power', type=int, default=0, help="show pattern for nth scattering")274 537 parser.add_argument('-p', '--probability', type=float, default=0.1, help="scattering probability") 275 538 parser.add_argument('-n', '--nq', type=int, default=1024, help='number of mesh points') … … 287 550 model = core.load_model(opts.model) 288 551 pars = parse_pars(model, opts) 289 res = MultipleScattering( opts.qmax, opts.nq, opts.probability, opts.is2d,290 window=opts.window, power=opts.power)552 res = MultipleScattering(qmax=opts.qmax, nq=opts.nq, window=opts.window, 553 probability=opts.probability, is2d=opts.is2d) 291 554 kernel = model.make_kernel(res.q_calc) 292 555 #print(pars) 293 556 bg = pars.get('background', 0.0) 294 557 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) 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.): 298 563 import pylab 299 if opts.power > 0: 300 pylab.title('scattering power %d'%opts.power) 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 301 570 else: 302 pylab.title('multiple scattering with fraction %g'%opts.probability) 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) 303 617 pylab.show() 304 618 305 def plotxy( Iq):619 def plotxy(q, Iq): 306 620 import pylab 307 if isinstance(Iq, tuple):308 q, Iq = Iq309 pylab.loglog(q , Iq)621 # q is a tuple of (q,) or (qx, qy) 622 if len(q) == 1: 623 pylab.loglog(q[0], Iq) 310 624 else: 311 data = Iq +0.312 data[Iq <= 0] = np.min(Iq[Iq >0])/2625 data = Iq.copy() 626 data[Iq <= 0] = np.min(Iq[Iq > 0])/2 313 627 pylab.imshow(np.log10(data)) 314 #import pylab; pylab.show()315 628 316 629 if __name__ == "__main__":
Note: See TracChangeset
for help on using the changeset viewer.