Changes in / [402c351:5442c77] in sasmodels

Location:

sasmodels

Files:

• compare.py (modified) (7 diffs)
• compare_many.py (modified) (5 diffs)
• generate.py (modified) (3 diffs)
• kernel_template.c (modified) (4 diffs)
• kernelcl.py (modified) (1 diff)
• kerneldll.py (modified) (8 diffs)
• kernelpy.py (modified) (3 diffs)
• resolution.py (modified) (1 diff)

sasmodels/compare.py

@@ -130 +130 @@
     if name == 'capped_cylinder' and 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']
+
+    # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff)
+    if name == 'guinier':
+        #q_max = 0.2  # mid q maximum
+        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)
 
     # These constraints are only needed for comparison to sasview
@@ -224 +233 @@
     model_definition = core.load_model_definition(name)
 
-    view = 'linear' if '-linear' in opts else 'log' if '-log' in opts else 'q4' if '-q4' in opts else 'log'
+    view = ('linear' if '-linear' in opts
+            else 'log' if '-log' in opts
+            else 'q4' if '-q4' in opts
+            else 'log')
 
     opt_values = dict(split
@@ -230 +242 @@
                       if len(split) == 2)
     # Sort out data
-    qmax = 10.0 if '-exq' in opts else 1.0 if '-highq' in opts else 0.2 if '-midq' in opts else 0.05
+    qmax = (10.0 if '-exq' in opts
+            else 1.0 if '-highq' in opts
+            else 0.2 if '-midq' in opts
+            else 0.05)
     Nq = int(opt_values.get('-Nq', '128'))
     res = float(opt_values.get('-res', '0'))
@@ -239 +254 @@
 
     # modelling accuracy is determined by dtype and cutoff
-    dtype = 'double' if '-double' in opts else 'single'
+    dtype = ('longdouble' if '-quad' in opts
+             else 'double' if '-double' in opts
+             else 'single')
     cutoff = float(opt_values.get('-cutoff','1e-5'))
 
@@ -247 +264 @@
         seed = int(opt_values['-random']) if '-random' in opt_values else None
         pars, seed = randomize_model(pars, seed=seed)
-        constrain_pars(model_definition, pars)
         print "Randomize using -random=%i"%seed
     pars.update(set_pars)  # set value after random to control value
+    constrain_pars(model_definition, pars)
 
     # parameter selection
@@ -362 +379 @@
 
     -plot*/-noplot plots or suppress the plot of the model
-    -single*/-double uses double precision for comparison
+    -single*/-double/-quad use single/double/quad precision for comparison
     -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
