# Changeset b3703f5 in sasmodels

Ignore:
Timestamp:
Mar 26, 2018 5:50:35 PM (4 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, c462169
Parents:
802c412
Message:

lint reduction

Location:
sasmodels
Files:
3 edited

### Legend:

Unmodified
 rd86f0fc import argparse import time import os.path import numpy as np from sasmodels import core from sasmodels import compare from sasmodels import resolution2d from sasmodels.resolution import Resolution, bin_edges from sasmodels.data import empty_data1D, empty_data2D, plot_data from sasmodels.direct_model import call_kernel import sasmodels.kernelcl USE_FAST = True  # OpenCL faster, less accurate math class NumpyCalculator: class ICalculator: """ Multiple scattering calculator """ def fft(self, Iq): """ Compute the forward FFT for an image, real -> complex. """ raise NotImplementedError() def ifft(self, Iq): """ Compute the inverse FFT for an image, complex -> complex. """ raise NotImplementedError() def mulitple_scattering(self, Iq): r""" Compute multiple scattering for I(q) given scattering probability p. Given a probability p of scattering with the thickness, the expected number of scattering events, $\lambda$ is $-\log(1 - p)$, giving a Poisson weighted sum of single, double, triple, etc. scattering patterns. The number of patterns used is based on coverage (default 99%). """ raise NotImplementedError() class NumpyCalculator(ICalculator): """ Multiple scattering calculator using numpy fft. """ def __init__(self, dims=None, dtype=PRECISION): self.dtype = dtype self.complex_dtype = np.dtype('F') if dtype == np.dtype('f') else np.dtype('D') pass def fft(self, Iq): def multiple_scattering(self, Iq, p, coverage=0.99): r""" Compute multiple scattering for I(q) given scattering probability p. Given a probability p of scattering with the thickness, the expected number of scattering events, $\lambda$ is $-\log(1 - p)$, giving a Poisson weighted sum of single, double, triple, etc. scattering patterns. The number of patterns used is based on coverage (default 99%). """ #t0 = time.time() coeffs = scattering_coeffs(p, coverage) scale = np.sum(Iq) frame = _forward_shift(Iq/scale, dtype=self.dtype) F = np.fft.fft2(frame) F_convolved = F * np.polyval(poly, F) frame = np.fft.ifft2(F_convolved) fourier_frame = np.fft.fft2(frame) convolved = fourier_frame * np.polyval(poly, fourier_frame) frame = np.fft.ifft2(convolved) result = scale * _inverse_shift(frame.real, dtype=self.dtype) #print("numpy multiscat time", time.time()-t0) """ class OpenclCalculator(NumpyCalculator): class OpenclCalculator(ICalculator): """ Multiple scattering calculator using OpenCL via pyfft. """ polyval1f = None polyval1d = None context = env.get_context(dtype) if dtype == np.dtype('f'): if self.polyval1f is None: if OpenclCalculator.polyval1f is None: program = sasmodels.kernelcl.compile_model( context, POLYVAL1_KERNEL, dtype, fast=USE_FAST) self.dtype = dtype self.complex_dtype = np.dtype('F') self.polyval1 = self.polyval1f self.polyval1 = OpenclCalculator.polyval1f else: if self.polyval1d is None: if OpenclCalculator.polyval1d is None: program = sasmodels.kernelcl.compile_model( context, POLYVAL1_KERNEL, dtype, fast=False) self.dtype = dtype self.complex_type = np.dtype('D') self.polyval1 = self.polyval1d self.polyval1 = OpenclCalculator.polyval1d self.queue = env.get_queue(dtype) self.plan = pyfft.cl.Plan(dims, queue=self.queue) gpu_poly = cl_array.to_device(self.queue, poly) self.plan.execute(gpu_data.data) degree, n = poly.shape[0], frame.shape[0]*frame.shape[1] degree, data_size= poly.shape[0], frame.shape[0]*frame.shape[1] self.polyval1( self.queue, [n], None, np.int32(degree), gpu_poly.data, np.int32(n), gpu_data.data) self.queue, [data_size], None, np.int32(degree), gpu_poly.data, np.int32(data_size), gpu_data.data) self.plan.execute(gpu_data.data, inverse=True) frame = gpu_data.get() """ if transform is None: nx, ny = Iq.shape transform = Calculator(dims=(nx*2, ny*2), dtype=dtype) n_x, n_y = Iq.shape transform = Calculator(dims=(n_x*2, n_y*2), dtype=dtype) scale = np.sum(Iq) frame = _forward_shift(Iq/scale, dtype=dtype) def parse_pars(model, opts): # type: (ModelInfo, argparse.Namespace) -> Dict[str, float] """ Parse par=val arguments from the command line. """ seed = np.random.randint(1000000) if opts.random and opts.seed < 0 else opts.seed 'is2d': opts.is2d, } pars, pars2 = compare.parse_pars(compare_opts) # Note: sascomp allows comparison on a pair of models, so ignore the second. pars, _ = compare.parse_pars(compare_opts) return pars formatter_class=argparse.ArgumentDefaultsHelpFormatter, ) parser.add_argument('-p', '--probability', type=float, default=0.1, help="scattering probability") parser.add_argument('-n', '--nq', type=int, default=1024, help='number of mesh points') parser.add_argument('-q', '--qmax', type=float, default=0.5, help='max q') parser.add_argument('-w', '--window', type=float, default=2.0, help='q calc = q max * window') parser.add_argument('-2', '--2d', dest='is2d', action='store_true', help='oriented sample') parser.add_argument('-s', '--seed', default=-1, help='random pars with given seed') parser.add_argument('-r', '--random', action='store_true', help='random pars with random seed') parser.add_argument('-o', '--outfile', type=str, default="", help='random pars with random seed') parser.add_argument('model', type=str, help='sas model name such as cylinder') parser.add_argument('pars', type=str, nargs='*', help='model parameters such as radius=30') parser.add_argument('-p', '--probability', type=float, default=0.1, help="scattering probability") parser.add_argument('-n', '--nq', type=int, default=1024, help='number of mesh points') parser.add_argument('-q', '--qmax', type=float, default=0.5, help='max q') parser.add_argument('-w', '--window', type=float, default=2.0, help='q calc = q max * window') parser.add_argument('-2', '--2d', dest='is2d', action='store_true', help='oriented sample') parser.add_argument('-s', '--seed', default=-1, help='random pars with given seed') parser.add_argument('-r', '--random', action='store_true', help='random pars with random seed') parser.add_argument('-o', '--outfile', type=str, default="", help='random pars with random seed') parser.add_argument('model', type=str, help='sas model name such as cylinder') parser.add_argument('pars', type=str, nargs='*', help='model parameters such as radius=30') opts = parser.parse_args() assert opts.nq%2 == 0, "require even # points" plotxy((res._q_steps, res._q_steps), res.Iqxy+background) pylab.title("total scattering for p=%g" % probability) if res.resolution is not None: pylab.figure() plotxy((res._q_steps, res._q_steps), result) pylab.title("total scattering with resolution") else: q = res._q # Plot 1D pattern for partial scattering pylab.loglog(q, res.Iq+background, label="total for p=%g"%probability) if res.resolution is not None: pylab.loglog(q, result, label="total with dQ") #new_annulus = annular_average(res._radius, res.Iqxy, res._edges) #pylab.loglog(q, new_annulus+background, label="new total for p=%g"%probability)