Changeset c0c3393 in sasmodels

Timestamp:

Mar 21, 2016 4:57:19 AM (9 years ago)

Author:

wojciech

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

9f51b0e

Parents:

d148fbc (diff), 7ff3cf3 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'polydisp' of ​https://github.com/SasView/sasmodels into polydisp

Location:

sasmodels

Files:

• core.py (modified) (2 diffs)
• generate.py (modified) (4 diffs)
• kernel_iq.c (modified) (2 diffs)
• kerneldll.py (modified) (2 diffs)
• kernelpy.py (modified) (3 diffs)
• models/src (added)
• src (added)

sasmodels/core.py

@@ -226 +226 @@
         weights = np.array([1.0])
         details = kernel.info['mono_details']
-        return kernel(pars, weights, details, cutoff)
+        return kernel(details, pars, weights, cutoff)
     else:
         pairs = [get_weights(p, values) for p in kernel.info['parameters']]
@@ -232 +232 @@
         details = generate.poly_details(kernel.info, weights, pars)
         weights, pars = [np.hstack(v) for v in (weights, pars)]
-        return kernel(pars, weights, details, cutoff)
+        return kernel(details, pars, weights, cutoff)
 
 def call_ER_VR(model_info, vol_pars):

sasmodels/generate.py

@@ -650 +650 @@
     model_info['has_2d'] = par_type['orientation'] or par_type['magnetic']
 
-def mono_details(max_pd, npars):
+def mono_details(model_info):
+    max_pd = model_info['max_pd']
+    npars = len(model_info['parameters']) - 2
     p = 5*max_pd
     c = p + 3*npars
 
-    mono = np.zeros(c + 3, 'i')
-    mono[0*max_pd:1*max_pd] = range(max_pd)       # pd_par: arbitrary order; use first
-    mono[1*max_pd:2*max_pd] = [1]*max_pd          # pd_length: only one element
-    mono[2*max_pd:3*max_pd] = range(2, max_pd+2)  # pd_offset: skip scale and background
-    mono[3*max_pd:4*max_pd] = [1]*max_pd          # pd_stride: vectors of length 1
-    mono[4*max_pd:5*max_pd] = [0]*max_pd          # pd_isvol: doens't matter if no norm
-    mono[p+0*npars:p+1*npars] = range(2, npars+2) # par_offset
-    mono[p+1*npars:p+2*npars] = [0]*npars         # no coordination
-    #mono[p+npars] = 1 # par_coord[0] is coordinated with the first par?
-    mono[p+2*npars:p+3*npars] = 0 # fast coord with 0
-    mono[c]   = 1     # fast_coord_count: one fast index
-    mono[c+1] = -1  # theta_var: None
-    mono[c+2] = 0   # fast_theta: False
-    return mono
-
-def poly_details(model_info, weights, pars, constraints=None):
+    details = np.zeros(c + 3, 'int32')
+    details[0*max_pd:1*max_pd] = range(max_pd)       # pd_par: arbitrary order; use first
+    details[1*max_pd:2*max_pd] = [1]*max_pd          # pd_length: only one element
+    details[2*max_pd:3*max_pd] = range(2, max_pd+2)  # pd_offset: skip scale and background
+    details[3*max_pd:4*max_pd] = [1]*max_pd          # pd_stride: vectors of length 1
+    details[4*max_pd:5*max_pd] = [0]*max_pd          # pd_isvol: doens't matter if no norm
+    details[p+0*npars:p+1*npars] = range(2, npars+2) # par_offset
+    details[p+1*npars:p+2*npars] = [0]*npars         # no coordination
+    #details[p+npars] = 1 # par_coord[0] is coordinated with the first par?
+    details[p+2*npars:p+3*npars] = 0 # fast coord with 0
+    details[c]   = 1     # fast_coord_count: one fast index
+    details[c+1] = -1  # theta_var: None
+    details[c+2] = 0   # fast_theta: False
+    return details
+
+def poly_details(model_info, weights, constraints=None):
     if constraints is not None:
         # Need to find the independently varying pars and sort them
@@ -676 +678 @@
         # and weights, returning par blocks
         raise ValueError("Can't handle constraints yet")
+
+    max_pd = model_info['max_pd']
+    npars = len(model_info['parameters']) - 2
+    p = 5*max_pd
+    c = p + 3*npars
+
+    # Decreasing list of polydispersity lengths
+    # Note: the reversing view, x[::-1], does not require a copy
+    idx = np.argsort([len(w) for w in weights])[::-1]
+    details = np.zeros(c + 3, 'int32')
+    details[0*max_pd:1*max_pd] = idx[:max_pd]        # pd_par: arbitrary order; use first
+    details[1*max_pd:2*max_pd] = [1]*max_pd          # pd_length: only one element
+    details[2*max_pd:3*max_pd] = range(2, max_pd+2)  # pd_offset: skip scale and background
+    details[3*max_pd:4*max_pd] = [1]*max_pd          # pd_stride: vectors of length 1
+    details[4*max_pd:5*max_pd] = [0]*max_pd          # pd_isvol: doens't matter if no norm
+    details[p+0*npars:p+1*npars] = range(2, npars+2) # par_offset
+    details[p+1*npars:p+2*npars] = [0]*npars         # no coordination
+    #details[p+npars] = 1 # par_coord[0] is coordinated with the first par?
+    details[p+2*npars:p+3*npars] = 0 # fast coord with 0
+    details[c]   = 1     # fast_coord_count: one fast index
+    details[c+1] = -1  # theta_var: None
+    details[c+2] = 0   # fast_theta: False
+    return details
+
+
 
     raise ValueError("polydispersity not supported")
@@ -768 +795 @@
         oldpars=getattr(kernel_module, 'oldpars', {}),
         tests=getattr(kernel_module, 'tests', []),
-        mono_details = mono_details(MAX_PD, len(kernel_module.parameters))
         )
     process_parameters(model_info)
@@ -775 +801 @@
     model_info.update((k, getattr(kernel_module, k, None)) for k in functions)
     create_default_functions(model_info)
+    # Precalculate the monodisperse parameters
+    # TODO: make this a lazy evaluator
+    # make_model_info is called for every model on sasview startup
+    model_info['mono_details'] = mono_details(model_info)
     return model_info
 

sasmodels/kernel_iq.c

@@ -50 +50 @@
     )
 {
-printf("aliasing\n");
   // Storage for the current parameter values.  These will be updated as we
   // walk the polydispersity cube.
@@ -56 +55 @@
   double *pvec = (double *)(&local_pars);  // Alias named parameters with a vector
 
-printf("allocating\n");
   local int offset[NPARS-2];
 
 #if 1 // defined(USE_SHORTCUT_OPTIMIZATION)
-printf("dereferencing %p as %d\n", problem, problem->pd_length[0]);
   if (problem->pd_length[0] == 1) {
     // Shouldn't need to copy!!

sasmodels/kerneldll.py

@@ -50 +50 @@
 import tempfile
 import ctypes as ct
-from ctypes import c_void_p, c_int, c_longdouble, c_double, c_float
+from ctypes import c_void_p, c_int32, c_longdouble, c_double, c_float
 
 import numpy as np
@@ -205 +205 @@
 
         # int, int, int, int*, double*, double*, double*, double*, double*, double
-        argtypes = [c_int]*3 + [c_void_p]*5 + [fp]
+        argtypes = [c_int32]*3 + [c_void_p]*5 + [fp]
         self.Iq = self.dll[generate.kernel_name(self.info, False)]
         self.Iqxy = self.dll[generate.kernel_name(self.info, True)]

sasmodels/kernelpy.py

@@ -128 +128 @@
 
     def __call__(self, fixed, pd, cutoff=1e-5):
-        print("fixed",fixed)
-        print("pd", pd)
+        #print("fixed",fixed)
+        #print("pd", pd)
         args = self.args[:]  # grab a copy of the args
         form, form_volume = self.kernel, self.info['form_volume']
@@ -187 +187 @@
     ################################################################
 
-    #TODO: Wojtek's notes
-    #TODO: The goal is to restructure polydispersity loop
-    #so it allows fitting arbitrary polydispersity parameters
-    #and it accepts cases like coupled parameters
     weight = np.empty(len(pd), 'd')
     if weight.size > 0:
@@ -264 +260 @@
     result = scale * ret / norm + background
     return result
-
-"""
-Python driver for python kernels
-
-Calls the kernel with a vector of $q$ values for a single parameter set.
-Polydispersity is supported by looping over different parameter sets and
-summing the results.  The interface to :class:`PyModel` matches those for
-:class:`kernelcl.GpuModel` and :class:`kerneldll.DllModel`.
-"""
-import numpy as np
-from numpy import pi, cos
-
-from .generate import F64
-
-class PyModel(object):
-    """
-    Wrapper for pure python models.
-    """
-    def __init__(self, model_info):
-        self.info = model_info
-
-    def __call__(self, q_vectors):
-        q_input = PyInput(q_vectors, dtype=F64)
-        kernel = self.info['Iqxy'] if q_input.is_2d else self.info['Iq']
-        return PyKernel(kernel, self.info, q_input)
-
-    def release(self):
-        """
-        Free resources associated with the model.
-        """
-        pass
-
-class PyInput(object):
-    """
-    Make q data available to the gpu.
-
-    *q_vectors* is a list of q vectors, which will be *[q]* for 1-D data,
-    and *[qx, qy]* for 2-D data.  Internally, the vectors will be reallocated
-    to get the best performance on OpenCL, which may involve shifting and
-    stretching the array to better match the memory architecture.  Additional
-    points will be evaluated with *q=1e-3*.
-
-    *dtype* is the data type for the q vectors. The data type should be
-    set to match that of the kernel, which is an attribute of
-    :class:`GpuProgram`.  Note that not all kernels support double
-    precision, so even if the program was created for double precision,
-    the *GpuProgram.dtype* may be single precision.
-
-    Call :meth:`release` when complete.  Even if not called directly, the
-    buffer will be released when the data object is freed.
-    """
-    def __init__(self, q_vectors, dtype):
-        self.nq = q_vectors[0].size
-        self.dtype = dtype
-        self.is_2d = (len(q_vectors) == 2)
-        self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors]
-        self.q_pointers = [q.ctypes.data for q in self.q_vectors]
-
-    def release(self):
-        """
-        Free resources associated with the model inputs.
-        """
-        self.q_vectors = []
-
-class PyKernel(object):
-    """
-    Callable SAS kernel.
-
-    *kernel* is the DllKernel object to call.
-
-    *model_info* is the module information
-
-    *q_input* is the DllInput q vectors at which the kernel should be
-    evaluated.
-
-    The resulting call method takes the *pars*, a list of values for
-    the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight)
-    vectors for the polydisperse parameters.  *cutoff* determines the
-    integration limits: any points with combined weight less than *cutoff*
-    will not be calculated.
-
-    Call :meth:`release` when done with the kernel instance.
-    """
-    def __init__(self, kernel, model_info, q_input):
-        self.info = model_info
-        self.q_input = q_input
-        self.res = np.empty(q_input.nq, q_input.dtype)
-        dim = '2d' if q_input.is_2d else '1d'
-        # Loop over q unless user promises that the kernel is vectorized by
-        # taggining it with vectorized=True
-        if not getattr(kernel, 'vectorized', False):
-            if dim == '2d':
-                def vector_kernel(qx, qy, *args):
-                    """
-                    Vectorized 2D kernel.
-                    """
-                    return np.array([kernel(qxi, qyi, *args)
-                                     for qxi, qyi in zip(qx, qy)])
-            else:
-                def vector_kernel(q, *args):
-                    """
-                    Vectorized 1D kernel.
-                    """
-                    return np.array([kernel(qi, *args)
-                                     for qi in q])
-            self.kernel = vector_kernel
-        else:
-            self.kernel = kernel
-        fixed_pars = model_info['partype']['fixed-' + dim]
-        pd_pars = model_info['partype']['pd-' + dim]
-        vol_pars = model_info['partype']['volume']
-
-        # First two fixed pars are scale and background
-        pars = [p.name for p in model_info['parameters'][2:]]
-        offset = len(self.q_input.q_vectors)
-        self.args = self.q_input.q_vectors + [None] * len(pars)
-        self.fixed_index = np.array([pars.index(p) + offset
-                                     for p in fixed_pars[2:]])
-        self.pd_index = np.array([pars.index(p) + offset
-                                  for p in pd_pars])
-        self.vol_index = np.array([pars.index(p) + offset
-                                   for p in vol_pars])
-        try: self.theta_index = pars.index('theta') + offset
-        except ValueError: self.theta_index = -1
-
-        # Caller needs fixed_pars and pd_pars
-        self.fixed_pars = fixed_pars
-        self.pd_pars = pd_pars
-
-    def __call__(self, fixed, pd, cutoff=1e-5):
-        #print("fixed",fixed)
-        #print("pd", pd)
-        args = self.args[:]  # grab a copy of the args
-        form, form_volume = self.kernel, self.info['form_volume']
-        # First two fixed
-        scale, background = fixed[:2]
-        for index, value in zip(self.fixed_index, fixed[2:]):
-            args[index] = float(value)
-        res = _loops(form, form_volume, cutoff, scale, background, args,
-                     pd, self.pd_index, self.vol_index, self.theta_index)
-
-        return res
-
-    def release(self):
-        """
-        Free resources associated with the kernel.
-        """
-        self.q_input = None
-
-def _loops(form, form_volume, cutoff, scale, background,
-           args, pd, pd_index, vol_index, theta_index):
-    """
-    *form* is the name of the form function, which should be vectorized over
-    q, but otherwise have an interface like the opencl kernels, with the
-    q parameters first followed by the individual parameters in the order
-    given in model.parameters (see :mod:`sasmodels.generate`).
-
-    *form_volume* calculates the volume of the shape.  *vol_index* gives
-    the list of volume parameters
-
-    *cutoff* ignores the corners of the dispersion hypercube
-
-    *scale*, *background* multiplies the resulting form and adds an offset
-
-    *args* is the prepopulated set of arguments to the form function, starting
-    with the q vectors, and including slots for all the parameters.  The
-    values for the parameters will be substituted with values from the
-    dispersion functions.
-
-    *pd* is the list of dispersion parameters
-
-    *pd_index* are the indices of the dispersion parameters in *args*
-
-    *vol_index* are the indices of the volume parameters in *args*
-
-    *theta_index* is the index of the theta parameter for the sasview
-    spherical correction, or -1 if there is no angular dispersion
-    """
-
-    ################################################################
-    #                                                              #
-    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
-    #   !!                                                    !!   #
-    #   !!  KEEP THIS CODE CONSISTENT WITH KERNEL_TEMPLATE.C  !!   #
-    #   !!                                                    !!   #
-    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
-    #                                                              #
-    ################################################################
-
-    weight = np.empty(len(pd), 'd')
-    if weight.size > 0:
-        # weight vector, to be populated by polydispersity loops
-
-        # identify which pd parameters are volume parameters
-        vol_weight_index = np.array([(index in vol_index) for index in pd_index])
-
-        # Sort parameters in decreasing order of pd length
-        Npd = np.array([len(pdi[0]) for pdi in pd], 'i')
-        order = np.argsort(Npd)[::-1]
-        stride = np.cumprod(Npd[order])
-        pd = [pd[index] for index in order]
-        pd_index = pd_index[order]
-        vol_weight_index = vol_weight_index[order]
-
-        fast_value = pd[0][0]
-        fast_weight = pd[0][1]
-    else:
-        stride = np.array([1])
-        vol_weight_index = slice(None, None)
-        # keep lint happy
-        fast_value = [None]
-        fast_weight = [None]
-
-    ret = np.zeros_like(args[0])
-    norm = np.zeros_like(ret)
-    vol = np.zeros_like(ret)
-    vol_norm = np.zeros_like(ret)
-    for k in range(stride[-1]):
-        # update polydispersity parameter values
-        fast_index = k % stride[0]
-        if fast_index == 0:  # bottom loop complete ... check all other loops
-            if weight.size > 0:
-                for i, index, in enumerate(k % stride):
-                    args[pd_index[i]] = pd[i][0][index]
-                    weight[i] = pd[i][1][index]
-        else:
-            args[pd_index[0]] = fast_value[fast_index]
-            weight[0] = fast_weight[fast_index]
-
-        # Computes the weight, and if it is not sufficient then ignore this
-        # parameter set.
-        # Note: could precompute w1*...*wn so we only need to multiply by w0
-        w = np.prod(weight)
-        if w > cutoff:
-            # Note: can precompute spherical correction if theta_index is not
-            # the fast index. Correction factor for spherical integration
-            #spherical_correction = abs(cos(pi*args[phi_index])) if phi_index>=0 else 1.0
-            spherical_correction = (abs(cos(pi * args[theta_index])) * pi / 2
-                                    if theta_index >= 0 else 1.0)
-            #spherical_correction = 1.0
-
-            # Call the scattering function and adds it to the total.
-            I = form(*args)
-            if np.isnan(I).any(): continue
-            ret += w * I * spherical_correction
-            norm += w
-
-            # Volume normalization.
-            # If there are "volume" polydispersity parameters, then these
-            # will be used to call the form_volume function from the user
-            # supplied kernel, and accumulate a normalized weight.
-            # Note: can precompute volume norm if fast index is not a volume
-            if form_volume:
-                vol_args = [args[index] for index in vol_index]
-                vol_weight = np.prod(weight[vol_weight_index])
-                vol += vol_weight * form_volume(*vol_args)
-                vol_norm += vol_weight
-
-    positive = (vol * vol_norm != 0.0)
-    ret[positive] *= vol_norm[positive] / vol[positive]
-    result = scale * ret / norm + background
-    return result