Changeset e44432d in sasmodels

Timestamp:

Oct 19, 2018 3:46:26 PM (6 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

be43e39

Parents:

2586ab72

Message:

support hollow models in structure factor calculations

Location:

sasmodels

Files:

• generate.py (modified) (4 diffs)
• kernel.py (modified) (3 diffs)
• kernel_iq.c (modified) (11 diffs)
• kernelcl.py (modified) (1 diff)
• kerneldll.py (modified) (3 diffs)
• kernelpy.py (modified) (4 diffs)
• modelinfo.py (modified) (4 diffs)
• models/hollow_cylinder.c (modified) (5 diffs)
• models/hollow_rectangular_prism.c (modified) (2 diffs)
• models/hollow_rectangular_prism_thin_walls.c (modified) (2 diffs)
• models/vesicle.c (modified) (2 diffs)
• product.py (modified) (5 diffs)
• sasview_model.py (modified) (1 diff)

sasmodels/generate.py

@@ -715 +715 @@
     return False
 
+_SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", flags=re.MULTILINE)
+def contains_shell_volume(source):
+    # type: (List[str]) -> bool
+    """
+    Return True if C source defines "void Fq(".
+    """
+    for code in source:
+        m = _SHELL_VOLUME_PATTERN.search(code)
+        if m is not None:
+            return True
+    return False
+
 def _add_source(source, code, path, lineno=1):
     """
@@ -767 +779 @@
         pars = partable.form_volume_parameters
         source.append(_gen_fn(model_info, 'form_volume', pars))
+    if isinstance(model_info.shell_volume, str):
+        pars = partable.form_volume_parameters
+        source.append(_gen_fn(model_info, 'shell_volume', pars))
     if isinstance(model_info.Iq, str):
         pars = [q] + partable.iq_parameters
@@ -779 +794 @@
         pars = [qa, qb, qc] + partable.iq_parameters
         source.append(_gen_fn(model_info, 'Iqabc', pars))
+
+    # Check for shell_volume in source
+    is_hollow = contains_shell_volume(source)
 
     # What kind of 2D model do we need?  Is it consistent with the parameters?
@@ -802 +820 @@
                               for p in partable.kernel_parameters))
     # Define the function calls
-    call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(mode, v) 0.0"
+    call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) 0.0"
     if partable.form_volume_parameters:
         refs = _call_pars("_v.", partable.form_volume_parameters)
-        call_volume = "#define CALL_VOLUME(_v) form_volume(%s)"%(",".join(refs))
+        if is_hollow:
+            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = form_volume(%s); _shell = shell_volume(%s); } while (0)"%((",".join(refs),)*2)
+        else:
+            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = form_volume(%s); } while (0)"%(",".join(refs))
         if model_info.effective_radius_type:
-            call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(mode, _v) effective_radius(mode, %s)"%(",".join(refs))
+            call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) effective_radius(_mode, %s)"%(",".join(refs))
     else:
         # Model doesn't have volume.  We could make the kernel run a little
         # faster by not using/transferring the volume normalizations, but
         # the ifdef's reduce readability more than is worthwhile.
-        call_volume = "#define CALL_VOLUME(v) 1.0"
+        call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = 1.0; } while (0)"
     source.append(call_volume)
     source.append(call_effective_radius)

sasmodels/kernel.py

@@ -43 +43 @@
     def Iq(self, call_details, values, cutoff, magnetic):
         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
-        Pq, Reff = self.Pq_Reff(call_details, values, cutoff, magnetic, effective_radius_type=0)
-        return Pq
+        r"""
+        Returns I(q) from the polydisperse average scattering.
+
+        .. math::
+
+            I(q) = \text{scale} \cdot P(q) + \text{background}
+
+        With the correct choice of model and contrast, setting *scale* to
+        the volume fraction $V_f$ of particles should match the measured
+        absolute scattering.  Some models (e.g., vesicle) have volume fraction
+        built into the model, and do not need an additional scale.
+        """
+        _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff, magnetic,
+                              effective_radius_type=0)
+        combined_scale = values[0]/shell_volume
+        background = values[1]
+        return combined_scale*F2 + background
     __call__ = Iq
 
-    def Pq_Reff(self, call_details, values, cutoff, magnetic, effective_radius_type):
+    def Fq(self, call_details, values, cutoff, magnetic, effective_radius_type=0):
         # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray
+        r"""
+        Returns <F(q)>, <F(q)^2>, effective radius, shell volume and
+        form:shell volume ratio.  The <F(q)> term may be None if the
+        form factor does not support direct computation of $F(q)$
+
+        $P(q) = <F^2(q)>/<V>$ is used for structure factor calculations,
+
+        .. math::
+
+            I(q) = \text{scale} \cdot P(q) \cdot S(q) + \text{background}
+
+        For the beta approximation, this becomes
+
+        .. math::
+
+            I(q) = \text{scale} * P (1 + <F>^2/<F^2> (S - 1)) + \text{background}
+                 = \text{scale}/<V> (<F^2> + <F>^2 (S - 1)) + \text{background}
+
+        $<F(q)>$ and $<F^2(q)>$ are averaged by polydispersity in shape
+        and orientation, with each configuration $x_k$ having form factor
+        $F(q, x_k)$, weight $w_k$ and volume $V_k$.  The result is:
+
+        .. math::
+
+            P(q) = \frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k}
+
+        The form factor itself is scaled by volume and contrast to compute the
+        total scattering.  This is then squared, and the volume weighted
+        F^2 is then normalized by volume F.  For a given density, the number
+        of scattering centers is assumed to scale linearly with volume.  Later
+        scaling the resulting $P(q)$ by the volume fraction of particles
+        gives the total scattering on an absolute scale. Most models
+        incorporate the volume fraction into the overall scale parameter.  An
+        exception is vesicle, which includes the volume fraction parameter in
+        the model itself, scaling $F$ by $\surd V_f$ so that the math for
+        the beta approximation works out.
+
+        By scaling $P(q)$ by total weight $\sum w_k$, there is no need to make
+        sure that the polydisperisity distributions normalize to one.  In
+        particular, any distibution values $x_k$ outside the valid domain
+        of $F$ will not be included, and the distribution will be implicitly
+        truncated.  This is controlled by the parameter limits defined in the
+        model (which truncate the distribution before calling the kernel) as
+        well as any region excluded using the *INVALID* macro defined within
+        the model itself.
+
+        The volume used in the polydispersity calculation is the form volume
+        for solid objects or the shell volume for hollow objects.  Shell
+        volume should be used within $F$ so that the normalizing scale
+        represents the volume fraction of the shell rather than the entire
+        form.  This corresponds to the volume fraction of shell-forming
+        material added to the solvent.
+
+        The calculation of $S$ requires the effective radius and the
+        volume fraction of the particles.  The model can have several
+        different ways to compute effective radius, with the
+        *effective_radius_type* parameter used to select amongst them.  The
+        volume fraction of particles should be determined from the total
+        volume fraction of the form, not just the shell volume fraction.
+        This makes a difference for hollow shapes, which need to scale
+        the volume fraction by the returned volume ratio when computing $S$.
+        For solid objects, the shell volume is set to the form volume so
+        this scale factor evaluates to one and so can be used for both
+        hollow and solid shapes.
+        """
         self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type)
         #print("returned",self.q_input.q, self.result)
@@ -55 +135 @@
         if total_weight == 0.:
             total_weight = 1.
-        weighted_volume = self.result[nout*self.q_input.nq + 1]
-        weighted_radius = self.result[nout*self.q_input.nq + 2]
-        effective_radius = weighted_radius/total_weight
-        # compute I = scale*P + background
-        #           = scale*(sum(w*F^2)/sum w)/(sum (w*V)/sum w) + background
-        #           = scale/sum (w*V) * sum(w*F^2) + background
-        F2 = self.result[0:nout*self.q_input.nq:nout]
-        scale = values[0]/(weighted_volume if weighted_volume != 0.0 else 1.0)
-        background = values[1]
-        Pq = scale*F2 + background
-        #print("scale",scale,background)
-        return Pq, effective_radius
-
-    def beta(self, call_details, values, cutoff, magnetic, effective_radius_type):
-        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray
-        if self.dim == '2d':
-            raise NotImplementedError("beta not yet supported for 2D")
-        if not self.info.have_Fq:
-            raise NotImplementedError("beta not yet supported for "+self.info.id)
-        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type)
-        total_weight = self.result[2*self.q_input.nq + 0]
-        if total_weight == 0.:
-            total_weight = 1.
-        weighted_volume = self.result[2*self.q_input.nq + 1]
-        weighted_radius = self.result[2*self.q_input.nq + 2]
-        volume_average = weighted_volume/total_weight
-        effective_radius = weighted_radius/total_weight
-        F2 = self.result[0:2*self.q_input.nq:2]/total_weight
-        F1 = self.result[1:2*self.q_input.nq:2]/total_weight
-        return F1, F2, volume_average, effective_radius
+        # Note: shell_volume == form_volume for solid objects
+        form_volume = self.result[nout*self.q_input.nq + 1]/total_weight
+        shell_volume = self.result[nout*self.q_input.nq + 2]/total_weight
+        effective_radius = self.result[nout*self.q_input.nq + 3]/total_weight
+        if shell_volume == 0.:
+            shell_volume = 1.
+        F1 = self.result[1:nout*self.q_input.nq:nout]/total_weight if nout == 2 else None
+        F2 = self.result[0:nout*self.q_input.nq:nout]/total_weight
+        return F1, F2, effective_radius, shell_volume, form_volume/shell_volume
 
     def release(self):
@@ -92 +151 @@
     def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type):
         # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray
+        """
+        Call the kernel.  Subclasses defining kernels for particular execution
+        engines need to provide an implementation for this.
+        """
         raise NotImplementedError()

sasmodels/kernel_iq.c

@@ -26 +26 @@
 //  MAGNETIC_PARS : a comma-separated list of indices to the sld
 //      parameters in the parameter table.
-//  CALL_VOLUME(table) : call the form volume function
+//  CALL_VOLUME(form, shell, table) : assign form and shell values
 //  CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function
 //  CALL_IQ(q, table) : call the Iq function for 1D calcs.
@@ -275 +275 @@
 
 // ==================== KERNEL CODE ========================
-#define COMPUTE_F1_F2 defined(CALL_FQ)
 kernel
 void KERNEL_NAME(
@@ -345 +344 @@
   // version must loop over all q.
   #ifdef USE_OPENCL
-    #if COMPUTE_F1_F2
+    #if defined(CALL_FQ)
       double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]);
-      double weighted_volume = (pd_start == 0 ? 0.0 : result[2*nq+1]);
-      double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+2]);
+      double weighted_form = (pd_start == 0 ? 0.0 : result[2*nq+1]);
+      double weighted_shell = (pd_start == 0 ? 0.0 : result[2*nq+2]);
+      double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+3]);
       double this_F2 = (pd_start == 0 ? 0.0 : result[2*q_index+0]);
       double this_F1 = (pd_start == 0 ? 0.0 : result[2*q_index+1]);
     #else
       double weight_norm = (pd_start == 0 ? 0.0 : result[nq]);
-      double weighted_volume = (pd_start == 0 ? 0.0 : result[nq+1]);
-      double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+2]);
+      double weighted_form = (pd_start == 0 ? 0.0 : result[nq+1]);
+      double weighted_shell = (pd_start == 0 ? 0.0 : result[nq+2]);
+      double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+3]);
       double this_result = (pd_start == 0 ? 0.0 : result[q_index]);
     #endif
   #else // !USE_OPENCL
-    #if COMPUTE_F1_F2
+    #if defined(CALL_FQ)
       double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]);
-      double weighted_volume = (pd_start == 0 ? 0.0 : result[2*nq+1]);
-      double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+2]);
+      double weighted_form = (pd_start == 0 ? 0.0 : result[2*nq+1]);
+      double weighted_shell = (pd_start == 0 ? 0.0 : result[2*nq+2]);
+      double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+3]);
     #else
       double weight_norm = (pd_start == 0 ? 0.0 : result[nq]);
-      double weighted_volume = (pd_start == 0 ? 0.0 : result[nq+1]);
-      double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+2]);
+      double weighted_form = (pd_start == 0 ? 0.0 : result[nq+1]);
+      double weighted_shell = (pd_start == 0 ? 0.0 : result[nq+2]);
+      double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+3]);
     #endif
     if (pd_start == 0) {
@@ -371 +374 @@
       #pragma omp parallel for
       #endif
-      #if COMPUTE_F1_F2
+      #if defined(CALL_FQ)
           // 2*nq for F^2,F pairs
           for (int q_index=0; q_index < 2*nq; q_index++) result[q_index] = 0.0;
@@ -378 +381 @@
       #endif
     }
-    //if (q_index==0) printf("start %d %g %g\n", pd_start, weighted_volume, result[0]);
+    //if (q_index==0) printf("start %d %g %g\n", pd_start, weighted_shell, result[0]);
 #endif // !USE_OPENCL
 
@@ -693 +696 @@
     // Note: weight==0 must always be excluded
     if (weight > cutoff) {
-      weighted_volume += weight * CALL_VOLUME(local_values.table);
-      #if COMPUTE_F1_F2
+      double form, shell;
+      CALL_VOLUME(form, shell, local_values.table);
       weight_norm += weight;
-      #endif
+      weighted_form += weight * form;
+      weighted_shell += weight * shell;
       if (effective_radius_type != 0) {
         weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table);
@@ -747 +751 @@
           }
         #else  // !MAGNETIC
-          #if COMPUTE_F1_F2
+          #if defined(CALL_FQ)
             CALL_KERNEL(); // sets F1 and F2 by reference
           #else
@@ -756 +760 @@
 
         #ifdef USE_OPENCL
-          #if COMPUTE_F1_F2
+          #if defined(CALL_FQ)
             this_F2 += weight * F2;
             this_F1 += weight * F1;
@@ -763 +767 @@
           #endif
         #else // !USE_OPENCL
-          #if COMPUTE_F1_F2
+          #if defined(CALL_FQ)
             result[2*q_index+0] += weight * F2;
             result[2*q_index+1] += weight * F1;
@@ -793 +797 @@
 // Remember the current result and the updated norm.
 #ifdef USE_OPENCL
-  #if COMPUTE_F1_F2
+  #if defined(CALL_FQ)
     result[2*q_index+0] = this_F2;
     result[2*q_index+1] = this_F1;
     if (q_index == 0) {
       result[2*nq+0] = weight_norm;
-      result[2*nq+1] = weighted_volume;
-      result[2*nq+2] = weighted_radius;
+      result[2*nq+1] = weighted_form;
+      result[2*nq+3] = weighted_shell;
+      result[2*nq+3] = weighted_radius;
     }
   #else
@@ -805 +810 @@
     if (q_index == 0) {
       result[nq+0] = weight_norm;
-      result[nq+1] = weighted_volume;
-      result[nq+2] = weighted_radius;
+      result[nq+1] = weighted_form;
+      result[nq+2] = weighted_shell;
+      result[nq+3] = weighted_radius;
     }
   #endif
 
-//if (q_index == 0) printf("res: %g/%g\n", result[0], weigthed_volume);
+//if (q_index == 0) printf("res: %g/%g\n", result[0], weighted_shell);
 #else // !USE_OPENCL
-  #if COMPUTE_F1_F2
+  #if defined(CALL_FQ)
     result[2*nq] = weight_norm;
-    result[2*nq+1] = weighted_volume;
-    result[2*nq+2] = weighted_radius;
+    result[2*nq+1] = weighted_form;
+    result[2*nq+2] = weighted_shell;
+    result[2*nq+3] = weighted_radius;
   #else
     result[nq] = weight_norm;
-    result[nq+1] = weighted_volume;
-    result[nq+2] = weighted_radius;
+    result[nq+1] = weighted_form;
+    result[nq+2] = weighted_shell;
+    result[nq+3] = weighted_radius;
   #endif
-//printf("res: %g/%g\n", result[0], weighted_volume);
+//printf("res: %g/%g\n", result[0], weighted_shell);
 #endif // !USE_OPENCL
 
 // ** clear the macros in preparation for the next kernel **
-#undef COMPUTE_F1_F2
 #undef PD_INIT
 #undef PD_OPEN

sasmodels/kernelcl.py

@@ -539 +539 @@
         # leave room for f1/f2 results in case we need to compute beta for 1d models
         nout = 2 if self.info.have_Fq and self.dim == '1d' else 1
-        # plus 3 weight, volume, radius
-        self.result = np.empty(q_input.nq*nout + 3, self.dtype)
+        # +4 for total weight, shell volume, effective radius, form volume
+        self.result = np.empty(q_input.nq*nout + 4, self.dtype)
 
         # Inputs and outputs for each kernel call

sasmodels/kerneldll.py

@@ -99 +99 @@
     pass
 # pylint: enable=unused-import
+
+# Compiler output is a byte stream that needs to be decode in python 3
+decode = (lambda s: s) if sys.version_info[0] < 3 else (lambda s: s.decode('utf8'))
 
 if "SAS_DLL_PATH" in os.environ:
@@ -184 +187 @@
         subprocess.check_output(command, shell=shell, stderr=subprocess.STDOUT)
     except subprocess.CalledProcessError as exc:
-        raise RuntimeError("compile failed.\n%s\n%s"%(command_str, exc.output))
+        output = decode(exc.output)
+        raise RuntimeError("compile failed.\n%s\n%s"%(command_str, output))
     if not os.path.exists(output):
         raise RuntimeError("compile failed.  File is in %r"%source)
@@ -383 +387 @@
         # leave room for f1/f2 results in case we need to compute beta for 1d models
         nout = 2 if self.info.have_Fq else 1
-        # plus 3 weight, volume, radius
-        self.result = np.empty(q_input.nq*nout + 3, self.dtype)
+        # +4 for total weight, shell volume, effective radius, form volume
+        self.result = np.empty(q_input.nq*nout + 4, self.dtype)
         self.real = (np.float32 if self.q_input.dtype == generate.F32
                      else np.float64 if self.q_input.dtype == generate.F64

sasmodels/kernelpy.py

@@ -164 +164 @@
         self._volume_args = volume_args
         volume = model_info.form_volume
+        shell = model_info.shell_volume
         radius = model_info.effective_radius
-        self._volume = ((lambda: volume(*volume_args)) if volume
-                        else (lambda: 1.0))
+        self._volume = ((lambda: (shell(*volume_args), volume(*volume_args))) if shell and volume
+                        else (lambda: [volume(*volume_args)]*2) if volume
+                        else (lambda: (1.0, 1.0)))
         self._radius = ((lambda mode: radius(mode, *volume_args)) if radius
                         else (lambda mode: cbrt(0.75/pi*volume(*volume_args))) if volume
@@ -225 +227 @@
         total = form()
         weight_norm = 1.0
-        weighted_volume = form_volume()
+        weighted_shell, weighted_form = form_volume()
         weighted_radius = form_radius()
 
@@ -233 +235 @@
 
         weight_norm = 0.0
-        weighted_volume = 0.0
+        weighted_form = 0.0
+        weighted_shell = 0.0
         weighted_radius = 0.0
         partial_weight = np.NaN
@@ -272 +275 @@
                 total += weight * Iq
                 weight_norm += weight
-                weighted_volume += weight * form_volume()
+                shell, form = form_volume()
+                weighted_form += weight * form
+                weighted_shell += weight * shell
                 weighted_radius += weight * form_radius()
 
-    result = np.hstack((total, weight_norm, weighted_volume, weighted_radius))
+    result = np.hstack((total, weight_norm, weighted_form, weighted_shell, weighted_radius))
     return result
 

sasmodels/modelinfo.py

@@ -718 +718 @@
 
 #: Set of variables defined in the model that might contain C code
-C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code']
+C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'shell_volume', 'c_code']
 
 def _find_source_lines(model_info, kernel_module):
@@ -803 +803 @@
     info.tests = getattr(kernel_module, 'tests', [])
     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.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore
@@ -921 +922 @@
     #: 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.
+    effective_radius_type = None   # type: List[str]
     #: List of C source files used to define the model.  The source files
     #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the
     #: model defines orientation parameters. Files containing the most basic
     #: functions must appear first in the list, followed by the files that
-    #: use those functions.  Form factors are indicated by providing
-    #: an :attr:`ER` function.
-    effective_radius_type = None   # type: List[str]
-    #: Returns the occupied volume and the total volume for each parameter set.
-    #: See :attr:`ER` for details on the parameters.
+    #: use those functions.
     source = None           # type: List[str]
     #: The set of tests that must pass.  The format of the tests is described
@@ -958 +958 @@
     #: C code, either defined as a string, or in the sources.
     form_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]]
+    #: Returns the shell volume for python-based models.  Form volume and
+    #: shell volume are needed for volume normalization in the polydispersity
+    #: integral and structure interactions for hollow shapes.  If no
+    #: parameters are *volume* parameters, then shell volume is not needed.
+    #: For C-based models, (with :attr:`sources` defined, or with :attr:`Iq`
+    #: defined using a string containing C code), shell_volume must also be
+    #: C code, either defined as a string, or in the sources.
+    shell_volume = None      # type: Union[None, str, Callable[[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

sasmodels/models/hollow_cylinder.c

@@ -2 +2 @@
 
 static double
-_fq(double qab, double qc,
-    double radius, double thickness, double length)
+shell_volume(double radius, double thickness, double length)
 {
-    const double lam1 = sas_2J1x_x((radius+thickness)*qab);
-    const double lam2 = sas_2J1x_x(radius*qab);
-    const double gamma_sq = square(radius/(radius+thickness));
-    //Note: lim_{thickness -> 0} psi = sas_J0(radius*qab)
-    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab)
-    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq);    //SRK 10/19/00
-    const double t2 = sas_sinx_x(0.5*length*qc);
-    return psi*t2;
-}
-
-// TODO: interface to form_volume/shell_volume not yet settled
-static double
-shell_volume(double *total, double radius, double thickness, double length)
-{
-    *total = M_PI*length*square(radius+thickness);
-    return *total - M_PI*length*radius*radius;
+    return M_PI*length*(square(radius+thickness) - radius*radius);
 }
 
@@ -26 +10 @@
 form_volume(double radius, double thickness, double length)
 {
-    double total;
-    return shell_volume(&total, radius, thickness, length);
+    return M_PI*length*square(radius+thickness);
 }
 
@@ -62 +45 @@
 }
 
+static double
+_fq(double qab, double qc,
+    double radius, double thickness, double length)
+{
+    const double lam1 = sas_2J1x_x((radius+thickness)*qab);
+    const double lam2 = sas_2J1x_x(radius*qab);
+    const double gamma_sq = square(radius/(radius+thickness));
+    //Note: lim_{thickness -> 0} psi = sas_J0(radius*qab)
+    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab)
+    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq);    //SRK 10/19/00
+    const double t2 = sas_sinx_x(0.5*length*qc);
+    return psi*t2;
+}
+
 static void
 Fq(double q, double *F1, double *F2, double radius, double thickness, double length,
@@ -81 +78 @@
     total_F1 *= 0.5*(upper-lower);
     total_F2 *= 0.5*(upper-lower);
-    const double s = (sld - solvent_sld) * form_volume(radius, thickness, length);
+    const double s = (sld - solvent_sld) * shell_volume(radius, thickness, length);
     *F1 = 1e-2 * s * total_F1;
     *F2 = 1e-4 * s*s * total_F2;
@@ -93 +90 @@
 {
     const double form = _fq(qab, qc, radius, thickness, length);
-    const double s = (sld - solvent_sld) * form_volume(radius, thickness, length);
+    const double s = (sld - solvent_sld) * shell_volume(radius, thickness, length);
     return 1.0e-4*square(s * form);
 }

sasmodels/models/hollow_rectangular_prism.c

@@ -1 @@
-// TODO: interface to form_volume/shell_volume not yet settled
 static double
-shell_volume(double *total, double length_a, double b2a_ratio, double c2a_ratio, double thickness)
+shell_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness)
 {
-    double length_b = length_a * b2a_ratio;
-    double length_c = length_a * c2a_ratio;
-    double a_core = length_a - 2.0*thickness;
-    double b_core = length_b - 2.0*thickness;
-    double c_core = length_c - 2.0*thickness;
-    double vol_core = a_core * b_core * c_core;
-    *total = length_a * length_b * length_c;
-    return *total - vol_core;
+    const double length_b = length_a * b2a_ratio;
+    const double length_c = length_a * c2a_ratio;
+    const double form_volume = length_a * length_b * length_c;
+    const double a_core = length_a - 2.0*thickness;
+    const double b_core = length_b - 2.0*thickness;
+    const double c_core = length_c - 2.0*thickness;
+    const double core_volume = a_core * b_core * c_core;
+    return form_volume - core_volume;
 }
 
@@ -16 +15 @@
 form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness)
 {
-    double total;
-    return shell_volume(&total, length_a, b2a_ratio, c2a_ratio, thickness);
+    const double length_b = length_a * b2a_ratio;
+    const double length_c = length_a * c2a_ratio;
+    const double form_volume = length_a * length_b * length_c;
+    return form_volume;
 }
 
 static double
 effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio, double thickness)
-//effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c",
-//                         "equivalent outer circular cross-section","half ab diagonal","half diagonal"]
 // NOTE length_a is external dimension!
 {

sasmodels/models/hollow_rectangular_prism_thin_walls.c

@@ -1 @@
-// TODO: interface to form_volume/shell_volume not yet settled
 static double
-shell_volume(double *total, double length_a, double b2a_ratio, double c2a_ratio)
+shell_volume(double length_a, double b2a_ratio, double c2a_ratio)
 {
-    double length_b = length_a * b2a_ratio;
-    double length_c = length_a * c2a_ratio;
-    double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c);
-    *total = length_a * length_b * length_c;
-    return vol_shell;
+    const double length_b = length_a * b2a_ratio;
+    const double length_c = length_a * c2a_ratio;
+    const double shell_volume = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c);
+    return shell_volume;
 }
 
@@ -13 +11 @@
 form_volume(double length_a, double b2a_ratio, double c2a_ratio)
 {
-    double total;
-    return shell_volume(&total, length_a, b2a_ratio, c2a_ratio);
+    const double length_b = length_a * b2a_ratio;
+    const double length_c = length_a * c2a_ratio;
+    const double form_volume = length_a * length_b * length_c;
+    return form_volume;
 }
-
 
 static double

sasmodels/models/vesicle.c

@@ -1 @@
-// TODO: interface to form_volume/shell_volume not yet settled
 static double
-shell_volume(double *total, double radius, double thickness)
+shell_volume(double radius, double thickness)
 {
-    //note that for the vesicle model, the volume is ONLY the shell volume
-    *total = M_4PI_3 * cube(radius+thickness);
-    return *total - M_4PI_3*cube(radius);
+    return M_4PI_3 * (cube(radius+thickness) - cube(radius));
 }
 
@@ -11 +8 @@
 form_volume(double radius, double thickness)
 {
-    //note that for the vesicle model, the volume is ONLY the shell volume
-    double total;
-    return shell_volume(&total, radius, thickness);
+    return M_4PI_3 * cube(radius+thickness);
 }
+
 
 static double

sasmodels/product.py

@@ -168 +168 @@
     return par
 
-def _intermediates(P, S, effective_radius):
-    # type: (np.ndarray, np.ndarray, float) -> OrderedDict[str, np.ndarray]
-    """
-    Returns intermediate results for standard product (P(Q)*S(Q))
-    """
-    return OrderedDict((
-        ("P(Q)", P),
-        ("S(Q)", S),
-        ("effective_radius", effective_radius),
-    ))
-
-def _intermediates_beta(F1,              # type: np.ndarray
-                        F2,              # type: np.ndarray
-                        S,               # type: np.ndarray
-                        scale,           # type: np.ndarray
-                        bg,              # type: np.ndarray
-                        effective_radius # type: float
-                        ):
+def _intermediates(
+        F1,               # type: np.ndarray
+        F2,               # type: np.ndarray
+        S,                # type: np.ndarray
+        scale,            # type: float
+        effective_radius, # type: float
+        beta_mode,        # type: bool
+        ):
     # type: (...) -> OrderedDict[str, Union[np.ndarray, float]]
     """
