Changes in / [edf06e1:ec87470] in sasmodels

Location:

sasmodels

Files:

• compare.py (modified) (2 diffs)
• conversion_table.py (modified) (1 diff)
• convert.py (modified) (1 diff)
• core.py (modified) (6 diffs)
• details.py (modified) (1 diff)
• generate.py (modified) (2 diffs)
• kernel_iq.c (modified) (2 diffs)
• kernel_iq.cl (modified) (12 diffs)
• kernelcl.py (modified) (11 diffs)
• kerneldll.py (modified) (3 diffs)
• models/core_shell_ellipsoid.py (modified) (1 diff)
• models/gaussian_peak.py (modified) (1 diff)
• models/hardsphere.py (modified) (1 diff)
• models/hollow_cylinder.c (modified) (2 diffs)
• models/hollow_cylinder.py (modified) (1 diff)
• models/lamellar_hg_stack_caille.py (modified) (1 diff)
• models/lamellar_stack_caille.py (modified) (1 diff)
• models/linear_pearls.py (modified) (1 diff)
• models/multilayer_vesicle.py (modified) (1 diff)
• models/onion.c (modified) (2 diffs)
• models/onion.py (modified) (1 diff)
• models/polymer_micelle.py (modified) (1 diff)
• models/polymer_micelle_kernel.c (modified) (1 diff)
• models/rpa.c (modified) (1 diff)
• models/surface_fractal.py (modified) (1 diff)
• models/unified_power_Rg.py (modified) (1 diff)
• sasview_model.py (modified) (3 diffs)

sasmodels/compare.py

@@ -451 +451 @@
         # paying for parameter conversion each time to keep life simple, if not fast
         oldpars = revert_pars(model_info, pars)
+        #print("sasview pars",oldpars)
         for k, v in oldpars.items():
             name_attr = k.split('.')  # polydispersity components
@@ -667 +668 @@
     #    h = plt.colorbar()
     #    h.ax.set_title(cbar_title)
+    fig = plt.gcf()
+    fig.suptitle(opts['name'])
 
     if Ncomp > 0 and Nbase > 0 and '-hist' in opts:

sasmodels/conversion_table.py

@@ -509 +509 @@
         "MultiShellModel",
         {
+            "radius": "core_radius",
+            "sld_solvent": "core_sld",
+            "n_pairs": "n_pairs",
+            "thick_shell": "s_thickness",
             "sld": "shell_sld",
             "thick_solvent": "w_thickness",
-            "thick_shell": "s_thickness",
-            "radius": "core_radius",
-            "sld_solvent": "core_sld"
         }
     ],

sasmodels/convert.py

@@ -316 +316 @@
         elif name == 'core_multi_shell':
             pars['n'] = min(math.ceil(pars['n']), 4)
+        elif name == 'onion':
+            pars['n_shells'] = math.ceil(pars['n_shells'])
         elif name == 'spherical_sld':
             for k in range(1, 11):

sasmodels/core.py

@@ -6 +6 @@
 __all__ = [
     "list_models", "load_model", "load_model_info",
-    "build_model", "precompile_dll",
+    "build_model", "precompile_dlls",
     ]
 
 import os
-from os.path import basename, dirname, join as joinpath, splitext
+from os.path import basename, dirname, join as joinpath
 from glob import glob
 
@@ -57 +57 @@
 #    build_model
 
-def list_models():
+KINDS = ("all", "py", "c", "double", "single", "1d", "2d",
+         "nonmagnetic", "magnetic")
+def list_models(kind=None):
     # type: () -> List[str]
     """
     Return the list of available models on the model path.
     """
+    if kind and kind not in KINDS:
+        raise ValueError("kind not in " + ", ".join(KINDS))
     root = dirname(__file__)
     files = sorted(glob(joinpath(root, 'models', "[a-zA-Z]*.py")))
     available_models = [basename(f)[:-3] for f in files]
-    return available_models
+    selected = [name for name in available_models if _matches(name, kind)]
+
+    return selected
+
+def _matches(name, kind):
+    if kind is None or kind == "all":
+        return True
+    info = load_model_info(name)
+    pars = info.parameters.kernel_parameters
+    if kind == "py" and callable(info.Iq):
+        return True
+    elif kind == "c" and not callable(info.Iq):
+        return True
+    elif kind == "double" and not info.single:
+        return True
+    elif kind == "single" and info.single:
+        return True
+    elif kind == "2d" and any(p.type == 'orientation' for p in pars):
+        return True
+    elif kind == "1d" and any(p.type != 'orientation' for p in pars):
+        return True
+    elif kind == "magnetic" and any(p.type == 'sld' for p in pars):
+        return True
+    elif kind == "nonmagnetic" and any(p.type != 'sld' for p in pars):
+        return True
+    return False
 
 def load_model(model_name, dtype=None, platform='ocl'):
@@ -129 +158 @@
             return mixture.MixtureModel(model_info, models)
         elif composition_type == 'product':
-            from . import product
             P, S = models
             return product.ProductModel(model_info, P, S)
@@ -189 +217 @@
     if platform is None:
         platform = "ocl"
-    if platform=="ocl" and not HAVE_OPENCL:
+    if platform == "ocl" and not HAVE_OPENCL:
         platform = "dll"
 
@@ -198 +226 @@
 
     # Convert special type names "half", "fast", and "quad"
-    fast = (dtype=="fast")
+    fast = (dtype == "fast")
     if fast:
         dtype = "single"
-    elif dtype=="quad":
+    elif dtype == "quad":
         dtype = "longdouble"
-    elif dtype=="half":
+    elif dtype == "half":
         dtype = "f16"
 
     # Convert dtype string to numpy dtype.
     if dtype is None:
