Changeset b3703f5 in sasmodels for sasmodels/multiscat.py

Timestamp:

Mar 26, 2018 5:50:35 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:

c11d09f, c462169

Parents:

802c412

Message:

lint reduction

File:

• sasmodels/multiscat.py (modified) (16 diffs)

sasmodels/multiscat.py

@@ -73 +73 @@
 import argparse
 import time
-import os.path
 
 import numpy as np
@@ -81 +80 @@
 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
@@ -106 +103 @@
 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):
@@ -127 +153 @@
 
     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)
@@ -140 +158 @@
         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)
@@ -173 +191 @@
 """
 
-class OpenclCalculator(NumpyCalculator):
+class OpenclCalculator(ICalculator):
+    """
+    Multiple scattering calculator using OpenCL via pyfft.
+    """
     polyval1f = None
     polyval1d = None
@@ -180 +201 @@
         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)
@@ -187 +208 @@
             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)
@@ -196 +217 @@
             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)
@@ -229 +250 @@
         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()
@@ -251 +272 @@
     """
     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)
@@ -528 +549 @@
 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
@@ -541 +565 @@
         '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
 
@@ -550 +575 @@
         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"
@@ -607 +642 @@
             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
@@ -624 +663 @@
             # 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)