-    Returns intermediate results for beta approximation-enabled product. The result may be an array or a float.
-    """
-    # TODO: 1. include calculated Q vector
-    # TODO: 2. consider implications if there are intermediate results in P(Q)
-    return OrderedDict((
-        ("P(Q)", scale*F2),
-        ("S(Q)", S),
-        ("beta(Q)", F1**2 / F2),
-        ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)),
-        ("effective_radius", effective_radius),
-        # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg),
-    ))
+    Returns intermediate results for beta approximation-enabled product.
+    The result may be an array or a float.
+    """
+    if beta_mode:
+        # TODO: 1. include calculated Q vector
+        # TODO: 2. consider implications if there are intermediate results in P(Q)
+        parts = OrderedDict((
+            ("P(Q)", scale*F2),
+            ("S(Q)", S),
+            ("beta(Q)", F1**2 / F2),
+            ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)),
+            ("effective_radius", effective_radius),
+            # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg),
+        ))
+    else:
+        parts = OrderedDict((
+            ("P(Q)", scale*F2),
+            ("S(Q)", S),
+            ("effective_radius", effective_radius),
+        ))
+    return parts
 
 class ProductModel(KernelModel):
@@ -253 +252 @@
     def __call__(self, call_details, values, cutoff, magnetic):
         # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray
+
         p_info, s_info = self.info.composition[1]
+        p_npars = p_info.parameters.npars
+        p_length = call_details.length[:p_npars]
+        p_offset = call_details.offset[:p_npars]
+        s_npars = s_info.parameters.npars
+        s_length = call_details.length[p_npars:p_npars+s_npars]
+        s_offset = call_details.offset[p_npars:p_npars+s_npars]
+
+        # Beta mode parameter is the first parameter after P and S parameters
+        have_beta_mode = p_info.have_Fq
+        beta_mode_offset = 2+p_npars+s_npars
+        beta_mode = (values[beta_mode_offset] > 0) if have_beta_mode else False
+        if beta_mode and self.p_kernel.dim== '2d':
+            raise NotImplementedError("beta not yet supported for 2D")
+
+        # R_eff type parameter is the second parameter after P and S parameters
+        # unless the model doesn't support beta mode, in which case it is first
+        have_radius_type = p_info.effective_radius_type is not None
+        radius_type_offset = 2+p_npars+s_npars + (1 if have_beta_mode else 0)
+        radius_type = int(values[radius_type_offset]) if have_radius_type else 0
+
+        # Retrieve the volume fraction, which is the second of the
+        # 'S' parameters in the parameter list, or 2+np in 0-origin,
+        # as well as the scale and background.
+        volfrac = values[3+p_npars]
+        scale, background = values[0], values[1]
 
         # if there are magnetic parameters, they will only be on the
