-                      reb46451
+                      rcaeb06d
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
+"""
+Program to compare models using different compute engines.
+This program lets you compare results between OpenCL and DLL versions
+of the code and between precision (half, fast, single, double, quad),
+where fast precision is single precision using native functions for
+trig, etc., and may not be completely IEEE 754 compliant.  This lets
+make sure that the model calculations are stable, or if you need to
+tag the model as double precision only.
+Run using ./compare.sh (Linux, Mac) or compare.bat (Windows) in the
+sasmodels root to see the command line options.
+Note that there is no way within sasmodels to select between an
+OpenCL CPU device and a GPU device, but you can do so by setting the
+PYOPENCL_CTX environment variable ahead of time.  Start a python
+interpreter and enter::
+    import pyopencl as cl
+    cl.create_some_context()
+This will prompt you to select from the available OpenCL devices
+and tell you which string to use for the PYOPENCL_CTX variable.
+On Windows you will need to remove the quotes.
+"""
+from __future__ import print_function
+USAGE = """
+usage: compare.py model N1 N2 [options...] [key=val]
+Compare the speed and value for a model between the SasView original and the
+sasmodels rewrite.
+model is the name of the model to compare (see below).
+N1 is the number of times to run sasmodels (default=1).
+N2 is the number times to run sasview (default=1).
+Options (* for default):
+    -plot*/-noplot plots or suppress the plot of the model
+    -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0
+    -nq=128 sets the number of Q points in the data set
+    -1d*/-2d computes 1d or 2d data
+    -preset*/-random[=seed] preset or random parameters
+    -mono/-poly* force monodisperse/polydisperse
+    -cutoff=1e-5* cutoff value for including a point in polydispersity
+    -pars/-nopars* prints the parameter set or not
+    -abs/-rel* plot relative or absolute error
+    -linear/-log*/-q4 intensity scaling
+    -hist/-nohist* plot histogram of relative error
+    -res=0 sets the resolution width dQ/Q if calculating with resolution
+    -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh
+    -edit starts the parameter explorer
+Any two calculation engines can be selected for comparison:
+    -single/-double/-half/-fast sets an OpenCL calculation engine
+    -single!/-double!/-quad! sets an OpenMP calculation engine
+    -sasview sets the sasview calculation engine
+The default is -single -sasview.  Note that the interpretation of quad
+precision depends on architecture, and may vary from 64-bit to 128-bit,
+with 80-bit floats being common (1e-19 precision).
+Key=value pairs allow you to set specific values for the model parameters.
+"""
+# Update docs with command line usage string.   This is separate from the usual
+# doc string so that we can display it at run time if there is an error.
+# lin
+__doc__ = __doc__ + """
+Program description
+-------------------
+""" + USAGE
 import sys
 …
 # List of available models
 MODELS = [basename(f)[:-3]
           for f in sorted(glob.glob(joinpath(ROOT,"models","[a-zA-Z]*.py")))]
+          for f in sorted(glob.glob(joinpath(ROOT, "models", "[a-zA-Z]*.py")))]
 # CRUFT python 2.6
 …
     elif p.endswith('_pd_type'):
         return v
     elif any(s in p for s in ('theta','phi','psi')):
