Changeset c1799d3 in sasmodels for sasmodels/generate.py

Timestamp:

Nov 9, 2018 2:16:05 PM (5 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:

c11d09f

Parents:

17695aa (diff), cf3d0ce (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 'beta_approx' into ticket-1157

File:

• sasmodels/generate.py (modified) (14 diffs)

sasmodels/generate.py

@@ -24 +24 @@
     dimension, or 1.0 if no volume normalization is required.
 
-    *ER(p1, p2, ...)* returns the effective radius of the form with
-    particular dimensions.
-
-    *VR(p1, p2, ...)* returns the volume ratio for core-shell style forms.
+    *shell_volume(p1, p2, ...)* returns the volume of the shell for forms
+    which are hollow.
+
+    *effective_radius(mode, p1, p2, ...)* returns the effective radius of
+    the form with particular dimensions.  Mode determines the type of
+    effective radius returned, with mode=1 for equivalent volume.
 
     #define INVALID(v) (expr)  returns False if v.parameter is invalid
@@ -72 +74 @@
 background.  This will work correctly even when polydispersity is off.
 
-*ER* and *VR* are python functions which operate on parameter vectors.
-The constructor code will generate the necessary vectors for computing
-them with the desired polydispersity.
 The kernel module must set variables defining the kernel meta data:
 
@@ -106 +105 @@
     create the OpenCL kernel functions.  The files defining the functions
     need to be listed before the files which use the functions.
-
-    *ER* is a python function defining the effective radius.  If it is
-    not present, the effective radius is 0.
-
-    *VR* is a python function defining the volume ratio.  If it is not
-    present, the volume ratio is 1.
 
     *form_volume*, *Iq*, *Iqac*, *Iqabc* are strings containing
@@ -671 +664 @@
 
 
-# type in IQXY pattern could be single, float, double, long double, ...
-_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]",
+_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ac|abc|xy))\s*[(]",
                            flags=re.MULTILINE)
 def find_xy_mode(source):
@@ -701 +693 @@
     return 'qa'
 
+# Note: not presently used.  Need to know whether Fq is available before
+# trying to compile the source.  Leave the code here in case we decide that
+# define have_Fq for each form factor is too tedious and error prone.
+_FQ_PATTERN = re.compile(r"(^|\s)void\s+Fq[(]", flags=re.MULTILINE)
+def contains_Fq(source):
+    # type: (List[str]) -> bool
+    """
+    Return True if C source defines "void Fq(".
+    """
+    for code in source:
+        m = _FQ_PATTERN.search(code)
+        if m is not None:
+            return True
+    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):
@@ -730 +748 @@
     # dispersion.  Need to be careful that necessary parameters are available
     # for computing volume even if we allow non-disperse volume parameters.
-
     partable = model_info.parameters
 
@@ -743 +760 @@
     for path, code in user_code:
         _add_source(source, code, path)
-
     if model_info.c_code:
         _add_source(source, model_info.c_code, model_info.filename,
@@ -751 +767 @@
     q, qx, qy, qab, qa, qb, qc \
         = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()]
+
     # Generate form_volume function, etc. from body only
     if isinstance(model_info.form_volume, str):
         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
@@ -767 +787 @@
         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?
@@ -789 +812 @@
     source.append("\\\n".join(p.as_definition()
                               for p in partable.kernel_parameters))
-
     # Define the function calls
+    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))
     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)
     model_refs = _call_pars("_v.", partable.iq_parameters)
-    pars = ",".join(["_q"] + model_refs)
-    call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars
+
+    if model_info.have_Fq:
+        pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs)
+        call_iq = "#define CALL_FQ(_q, _F1, _F2, _v) Fq(%s)" % pars
+        clear_iq = "#undef CALL_FQ"
+    else:
+        pars = ",".join(["_q"] + model_refs)
+        call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars
+        clear_iq = "#undef CALL_IQ"
     if xy_mode == 'qabc':
         pars = ",".join(["_qa", "_qb", "_qc"] + model_refs)
@@ -812 +847 @@
         call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars
         clear_iqxy = "#undef CALL_IQ_AC"
-    elif xy_mode == 'qa':
+    elif xy_mode == 'qa' and not model_info.have_Fq:
         pars = ",".join(["_qa"] + model_refs)
         call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars
         clear_iqxy = "#undef CALL_IQ_A"
+    elif xy_mode == 'qa' and model_info.have_Fq:
+        pars = ",".join(["_qa", "&_F1", "&_F2",] + model_refs)
+        # Note: uses rare C construction (expr1, expr2) which computes
+        # expr1 then expr2 and evaluates to expr2.  This allows us to
+        # leave it looking like a function even though it is returning
+        # its values by reference.
+        call_iqxy = "#define CALL_FQ_A(_qa,_F1,_F2,_v) (Fq(%s),_F2)" % pars
+        clear_iqxy = "#undef CALL_FQ_A"
     elif xy_mode == 'qxy':
         orientation_refs = _call_pars("_v.", partable.orientation_parameters)
@@ -831 +874 @@
     magpars = [k-2 for k, p in enumerate(partable.call_parameters)
                if p.type == 'sld']
-
     # Fill in definitions for numbers of parameters
     source.append("#define MAX_PD %s"%partable.max_pd)
@@ -839 +881 @@
     source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars))
     source.append("#define PROJECTION %d"%PROJECTION)
-
     # TODO: allow mixed python/opencl kernels?
-
-    ocl = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name)
-    dll = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name)
+    ocl = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name)
+    dll = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name)
+
     result = {
         'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]),
         'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]),
     }
-
     return result
 
 
-def _kernels(kernel, call_iq, call_iqxy, clear_iqxy, name):
+def _kernels(kernel, call_iq, clear_iq, call_iqxy, clear_iqxy, name):
     # type: ([str,str], str, str, str) -> List[str]
     code = kernel[0]
@@ -862 +902 @@
         '#line 1 "%s Iq"' % path,
         code,
-        "#undef CALL_IQ",
+        clear_iq,
         "#undef KERNEL_NAME",
         ]