source: sasmodels/sasmodels/sasview_model.py @ 3221de0

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

restructure handling of opencl flags so it works with sasview

  • Property mode set to 100644
File size: 30.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#: 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            if not self._model_info.parameters.has_2d:
596                return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2))
597            else:
598                return self.calculate_Iq(qx, qy)
599
600        elif isinstance(qdist, np.ndarray):
601            # We have a simple 1D distribution of q-values
602            return self.calculate_Iq(qdist)
603
604        else:
605            raise TypeError("evalDistribution expects q or [qx, qy], not %r"
606                            % type(qdist))
607
608    def calc_composition_models(self, qx):
609        """
610        returns parts of the composition model or None if not a composition
611        model.
612        """
613        # TODO: have calculate_Iq return the intermediates.
614        #
615        # The current interface causes calculate_Iq() to be called twice,
616        # once to get the combined result and again to get the intermediate
617        # results.  This is necessary for now.
618        # Long term, the solution is to change the interface to calculate_Iq
619        # so that it returns a results object containing all the bits:
620        #     the A, B, C, ... of the composition model (and any subcomponents?)
621        #     the P and S of the product model,
622        #     the combined model before resolution smearing,
623        #     the sasmodel before sesans conversion,
624        #     the oriented 2D model used to fit oriented usans data,
625        #     the final I(q),
626        #     ...
627        #
628        # Have the model calculator add all of these blindly to the data
629        # tree, and update the graphs which contain them.  The fitter
630        # needs to be updated to use the I(q) value only, ignoring the rest.
631        #
632        # The simple fix of returning the existing intermediate results
633        # will not work for a couple of reasons: (1) another thread may
634        # sneak in to compute its own results before calc_composition_models
635        # is called, and (2) calculate_Iq is currently called three times:
636        # once with q, once with q values before qmin and once with q values
637        # after q max.  Both of these should be addressed before
638        # replacing this code.
639        composition = self._model_info.composition
640        if composition and composition[0] == 'product': # only P*S for now
641            with calculation_lock:
642                self._calculate_Iq(qx)
643                return self._intermediate_results
644        else:
645            return None
646
647    def calculate_Iq(self, qx, qy=None):
648        # type: (Sequence[float], Optional[Sequence[float]]) -> np.ndarray
649        """
650        Calculate Iq for one set of q with the current parameters.
651
652        If the model is 1D, use *q*.  If 2D, use *qx*, *qy*.
653
654        This should NOT be used for fitting since it copies the *q* vectors
655        to the card for each evaluation.
656        """
657        ## uncomment the following when trying to debug the uncoordinated calls
658        ## to calculate_Iq
659        #if calculation_lock.locked():
660        #    logger.info("calculation waiting for another thread to complete")
661        #    logger.info("\n".join(traceback.format_stack()))
662
663        with calculation_lock:
664            return self._calculate_Iq(qx, qy)
665
666    def _calculate_Iq(self, qx, qy=None):
667        if self._model is None:
668            self._model = core.build_model(self._model_info)
669        if qy is not None:
670            q_vectors = [np.asarray(qx), np.asarray(qy)]
671        else:
672            q_vectors = [np.asarray(qx)]
673        calculator = self._model.make_kernel(q_vectors)
674        parameters = self._model_info.parameters
675        pairs = [self._get_weights(p) for p in parameters.call_parameters]
676        #weights.plot_weights(self._model_info, pairs)
677        call_details, values, is_magnetic = make_kernel_args(calculator, pairs)
678        #call_details.show()
679        #print("pairs", pairs)
680        #for k, p in enumerate(self._model_info.parameters.call_parameters):
681        #    print(k, p.name, *pairs[k])
682        #print("params", self.params)
683        #print("values", values)
684        #print("is_mag", is_magnetic)
685        result = calculator(call_details, values, cutoff=self.cutoff,
686                            magnetic=is_magnetic)
687        #print("result", result)
688        self._intermediate_results = getattr(calculator, 'results', None)
689        calculator.release()
690        self._model.release()
691        return result
692
693    def calculate_ER(self):
694        # type: () -> float
695        """
696        Calculate the effective radius for P(q)*S(q)
697
698        :return: the value of the effective radius
699        """
700        if self._model_info.ER is None:
701            return 1.0
702        else:
703            value, weight = self._dispersion_mesh()
704            fv = self._model_info.ER(*value)
705            #print(values[0].shape, weights.shape, fv.shape)
706            return np.sum(weight * fv) / np.sum(weight)
707
708    def calculate_VR(self):
709        # type: () -> float
710        """
711        Calculate the volf ratio for P(q)*S(q)
712
713        :return: the value of the volf ratio
714        """
715        if self._model_info.VR is None:
716            return 1.0
717        else:
718            value, weight = self._dispersion_mesh()
719            whole, part = self._model_info.VR(*value)
720            return np.sum(weight * part) / np.sum(weight * whole)
721
722    def set_dispersion(self, parameter, dispersion):
723        # type: (str, weights.Dispersion) -> Dict[str, Any]
724        """
725        Set the dispersion object for a model parameter
726
727        :param parameter: name of the parameter [string]
728        :param dispersion: dispersion object of type Dispersion
729        """
730        if parameter in self.params:
731            # TODO: Store the disperser object directly in the model.
732            # The current method of relying on the sasview GUI to
733            # remember them is kind of funky.
734            # Note: can't seem to get disperser parameters from sasview
735            # (1) Could create a sasview model that has not yet been
736            # converted, assign the disperser to one of its polydisperse
737            # parameters, then retrieve the disperser parameters from the
738            # sasview model.
739            # (2) Could write a disperser parameter retriever in sasview.
740            # (3) Could modify sasview to use sasmodels.weights dispersers.
741            # For now, rely on the fact that the sasview only ever uses
742            # new dispersers in the set_dispersion call and create a new
743            # one instead of trying to assign parameters.
744            self.dispersion[parameter] = dispersion.get_pars()
745        else:
746            raise ValueError("%r is not a dispersity or orientation parameter")
747
748    def _dispersion_mesh(self):
749        # type: () -> List[Tuple[np.ndarray, np.ndarray]]
750        """
751        Create a mesh grid of dispersion parameters and weights.
752
753        Returns [p1,p2,...],w where pj is a vector of values for parameter j
754        and w is a vector containing the products for weights for each
755        parameter set in the vector.
756        """
757        pars = [self._get_weights(p)
758                for p in self._model_info.parameters.call_parameters
759                if p.type == 'volume']
760        return dispersion_mesh(self._model_info, pars)
761
762    def _get_weights(self, par):
763        # type: (Parameter) -> Tuple[np.ndarray, np.ndarray]
764        """
765        Return dispersion weights for parameter
766        """
767        if par.name not in self.params:
768            if par.name == self.multiplicity_info.control:
769                return self.multiplicity, [self.multiplicity], [1.0]
770            else:
771                # For hidden parameters use default values.  This sets
772                # scale=1 and background=0 for structure factors
773                default = self._model_info.parameters.defaults.get(par.name, np.NaN)
774                return default, [default], [1.0]
775        elif par.polydisperse:
776            value = self.params[par.name]
777            dis = self.dispersion[par.name]
778            if dis['type'] == 'array':
779                dispersity, weight = dis['values'], dis['weights']
780            else:
781                dispersity, weight = weights.get_weights(
782                    dis['type'], dis['npts'], dis['width'], dis['nsigmas'],
783                    value, par.limits, par.relative_pd)
784            return value, dispersity, weight
785        else:
786            value = self.params[par.name]
787            return value, [value], [1.0]
788
789def test_cylinder():
790    # type: () -> float
791    """
792    Test that the cylinder model runs, returning the value at [0.1,0.1].
793    """
794    Cylinder = _make_standard_model('cylinder')
795    cylinder = Cylinder()
796    return cylinder.evalDistribution([0.1, 0.1])
797
798def test_structure_factor():
799    # type: () -> float
800    """
801    Test that 2-D hardsphere model runs and doesn't produce NaN.
802    """
803    Model = _make_standard_model('hardsphere')
804    model = Model()
805    value2d = model.evalDistribution([0.1, 0.1])
806    value1d = model.evalDistribution(np.array([0.1*np.sqrt(2)]))
807    #print("hardsphere", value1d, value2d)
808    if np.isnan(value1d) or np.isnan(value2d):
809        raise ValueError("hardsphere returns nan")
810
811def test_product():
812    # type: () -> float
813    """
814    Test that 2-D hardsphere model runs and doesn't produce NaN.
815    """
816    S = _make_standard_model('hayter_msa')()
817    P = _make_standard_model('cylinder')()
818    model = MultiplicationModel(P, S)
819    value = model.evalDistribution([0.1, 0.1])
820    if np.isnan(value):
821        raise ValueError("cylinder*hatyer_msa returns null")
822
823def test_rpa():
824    # type: () -> float
825    """
826    Test that the 2-D RPA model runs
827    """
828    RPA = _make_standard_model('rpa')
829    rpa = RPA(3)
830    return rpa.evalDistribution([0.1, 0.1])
831
832def test_empty_distribution():
833    # type: () -> None
834    """
835    Make sure that sasmodels returns NaN when there are no polydispersity points
836    """
837    Cylinder = _make_standard_model('cylinder')
838    cylinder = Cylinder()
839    cylinder.setParam('radius', -1.0)
840    cylinder.setParam('background', 0.)
841    Iq = cylinder.evalDistribution(np.asarray([0.1]))
842    assert Iq[0] == 0., "empty distribution fails"
843
844def test_model_list():
845    # type: () -> None
846    """
847    Make sure that all models build as sasview models
848    """
849    from .exception import annotate_exception
850    for name in core.list_models():
851        try:
852            _make_standard_model(name)
853        except:
854            annotate_exception("when loading "+name)
855            raise
856
857def test_old_name():
858    # type: () -> None
859    """
860    Load and run cylinder model as sas-models-CylinderModel
861    """
862    if not SUPPORT_OLD_STYLE_PLUGINS:
863        return
864    try:
865        # if sasview is not on the path then don't try to test it
866        import sas
867    except ImportError:
868        return
869    load_standard_models()
870    from sas.models.CylinderModel import CylinderModel
871    CylinderModel().evalDistribution([0.1, 0.1])
872
873if __name__ == "__main__":
874    print("cylinder(0.1,0.1)=%g"%test_cylinder())
875    #test_product()
876    #test_structure_factor()
877    #print("rpa:", test_rpa())
878    #test_empty_distribution()
Note: See TracBrowser for help on using the repository browser.