Changeset 5399809 in sasmodels

Timestamp:

Aug 21, 2018 1:32:40 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:

2a12351b

Parents:

c57ee9e

Message:

fix R_eff support infrastructure so more tests pass

Location:

sasmodels

Files:

• core.py (modified) (1 diff)
• generate.py (modified) (2 diffs)
• kernel.py (modified) (3 diffs)
• kernelcl.py (modified) (3 diffs)
• kernelpy.py (modified) (2 diffs)
• modelinfo.py (modified) (1 diff)
• product.py (modified) (2 diffs)

sasmodels/core.py

@@ -363 +363 @@
     model = load_model("cylinder@hardsphere*sphere")
     actual = [p.name for p in model.info.parameters.kernel_parameters]
-    target = ("sld sld_solvent radius length theta phi volfraction"
-              " beta_mode A_sld A_sld_solvent A_radius").split()
+    target = ("sld sld_solvent radius length theta phi"
+              " radius_effective volfraction structure_factor_mode"
+              " A_sld A_sld_solvent A_radius").split()
     assert target == actual, "%s != %s"%(target, actual)
 

sasmodels/generate.py

@@ -802 +802 @@
                               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))
-        call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(mode, _v) effective_radius(mode, %s)"%(",".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
@@ -811 +813 @@
         # the ifdef's reduce readability more than is worthwhile.
         call_volume = "#define CALL_VOLUME(v) 1.0"
-        call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(mode, v) 0.0"
     source.append(call_volume)
     source.append(call_effective_radius)

sasmodels/kernel.py

@@ -44 +44 @@
         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
         Pq, Reff = self.Pq_Reff(call_details, values, cutoff, magnetic, effective_radius_type=0)
-        scale = values[0]
-        background = values[1]
-        return scale*Pq + background
+        return Pq
     __call__ = Iq
 
@@ -53 +51 @@
         self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type)
         #print("returned",self.q_input.q, self.result)
-        nout = 2 if self.info.have_Fq else 1
+        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1
         total_weight = self.result[nout*self.q_input.nq + 0]
         if total_weight == 0.:
@@ -63 +61 @@
         #           = 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 self.result[0:nout*self.q_input.nq:nout], effective_radius
+        return Pq, effective_radius
 
     def beta(self, call_details, values, cutoff, magnetic, effective_radius_type):

sasmodels/kernelcl.py

@@ -481 +481 @@
         # at this point, so instead using 32, which is good on the set of
         # architectures tested so far.
+        extra_q = 3  # total weight, weighted volume and weighted radius
         if self.is_2d:
