Changeset dd7fc12 in sasmodels for sasmodels/compare.py

Timestamp:

Apr 15, 2016 11:11:43 AM (8 years ago)

Author:

Paul Kienzle <pkienzle@…>

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:

3599d36

Parents:

b151003

Message:

fix kerneldll dtype problem; more type hinting

File:

• sasmodels/compare.py (modified) (41 diffs)

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.