Changeset 6dc78e4 in sasmodels

Timestamp:

Aug 17, 2016 5:12:51 PM (8 years ago)

Author:

Paul Kienzle <pkienzle@…>

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:

deb854f, 300a2f7

Parents:

6b7e2f7

Message:

re-enable structure factor product

Location:

sasmodels

Files:

• convert.py (modified) (1 diff)
• details.py (modified) (3 diffs)
• mixture.py (modified) (2 diffs)
• product.py (modified) (9 diffs)

sasmodels/convert.py

@@ -180 +180 @@
         composition_type, parts = model_info.composition
         if composition_type == 'product':
-            oldpars = {'scale':'scale_factor'}
-            oldpars.update(_get_old_pars(parts[0]))
-            oldpars.update(_get_old_pars(parts[1]))
+            translation = {'scale':'scale_factor'}
+            translation.update(_get_translation_table(parts[0]))
+            translation.update(_get_translation_table(parts[1]))
         else:
             raise NotImplementedError("cannot convert to sasview sum")
     else:
         translation = _get_translation_table(model_info)
-        oldpars = _revert_pars(_unscale_sld(model_info, pars), translation)
-        oldpars = _trim_vectors(model_info, pars, oldpars)
+    oldpars = _revert_pars(_unscale_sld(model_info, pars), translation)
+    oldpars = _trim_vectors(model_info, pars, oldpars)
 
     # Make sure the control parameter is an integer

sasmodels/details.py

@@ -177 +177 @@
 
 def make_details(model_info, length, offset, num_weights):
+    # type: (ModelInfo, np.ndarray, np.ndarray, int) -> CallDetails
     """
     Return a :class:`CallDetails` object for a polydisperse calculation
@@ -186 +187 @@
     array.
     """
-    # type: (ModelInfo, np.ndarray, np.ndarray, int) -> CallDetails
     #pars = model_info.parameters.call_parameters[2:model_info.parameters.npars+2]
     #print(", ".join(str(i)+"-"+p.id for i,p in enumerate(pars)))
@@ -275 +275 @@
     """
     value, weight = zip(*pars)
-    weight = [w if w else [1.] for w in weight]
+    #weight = [w if len(w)>0 else [1.] for w in weight]
     weight = np.vstack([v.flatten() for v in meshgrid(*weight)])
     weight = np.prod(weight, axis=0)

sasmodels/mixture.py

@@ -116 +116 @@
         self.kernels = kernels
         self.dtype = self.kernels[0].dtype
+        self.results = []  # type: List[np.ndarray]
 
     def __call__(self, call_details, values, cutoff, magnetic):
@@ -122 +123 @@
         total = 0.0
         # remember the parts for plotting later
-        self.results = []
+        self.results = []  # type: List[np.ndarray]
         offset = 2 # skip scale & background
         parts = MixtureParts(self.info, self.kernels, call_details, values)

sasmodels/product.py

@@ -9 +9 @@
 
 To use it, first load form factor P and structure factor S, then create
-*ProductModel(P, S)*.
+*make_product_info(P, S)*.
 """
+from __future__ import print_function, division
+
+from copy import copy
 import numpy as np  # type: ignore
 
-from .modelinfo import suffix_parameter, ParameterTable, ModelInfo
+from .modelinfo import Parameter, ParameterTable, ModelInfo
 from .kernel import KernelModel, Kernel
-from .details import dispersion_mesh
+from .details import make_details, dispersion_mesh
 
 try:
     from typing import Tuple
-    from .modelinfo import ParameterSet
-    from .details import CallDetails
 except ImportError:
     pass
@@ -26 +27 @@
 # TODO: make estimates available to constraints
 #ESTIMATED_PARAMETERS = [
-#    ["est_effect_radius", "A", 0.0, [0, np.inf], "", "Estimated effective radius"],
+#    ["est_radius_effective", "A", 0.0, [0, np.inf], "", "Estimated effective radius"],
 #    ["est_volume_ratio", "", 1.0, [0, np.inf], "", "Estimated volume ratio"],
 #]
+
+ER_ID = "radius_effective"
+VF_ID = "volfraction"
 
 # TODO: core_shell_sphere model has suppressed the volume ratio calculation
