source: sasmodels/sasmodels/sasview_model.py @ a8a1d48

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since a8a1d48 was a8a1d48, checked in by GitHub <noreply@…>, 4 months ago

Set structure factor background default to 0.0. Closes #1157.

  • Property mode set to 100644
File size: 34.9 KB
Line 
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
7    from sasmodels.sasview_model import load_custom_model
8    CylinderModel = load_custom_model('sasmodels/models/cylinder.py')
9"""
10from __future__ import print_function
11
12import math
13from copy import deepcopy
14import collections
15import traceback
16import logging
17from os.path import basename, splitext, abspath, getmtime
18try:
19    import _thread as thread
20except ImportError:
21    import thread
22
23import numpy as np  # type: ignore
24
25from . import core
26from . import custom
27from . import product
28from . import generate
29from . import weights
30from . import modelinfo
31from .details import make_kernel_args, dispersion_mesh
32
33# pylint: disable=unused-import
34try:
35    from typing import (Dict, Mapping, Any, Sequence, Tuple, NamedTuple,
36                        List, Optional, Union, Callable)
37    from .modelinfo import ModelInfo, Parameter
38    from .kernel import KernelModel
39    MultiplicityInfoType = NamedTuple(
40        'MultiplicityInfo',
41        [("number", int), ("control", str), ("choices", List[str]),
42         ("x_axis_label", str)])
43    SasviewModelType = Callable[[int], "SasviewModel"]
44except ImportError:
45    pass
46# pylint: enable=unused-import
47
48logger = logging.getLogger(__name__)
49
50calculation_lock = thread.allocate_lock()
51
52#: True if pre-existing plugins, with the old names and parameters, should
53#: continue to be supported.
54SUPPORT_OLD_STYLE_PLUGINS = True
55
56# TODO: separate x_axis_label from multiplicity info
57MultiplicityInfo = collections.namedtuple(
58    'MultiplicityInfo',
59    ["number", "control", "choices", "x_axis_label"],
60)
61
62#: set of defined models (standard and custom)
63MODELS = {}  # type: Dict[str, SasviewModelType]
64# TODO: remove unused MODEL_BY_PATH cache once sasview no longer references it
65#: custom model {path: model} mapping so we can check timestamps
66MODEL_BY_PATH = {}  # type: Dict[str, SasviewModelType]
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"]
70
71def find_model(modelname):
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    """
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
86
87# TODO: figure out how to say that the return type is a subclass
88def load_standard_models():
89    # type: () -> List[SasviewModelType]
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:
98            MODELS[name] = _make_standard_model(name)
99        except Exception:
100            logger.error(traceback.format_exc())
101    if SUPPORT_OLD_STYLE_PLUGINS:
102        _register_old_models()
103
104    return list(MODELS.values())
105
106
107def load_custom_model(path):
108    # type: (str) -> SasviewModelType
109    """
110    Load a custom model given the model path.
111    """
112    #logger.info("Loading model %s", path)
113
114    # Load the kernel module.  This may already be cached by the loader, so
115    # only requires checking the timestamps of the dependents.
116    kernel_module = custom.load_custom_kernel_module(path)
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:
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]
132        if not hasattr(model, 'filename'):
133            model.filename = abspath(kernel_module.__file__).replace('.pyc', '.py')
134        if not hasattr(model, 'id'):
135            model.id = splitext(basename(model.filename))[0]
136    else:
137        model_info = modelinfo.make_model_info(kernel_module)
138        model = make_model_from_info(model_info)
139
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
145
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'
151        logger.info("Model %s already exists: using %s [%s]",
152                    _previous_name, model.name, model.filename)
153
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]
159
160
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
175def _make_standard_model(name):
176    # type: (str) -> SasviewModelType
177    """
178    Load the sasview model defined by *name*.
179
180    *name* can be a standard model name or a path to a custom model.
181
182    Returns a class that can be used directly as a sasview model.
183    """
184    kernel_module = generate.load_kernel_module(name)
185    model_info = modelinfo.make_model_info(kernel_module)
186    return make_model_from_info(model_info)
187
188
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
204
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
216def MultiplicationModel(form_factor, structure_factor):
217    # type: ("SasviewModel", "SasviewModel") -> "SasviewModel"
218    """
219    Returns a constructed product model from form_factor and structure_factor.
220    """
221    model_info = product.make_product_info(form_factor._model_info,
222                                           structure_factor._model_info)
223    ConstructedModel = make_model_from_info(model_info)
224    return ConstructedModel(form_factor.multiplicity)
225
226
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
240    # Process multiplicity
241    non_fittable = []  # type: List[str]
242    xlabel = model_info.profile_axes[0] if model_info.profile is not None else ""
243    variants = MultiplicityInfo(0, "", [], xlabel)
244    for p in model_info.parameters.kernel_parameters:
245        if p.name == model_info.control:
246            non_fittable.append(p.name)
247            variants = MultiplicityInfo(
248                len(p.choices) if p.choices else int(p.limits[1]),
249                p.name, p.choices, xlabel
250            )
251            break
252
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
262    # Organize parameter sets
263    orientation_params = []
264    magnetic_params = []
265    fixed = []
266    for p in model_info.parameters.user_parameters({}, is2d=True):
267        if p.type == 'orientation':
268            orientation_params.append(p.name)
269            orientation_params.append(p.name+".width")
270            fixed.append(p.name+".width")
271        elif p.type == 'magnetic':
272            orientation_params.append(p.name)
273            magnetic_params.append(p.name)
274            fixed.append(p.name+".width")
275
276
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
285    attrs['is_form_factor'] = model_info.effective_radius_type is not None
286    attrs['is_multiplicity_model'] = variants[0] > 1
287    attrs['multiplicity_info'] = variants
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)
292    attrs['fun_list'] = tuple(fun_list)
293
294    return attrs
295
296class SasviewModel(object):
297    """
298    Sasview wrapper for opencl/ctypes model.
299    """
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
316    orientation_params = None # type: List[str]
317    #: names of the magnetic parameters in the order they appear
318    magnetic_params = None    # type: List[str]
319    #: names of the fittable parameters
320    fixed = None              # type: List[str]
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
342    #: Multiplicity information
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
351    details = None     # type: Dict[str, Sequence[Any]]
352    #                  # actual type is Dict[str, List[str, float, float]]
353    #: multiplicity value, or None if no multiplicity on the model
354    multiplicity = None     # type: Optional[int]
355    #: memory for polydispersity array if using ArrayDispersion (used by sasview).
356    _persistency_dict = None # type: Dict[str, Tuple[np.ndarray, np.ndarray]]
357
358    def __init__(self, multiplicity=None):
359        # type: (Optional[int]) -> None
360
361        # TODO: _persistency_dict to persistency_dict throughout sasview
362        # TODO: refactor multiplicity to encompass variants
363        # TODO: dispersion should be a class
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
373        # Get the list of hidden parameters given the multiplicity
374        # Don't include multiplicity in the list of parameters
375        self.multiplicity = multiplicity
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()
381        if self._model_info.structure_factor:
382            hidden.add('scale')
383            hidden.add('background')
384
385        # Update the parameter lists to exclude any hidden parameters
386        self.magnetic_params = tuple(pname for pname in self.magnetic_params
387                                     if pname not in hidden)
388        self.orientation_params = tuple(pname for pname in self.orientation_params
389                                        if pname not in hidden)
390
391        self._persistency_dict = {}
392        self.params = collections.OrderedDict()
393        self.dispersion = collections.OrderedDict()
394        self.details = {}
395        for p in self._model_info.parameters.user_parameters({}, is2d=True):
396            if p.name in hidden:
397                continue
398            self.params[p.name] = p.default
399            self.details[p.id] = [p.units, p.limits[0], p.limits[1]]
400            if p.polydisperse:
401                self.details[p.id+".width"] = [
402                    "", 0.0, 1.0 if p.relative_pd else np.inf
403                ]
404                self.dispersion[p.name] = {
405                    'width': 0,
406                    'npts': 35,
407                    'nsigmas': 3,
408                    'type': 'gaussian',
409                }
410
411    def __get_state__(self):
412        # type: () -> Dict[str, Any]
413        state = self.__dict__.copy()
414        state.pop('_model')
415        # May need to reload model info on set state since it has pointers
416        # to python implementations of Iq, etc.
417        #state.pop('_model_info')
418        return state
419
420    def __set_state__(self, state):
421        # type: (Dict[str, Any]) -> None
422        self.__dict__ = state
423        self._model = None
424
425    def __str__(self):
426        # type: () -> str
427        """
428        :return: string representation
429        """
430        return self.name
431
432    def is_fittable(self, par_name):
433        # type: (str) -> bool
434        """
435        Check if a given parameter is fittable or not
436
437        :param par_name: the parameter name to check
438        """
439        return par_name in self.fixed
440        #For the future
441        #return self.params[str(par_name)].is_fittable()
442
443
444    def getProfile(self):
445        # type: () -> (np.ndarray, np.ndarray)
446        """
447        Get SLD profile
448
449        : return: (z, beta) where z is a list of depth of the transition points
450                beta is a list of the corresponding SLD values
451        """
452        args = {} # type: Dict[str, Any]
453        for p in self._model_info.parameters.kernel_parameters:
454            if p.id == self.multiplicity_info.control:
455                value = float(self.multiplicity)
456            elif p.length == 1:
457                value = self.params.get(p.id, np.NaN)
458            else:
459                value = np.array([self.params.get(p.id+str(k), np.NaN)
460                                  for k in range(1, p.length+1)])
461            args[p.id] = value
462
463        x, y = self._model_info.profile(**args)
464        return x, 1e-6*y
465
466    def setParam(self, name, value):
467        # type: (str, float) -> None
468        """
469        Set the value of a model parameter
470
471        :param name: name of the parameter
472        :param value: value of the parameter
473
474        """
475        # Look for dispersion parameters
476        toks = name.split('.')
477        if len(toks) == 2:
478            for item in self.dispersion.keys():
479                if item == toks[0]:
480                    for par in self.dispersion[item]:
481                        if par == toks[1]:
482                            self.dispersion[item][par] = value
483                            return
484        else:
485            # Look for standard parameter
486            for item in self.params.keys():
487                if item == name:
488                    self.params[item] = value
489                    return
490
491        raise ValueError("Model does not contain parameter %s" % name)
492
493    def getParam(self, name):
494        # type: (str) -> float
495        """
496        Set the value of a model parameter
497
498        :param name: name of the parameter
499
500        """
501        # Look for dispersion parameters
502        toks = name.split('.')
503        if len(toks) == 2:
504            for item in self.dispersion.keys():
505                if item == toks[0]:
506                    for par in self.dispersion[item]:
507                        if par == toks[1]:
508                            return self.dispersion[item][par]
509        else:
510            # Look for standard parameter
511            for item in self.params.keys():
512                if item == name:
513                    return self.params[item]
514
515        raise ValueError("Model does not contain parameter %s" % name)
516
517    def getParamList(self):
518        # type: () -> Sequence[str]
519        """
520        Return a list of all available parameters for the model
521        """
522        param_list = list(self.params.keys())
523        # WARNING: Extending the list with the dispersion parameters
524        param_list.extend(self.getDispParamList())
525        return param_list
526
527    def getDispParamList(self):
528        # type: () -> Sequence[str]
529        """
530        Return a list of polydispersity parameters for the model
531        """
532        # TODO: fix test so that parameter order doesn't matter
533        ret = ['%s.%s' % (p_name, ext)
534               for p_name in self.dispersion.keys()
535               for ext in ('npts', 'nsigmas', 'width')]
536        #print(ret)
537        return ret
538
539    def clone(self):
540        # type: () -> "SasviewModel"
541        """ Return a identical copy of self """
542        return deepcopy(self)
543
544    def run(self, x=0.0):
545        # type: (Union[float, (float, float), List[float]]) -> float
546        """
547        Evaluate the model
548
549        :param x: input q, or [q,phi]
550
551        :return: scattering function P(q)
552
553        **DEPRECATED**: use calculate_Iq instead
554        """
555        if isinstance(x, (list, tuple)):
556            # pylint: disable=unpacking-non-sequence
557            q, phi = x
558            result, _ = self.calculate_Iq([q*math.cos(phi)], [q*math.sin(phi)])
559            return result[0]
560        else:
561            result, _ = self.calculate_Iq([x])
562            return result[0]
563
564
565    def runXY(self, x=0.0):
566        # type: (Union[float, (float, float), List[float]]) -> float
567        """
568        Evaluate the model in cartesian coordinates
569
570        :param x: input q, or [qx, qy]
571
572        :return: scattering function P(q)
573
574        **DEPRECATED**: use calculate_Iq instead
575        """
576        if isinstance(x, (list, tuple)):
577            result, _ = self.calculate_Iq([x[0]], [x[1]])
578            return result[0]
579        else:
580            result, _ = self.calculate_Iq([x])
581            return result[0]
582
583    def evalDistribution(self, qdist):
584        # type: (Union[np.ndarray, Tuple[np.ndarray, np.ndarray], List[np.ndarray]]) -> np.ndarray
585        r"""
586        Evaluate a distribution of q-values.
587
588        :param qdist: array of q or a list of arrays [qx,qy]
589
590        * For 1D, a numpy array is expected as input
591
592        ::
593
594            evalDistribution(q)
595
596          where *q* is a numpy array.
597
598        * For 2D, a list of *[qx,qy]* is expected with 1D arrays as input
599
600        ::
601
602              qx = [ qx[0], qx[1], qx[2], ....]
603              qy = [ qy[0], qy[1], qy[2], ....]
604
605        If the model is 1D only, then
606
607        .. math::
608
609            q = \sqrt{q_x^2+q_y^2}
610
611        """
612        if isinstance(qdist, (list, tuple)):
613            # Check whether we have a list of ndarrays [qx,qy]
614            qx, qy = qdist
615            result, _ = self.calculate_Iq(qx, qy)
616            return result
617
618        elif isinstance(qdist, np.ndarray):
619            # We have a simple 1D distribution of q-values
620            result, _ = self.calculate_Iq(qdist)
621            return result
622
623        else:
624            raise TypeError("evalDistribution expects q or [qx, qy], not %r"
625                            % type(qdist))
626
627    def calc_composition_models(self, qx):
628        """
629        returns parts of the composition model or None if not a composition
630        model.
631        """
632        # TODO: have calculate_Iq return the intermediates.
633        #
634        # The current interface causes calculate_Iq() to be called twice,
635        # once to get the combined result and again to get the intermediate
636        # results.  This is necessary for now.
637        # Long term, the solution is to change the interface to calculate_Iq
638        # so that it returns a results object containing all the bits:
639        #     the A, B, C, ... of the composition model (and any subcomponents?)
640        #     the P and S of the product model
641        #     the combined model before resolution smearing,
642        #     the sasmodel before sesans conversion,
643        #     the oriented 2D model used to fit oriented usans data,
644        #     the final I(q),
645        #     ...
646        #
647        # Have the model calculator add all of these blindly to the data
648        # tree, and update the graphs which contain them.  The fitter
649        # needs to be updated to use the I(q) value only, ignoring the rest.
650        #
651        # The simple fix of returning the existing intermediate results
652        # will not work for a couple of reasons: (1) another thread may
653        # sneak in to compute its own results before calc_composition_models
654        # is called, and (2) calculate_Iq is currently called three times:
655        # once with q, once with q values before qmin and once with q values
656        # after q max.  Both of these should be addressed before
657        # replacing this code.
658        composition = self._model_info.composition
659        if composition and composition[0] == 'product': # only P*S for now
660            with calculation_lock:
661                _, lazy_results = self._calculate_Iq(qx)
662                # for compatibility with sasview 4.x
663                results = lazy_results()
664                pq_data = results.get("P(Q)")
665                sq_data = results.get("S(Q)")
666                return pq_data, sq_data
667        else:
668            return None
669
670    def calculate_Iq(self,
671                     qx,     # type: Sequence[float]
672                     qy=None # type: Optional[Sequence[float]]
673                     ):
674        # type: (...) -> Tuple[np.ndarray, Callable[[], collections.OrderedDict[str, np.ndarray]]]
675        """
676        Calculate Iq for one set of q with the current parameters.
677
678        If the model is 1D, use *q*.  If 2D, use *qx*, *qy*.
679
680        This should NOT be used for fitting since it copies the *q* vectors
681        to the card for each evaluation.
682
683        The returned tuple contains the scattering intensity followed by a
684        callable which returns a dictionary of intermediate data from
685        ProductKernel.
686        """
687        ## uncomment the following when trying to debug the uncoordinated calls
688        ## to calculate_Iq
689        #if calculation_lock.locked():
690        #    logger.info("calculation waiting for another thread to complete")
691        #    logger.info("\n".join(traceback.format_stack()))
692
693        with calculation_lock:
694            return self._calculate_Iq(qx, qy)
695
696    def _calculate_Iq(self, qx, qy=None):
697        if self._model is None:
698            self._model = core.build_model(self._model_info)
699        if qy is not None:
700            q_vectors = [np.asarray(qx), np.asarray(qy)]
701        else:
702            q_vectors = [np.asarray(qx)]
703        calculator = self._model.make_kernel(q_vectors)
704        parameters = self._model_info.parameters
705        pairs = [self._get_weights(p) for p in parameters.call_parameters]
706        #weights.plot_weights(self._model_info, pairs)
707        call_details, values, is_magnetic = make_kernel_args(calculator, pairs)
708        #call_details.show()
709        #print("================ parameters ==================")
710        #for p, v in zip(parameters.call_parameters, pairs): print(p.name, v[0])
711        #for k, p in enumerate(self._model_info.parameters.call_parameters):
712        #    print(k, p.name, *pairs[k])
713        #print("params", self.params)
714        #print("values", values)
715        #print("is_mag", is_magnetic)
716        result = calculator(call_details, values, cutoff=self.cutoff,
717                            magnetic=is_magnetic)
718        lazy_results = getattr(calculator, 'results',
719                               lambda: collections.OrderedDict())
720        #print("result", result)
721
722        calculator.release()
723        #self._model.release()
724
725        return result, lazy_results
726
727
728    def calculate_ER(self, mode=1):
729        # type: () -> float
730        """
731        Calculate the effective radius for P(q)*S(q)
732
733        *mode* is the R_eff type, which defaults to 1 to match the ER
734        calculation for sasview models from version 3.x.
735
736        :return: the value of the effective radius
737        """
738        # ER and VR are only needed for old multiplication models, based on
739        # sas.sascalc.fit.MultiplicationModel.  Fail for now.  If we want to
740        # continue supporting them then add some test cases so that the code
741        # is exercised.  We can access ER/VR using the kernel Fq function by
742        # extending _calculate_Iq so that it calls:
743        #    if er_mode > 0:
744        #        res = calculator.Fq(call_details, values, cutoff=self.cutoff,
745        #                            magnetic=False, effective_radius_type=mode)
746        #        R_eff, form_shell_ratio = res[2], res[4]
747        #        return R_eff, form_shell_ratio
748        # Then use the following in calculate_ER:
749        #    ER, VR = self._calculate_Iq(q=[0.1], er_mode=mode)
750        #    return ER
751        # Similarly, for calculate_VR:
752        #    ER, VR = self._calculate_Iq(q=[0.1], er_mode=1)
753        #    return VR
754        # Obviously a combined calculate_ER_VR method would be better, but
755        # we only need them to support very old models, so ignore the 2x
756        # performance hit.
757        raise NotImplementedError("ER function is no longer available.")
758
759    def calculate_VR(self):
760        # type: () -> float
761        """
762        Calculate the volf ratio for P(q)*S(q)
763
764        :return: the value of the form:shell volume ratio
765        """
766        # See comments in calculate_ER.
767        raise NotImplementedError("VR function is no longer available.")
768
769    def set_dispersion(self, parameter, dispersion):
770        # type: (str, weights.Dispersion) -> None
771        """
772        Set the dispersion object for a model parameter
773
774        :param parameter: name of the parameter [string]
775        :param dispersion: dispersion object of type Dispersion
776        """
777        if parameter in self.params:
778            # TODO: Store the disperser object directly in the model.
779            # The current method of relying on the sasview GUI to
780            # remember them is kind of funky.
781            # Note: can't seem to get disperser parameters from sasview
782            # (1) Could create a sasview model that has not yet been
783            # converted, assign the disperser to one of its polydisperse
784            # parameters, then retrieve the disperser parameters from the
785            # sasview model.
786            # (2) Could write a disperser parameter retriever in sasview.
787            # (3) Could modify sasview to use sasmodels.weights dispersers.
788            # For now, rely on the fact that the sasview only ever uses
789            # new dispersers in the set_dispersion call and create a new
790            # one instead of trying to assign parameters.
791            self.dispersion[parameter] = dispersion.get_pars()
792        else:
793            raise ValueError("%r is not a dispersity or orientation parameter"
794                             % parameter)
795
796    def _dispersion_mesh(self):
797        # type: () -> List[Tuple[np.ndarray, np.ndarray]]
798        """
799        Create a mesh grid of dispersion parameters and weights.
800
801        Returns [p1,p2,...],w where pj is a vector of values for parameter j
802        and w is a vector containing the products for weights for each
803        parameter set in the vector.
804        """
805        pars = [self._get_weights(p)
806                for p in self._model_info.parameters.call_parameters
807                if p.type == 'volume']
808        return dispersion_mesh(self._model_info, pars)
809
810    def _get_weights(self, par):
811        # type: (Parameter) -> Tuple[np.ndarray, np.ndarray]
812        """
813        Return dispersion weights for parameter
814        """
815        if par.name not in self.params:
816            if par.name == self.multiplicity_info.control:
817                return self.multiplicity, [self.multiplicity], [1.0]
818            else:
819                # For hidden parameters use default values.  This sets
820                # scale=1 and background=0 for structure factors
821                default = self._model_info.parameters.defaults.get(par.name, np.NaN)
822                return default, [default], [1.0]
823        elif par.polydisperse:
824            value = self.params[par.name]
825            dis = self.dispersion[par.name]
826            if dis['type'] == 'array':
827                dispersity, weight = dis['values'], dis['weights']
828            else:
829                dispersity, weight = weights.get_weights(
830                    dis['type'], dis['npts'], dis['width'], dis['nsigmas'],
831                    value, par.limits, par.relative_pd)
832            return value, dispersity, weight
833        else:
834            value = self.params[par.name]
835            return value, [value], [1.0]
836
837    @classmethod
838    def runTests(cls):
839        """
840        Run any tests built into the model and captures the test output.
841
842        Returns success flag and output
843        """
844        from .model_test import check_model
845        return check_model(cls._model_info)
846
847def test_cylinder():
848    # type: () -> float
849    """
850    Test that the cylinder model runs, returning the value at [0.1,0.1].
851    """
852    Cylinder = _make_standard_model('cylinder')
853    cylinder = Cylinder()
854    return cylinder.evalDistribution([0.1, 0.1])
855
856def test_structure_factor():
857    # type: () -> float
858    """
859    Test that 2-D hardsphere model runs and doesn't produce NaN.
860    """
861    Model = _make_standard_model('hardsphere')
862    model = Model()
863    value2d = model.evalDistribution([0.1, 0.1])
864    value1d = model.evalDistribution(np.array([0.1*np.sqrt(2)]))
865    #print("hardsphere", value1d, value2d)
866    if np.isnan(value1d) or np.isnan(value2d):
867        raise ValueError("hardsphere returns nan")
868
869def test_product():
870    # type: () -> float
871    """
872    Test that 2-D hardsphere model runs and doesn't produce NaN.
873    """
874    S = _make_standard_model('hayter_msa')()
875    P = _make_standard_model('cylinder')()
876    model = MultiplicationModel(P, S)
877    model.setParam(product.RADIUS_MODE_ID, 1.0)
878    value = model.evalDistribution([0.1, 0.1])
879    if np.isnan(value):
880        raise ValueError("cylinder*hatyer_msa returns null")
881
882def test_rpa():
883    # type: () -> float
884    """
885    Test that the 2-D RPA model runs
886    """
887    RPA = _make_standard_model('rpa')
888    rpa = RPA(3)
889    return rpa.evalDistribution([0.1, 0.1])
890
891def test_empty_distribution():
892    # type: () -> None
893    """
894    Make sure that sasmodels returns NaN when there are no polydispersity points
895    """
896    Cylinder = _make_standard_model('cylinder')
897    cylinder = Cylinder()
898    cylinder.setParam('radius', -1.0)
899    cylinder.setParam('background', 0.)
900    Iq = cylinder.evalDistribution(np.asarray([0.1]))
901    assert Iq[0] == 0., "empty distribution fails"
902
903def test_model_list():
904    # type: () -> None
905    """
906    Make sure that all models build as sasview models
907    """
908    from .exception import annotate_exception
909    for name in core.list_models():
910        try:
911            _make_standard_model(name)
912        except:
913            annotate_exception("when loading "+name)
914            raise
915
916def test_old_name():
917    # type: () -> None
918    """
919    Load and run cylinder model as sas-models-CylinderModel
920    """
921    if not SUPPORT_OLD_STYLE_PLUGINS:
922        return
923    try:
924        # if sasview is not on the path then don't try to test it
925        import sas
926    except ImportError:
927        return
928    load_standard_models()
929    from sas.models.CylinderModel import CylinderModel
930    CylinderModel().evalDistribution([0.1, 0.1])
931
932def test_structure_factor_background():
933    # type: () -> None
934    """
935    Check that sasview model and direct model match, with background=0.
936    """
937    from .data import empty_data1D
938    from .core import load_model_info, build_model
939    from .direct_model import DirectModel
940
941    model_name = "hardsphere"
942    q = [0.0]
943
944    sasview_model = _make_standard_model(model_name)()
945    sasview_value = sasview_model.evalDistribution(np.array(q))[0]
946
947    data = empty_data1D(q)
948    model_info = load_model_info(model_name)
949    model = build_model(model_info)
950    direct_model = DirectModel(data, model)
951    direct_value_zero_background = direct_model(background=0.0)
952
953    assert sasview_value == direct_value_zero_background
954
955    # Additionally check that direct value background defaults to zero
956    direct_value_default = direct_model()
957    assert sasview_value == direct_value_default
958
959
960def magnetic_demo():
961    Model = _make_standard_model('sphere')
962    model = Model()
963    model.setParam('sld_M0', 8)
964    q = np.linspace(-0.35, 0.35, 500)
965    qx, qy = np.meshgrid(q, q)
966    result, _ = model.calculate_Iq(qx.flatten(), qy.flatten())
967    result = result.reshape(qx.shape)
968
969    import pylab
970    pylab.imshow(np.log(result + 0.001))
971    pylab.show()
972
973if __name__ == "__main__":
974    print("cylinder(0.1,0.1)=%g"%test_cylinder())
975    #magnetic_demo()
976    #test_product()
977    #test_structure_factor()
978    #print("rpa:", test_rpa())
979    #test_empty_distribution()
980    #test_structure_factor_background()
Note: See TracBrowser for help on using the repository browser.