Changeset 0ceca73 in sasmodels

Timestamp:

Apr 15, 2016 11:24:16 AM (9 years ago)

Author:

wojciech

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

299dcce

Parents:

391ea92 (diff), 3599d36 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'polydisp' of ​https://github.com/SasView/sasmodels into polydisp

Location:

sasmodels

Files:

• compare.py (modified) (41 diffs)
• core.py (modified) (6 diffs)
• kernelcl.py (modified) (22 diffs)
• models/spherical_sld.c (added)

sasmodels/compare.py

@@ -41 +41 @@
 from .direct_model import DirectModel
 from .convert import revert_name, revert_pars, constrain_new_to_old
+
+try:
+    from typing import Optional, Dict, Any, Callable, Tuple
+except:
+    pass
+else:
+    from .modelinfo import ModelInfo, Parameter, ParameterSet
+    from .data import Data
+    Calculator = Callable[[float, ...], np.ndarray]
 
 USAGE = """
@@ -97 +106 @@
 kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True
 
-MODELS = core.list_models()
-
 # CRUFT python 2.6
 if not hasattr(datetime.timedelta, 'total_seconds'):
@@ -160 +167 @@
         ...            print(randint(0,1000000,3))
         ...            raise Exception()
-        ...    except:
+        ...    except Exception:
         ...        print("Exception raised")
         ...    print(randint(0,1000000))
@@ -169 +176 @@
     """
     def __init__(self, seed=None):
+        # type: (Optional[int]) -> None
         self._state = np.random.get_state()
         np.random.seed(seed)
 
     def __enter__(self):
-        return None
-
-    def __exit__(self, *args):
+        # type: () -> None
+        pass
+
+    def __exit__(self, type, value, traceback):
+        # type: (Any, BaseException, Any) -> None
+        # TODO: better typing for __exit__ method
         np.random.set_state(self._state)
 
 def tic():
+    # type: () -> Callable[[], float]
     """
     Timer function.
@@ -190 +202 @@
 
 def set_beam_stop(data, radius, outer=None):
+    # type: (Data, float, float) -> None
     """
     Add a beam stop of the given *radius*.  If *outer*, make an annulus.
 
-    Note: this function does not use the sasview package
+    Note: this function does not require sasview
     """
     if hasattr(data, 'qx_data'):
@@ -207 +220 @@
 
 def parameter_range(p, v):
+    # type: (str, float) -> Tuple[float, float]
     """
     Choose a parameter range based on parameter name and initial value.
@@ -212 +226 @@
     # process the polydispersity options
     if p.endswith('_pd_n'):
-        return [0, 100]
+        return 0., 100.
     elif p.endswith('_pd_nsigma'):
-        return [0, 5]
+        return 0., 5.
     elif p.endswith('_pd_type'):
-        return v
+        raise ValueError("Cannot return a range for a string value")
     elif any(s in p for s in ('theta', 'phi', 'psi')):
         # orientation in [-180,180], orientation pd in [0,45]
         if p.endswith('_pd'):
-            return [0, 45]
+            return 0., 45.
         else:
-            return [-180, 180]
+            return -180., 180.
     elif p.endswith('_pd'):
-        return [0, 1]
+        return 0., 1.
     elif 'sld' in p:
-        return [-0.5, 10]
+        return -0.5, 10.
     elif p == 'background':
-        return [0, 10]
+        return 0., 10.
     elif p == 'scale':
-        return [0, 1e3]
-    elif v < 0:
-        return [2*v, -2*v]
+        return 0., 1.e3
+    elif v < 0.:
+        return 2.*v, -2.*v
     else:
-        return [0, (2*v if v > 0 else 1)]
+        return 0., (2.*v if v > 0. else 1.)
 
 
 def _randomize_one(model_info, p, v):
+    # type: (ModelInfo, str, float) -> float
+    # type: (ModelInfo, str, str) -> str
     """
     Randomize a single parameter.
@@ -263 +279 @@
 
 def randomize_pars(model_info, pars, seed=None):
+    # type: (ModelInfo, ParameterSet, int) -> ParameterSet
     """
     Generate random values for all of the parameters.