-            # Note: 16 rather than 15 because result is 1 longer than input.
-            width = ((self.nq+16)//16)*16
+            width = ((self.nq+15+extra_q)//16)*16
             self.q = np.empty((width, 2), dtype=dtype)
             self.q[:self.nq, 0] = q_vectors[0]
             self.q[:self.nq, 1] = q_vectors[1]
         else:
-            # Note: 32 rather than 31 because result is 1 longer than input.
-            width = ((self.nq+32)//32)*32
+            width = ((self.nq+31+extra_q)//32)*32
             self.q = np.empty(width, dtype=dtype)
             self.q[:self.nq] = q_vectors[0]
@@ -539 +538 @@
         self.dim = '2d' if q_input.is_2d else '1d'
         # leave room for f1/f2 results in case we need to compute beta for 1d models
-        num_returns = 1 if self.dim == '2d' else 2  #
-        # plus 1 for the normalization value, plus another for R_eff
-        self.result = np.empty((q_input.nq+2)*num_returns, dtype)
+        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)
 
         # Inputs and outputs for each kernel call
@@ -549 +548 @@
 
         self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
-                                  q_input.global_size[0] * num_returns * dtype.itemsize)
+                                  q_input.global_size[0] * nout * dtype.itemsize)
         self.q_input = q_input # allocated by GpuInput above
 

sasmodels/kernelpy.py

@@ -182 +182 @@
         radius = ((lambda: 0.0) if effective_radius_type == 0
                   else (lambda: self._radius(effective_radius_type)))
-        res = _loops(self._parameter_vector, self._form, self._volume, radius,
-                     self.q_input.nq, call_details, values, cutoff)
-        return res
+        self.result = _loops(
+            self._parameter_vector, self._form, self._volume, radius,
+            self.q_input.nq, call_details, values, cutoff)
 
     def release(self):
@@ -213 +213 @@
     #                                                              #
     ################################################################
+
+    # WARNING: Trickery ahead
+    # The parameters[] vector is embedded in the closures for form(),
+    # form_volume() and form_radius().  We set the initial vector from
+    # the values for the model parameters. As we loop through the polydispesity
+    # mesh, we update the components with the polydispersity values before
+    # calling the respective functions.
     n_pars = len(parameters)
     parameters[:] = values[2:n_pars+2]
+
     if call_details.num_active == 0:
-        pd_norm = float(form_volume())
-        scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0)
-        background = values[1]
-        return scale*form() + background
-
-    pd_value = values[2+n_pars:2+n_pars + call_details.num_weights]
-    pd_weight = values[2+n_pars + call_details.num_weights:]
-
-    weight_norm = 0.0
-    weighted_volume = 0.0
-    weighted_radius = 0.0
-    partial_weight = np.NaN
-    weight = np.NaN
-
-    p0_par = call_details.pd_par[0]
-    p0_length = call_details.pd_length[0]
-    p0_index = p0_length
-    p0_offset = call_details.pd_offset[0]
-
-    pd_par = call_details.pd_par[:call_details.num_active]
-    pd_offset = call_details.pd_offset[:call_details.num_active]
-    pd_stride = call_details.pd_stride[:call_details.num_active]
-    pd_length = call_details.pd_length[:call_details.num_active]
-
-    total = np.zeros(nq, 'd')
-    for loop_index in range(call_details.num_eval):
-        # update polydispersity parameter values
-        if p0_index == p0_length:
-            pd_index = (loop_index//pd_stride)%pd_length
-            parameters[pd_par] = pd_value[pd_offset+pd_index]
-            partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:])
-            p0_index = loop_index%p0_length
-
-        weight = partial_weight * pd_weight[p0_offset + p0_index]
-        parameters[p0_par] = pd_value[p0_offset + p0_index]
-        p0_index += 1
-        if weight > cutoff:
-            # Call the scattering function
-            # Assume that NaNs are only generated if the parameters are bad;
-            # exclude all q for that NaN.  Even better would be to have an
-            # INVALID expression like the C models, but that is too expensive.
-            Iq = np.asarray(form(), 'd')
-            if np.isnan(Iq).any():
-                continue
-
-            # update value and norm
-            total += weight * Iq
-            weight_norm += weight
-            weighted_volume += weight * form_volume()
-            weighted_radius += weight * form_radius()
+        total = form()
+        weight_norm = 1.0
+        weighted_volume = form_volume()
+        weighted_radius = form_radius()
+
+    else:
+        pd_value = values[2+n_pars:2+n_pars + call_details.num_weights]
+        pd_weight = values[2+n_pars + call_details.num_weights:]
+
+        weight_norm = 0.0
+        weighted_volume = 0.0
+        weighted_radius = 0.0
+        partial_weight = np.NaN
+        weight = np.NaN
+
+        p0_par = call_details.pd_par[0]
+        p0_length = call_details.pd_length[0]
+        p0_index = p0_length
+        p0_offset = call_details.pd_offset[0]
+
+        pd_par = call_details.pd_par[:call_details.num_active]
+        pd_offset = call_details.pd_offset[:call_details.num_active]
+        pd_stride = call_details.pd_stride[:call_details.num_active]
+        pd_length = call_details.pd_length[:call_details.num_active]
+
+        total = np.zeros(nq, 'd')
+        for loop_index in range(call_details.num_eval):
+            # update polydispersity parameter values
+            if p0_index == p0_length:
+                pd_index = (loop_index//pd_stride)%pd_length
+                parameters[pd_par] = pd_value[pd_offset+pd_index]
+                partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:])
+                p0_index = loop_index%p0_length
+
+            weight = partial_weight * pd_weight[p0_offset + p0_index]
+            parameters[p0_par] = pd_value[p0_offset + p0_index]
+            p0_index += 1
+            if weight > cutoff:
+                # Call the scattering function
+                # Assume that NaNs are only generated if the parameters are bad;
+                # exclude all q for that NaN.  Even better would be to have an
+                # INVALID expression like the C models, but that is expensive.
+                Iq = np.asarray(form(), 'd')
+                if np.isnan(Iq).any():
+                    continue
+
+                # update value and norm
+                total += weight * Iq
+                weight_norm += weight
+                weighted_volume += weight * form_volume()
+                weighted_radius += weight * form_radius()
 
     result = np.hstack((total, weight_norm, weighted_volume, weighted_radius))

sasmodels/modelinfo.py

@@ -794 +794 @@
     info.source = getattr(kernel_module, 'source', [])
     info.c_code = getattr(kernel_module, 'c_code', None)
+    info.effective_radius = getattr(kernel_module, 'effective_radius', None)
     info.ER = None  # CRUFT
     info.VR = None  # CRUFT

sasmodels/product.py

@@ -301 +301 @@
                 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)
@@ -312 +313 @@
                 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)