source: sasmodels/sasmodels/sasview_model.py @ 3a1afed

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

remove calculate_ER/calculate_VR since it was broken and only needed for 3.x multiplication models

  • Property mode set to 100644
File size: 34.1 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            self._model_info.parameters.defaults['background'] = 0.
385
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
392        self._persistency_dict = {}
393        self.params = collections.OrderedDict()
394        self.dispersion = collections.OrderedDict()
395        self.details = {}
396        for p in self._model_info.parameters.user_parameters({}, is2d=True):
397            if p.name in hidden:
398                continue
399            self.params[p.name] = p.default
400            self.details[p.id] = [p.units, p.limits[0], p.limits[1]]
401            if p.polydisperse:
402                self.details[p.id+".width"] = [
403                    "", 0.0, 1.0 if p.relative_pd else np.inf
404                ]
405                self.dispersion[p.name] = {
406                    'width': 0,
407                    'npts': 35,
408                    'nsigmas': 3,
409                    'type': 'gaussian',
410                }
411
412    def __get_state__(self):
413        # type: () -> Dict[str, Any]
414        state = self.__dict__.copy()
415        state.pop('_model')
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):
422        # type: (Dict[str, Any]) -> None
423        self.__dict__ = state
424        self._model = None
425
426    def __str__(self):
427        # type: () -> str
428        """
429        :return: string representation
430        """
431        return self.name
432
433    def is_fittable(self, par_name):
434        # type: (str) -> bool
435        """
436        Check if a given parameter is fittable or not
437
438        :param par_name: the parameter name to check
439        """
440        return par_name in self.fixed
441        #For the future
442        #return self.params[str(par_name)].is_fittable()
443
444
445    def getProfile(self):
446        # type: () -> (np.ndarray, np.ndarray)
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        """
453        args = {} # type: Dict[str, Any]
454        for p in self._model_info.parameters.kernel_parameters:
455            if p.id == self.multiplicity_info.control:
456                value = float(self.multiplicity)
457            elif p.length == 1:
458                value = self.params.get(p.id, np.NaN)
459            else:
460                value = np.array([self.params.get(p.id+str(k), np.NaN)
461                                  for k in range(1, p.length+1)])
462            args[p.id] = value
463
464        x, y = self._model_info.profile(**args)
465        return x, 1e-6*y
466
467    def setParam(self, name, value):
468        # type: (str, float) -> None
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('.')
478        if len(toks) == 2:
479            for item in self.dispersion.keys():
480                if item == toks[0]:
481                    for par in self.dispersion[item]:
482                        if par == toks[1]:
483                            self.dispersion[item][par] = value
484                            return
485        else:
486            # Look for standard parameter
487            for item in self.params.keys():
488                if item == name:
489                    self.params[item] = value
490                    return
491
492        raise ValueError("Model does not contain parameter %s" % name)
493
494    def getParam(self, name):
495        # type: (str) -> float
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('.')
504        if len(toks) == 2:
505            for item in self.dispersion.keys():
506                if item == toks[0]:
507                    for par in self.dispersion[item]:
508                        if par == toks[1]:
509                            return self.dispersion[item][par]
510        else:
511            # Look for standard parameter
512            for item in self.params.keys():
513                if item == name:
514                    return self.params[item]
515
516        raise ValueError("Model does not contain parameter %s" % name)
517
518    def getParamList(self):
519        # type: () -> Sequence[str]
520        """
521        Return a list of all available parameters for the model
522        """
523        param_list = list(self.params.keys())
524        # WARNING: Extending the list with the dispersion parameters
525        param_list.extend(self.getDispParamList())
526        return param_list
527
528    def getDispParamList(self):
529        # type: () -> Sequence[str]
530        """
531        Return a list of polydispersity parameters for the model
532        """
533        # TODO: fix test so that parameter order doesn't matter
534        ret = ['%s.%s' % (p_name, ext)
535               for p_name in self.dispersion.keys()
536               for ext in ('npts', 'nsigmas', 'width')]
537        #print(ret)
538        return ret
539
540    def clone(self):
541        # type: () -> "SasviewModel"
542        """ Return a identical copy of self """
543        return deepcopy(self)
544
545    def run(self, x=0.0):
546        # type: (Union[float, (float, float), List[float]]) -> float
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        """
556        if isinstance(x, (list, tuple)):
557            # pylint: disable=unpacking-non-sequence
558            q, phi = x
559            result, _ = self.calculate_Iq([q*math.cos(phi)], [q*math.sin(phi)])
560            return result[0]
561        else:
562            result, _ = self.calculate_Iq([x])
563            return result[0]
564
565
566    def runXY(self, x=0.0):
567        # type: (Union[float, (float, float), List[float]]) -> float
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        """
577        if isinstance(x, (list, tuple)):
578            result, _ = self.calculate_Iq([x[0]], [x[1]])
579            return result[0]
580        else:
581            result, _ = self.calculate_Iq([x])
582            return result[0]
583
584    def evalDistribution(self, qdist):
585        # type: (Union[np.ndarray, Tuple[np.ndarray, np.ndarray], List[np.ndarray]]) -> np.ndarray
586        r"""
587        Evaluate a distribution of q-values.
588
589        :param qdist: array of q or a list of arrays [qx,qy]
590
591        * For 1D, a numpy array is expected as input
592
593        ::
594
595            evalDistribution(q)
596
597          where *q* is a numpy array.
598
599        * For 2D, a list of *[qx,qy]* is expected with 1D arrays as input
600
601        ::
602
603              qx = [ qx[0], qx[1], qx[2], ....]
604              qy = [ qy[0], qy[1], qy[2], ....]
605
606        If the model is 1D only, then
607
608        .. math::
609
610            q = \sqrt{q_x^2+q_y^2}
611
612        """
613        if isinstance(qdist, (list, tuple)):
614            # Check whether we have a list of ndarrays [qx,qy]
615            qx, qy = qdist
616            result, _ = self.calculate_Iq(qx, qy)
617            return result
618
619        elif isinstance(qdist, np.ndarray):
620            # We have a simple 1D distribution of q-values
621            result, _ = self.calculate_Iq(qdist)
622            return result
623
624        else:
625            raise TypeError("evalDistribution expects q or [qx, qy], not %r"
626                            % type(qdist))
627
628    def calc_composition_models(self, qx):
629        """
630        returns parts of the composition model or None if not a composition
631        model.
632        """
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:
640        #     the A, B, C, ... of the composition model (and any subcomponents?)
641        #     the P and S of the product model
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        #     ...
647        #
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.
659        composition = self._model_info.composition
660        if composition and composition[0] == 'product': # only P*S for now
661            with calculation_lock:
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
668        else:
669            return None
670
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]]]
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.
683
684        The returned tuple contains the scattering intensity followed by a
685        callable which returns a dictionary of intermediate data from
686        ProductKernel.
687        """
688        ## uncomment the following when trying to debug the uncoordinated calls
689        ## to calculate_Iq
690        #if calculation_lock.locked():
691        #    logger.info("calculation waiting for another thread to complete")
692        #    logger.info("\n".join(traceback.format_stack()))
693
694        with calculation_lock:
695            return self._calculate_Iq(qx, qy)
696
697    def _calculate_Iq(self, qx, qy=None):
698        if self._model is None:
699            self._model = core.build_model(self._model_info)
700        if qy is not None:
701            q_vectors = [np.asarray(qx), np.asarray(qy)]
702        else:
703            q_vectors = [np.asarray(qx)]
704        calculator = self._model.make_kernel(q_vectors)
705        parameters = self._model_info.parameters
706        pairs = [self._get_weights(p) for p in parameters.call_parameters]
707        #weights.plot_weights(self._model_info, pairs)
708        call_details, values, is_magnetic = make_kernel_args(calculator, pairs)
709        #call_details.show()
710        #print("================ parameters ==================")
711        #for p, v in zip(parameters.call_parameters, pairs): print(p.name, v[0])
712        #for k, p in enumerate(self._model_info.parameters.call_parameters):
713        #    print(k, p.name, *pairs[k])
714        #print("params", self.params)
715        #print("values", values)
716        #print("is_mag", is_magnetic)
717        result = calculator(call_details, values, cutoff=self.cutoff,
718                            magnetic=is_magnetic)
719        lazy_results = getattr(calculator, 'results',
720                               lambda: collections.OrderedDict())
721        #print("result", result)
722
723        calculator.release()
724        #self._model.release()
725
726        return result, lazy_results
727
728
729    def calculate_ER(self, mode=1):
730        # type: () -> float
731        """
732        Calculate the effective radius for P(q)*S(q)
733
734        *mode* is the R_eff type, which defaults to 1 to match the ER
735        calculation for sasview models from version 3.x.
736
737        :return: the value of the effective radius
738        """
739        # ER and VR are only needed for old multiplication models, based on
740        # sas.sascalc.fit.MultiplicationModel.  Fail for now.  If we want to
741        # continue supporting them then add some test cases so that the code
742        # is exercised.  We can access ER/VR using the kernel Fq function by
743        # extending _calculate_Iq so that it calls:
744        #    if er_mode > 0:
745        #        res = calculator.Fq(call_details, values, cutoff=self.cutoff,
746        #                            magnetic=False, effective_radius_type=mode)
747        #        R_eff, form_shell_ratio = res[2], res[4]
748        #        return R_eff, form_shell_ratio
749        # Then use the following in calculate_ER:
750        #    ER, VR = self._calculate_Iq(q=[0.1], er_mode=mode)
751        #    return ER
752        # Similarly, for calculate_VR:
753        #    ER, VR = self._calculate_Iq(q=[0.1], er_mode=1)
754        #    return VR
755        # Obviously a combined calculate_ER_VR method would be better, but
756        # we only need them to support very old models, so ignore the 2x
757        # performance hit.
758        raise NotImplementedError("ER function is no longer available.")
759
760    def calculate_VR(self):
761        # type: () -> float
762        """
763        Calculate the volf ratio for P(q)*S(q)
764
765        :return: the value of the form:shell volume ratio
766        """
767        # See comments in calculate_ER.
768        raise NotImplementedError("VR function is no longer available.")
769
770    def set_dispersion(self, parameter, dispersion):
771        # type: (str, weights.Dispersion) -> None
772        """
773        Set the dispersion object for a model parameter
774
775        :param parameter: name of the parameter [string]
776        :param dispersion: dispersion object of type Dispersion
777        """
778        if parameter in self.params:
779            # TODO: Store the disperser object directly in the model.
780            # The current method of relying on the sasview GUI to
781            # remember them is kind of funky.
782            # Note: can't seem to get disperser parameters from sasview
783            # (1) Could create a sasview model that has not yet been
784            # converted, assign the disperser to one of its polydisperse
785            # parameters, then retrieve the disperser parameters from the
786            # sasview model.
787            # (2) Could write a disperser parameter retriever in sasview.
788            # (3) Could modify sasview to use sasmodels.weights dispersers.
789            # For now, rely on the fact that the sasview only ever uses
790            # new dispersers in the set_dispersion call and create a new
791            # one instead of trying to assign parameters.
792            self.dispersion[parameter] = dispersion.get_pars()
793        else:
794            raise ValueError("%r is not a dispersity or orientation parameter"
795                             % parameter)
796
797    def _dispersion_mesh(self):
798        # type: () -> List[Tuple[np.ndarray, np.ndarray]]
799        """
800        Create a mesh grid of dispersion parameters and weights.
801
802        Returns [p1,p2,...],w where pj is a vector of values for parameter j
803        and w is a vector containing the products for weights for each
804        parameter set in the vector.
805        """
806        pars = [self._get_weights(p)
807                for p in self._model_info.parameters.call_parameters
808                if p.type == 'volume']
809        return dispersion_mesh(self._model_info, pars)
810
811    def _get_weights(self, par):
812        # type: (Parameter) -> Tuple[np.ndarray, np.ndarray]
813        """
814        Return dispersion weights for parameter
815        """
816        if par.name not in self.params:
817            if par.name == self.multiplicity_info.control:
818                return self.multiplicity, [self.multiplicity], [1.0]
819            else:
820                # For hidden parameters use default values.  This sets
821                # scale=1 and background=0 for structure factors
822                default = self._model_info.parameters.defaults.get(par.name, np.NaN)
823                return default, [default], [1.0]
824        elif par.polydisperse:
825            value = self.params[par.name]
826            dis = self.dispersion[par.name]
827            if dis['type'] == 'array':
828                dispersity, weight = dis['values'], dis['weights']
829            else:
830                dispersity, weight = weights.get_weights(
831                    dis['type'], dis['npts'], dis['width'], dis['nsigmas'],
832                    value, par.limits, par.relative_pd)
833            return value, dispersity, weight
834        else:
835            value = self.params[par.name]
836            return value, [value], [1.0]
837
838    @classmethod
839    def runTests(cls):
840        """
841        Run any tests built into the model and captures the test output.
842
843        Returns success flag and output
844        """
845        from .model_test import check_model
846        return check_model(cls._model_info)
847
848def test_cylinder():
849    # type: () -> float
850    """
851    Test that the cylinder model runs, returning the value at [0.1,0.1].
852    """
853    Cylinder = _make_standard_model('cylinder')
854    cylinder = Cylinder()
855    return cylinder.evalDistribution([0.1, 0.1])
856
857def test_structure_factor():
858    # type: () -> float
859    """
860    Test that 2-D hardsphere model runs and doesn't produce NaN.
861    """
862    Model = _make_standard_model('hardsphere')
863    model = Model()
864    value2d = model.evalDistribution([0.1, 0.1])
865    value1d = model.evalDistribution(np.array([0.1*np.sqrt(2)]))
866    #print("hardsphere", value1d, value2d)
867    if np.isnan(value1d) or np.isnan(value2d):
868        raise ValueError("hardsphere returns nan")
869
870def test_product():
871    # type: () -> float
872    """
873    Test that 2-D hardsphere model runs and doesn't produce NaN.
874    """
875    S = _make_standard_model('hayter_msa')()
876    P = _make_standard_model('cylinder')()
877    model = MultiplicationModel(P, S)
878    model.setParam(product.RADIUS_MODE_ID, 1.0)
879    value = model.evalDistribution([0.1, 0.1])
880    if np.isnan(value):
881        raise ValueError("cylinder*hatyer_msa returns null")
882
883def test_rpa():
884    # type: () -> float
885    """
886    Test that the 2-D RPA model runs
887    """
888    RPA = _make_standard_model('rpa')
889    rpa = RPA(3)
890    return rpa.evalDistribution([0.1, 0.1])
891
892def test_empty_distribution():
893    # type: () -> None
894    """
895    Make sure that sasmodels returns NaN when there are no polydispersity points
896    """
897    Cylinder = _make_standard_model('cylinder')
898    cylinder = Cylinder()
899    cylinder.setParam('radius', -1.0)
900    cylinder.setParam('background', 0.)
901    Iq = cylinder.evalDistribution(np.asarray([0.1]))
902    assert Iq[0] == 0., "empty distribution fails"
903
904def test_model_list():
905    # type: () -> None
906    """
907    Make sure that all models build as sasview models
908    """
909    from .exception import annotate_exception
910    for name in core.list_models():
911        try:
912            _make_standard_model(name)
913        except:
914            annotate_exception("when loading "+name)
915            raise
916
917def test_old_name():
918    # type: () -> None
919    """
920    Load and run cylinder model as sas-models-CylinderModel
921    """
922    if not SUPPORT_OLD_STYLE_PLUGINS:
923        return
924    try:
925        # if sasview is not on the path then don't try to test it
926        import sas
927    except ImportError:
928        return
929    load_standard_models()
930    from sas.models.CylinderModel import CylinderModel
931    CylinderModel().evalDistribution([0.1, 0.1])
932
933def magnetic_demo():
934    Model = _make_standard_model('sphere')
935    model = Model()
936    model.setParam('sld_M0', 8)
937    q = np.linspace(-0.35, 0.35, 500)
938    qx, qy = np.meshgrid(q, q)
939    result, _ = model.calculate_Iq(qx.flatten(), qy.flatten())
940    result = result.reshape(qx.shape)
941
942    import pylab
943    pylab.imshow(np.log(result + 0.001))
944    pylab.show()
945
946if __name__ == "__main__":
947    print("cylinder(0.1,0.1)=%g"%test_cylinder())
948    #magnetic_demo()
949    #test_product()
950    #test_structure_factor()
951    #print("rpa:", test_rpa())
952    #test_empty_distribution()
Note: See TracBrowser for help on using the repository browser.