@@ -273 +290 @@
     with push_seed(seed):
         # Note: the sort guarantees order `of calls to random number generator
-        pars = dict((p, _randomize_one(model_info, p, v))
-                    for p, v in sorted(pars.items()))
-    return pars
+        random_pars = dict((p, _randomize_one(model_info, p, v))
+                           for p, v in sorted(pars.items()))
+    return random_pars
 
 def constrain_pars(model_info, pars):
+    # type: (ModelInfo, ParameterSet) -> None
     """
     Restrict parameters to valid values.
@@ -284 +302 @@
     which need to support within model constraints (cap radius more than
     cylinder radius in this case).
+
+    Warning: this updates the *pars* dictionary in place.
     """
     name = model_info.id
@@ -315 +335 @@
 
 def parlist(model_info, pars, is2d):
+    # type: (ModelInfo, ParameterSet, bool) -> str
     """
     Format the parameter list for printing.
@@ -326 +347 @@
             n=int(pars.get(p.id+"_pd_n", 0)),
             nsigma=pars.get(p.id+"_pd_nsgima", 3.),
-            type=pars.get(p.id+"_pd_type", 'gaussian'))
+            pdtype=pars.get(p.id+"_pd_type", 'gaussian'),
+        )
         lines.append(_format_par(p.name, **fields))
     return "\n".join(lines)
@@ -332 +354 @@
     #return "\n".join("%s: %s"%(p, v) for p, v in sorted(pars.items()))
 
-def _format_par(name, value=0., pd=0., n=0, nsigma=3., type='gaussian'):
+def _format_par(name, value=0., pd=0., n=0, nsigma=3., pdtype='gaussian'):
+    # type: (str, float, float, int, float, str) -> str
     line = "%s: %g"%(name, value)
     if pd != 0.  and n != 0:
         line += " +/- %g  (%d points in [-%g,%g] sigma %s)"\
-                % (pd, n, nsigma, nsigma, type)
+                % (pd, n, nsigma, nsigma, pdtype)
     return line
 
 def suppress_pd(pars):
+    # type: (ParameterSet) -> ParameterSet
     """
     Suppress theta_pd for now until the normalization is resolved.
@@ -352 +376 @@
 
 def eval_sasview(model_info, data):
+    # type: (Modelinfo, Data) -> Calculator
     """
     Return a model calculator using the pre-4.0 SasView models.
@@ -359 +384 @@
     import sas
     from sas.models.qsmearing import smear_selection
+    import sas.models
 
     def get_model(name):
+        # type: (str) -> "sas.models.BaseComponent"
         #print("new",sorted(_pars.items()))
-        sas = __import__('sas.models.' + name)
+        __import__('sas.models.' + name)
         ModelClass = getattr(getattr(sas.models, name, None), name, None)
         if ModelClass is None:
@@ -400 +427 @@
 
     def calculator(**pars):
+        # type: (float, ...) -> np.ndarray
         """
         Sasview calculator for model.
@@ -406 +434 @@
         pars = revert_pars(model_info, pars)
         for k, v in pars.items():
-            parts = k.split('.')  # polydispersity components
-            if len(parts) == 2:
-                model.dispersion[parts[0]][parts[1]] = v
+            name_attr = k.split('.')  # polydispersity components
+            if len(name_attr) == 2:
+                model.dispersion[name_attr[0]][name_attr[1]] = v
             else:
                 model.setParam(k, v)
@@ -428 +456 @@
 }
 def eval_opencl(model_info, data, dtype='single', cutoff=0.):
