source: sasmodels/sasmodels/sasview_model.py @ 5024a56

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 5024a56 was 5024a56, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

Make sure that the label radius_effective_mode is used throughout

  • Property mode set to 100644
File size: 33.1 KB
RevLine 
[87985ca]1"""
2Sasview model constructor.
3
4Given a module defining an OpenCL kernel such as sasmodels.models.cylinder,
5create a sasview model class to run that kernel as follows::
6
[92d38285]7    from sasmodels.sasview_model import load_custom_model
8    CylinderModel = load_custom_model('sasmodels/models/cylinder.py')
[87985ca]9"""
[4d76711]10from __future__ import print_function
[87985ca]11
[ce27e21]12import math
13from copy import deepcopy
[2622b3f]14import collections
[4d76711]15import traceback
16import logging
[724257c]17from os.path import basename, splitext, abspath, getmtime
[9f8ade1]18try:
19    import _thread as thread
20except ImportError:
21    import thread
[ce27e21]22
[7ae2b7f]23import numpy as np  # type: ignore
[ce27e21]24
[aa4946b]25from . import core
[4d76711]26from . import custom
[a80e64c]27from . import product
[72a081d]28from . import generate
[fb5914f]29from . import weights
[6d6508e]30from . import modelinfo
[bde38b5]31from .details import make_kernel_args, dispersion_mesh
[ff7119b]32
[2d81cfe]33# pylint: disable=unused-import
[fa5fd8d]34try:
[2d81cfe]35    from typing import (Dict, Mapping, Any, Sequence, Tuple, NamedTuple,
36                        List, Optional, Union, Callable)
[fa5fd8d]37    from .modelinfo import ModelInfo, Parameter
38    from .kernel import KernelModel
39    MultiplicityInfoType = NamedTuple(
[a9bc435]40        'MultiplicityInfo',
[fa5fd8d]41        [("number", int), ("control", str), ("choices", List[str]),
42         ("x_axis_label", str)])
[60f03de]43    SasviewModelType = Callable[[int], "SasviewModel"]
[fa5fd8d]44except ImportError:
45    pass
[2d81cfe]46# pylint: enable=unused-import
[fa5fd8d]47
[724257c]48logger = logging.getLogger(__name__)
49
[a38b065]50calculation_lock = thread.allocate_lock()
51
[724257c]52#: True if pre-existing plugins, with the old names and parameters, should
53#: continue to be supported.
[c95dfc63]54SUPPORT_OLD_STYLE_PLUGINS = True
55
[fa5fd8d]56# TODO: separate x_axis_label from multiplicity info
57MultiplicityInfo = collections.namedtuple(
58    'MultiplicityInfo',
59    ["number", "control", "choices", "x_axis_label"],
60)
61
[724257c]62#: set of defined models (standard and custom)
63MODELS = {}  # type: Dict[str, SasviewModelType]
[839fd68]64# TODO: remove unused MODEL_BY_PATH cache once sasview no longer references it
[724257c]65#: custom model {path: model} mapping so we can check timestamps
66MODEL_BY_PATH = {}  # type: Dict[str, SasviewModelType]
[d321747]67#: Track modules that we have loaded so we can determine whether the model
68#: has changed since we last reloaded.
69_CACHED_MODULE = {}  # type: Dict[str, "module"]
[724257c]70
[92d38285]71def find_model(modelname):
[b32dafd]72    # type: (str) -> SasviewModelType
73    """
74    Find a model by name.  If the model name ends in py, try loading it from
75    custom models, otherwise look for it in the list of builtin models.
76    """
[92d38285]77    # TODO: used by sum/product model to load an existing model
78    # TODO: doesn't handle custom models properly
79    if modelname.endswith('.py'):
80        return load_custom_model(modelname)
81    elif modelname in MODELS:
82        return MODELS[modelname]
83    else:
84        raise ValueError("unknown model %r"%modelname)
85
[56b2687]86
[fa5fd8d]87# TODO: figure out how to say that the return type is a subclass
[4d76711]88def load_standard_models():
[60f03de]89    # type: () -> List[SasviewModelType]
[4d76711]90    """
91    Load and return the list of predefined models.
92
93    If there is an error loading a model, then a traceback is logged and the
94    model is not returned.
95    """
96    for name in core.list_models():
97        try:
[92d38285]98            MODELS[name] = _make_standard_model(name)
[ee8f734]99        except Exception:
[724257c]100            logger.error(traceback.format_exc())
[c95dfc63]101    if SUPPORT_OLD_STYLE_PLUGINS:
102        _register_old_models()
103
[724257c]104    return list(MODELS.values())
[de97440]105
[4d76711]106
107def load_custom_model(path):
[60f03de]108    # type: (str) -> SasviewModelType
[4d76711]109    """
110    Load a custom model given the model path.
[ff7119b]111    """
[724257c]112    #logger.info("Loading model %s", path)
[d321747]113
114    # Load the kernel module.  This may already be cached by the loader, so
115    # only requires checking the timestamps of the dependents.
[4d76711]116    kernel_module = custom.load_custom_kernel_module(path)
[d321747]117
118    # Check if the module has changed since we last looked.
119    reloaded = kernel_module != _CACHED_MODULE.get(path, None)
120    _CACHED_MODULE[path] = kernel_module
121
122    # Turn the module into a model.  We need to do this in even if the
123    # model has already been loaded so that we can determine the model
124    # name and retrieve it from the MODELS cache.
125    model = getattr(kernel_module, 'Model', None)
126    if model is not None:
[9457498]127        # Old style models do not set the name in the class attributes, so
128        # set it here; this name will be overridden when the object is created
129        # with an instance variable that has the same value.
130        if model.name == "":
131            model.name = splitext(basename(path))[0]
[20a70bc]132        if not hasattr(model, 'filename'):
[724257c]133            model.filename = abspath(kernel_module.__file__).replace('.pyc', '.py')
[e4bf271]134        if not hasattr(model, 'id'):
135            model.id = splitext(basename(model.filename))[0]
[724257c]136    else:
[56b2687]137        model_info = modelinfo.make_model_info(kernel_module)
[bcdd6c9]138        model = make_model_from_info(model_info)
[ed10b57]139
[2f2c70c]140    # If a model name already exists and we are loading a different model,
141    # use the model file name as the model name.
142    if model.name in MODELS and not model.filename == MODELS[model.name].filename:
143        _previous_name = model.name
144        model.name = model.id
[bf8c271]145
[2f2c70c]146        # If the new model name is still in the model list (for instance,
147        # if we put a cylinder.py in our plug-in directory), then append
148        # an identifier.
149        if model.name in MODELS and not model.filename == MODELS[model.name].filename:
150            model.name = model.id + '_user'
[724257c]151        logger.info("Model %s already exists: using %s [%s]",
152                    _previous_name, model.name, model.filename)
[ed10b57]153
[d321747]154    # Only update the model if the module has changed
155    if reloaded or model.name not in MODELS:
156        MODELS[model.name] = model
157
158    return MODELS[model.name]
[4d76711]159
[87985ca]160
[bcdd6c9]161def make_model_from_info(model_info):
162    # type: (ModelInfo) -> SasviewModelType
163    """
164    Convert *model_info* into a SasView model wrapper.
165    """
166    def __init__(self, multiplicity=None):
167        SasviewModel.__init__(self, multiplicity=multiplicity)
168    attrs = _generate_model_attributes(model_info)
169    attrs['__init__'] = __init__
170    attrs['filename'] = model_info.filename
171    ConstructedModel = type(model_info.name, (SasviewModel,), attrs) # type: SasviewModelType
172    return ConstructedModel
173
174
[4d76711]175def _make_standard_model(name):
[60f03de]176    # type: (str) -> SasviewModelType
[ff7119b]177    """
[4d76711]178    Load the sasview model defined by *name*.
[72a081d]179
[4d76711]180    *name* can be a standard model name or a path to a custom model.
[87985ca]181
[4d76711]182    Returns a class that can be used directly as a sasview model.
[ff7119b]183    """
[4d76711]184    kernel_module = generate.load_kernel_module(name)
[fa5fd8d]185    model_info = modelinfo.make_model_info(kernel_module)
[bcdd6c9]186    return make_model_from_info(model_info)
[72a081d]187
188
[724257c]189def _register_old_models():
190    # type: () -> None
191    """
192    Place the new models into sasview under the old names.
193
194    Monkey patch sas.sascalc.fit as sas.models so that sas.models.pluginmodel
195    is available to the plugin modules.
196    """
197    import sys
198    import sas   # needed in order to set sas.models
199    import sas.sascalc.fit
200    sys.modules['sas.models'] = sas.sascalc.fit
201    sas.models = sas.sascalc.fit
202    import sas.models
203    from sasmodels.conversion_table import CONVERSION_TABLE
[e65c3ba]204
[724257c]205    for new_name, conversion in CONVERSION_TABLE.get((3, 1, 2), {}).items():
206        # CoreShellEllipsoidModel => core_shell_ellipsoid:1
207        new_name = new_name.split(':')[0]
208        old_name = conversion[0] if len(conversion) < 3 else conversion[2]
209        module_attrs = {old_name: find_model(new_name)}
210        ConstructedModule = type(old_name, (), module_attrs)
211        old_path = 'sas.models.' + old_name
212        setattr(sas.models, old_path, ConstructedModule)
213        sys.modules[old_path] = ConstructedModule
214
215
[a80e64c]216def MultiplicationModel(form_factor, structure_factor):
217    # type: ("SasviewModel", "SasviewModel") -> "SasviewModel"
[e65c3ba]218    """
219    Returns a constructed product model from form_factor and structure_factor.
220    """
[a80e64c]221    model_info = product.make_product_info(form_factor._model_info,
222                                           structure_factor._model_info)
[bcdd6c9]223    ConstructedModel = make_model_from_info(model_info)
[a06af5d]224    return ConstructedModel(form_factor.multiplicity)
[a80e64c]225
[ce27e21]226
[fa5fd8d]227def _generate_model_attributes(model_info):
228    # type: (ModelInfo) -> Dict[str, Any]
229    """
230    Generate the class attributes for the model.
231
232    This should include all the information necessary to query the model
233    details so that you do not need to instantiate a model to query it.
234
235    All the attributes should be immutable to avoid accidents.
236    """
237
238    # TODO: allow model to override axis labels input/output name/unit
239
[a18c5b3]240    # Process multiplicity
[fa5fd8d]241    non_fittable = []  # type: List[str]
[04045f4]242    xlabel = model_info.profile_axes[0] if model_info.profile is not None else ""
243    variants = MultiplicityInfo(0, "", [], xlabel)
[a18c5b3]244    for p in model_info.parameters.kernel_parameters:
[04045f4]245        if p.name == model_info.control:
[fa5fd8d]246            non_fittable.append(p.name)
[04045f4]247            variants = MultiplicityInfo(
[ce176ca]248                len(p.choices) if p.choices else int(p.limits[1]),
249                p.name, p.choices, xlabel
[fa5fd8d]250            )
251            break
252
[50ec515]253    # Only a single drop-down list parameter available
254    fun_list = []
255    for p in model_info.parameters.kernel_parameters:
256        if p.choices:
257            fun_list = p.choices
258            if p.length > 1:
259                non_fittable.extend(p.id+str(k) for k in range(1, p.length+1))
260            break
261
[a18c5b3]262    # Organize parameter sets
[fa5fd8d]263    orientation_params = []
264    magnetic_params = []
265    fixed = []
[85fe7f8]266    for p in model_info.parameters.user_parameters({}, is2d=True):
[fa5fd8d]267        if p.type == 'orientation':
268            orientation_params.append(p.name)
269            orientation_params.append(p.name+".width")
270            fixed.append(p.name+".width")
[32e3c9b]271        elif p.type == 'magnetic':
[fa5fd8d]272            orientation_params.append(p.name)
273            magnetic_params.append(p.name)
274            fixed.append(p.name+".width")
[a18c5b3]275
[32e3c9b]276
[a18c5b3]277    # Build class dictionary
278    attrs = {}  # type: Dict[str, Any]
279    attrs['_model_info'] = model_info
280    attrs['name'] = model_info.name
281    attrs['id'] = model_info.id
282    attrs['description'] = model_info.description
283    attrs['category'] = model_info.category
284    attrs['is_structure_factor'] = model_info.structure_factor
[6e7ba14]285    attrs['is_form_factor'] = model_info.effective_radius_type is not None
[a18c5b3]286    attrs['is_multiplicity_model'] = variants[0] > 1
287    attrs['multiplicity_info'] = variants
[fa5fd8d]288    attrs['orientation_params'] = tuple(orientation_params)
289    attrs['magnetic_params'] = tuple(magnetic_params)
290    attrs['fixed'] = tuple(fixed)
291    attrs['non_fittable'] = tuple(non_fittable)
[50ec515]292    attrs['fun_list'] = tuple(fun_list)
[fa5fd8d]293
294    return attrs
[4d76711]295
[ce27e21]296class SasviewModel(object):
297    """
298    Sasview wrapper for opencl/ctypes model.
299    """
[fa5fd8d]300    # Model parameters for the specific model are set in the class constructor
301    # via the _generate_model_attributes function, which subclasses
302    # SasviewModel.  They are included here for typing and documentation
303    # purposes.
304    _model = None       # type: KernelModel
305    _model_info = None  # type: ModelInfo
306    #: load/save name for the model
307    id = None           # type: str
308    #: display name for the model
309    name = None         # type: str
310    #: short model description
311    description = None  # type: str
312    #: default model category
313    category = None     # type: str
314
315    #: names of the orientation parameters in the order they appear
[724257c]316    orientation_params = None # type: List[str]
[fa5fd8d]317    #: names of the magnetic parameters in the order they appear
[724257c]318    magnetic_params = None    # type: List[str]
[fa5fd8d]319    #: names of the fittable parameters
[724257c]320    fixed = None              # type: List[str]
[fa5fd8d]321    # TODO: the attribute fixed is ill-named
322
323    # Axis labels
324    input_name = "Q"
325    input_unit = "A^{-1}"
326    output_name = "Intensity"
327    output_unit = "cm^{-1}"
328
329    #: default cutoff for polydispersity
330    cutoff = 1e-5
331
332    # Note: Use non-mutable values for class attributes to avoid errors
333    #: parameters that are not fitted
334    non_fittable = ()        # type: Sequence[str]
335
336    #: True if model should appear as a structure factor
337    is_structure_factor = False
338    #: True if model should appear as a form factor
339    is_form_factor = False
340    #: True if model has multiplicity
341    is_multiplicity_model = False
[1f35235]342    #: Multiplicity information
[fa5fd8d]343    multiplicity_info = None # type: MultiplicityInfoType
344
345    # Per-instance variables
346    #: parameter {name: value} mapping
347    params = None      # type: Dict[str, float]
348    #: values for dispersion width, npts, nsigmas and type
349    dispersion = None  # type: Dict[str, Any]
350    #: units and limits for each parameter
[60f03de]351    details = None     # type: Dict[str, Sequence[Any]]
352    #                  # actual type is Dict[str, List[str, float, float]]
[04dc697]353    #: multiplicity value, or None if no multiplicity on the model
[fa5fd8d]354    multiplicity = None     # type: Optional[int]
[04dc697]355    #: memory for polydispersity array if using ArrayDispersion (used by sasview).
356    _persistency_dict = None # type: Dict[str, Tuple[np.ndarray, np.ndarray]]
[fa5fd8d]357
358    def __init__(self, multiplicity=None):
[04dc697]359        # type: (Optional[int]) -> None
[2622b3f]360
[04045f4]361        # TODO: _persistency_dict to persistency_dict throughout sasview
362        # TODO: refactor multiplicity to encompass variants
363        # TODO: dispersion should be a class
[fa5fd8d]364        # TODO: refactor multiplicity info
365        # TODO: separate profile view from multiplicity
366        # The button label, x and y axis labels and scale need to be under
367        # the control of the model, not the fit page.  Maximum flexibility,
368        # the fit page would supply the canvas and the profile could plot
369        # how it wants, but this assumes matplotlib.  Next level is that
370        # we provide some sort of data description including title, labels
371        # and lines to plot.
372
[1f35235]373        # Get the list of hidden parameters given the multiplicity
[04045f4]374        # Don't include multiplicity in the list of parameters
[fa5fd8d]375        self.multiplicity = multiplicity
[04045f4]376        if multiplicity is not None:
377            hidden = self._model_info.get_hidden_parameters(multiplicity)
378            hidden |= set([self.multiplicity_info.control])
379        else:
380            hidden = set()
[8f93522]381        if self._model_info.structure_factor:
382            hidden.add('scale')
383            hidden.add('background')
384            self._model_info.parameters.defaults['background'] = 0.
[04045f4]385
[bd547d0]386        # Update the parameter lists to exclude any hidden parameters
387        self.magnetic_params = tuple(pname for pname in self.magnetic_params
388                                     if pname not in hidden)
389        self.orientation_params = tuple(pname for pname in self.orientation_params
390                                        if pname not in hidden)
391
[04dc697]392        self._persistency_dict = {}
[fa5fd8d]393        self.params = collections.OrderedDict()
[b3a85cd]394        self.dispersion = collections.OrderedDict()
[fa5fd8d]395        self.details = {}
[8977226]396        for p in self._model_info.parameters.user_parameters({}, is2d=True):
[04045f4]397            if p.name in hidden:
[fa5fd8d]398                continue
[fcd7bbd]399            self.params[p.name] = p.default
[c952e59]400            self.details[p.id] = [p.units, p.limits[0], p.limits[1]]
[fb5914f]401            if p.polydisperse:
[fa5fd8d]402                self.details[p.id+".width"] = [
403                    "", 0.0, 1.0 if p.relative_pd else np.inf
404                ]
[fb5914f]405                self.dispersion[p.name] = {
406                    'width': 0,
407                    'npts': 35,
408                    'nsigmas': 3,
409                    'type': 'gaussian',
410                }
[ce27e21]411
[de97440]412    def __get_state__(self):
[fa5fd8d]413        # type: () -> Dict[str, Any]
[de97440]414        state = self.__dict__.copy()
[4d76711]415        state.pop('_model')
[de97440]416        # May need to reload model info on set state since it has pointers
417        # to python implementations of Iq, etc.
418        #state.pop('_model_info')
419        return state
420
421    def __set_state__(self, state):
[fa5fd8d]422        # type: (Dict[str, Any]) -> None
[de97440]423        self.__dict__ = state
[fb5914f]424        self._model = None
[de97440]425
[ce27e21]426    def __str__(self):
[fa5fd8d]427        # type: () -> str
[ce27e21]428        """
429        :return: string representation
430        """
431        return self.name
432
433    def is_fittable(self, par_name):
[fa5fd8d]434        # type: (str) -> bool
[ce27e21]435        """
436        Check if a given parameter is fittable or not
437
438        :param par_name: the parameter name to check
439        """
[e758662]440        return par_name in self.fixed
[ce27e21]441        #For the future
442        #return self.params[str(par_name)].is_fittable()
443
444
445    def getProfile(self):
[fa5fd8d]446        # type: () -> (np.ndarray, np.ndarray)
[ce27e21]447        """
448        Get SLD profile
449
450        : return: (z, beta) where z is a list of depth of the transition points
451                beta is a list of the corresponding SLD values
452        """
[745b7bb]453        args = {} # type: Dict[str, Any]
[fa5fd8d]454        for p in self._model_info.parameters.kernel_parameters:
455            if p.id == self.multiplicity_info.control:
[745b7bb]456                value = float(self.multiplicity)
[fa5fd8d]457            elif p.length == 1:
[745b7bb]458                value = self.params.get(p.id, np.NaN)
[fa5fd8d]459            else:
[745b7bb]460                value = np.array([self.params.get(p.id+str(k), np.NaN)
[b32dafd]461                                  for k in range(1, p.length+1)])
[745b7bb]462            args[p.id] = value
463
[e7fe459]464        x, y = self._model_info.profile(**args)
465        return x, 1e-6*y
[ce27e21]466
467    def setParam(self, name, value):
[fa5fd8d]468        # type: (str, float) -> None
[ce27e21]469        """
470        Set the value of a model parameter
471
472        :param name: name of the parameter
473        :param value: value of the parameter
474
475        """
476        # Look for dispersion parameters
477        toks = name.split('.')
[de0c4ba]478        if len(toks) == 2:
[ce27e21]479            for item in self.dispersion.keys():
[e758662]480                if item == toks[0]:
[ce27e21]481                    for par in self.dispersion[item]:
[e758662]482                        if par == toks[1]:
[ce27e21]483                            self.dispersion[item][par] = value
484                            return
485        else:
486            # Look for standard parameter
487            for item in self.params.keys():
[e758662]488                if item == name:
[ce27e21]489                    self.params[item] = value
490                    return
491
[63b32bb]492        raise ValueError("Model does not contain parameter %s" % name)
[ce27e21]493
494    def getParam(self, name):
[fa5fd8d]495        # type: (str) -> float
[ce27e21]496        """
497        Set the value of a model parameter
498
499        :param name: name of the parameter
500
501        """
502        # Look for dispersion parameters
503        toks = name.split('.')
[de0c4ba]504        if len(toks) == 2:
[ce27e21]505            for item in self.dispersion.keys():
[e758662]506                if item == toks[0]:
[ce27e21]507                    for par in self.dispersion[item]:
[e758662]508                        if par == toks[1]:
[ce27e21]509                            return self.dispersion[item][par]
510        else:
511            # Look for standard parameter
512            for item in self.params.keys():
[e758662]513                if item == name:
[ce27e21]514                    return self.params[item]
515
[63b32bb]516        raise ValueError("Model does not contain parameter %s" % name)
[ce27e21]517
518    def getParamList(self):
[04dc697]519        # type: () -> Sequence[str]
[ce27e21]520        """
521        Return a list of all available parameters for the model
522        """
[04dc697]523        param_list = list(self.params.keys())
[ce27e21]524        # WARNING: Extending the list with the dispersion parameters
[de0c4ba]525        param_list.extend(self.getDispParamList())
526        return param_list
[ce27e21]527
528    def getDispParamList(self):
[04dc697]529        # type: () -> Sequence[str]
[ce27e21]530        """
[fb5914f]531        Return a list of polydispersity parameters for the model
[ce27e21]532        """
[1780d59]533        # TODO: fix test so that parameter order doesn't matter
[3bcb88c]534        ret = ['%s.%s' % (p_name, ext)
535               for p_name in self.dispersion.keys()
536               for ext in ('npts', 'nsigmas', 'width')]
[9404dd3]537        #print(ret)
[1780d59]538        return ret
[ce27e21]539
540    def clone(self):
[04dc697]541        # type: () -> "SasviewModel"
[ce27e21]542        """ Return a identical copy of self """
543        return deepcopy(self)
544
545    def run(self, x=0.0):
[fa5fd8d]546        # type: (Union[float, (float, float), List[float]]) -> float
[ce27e21]547        """
548        Evaluate the model
549
550        :param x: input q, or [q,phi]
551
552        :return: scattering function P(q)
553
554        **DEPRECATED**: use calculate_Iq instead
555        """
[de0c4ba]556        if isinstance(x, (list, tuple)):
[3c56da87]557            # pylint: disable=unpacking-non-sequence
[ce27e21]558            q, phi = x
[84f2962]559            result, _ = self.calculate_Iq([q*math.cos(phi)], [q*math.sin(phi)])
560            return result[0]
[ce27e21]561        else:
[84f2962]562            result, _ = self.calculate_Iq([x])
563            return result[0]
[ce27e21]564
565
566    def runXY(self, x=0.0):
[fa5fd8d]567        # type: (Union[float, (float, float), List[float]]) -> float
[ce27e21]568        """
569        Evaluate the model in cartesian coordinates
570
571        :param x: input q, or [qx, qy]
572
573        :return: scattering function P(q)
574
575        **DEPRECATED**: use calculate_Iq instead
576        """
[de0c4ba]577        if isinstance(x, (list, tuple)):
[84f2962]578            result, _ = self.calculate_Iq([x[0]], [x[1]])
579            return result[0]
[ce27e21]580        else:
[84f2962]581            result, _ = self.calculate_Iq([x])
582            return result[0]
[ce27e21]583
584    def evalDistribution(self, qdist):
[04dc697]585        # type: (Union[np.ndarray, Tuple[np.ndarray, np.ndarray], List[np.ndarray]]) -> np.ndarray
[d138d43]586        r"""
[ce27e21]587        Evaluate a distribution of q-values.
588
[d138d43]589        :param qdist: array of q or a list of arrays [qx,qy]
[ce27e21]590
[d138d43]591        * For 1D, a numpy array is expected as input
[ce27e21]592
[d138d43]593        ::
[ce27e21]594
[d138d43]595            evalDistribution(q)
[ce27e21]596
[d138d43]597          where *q* is a numpy array.
[ce27e21]598
[d138d43]599        * For 2D, a list of *[qx,qy]* is expected with 1D arrays as input
[ce27e21]600
[d138d43]601        ::
[ce27e21]602
[d138d43]603              qx = [ qx[0], qx[1], qx[2], ....]
604              qy = [ qy[0], qy[1], qy[2], ....]
[ce27e21]605
[d138d43]606        If the model is 1D only, then
[ce27e21]607
[d138d43]608        .. math::
[ce27e21]609
[d138d43]610            q = \sqrt{q_x^2+q_y^2}
[ce27e21]611
612        """
[de0c4ba]613        if isinstance(qdist, (list, tuple)):
[ce27e21]614            # Check whether we have a list of ndarrays [qx,qy]
615            qx, qy = qdist
[84f2962]616            result, _ = self.calculate_Iq(qx, qy)
617            return result
[ce27e21]618
619        elif isinstance(qdist, np.ndarray):
620            # We have a simple 1D distribution of q-values
[84f2962]621            result, _ = self.calculate_Iq(qdist)
622            return result
[ce27e21]623
624        else:
[3c56da87]625            raise TypeError("evalDistribution expects q or [qx, qy], not %r"
626                            % type(qdist))
[ce27e21]627
[9dcb21d]628    def calc_composition_models(self, qx):
[64614ad]629        """
[9dcb21d]630        returns parts of the composition model or None if not a composition
631        model.
[64614ad]632        """
[946c8d27]633        # TODO: have calculate_Iq return the intermediates.
634        #
635        # The current interface causes calculate_Iq() to be called twice,
636        # once to get the combined result and again to get the intermediate
637        # results.  This is necessary for now.
638        # Long term, the solution is to change the interface to calculate_Iq
639        # so that it returns a results object containing all the bits:
[9644b5a]640        #     the A, B, C, ... of the composition model (and any subcomponents?)
[d32de68]641        #     the P and S of the product model
[946c8d27]642        #     the combined model before resolution smearing,
643        #     the sasmodel before sesans conversion,
644        #     the oriented 2D model used to fit oriented usans data,
645        #     the final I(q),
646        #     ...
[9644b5a]647        #
[946c8d27]648        # Have the model calculator add all of these blindly to the data
649        # tree, and update the graphs which contain them.  The fitter
650        # needs to be updated to use the I(q) value only, ignoring the rest.
651        #
652        # The simple fix of returning the existing intermediate results
653        # will not work for a couple of reasons: (1) another thread may
654        # sneak in to compute its own results before calc_composition_models
655        # is called, and (2) calculate_Iq is currently called three times:
656        # once with q, once with q values before qmin and once with q values
657        # after q max.  Both of these should be addressed before
658        # replacing this code.
[9644b5a]659        composition = self._model_info.composition
660        if composition and composition[0] == 'product': # only P*S for now
661            with calculation_lock:
[84f2962]662                _, lazy_results = self._calculate_Iq(qx)
663                # for compatibility with sasview 4.x
664                results = lazy_results()
665                pq_data = results.get("P(Q)")
666                sq_data = results.get("S(Q)")
667                return pq_data, sq_data
[9644b5a]668        else:
669            return None
[bf8c271]670
[84f2962]671    def calculate_Iq(self,
672                     qx,     # type: Sequence[float]
673                     qy=None # type: Optional[Sequence[float]]
674                     ):
675        # type: (...) -> Tuple[np.ndarray, Callable[[], collections.OrderedDict[str, np.ndarray]]]
[ff7119b]676        """
677        Calculate Iq for one set of q with the current parameters.
678
679        If the model is 1D, use *q*.  If 2D, use *qx*, *qy*.
680
681        This should NOT be used for fitting since it copies the *q* vectors
682        to the card for each evaluation.
[84f2962]683
684        The returned tuple contains the scattering intensity followed by a
685        callable which returns a dictionary of intermediate data from
686        ProductKernel.
[ff7119b]687        """
[a38b065]688        ## uncomment the following when trying to debug the uncoordinated calls
689        ## to calculate_Iq
690        #if calculation_lock.locked():
[724257c]691        #    logger.info("calculation waiting for another thread to complete")
692        #    logger.info("\n".join(traceback.format_stack()))
[a38b065]693
694        with calculation_lock:
695            return self._calculate_Iq(qx, qy)
696
[39a06c9]697    def _calculate_Iq(self, qx, qy=None, Fq=False, effective_radius_type=1):
[fb5914f]698        if self._model is None:
[d2bb604]699            self._model = core.build_model(self._model_info)
[fa5fd8d]700        if qy is not None:
701            q_vectors = [np.asarray(qx), np.asarray(qy)]
702        else:
703            q_vectors = [np.asarray(qx)]
[a738209]704        calculator = self._model.make_kernel(q_vectors)
[6a0d6aa]705        parameters = self._model_info.parameters
706        pairs = [self._get_weights(p) for p in parameters.call_parameters]
[9c1a59c]707        #weights.plot_weights(self._model_info, pairs)
[bde38b5]708        call_details, values, is_magnetic = make_kernel_args(calculator, pairs)
[4edec6f]709        #call_details.show()
[05df1de]710        #print("================ parameters ==================")
711        #for p, v in zip(parameters.call_parameters, pairs): print(p.name, v[0])
[ce99754]712        #for k, p in enumerate(self._model_info.parameters.call_parameters):
713        #    print(k, p.name, *pairs[k])
[4edec6f]714        #print("params", self.params)
715        #print("values", values)
716        #print("is_mag", is_magnetic)
[39a06c9]717        if Fq:
718            result = calculator.Fq(call_details, values, cutoff=self.cutoff,
719                                   magnetic=is_magnetic,
720                                   effective_radius_type=effective_radius_type)
[6a0d6aa]721        result = calculator(call_details, values, cutoff=self.cutoff,
[9eb3632]722                            magnetic=is_magnetic)
[84f2962]723        lazy_results = getattr(calculator, 'results',
724                               lambda: collections.OrderedDict())
[ce99754]725        #print("result", result)
[84f2962]726
[a738209]727        calculator.release()
[d533590]728        #self._model.release()
[84f2962]729
730        return result, lazy_results
[ce27e21]731
[39a06c9]732
733    def calculate_ER(self, mode=1):
[fa5fd8d]734        # type: () -> float
[ce27e21]735        """
736        Calculate the effective radius for P(q)*S(q)
737
738        :return: the value of the effective radius
739        """
[39a06c9]740        Fq = self._calculate_Iq([0.1], True, mode)
741        return Fq[2]
[ce27e21]742
743    def calculate_VR(self):
[fa5fd8d]744        # type: () -> float
[ce27e21]745        """
746        Calculate the volf ratio for P(q)*S(q)
747
[39a06c9]748        :return: the value of the form:shell volume ratio
[ce27e21]749        """
[39a06c9]750        Fq = self._calculate_Iq([0.1], True, mode)
751        return Fq[4]
[ce27e21]752
753    def set_dispersion(self, parameter, dispersion):
[7c3fb15]754        # type: (str, weights.Dispersion) -> None
[ce27e21]755        """
756        Set the dispersion object for a model parameter
757
758        :param parameter: name of the parameter [string]
759        :param dispersion: dispersion object of type Dispersion
760        """
[fa800e72]761        if parameter in self.params:
[1780d59]762            # TODO: Store the disperser object directly in the model.
[56b2687]763            # The current method of relying on the sasview GUI to
[fa800e72]764            # remember them is kind of funky.
[1780d59]765            # Note: can't seem to get disperser parameters from sasview
[9c1a59c]766            # (1) Could create a sasview model that has not yet been
[1780d59]767            # converted, assign the disperser to one of its polydisperse
768            # parameters, then retrieve the disperser parameters from the
[9c1a59c]769            # sasview model.
770            # (2) Could write a disperser parameter retriever in sasview.
771            # (3) Could modify sasview to use sasmodels.weights dispersers.
[1780d59]772            # For now, rely on the fact that the sasview only ever uses
773            # new dispersers in the set_dispersion call and create a new
774            # one instead of trying to assign parameters.
[ce27e21]775            self.dispersion[parameter] = dispersion.get_pars()
776        else:
[7c3fb15]777            raise ValueError("%r is not a dispersity or orientation parameter"
778                             % parameter)
[ce27e21]779
[aa4946b]780    def _dispersion_mesh(self):
[fa5fd8d]781        # type: () -> List[Tuple[np.ndarray, np.ndarray]]
[ce27e21]782        """
783        Create a mesh grid of dispersion parameters and weights.
784
785        Returns [p1,p2,...],w where pj is a vector of values for parameter j
786        and w is a vector containing the products for weights for each
787        parameter set in the vector.
788        """
[4bfd277]789        pars = [self._get_weights(p)
790                for p in self._model_info.parameters.call_parameters
791                if p.type == 'volume']
[9eb3632]792        return dispersion_mesh(self._model_info, pars)
[ce27e21]793
794    def _get_weights(self, par):
[fa5fd8d]795        # type: (Parameter) -> Tuple[np.ndarray, np.ndarray]
[de0c4ba]796        """
[fb5914f]797        Return dispersion weights for parameter
[de0c4ba]798        """
[fa5fd8d]799        if par.name not in self.params:
800            if par.name == self.multiplicity_info.control:
[32f87a5]801                return self.multiplicity, [self.multiplicity], [1.0]
[fa5fd8d]802            else:
[17db833]803                # For hidden parameters use default values.  This sets
804                # scale=1 and background=0 for structure factors
805                default = self._model_info.parameters.defaults.get(par.name, np.NaN)
806                return default, [default], [1.0]
[fa5fd8d]807        elif par.polydisperse:
[32f87a5]808            value = self.params[par.name]
[fb5914f]809            dis = self.dispersion[par.name]
[9c1a59c]810            if dis['type'] == 'array':
[32f87a5]811                dispersity, weight = dis['values'], dis['weights']
[9c1a59c]812            else:
[32f87a5]813                dispersity, weight = weights.get_weights(
[9c1a59c]814                    dis['type'], dis['npts'], dis['width'], dis['nsigmas'],
[32f87a5]815                    value, par.limits, par.relative_pd)
816            return value, dispersity, weight
[fb5914f]817        else:
[32f87a5]818            value = self.params[par.name]
[ce99754]819            return value, [value], [1.0]
[ce27e21]820
[12eec1e]821    @classmethod
822    def runTests(cls):
823        """
824        Run any tests built into the model and captures the test output.
825
826        Returns success flag and output
827        """
828        from .model_test import check_model
829        return check_model(cls._model_info)
830
[749a7d4]831def test_cylinder():
[fa5fd8d]832    # type: () -> float
[4d76711]833    """
[749a7d4]834    Test that the cylinder model runs, returning the value at [0.1,0.1].
[4d76711]835    """
836    Cylinder = _make_standard_model('cylinder')
[fb5914f]837    cylinder = Cylinder()
[b32dafd]838    return cylinder.evalDistribution([0.1, 0.1])
[de97440]839
[8f93522]840def test_structure_factor():
841    # type: () -> float
842    """
[749a7d4]843    Test that 2-D hardsphere model runs and doesn't produce NaN.
[8f93522]844    """
845    Model = _make_standard_model('hardsphere')
846    model = Model()
[17db833]847    value2d = model.evalDistribution([0.1, 0.1])
848    value1d = model.evalDistribution(np.array([0.1*np.sqrt(2)]))
849    #print("hardsphere", value1d, value2d)
850    if np.isnan(value1d) or np.isnan(value2d):
851        raise ValueError("hardsphere returns nan")
[8f93522]852
[ce99754]853def test_product():
854    # type: () -> float
855    """
856    Test that 2-D hardsphere model runs and doesn't produce NaN.
857    """
858    S = _make_standard_model('hayter_msa')()
859    P = _make_standard_model('cylinder')()
860    model = MultiplicationModel(P, S)
[5024a56]861    model.setParam(product.RADIUS_MODE_ID, 1.0)
[ce99754]862    value = model.evalDistribution([0.1, 0.1])
863    if np.isnan(value):
864        raise ValueError("cylinder*hatyer_msa returns null")
865
[04045f4]866def test_rpa():
867    # type: () -> float
868    """
[749a7d4]869    Test that the 2-D RPA model runs
[04045f4]870    """
871    RPA = _make_standard_model('rpa')
872    rpa = RPA(3)
[b32dafd]873    return rpa.evalDistribution([0.1, 0.1])
[04045f4]874
[749a7d4]875def test_empty_distribution():
876    # type: () -> None
877    """
878    Make sure that sasmodels returns NaN when there are no polydispersity points
879    """
880    Cylinder = _make_standard_model('cylinder')
881    cylinder = Cylinder()
882    cylinder.setParam('radius', -1.0)
883    cylinder.setParam('background', 0.)
884    Iq = cylinder.evalDistribution(np.asarray([0.1]))
[2d81cfe]885    assert Iq[0] == 0., "empty distribution fails"
[4d76711]886
887def test_model_list():
[fa5fd8d]888    # type: () -> None
[4d76711]889    """
[749a7d4]890    Make sure that all models build as sasview models
[4d76711]891    """
892    from .exception import annotate_exception
893    for name in core.list_models():
894        try:
895            _make_standard_model(name)
896        except:
897            annotate_exception("when loading "+name)
898            raise
899
[c95dfc63]900def test_old_name():
901    # type: () -> None
902    """
[a69d8cd]903    Load and run cylinder model as sas-models-CylinderModel
[c95dfc63]904    """
905    if not SUPPORT_OLD_STYLE_PLUGINS:
906        return
907    try:
908        # if sasview is not on the path then don't try to test it
909        import sas
910    except ImportError:
911        return
912    load_standard_models()
913    from sas.models.CylinderModel import CylinderModel
914    CylinderModel().evalDistribution([0.1, 0.1])
915
[05df1de]916def magnetic_demo():
917    Model = _make_standard_model('sphere')
918    model = Model()
[610ef23]919    model.setParam('sld_M0', 8)
[05df1de]920    q = np.linspace(-0.35, 0.35, 500)
921    qx, qy = np.meshgrid(q, q)
[84f2962]922    result, _ = model.calculate_Iq(qx.flatten(), qy.flatten())
[05df1de]923    result = result.reshape(qx.shape)
924
925    import pylab
926    pylab.imshow(np.log(result + 0.001))
927    pylab.show()
928
[fb5914f]929if __name__ == "__main__":
[749a7d4]930    print("cylinder(0.1,0.1)=%g"%test_cylinder())
[05df1de]931    #magnetic_demo()
[ce99754]932    #test_product()
[17db833]933    #test_structure_factor()
934    #print("rpa:", test_rpa())
[749a7d4]935    #test_empty_distribution()
Note: See TracBrowser for help on using the repository browser.