source: sasmodels/sasmodels/sasview_model.py @ c64a68e

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

do not release kernel in sasview_model wrapper; this avoids error for sasview on linux when dll handle is released

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