source: sasmodels/sasmodels/sasview_model.py @ c952e59

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since c952e59 was c952e59, checked in by Torin Cooper-Bennun <torin.cooper-bennun@…>, 6 years ago

make SasviewModel?'s parsing of par units, bounds sane again

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