Changeset bb39b4a in sasmodels

Timestamp:

Aug 6, 2017 11:36:58 AM (7 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

97d89af

Parents:

232bb12

Message:

make the sascomp command line more regular

Location:

sasmodels

Files:

• compare.py (modified) (20 diffs)
• core.py (modified) (1 diff)
• models/ellipsoid.py (modified) (2 diffs)

sasmodels/compare.py

@@ -56 +56 @@
 
 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.
+usage: sascomp model [options...] [key=val]
+
+Generate and compare SAS models.  If a single model is specified it shows
+a plot of that model.  Different models can be compared, or the same model
+with different parameters.  The same model with the same parameters can
+be compared with different calculation engines to see the effects of precision
+on the resultant values.
 
 model or model1,model2 are the names of the models 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
+    === data generation ===
+    -data="path" uses q, dq from the data file
+    -noise=0 sets the measurement error dI/I
+    -res=0 sets the resolution width dQ/Q if calculating with resolution
     -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
     -zero indicates that q=0 should be included
-    -1d*/-2d computes 1d or 2d data
+
+    === model parameters ===
     -preset*/-random[=seed] preset or random parameters
+    -sets=n generates n random datasets, with the seed given by -random=seed
+    -pars/-nopars* prints the parameter set or not
+    -default/-demo* use demo vs default parameters
+
+    === calculation options ===
     -mono*/-poly force monodisperse or allow polydisperse demo parameters
+    -cutoff=1e-5* cutoff value for including a point in polydispersity
     -magnetic/-nonmagnetic* suppress magnetism
-    -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
-    -default/-demo* use demo vs default parameters
-    -help/-html shows the model docs instead of running the model
-    -title="note" adds note to the plot title, after the model name
-    -data="path" uses q, dq from the data file
-
-Any two calculation engines can be selected for comparison:
-
+
+    === precision options ===
+    -calc=default uses the default calcution precision
     -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 -double.  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).
+    === plotting ===
+    -plot*/-noplot plots or suppress the plot of the model
+    -linear/-log*/-q4 intensity scaling on plots
+    -hist/-nohist* plot histogram of relative error
+    -abs/-rel* plot relative or absolute error
+    -title="note" adds note to the plot title, after the model name
+
+    === output options ===
+    -edit starts the parameter explorer
+    -help/-html shows the model docs instead of running the model
+
+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).
+On unix and mac you may need single quotes around the DLL computation
+engines, such as -calc='single!,double!' since !, is treated as a history
+expansion request in the shell.
 
 Key=value pairs allow you to set specific values for the model parameters.