+    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 [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(p, v):
     """
     Randomizing parameter.
     """
     if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')):
+    Randomize a single parameter.
+    """
+    if any(p.endswith(s) for s in ('_pd_n', '_pd_nsigma', '_pd_type')):
         return v
     else:
 …
 def randomize_pars(pars, seed=None):
+    """
+    Generate random values for all of the parameters.
+    Valid ranges for the random number generator are guessed from the name of
+    the parameter; this will not account for constraints such as cap radius
+    greater than cylinder radius in the capped_cylinder model, so
+    :func:`constrain_pars` needs to be called afterward..
+    """
     np.random.seed(seed)
     # Note: the sort guarantees order `of calls to random number generator
     pars = dict((p,_randomize_one(p,v))
                 for p,v in sorted(pars.items()))
+    pars = dict((p, _randomize_one(p, v))
+                for p, v in sorted(pars.items()))
     return pars
 …
     """
     Restrict parameters to valid values.
+    This includes model specific code for models such as capped_cylinder
+    which need to support within model constraints (cap radius more than
+    cylinder radius in this case).
     """
     name = model_definition.name
     if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']:
         pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius']
+        pars['radius'], pars['cap_radius'] = pars['cap_radius'], pars['radius']
     if name == 'barbell' and pars['bell_radius'] < pars['radius']:
         pars['radius'],pars['bell_radius'] = pars['bell_radius'],pars['radius']
+        pars['radius'], pars['bell_radius'] = pars['bell_radius'], pars['radius']
     # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff)
 …
         q_max = 1.0  # high q maximum
         rg_max = np.sqrt(90*np.log(10) + 3*np.log(pars['scale']))/q_max
         pars['rg'] = min(pars['rg'],rg_max)
+        pars['rg'] = min(pars['rg'], rg_max)
     if name == 'rpa':
 …
 def parlist(pars):
+    return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items()))
+    """
+    Format the parameter list for printing.
+    """
+    return "\n".join("%s: %s"%(p, v) for p, v in sorted(pars.items()))
 def suppress_pd(pars):
 …
 def eval_sasview(model_definition, data):
+    """
+    Return a model calculator using the SasView fitting engine.
+    """
     # importing sas here so that the error message will be that sas failed to
     # import rather than the more obscure smear_selection not imported error
 …
     # convert model parameters from sasmodel form to sasview form
     #print("old",sorted(pars.items()))
     modelname, pars = revert_model(model_definition, {})
     #print("new",sorted(pars.items()))
+    modelname, _pars = revert_model(model_definition, {})
+    #print("new",sorted(_pars.items()))
     sas = __import__('sas.models.'+modelname)
     ModelClass = getattr(getattr(sas.models,modelname,None),modelname,None)
+    ModelClass = getattr(getattr(sas.models, modelname, None), modelname, None)
     if ModelClass is None:
         raise ValueError("could not find model %r in sas.models"%modelname)
 …
     def calculator(**pars):
+        """
+        Sasview calculator for model.
+        """
         # paying for parameter conversion each time to keep life simple, if not fast
         _, pars = revert_model(model_definition, pars)
         for k,v in pars.items():
+        for k, v in pars.items():
             parts = k.split('.')  # polydispersity components
             if len(parts) == 2:
 …
+}
 def eval_opencl(model_definition, data, dtype='single', cutoff=0.):
+    """
+    Return a model calculator using the OpenCL calculation engine.
+    """
     try:
         model = core.load_model(model_definition, dtype=dtype, platform="ocl")
 …
 def eval_ctypes(model_definition, data, dtype='double', cutoff=0.):
     if dtype=='quad':
+    if dtype == 'quad':
         dtype = 'longdouble'
     model = core.load_model(model_definition, dtype=dtype, platform="dll")
 …
 def time_calculation(calculator, pars, Nevals=1):
+    """
+    Compute the average calculation time over N evaluations.
+    An additional call is generated without polydispersity in order to
+    initialize the calculation engine, and make the average more stable.
+    """
     # initialize the code so time is more accurate
     value = calculator(**suppress_pd(pars))
 …
 def make_data(opts):
+    """
+    Generate an empty dataset, used with the model to set Q points
+    and resolution.
+    *opts* contains the options, with 'qmax', 'nq', 'res',
+    'accuracy', 'is2d' and 'view' parsed from the command line.
+    """
     qmax, nq, res = opts['qmax'], opts['nq'], opts['res']
     if opts['is2d']:
 …
 def make_engine(model_definition, data, dtype, cutoff):
+    """
+    Generate the appropriate calculation engine for the given datatype.
+    Datatypes with '!' appended are evaluated using external C DLLs rather
+    than OpenCL.
+    """
     if dtype == 'sasview':
         return eval_sasview(model_definition, data)
 …
 def compare(opts, limits=None):
+    """
+    Preform a comparison using options from the command line.
+    *limits* are the limits on the values to use, either to set the y-axis
+    for 1D or to set the colormap scale for 2D.  If None, then they are
+    inferred from the data and returned. When exploring using Bumps,
+    the limits are set when the model is initially called, and maintained
+    as the values are adjusted, making it easier to see the effects of the
+    parameters.
+    """
     Nbase, Ncomp = opts['N1'], opts['N2']
     pars = opts['pars']
 …
         resid = (base_value - comp_value)
         relerr = resid/comp_value
+        _print_stats("|%s - %s|"%(base.engine,comp.engine)+(" "*(3+len(comp.engine))), resid)
+        _print_stats("|(%s - %s) / %s|"%(base.engine,comp.engine,comp.engine), relerr)
+        _print_stats("|%s-%s|"%(base.engine, comp.engine) + (" "*(3+len(comp.engine))),
+                     resid)
+        _print_stats("|(%s-%s)/%s|"%(base.engine, comp.engine, comp.engine),
+                     relerr)
     # Plot if requested
 …
         if Nbase > 0: plt.subplot(132)
         plot_theory(data, comp_value, view=view, plot_data=False, limits=limits)
         plt.title("%s t=%.1f ms"%(comp.engine,comp_time))
+        plt.title("%s t=%.1f ms"%(comp.engine, comp_time))
         #cbar_title = "log I"
     if Ncomp > 0 and Nbase > 0:
         plt.subplot(133)
         if '-abs' in opts:
             err,errstr,errview = resid, "abs err", "linear"
+            err, errstr, errview = resid, "abs err", "linear"
         else:
             err,errstr,errview = abs(relerr), "rel err", "log"
+            err, errstr, errview = abs(relerr), "rel err", "log"
         #err,errstr = base/comp,"ratio"
         plot_theory(data, None, resid=err, view=errview, plot_data=False)
 …
         plt.figure()
         v = relerr
+        v[v==0] = 0.5*np.min(np.abs(v[v!=0]))
+        plt.hist(np.log10(np.abs(v)), normed=1, bins=50);
+        plt.xlabel('log10(err), err = |(%s - %s) / %s|'%(base.engine, comp.engine, comp.engine));
+        v[v == 0] = 0.5*np.min(np.abs(v[v != 0]))
+        plt.hist(np.log10(np.abs(v)), normed=1, bins=50)
+        plt.xlabel('log10(err), err = |(%s - %s) / %s|'
+                   % (base.engine, comp.engine, comp.engine))
         plt.ylabel('P(err)')
         plt.title('Distribution of relative error between calculation engines')
 …
         "zero-offset:%+.3e"%np.mean(err),
+        ]
     print(label+"  ".join(data))
+    print(label+"  "+"  ".join(data))
 …
 # ===========================================================================
+#
-USAGE="""
-usage: compare.py model N1 N2 [options...] [key=val]
-Compare the speed and value for a model between the SasView original and the
-sasmodels rewrite.
-model is the name of the model to compare (see below).
-N1 is the number of times to run sasmodels (default=1).
-N2 is the number times to run sasview (default=1).
-Options (* for default):
-    -plot*/-noplot plots or suppress the plot of the model
-    -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0
-    -nq=128 sets the number of Q points in the data set
-    -1d*/-2d computes 1d or 2d data
-    -preset*/-random[=seed] preset or random parameters
-    -mono/-poly* force monodisperse/polydisperse
-    -cutoff=1e-5* cutoff value for including a point in polydispersity
-    -pars/-nopars* prints the parameter set or not
-    -abs/-rel* plot relative or absolute error
-    -linear/-log*/-q4 intensity scaling
-    -hist/-nohist* plot histogram of relative error
-    -res=0 sets the resolution width dQ/Q if calculating with resolution
-    -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh
-    -edit starts the parameter explorer
-Any two calculation engines can be selected for comparison:
-    -single/-double/-half/-fast sets an OpenCL calculation engine
-    -single!/-double!/-quad! sets an OpenMP calculation engine
-    -sasview sets the sasview calculation engine
-The default is -single -sasview.  Note that the interpretation of quad
-precision depends on architecture, and may vary from 64-bit to 128-bit,
-with 80-bit floats being common (1e-19 precision).
-Key=value pairs allow you to set specific values for the model parameters.
-Available models:
-"""
 NAME_OPTIONS = set([
     'plot', 'noplot',
 …
 def columnize(L, indent="", width=79):
+    """
+    Format a list of strings into columns for printing.
+    """
     column_width = max(len(w) for w in L) + 1
     num_columns = (width - len(indent)) // column_width
 …
 def get_demo_pars(model_definition):
+    """
+    Extract demo parameters from the model definition.
+    """
     info = generate.make_info(model_definition)
     # Get the default values for the parameters
 …
 def parse_opts():
+    flags = [arg for arg in sys.argv[1:] if arg.startswith('-')]
+    values = [arg for arg in sys.argv[1:] if not arg.startswith('-') and '=' in arg]
+    args = [arg for arg in sys.argv[1:] if not arg.startswith('-') and '=' not in arg]
+    """
+    Parse command line options.
+    """
+    flags = [arg for arg in sys.argv[1:]
+             if arg.startswith('-')]
+    values = [arg for arg in sys.argv[1:]
+              if not arg.startswith('-') and '=' in arg]
+    args = [arg for arg in sys.argv[1:]
+            if not arg.startswith('-') and '=' not in arg]
     models = "\n    ".join("%-15s"%v for v in MODELS)
     if len(args) == 0:
         print(USAGE)
+        print("\nAvailable models:")
         print(columnize(MODELS, indent="  "))
         sys.exit(1)
     if args[0] not in MODELS:
         print("Model %r not available. Use one of:\n    %s"%(args[0],models))
+        print("Model %r not available. Use one of:\n    %s"%(args[0], models))
         sys.exit(1)
     if len(args) > 3:
 …
     invalid = [o[1:] for o in flags
                if o[1:] not in NAME_OPTIONS
                   and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)]
+                   and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)]
     if invalid:
         print("Invalid options: %s"%(", ".join(invalid)))

Note: See TracChangeset for help on using the changeset viewer.

SasView

Changeset caeb06d in sasmodels

Legend:

sasmodels/compare.py

Download in other formats: