Allow pure python Fq for beta approximation.

sasmodels/kernelpy.py

@@ -161,10 +161,17 @@ def __init__(self, model_info, q_input):
             form = model_info.Iqxy
             qx, qy = q_input.q[:, 0], q_input.q[:, 1]
             self._form = lambda: form(qx, qy, *kernel_args)
+            self._nq_copies = 1
         else:
-            form = model_info.Iq
             q = q_input.q
-            self._form = lambda: form(q, *kernel_args)
+            if model_info.have_Fq:
+                form = model_info.Fq
+                self._form = lambda: _pack_fq(*form(q, *kernel_args))
+                self._nq_copies = 2
+            else:
+                form = model_info.Iq
+                self._form = lambda: form(q, *kernel_args)
+                self._nq_copies = 1
 
         # Generate a closure which calls the form_volume if it exists.
         self._volume_args = volume_args
@@ -188,7 +195,7 @@ def _call_kernel(self, call_details, values, cutoff, magnetic, radius_effective_
                   else (lambda: self._radius(radius_effective_mode)))
         self.result = _loops(
             self._parameter_vector, self._form, self._volume, radius,
-            self.q_input.nq, call_details, values, cutoff)
+            self.q_input.nq*self._nq_copies, call_details, values, cutoff)
 
     def release(self):
         # type: () -> None
@@ -198,6 +205,13 @@ def release(self):
         self.q_input.release()
         self.q_input = None
 
+def _pack_fq(fq, fqsq):
+    result = np.empty(2*fq.shape[0], dtype=fq.dtype)
+    # PAK: For reasons unknown, the original author (me) decided to pack the
+    # values into the returned array as (f^2, f) pairs.
+    result[0::2] = fqsq
+    result[1::2] = fq
+    return result
 
 def _loops(parameters, form, form_volume, form_radius, nq, call_details,
            values, cutoff):
@@ -291,6 +305,7 @@ def _create_default_functions(model_info):
     """
     # Note: Must call create_vector_Iq before create_vector_Iqxy.
     _create_vector_Iq(model_info)
+    _create_vector_Fq(model_info)
     _create_vector_Iqxy(model_info)
 
 
@@ -309,6 +324,22 @@ def vector_Iq(q, *args):
         vector_Iq.vectorized = True
         model_info.Iq = vector_Iq
 
+def _create_vector_Fq(model_info):
+    """
+    Define Fq as a vector function if it exists.
+    """
+    # Note that this is doing slightly too much work since we are composing
+    # two vector results which are then interlaced via _pack_fq above.
+    Fq = model_info.Fq
+    if callable(Fq) and not getattr(Fq, 'vectorized', False):
+        def vector_Fq(q, *args):
+            """
+            Vectorized 1D kernel returning Fq, Fq**2
+            """
+            fq, fqsq = zip(*(Fq(qi, *args) for qi in q))
+            return np.array(fq), np.array(fqsq)
+        vector_Fq.vectorized = True
+        model_info.Fq = vector_Fq
 
 def _create_vector_Iqxy(model_info):
     """

sasmodels/modelinfo.py

@@ -943,9 +943,7 @@ def make_model_info(kernel_module):
     info.docs = kernel_module.__doc__
     info.category = getattr(kernel_module, 'category', None)
     info.structure_factor = getattr(kernel_module, 'structure_factor', False)
-    # TODO: find Fq by inspection
     info.radius_effective_modes = getattr(kernel_module, 'radius_effective_modes', None)
-    info.have_Fq = getattr(kernel_module, 'have_Fq', False)
     info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y'])
     # Note: custom.load_custom_kernel_module assumes the C sources are defined
     # by this attribute.
@@ -958,6 +956,9 @@ def make_model_info(kernel_module):
     info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore
     info.shell_volume = getattr(kernel_module, 'shell_volume', None) # type: ignore
     info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore
+    info.Fq = getattr(kernel_module, 'Fq', None) # type: ignore
+    # TODO: We should be able to find Fq in C code by inspection.
+    info.have_Fq = getattr(kernel_module, 'have_Fq', (info.Fq is not None))
     info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore
     info.Iqac = getattr(kernel_module, 'Iqac', None) # type: ignore
     info.Iqabc = getattr(kernel_module, 'Iqabc', None) # type: ignore
@@ -1096,9 +1097,6 @@ class ModelInfo(object):
     #: between form factor models.  This will default to False if it is not
     #: provided in the file.
     structure_factor = None # type: bool
-    #: True if the model defines an Fq function with signature
-    #: ``void Fq(double q, double *F1, double *F2, ...)``
-    have_Fq = False
     #: List of options for computing the effective radius of the shape,
     #: or None if the model is not usable as a form factor model.
     radius_effective_modes = None   # type: List[str]
@@ -1136,20 +1134,16 @@ class ModelInfo(object):
     #: monodisperse approximation for non-dilute solutions, P@S.  The first
     #: argument is the integer effective radius mode, with default 0.
     radius_effective = None  # type: Union[None, Callable[[int, np.ndarray], float]]
-    #: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined
-    #: by the parameter table.  *Iq* can be defined as a python function, or
-    #: as a C function.  If it is defined in C, then set *Iq* to the body of
-    #: the C function, including the return statement.  This function takes
-    #: values for *q* and each of the parameters as separate *double* values
-    #: (which may be converted to float or long double by sasmodels).  All
-    #: source code files listed in :attr:`source` will be loaded before the
-    #: *Iq* function is defined.  If *Iq* is not present, then sources should
-    #: define *static double Iq(double q, double a, double b, ...)* which
-    #: will return *I(q, a, b, ...)*.  Multiplicity parameters are sent as
-    #: pointers to doubles.  Constants in floating point expressions should
-    #: include the decimal point. See :mod:`.generate` for more details. If
-    #: *have_Fq* is True, then Iq should return an interleaved array of
-    #: $[\sum F(q_1), \sum F^2(q_1), \ldots, \sum F(q_n), \sum F^2(q_n)]$.
+    #: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined by the
+    #: parameter table. Multiplicity parameters such as the number of shells are
+    #: sent as floating point values. If the function can operate with a vector
+    #: of *q* values (that is, you aren't doing simple comparisons such as
+    #: *q > 0*), then set *Iq.vectorized = True* to make your code run faster.
+    #: You can also set *Iq* to a string containing the body of a C function.
+    #: For example, ``Iq = return (a*q + b)*q + c;`` will generate a quadratic.
+    #: All source code files listed in :attr:`source` will be loaded before the
+    #: C function is defined. Constants in floating point expressions should
+    #: include the decimal point. See :mod:`.generate` for more details.
     Iq = None               # type: Union[None, str, Callable[[...], np.ndarray]]
     #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`.
     Iqxy = None             # type: Union[None, str, Callable[[...], np.ndarray]]
@@ -1159,6 +1153,13 @@ class ModelInfo(object):
     Iqabc = None            # type: Union[None, str, Callable[[...], np.ndarray]]
     #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`.
     Imagnetic = None        # type: Union[None, str, Callable[[...], np.ndarray]]
+    #: Returns *F(q, a, b, ...), F^2(q, a, b, c, ...)*. Note: you cannot assign
+    #: a C source code body to *Fq*.
+    Fq = None               # type: Union[None, Callable[[...], Tuple[np.ndarray,np.ndarray]]]
+    #: True if the model defines an Fq function in C with signature
+    #: ``void Fq(double q, double *F1, double *F2, ...)``
+    #: or in python as ``Fq(q, ...) -> (fq, fq^2)``.
+    have_Fq = False
     #: Returns a model profile curve *x, y*.  If *profile* is defined, this
     #: curve will appear in response to the *Show* button in SasView.  Use
     #: :attr:`profile_axes` to set the axis labels.  Note that *y* values

sasmodels/models/_spherepy.py

@@ -79,23 +79,29 @@ def radius_effective(mode, radius):
     """Calculate R_eff for sphere"""
     return radius if mode else 0.
 
-def Iq(q, sld, sld_solvent, radius):
-    """Calculate I(q) for sphere"""
+# Variants for testing purposes: toggle True/False
+vectorized = True  # Whether to call a q vector or loop over q scalars.
+have_Fq = True  # Whether to use Fq or Iq for evaluation.
+def Fq(q, sld, sld_solvent, radius):
+    """Calculate F(q), F^2(q) for sphere"""
     #print "q",q
     #print "sld,r",sld,sld_solvent,radius
     qr = q * radius
     sn, cn = sin(qr), cos(qr)
-    ## The natural expression for the Bessel function is the following:
-    ##     bes = 3 * (sn-qr*cn)/qr**3 if qr>0 else 1
-    ## however, to support vector q values we need to handle the conditional
-    ## as a vector, which we do by first evaluating the full expression
-    ## everywhere, then fixing it up where it is broken.   We should probably
-    ## set numpy to ignore the 0/0 error before we do though...
-    bes = 3 * (sn - qr * cn) / qr ** 3 # may be 0/0 but we fix that next line
-    bes[qr == 0] = 1
-    fq = bes * (sld - sld_solvent) * form_volume(radius)
-    return 1.0e-4 * fq ** 2
-Iq.vectorized = True  # Iq accepts an array of q values
+    if vectorized:
+        with np.errstate(all='ignore'):
+            bes = 3 * (sn - qr * cn) / qr ** 3 # may be 0/0 but we fix that next line
+            bes[qr == 0] = 1
+    else:
+        bes = 3 * (sn-qr*cn)/qr**3 if qr != 0 else 1
+    fq = bes * (1e-2 * (sld - sld_solvent) * form_volume(radius))
+    return fq, fq**2
+Fq.vectorized = vectorized  # Fq accepts an array of q value
+
+def Iq(q, sld, sld_solvent, radius):
+    """Calculate I(q) for sphere"""
+    return Fq(q, sld, sld_solvent, radius)[1]
+Iq.vectorized = vectorized  # Iq accepts an array of q value
 
 def sesans(z, sld, sld_solvent, radius):
     """