+    # type: (ModelInfo, Data, str, float) -> Calculator
     """
     Return a model calculator using the OpenCL calculation engine.
@@ -442 +471 @@
 
 def eval_ctypes(model_info, data, dtype='double', cutoff=0.):
+    # type: (ModelInfo, Data, str, float) -> Calculator
     """
     Return a model calculator using the DLL calculation engine.
     """
-    if dtype == 'quad':
-        dtype = 'longdouble'
     model = core.build_model(model_info, dtype=dtype, platform="dll")
     calculator = DirectModel(data, model, cutoff=cutoff)
@@ -453 +481 @@
 
 def time_calculation(calculator, pars, Nevals=1):
+    # type: (Calculator, ParameterSet, int) -> Tuple[np.ndarray, float]
     """
     Compute the average calculation time over N evaluations.
@@ -461 +490 @@
     # initialize the code so time is more accurate
     if Nevals > 1:
-        value = calculator(**suppress_pd(pars))
+        calculator(**suppress_pd(pars))
     toc = tic()
-    for _ in range(max(Nevals, 1)):  # make sure there is at least one eval
+    # make sure there is at least one eval
+    value = calculator(**pars)
+    for _ in range(Nevals-1):
         value = calculator(**pars)
     average_time = toc()*1000./Nevals
@@ -469 +500 @@
 
 def make_data(opts):
+    # type: (Dict[str, Any]) -> Tuple[Data, np.ndarray]
     """
     Generate an empty dataset, used with the model to set Q points
@@ -478 +510 @@
     qmax, nq, res = opts['qmax'], opts['nq'], opts['res']
     if opts['is2d']:
-        data = empty_data2D(np.linspace(-qmax, qmax, nq), resolution=res)
+        q = np.linspace(-qmax, qmax, nq)  # type: np.ndarray
+        data = empty_data2D(q, resolution=res)
         data.accuracy = opts['accuracy']
         set_beam_stop(data, 0.0004)
@@ -495 +528 @@
 
 def make_engine(model_info, data, dtype, cutoff):
+    # type: (ModelInfo, Data, str, float) -> Calculator
     """
     Generate the appropriate calculation engine for the given datatype.
@@ -509 +543 @@
 
 def _show_invalid(data, theory):
+    # type: (Data, np.ma.ndarray) -> None
+    """
+    Display a list of the non-finite values in theory.
+    """
     if not theory.mask.any():
         return
@@ -514 +552 @@
     if hasattr(data, 'x'):
         bad = zip(data.x[theory.mask], theory[theory.mask])
-        print("   *** ", ", ".join("I(%g)=%g"%(x, y) for x,y in bad))
+        print("   *** ", ", ".join("I(%g)=%g"%(x, y) for x, y in bad))
 
 
 def compare(opts, limits=None):
+    # type: (Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float]
     """
     Preform a comparison using options from the command line.
@@ -532 +571 @@
     data = opts['data']
 
+    # silence the linter
+    base = opts['engines'][0] if Nbase else None
+    comp = opts['engines'][1] if Ncomp else None
+    base_time = comp_time = None
+    base_value = comp_value = resid = relerr = None
+
     # Base calculation
     if Nbase > 0:
-        base = opts['engines'][0]
         try:
-            base_value, base_time = time_calculation(base, pars, Nbase)
-            base_value = np.ma.masked_invalid(base_value)
+            base_raw, base_time = time_calculation(base, pars, Nbase)
+            base_value = np.ma.masked_invalid(base_raw)
             print("%s t=%.2f ms, intensity=%.0f"
                   % (base.engine, base_time, base_value.sum()))
@@ -547 +591 @@
     # Comparison calculation
     if Ncomp > 0:
-        comp = opts['engines'][1]
         try:
-            comp_value, comp_time = time_calculation(comp, pars, Ncomp)
-            comp_value = np.ma.masked_invalid(comp_value)
+            comp_raw, comp_time = time_calculation(comp, pars, Ncomp)
+            comp_value = np.ma.masked_invalid(comp_raw)
             print("%s t=%.2f ms, intensity=%.0f"
                   % (comp.engine, comp_time, comp_value.sum()))
@@ -625 +668 @@
 
 def _print_stats(label, err):
+    # type: (str, np.ma.ndarray) -> None
+    # work with trimmed data, not the full set
     sorted_err = np.sort(abs(err.compressed()))
-    p50 = int((len(err)-1)*0.50)
-    p98 = int((len(err)-1)*0.98)
+    p50 = int((len(sorted_err)-1)*0.50)
+    p98 = int((len(sorted_err)-1)*0.98)
     data = [
         "max:%.3e"%sorted_err[-1],
         "median:%.3e"%sorted_err[p50],
         "98%%:%.3e"%sorted_err[p98],
-        "rms:%.3e"%np.sqrt(np.mean(err**2)),
-        "zero-offset:%+.3e"%np.mean(err),
+        "rms:%.3e"%np.sqrt(np.mean(sorted_err**2)),
+        "zero-offset:%+.3e"%np.mean(sorted_err),
         ]
     print(label+"  "+"  ".join(data))
@@ -662 +707 @@
 
 def columnize(L, indent="", width=79):
+    # type: (List[str], str, int) -> str
     """
     Format a list of strings into columns.
