Changeset 0edb6c1 in sasmodels

Timestamp:

Sep 12, 2017 5:12:03 AM (7 years ago)

Author:

lewis

Branches:

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

Children:

9d1e3e4

Parents:

2ad5d30 (diff), 22e922e (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 'ticket-767' of ​https://github.com/lewisodriscoll/sasmodels into ticket-767

Location:

sasmodels

Files:

• core.py (modified) (4 diffs)
• mixture.py (modified) (7 diffs)
• model_test.py (modified) (2 diffs)
• modelinfo.py (modified) (1 diff)
• product.py (modified) (1 diff)
• sasview_model.py (modified) (5 diffs)
• compare.py (modified) (4 diffs)
• models/spherical_sld.py (modified) (1 diff)

sasmodels/core.py

@@ -10 +10 @@
 
 import os
+import re
 from os.path import basename, dirname, join as joinpath
 from glob import glob
@@ -21 +22 @@
 from . import kernelpy
 from . import kerneldll
+from . import custom
 
 if os.environ.get("SAS_OPENCL", "").lower() == "none":
@@ -30 +32 @@
     except Exception:
         HAVE_OPENCL = False
+
+CUSTOM_MODEL_PATH = os.environ.get('SAS_MODELPATH', "")
+if CUSTOM_MODEL_PATH == "":
+    path = joinpath(os.path.expanduser("~"), ".sasmodels", "custom_models")
+    if not os.path.isdir(path):
+        os.makedirs(path)
+    CUSTOM_MODEL_PATH = path
 
 try:
@@ -125 +134 @@
                        dtype=dtype, platform=platform)
 
-
-def load_model_info(model_name):
+def load_model_info(model_string):
     # type: (str) -> modelinfo.ModelInfo
     """
     Load a model definition given the model name.
 
-    *model_name* is the name of the model, or perhaps a model expression
-    such as sphere*hardsphere or sphere+cylinder.
+    *model_string* is the name of the model, or perhaps a model expression
+    such as sphere*cylinder or sphere+cylinder. Use '@' for a structure
+    factor product, e.g. sphere@hardsphere. Custom models can be specified by
+    prefixing the model name with 'custom.', e.g. 'custom.MyModel+sphere'.
 
     This returns a handle to the module defining the model.  This can be
     used with functions in generate to build the docs or extract model info.
     """
-    parts = model_name.split('+')
-    if len(parts) > 1:
-        model_info_list = [load_model_info(p) for p in parts]
-        return mixture.make_mixture_info(model_info_list)
-
-    parts = model_name.split('*')
-    if len(parts) > 1:
-        if len(parts) > 2:
-            raise ValueError("use P*S to apply structure factor S to model P")
-        P_info, Q_info = [load_model_info(p) for p in parts]
+    if '@' in model_string:
+        parts = model_string.split('@')
+        if len(parts) != 2:
+            raise ValueError("Use P@S to apply a structure factor S to model P")
+        P_info, Q_info = [load_model_info(part) for part in parts]
         return product.make_product_info(P_info, Q_info)
 
-    kernel_module = generate.load_kernel_module(model_name)
-    return modelinfo.make_model_info(kernel_module)
+    product_parts = []
+    addition_parts = []
+
+    addition_parts_names = model_string.split('+')
+    if len(addition_parts_names) >= 2:
+        addition_parts = [load_model_info(part) for part in addition_parts_names]
+    elif len(addition_parts_names) == 1:
+        product_parts_names = model_string.split('*')
+        if len(product_parts_names) >= 2:
+            product_parts = [load_model_info(part) for part in product_parts_names]
+        elif len(product_parts_names) == 1:
+            if "custom." in product_parts_names[0]:
+                # Extract ModelName from "custom.ModelName"
+                pattern = "custom.([A-Za-z0-9_-]+)"
+                result = re.match(pattern, product_parts_names[0])
+                if result is None:
+                    raise ValueError("Model name in invalid format: " + product_parts_names[0])
+                model_name = result.group(1)
+                # Use ModelName to find the path to the custom model file
+                model_path = joinpath(CUSTOM_MODEL_PATH, model_name + ".py")
+                if not os.path.isfile(model_path):
+                    raise ValueError("The model file {} doesn't exist".format(model_path))
+                kernel_module = custom.load_custom_kernel_module(model_path)
+                return modelinfo.make_model_info(kernel_module)
+            # Model is a core model
+            kernel_module = generate.load_kernel_module(product_parts_names[0])
+            return modelinfo.make_model_info(kernel_module)
+
+    model = None
+    if len(product_parts) > 1:
+        model = mixture.make_mixture_info(product_parts, operation='*')
+    if len(addition_parts) > 1:
+        if model is not None:
+            addition_parts.append(model)
+        model = mixture.make_mixture_info(addition_parts, operation='+')
+    return model
 
 

