source: sasmodels/sasmodels/sasview_model.py @ e2592f0

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

restrict magnetic parameters to visible shells. Refs #1188

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