@@ -679 +725 @@
 
 def get_pars(model_info, use_demo=False):
+    # type: (ModelInfo, bool) -> ParameterSet
     """
     Extract demo parameters from the model definition.
@@ -704 +751 @@
 
 def parse_opts():
+    # type: () -> Dict[str, Any]
     """
     Parse command line options.
@@ -757 +805 @@
         'explore'   : False,
         'use_demo'  : True,
+        'zero'      : False,
     }
     engines = []
@@ -777 +826 @@
         elif arg.startswith('-cutoff='):   opts['cutoff'] = float(arg[8:])
         elif arg.startswith('-random='):   opts['seed'] = int(arg[8:])
-        elif arg == '-random':  opts['seed'] = np.random.randint(1e6)
+        elif arg == '-random':  opts['seed'] = np.random.randint(1000000)
         elif arg == '-preset':  opts['seed'] = -1
         elif arg == '-mono':    opts['mono'] = True
@@ -874 +923 @@
 
 def explore(opts):
+    # type: (Dict[str, Any]) -> None
     """
     Explore the model using the Bumps GUI.
@@ -900 +950 @@
     """
     def __init__(self, opts):
+        # type: (Dict[str, Any]) -> None
         from bumps.cli import config_matplotlib  # type: ignore
         from . import bumps_model
@@ -923 +974 @@
 
     def numpoints(self):
+        # type: () -> int
         """
         Return the number of points.
@@ -929 +981 @@
 
     def parameters(self):
+        # type: () -> Any   # Dict/List hierarchy of parameters
         """
         Return a dictionary of parameters.
@@ -935 +988 @@
 
     def nllf(self):
+        # type: () -> float
         """
         Return cost.
@@ -942 +996 @@
 
     def plot(self, view='log'):
+        # type: (str) -> None
         """
         Plot the data and residuals.
@@ -951 +1006 @@
         if self.limits is None:
             vmin, vmax = limits
-            vmax = 1.3*vmax
-            vmin = vmax*1e-7
-            self.limits = vmin, vmax
+            self.limits = vmax*1e-7, 1.3*vmax
 
 
 def main():
+    # type: () -> None
     """
     Main program.

sasmodels/core.py

@@ -28 +28 @@
 try:
     from typing import List, Union, Optional, Any
-    DType = Union[None, str, np.dtype]
     from .kernel import KernelModel
+    from .modelinfo import ModelInfo
 except ImportError:
     pass
@@ -56 +56 @@
     return available_models
 
-def isstr(s):
-    # type: (Any) -> bool
-    """
-    Return True if *s* is a string-like object.
-    """
-    try: s + ''
-    except Exception: return False
-    return True
-
 def load_model(model_name, dtype=None, platform='ocl'):
-    # type: (str, DType, str) -> KernelModel
+    # type: (str, str, str) -> KernelModel
     """
     Load model info and build model.
@@ -102 +93 @@
 
 def build_model(model_info, dtype=None, platform="ocl"):