-        numpy_dtype = generate.F32 if platform=="ocl" and model_info.single else generate.F64
+        numpy_dtype = (generate.F32 if platform == "ocl" and model_info.single
+                       else generate.F64)
     else:
         numpy_dtype = np.dtype(dtype)
 
     # Make sure that the type is supported by opencl, otherwise use dll
-    if platform=="ocl":
+    if platform == "ocl":
         env = kernelcl.environment()
         if not env.has_type(numpy_dtype):
@@ -221 +250 @@
 
     return numpy_dtype, fast, platform
+
+def list_models_main():
+    import sys
+    kind = sys.argv[1] if len(sys.argv) > 1 else "all"
+    print("\n".join(list_models(kind)))
+
+if __name__ == "__main__":
+    list_models_main()

sasmodels/details.py

@@ -176 +176 @@
     if all(len(w)==1 for w in weights):
         call_details = mono_details(kernel.info)
-        data = np.array(scalars+scalars+[1]*len(scalars), dtype=kernel.dtype)
+        # Pad value array to a 32 value boundary
+        data_len = 3*len(scalars)
+        extra = ((data_len+31)//32)*32 - data_len
+        data = np.array(scalars+scalars+[1.]*len(scalars)+[0.]*extra, dtype=kernel.dtype)
     else:
         call_details = poly_details(kernel.info, weights)
-        data = np.hstack(scalars+list(values)+list(weights)).astype(kernel.dtype)
+        # Pad value array to a 32 value boundary
+        data_len = len(scalars) + 2*sum(len(v) for v in values)
+        extra = ((data_len+31)//32)*32 - data_len
+        data = np.hstack(scalars+list(values)+list(weights)+[0.]*extra).astype(kernel.dtype)
     is_magnetic = convert_magnetism(kernel.info.parameters, data)
     #call_details.show()

sasmodels/generate.py

@@ -166 +166 @@
 import re
 import string
-import warnings
 
 import numpy as np  # type: ignore
@@ -497 +496 @@
     if callable(model_info.Iq):
         raise ValueError("can't compile python model")
+        #return None
 
     # TODO: need something other than volume to indicate dispersion parameters

sasmodels/kernel_iq.c

@@ -173 +173 @@
 
 
-  double spherical_correction=1.0;
+#if MAX_PD>0
   const int theta_par = details->theta_par;
-#if MAX_PD>0
   const int fast_theta = (theta_par == p0);
   const int slow_theta = (theta_par >= 0 && !fast_theta);
+  double spherical_correction = 1.0;
 #else
-  const int slow_theta = (theta_par >= 0);
+  // Note: if not polydisperse the weights cancel and we don't need the
+  // spherical correction.
+  const double spherical_correction = 1.0;
 #endif
 
@@ -218 +220 @@
     const double weight1 = 1.0;
 #endif
-    if (slow_theta) { // Theta is not in inner loop
-      spherical_correction = fmax(fabs(cos(M_PI_180*pvec[theta_par])), 1.e-6);
-    }
-#if MAX_PD>0
+#if MAX_PD>0
+  if (slow_theta) { // Theta is not in inner loop
+    spherical_correction = fmax(fabs(cos(M_PI_180*pvec[theta_par])), 1.e-6);
+  }
   while(i0 < n0) {
     pvec[p0] = v0[i0];

sasmodels/kernel_iq.cl

@@ -28 +28 @@
 } ProblemDetails;
 
+// Intel HD 4000 needs private arrays to be a multiple of 4 long
 typedef struct {
     PARAMETER_TABLE
+} ParameterTable;
+typedef union {
+    ParameterTable table;
+    double vector[4*((NUM_PARS+3)/4)];
 } ParameterBlock;
 #endif // _PAR_BLOCK_
@@ -87 +92 @@
   // walk the polydispersity cube.  local_values will be aliased to pvec.
   ParameterBlock local_values;
-  double *pvec = (double *)&local_values;
 
   // Fill in the initial variables
   for (int i=0; i < NUM_PARS; i++) {
-    pvec[i] = values[2+i];
-//if (q_index==0) printf("p%d = %g\n",i, pvec[i]);
+    local_values.vector[i] = values[2+i];
+//if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]);
   }
 
 #if defined(MAGNETIC) && NUM_MAGNETIC>0
-  // Location of the sld parameters in the parameter pvec.
+  // Location of the sld parameters in the parameter vector.
   // These parameters are updated with the effective sld due to magnetism.
   #if NUM_MAGNETIC > 3
@@ -166 +170 @@
 
 
-  double spherical_correction=1.0;
+#if MAX_PD>0
   const int theta_par = details->theta_par;
-#if MAX_PD>0
-  const int fast_theta = (theta_par == p0);
-  const int slow_theta = (theta_par >= 0 && !fast_theta);
+  const bool fast_theta = (theta_par == p0);
+  const bool slow_theta = (theta_par >= 0 && !fast_theta);
+  double spherical_correction = 1.0;
 #else
-  const int slow_theta = (theta_par >= 0);
+  // Note: if not polydisperse the weights cancel and we don't need the
+  // spherical correction.
+  const double spherical_correction = 1.0;
 #endif
 
@@ -181 +187 @@
   const double weight5 = 1.0;
   while (i4 < n4) {
-    pvec[p4] = v4[i4];
+    local_values.vector[p4] = v4[i4];
     double weight4 = w4[i4] * weight5;
-//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 4, p4, i4, n4, pvec[p4], weight4);
+//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 4, p4, i4, n4, local_values.vector[p4], weight4);
 #elif MAX_PD>3
     const double weight4 = 1.0;
@@ -189 +195 @@
 #if MAX_PD>3
   while (i3 < n3) {
-    pvec[p3] = v3[i3];
+    local_values.vector[p3] = v3[i3];
     double weight3 = w3[i3] * weight4;
-//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 3, p3, i3, n3, pvec[p3], weight3);
+//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 3, p3, i3, n3, local_values.vector[p3], weight3);
 #elif MAX_PD>2
     const double weight3 = 1.0;
@@ -197 +203 @@
 #if MAX_PD>2
   while (i2 < n2) {
-    pvec[p2] = v2[i2];
+    local_values.vector[p2] = v2[i2];
     double weight2 = w2[i2] * weight3;
-//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 2, p2, i2, n2, pvec[p2], weight2);
+//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 2, p2, i2, n2, local_values.vector[p2], weight2);
 #elif MAX_PD>1
     const double weight2 = 1.0;
@@ -205 +211 @@
 #if MAX_PD>1
   while (i1 < n1) {
-    pvec[p1] = v1[i1];
+    local_values.vector[p1] = v1[i1];
     double weight1 = w1[i1] * weight2;
-//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 1, p1, i1, n1, pvec[p1], weight1);
+//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 1, p1, i1, n1, local_values.vector[p1], weight1);
 #elif MAX_PD>0
     const double weight1 = 1.0;
 #endif
-    if (slow_theta) { // Theta is not in inner loop
-      spherical_correction = fmax(fabs(cos(M_PI_180*pvec[theta_par])), 1.e-6);
-    }
-#if MAX_PD>0
+#if MAX_PD>0
+  if (slow_theta) { // Theta is not in inner loop
+    spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6);
+  }
   while(i0 < n0) {
-    pvec[p0] = v0[i0];
+    local_values.vector[p0] = v0[i0];
     double weight0 = w0[i0] * weight1;
-//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, pvec[p0], weight0);
+//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0);
     if (fast_theta) { // Theta is in inner loop
-      spherical_correction = fmax(fabs(cos(M_PI_180*pvec[p0])), 1.e-6);
+      spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6);
     }
 #else
@@ -226 +232 @@
 #endif
 
-//if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, pvec[i]); printf("\n"); }
+//if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, local_values.vector[i]); printf("\n"); }
 //if (q_index == 0) printf("sphcor: %g\n", spherical_correction);
 
     #ifdef INVALID
-    if (!INVALID(local_values))
+    if (!INVALID(local_values.table))
     #endif
     {
@@ -239 +245 @@
         // would be problems looking at models with theta=90.
         const double weight = weight0 * spherical_correction;
-        pd_norm += weight * CALL_VOLUME(local_values);
+        pd_norm += weight * CALL_VOLUME(local_values.table);
 
 #if defined(MAGNETIC) && NUM_MAGNETIC > 0
@@ -265 +271 @@
                 #define M3 NUM_PARS+13
                 #define SLD(_M_offset, _sld_offset) \
-                    pvec[_sld_offset] = xs * (axis \
+                    local_values.vector[_sld_offset] = xs * (axis \
                     ? (index==1 ? -values[_M_offset+2] : values[_M_offset+2]) \
                     : mag_sld(qx, qy, pk, values[_M_offset], values[_M_offset+1], \
@@ -283 +289 @@
                 }
                 #endif
-                scattering += CALL_IQ(q, q_index, local_values);
+                scattering += CALL_IQ(q, q_index, local_values.table);
               }
             }
@@ -289 +295 @@
         }
 #else  // !MAGNETIC
-        const double scattering = CALL_IQ(q, q_index, local_values);
+        const double scattering = CALL_IQ(q, q_index, local_values.table);
 #endif // !MAGNETIC
         this_result += weight * scattering;

sasmodels/kernelcl.py

@@ -65 +65 @@
     raise RuntimeError("OpenCL not available")
 
-from pyopencl import mem_flags as mf  # type: ignore
-from pyopencl.characterize import get_fast_inaccurate_build_options  # type: ignore
-# CRUFT: pyopencl < 2017.1  (as of June 2016 needs quotes around include path)
-def _quote_path(v):
-    """
-    Quote the path if it is not already quoted.
-
-    If v starts with '-', then assume that it is a -I option or similar
-    and do not quote it.  This is fragile:  -Ipath with space needs to
-    be quoted.
-    """
-    return '"'+v+'"' if v and ' ' in v and not v[0] in "\"'-" else v
-
-if hasattr(cl, '_DEFAULT_INCLUDE_OPTIONS'):
-    cl._DEFAULT_INCLUDE_OPTIONS = [_quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS]
-
 from pyopencl import mem_flags as mf
 from pyopencl.characterize import get_fast_inaccurate_build_options
@@ -93 +77 @@
 except ImportError:
     pass
+
+# CRUFT: pyopencl < 2017.1  (as of June 2016 needs quotes around include path)
+def quote_path(v):
+    """
+    Quote the path if it is not already quoted.
+
+    If v starts with '-', then assume that it is a -I option or similar
+    and do not quote it.  This is fragile:  -Ipath with space needs to
+    be quoted.
+    """
+    return '"'+v+'"' if v and ' ' in v and not v[0] in "\"'-" else v
+
+def fix_pyopencl_include():
+    import pyopencl as cl
+    if hasattr(cl, '_DEFAULT_INCLUDE_OPTIONS'):
+        cl._DEFAULT_INCLUDE_OPTIONS = [quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS]
+
+fix_pyopencl_include()
+
 
 # The max loops number is limited by the amount of local memory available
@@ -256 +259 @@
         Return a OpenCL context for the kernels of type dtype.
         """
-        for context, queue in zip(self.context, self.queues):
+        for context in self.context:
             if all(has_type(d, dtype) for d in context.devices):
                 return context
@@ -282 +285 @@
         if key not in self.compiled:
             context = self.get_context(dtype)
-            logging.info("building %s for OpenCL %s"
-                         % (key, context.devices[0].name.strip()))
+            logging.info("building %s for OpenCL %s", key,
+                         context.devices[0].name.strip())
             program = compile_model(self.get_context(dtype),
                                     str(source), dtype, fast)
@@ -299 +302 @@
 
 def _get_default_context():
-    # type: () -> cl.Context
+    # type: () -> List[cl.Context]
     """
     Get an OpenCL context, preferring GPU over CPU, and preferring Intel
@@ -318 +321 @@
     for platform in cl.get_platforms():
         # AMD provides a much weaker CPU driver than Intel/Apple, so avoid it.
-        # If someone has bothered to install the AMD/NVIDIA drivers, prefer them over the integrated
-        # graphics driver that may have been supplied with the CPU chipset.
-        preferred_cpu = platform.vendor.startswith('Intel') or platform.vendor.startswith('Apple')
-        preferred_gpu = platform.vendor.startswith('Advanced') or platform.vendor.startswith('NVIDIA')
+        # If someone has bothered to install the AMD/NVIDIA drivers, prefer
+        # them over the integrated graphics driver that may have been supplied
+        # with the CPU chipset.
+        preferred_cpu = (platform.vendor.startswith('Intel')
+                         or platform.vendor.startswith('Apple'))
+        preferred_gpu = (platform.vendor.startswith('Advanced')
+                         or platform.vendor.startswith('NVIDIA'))
         for device in platform.get_devices():
             if device.type == cl.device_type.GPU:
-                # If the existing type is not GPU then it will be CUSTOM or ACCELERATOR,
-                # so don't override it.
+                # If the existing type is not GPU then it will be CUSTOM
+                # or ACCELERATOR so don't override it.
                 if gpu is None or (preferred_gpu and gpu.type == cl.device_type.GPU):
                     gpu = device
@@ -334 +340 @@
                 # System has cl.device_type.ACCELERATOR or cl.device_type.CUSTOM
                 # Intel Phi for example registers as an accelerator
-                # Since the user installed a custom device on their system and went through the
-                # pain of sorting out OpenCL drivers for it, lets assume they really do want to
-                # use it as their primary compute device.
+                # Since the user installed a custom device on their system
+                # and went through the pain of sorting out OpenCL drivers for
+                # it, lets assume they really do want to use it as their
+                # primary compute device.
                 gpu = device
 
-    # order the devices by gpu then by cpu; when searching for an available device by data type they
-    # will be checked in this order, which means that if the gpu supports double then the cpu will never
-    # be used (though we may make it possible to explicitly request the cpu at some point).
+    # order the devices by gpu then by cpu; when searching for an available
+    # device by data type they will be checked in this order, which means
+    # that if the gpu supports double then the cpu will never be used (though
+    # we may make it possible to explicitly request the cpu at some point).
     devices = []
     if gpu is not None:
@@ -372 +380 @@
         self.fast = fast
         self.program = None # delay program creation
+        self._kernels = None
 
     def __getstate__(self):
@@ -394 +403 @@
             names = [generate.kernel_name(self.info, k) for k in variants]
             kernels = [getattr(self.program, k) for k in names]
-            self._kernels = dict((k,v) for k,v in zip(variants, kernels))
+            self._kernels = dict((k, v) for k, v in zip(variants, kernels))
         is_2d = len(q_vectors) == 2
         if is_2d:
@@ -500 +509 @@
     def __init__(self, kernel, dtype, model_info, q_vectors):
         # type: (cl.Kernel, np.dtype, ModelInfo, List[np.ndarray]) -> None
-        max_pd = model_info.parameters.max_pd
-        npars = len(model_info.parameters.kernel_parameters)-2
         q_input = GpuInput(q_vectors, dtype)
         self.kernel = kernel
@@ -516 +523 @@
 
         self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
-                               q_input.global_size[0] * dtype.itemsize)
+                                  q_input.global_size[0] * dtype.itemsize)
         self.q_input = q_input # allocated by GpuInput above
 
-        self._need_release = [ self.result_b, self.q_input ]
+        self._need_release = [self.result_b, self.q_input]
         self.real = (np.float32 if dtype == generate.F32
                      else np.float64 if dtype == generate.F64

sasmodels/kerneldll.py

@@ -57 +57 @@
 import numpy as np  # type: ignore
 
+try:
+    import tinycc
+except ImportError:
+    tinycc = None
+
 from . import generate
 from .kernel import KernelModel, Kernel
@@ -70 +75 @@
     pass
 
-if os.name == 'nt':
-    ARCH = "" if sys.maxint > 2**32 else "x86"  # maxint=2**31-1 on 32 bit
-    # Windows compiler; check if TinyCC is available
-    try:
-        import tinycc
-    except ImportError:
-        tinycc = None
-    # call vcvarsall.bat before compiling to set path, headers, libs, etc.
+if "SAS_COMPILER" in os.environ:
+    compiler = os.environ["SAS_COMPILER"]
+elif os.name == 'nt':
+    # If vcvarsall.bat has been called, then VCINSTALLDIR is in the environment
+    # and we can use the MSVC compiler.  Otherwise, if tinycc is available
+    # the use it.  Otherwise, hope that mingw is available.
     if "VCINSTALLDIR" in os.environ:
-        # MSVC compiler is available, so use it.  OpenMP requires a copy of
-        # vcomp90.dll on the path.  One may be found here:
-        #       C:/Windows/winsxs/x86_microsoft.vc90.openmp*/vcomp90.dll
-        # Copy this to the python directory and uncomment the OpenMP COMPILE
-        # TODO: remove intermediate OBJ file created in the directory
-        # TODO: maybe don't use randomized name for the c file
-        # TODO: maybe ask distutils to find MSVC
-        CC = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG".split()
-        if "SAS_OPENMP" in os.environ:
-            CC.append("/openmp")
-        LN = "/link /DLL /INCREMENTAL:NO /MANIFEST".split()
-        def compile_command(source, output):
-            return CC + ["/Tp%s"%source] + LN + ["/OUT:%s"%output]
+        compiler = "msvc"
     elif tinycc:
-        # TinyCC compiler.
-        CC = [tinycc.TCC] + "-shared -rdynamic -Wall".split()
-        def compile_command(source, output):
-            return CC + [source, "-o", output]
+        compiler = "tinycc"
     else:
-        # MinGW compiler.
-        CC = "gcc -shared -std=c99 -O2 -Wall".split()
-        if "SAS_OPENMP" in os.environ:
-            CC.append("-fopenmp")
-        def compile_command(source, output):
-            return CC + [source, "-o", output, "-lm"]
+        compiler = "mingw"
 else:
-    ARCH = ""
+    compiler = "unix"
+
+ARCH = "" if sys.maxint > 2**32 else "x86"  # maxint=2**31-1 on 32 bit
+if compiler == "unix":
     # Generic unix compile
     # On mac users will need the X code command line tools installed
@@ -112 +98 @@
     # add openmp support if not running on a mac
     if sys.platform != "darwin":
+        CC.append("-fopenmp")
+    def compile_command(source, output):
+        return CC + [source, "-o", output, "-lm"]
+elif compiler == "msvc":
+    # Call vcvarsall.bat before compiling to set path, headers, libs, etc.
+    # MSVC compiler is available, so use it.  OpenMP requires a copy of
+    # vcomp90.dll on the path.  One may be found here:
+    #       C:/Windows/winsxs/x86_microsoft.vc90.openmp*/vcomp90.dll
+    # Copy this to the python directory and uncomment the OpenMP COMPILE
+    # TODO: remove intermediate OBJ file created in the directory
+    # TODO: maybe don't use randomized name for the c file
+    # TODO: maybe ask distutils to find MSVC
+    CC = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG".split()
+    if "SAS_OPENMP" in os.environ:
+        CC.append("/openmp")
+    LN = "/link /DLL /INCREMENTAL:NO /MANIFEST".split()
+    def compile_command(source, output):
+        return CC + ["/Tp%s"%source] + LN + ["/OUT:%s"%output]
+elif compiler == "tinycc":
+    # TinyCC compiler.
+    CC = [tinycc.TCC] + "-shared -rdynamic -Wall".split()
+    def compile_command(source, output):
+        return CC + [source, "-o", output]
+elif compiler == "mingw":
+    # MinGW compiler.
+    CC = "gcc -shared -std=c99 -O2 -Wall".split()
+    if "SAS_OPENMP" in os.environ:
         CC.append("-fopenmp")
     def compile_command(source, output):

sasmodels/models/core_shell_ellipsoid.py

@@ -99 +99 @@
 category = "shape:ellipsoid"
 
-single = False  # TODO: maybe using sph_j1c inside gfn would help?
 # pylint: disable=bad-whitespace, line-too-long
 #             ["name", "units", default, [lower, upper], "type", "description"],

sasmodels/models/gaussian_peak.py

@@ -38 +38 @@
 category = "shape-independent"
 
-single = False
 #             ["name", "units", default, [lower, upper], "type","description"],
 parameters = [["q0", "1/Ang", 0.05, [-inf, inf], "", "Peak position"],

sasmodels/models/hardsphere.py

@@ -58 +58 @@
 category = "structure-factor"
 structure_factor = True
-single = False
+single = False # TODO: check
 
 #             ["name", "units", default, [lower, upper], "type","description"],

sasmodels/models/hollow_cylinder.c

@@ -6 +6 @@
         double solvent_sld, double theta, double phi);
 
-#define INVALID(v) (v.core_radius >= v.radius || v.radius >= v.length)
+#define INVALID(v) (v.core_radius >= v.radius)
 
 // From Igor library
-static double hollow_cylinder_scaling(double integrand, double delrho, double volume)
+static double hollow_cylinder_scaling(
+    double integrand, double delrho, double volume)
 {
-        double answer;
-        // Multiply by contrast^2
-        answer = integrand*delrho*delrho;
+    double answer;
+    // Multiply by contrast^2
+    answer = integrand*delrho*delrho;
 
-        //normalize by cylinder volume
-        answer *= volume*volume;
+    //normalize by cylinder volume
+    answer *= volume*volume;
 
-        //convert to [cm-1]
-        answer *= 1.0e-4;
+    //convert to [cm-1]
+    answer *= 1.0e-4;
 
-        return answer;
+    return answer;
 }
 
-static double _hollow_cylinder_kernel(double q, double core_radius, double radius,
-        double length, double dum)
+
+static double _hollow_cylinder_kernel(
+    double q, double core_radius, double radius, double length, double dum)
 {
-    double gamma,arg1,arg2,lam1,lam2,psi,sinarg,t2,retval;              //local variables
-    
-    gamma = core_radius/radius;
-    arg1 = q*radius*sqrt(1.0-dum*dum);          //1=shell (outer radius)
-    arg2 = q*core_radius*sqrt(1.0-dum*dum);                     //2=core (inner radius)
-    if (arg1 == 0.0){
-        lam1 = 1.0;
-    }else{
-        lam1 = sas_J1c(arg1);
-    }
-    if (arg2 == 0.0){
-        lam2 = 1.0;
-    }else{
-        lam2 = sas_J1c(arg2);
-    }
-    //Todo: Need to check psi behavior as gamma goes to 1.
-    psi = (lam1 -  gamma*gamma*lam2)/(1.0-gamma*gamma);         //SRK 10/19/00
-    sinarg = q*length*dum/2.0;
-    if (sinarg == 0.0){
-        t2 = 1.0;
-    }else{
-        t2 = sin(sinarg)/sinarg;
-    }
+    const double qs = q*sqrt(1.0-dum*dum);
+    const double lam1 = sas_J1c(radius*qs);
+    const double lam2 = sas_J1c(core_radius*qs);
+    const double gamma_sq = square(core_radius/radius);
+    //Note: lim_{r -> r_c} psi = J0(core_radius*qs)
+    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00
+    const double t2 = sinc(q*length*dum/2.0);
+    return square(psi*t2);
+}
 
-    retval = psi*psi*t2*t2;
-    
-    return(retval);
-}
-static double hollow_cylinder_analytical_2D_scaled(double q, double q_x, double q_y, double radius, double core_radius, double length, double sld,
-        double solvent_sld, double theta, double phi) {
-        double cyl_x, cyl_y; //, cyl_z
-        //double q_z;
-        double vol, cos_val, delrho;
-        double answer;
-        //convert angle degree to radian
-        double pi = 4.0*atan(1.0);
-        theta = theta * pi/180.0;
-        phi = phi * pi/180.0;
-        delrho = solvent_sld - sld;
 
-        // Cylinder orientation
-        cyl_x = cos(theta) * cos(phi);
-        cyl_y = sin(theta);
-        //cyl_z = -cos(theta) * sin(phi);
+static double hollow_cylinder_analytical_2D_scaled(
+    double q, double q_x, double q_y, double radius, double core_radius,
+    double length, double sld, double solvent_sld, double theta, double phi)
+{
+    double cyl_x, cyl_y; //, cyl_z
+    //double q_z;
+    double vol, cos_val, delrho;
+    double answer;
+    //convert angle degree to radian
+    theta = theta * M_PI_180;
+    phi = phi * M_PI_180;
+    delrho = solvent_sld - sld;
 
-        // q vector
-        //q_z = 0;
+    // Cylinder orientation
+    cyl_x = cos(theta) * cos(phi);
+    cyl_y = sin(theta);
+    //cyl_z = -cos(theta) * sin(phi);
 
-        // Compute the angle btw vector q and the
-        // axis of the cylinder
-        cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
+    // q vector
+    //q_z = 0;
 
-        answer = _hollow_cylinder_kernel(q, core_radius, radius, length, cos_val);
+    // Compute the angle btw vector q and the
+    // axis of the cylinder
+    cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
 
-        vol = form_volume(radius, core_radius, length);
-        answer = hollow_cylinder_scaling(answer, delrho, vol);
+    answer = _hollow_cylinder_kernel(q, core_radius, radius, length, cos_val);
 
-        return answer;
+    vol = form_volume(radius, core_radius, length);
+    answer = hollow_cylinder_scaling(answer, delrho, vol);
+
+    return answer;
 }
 
@@ -90 +76 @@
 double form_volume(double radius, double core_radius, double length)
 {
-        double pi = 4.0*atan(1.0);
-        double v_shell = pi*length*(radius*radius-core_radius*core_radius);
-        return(v_shell);
+    double v_shell = M_PI*length*(radius*radius-core_radius*core_radius);
+    return(v_shell);
 }
 
 
-double Iq(double q, double radius, double core_radius, double length, double sld,
-        double solvent_sld)
+double Iq(double q, double radius, double core_radius, double length,
+    double sld, double solvent_sld)
 {
     int i;
-        int nord=76;                    //order of integration
-        double lower,upper,zi, inter;           //upper and lower integration limits
-        double summ,answer,delrho;                      //running tally of integration
-        double norm,volume;     //final calculation variables
-        
-        delrho = solvent_sld - sld;
-        lower = 0.0;
-        upper = 1.0;            //limits of numerical integral
+    double lower,upper,zi, inter;               //upper and lower integration limits
+    double summ,answer,delrho;                  //running tally of integration
+    double norm,volume; //final calculation variables
 
-        summ = 0.0;                     //initialize intergral
-        for(i=0;i<nord;i++) {
-                zi = ( Gauss76Z[i] * (upper-lower) + lower + upper )/2.0;
-                inter = Gauss76Wt[i] * _hollow_cylinder_kernel(q, core_radius, radius, length, zi);
-                summ += inter;
-        }
-        
-        norm = summ*(upper-lower)/2.0;
-        volume = form_volume(radius, core_radius, length);
-        answer = hollow_cylinder_scaling(norm, delrho, volume);
-        
-        return(answer);
+    lower = 0.0;
+    upper = 1.0;                //limits of numerical integral
+
+    summ = 0.0;                 //initialize intergral
+    for (i=0;i<76;i++) {
+        zi = ( Gauss76Z[i] * (upper-lower) + lower + upper )/2.0;
+        inter = Gauss76Wt[i] * _hollow_cylinder_kernel(q, core_radius, radius, length, zi);
+        summ += inter;
+    }
+
+    norm = summ*(upper-lower)/2.0;
+    volume = form_volume(radius, core_radius, length);
+    delrho = solvent_sld - sld;
+    answer = hollow_cylinder_scaling(norm, delrho, volume);
+
+    return(answer);
 }
 
-double Iqxy(double qx, double qy, double radius, double core_radius, double length, double sld,
-        double solvent_sld, double theta, double phi)
+
+double Iqxy(double qx, double qy, double radius, double core_radius,
+    double length, double sld, double solvent_sld, double theta, double phi)
 {
-        double q;
-        q = sqrt(qx*qx+qy*qy);
-        return hollow_cylinder_analytical_2D_scaled(q, qx/q, qy/q, radius, core_radius, length, sld, solvent_sld, theta, phi);
+    const double q = sqrt(qx*qx+qy*qy);
+    return hollow_cylinder_analytical_2D_scaled(q, qx/q, qy/q, radius, core_radius, length, sld, solvent_sld, theta, phi);
 }

sasmodels/models/hollow_cylinder.py

@@ -79 +79 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/polevl.c","lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"]
+source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"]
 
 # pylint: disable=W0613

sasmodels/models/lamellar_hg_stack_caille.py

@@ -90 +90 @@
 category = "shape:lamellae"
 
-single = False
+single = False  # TODO: check
 parameters = [
     #   [ "name", "units", default, [lower, upper], "type",

sasmodels/models/lamellar_stack_caille.py

@@ -83 +83 @@
 category = "shape:lamellae"
 
-single = False
+single = False  # TODO: check
 # pylint: disable=bad-whitespace, line-too-long
 #             ["name", "units", default, [lower, upper], "type","description"],

sasmodels/models/linear_pearls.py

@@ -61 +61 @@
     ]
 # pylint: enable=bad-whitespace, line-too-long
+single=False
 
 source = ["linear_pearls.c"]

sasmodels/models/multilayer_vesicle.py

@@ -74 +74 @@
     ["volfraction", "",  0.05, [0.0, 1],  "", "volume fraction of vesicles"],
     ["radius", "Ang", 60.0, [0.0, inf],  "", "Core radius of the multishell"],
-    ["thick_shell", "Ang",        10.0, [0.0, inf],  "sld", "Shell thickness"],
+    ["thick_shell", "Ang",        10.0, [0.0, inf],  "", "Shell thickness"],
     ["thick_solvent", "Ang",        10.0, [0.0, inf],  "", "Water thickness"],
     ["sld_solvent",    "1e-6/Ang^2",  6.4, [-inf, inf], "sld", "Core scattering length density"],

sasmodels/models/onion.c

@@ -2 +2 @@
 static double
 f_exp(double q, double r, double sld_in, double sld_out,
-    double thickness, double A)
-{
-  const double vol = M_4PI_3 * cube(r);
-  const double qr = q * r;
-  const double alpha = A * r/thickness;
-  const double bes = sph_j1c(qr);
-  const double B = (sld_out - sld_in)/expm1(A);
-  const double C = sld_in - B;
-  double fun;
-  if (qr == 0.0) {
-    fun = 1.0;
-  } else if (fabs(A) > 0.0) {
-    const double qrsq = qr*qr;
-    const double alphasq = alpha*alpha;
-    const double sumsq = alphasq + qrsq;
-    double sinqr, cosqr;
-    SINCOS(qr, sinqr, cosqr);
-    fun = -3.0*(
-            ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq
-                - (alpha*sinqr/qr - cosqr)
-        ) / sumsq;
-  } else {
-    fun = bes;
-  }
-  return vol * (B*fun + C*bes);
-}
-
-static double
-f_linear(double q, double r, double sld, double slope)
+    double thickness, double A, double side)
 {
   const double vol = M_4PI_3 * cube(r);
   const double qr = q * r;
   const double bes = sph_j1c(qr);
-  double fun = 0.0;
-  if (qr > 0.0) {
-    const double qrsq = qr*qr;
+  const double alpha = A * r/thickness;
+  double result;
+  if (qr == 0.0) {
+    result = 1.0;
+  } else if (fabs(A) > 0.0) {
+    const double qrsq = qr * qr;
+    const double alphasq = alpha * alpha;
+    const double sumsq = alphasq + qrsq;
     double sinqr, cosqr;
     SINCOS(qr, sinqr, cosqr);
-    // Jae-He's code seems to simplify to this
-    //     fun = 3.0 * slope * r * (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq);
-    // Rederiving the math, we get the following instead:
-    fun = 3.0 * slope * r * (2.0*cosqr + qr*sinqr)/(qrsq*qrsq);
+    const double t1 = (alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr;
+    const double t2 = alpha*sinqr/qr - cosqr;
+    const double fun = -3.0*(t1/sumsq - t2)/sumsq;
+    const double slope = (sld_out - sld_in)/expm1(A);
+    const double contrast = slope*exp(A*side);
+    result = contrast*fun + (sld_in-slope)*bes;
+  } else {
+    result = sld_in*bes;
   }
-  return vol * (sld*bes + fun);
-}
-
-static double
-f_constant(double q, double r, double sld)
-{
-  const double bes = sph_j1c(q * r);
-  const double vol = M_4PI_3 * cube(r);
-  return sld * vol * bes;
+  return vol * result;
 }
 
@@ -69 +42 @@
 static double
 Iq(double q, double sld_core, double core_radius, double sld_solvent,
-    double n, double sld_in[], double sld_out[], double thickness[],
+    double n_shells, double sld_in[], double sld_out[], double thickness[],
     double A[])
 {
-  int i;
-  double r = core_radius;
-  double f = f_constant(q, r, sld_core);
-  for (i=0; i<n; i++){
-    const double r0 = r;
-    r += thickness[i];
-    if (r == r0) {
-      // no thickness, so nothing to add
-    } else if (fabs(A[i]) < 1.0e-16 || sld_out[i] == sld_in[i]) {
-      f -= f_constant(q, r0, sld_in[i]);
-      f += f_constant(q, r, sld_in[i]);
-    } else if (fabs(A[i]) < 1.0e-4) {
-      const double slope = (sld_out[i] - sld_in[i])/thickness[i];
-      f -= f_linear(q, r0, sld_in[i], slope);
-      f += f_linear(q, r, sld_out[i], slope);
-    } else {
-      f -= f_exp(q, r0, sld_in[i], sld_out[i], thickness[i], A[i]);
-      f += f_exp(q, r, sld_in[i], sld_out[i], thickness[i], A[i]);
-    }
+  int n = (int)(n_shells+0.5);
+  double r_out = core_radius;
+  double f = f_exp(q, r_out, sld_core, 0.0, 0.0, 0.0, 0.0);
+  for (int i=0; i < n; i++){
+    const double r_in = r_out;
+    r_out += thickness[i];
+    f -= f_exp(q, r_in, sld_in[i], sld_out[i], thickness[i], A[i], 0.0);
+    f += f_exp(q, r_out, sld_in[i], sld_out[i], thickness[i], A[i], 1.0);
   }
-  f -= f_constant(q, r, sld_solvent);
+  f -= f_exp(q, r_out, sld_solvent, 0.0, 0.0, 0.0, 0.0);
   const double f2 = f * f * 1.0e-4;
 

sasmodels/models/onion.py

@@ -394 +394 @@
     # Could also specify them individually as
     # "A1": 0, "A2": -1, "A3": 1e-4, "A4": 1,
+    #"core_radius_pd_n": 10,
+    #"core_radius_pd": 0.4,
+    #"thickness4_pd_n": 10,
+    #"thickness4_pd": 0.4,
     }
 

sasmodels/models/polymer_micelle.py

@@ -51 +51 @@
 # pylint: enable=bad-whitespace, line-too-long
 
+single = False
+
 source = ["lib/sph_j1c.c", "polymer_micelle_kernel.c"]
 

sasmodels/models/polymer_micelle_kernel.c

@@ -36 +36 @@
     // Self-correlation term of the chains
     const double qrg2 = q*radius_gyr*q*radius_gyr;
-    const double debye_chain = (qrg2 == 0.0) ? 1.0 : 2.0*(exp(-qrg2)-1+qrg2)/(qrg2*qrg2);
+    const double debye_chain = (qrg2 == 0.0) ? 1.0 : 2.0*(expm1(-qrg2)+qrg2)/(qrg2*qrg2);
     const double term2 = n_aggreg * beta_corona * beta_corona * debye_chain;
 
     // Interference cross-term between core and chains
-    const double chain_ampl = (qrg2 == 0.0) ? 1.0 : (1-exp(-qrg2))/qrg2;
+    const double chain_ampl = (qrg2 == 0.0) ? 1.0 : -expm1(-qrg2)/qrg2;
     const double bes_corona = sinc(q*(radius_core + d_penetration * radius_gyr));
     const double term3 = 2 * n_aggreg * n_aggreg * beta_core * beta_corona *

sasmodels/models/rpa.c

@@ -25 +25 @@
   double S0ba,Pbb,S0bb,Pbc,S0bc,Pbd,S0bd;
   double S0ca,S0cb,Pcc,S0cc,Pcd,S0cd;
-  double S0da,S0db,S0dc,Pdd,S0dd;
-  double Kaa,Kbb,Kcc,Kdd;
-  double Kba,Kca,Kda,Kcb,Kdb,Kdc;
+  double S0da,S0db,S0dc;
+  double Pdd,S0dd;
+  double Kaa,Kbb,Kcc;
+  double Kba,Kca,Kcb;
+  double Kda,Kdb,Kdc,Kdd;
   double Zaa,Zab,Zac,Zba,Zbb,Zbc,Zca,Zcb,Zcc;
   double DenT,T11,T12,T13,T21,T22,T23,T31,T32,T33;

sasmodels/models/surface_fractal.py

@@ -80 +80 @@
 parameters = [["radius",        "Ang", 10.0, [0, inf],   "",
                "Particle radius"],
-              ["surface_dim",   "",    2.0,  [0, inf],   "",
+              ["surface_dim",   "",    2.0,  [1, 3],   "",
                "Surface fractal dimension"],
               ["cutoff_length", "Ang", 500., [0.0, inf], "",

sasmodels/models/unified_power_Rg.py

@@ -80 +80 @@
 
     with errstate(divide='ignore', invalid='ignore'):
-        result = np.empty(q.shape, 'd')
+        result = np.zeros(q.shape, 'd')
         for i in range(ilevel):
             exp_now = exp(-(q*rg[i])**2/3.)

sasmodels/sasview_model.py

@@ -514 +514 @@
         pairs = [self._get_weights(p) for p in parameters.call_parameters]
         call_details, values, is_magnetic = build_details(calculator, pairs)
+        #call_details.show()
+        #print("pairs", pairs)
+        #print("params", self.params)
+        #print("values", values)
+        #print("is_mag", is_magnetic)
         result = calculator(call_details, values, cutoff=self.cutoff,
                             magnetic=is_magnetic)
@@ -598 +603 @@
         if par.name not in self.params:
             if par.name == self.multiplicity_info.control:
-                return [self.multiplicity], []
+                return [self.multiplicity], [1.0]
             else:
-                return [np.NaN], []
+                return [np.NaN], [1.0]
         elif par.polydisperse:
             dis = self.dispersion[par.name]
@@ -608 +613 @@
             return value, weight / np.sum(weight)
         else:
-            return [self.params[par.name]], []
+            return [self.params[par.name]], [1.0]
 
 def test_model():