@@ -37 +41 @@
     Create info block for product model.
     """
+    # Make sure effective radius is the first parameter and
+    # make sure volume fraction is the second parameter of the
+    # structure factor calculator.  Structure factors should not
+    # have any magnetic parameters
+    assert(s_info.parameters.kernel_parameters[0].id == ER_ID)
+    assert(s_info.parameters.kernel_parameters[1].id == VF_ID)
+    assert(s_info.parameters.magnetism_index == [])
     p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters
     s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters
@@ -44 +55 @@
     if p_set & s_set:
         # there is some overlap between the parameter names; tag the
-        # overlapping S parameters with name_S
-        s_list = [(suffix_parameter(par, "_S") if par.id in p_set else par)
-                  for par in s_pars.kernel_parameters]
-        combined_pars = p_pars.kernel_parameters + s_list
+        # overlapping S parameters with name_S.
+        # Skip the first parameter of s, which is effective radius
+        s_list = [(suffix_parameter(par) if par.id in p_set else par)
+                  for par in s_pars.kernel_parameters[1:]]
     else:
-        combined_pars = p_pars.kernel_parameters + s_pars.kernel_parameters
+        # Skip the first parameter of s, which is effective radius
+        s_list = s_pars.kernel_parameters[1:]
+    translate_name = dict((old.id, new.id) for old, new
+                          in zip(s_pars.kernel_parameters[1:], s_list))
+    demo = {}
+    demo.update((k, v) for k, v in p_info.demo.items()
+                if k not in ("background", "scale"))
+    demo.update((translate_name[k], v) for k, v in s_info.demo.items()
+                if k not in ("background", "scale") and not k.startswith(ER_ID))
+    combined_pars = p_pars.kernel_parameters + s_list
     parameters = ParameterTable(combined_pars)
+    parameters.max_pd = p_pars.max_pd + s_pars.max_pd
 
     model_info = ModelInfo()
     model_info.id = '*'.join((p_id, s_id))
-    model_info.name = ' X '.join((p_name, s_name))
+    model_info.name = '*'.join((p_name, s_name))
     model_info.filename = None
     model_info.title = 'Product of %s and %s'%(p_name, s_name)
@@ -67 +88 @@
     #model_info.source = []
     # Iq, Iqxy, form_volume, ER, VR and sesans
+    # Remember the component info blocks so we can build the model
     model_info.composition = ('product', [p_info, s_info])
+    model_info.demo = {}
     return model_info
+
+def suffix_parameter(par, suffix):
+    par = copy(par)
+    par.name = par.name + " S"
+    par.id = par.id + "_S"
+    return par
 
 class ProductModel(KernelModel):
@@ -77 +106 @@
         self.S = S
 
-    def __call__(self, q_vectors):
+    def make_kernel(self, q_vectors):
         # type: (List[np.ndarray]) -> Kernel
         # Note: may be sending the q_vectors to the GPU twice even though they
@@ -104 +133 @@
         self.p_kernel = p_kernel
         self.s_kernel = s_kernel
-
-    def __call__(self, details, weights, values, cutoff):
-        # type: (CallDetails, np.ndarray, np.ndarray, float) -> np.ndarray
-        effect_radius, vol_ratio = call_ER_VR(self.p_kernel.info, vol_pars)
-
-        p_details, s_details = details.parts
-        p_res = self.p_kernel(p_details, weights, values, cutoff)
-        s_res = self.s_kernel(s_details, weights, values, cutoff)
-        #print s_fixed, s_pd, p_fixed, p_pd
-
-        return values[0]*(p_res*s_res) + values[1]
+        self.dtype = p_kernel.dtype
+        self.results = []  # type: List[np.ndarray]
+
+    def __call__(self, call_details, values, cutoff, magnetic):
+        # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray
+        p_info, s_info = self.info.composition[1]
+
+        # if there are magnetic parameters, they will only be on the
+        # form factor P, not the structure factor S.
+        nmagnetic = len(self.info.parameters.magnetism_index)
+        if nmagnetic:
+            spin_index = self.info.parameters.npars + 2
+            magnetism = values[spin_index: spin_index+3+3*nmagnetic]
+        else:
+            magnetism = []
+        nvalues = self.info.parameters.nvalues
+        nweights = call_details.num_weights
+        weights = values[nvalues:nvalues + 2*nweights]
+
+        # 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]
+        spacer = (32 - sum(len(v) for v in p_values)%32)%32
+        p_values.append([0.]*spacer)
+        p_values = np.hstack(p_values).astype(self.p_kernel.dtype)
+
+        # Call ER and VR for P since these are needed for S.
+        p_er, p_vr = calc_er_vr(p_info, p_details, p_values)
+        s_vr = (volfrac/p_vr if p_vr != 0. else volfrac)
+        #print("volfrac:%g p_er:%g p_vr:%g s_vr:%g"%(volfrac,p_er,p_vr,s_vr))
+
+        # Construct the calling parameters for S.
+        # The  effective radius is not in the combined parameter list, so
+        # the number of 'S' parameters is one less than expected.  The
+        # computed effective radius needs to be added into the weights
+        # vector, especially since it is a polydisperse parameter in the
+        # stand-alone structure factor models.  We will added it at the
+        # end so the remaining offsets don't need to change.
+        s_npars = s_info.parameters.npars-1
+        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)
+        v, w = weights[:nweights], weights[nweights:]
+        s_values = [[1., 0., p_er, s_vr],
+                    # er and vf already included, so start one past the
+                    # volfrac parameter, and go two less than the number
+                    # of S parameters.
+                    values[2+p_npars+2:2+p_npars+s_npars-1],
+                    # no magnetism parameters to include for S
+                    # add er into the (value, weights) pairs
+                    v, [p_er], w, [1.0]]
+        spacer = (32 - sum(len(v) for v in s_values)%32)%32
+        s_values.append([0.]*spacer)
+        s_values = np.hstack(s_values).astype(self.s_kernel.dtype)
+
+        # Call the kernels
+        p_result = self.p_kernel(p_details, p_values, cutoff, magnetic)
+        s_result = self.s_kernel(s_details, s_values, cutoff, False)
+
+        #call_details.show(values)
+        #print("values", values)
+        #p_details.show(p_values)
+        #print("=>", p_result)
+        #print("p val", s_values)
+        #s_details.show(s_values)
+        #print("=>", s_result)
+
+        # remember the parts for plotting later
+        self.results = [p_result, 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()
+
+        return values[0]*(p_result*s_result) + values[1]
 
     def release(self):
@@ -121 +223 @@
         self.s_kernel.release()
 
-def call_ER_VR(model_info, pars):
-    """
-    Return effect radius and volume ratio for the model.
-    """
+
+def calc_er_vr(model_info, call_details, values):
+    # type: (ModelInfo, ParameterSet) -> Tuple[float, float]
+
     if model_info.ER is None and model_info.VR is None:
         return 1.0, 1.0
 
-    value, weight = _vol_pars(model_info, pars)
+    nvalues = model_info.parameters.nvalues
+    value = values[nvalues:nvalues + call_details.num_weights]
+    weight = values[nvalues + call_details.num_weights: nvalues + 2*call_details.num_weights]
+    npars = model_info.parameters.npars
+    pairs = [(value[offset:offset+length], weight[offset:offset+length])
+             for p, offset, length
+             in zip(model_info.parameters.call_parameters[2:2+npars],
+                    call_details.offset,
+                    call_details.length)
+             if p.type == 'volume']
+    value, weight = dispersion_mesh(model_info, pairs)
 
     if model_info.ER is not None:
         individual_radii = model_info.ER(*value)
-        effect_radius = np.sum(weight*individual_radii) / np.sum(weight)
+        radius_effective = np.sum(weight*individual_radii) / np.sum(weight)
     else:
-        effect_radius = 1.0
+        radius_effective = 1.0
 
     if model_info.VR is not None:
@@ -142 +254 @@
         volume_ratio = 1.0
 
-    return effect_radius, volume_ratio
-
-def _vol_pars(model_info, pars):
-    # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray]
-    vol_pars = [get_weights(p, pars)
-                for p in model_info.parameters.call_parameters
-                if p.type == 'volume']
-    value, weight = dispersion_mesh(model_info, vol_pars)
-    return value, weight
-
+    return radius_effective, volume_ratio