-    # type: (modelinfo.ModelInfo, np.dtype, str) -> KernelModel
+    # type: (modelinfo.ModelInfo, str, str) -> KernelModel
     """
     Prepare the model for the default execution platform.
@@ -113 +104 @@
 
     *dtype* indicates whether the model should use single or double precision
-    for the calculation. Any valid numpy single or double precision identifier
-    is valid, such as 'single', 'f', 'f32', or np.float32 for single, or
-    'double', 'd', 'f64'  and np.float64 for double.  If *None*, then use
-    'single' unless the model defines single=False.
+    for the calculation.  Choices are 'single', 'double', 'quad', 'half',
+    or 'fast'.  If *dtype* ends with '!', then force the use of the DLL rather
+    than OpenCL for the calculation.
 
     *platform* should be "dll" to force the dll to be used for C models,
@@ -147 +137 @@
     # source = open(model_info.name+'.cl','r').read()
     source = generate.make_source(model_info)
-    if dtype is None:
-        dtype = generate.F32 if model_info.single else generate.F64
+    numpy_dtype, fast = parse_dtype(model_info, dtype)
     if (platform == "dll"
+            or dtype.endswith('!')
             or not HAVE_OPENCL
-            or not kernelcl.environment().has_type(dtype)):
-        return kerneldll.load_dll(source, model_info, dtype)
+            or not kernelcl.environment().has_type(numpy_dtype)):
+        return kerneldll.load_dll(source, model_info, numpy_dtype)
     else:
-        return kernelcl.GpuModel(source, model_info, dtype)
+        return kernelcl.GpuModel(source, model_info, numpy_dtype, fast=fast)
 
 def precompile_dll(model_name, dtype="double"):
-    # type: (str, DType) -> Optional[str]
+    # type: (str, str) -> Optional[str]
     """
     Precompile the dll for a model.
@@ -172 +162 @@
     """
     model_info = load_model_info(model_name)
+    numpy_dtype, fast = parse_dtype(model_info, dtype)
     source = generate.make_source(model_info)
-    return kerneldll.make_dll(source, model_info, dtype=dtype) if source else None
+    return kerneldll.make_dll(source, model_info, dtype=numpy_dtype) if source else None
+
+def parse_dtype(model_info, dtype):
+    # type: (ModelInfo, str) -> Tuple[np.dtype, bool]
+    """
+    Interpret dtype string, returning np.dtype and fast flag.
+
+    Possible types include 'half', 'single', 'double' and 'quad'.  If the
+    type is 'fast', then this is equivalent to dtype 'single' with the
+    fast flag set to True.
+    """
+    # Fill in default type based on required precision in the model
+    if dtype is None:
+        dtype = 'single' if model_info.single else 'double'
+
+    # Ignore platform indicator
+    if dtype.endswith('!'):
+        dtype = dtype[:-1]
+
+    # Convert type string to type
+    if dtype == 'quad':
+        return generate.F128, False
+    elif dtype == 'half':
+        return generate.F16, False
+    elif dtype == 'fast':
+        return generate.F32, True
+    else:
+        return np.dtype(dtype), False
+

sasmodels/kernelcl.py

@@ -103 +103 @@
 ENV = None
 def environment():
+    # type: () -> "GpuEnvironment"
     """
     Returns a singleton :class:`GpuEnvironment`.
@@ -114 +115 @@
 
 def has_type(device, dtype):
+    # type: (cl.Device, np.dtype) -> bool
     """
     Return true if device supports the requested precision.
@@ -127 +129 @@
 
 def get_warp(kernel, queue):
+    # type: (cl.Kernel, cl.CommandQueue) -> int
     """
     Return the size of an execution batch for *kernel* running on *queue*.
@@ -135 +138 @@
 
 def _stretch_input(vector, dtype, extra=1e-3, boundary=32):
+    # type: (np.ndarray, np.dtype, float, int) -> np.ndarray
     """
     Stretch an input vector to the correct boundary.
@@ -157 +161 @@
 
 def compile_model(context, source, dtype, fast=False):
+    # type: (cl.Context, str, np.dtype, bool) -> cl.Program
     """
     Build a model to run on the gpu.
