Changes in / [55b283e8:b968fe8] in sasmodels

Location:

sasmodels

Files:

• compare.py (modified) (34 diffs)
• models/hollow_cylinder.c (modified) (2 diffs)
• models/hollow_cylinder.py (modified) (4 diffs)
• models/lib/Si.c (modified) (3 diffs)

sasmodels/compare.py

@@ -1 +1 @@
 #!/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
@@ -25 +103 @@
 # 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
@@ -76 +154 @@
     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]
@@ -86 +164 @@
     elif p.endswith('_pd'):
         return [0, 1]
-    elif p in ['background', 'scale']:
+    elif p == 'background':
+        return [0, 10]
+    elif p == 'scale':
         return [0, 1e3]
+    elif p == 'case_num':
+        # RPA hack
+        return [0, 10]
+    elif v < 0:
+        # Kxy parameters in rpa model can be negative
+        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:
@@ -101 +187 @@
 
 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
 
@@ -110 +204 @@
     """
     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)
@@ -122 +220 @@
         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':
@@ -136 +234 @@
 
 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):
@@ -151 +252 @@
 
 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
@@ -158 +262 @@
     # 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, _ = 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)
@@ -184 +288 @@
 
     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:
@@ -209 +316 @@
 }
 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")
@@ -221 +331 @@
 
 def eval_ctypes(model_definition, data, dtype='double', cutoff=0.):
-    if dtype=='quad':
+    """
+    Return a model calculator using the DLL calculation engine.
+    """
+    if dtype == 'quad':
         dtype = 'longdouble'
     model = core.load_model(model_definition, dtype=dtype, platform="dll")
@@ -229 +342 @@
 
 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))
@@ -238 +357 @@
 
 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']:
@@ -255 +381 @@
 
 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)
@@ -265 +397 @@
 
 def compare(opts, limits=None):
-    Nbase, Ncomp = opts['N1'], opts['N2']
+    """
+    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']
     data = opts['data']
@@ -296 +438 @@
         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
@@ -321 +465 @@
         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)
@@ -340 +484 @@
         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')
@@ -362 +507 @@
         "zero-offset:%+.3e"%np.mean(err),
         ]
-    print(label+"  ".join(data))
+    print(label+"  "+"  ".join(data))
 
 
@@ -368 +513 @@
 # ===========================================================================
 #
-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',
@@ -431 +533 @@
 
 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
@@ -443 +548 @@
 
 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
-    pars = dict((p[0],p[2]) for p in info['parameters'])
+    pars = dict((p[0], p[2]) for p in info['parameters'])
 
     # Fill in default values for the polydispersity parameters
@@ -460 +568 @@
 
 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:
-        print("expected parameters: model Nopencl Nsasview")
+        print("expected parameters: model N1 N2")
 
     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)))
@@ -538 +653 @@
 
     if len(engines) == 0:
-        engines.extend(['single','sasview'])
+        engines.extend(['single', 'sasview'])
     elif len(engines) == 1:
         if engines[0][0] != 'sasview':
@@ -550 +665 @@
     model_definition = core.load_model_definition(name)
 
-    N1 = int(args[1]) if len(args) > 1 else 1
-    N2 = int(args[2]) if len(args) > 2 else 1
+    n1 = int(args[1]) if len(args) > 1 else 1
+    n2 = int(args[2]) if len(args) > 2 else 1
 
     # Get demo parameters from model definition, or use default parameters
@@ -583 +698 @@
     # Create the computational engines
     data, _index = make_data(opts)
-    if N1:
+    if n1:
         base = make_engine(model_definition, data, engines[0], opts['cutoff'])
     else:
         base = None
-    if N2:
+    if n2:
         comp = make_engine(model_definition, data, engines[1], opts['cutoff'])
     else:
@@ -596 +711 @@
         'name'      : name,
         'def'       : model_definition,
-        'N1'        : N1,
-        'N2'        : N2,
+        'n1'        : n1,
+        'n2'        : n2,
         'presets'   : presets,
         'pars'      : pars,
@@ -635 +750 @@
     def __init__(self, opts):
         from bumps.cli import config_matplotlib
-        import bumps_model
+        from . import bumps_model
         config_matplotlib()
         self.opts = opts
@@ -643 +758 @@
             active = [base + ext
                       for base in info['partype']['pd-1d']
-                      for ext in ['','_pd','_pd_n','_pd_nsigma']]
+                      for ext in ['', '_pd', '_pd_n', '_pd_nsigma']]
             active.extend(info['partype']['fixed-1d'])
             for k in active:
@@ -658 +773 @@
     def numpoints(self):
         """
-        Return the number of points
+        Return the number of points.
         """
         return len(self.pars) + 1  # so dof is 1