sasmodels/mixture.py

@@ -25 +25 @@
     pass
 
-def make_mixture_info(parts):
+def make_mixture_info(parts, operation='+'):
     # type: (List[ModelInfo]) -> ModelInfo
     """
     Create info block for mixture model.
     """
-    flatten = []
-    for part in parts:
-        if part.composition and part.composition[0] == 'mixture':
-            flatten.extend(part.composition[1])
-        else:
-            flatten.append(part)
-    parts = flatten
-
     # Build new parameter list
     combined_pars = []
-    for k, part in enumerate(parts):
+
+    model_num = 0
+    all_parts = copy(parts)
+    is_flat = False
+    while not is_flat:
+        is_flat = True
+        for part in all_parts:
+            if part.composition and part.composition[0] == 'mixture' and \
+                len(part.composition[1]) > 1:
+                all_parts += part.composition[1]
+                all_parts.remove(part)
+                is_flat = False
+
+    # When creating a mixture model that is a sum of product models (ie (1*2)+(3*4))
+    # the parameters for models 1 & 2 will be prefixed with A & B respectively,
+    # but so will the parameters for models 3 & 4. We need to rename models 3 & 4
+    # so that they are prefixed with C & D to avoid overlap of parameter names.
+    used_prefixes = []
+    for part in parts:
+        i = 0
+        if part.composition and part.composition[0] == 'mixture':
+            npars_list = [info.parameters.npars for info in part.composition[1]]
+            for npars in npars_list:
+                # List of params of one of the constituent models of part
+                submodel_pars = part.parameters.kernel_parameters[i:i+npars]
+                # Prefix of the constituent model
+                prefix = submodel_pars[0].name[0]
+                if prefix not in used_prefixes: # Haven't seen this prefix so far
+                    used_prefixes.append(prefix)
+                    i += npars
+                    continue
+                while prefix in used_prefixes:
+                    # This prefix has been already used, so change it to the
+                    # next letter that hasn't been used
+                    prefix = chr(ord(prefix) + 1)
+                used_prefixes.append(prefix)
+                prefix += "_"
+                # Update the parameters of this constituent model to use the
+                # new prefix
+                for par in submodel_pars:
+                    par.id = prefix + par.id[2:]
+                    par.name = prefix + par.name[2:]
+                    if par.length_control is not None:
+                        par.length_control = prefix + par.length_control[2:]
+                i += npars
+
+    for part in parts:
         # Parameter prefix per model, A_, B_, ...
         # Note that prefix must also be applied to id and length_control
         # to support vector parameters
-        prefix = chr(ord('A')+k) + '_'
-        scale =  Parameter(prefix+'scale', default=1.0,
-                           description="model intensity for " + part.name)
-        combined_pars.append(scale)
+        prefix = ''
+        if not part.composition:
+            # Model isn't a composition model, so it's parameters don't have a
+            # a prefix. Add the next available prefix
+            prefix = chr(ord('A')+len(used_prefixes))
+            used_prefixes.append(prefix)
+            prefix += '_'
+            
+        if operation == '+':
+            # If model is a sum model, each constituent model gets its own scale parameter
+            scale_prefix = prefix
+            if prefix == '' and part.operation == '*':
+                # `part` is a composition product model. Find the prefixes of
+                # it's parameters to form a new prefix for the scale, eg:
+                # a model with A*B*C will have ABC_scale
+                sub_prefixes = []
+                for param in part.parameters.kernel_parameters:
+                    # Prefix of constituent model
+                    sub_prefix = param.id.split('_')[0]
+                    if sub_prefix not in sub_prefixes:
+                        sub_prefixes.append(sub_prefix)
+                # Concatenate sub_prefixes to form prefix for the scale
+                scale_prefix = ''.join(sub_prefixes) + '_'
+            scale =  Parameter(scale_prefix + 'scale', default=1.0,
+                            description="model intensity for " + part.name)
+            combined_pars.append(scale)
         for p in part.parameters.kernel_parameters:
             p = copy(p)