@@ -192 +197 @@
     """
     def __init__(self):
+        # type: () -> None
         # find gpu context
         #self.context = cl.create_some_context()
@@ -210 +216 @@
 
     def has_type(self, dtype):
+        # type: (np.dtype) -> bool
         """
         Return True if all devices support a given type.
         """
-        dtype = generate.F32 if dtype == 'fast' else np.dtype(dtype)
         return any(has_type(d, dtype)
                    for context in self.context
@@ -219 +225 @@
 
     def get_queue(self, dtype):
+        # type: (np.dtype) -> cl.CommandQueue
         """
         Return a command queue for the kernels of type dtype.
@@ -227 +234 @@
 
     def get_context(self, dtype):
+        # type: (np.dtype) -> cl.Context
         """
         Return a OpenCL context for the kernels of type dtype.
@@ -235 +243 @@
 
     def _create_some_context(self):
+        # type: () -> cl.Context
         """
         Protected call to cl.create_some_context without interactivity.  Use
@@ -248 +257 @@
 
     def compile_program(self, name, source, dtype, fast=False):
+        # type: (str, str, np.dtype, bool) -> cl.Program
         """
         Compile the program for the device in the given context.
@@ -261 +271 @@
 
     def release_program(self, name):
+        # type: (str) -> None
         """
         Free memory associated with the program on the device.
@@ -269 +280 @@
 
 def _get_default_context():
+    # type: () -> cl.Context
     """
     Get an OpenCL context, preferring GPU over CPU, and preferring Intel
@@ -334 +346 @@
     that the compiler is allowed to take shortcuts.
     """
-    def __init__(self, source, model_info, dtype=generate.F32):
+    def __init__(self, source, model_info, dtype=generate.F32, fast=False):
+        # type: (str, ModelInfo, np.dtype, bool) -> None
         self.info = model_info
         self.source = source
-        self.dtype = generate.F32 if dtype == 'fast' else np.dtype(dtype)
-        self.fast = (dtype == 'fast')
+        self.dtype = dtype
+        self.fast = fast
         self.program = None # delay program creation
 
     def __getstate__(self):
+        # type: () -> Tuple[ModelInfo, str, np.dtype, bool]
         return self.info, self.source, self.dtype, self.fast
 
     def __setstate__(self, state):
+        # type: (Tuple[ModelInfo, str, np.dtype, bool]) -> None
         self.info, self.source, self.dtype, self.fast = state
         self.program = None
 
     def make_kernel(self, q_vectors):
+        # type: (List[np.ndarray]) -> "GpuKernel"
         if self.program is None:
             compiler = environment().compile_program
@@ -356 +372 @@
         kernel_name = generate.kernel_name(self.info, is_2d)
         kernel = getattr(self.program, kernel_name)
-        return GpuKernel(kernel, self.info, q_vectors, self.dtype)
+        return GpuKernel(kernel, self.info, q_vectors)
 
     def release(self):
+        # type: () -> None
         """
         Free the resources associated with the model.
@@ -367 +384 @@
 
     def __del__(self):
+        # type: () -> None
         self.release()
 
@@ -390 +408 @@
     """
     def __init__(self, q_vectors, dtype=generate.F32):
+        # type: (List[np.ndarray], np.dtype) -> None
         # TODO: do we ever need double precision q?
         env = environment()
@@ -419 +438 @@
 
     def release(self):
+        # type: () -> None
         """
         Free the memory.
@@ -427 +447 @@
 
     def __del__(self):
+        # type: () -> None
         self.release()
 
@@ -450 +471 @@
     """
     def __init__(self, kernel, model_info, q_vectors):
-        # type: (KernelModel, ModelInfo, List[np.ndarray]) -> None
+        # type: (cl.Kernel, ModelInfo, List[np.ndarray]) -> None
         max_pd = model_info.parameters.max_pd
         npars = len(model_info.parameters.kernel_parameters)-2
@@ -505 +526 @@
 
     def release(self):
+        # type: () -> None
         """
         Release resources associated with the kernel.
@@ -513 +535 @@
 
     def __del__(self):
+        # type: () -> None
         self.release()