@@ -664 +779 @@
     def parameters(self):
         """
-        Return a dictionary of parameters
+        Return a dictionary of parameters.
         """
         return self.pars
 
     def nllf(self):
+        """
+        Return cost.
+        """
         return 0.  # No nllf
 
@@ -675 +793 @@
         Plot the data and residuals.
         """
-        pars = dict((k, v.value) for k,v in self.pars.items())
+        pars = dict((k, v.value) for k, v in self.pars.items())
         pars.update(self.pd_types)
         self.opts['pars'] = pars

sasmodels/models/hollow_cylinder.c

@@ -114 +114 @@
         double norm,volume;     //final calculation variables
         
+        if (core_radius >= radius || radius >= length) {
+                return NAN;
+        }
+        
         delrho = solvent_sld - sld;
         lower = 0.0;
@@ -132 +136 @@
 }
 
-//FIXME: Factor of two difference
 double Iqxy(double qx, double qy, double radius, double core_radius, double length, double sld,
         double solvent_sld, double theta, double phi)

sasmodels/models/hollow_cylinder.py

@@ -64 +64 @@
 """
 
-from numpy import inf
+from numpy import pi, inf
 
 name = "hollow_cylinder"
@@ -92 +92 @@
 source = ["lib/J1.c", "lib/gauss76.c", "hollow_cylinder.c"]
 
+def ER(radius, core_radius, length):
+    if radius == 0 or length == 0:
+        return 0.0
+    len1 = radius
+    len2 = length/2.0
+    term1 = len1*len1*2.0*len2/2.0
+    term2 = 1.0 + (len2/len1)*(1.0 + 1/len2/2.0)*(1.0 + pi*len1/len2/2.0)
+    ddd = 3.0*term1*term2
+    diam = pow(ddd, (1.0/3.0))
+    return diam
+
+def VR(radius, core_radius, length):
+    vol_core = pi*core_radius*core_radius*length
+    vol_total = pi*radius*radius*length
+    vol_shell = vol_total - vol_core
+    return vol_shell, vol_total
+
 # parameters for demo
 demo = dict(scale=1.0,background=0.0,length=400.0,radius=30.0,core_radius=20.0,
@@ -97 +114 @@
             radius_pd=.2, radius_pd_n=9,
             length_pd=.2, length_pd_n=10,
+            core_radius_pd=.2, core_radius_pd_n=9,
             theta_pd=10, theta_pd_n=5,
             )
@@ -109 +127 @@
 # Parameters for unit tests
 tests = [
-         [{"radius" : 30.0},0.00005,1764.926]
+         [{"radius" : 30.0},0.00005,1764.926],
+         [{},'VR',1.8],
+         [{},0.001,1756.76]
          ]

sasmodels/models/lib/Si.c

@@ -1 @@
-int factorial(int f);
 double Si(double x);
 
-// integral of sin(x)/x: approximated to w/i 1%
+// integral of sin(x)/x Taylor series approximated to w/i 0.1%
 double Si(double x)
 {
-        int i;
-        int nmax=6;
         double out;
-        long power;
         double pi = 4.0*atan(1.0);
 
@@ -16 +12 @@
                 out = pi/2.0;
 
-                for (i=0; i<nmax-2; i+=1) {
-                        out_cos += pow(-1.0, i) * (double)factorial(2*i) / pow(x, 2*i+1);
-                        out_sin += pow(-1.0, i) * (double)factorial(2*i+1) / pow(x, 2*i+2);
-                }
+                // Explicitly writing factorial values triples the speed of the calculation
+                out_cos = 1/x - 2/pow(x,3) + 24/pow(x,5) - 720/pow(x,7) + 40320/pow(x,9);
+                out_sin = 1/x - 6/pow(x,4) + 120/pow(x,6) - 5040/pow(x,8) + 362880/pow(x,10);
 
                 out -= cos(x) * out_cos;
@@ -26 +21 @@
         }
 
-        out = 0.0;
+        // Explicitly writing factorial values triples the speed of the calculation
+        out = x - pow(x, 3)/18 + pow(x,5)/600 - pow(x,7)/35280 + pow(x,9)/3265920;
 
-        for (i=0; i<nmax; i+=1) {
-                if (i==0) {
-                        out += x;
-                        continue;
-                }
-
-                power = pow(x,(2 * i + 1));
-                out += pow(-1.0, i) * power / ((2.0 * (double)i + 1.0) * (double)factorial(2 * i + 1));
-
-                //printf ("Si=%g %g %d\n", x, out, i);
-        }
-
+        //printf ("Si=%g %g\n", x, out);
         return out;
 }
-
-int factorial(int f)
-{
-    if ( f == 0 ) 
-        return 1;
-    return(f * factorial(f - 1));
-}