-Key=value1,value2 to compare different values of the same parameter.
-value can be an expression including other parameters
+Key=value1,value2 to compare different values of the same parameter. The
+value can be an expression including other parameters.
+
+Items later on the command line override those that appear earlier.
+
+Examples:
+
+    # compare single and double precision calculation for a barbell
+    sascomp barbell -calc=single,double
+
+    # generate 10 random lorentz models, with seed=27
+    sascomp lorentz -sets=10 -seed=27
+
+    # compare ellipsoid with R = R_polar = R_equatorial to sphere of radius R
+    sascomp sphere,ellipsoid radius_polar=radius radius_equatorial=radius
+
+    # model timing test requires multiple evals to perform the estimate
+    sascomp pringle -calc=single,double -timing=100,100 -noplot
 """
 
@@ -111 +141 @@
 -------------------
 
-"""
-           + USAGE)
+""" + USAGE)
 
 kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True
@@ -775 +804 @@
     # type: (Dict[str, Any]) -> Dict[str, Any]
 
-    n_base, n_comp = opts['count']
-    pars, pars2 = opts['pars']
+    base, comp = opts['engines']
+    base_n, comp_n = opts['count']
+    base_pars, comp_pars = opts['pars']
     data = opts['data']
 
-    # silence the linter
-    base = opts['engines'][0] if n_base else None
-    comp = opts['engines'][1] if n_comp else None
+    comparison = comp is not None
 
     base_time = comp_time = None
@@ -787 +815 @@
 
     # Base calculation
-    if n_base > 0:
+    try:
+        base_raw, base_time = time_calculation(base, base_pars, base_n)
+        base_value = np.ma.masked_invalid(base_raw)
+        if verbose:
+            print("%s t=%.2f ms, intensity=%.0f"
+                  % (base.engine, base_time, base_value.sum()))
+        _show_invalid(data, base_value)
+    except ImportError:
+        traceback.print_exc()
+
+    # Comparison calculation
+    if comparison:
         try:
-            base_raw, base_time = time_calculation(base, pars, n_base)
-            base_value = np.ma.masked_invalid(base_raw)
-            if verbose:
-                print("%s t=%.2f ms, intensity=%.0f"
-                      % (base.engine, base_time, base_value.sum()))
-            _show_invalid(data, base_value)
-        except ImportError:
-            traceback.print_exc()
-            n_base = 0
-
-    # Comparison calculation
-    if n_comp > 0:
-        try:
-            comp_raw, comp_time = time_calculation(comp, pars2, n_comp)
+            comp_raw, comp_time = time_calculation(comp, comp_pars, comp_n)
             comp_value = np.ma.masked_invalid(comp_raw)
             if verbose:
@@ -810 +836 @@
         except ImportError:
             traceback.print_exc()
-            n_comp = 0
 
     # Compare, but only if computing both forms
-    if n_base > 0 and n_comp > 0:
+    if comparison:
         resid = (base_value - comp_value)
         relerr = resid/np.where(comp_value != 0., abs(comp_value), 1.0)
@@ -856 +881 @@
 
     have_base, have_comp = (base_value is not None), (comp_value is not None)
-    base = opts['engines'][0] if have_base else None
-    comp = opts['engines'][1] if have_comp else None
+    base, comp = opts['engines']
     data = opts['data']
     use_data = (opts['datafile'] is not None) and (have_base ^ have_comp)
@@ -874 +898 @@
 
     if have_base:
-        if have_comp: plt.subplot(131)
+        if have_comp:
+            plt.subplot(131)
         plot_theory(data, base_value, view=view, use_data=use_data, limits=limits)
         plt.title("%s t=%.2f ms"%(base.engine, base_time))
         #cbar_title = "log I"
     if have_comp:
-        if have_base: plt.subplot(132)
+        if have_base:
+            plt.subplot(132)
         if not opts['is2d'] and have_base:
             plot_theory(data, base_value, view=view, use_data=use_data, limits=limits)
@@ -894 +920 @@
             sorted = np.sort(err.flatten())
             cutoff = sorted[int(sorted.size*0.95)]
-            err[err>cutoff] = cutoff
+            err[err > cutoff] = cutoff
         #err,errstr = base/comp,"ratio"
         plot_theory(data, None, resid=err, view=errview, use_data=use_data)
@@ -923 +949 @@
 # ===========================================================================
 #
-NAME_OPTIONS = set([
+
+# Set of command line options.
+# Normal options such as -plot/-noplot are specified as 'name'.
+# For options such as -nq=500 which require a value use 'name='.
+#
+OPTIONS = [
+    # Plotting
     'plot', 'noplot',
-    'half', 'fast', 'single', 'double',
-    'single!', 'double!', 'quad!', 'sasview',
-    'lowq', 'midq', 'highq', 'exq', 'zero',
+    'linear', 'log', 'q4',
+    'rel', 'abs',
+    'hist', 'nohist',
+    'title=',
+
+    # Data generation
+    'data=', 'noise=', 'res=',
+    'nq=', 'lowq', 'midq', 'highq', 'exq', 'zero',
     '2d', '1d',
-    'preset', 'random',
-    'poly', 'mono',
+
+    # Parameter set
+    'preset', 'random', 'random=', 'sets=',
+    'demo', 'default',  # TODO: remove demo/default
+    'nopars', 'pars',
+
+    # Calculation options
+    'poly', 'mono', 'cutoff=',
     'magnetic', 'nonmagnetic',
-    'nopars', 'pars',
-    'rel', 'abs',
-    'linear', 'log', 'q4',
-    'hist', 'nohist',
-    'edit', 'html', 'help',
-    'demo', 'default',
-    ])
-VALUE_OPTIONS = [
-    # Note: random is both a name option and a value option
-    'cutoff', 'random', 'nq', 'res', 'accuracy', 'title', 'data', 'sets'
+    'accuracy=',
+
+    # Precision options
+    'calc=',
+    'half', 'fast', 'single', 'double', 'single!', 'double!', 'quad!',
+    'sasview',  # TODO: remove sasview 3.x support
+    'timing=',
+
+    # Output options
+    'help', 'html', 'edit',
     ]
+
+NAME_OPTIONS = set(k for k in OPTIONS if not k.endswith('='))
+VALUE_OPTIONS = [k[:-1] for k in OPTIONS if k.endswith('=')]
+
 
 def columnize(items, indent="", width=79):
@@ -1002 +1049 @@
 # not sure about brackets [] {}
 # maybe one of the following @ ? ^ ! ,
-MODEL_SPLIT = ','
+PAR_SPLIT = ','
 def parse_opts(argv):
     # type: (List[str]) -> Dict[str, Any]
@@ -1021 +1068 @@
         print(columnize(MODELS, indent="  "))
         return None
-    if len(positional_args) > 3:
-        print("expected parameters: model N1 N2")
 
     invalid = [o[1:] for o in flags
@@ -1031 +1076 @@
         return None
 
-    name = positional_args[0]
-    n1 = int(positional_args[1]) if len(positional_args) > 1 else 1
-    n2 = int(positional_args[2]) if len(positional_args) > 2 else 1
+    name = positional_args[-1]
 
     # pylint: disable=bad-whitespace
@@ -1044 +1087 @@
         'nq'        : 128,
         'res'       : 0.0,
+        'noise'     : 0.0,
         'accuracy'  : 'Low',
-        'cutoff'    : 0.0,
+        'cutoff'    : '0.0',
         'seed'      : -1,  # default to preset
         'mono'      : True,
@@ -1060 +1104 @@
         'datafile'  : None,
         'sets'      : 1,
+        'engine'    : 'default',
+        'evals'     : '1',
     }
-    engines = []
     for arg in flags:
         if arg == '-noplot':    opts['plot'] = False
@@ -1077 +1122 @@
         elif arg.startswith('-nq='):       opts['nq'] = int(arg[4:])
         elif arg.startswith('-res='):      opts['res'] = float(arg[5:])
+        elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:])
         elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:])
         elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:]
-        elif arg.startswith('-cutoff='):   opts['cutoff'] = float(arg[8:])
+        elif arg.startswith('-cutoff='):   opts['cutoff'] = arg[8:]
         elif arg.startswith('-random='):   opts['seed'] = int(arg[8:])
         elif arg.startswith('-title='):    opts['title'] = arg[7:]
         elif arg.startswith('-data='):     opts['datafile'] = arg[6:]
+        elif arg.startswith('-calc='):     opts['engine'] = arg[6:]
+        elif arg.startswith('-neval='):    opts['evals'] = arg[7:]
         elif arg == '-random':  opts['seed'] = np.random.randint(1000000)
         elif arg == '-preset':  opts['seed'] = -1
@@ -1095 +1143 @@
         elif arg == '-rel':     opts['rel_err'] = True
         elif arg == '-abs':     opts['rel_err'] = False
-        elif arg == '-half':    engines.append(arg[1:])
-        elif arg == '-fast':    engines.append(arg[1:])
-        elif arg == '-single':  engines.append(arg[1:])
-        elif arg == '-double':  engines.append(arg[1:])
-        elif arg == '-single!': engines.append(arg[1:])
-        elif arg == '-double!': engines.append(arg[1:])
-        elif arg == '-quad!':   engines.append(arg[1:])
-        elif arg == '-sasview': engines.append(arg[1:])
+        elif arg == '-half':    opts['engine'] = 'half'
+        elif arg == '-fast':    opts['engine'] = 'fast'
+        elif arg == '-single':  opts['engine'] = 'single'
+        elif arg == '-double':  opts['engine'] = 'double'
+        elif arg == '-single!': opts['engine'] = 'single!'
+        elif arg == '-double!': opts['engine'] = 'double!'
+        elif arg == '-quad!':   opts['engine'] = 'quad!'
+        elif arg == '-sasview': opts['engine'] = 'sasview'
         elif arg == '-edit':    opts['explore'] = True
         elif arg == '-demo':    opts['use_demo'] = True
@@ -1114 +1162 @@
         opts['seed'] = np.random.randint(1000000)
 
-    if MODEL_SPLIT in name:
-        name, name2 = name.split(MODEL_SPLIT, 2)
-    else:
-        name2 = name
-    try:
-        model_info = core.load_model_info(name)
-        model_info2 = core.load_model_info(name2) if name2 != name else model_info
-    except ImportError as exc:
-        print(str(exc))
-        print("Could not find model; use one of:\n    " + models)
-        return None
-
-    # TODO: check if presets are different when deciding if models are same
-    same_model = name == name2
-    if len(engines) == 0:
-        if same_model:
-            engines.extend(['single', 'double'])
-        else:
-            engines.extend(['single', 'single'])
-    elif len(engines) == 1:
-        if not same_model:
-            engines.append(engines[0])
-        elif engines[0] == 'double':
-            engines.append('single')
-        else:
-            engines.append('double')
-    elif len(engines) > 2:
-        del engines[2:]
-
     # Create the computational engines
     if opts['datafile'] is not None:
@@ -1148 +1167 @@
     else:
         data, _ = make_data(opts)
-    if n1:
-        base = make_engine(model_info, data, engines[0], opts['cutoff'])
+
+    comparison = any(PAR_SPLIT in v for v in values)
+    if PAR_SPLIT in name:
+        names = name.split(PAR_SPLIT, 2)
+        comparison = True
     else:
-        base = None
-    if n2:
-        comp = make_engine(model_info2, data, engines[1], opts['cutoff'])
+        names = [name]*2
+    try:
+        model_info = [core.load_model_info(k) for k in names]
+    except ImportError as exc:
+        print(str(exc))
+        print("Could not find model; use one of:\n    " + models)
+        return None
+
+    if PAR_SPLIT in opts['engine']:
+        engine_types = opts['engine'].split(PAR_SPLIT, 2)
+        comparison = True
+    else:
+        engine_types = [opts['engine']]*2
+
+    if PAR_SPLIT in opts['evals']:
+        evals = [int(k) for k in opts['evals'].split(PAR_SPLIT, 2)]
+        comparison = True
+    else:
+        evals = [int(opts['evals'])]*2
+
+    if PAR_SPLIT in opts['cutoff']:
+        cutoff = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)]
+        comparison = True
+    else:
+        cutoff = [float(opts['cutoff'])]*2
+
+    base = make_engine(model_info[0], data, engine_types[0], cutoff[0])
+    if comparison:
+        comp = make_engine(model_info[1], data, engine_types[1], cutoff[1])
     else:
         comp = None
@@ -1161 +1209 @@
     opts.update({
         'data'      : data,
-        'name'      : [name, name2],
-        'def'       : [model_info, model_info2],
-        'count'     : [n1, n2],
+        'name'      : names,
+        'def'       : model_info,
+        'count'     : evals,
         'engines'   : [base, comp],
         'values'    : values,
@@ -1210 +1258 @@
             print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s))))
             return None
-        v1, v2 = v.split(MODEL_SPLIT, 2) if MODEL_SPLIT in v else (v,v)
+        v1, v2 = v.split(PAR_SPLIT, 2) if PAR_SPLIT in v else (v,v)
         if v1 and k in pars:
             presets[k] = float(v1) if isnumber(v1) else v1

sasmodels/core.py

@@ -253 +253 @@
 
     # Convert dtype string to numpy dtype.
-    if dtype is None:
+    if dtype is None or dtype == "default":
         numpy_dtype = (generate.F32 if platform == "ocl" and model_info.single
                        else generate.F64)

sasmodels/models/ellipsoid.py

@@ -111 +111 @@
 *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*,
 Plenum Press, New York, 1987.
+
+A. Isihara. J. Chem. Phys. 18(1950) 1446-1449
 
 Authorship and Verification
@@ -162 +164 @@
 def ER(radius_polar, radius_equatorial):
     import numpy as np
-# see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
+    # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
     ee = np.empty_like(radius_polar)
     idx = radius_polar > radius_equatorial