@@ -67 +127 @@
 
     model_info = ModelInfo()
-    model_info.id = '+'.join(part.id for part in parts)
-    model_info.name = ' + '.join(part.name for part in parts)
+    model_info.id = operation.join(part.id for part in parts)
+    model_info.operation = operation
+    model_info.name = '(' + operation.join(part.name for part in parts) + ')'
     model_info.filename = None
     model_info.title = 'Mixture model with ' + model_info.name
@@ -121 +182 @@
         self.kernels = kernels
         self.dtype = self.kernels[0].dtype
+        self.operation = model_info.operation
         self.results = []  # type: List[np.ndarray]
 
@@ -129 +191 @@
         # remember the parts for plotting later
         self.results = []  # type: List[np.ndarray]
-        offset = 2 # skip scale & background
         parts = MixtureParts(self.info, self.kernels, call_details, values)
         for kernel, kernel_details, kernel_values in parts:
             #print("calling kernel", kernel.info.name)
             result = kernel(kernel_details, kernel_values, cutoff, magnetic)
-            #print(kernel.info.name, result)
-            total += result
+            result = np.array(result).astype(kernel.dtype)
+            # print(kernel.info.name, result)
+            if self.operation == '+':
+                total += result
+            elif self.operation == '*':
+                if np.all(total) == 0.0:
+                    total = result
+                else:
+                    total *= result
             self.results.append(result)
 
@@ -176 +244 @@
 
         self.part_num += 1
-        self.par_index += info.parameters.npars + 1
+        self.par_index += info.parameters.npars
+        if self.model_info.operation == '+':
+            self.par_index += 1 # Account for each constituent model's scale param
         self.mag_index += 3 * len(info.parameters.magnetism_index)
 
@@ -187 +257 @@
         # which includes the initial scale and background parameters.
         # We want the index into the weight length/offset for each parameter.
-        # Exclude the initial scale and background, so subtract two, but each
-        # component has its own scale factor which we need to skip when
-        # constructing the details for the kernel, so add one, giving a
-        # net subtract one.
-        index = slice(par_index - 1, par_index - 1 + info.parameters.npars)
+        # Exclude the initial scale and background, so subtract two. If we're
+        # building an addition model, each component has its own scale factor
+        # which we need to skip when constructing the details for the kernel, so
+        # add one, giving a net subtract one.
+        diff = 1 if self.model_info.operation == '+' else 2
+        index = slice(par_index - diff, par_index - diff + info.parameters.npars)
         length = full.length[index]
         offset = full.offset[index]
@@ -201 +272 @@
     def _part_values(self, info, par_index, mag_index):
         # type: (ModelInfo, int, int) -> np.ndarray
-        #print(info.name, par_index, self.values[par_index:par_index + info.parameters.npars + 1])
-        scale = self.values[par_index]
-        pars = self.values[par_index + 1:par_index + info.parameters.npars + 1]
+        # Set each constituent model's scale to 1 if this is a multiplication model
+        scale = self.values[par_index] if self.model_info.operation == '+' else 1.0
+        diff = 1 if self.model_info.operation == '+' else 0 # Skip scale if addition model
+        pars = self.values[par_index + diff:par_index + info.parameters.npars + diff]
         nmagnetic = len(info.parameters.magnetism_index)
         if nmagnetic:

sasmodels/model_test.py

@@ -201 +201 @@
                 ({}, 'VR', None),
                 ]