@@ -268 +293 @@
 
         # Construct the calling parameters for P.
-        p_npars = p_info.parameters.npars
-        p_length = call_details.length[:p_npars]
-        p_offset = call_details.offset[:p_npars]
         p_details = make_details(p_info, p_length, p_offset, nweights)
-
-        # Set p scale to the volume fraction in s, which is the first of the
-        # 'S' parameters in the parameter list, or 2+np in 0-origin.
-        volfrac = values[2+p_npars]
-        p_values = [[volfrac, 0.0], values[2:2+p_npars], magnetism, weights]
+        p_values = [
+            [1., 0.], # scale=1, background=0,
+            values[2:2+p_npars],
+            magnetism,
+            weights]
         spacer = (32 - sum(len(v) for v in p_values)%32)%32
         p_values.append([0.]*spacer)
@@ -282 +304 @@
 
         # Construct the calling parameters for S.
-        s_npars = s_info.parameters.npars
-        s_length = call_details.length[p_npars:p_npars+s_npars]
-        s_offset = call_details.offset[p_npars:p_npars+s_npars]
-        s_length = np.hstack((1, s_length))
-        s_offset = np.hstack((nweights, s_offset))
-        s_details = make_details(s_info, s_length, s_offset, nweights+1)
+        if radius_type > 0:
+            # If R_eff comes from form factor, make sure it is monodisperse.
+            # weight is set to 1 later, after the value array is created
+            s_length[0] = 1
+        s_details = make_details(s_info, s_length, s_offset, nweights)
         s_values = [
-            # scale=1, background=0,
-            [1., 0.],
+            [1., 0.], # scale=1, background=0,
             values[2+p_npars:2+p_npars+s_npars],
             weights,
@@ -298 +318 @@
         s_values = np.hstack(s_values).astype(self.s_kernel.dtype)
 
-        # beta mode is the first parameter after the structure factor pars
-        extra_offset = 2+p_npars+s_npars
-        if p_info.have_Fq:
-            beta_mode = values[extra_offset]
-            extra_offset += 1
-        else:
-            beta_mode = 0
-        if p_info.effective_radius_type is not None:
-            effective_radius_type = int(values[extra_offset])
-            extra_offset += 1
-        else:
-            effective_radius_type = 0
-
-        # Call the kernels
-        scale, background = values[0], values[1]
-        if beta_mode:
-            F1, F2, volume_avg, effective_radius = self.p_kernel.beta(
-                p_details, p_values, cutoff, magnetic, effective_radius_type)
-            if effective_radius_type > 0:
-                # Plug R_eff from p into S model (initial value and pd value)
-                s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius
-            s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False)
-            combined_scale = scale*volfrac/volume_avg
-
-            self.results = lambda: _intermediates_beta(F1, F2, s_result, volfrac/volume_avg, background, effective_radius)
-            final_result = combined_scale*(F2 + (F1**2)*(s_result - 1)) + background
-
-        else:
-            p_result, effective_radius = self.p_kernel.Pq_Reff(
-                p_details, p_values, cutoff, magnetic, effective_radius_type)
-            if effective_radius_type > 0:
-                # Plug R_eff from p into S model (initial value and pd value)
-                s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius
-            s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False)
-            # remember the parts for plotting later
-            self.results = lambda: _intermediates(p_result, s_result, effective_radius)
-            final_result = scale*(p_result*s_result) + background
-
-        #call_details.show(values)
-        #print("values", values)
-        #p_details.show(p_values)
-        #print("=>", p_result)
-        #s_details.show(s_values)
-        #print("=>", s_result)
-        #import pylab as plt
-        #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-')
-        #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-')
-        #plt.figure()
+        # Call the form factor kernel to compute <F> and <F^2>.
+        # If the model doesn't support Fq the returned <F> will be None.
+        F1, F2, effective_radius, shell_volume, volume_ratio = self.p_kernel.Fq(
+            p_details, p_values, cutoff, magnetic, radius_type)
+
+        # Call the structure factor kernel to compute S.
+        # Plug R_eff from the form factor into structure factor parameters
+        # and scale volume fraction by form:shell volume ratio. These changes
+        # needs to be both in the initial value slot as well as the
+        # polydispersity distribution slot in the values array due to
+        # implementation details in kernel_iq.c.
+        #print("R_eff=%d:%g, volfrac=%g, volume ratio=%g"%(radius_type, effective_radius, volfrac, volume_ratio))
+        if radius_type > 0:
+            # set the value to the model ER and set the weight to 1
+            s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius
+            s_values[2+s_npars+s_offset[0]+nweights] = 1.0
+        s_values[3] = s_values[2+s_npars+s_offset[1]] = volfrac*volume_ratio
+        S = self.s_kernel.Iq(s_details, s_values, cutoff, False)
+
+        # Determine overall scale factor. Hollow shapes are weighted by
+        # shell_volume, so that is needed for volume normalization.  For
+        # solid shapes we can use shell_volume as well since it is equal
+        # to form volume.
+        combined_scale = scale*volfrac/shell_volume
+
+        # Combine form factor and structure factor
+        #print("beta", beta_mode, F1, F2, S)
+        PS = F2 + F1**2*(S-1) if beta_mode else F2*S
+        final_result = combined_scale*PS + background
+
+        # Capture intermediate values so user can see them.  These are
+        # returned as a lazy evaluator since they are only needed in the
+        # GUI, and not for each evaluation during a fit.
+        # TODO: return the results structure with the final results
+        # That way the model calcs are idempotent. Further, we can
+        # generalize intermediates to various other model types if we put it
+        # kernel calling interface.  Could do this as an "optional"
+        # return value in the caller, though in that case we could return
+        # the results directly rather than through a lazy evaluator.
+        self.results = lambda: _intermediates(
+            F1, F2, S, combined_scale, effective_radius, beta_mode)
+
         return final_result
 

sasmodels/sasview_model.py

@@ -847 +847 @@
     P = _make_standard_model('cylinder')()
     model = MultiplicationModel(P, S)
+    model.setParam('radius_effective_mode', 1.0)
     value = model.evalDistribution([0.1, 0.1])
     if np.isnan(value):