@@ -386 +403 @@
 NAME_OPTIONS = set([
     'plot','noplot',
-    'single','double',
+    'single','double','longdouble',
     'lowq','midq','highq','exq',
     '2d','1d',

sasmodels/compare_many.py

@@ -57 +57 @@
     num_good = 0
     first = True
+    max_diff = 0
     for k in range(N):
         print >>sys.stderr, name, k
@@ -73 +74 @@
             columns.extend(stats)
             labels.append('GPU single')
-            good = good and (stats[0] < 1e-14)
+            max_diff = max(max_diff, stats[0])
+            good = good and (stats[0] < 5e-5)
         if 0 and environment().has_double:
             gpu_double_value = trymodel(eval_opencl, dtype='double', cutoff=cutoff)
@@ -79 +81 @@
             columns.extend(stats)
             labels.append('GPU double')
-            good = good and (stats[0] < 1e-14)
+            max_diff = max(max_diff, stats[0])
+            good = good and (stats[0] < 1e-12)
         if 1:
             cpu_double_value = trymodel(eval_ctypes, dtype='double', cutoff=cutoff)
@@ -85 +88 @@
             columns.extend(stats)
             labels.append('CPU double')
-            good = good and (stats[0] < 1e-14)
+            max_diff = max(max_diff, stats[0])
+            good = good and (stats[0] < 1e-12)
         if 0:
             stats = get_stats(cpu_double_value, gpu_single_value, index)
             columns.extend(stats)
             labels.append('single/double')
-            good = good and (stats[0] < 1e-14)
+            max_diff = max(max_diff, stats[0])
+            good = good and (stats[0] < 5e-5)
 
         columns += [v for _,v in sorted(pars_i.items())]
@@ -100 +105 @@
         else:
             print(("%d,"%seed)+','.join("%g"%v for v in columns))
-    print '"%d/%d good"'%(num_good, N)
+    print '"good","%d/%d","max diff",%g'%(num_good, N, max_diff)
 
 

sasmodels/generate.py

@@ -197 +197 @@
 # TODO: identify model files which have changed since loading and reload them.
 
-__all__ = ["make", "doc", "sources", "use_single"]
+__all__ = ["make", "doc", "sources", "use_single", "use_long_double"]
 
 import sys
@@ -206 +206 @@
 C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c')
 
+F32 = np.dtype('float32')
 F64 = np.dtype('float64')
-F32 = np.dtype('float32')
+try:  # CRUFT: older numpy does not support float128
+    F128 = np.dtype('float128')
+except TypeError:
+    F128 = None
+
 
 # Scale and background, which are parameters common to every form factor
@@ -322 +327 @@
     source = re.sub(r'[^a-zA-Z_](\d*[.]\d+|\d+[.]\d*)([eE][+-]?\d+)?',
                     r'\g<0>f', source)
+    return source
+
+def use_long_double(source):
+    """
+    Convert code from double precision to long double precision.
+    """
+    # Convert double keyword to float.  Accept an 'n' parameter for vector
+    # values, where n is 2, 4, 8 or 16. Assume complex numbers are represented
+    # as cdouble which is typedef'd to double2.
+    source = re.sub(r'(^|[^a-zA-Z0-9_]c?)double(([248]|16)?($|[^a-zA-Z0-9_]))',
+                    r'\1long double\2', source)
+    # Convert floating point constants to single by adding 'f' to the end.
+    # OS/X driver complains if you don't do this.
+    source = re.sub(r'[^a-zA-Z_](\d*[.]\d+|\d+[.]\d*)([eE][+-]?\d+)?',
+                    r'\g<0>L', source)
     return source
 

sasmodels/kernel_template.c

@@ -125 +125 @@
     if (weight > cutoff) {
       const double scattering = Iq(qi, IQ_PARAMETERS);
-      if (scattering >= 0.0) { // scattering cannot be negative
+      //if (scattering >= 0.0) { // scattering cannot be negative
         ret += weight*scattering;
         norm += weight;
@@ -133 +133 @@
         norm_vol += vol_weight;
       #endif
-      }
+      //}
     //else { printf("exclude qx,qy,I:%%g,%%g,%%g\n",qi,scattering); }
     }
@@ -198 +198 @@
 
       const double scattering = Iqxy(qxi, qyi, IQXY_PARAMETERS);
-      if (scattering >= 0.0) { // scattering cannot be negative
+      //if (scattering >= 0.0) { // scattering cannot be negative
         // TODO: use correct angle for spherical correction
         // Definition of theta and phi are probably reversed relative to the
@@ -219 +219 @@
       #endif
         norm_vol += vol_weight;
-      }
+      //}
       //else { printf("exclude qx,qy,I:%%g,%%g,%%g\n",qi,scattering); }
     }

sasmodels/kernelcl.py

@@ -22 +22 @@
 devices, where it can be combined with other structure factors and form
 factors and have instrumental resolution effects applied.
+
+In order to use OpenCL for your models, you will need OpenCL drivers for
+your machine.  These should be available from your graphics card vendor.
+Intel provides OpenCL drivers for CPUs as well as their integrated HD
+graphics chipsets.  AMD also provides drivers for Intel CPUs, but as of
+this writing the performance is lacking compared to the Intel drivers.
+NVidia combines drivers for CUDA and OpenCL in one package.  The result
+is a bit messy if you have multiple drivers installed.  You can see which
+drivers are available by starting python and running:
+
+    import pyopencl as cl
+    cl.create_some_context(interactive=True)
+
+Once you have done that, it will show the available drivers which you
+can select.  It will then tell you that you can use these drivers
+automatically by setting the PYOPENCL_CTX environment variable.
+
+Some graphics cards have multiple devices on the same card.  You cannot
+yet use both of them concurrently to evaluate models, but you can run
+the program twice using a different device for each session.
+
+OpenCL kernels are compiled when needed by the device driver.  Some
+drivers produce compiler output even when there is no error.  You
+can see the output by setting PYOPENCL_COMPILER_OUTPUT=1.  It should be
+harmless, albeit annoying.
 """
 import os

sasmodels/kerneldll.py

@@ -1 @@
-"""
+r"""
 C types wrapper for sasview models.
 
@@ -5 +5 @@
 you wish to allow single precision floating point evaluation for the compiled
 models, otherwise it defaults to *False*.
+
+The compiler command line is stored in the attribute *COMPILE*, with string
+substitutions for %(source)s and %(output)s indicating what to compile and
+where to store it.  The actual command is system dependent.
+
+On windows systems, you have a choice of compilers.  *MinGW* is the GNU
+compiler toolchain, available in packages such as anaconda and PythonXY,
+or available stand alone. This toolchain has had difficulties on some
+systems, and may or may not work for you.  In order to build DLLs, *gcc*
+must be on your path.  If the environment variable *SAS_OPENMP* is given
+then -fopenmp is added to the compiler flags.  This requires a version
+of MinGW compiled with OpenMP support.
+
+An alternative toolchain uses the Microsoft Visual C++ compiler, available
+free from microsoft:
+
+    `http://www.microsoft.com/en-us/download/details.aspx?id=44266`_
+
+Again, this requires that the compiler is available on your path.  This is
+done by running vcvarsall.bat in a windows terminal.  Install locations are
+system dependent, such as:
+
+    C:\Program Files (x86)\Common Files\Microsoft\Visual C++ for Python\9.0\vcvarsall.bat
+
+or maybe
+
+    C:\Users\yourname\AppData\Local\Programs\Common\Microsoft\Visual C++ for Python\9.0\vcvarsall.bat
+
+And again, the environment variable *SAS_OPENMP* controls whether OpenMP is
+used to compile the C code.  This requires the Microsoft vcomp90.dll library,
+which doesn't seem to be included with the compiler, nor does there appear
+to be a public download location.  There may be one on your machine already
+in a location such as:
+
+    C:\Windows\winsxs\x86_microsoft.vc90.openmp*\vcomp90.dll
+
+If you copy this onto your path, such as the python directory or the install
+directory for this application, then OpenMP should be supported.
 """
 
@@ -11 +49 @@
 import tempfile
 import ctypes as ct
-from ctypes import c_void_p, c_int, c_double, c_float
+from ctypes import c_void_p, c_int, c_longdouble, c_double, c_float
 
 import numpy as np
@@ -58 +96 @@
     if np.dtype(dtype) == generate.F32:
         basename += "32"
+    elif np.dtype(dtype) == generate.F64:
+        basename += "64"
+    else:
+        basename += "128"
     return joinpath(DLL_PATH, basename+'.so')
 
@@ -80 +122 @@
     models are allowed as DLLs.
     """
-    if not ALLOW_SINGLE_PRECISION_DLLS: dtype = "double"   # Force 64-bit dll
     dtype = np.dtype(dtype)
+    if dtype == generate.F32 and not ALLOW_SINGLE_PRECISION_DLLS:
+        dtype = generate.F64  # Force 64-bit dll
 
     if callable(info.get('Iq',None)):
@@ -89 +132 @@
         source = generate.use_single(source)
         tempfile_prefix = 'sas_'+info['name']+'32_'
+    elif dtype == generate.F64:
+        tempfile_prefix = 'sas_'+info['name']+'64_'
     else:
-        tempfile_prefix = 'sas_'+info['name']+'_'
+        source = generate.use_long_double(source)
+        tempfile_prefix = 'sas_'+info['name']+'128_'
 
     source_files = generate.sources(info) + [info['filename']]
@@ -158 +204 @@
             raise
 
-        fp = c_float if self.dtype == generate.F32 else c_double
+        fp = (c_float if self.dtype == generate.F32
+              else c_double if self.dtype == generate.F64
+              else c_longdouble)
         pd_args_1d = [c_void_p, fp] + [c_int]*Npd1d if Npd1d else []
         pd_args_2d= [c_void_p, fp] + [c_int]*Npd2d if Npd2d else []
@@ -227 +275 @@
 
     def __call__(self, fixed_pars, pd_pars, cutoff):
-        real = np.float32 if self.q_input.dtype == generate.F32 else np.float64
+        real = (np.float32 if self.q_input.dtype == generate.F32
+                else np.float64 if self.q_input.dtype == generate.F64
+                else np.float128)
 
         nq = c_int(self.q_input.nq)

sasmodels/kernelpy.py

@@ -211 +211 @@
         if w > cutoff:
             I = form(*args)
-            positive = (I >= 0.0)
+            #positive = (I >= 0.0)
 
             # Note: can precompute spherical correction if theta_index is not
@@ -219 +219 @@
                                     if theta_index >= 0 else 1.0)
             #spherical_correction = 1.0
-            ret += w * I * spherical_correction * positive
-            norm += w * positive
+            ret += w * I * spherical_correction #* positive
+            norm += w #* positive
 
             # Volume normalization.
@@ -230 +230 @@
                 vol_args = [args[index] for index in vol_index]
                 vol_weight = np.prod(weight[vol_weight_index])
-                vol += vol_weight * form_volume(*vol_args) * positive
-                vol_norm += vol_weight * positive
+                vol += vol_weight * form_volume(*vol_args) #* positive
+                vol_norm += vol_weight #* positive
 
     positive = (vol * vol_norm != 0.0)

sasmodels/resolution.py

@@ -1074 +1074 @@
 
 if __name__ == "__main__":
-    demo()
-    #main()
+    #demo()
+    main()