-
-            tests = smoke_tests + self.info.tests
+            tests = smoke_tests
+            if self.info.tests is not None:
+                tests += self.info.tests
             try:
                 model = build_model(self.info, dtype=self.dtype,
@@ -371 +372 @@
         stream.writeln(traceback.format_exc())
         return
-
     # Run the test suite
     suite.run(result)

sasmodels/modelinfo.py

@@ -727 +727 @@
     models when the model is first called, not when the model is loaded.
     """
+    if hasattr(kernel_module, "model_info"):
+        # Custom sum/multi models
+        return kernel_module.model_info
     info = ModelInfo()
     #print("make parameter table", kernel_module.parameters)

sasmodels/product.py

@@ -81 +81 @@
 
     model_info = ModelInfo()
-    model_info.id = '*'.join((p_id, s_id))
-    model_info.name = '*'.join((p_name, s_name))
+    model_info.id = '@'.join((p_id, s_id))
+    model_info.name = '@'.join((p_name, s_name))
     model_info.filename = None
     model_info.title = 'Product of %s and %s'%(p_name, s_name)

sasmodels/sasview_model.py

@@ -120 +120 @@
     else:
         model_info = modelinfo.make_model_info(kernel_module)
-        model = _make_model_from_info(model_info)
+        model = make_model_from_info(model_info)
     model.timestamp = getmtime(path)
 
@@ -142 +142 @@
 
 
+def make_model_from_info(model_info):
+    # type: (ModelInfo) -> SasviewModelType
+    """
+    Convert *model_info* into a SasView model wrapper.
+    """
+    def __init__(self, multiplicity=None):
+        SasviewModel.__init__(self, multiplicity=multiplicity)
+    attrs = _generate_model_attributes(model_info)
+    attrs['__init__'] = __init__
+    attrs['filename'] = model_info.filename
+    ConstructedModel = type(model_info.name, (SasviewModel,), attrs) # type: SasviewModelType
+    return ConstructedModel
+
+
 def _make_standard_model(name):
     # type: (str) -> SasviewModelType
@@ -153 +167 @@
     kernel_module = generate.load_kernel_module(name)
     model_info = modelinfo.make_model_info(kernel_module)
-    return _make_model_from_info(model_info)
-
-
+    return make_model_from_info(model_info)
+
+    
 def _register_old_models():
     # type: () -> None
@@ -187 +201 @@
     model_info = product.make_product_info(form_factor._model_info,
                                            structure_factor._model_info)
-    ConstructedModel = _make_model_from_info(model_info)
+    ConstructedModel = make_model_from_info(model_info)
     return ConstructedModel()
 
-def _make_model_from_info(model_info):
-    # type: (ModelInfo) -> SasviewModelType
-    """
-    Convert *model_info* into a SasView model wrapper.
-    """
-    def __init__(self, multiplicity=None):
-        SasviewModel.__init__(self, multiplicity=multiplicity)
-    attrs = _generate_model_attributes(model_info)
-    attrs['__init__'] = __init__
-    attrs['filename'] = model_info.filename
-    ConstructedModel = type(model_info.name, (SasviewModel,), attrs) # type: SasviewModelType
-    return ConstructedModel
 
 def _generate_model_attributes(model_info):
@@ -603 +605 @@
         if hasattr(self._model_info, "composition") \
            and self._model_info.composition is not None:
-            p_model = _make_model_from_info(self._model_info.composition[1][0])()
-            s_model = _make_model_from_info(self._model_info.composition[1][1])()
+            p_model = make_model_from_info(self._model_info.composition[1][0])()
+            s_model = make_model_from_info(self._model_info.composition[1][1])()
         return p_model, s_model
 

sasmodels/compare.py

@@ -88 +88 @@
     -magnetic/-nonmagnetic* suppress magnetism
     -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh
+    -neval=1 sets the number of evals for more accurate timing
 
     === precision options ===
@@ -1013 +1014 @@
     'magnetic', 'nonmagnetic',
     'accuracy=',
+    'neval=',  # for timing...
 
     # Precision options
@@ -1018 +1020 @@
     'half', 'fast', 'single', 'double', 'single!', 'double!', 'quad!',
     'sasview',  # TODO: remove sasview 3.x support
-    'timing=',
 
     # Output options
@@ -1068 +1069 @@
 
     # Plug in values given in demo
-    if use_demo:
+    if use_demo and model_info.demo:
         pars.update(model_info.demo)
     return pars

sasmodels/models/spherical_sld.py

@@ -18 +18 @@
 Interface shapes are as follows::
 
-    0: erf(|nu|*z)
-    1: Rpow(z^|nu|)
-    2: Lpow(z^|nu|)
-    3: Rexp(-|nu|z)
-    4: Lexp(-|nu|z)
+    0: erf($\nu z$)
+    1: Rpow($z^\nu$)
+    2: Lpow($z^\nu$)
+    3: Rexp($-\nu z$)
+    4: Lexp($-\nu z$)
 
 Definition