Changeset e2da671 in sasmodels


Ignore:
Timestamp:
Mar 6, 2019 6:10:02 PM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master
Children:
f64b154
Parents:
4e96703 (diff), 9150036 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'beta_approx' into ticket-608-user-defined-weights

Files:
5 added
21 edited
6 moved

Legend:

Unmodified
Added
Removed

  • README.rst

    re30d645 r2a64722  
    1010is available. 
    1111 
    12 Example 
     12Install 
    1313------- 
     14 
     15The easiest way to use sasmodels is from `SasView <http://www.sasview.org/>`_. 
     16 
     17You can also install sasmodels as a standalone package in python. Use 
     18`miniconda <https://docs.conda.io/en/latest/miniconda.html>`_ 
     19or `anaconda <https://www.anaconda.com/>`_ 
     20to create a python environment with the sasmodels dependencies:: 
     21 
     22    $ conda create -n sasmodels -c conda-forge numpy scipy matplotlib pyopencl 
     23 
     24The option ``-n sasmodels`` names the environment sasmodels, and the option 
     25``-c conda-forge`` selects the conda-forge package channel because pyopencl 
     26is not part of the base anaconda distribution. 
     27 
     28Activate the environment and install sasmodels:: 
     29 
     30    $ conda activate sasmodels 
     31    (sasmodels) $ pip install sasmodels 
     32 
     33Install `bumps <https://github.com/bumps/bumps>`_ if you want to use it to fit 
     34your data:: 
     35 
     36    (sasmodels) $ pip install bumps 
     37 
     38Usage 
     39----- 
     40 
     41Check that the works:: 
     42 
     43    (sasmodels) $ python -m sasmodels.compare cylinder 
     44 
     45To show the orientation explorer:: 
     46 
     47    (sasmodels) $ python -m sasmodels.jitter 
     48 
     49Documentation is available online as part of the SasView 
     50`fitting perspective <http://www.sasview.org/docs/index.html>`_ 
     51as well as separate pages for 
     52`individual models <http://www.sasview.org/docs/user/sasgui/perspectives/fitting/models/index.html>`_. 
     53Programming details for sasmodels are available in the 
     54`developer documentation <http://www.sasview.org/docs/dev/dev.html>`_. 
     55 
     56 
     57Fitting Example 
     58--------------- 
    1459 
    1560The example directory contains a radial+tangential data set for an oriented 
    1661rod-like shape. 
    1762 
    18 The data is loaded by sas.dataloader from the sasview package, so sasview 
    19 is needed to run the example. 
     63To load the example data, you will need the SAS data loader from the sasview 
     64package. This is not yet available on PyPI, so you will need a copy of the 
     65SasView source code to run it.  Create a directory somewhere to hold the 
     66sasview and sasmodels source code, which we will refer to as $SOURCE. 
    2067 
    21 To run the example, you need sasview, sasmodels and bumps.  Assuming these 
    22 repositories are installed side by side, change to the sasmodels/example 
    23 directory and enter:: 
     68Use the following to install sasview, and the sasmodels examples:: 
    2469 
    25     PYTHONPATH=..:../../sasview/src ../../bumps/run.py fit.py \ 
    26         cylinder --preview 
     70    (sasmodels) $ cd $SOURCE 
     71    (sasmodels) $ conda install git 
     72    (sasmodels) $ git clone https://github.com/sasview/sasview.git 
     73    (sasmodels) $ git clone https://github.com/sasview/sasmodels.git 
    2774 
    28 See bumps documentation for instructions on running the fit.  With the 
    29 python packages installed, e.g., into a virtual environment, then the 
    30 python path need not be set, and the command would be:: 
     75Set the path to the sasview source on your python path within the sasmodels 
     76environment.  On Windows, this will be:: 
    3177 
    32     bumps fit.py cylinder --preview 
     78    (sasmodels)> set PYTHONPATH="$SOURCE\sasview\src" 
     79    (sasmodels)> cd $SOURCE/sasmodels/example 
     80    (sasmodels)> python -m bumps.cli fit.py cylinder --preview 
     81 
     82On Mac/Linux with the standard shell this will be:: 
     83 
     84    (sasmodels) $ export PYTHONPATH="$SOURCE/sasview/src" 
     85    (sasmodels) $ cd $SOURCE/sasmodels/example 
     86    (sasmodels) $ bumps fit.py cylinder --preview 
    3387 
    3488The fit.py model accepts up to two arguments.  The first argument is the 
    … …  
    3892both radial and tangential simultaneously, use the word "both". 
    3993 
    40 Notes 
    41 ----- 
    42  
    43 cylinder.c + cylinder.py is the cylinder model with renamed variables and 
    44 sld scaled by 1e6 so the numbers are nicer.  The model name is "cylinder" 
    45  
    46 lamellar.py is an example of a single file model with embedded C code. 
     94See `bumps documentation <https://bumps.readthedocs.io/>`_ for detailed 
     95instructions on running the fit. 
    4796 
    4897|TravisStatus|_ 

  • doc/guide/plugin.rst

    raa8c6e0 r9150036  
    272272structure factor to account for interactions between particles.  See 
    273273`Form_Factors`_ for more details. 
     274 
     275**model_info = ...** lets you define a model directly, for example, by 
     276loading and modifying existing models.  This is done implicitly by 
     277:func:`sasmodels.core.load_model_info`, which can create a mixture model 
     278from a pair of existing models.  For example:: 
     279 
     280    from sasmodels.core import load_model_info 
     281    model_info = load_model_info('sphere+cylinder') 
     282 
     283See :class:`sasmodels.modelinfo.ModelInfo` for details about the model 
     284attributes that are defined. 
    274285 
    275286Model Parameters 
    … …  
    894905             - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right) 
    895906 
    896         For small arguments , 
     907        For small arguments, 
    897908 
    898909        .. math:: 

  • example/multiscatfit.py

    r49d1f8b8 r2c4a190  
    1515 
    1616    # Show the model without fitting 
    17     PYTHONPATH=..:../explore:../../bumps:../../sasview/src python multiscatfit.py 
     17    PYTHONPATH=..:../../bumps:../../sasview/src python multiscatfit.py 
    1818 
    1919    # Run the fit 
    20     PYTHONPATH=..:../explore:../../bumps:../../sasview/src ../../bumps/run.py \ 
     20    PYTHONPATH=..:../../bumps:../../sasview/src ../../bumps/run.py \ 
    2121    multiscatfit.py --store=/tmp/t1 
    2222 
    … …  
    5555    ) 
    5656 
     57# Tie the model to the data 
     58M = Experiment(data=data, model=model) 
     59 
     60# Stack mulitple scattering on top of the existing resolution function. 
     61M.resolution = MultipleScattering(resolution=M.resolution, probability=0.) 
     62 
    5763# SET THE FITTING PARAMETERS 
    5864model.radius_polar.range(15, 3000) 
    … …  
    6571model.scale.range(0, 0.1) 
    6672 
    67 # Mulitple scattering probability parameter 
    68 # HACK: the probability is stuffed in as an extra parameter to the experiment. 
    69 probability = Parameter(name="probability", value=0.0) 
    70 probability.range(0.0, 0.9) 
     73# The multiple scattering probability parameter is in the resolution function 
     74# instead of the scattering function, so access it through M.resolution 
     75M.scattering_probability.range(0.0, 0.9) 
    7176 
    72 M = Experiment(data=data, model=model, extra_pars={'probability': probability}) 
    73  
    74 # Stack mulitple scattering on top of the existing resolution function. 
    75 # Because resolution functions in sasview don't have fitting parameters, 
    76 # we instead allow the multiple scattering calculator to take a function 
    77 # instead of a probability.  This function returns the current value of 
    78 # the parameter. ** THIS IS TEMPORARY ** when multiple scattering is 
    79 # properly integrated into sasmodels and sasview, its fittable parameter 
    80 # will be treated like the model parameters. 
    81 M.resolution = MultipleScattering(resolution=M.resolution, 
    82                                   probability=lambda: probability.value, 
    83                                   ) 
    84 M._kernel_inputs = M.resolution.q_calc 
     77# Let bumps know that we are fitting this experiment 
    8578problem = FitProblem(M) 
    8679 

  • explore/beta/sasfit_compare.py

    r2a12351b r119073a  
    505505    } 
    506506 
    507     Q, IQ = load_sasfit(data_file('richard_test.txt')) 
    508     Q, IQSD = load_sasfit(data_file('richard_test2.txt')) 
    509     Q, IQBD = load_sasfit(data_file('richard_test3.txt')) 
     507    Q, IQ = load_sasfit(data_file('sasfit_sphere_schulz_IQD.txt')) 
     508    Q, IQSD = load_sasfit(data_file('sasfit_sphere_schulz_IQSD.txt')) 
     509    Q, IQBD = load_sasfit(data_file('sasfit_sphere_schulz_IQBD.txt')) 
    510510    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 
    511511    actual = sphere_r(Q, norm="sasfit", **pars) 
    … …  
    526526    } 
    527527 
    528     Q, IQ = load_sasfit(data_file('richard_test4.txt')) 
    529     Q, IQSD = load_sasfit(data_file('richard_test5.txt')) 
    530     Q, IQBD = load_sasfit(data_file('richard_test6.txt')) 
     528    Q, IQ = load_sasfit(data_file('sasfit_ellipsoid_shulz_IQD.txt')) 
     529    Q, IQSD = load_sasfit(data_file('sasfit_ellipsoid_shulz_IQSD.txt')) 
     530    Q, IQBD = load_sasfit(data_file('sasfit_ellipsoid_shulz_IQBD.txt')) 
    531531    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 
    532532    actual = ellipsoid_pe(Q, norm="sasfit", **pars) 

  • explore/precision.py

    raa8c6e0 rcd28947  
    207207    return model_info 
    208208 
    209 # Hack to allow second parameter A in two parameter functions 
     209# Hack to allow second parameter A in the gammainc and gammaincc functions. 
     210# Create a 2-D variant of the precision test if we need to handle other two 
     211# parameter functions. 
    210212A = 1 
    211213def parse_extra_pars(): 
     214    """ 
     215    Parse the command line looking for the second parameter "A=..." for the 
     216    gammainc/gammaincc functions. 
     217    """ 
    212218    global A 
    213219 
    … …  
    333339) 
    334340add_function( 
     341    # Note: "a" is given as A=... on the command line via parse_extra_pars 
    335342    name="gammainc(x)", 
    336343    mp_function=lambda x, a=A: mp.gammainc(a, a=0, b=x)/mp.gamma(a), 
    … …  
    339346) 
    340347add_function( 
     348    # Note: "a" is given as A=... on the command line via parse_extra_pars 
    341349    name="gammaincc(x)", 
    342350    mp_function=lambda x, a=A: mp.gammainc(a, a=x, b=mp.inf)/mp.gamma(a), 

  • sasmodels/__init__.py

    ra1ec908 r37f38ff  
    1414defining new models. 
    1515""" 
    16 __version__ = "0.98" 
     16__version__ = "0.99" 
    1717 
    1818def data_files(): 

  • sasmodels/bumps_model.py

    r49d1f8b8 r2c4a190  
    3535    # when bumps is not on the path. 
    3636    from bumps.names import Parameter # type: ignore 
     37    from bumps.parameter import Reference # type: ignore 
    3738except ImportError: 
    3839    pass 
    … …  
    139140    def __init__(self, data, model, cutoff=1e-5, name=None, extra_pars=None): 
    140141        # type: (Data, Model, float) -> None 
     142        # Allow resolution function to define fittable parameters.  We do this 
     143        # by creating reference parameters within the resolution object rather 
     144        # than modifying the object itself to use bumps parameters.  We need 
     145        # to reset the parameters each time the object has changed.  These 
     146        # additional parameters need to be returned from the fitting engine. 
     147        # To make them available to the user, they are added as top-level 
     148        # attributes to the experiment object.  The only change to the 
     149        # resolution function is that it needs an optional 'fittable' attribute 
     150        # which maps the internal name to the user visible name for the 
     151        # for the parameter. 
     152        self._resolution = None 
     153        self._resolution_pars = {} 
    141154        # remember inputs so we can inspect from outside 
    142155        self.name = data.filename if name is None else name 
    … …  
    145158        self._interpret_data(data, model.sasmodel) 
    146159        self._cache = {} 
     160        # CRUFT: no longer need extra parameters 
     161        # Multiple scattering probability is now retrieved directly from the 
     162        # multiple scattering resolution function. 
    147163        self.extra_pars = extra_pars 
    148164 
    … …  
    162178        return len(self.Iq) 
    163179 
     180    @property 
     181    def resolution(self): 
     182        return self._resolution 
     183 
     184    @resolution.setter 
     185    def resolution(self, value): 
     186        self._resolution = value 
     187 
     188        # Remove old resolution fitting parameters from experiment 
     189        for name in self._resolution_pars: 
     190            delattr(self, name) 
     191 
     192        # Create new resolution fitting parameters 
     193        res_pars = getattr(self._resolution, 'fittable', {}) 
     194        self._resolution_pars = { 
     195            name: Reference(self._resolution, refname, name=name) 
     196            for refname, name in res_pars.items() 
     197        } 
     198 
     199        # Add new resolution fitting parameters as experiment attributes 
     200        for name, ref in self._resolution_pars.items(): 
     201            setattr(self, name, ref) 
     202 
    164203    def parameters(self): 
    165204        # type: () -> Dict[str, Parameter] 
    … …  
    168207        """ 
    169208        pars = self.model.parameters() 
    170         if self.extra_pars: 
     209        if self.extra_pars is not None: 
    171210            pars.update(self.extra_pars) 
     211        pars.update(self._resolution_pars) 
    172212        return pars 
    173213 

  • sasmodels/compare.py

    r4e96703 re2da671  
    11531153        'rel_err'   : True, 
    11541154        'explore'   : False, 
    1155         'use_demo'  : True, 
     1155        'use_demo'  : False, 
    11561156        'zero'      : False, 
    11571157        'html'      : False, 

  • sasmodels/direct_model.py

    r5024a56 r9150036  
    224224            else: 
    225225                Iq, dIq = None, None 
    226             #self._theory = np.zeros_like(q) 
    227             q_vectors = [res.q_calc] 
    228226        elif self.data_type == 'Iqxy': 
    229227            #if not model.info.parameters.has_2d: 
    … …  
    242240            res = resolution2d.Pinhole2D(data=data, index=index, 
    243241                                         nsigma=3.0, accuracy=accuracy) 
    244             #self._theory = np.zeros_like(self.Iq) 
    245             q_vectors = res.q_calc 
    246242        elif self.data_type == 'Iq': 
    247243            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
    … …  
    268264            else: 
    269265                res = resolution.Perfect1D(data.x[index]) 
    270  
    271             #self._theory = np.zeros_like(self.Iq) 
    272             q_vectors = [res.q_calc] 
    273266        elif self.data_type == 'Iq-oriented': 
    274267            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
    … …  
    286279                                      qx_width=data.dxw[index], 
    287280                                      qy_width=data.dxl[index]) 
    288             q_vectors = res.q_calc 
    289281        else: 
    290282            raise ValueError("Unknown data type") # never gets here 
    … …  
    292284        # Remember function inputs so we can delay loading the function and 
    293285        # so we can save/restore state 
    294         self._kernel_inputs = q_vectors 
    295286        self._kernel = None 
    296287        self.Iq, self.dIq, self.index = Iq, dIq, index 
    … …  
    329320        # type: (ParameterSet, float) -> np.ndarray 
    330321        if self._kernel is None: 
    331             self._kernel = self._model.make_kernel(self._kernel_inputs) 
     322            # TODO: change interfaces so that resolution returns kernel inputs 
     323            # Maybe have resolution always return a tuple, or maybe have 
     324            # make_kernel accept either an ndarray or a pair of ndarrays. 
     325            kernel_inputs = self.resolution.q_calc 
     326            if isinstance(kernel_inputs, np.ndarray): 
     327                kernel_inputs = (kernel_inputs,) 
     328            self._kernel = self._model.make_kernel(kernel_inputs) 
    332329 
    333330        # Need to pull background out of resolution for multiple scattering 
    334         background = pars.get('background', DEFAULT_BACKGROUND) 
     331        default_background = self._model.info.parameters.common_parameters[1].default 
     332        background = pars.get('background', default_background) 
    335333        pars = pars.copy() 
    336334        pars['background'] = 0. 

  • sasmodels/generate.py

    r39a06c9 ra8a1d48  
    703703    """ 
    704704    for code in source: 
    705         m = _FQ_PATTERN.search(code) 
    706         if m is not None: 
     705        if _FQ_PATTERN.search(code) is not None: 
    707706            return True 
    708707    return False 
    … …  
    712711    # type: (List[str]) -> bool 
    713712    """ 
    714     Return True if C source defines "void Fq(". 
     713    Return True if C source defines "double shell_volume(". 
    715714    """ 
    716715    for code in source: 
    717         m = _SHELL_VOLUME_PATTERN.search(code) 
    718         if m is not None: 
     716        if _SHELL_VOLUME_PATTERN.search(code) is not None: 
    719717            return True 
    720718    return False 
    … …  
    10081006        pars = model_info.parameters.kernel_parameters 
    10091007    else: 
    1010         pars = model_info.parameters.COMMON + model_info.parameters.kernel_parameters 
     1008        pars = (model_info.parameters.common_parameters 
     1009                + model_info.parameters.kernel_parameters) 
    10111010    partable = make_partable(pars) 
    10121011    subst = dict(id=model_info.id.replace('_', '-'), 

  • sasmodels/jitter.py

    r1198f90 r7d97437  
    1515    pass 
    1616 
     17import matplotlib as mpl 
    1718import matplotlib.pyplot as plt 
    1819from matplotlib.widgets import Slider 
    … …  
    746747        pass 
    747748 
    748     axcolor = 'lightgoldenrodyellow' 
     749    # CRUFT: use axisbg instead of facecolor for matplotlib<2 
     750    facecolor_prop = 'facecolor' if mpl.__version__ > '2' else 'axisbg' 
     751    props = {facecolor_prop: 'lightgoldenrodyellow'} 
    749752 
    750753    ## add control widgets to plot 
    751     axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    752     axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
    753     axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor) 
     754    axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], **props) 
     755    axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], **props) 
     756    axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], **props) 
    754757    stheta = Slider(axes_theta, 'Theta', -90, 90, valinit=theta) 
    755758    sphi = Slider(axes_phi, 'Phi', -180, 180, valinit=phi) 
    756759    spsi = Slider(axes_psi, 'Psi', -180, 180, valinit=psi) 
    757760 
    758     axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    759     axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
    760     axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) 
     761    axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], **props) 
     762    axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], **props) 
     763    axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], **props) 
    761764    # Note: using ridiculous definition of rectangle distribution, whose width 
    762765    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 

  • sasmodels/kernel.py

    re44432d rcd28947  
    133133        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
    134134        total_weight = self.result[nout*self.q_input.nq + 0] 
     135        # Note: total_weight = sum(weight > cutoff), with cutoff >= 0, so it 
     136        # is okay to test directly against zero.  If weight is zero then I(q), 
     137        # etc. must also be zero. 
    135138        if total_weight == 0.: 
    136139            total_weight = 1. 

  • sasmodels/kernelcl.py

    rf872fd1 r8064d5e  
    5858import time 
    5959 
     60try: 
     61    from time import perf_counter as clock 
     62except ImportError: # CRUFT: python < 3.3 
     63    import sys 
     64    if sys.platform.count("darwin") > 0: 
     65        from time import time as clock 
     66    else: 
     67        from time import clock 
     68 
    6069import numpy as np  # type: ignore 
    61  
    6270 
    6371# Attempt to setup opencl. This may fail if the pyopencl package is not 
    … …  
    611619        #Call kernel and retrieve results 
    612620        wait_for = None 
    613         last_nap = time.clock() 
     621        last_nap = clock() 
    614622        step = 1000000//self.q_input.nq + 1 
    615623        for start in range(0, call_details.num_eval, step): 
    … …  
    622630                # Allow other processes to run 
    623631                wait_for[0].wait() 
    624                 current_time = time.clock() 
     632                current_time = clock() 
    625633                if current_time - last_nap > 0.5: 
    626634                    time.sleep(0.001) 

  • sasmodels/modelinfo.py

    r39a06c9 rc1799d3  
    404404      parameters counted as n individual parameters p1, p2, ... 
    405405 
     406    * *common_parameters* is the list of common parameters, with a unique 
     407      copy for each model so that structure factors can have a default 
     408      background of 0.0. 
     409 
    406410    * *call_parameters* is the complete list of parameters to the kernel, 
    407411      including scale and background, with vector parameters recorded as 
    … …  
    416420    parameters don't use vector notation, and instead use p1, p2, ... 
    417421    """ 
    418     # scale and background are implicit parameters 
    419     COMMON = [Parameter(*p) for p in COMMON_PARAMETERS] 
    420  
    421422    def __init__(self, parameters): 
    422423        # type: (List[Parameter]) -> None 
     424 
     425        # scale and background are implicit parameters 
     426        # Need them to be unique to each model in case they have different 
     427        # properties, such as default=0.0 for structure factor backgrounds. 
     428        self.common_parameters = [Parameter(*p) for p in COMMON_PARAMETERS] 
     429 
    423430        self.kernel_parameters = parameters 
    424431        self._set_vector_lengths() 
    … …  
    468475                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
    469476        self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse) 
     477 
     478    def set_zero_background(self): 
     479        """ 
     480        Set the default background to zero for this model.  This is done for 
     481        structure factor models. 
     482        """ 
     483        # type: () -> None 
     484        # Make sure background is the second common parameter. 
     485        assert self.common_parameters[1].id == "background" 
     486        self.common_parameters[1].default = 0.0 
     487        self.defaults = self._get_defaults() 
    470488 
    471489    def check_angles(self): 
    … …  
    567585    def _get_call_parameters(self): 
    568586        # type: () -> List[Parameter] 
    569         full_list = self.COMMON[:] 
     587        full_list = self.common_parameters[:] 
    570588        for p in self.kernel_parameters: 
    571589            if p.length == 1: 
    … …  
    670688 
    671689        # Gather the user parameters in order 
    672         result = control + self.COMMON 
     690        result = control + self.common_parameters 
    673691        for p in self.kernel_parameters: 
    674692            if not is2d and p.type in ('orientation', 'magnetic'): 
    … …  
    770788 
    771789    info = ModelInfo() 
     790 
     791    # Build the parameter table 
    772792    #print("make parameter table", kernel_module.parameters) 
    773793    parameters = make_parameter_table(getattr(kernel_module, 'parameters', [])) 
     794 
     795    # background defaults to zero for structure factor models 
     796    structure_factor = getattr(kernel_module, 'structure_factor', False) 
     797    if structure_factor: 
     798        parameters.set_zero_background() 
     799 
     800    # TODO: remove demo parameters 
     801    # The plots in the docs are generated from the model default values. 
     802    # Sascomp set parameters from the command line, and so doesn't need 
     803    # demo values for testing. 
    774804    demo = expand_pars(parameters, getattr(kernel_module, 'demo', None)) 
     805 
    775806    filename = abspath(kernel_module.__file__).replace('.pyc', '.py') 
    776807    kernel_id = splitext(basename(filename))[0] 

  • sasmodels/models/hardsphere.py

    r304c775 rc1799d3  
    162162    return pars 
    163163 
    164 demo = dict(radius_effective=200, volfraction=0.2, 
    165             radius_effective_pd=0.1, radius_effective_pd_n=40) 
    166164# Q=0.001 is in the Taylor series, low Q part, so add Q=0.1, 
    167165# assuming double precision sasview is correct 

  • sasmodels/models/pearl_necklace.c

    r99658f6 r9b5fd42  
    4040    const double si = sas_sinx_x(q*A_s); 
    4141    const double omsi = 1.0 - si; 
    42     const double pow_si = pow(si, num_pearls); 
     42    const double pow_si = pown(si, num_pearls); 
    4343 
    4444    // form factor for num_pearls 
    … …  
    8181radius_from_volume(double radius, double edge_sep, double thick_string, double fp_num_pearls) 
    8282{ 
    83     const int num_pearls = (int) fp_num_pearls +0.5; 
    8483    const double vol_tot = form_volume(radius, edge_sep, thick_string, fp_num_pearls); 
    8584    return cbrt(vol_tot/M_4PI_3); 

  • sasmodels/multiscat.py

    rb3703f5 r2c4a190  
    342342 
    343343    *probability* is related to the expected number of scattering 
    344     events in the sample $\lambda$ as $p = 1 = e^{-\lambda}$.  As a 
    345     hack to allow probability to be a fitted parameter, the "value" 
    346     can be a function that takes no parameters and returns the current 
    347     value of the probability.  *coverage* determines how many scattering 
    348     steps to consider.  The default is 0.99, which sets $n$ such that 
    349     $1 \ldots n$ covers 99% of the Poisson probability mass function. 
     344    events in the sample $\lambda$ as $p = 1 - e^{-\lambda}$. 
     345    *coverage* determines how many scattering steps to consider.  The 
     346    default is 0.99, which sets $n$ such that $1 \ldots n$ covers 99% 
     347    of the Poisson probability mass function. 
    350348 
    351349    *is2d* is True then 2D scattering is used, otherwise it accepts 
    … …  
    399397        self.qmin = qmin 
    400398        self.nq = nq 
    401         self.probability = probability 
     399        self.probability = 0. if probability is None else probability 
    402400        self.coverage = coverage 
    403401        self.is2d = is2d 
    … …  
    456454        self.Iqxy = None # type: np.ndarray 
    457455 
     456        # Label probability as a fittable parameter, and give its external name 
     457        # Note that the external name must be a valid python identifier, since 
     458        # is will be set as an experiment attribute. 
     459        self.fittable = {'probability': 'scattering_probability'} 
     460 
    458461    def apply(self, theory): 
    459462        if self.is2d: 
    … …  
    463466        Iq_calc = Iq_calc.reshape(self.nq, self.nq) 
    464467 
     468        # CRUFT: don't need probability as a function anymore 
    465469        probability = self.probability() if callable(self.probability) else self.probability 
    466470        coverage = self.coverage 

  • sasmodels/resolution.py

    r9e7837a re2592f0  
    445445    q = np.sort(q) 
    446446    if q_min + 2*MINIMUM_RESOLUTION < q[0]: 
    447         n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15 
     447        n_low = int(np.ceil((q[0]-q_min) / (q[1]-q[0]))) if q[1] > q[0] else 15 
    448448        q_low = np.linspace(q_min, q[0], n_low+1)[:-1] 
    449449    else: 
    450450        q_low = [] 
    451451    if q_max - 2*MINIMUM_RESOLUTION > q[-1]: 
    452         n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1] > q[-2] else 15 
     452        n_high = int(np.ceil((q_max-q[-1]) / (q[-1]-q[-2]))) if q[-1] > q[-2] else 15 
    453453        q_high = np.linspace(q[-1], q_max, n_high+1)[1:] 
    454454    else: 
    … …  
    499499            q_min = q[0]*MINIMUM_ABSOLUTE_Q 
    500500        n_low = log_delta_q * (log(q[0])-log(q_min)) 
    501         q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1] 
     501        q_low = np.logspace(log10(q_min), log10(q[0]), int(np.ceil(n_low))+1)[:-1] 
    502502    else: 
    503503        q_low = [] 
    504504    if q_max > q[-1]: 
    505505        n_high = log_delta_q * (log(q_max)-log(q[-1])) 
    506         q_high = np.logspace(log10(q[-1]), log10(q_max), np.ceil(n_high)+1)[1:] 
     506        q_high = np.logspace(log10(q[-1]), log10(q_max), int(np.ceil(n_high))+1)[1:] 
    507507    else: 
    508508        q_high = [] 

  • sasmodels/sasview_model.py

    r4e96703 re2da671  
    2525from . import core 
    2626from . import custom 
     27from . import kernelcl 
    2728from . import product 
    2829from . import generate 
    … …  
    3031from . import modelinfo 
    3132from .details import make_kernel_args, dispersion_mesh 
     33from .kernelcl import reset_environment 
    3234 
    3335# Hack: load in any custom distributions 
    … …  
    7375#: has changed since we last reloaded. 
    7476_CACHED_MODULE = {}  # type: Dict[str, "module"] 
     77 
     78def reset_environment(): 
     79    # type: () -> None 
     80    """ 
     81    Clear the compute engine context so that the GUI can change devices. 
     82 
     83    This removes all compiled kernels, even those that are active on fit 
     84    pages, but they will be restored the next time they are needed. 
     85    """ 
     86    kernelcl.reset_environment() 
     87    for model in MODELS.values(): 
     88        model._model = None 
    7589 
    7690def find_model(modelname): 
    … …  
    387401            hidden.add('scale') 
    388402            hidden.add('background') 
    389             self._model_info.parameters.defaults['background'] = 0. 
    390403 
    391404        # Update the parameter lists to exclude any hidden parameters 
    … …  
    700713            return self._calculate_Iq(qx, qy) 
    701714 
    702     def _calculate_Iq(self, qx, qy=None, Fq=False, effective_radius_type=1): 
     715    def _calculate_Iq(self, qx, qy=None): 
    703716        if self._model is None: 
    704             self._model = core.build_model(self._model_info) 
     717            # Only need one copy of the compiled kernel regardless of how many 
     718            # times it is used, so store it in the class.  Also, to reset the 
     719            # compute engine, need to clear out all existing compiled kernels, 
     720            # which is much easier to do if we store them in the class. 
     721            self.__class__._model = core.build_model(self._model_info) 
    705722        if qy is not None: 
    706723            q_vectors = [np.asarray(qx), np.asarray(qy)] 
    … …  
    720737        #print("values", values) 
    721738        #print("is_mag", is_magnetic) 
    722         if Fq: 
    723             result = calculator.Fq(call_details, values, cutoff=self.cutoff, 
    724                                    magnetic=is_magnetic, 
    725                                    effective_radius_type=effective_radius_type) 
    726739        result = calculator(call_details, values, cutoff=self.cutoff, 
    727740                            magnetic=is_magnetic) 
    … …  
    741754        Calculate the effective radius for P(q)*S(q) 
    742755 
     756        *mode* is the R_eff type, which defaults to 1 to match the ER 
     757        calculation for sasview models from version 3.x. 
     758 
    743759        :return: the value of the effective radius 
    744760        """ 
    745         Fq = self._calculate_Iq([0.1], True, mode) 
    746         return Fq[2] 
     761        # ER and VR are only needed for old multiplication models, based on 
     762        # sas.sascalc.fit.MultiplicationModel.  Fail for now.  If we want to 
     763        # continue supporting them then add some test cases so that the code 
     764        # is exercised.  We can access ER/VR using the kernel Fq function by 
     765        # extending _calculate_Iq so that it calls: 
     766        #    if er_mode > 0: 
     767        #        res = calculator.Fq(call_details, values, cutoff=self.cutoff, 
     768        #                            magnetic=False, effective_radius_type=mode) 
     769        #        R_eff, form_shell_ratio = res[2], res[4] 
     770        #        return R_eff, form_shell_ratio 
     771        # Then use the following in calculate_ER: 
     772        #    ER, VR = self._calculate_Iq(q=[0.1], er_mode=mode) 
     773        #    return ER 
     774        # Similarly, for calculate_VR: 
     775        #    ER, VR = self._calculate_Iq(q=[0.1], er_mode=1) 
     776        #    return VR 
     777        # Obviously a combined calculate_ER_VR method would be better, but 
     778        # we only need them to support very old models, so ignore the 2x 
     779        # performance hit. 
     780        raise NotImplementedError("ER function is no longer available.") 
    747781 
    748782    def calculate_VR(self): 
    … …  
    753787        :return: the value of the form:shell volume ratio 
    754788        """ 
    755         Fq = self._calculate_Iq([0.1], True, mode) 
    756         return Fq[4] 
     789        # See comments in calculate_ER. 
     790        raise NotImplementedError("VR function is no longer available.") 
    757791 
    758792    def set_dispersion(self, parameter, dispersion): 
    … …  
    919953    CylinderModel().evalDistribution([0.1, 0.1]) 
    920954 
     955def test_structure_factor_background(): 
     956    # type: () -> None 
     957    """ 
     958    Check that sasview model and direct model match, with background=0. 
     959    """ 
     960    from .data import empty_data1D 
     961    from .core import load_model_info, build_model 
     962    from .direct_model import DirectModel 
     963 
     964    model_name = "hardsphere" 
     965    q = [0.0] 
     966 
     967    sasview_model = _make_standard_model(model_name)() 
     968    sasview_value = sasview_model.evalDistribution(np.array(q))[0] 
     969 
     970    data = empty_data1D(q) 
     971    model_info = load_model_info(model_name) 
     972    model = build_model(model_info) 
     973    direct_model = DirectModel(data, model) 
     974    direct_value_zero_background = direct_model(background=0.0) 
     975 
     976    assert sasview_value == direct_value_zero_background 
     977 
     978    # Additionally check that direct value background defaults to zero 
     979    direct_value_default = direct_model() 
     980    assert sasview_value == direct_value_default 
     981 
     982 
    921983def magnetic_demo(): 
    922984    Model = _make_standard_model('sphere') 
    … …  
    9391001    #print("rpa:", test_rpa()) 
    9401002    #test_empty_distribution() 
     1003    #test_structure_factor_background() 

  • sasmodels/weights.py

    rf41027b re2da671  
    2323    default = dict(npts=35, width=0, nsigmas=3) 
    2424    def __init__(self, npts=None, width=None, nsigmas=None): 
    25         self.npts = self.default['npts'] if npts is None else npts 
     25        self.npts = self.default['npts'] if npts is None else int(npts) 
    2626        self.width = self.default['width'] if width is None else width 
    2727        self.nsigmas = self.default['nsigmas'] if nsigmas is None else nsigmas 

  • doc/guide/pd/polydispersity.rst

    rd089a00 ra5cb9bc  
    1111-------------------------------------------- 
    1212 
    13 For some models we can calculate the average intensity for a population of  
    14 particles that possess size and/or orientational (ie, angular) distributions.  
    15 In SasView we call the former *polydispersity* but use the parameter *PD* to  
    16 parameterise both. In other words, the meaning of *PD* in a model depends on  
     13For some models we can calculate the average intensity for a population of 
     14particles that possess size and/or orientational (ie, angular) distributions. 
     15In SasView we call the former *polydispersity* but use the parameter *PD* to 
     16parameterise both. In other words, the meaning of *PD* in a model depends on 
    1717the actual parameter it is being applied too. 
    1818 
    19 The resultant intensity is then normalized by the average particle volume such  
     19The resultant intensity is then normalized by the average particle volume such 
    2020that 
    2121 
    … …  
    2424  P(q) = \text{scale} \langle F^* F \rangle / V + \text{background} 
    2525 
    26 where $F$ is the scattering amplitude and $\langle\cdot\rangle$ denotes an  
     26where $F$ is the scattering amplitude and $\langle\cdot\rangle$ denotes an 
    2727average over the distribution $f(x; \bar x, \sigma)$, giving 
    2828 
    2929.. math:: 
    3030 
    31   P(q) = \frac{\text{scale}}{V} \int_\mathbb{R}  
     31  P(q) = \frac{\text{scale}}{V} \int_\mathbb{R} 
    3232  f(x; \bar x, \sigma) F^2(q, x)\, dx + \text{background} 
    3333 
    3434Each distribution is characterized by a center value $\bar x$ or 
    3535$x_\text{med}$, a width parameter $\sigma$ (note this is *not necessarily* 
    36 the standard deviation, so read the description of the distribution carefully),  
    37 the number of sigmas $N_\sigma$ to include from the tails of the distribution,  
    38 and the number of points used to compute the average. The center of the  
    39 distribution is set by the value of the model parameter. 
    40  
    41 The distribution width applied to *volume* (ie, shape-describing) parameters  
    42 is relative to the center value such that $\sigma = \mathrm{PD} \cdot \bar x$.  
    43 However, the distribution width applied to *orientation* parameters is just  
    44 $\sigma = \mathrm{PD}$. 
     36the standard deviation, so read the description carefully), the number of 
     37sigmas $N_\sigma$ to include from the tails of the distribution, and the 
     38number of points used to compute the average. The center of the distribution 
     39is set by the value of the model parameter. The meaning of a polydispersity 
     40parameter *PD* (not to be confused with a molecular weight distributions 
     41in polymer science) in a model depends on the type of parameter it is being 
     42applied too. 
     43 
     44The distribution width applied to *volume* (ie, shape-describing) parameters 
     45is relative to the center value such that $\sigma = \mathrm{PD} \cdot \bar x$. 
     46However, the distribution width applied to *orientation* (ie, angle-describing) 
     47parameters is just $\sigma = \mathrm{PD}$. 
    4548 
    4649$N_\sigma$ determines how far into the tails to evaluate the distribution, 
    … …  
    5255 
    5356Users should note that the averaging computation is very intensive. Applying 
    54 polydispersion and/or orientational distributions to multiple parameters at  
    55 the same time, or increasing the number of points in the distribution, will  
    56 require patience! However, the calculations are generally more robust with  
     57polydispersion and/or orientational distributions to multiple parameters at 
     58the same time, or increasing the number of points in the distribution, will 
     59require patience! However, the calculations are generally more robust with 
    5760more data points or more angles. 
    5861 
    … …  
    6669*  *Schulz Distribution* 
    6770*  *Array Distribution* 
     71*  *User-defined Distributions* 
    6872 
    6973These are all implemented as *number-average* distributions. 
    7074 
    71 Additional distributions are under consideration. 
    7275 
    7376**Beware: when the Polydispersity & Orientational Distribution panel in SasView is** 
    … …  
    7578**This may not be suitable. See Suggested Applications below.** 
    7679 
    77 .. note:: In 2009 IUPAC decided to introduce the new term 'dispersity' to replace  
    78            the term 'polydispersity' (see `Pure Appl. Chem., (2009), 81(2),  
    79            351-353 <http://media.iupac.org/publications/pac/2009/pdf/8102x0351.pdf>`_  
    80            in order to make the terminology describing distributions of chemical  
    81            properties unambiguous. However, these terms are unrelated to the  
    82            proportional size distributions and orientational distributions used in  
     80.. note:: In 2009 IUPAC decided to introduce the new term 'dispersity' to replace 
     81           the term 'polydispersity' (see `Pure Appl. Chem., (2009), 81(2), 
     82           351-353 <http://media.iupac.org/publications/pac/2009/pdf/8102x0351.pdf>`_ 
     83           in order to make the terminology describing distributions of chemical 
     84           properties unambiguous. However, these terms are unrelated to the 
     85           proportional size distributions and orientational distributions used in 
    8386           SasView models. 
    8487 
    … …  
    9295or angular orientations, consider using the Gaussian or Boltzmann distributions. 
    9396 
    94 If applying polydispersion to parameters describing angles, use the Uniform  
    95 distribution. Beware of using distributions that are always positive (eg, the  
     97If applying polydispersion to parameters describing angles, use the Uniform 
     98distribution. Beware of using distributions that are always positive (eg, the 
    9699Lognormal) because angles can be negative! 
    97100 
    98 The array distribution allows a user-defined distribution to be applied. 
     101The array distribution provides a very simple means of implementing a user- 
     102defined distribution, but without any fittable parameters. Greater flexibility 
     103is conferred by the user-defined distribution. 
    99104 
    100105.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
    … …  
    334339.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
    335340 
     341User-defined Distributions 
     342^^^^^^^^^^^^^^^^^^^^^^^^^^ 
     343 
     344You can also define your own distribution by creating a python file defining a 
     345*Distribution* object with a *_weights* method.  The *_weights* method takes 
     346*center*, *sigma*, *lb* and *ub* as arguments, and can access *self.npts* 
     347and *self.nsigmas* from the distribution.  They are interpreted as follows: 
     348 
     349* *center* the value of the shape parameter (for size dispersity) or zero 
     350  if it is an angular dispersity.  This parameter may be fitted. 
     351 
     352* *sigma* the width of the distribution, which is the polydispersity parameter 
     353  times the center for size dispersity, or the polydispersity parameter alone 
     354  for angular dispersity.  This parameter may be fitted. 
     355 
     356* *lb*, *ub* are the parameter limits (lower & upper bounds) given in the model 
     357  definition file.  For example, a radius parameter has *lb* equal to zero.  A 
     358  volume fraction parameter would have *lb* equal to zero and *ub* equal to one. 
     359 
     360* *self.nsigmas* the distance to go into the tails when evaluating the 
     361  distribution.  For a two parameter distribution, this value could be 
     362  co-opted to use for the second parameter, though it will not be available 
     363  for fitting. 
     364 
     365* *self.npts* the number of points to use when evaluating the distribution. 
     366  The user will adjust this to trade calculation time for accuracy, but the 
     367  distribution code is free to return more or fewer, or use it for the third 
     368  parameter in a three parameter distribution. 
     369 
     370As an example, the code following wraps the Laplace distribution from scipy stats:: 
     371 
     372    import numpy as np 
     373    from scipy.stats import laplace 
     374 
     375    from sasmodels import weights 
     376 
     377    class Dispersion(weights.Dispersion): 
     378        r""" 
     379        Laplace distribution 
     380 
     381        .. math:: 
     382 
     383            w(x) = e^{-\sigma |x - \mu|} 
     384        """ 
     385        type = "laplace" 
     386        default = dict(npts=35, width=0, nsigmas=3)  # default values 
     387        def _weights(self, center, sigma, lb, ub): 
     388            x = self._linspace(center, sigma, lb, ub) 
     389            wx = laplace.pdf(x, center, sigma) 
     390            return x, wx 
     391 
     392You can plot the weights for a given value and width using the following:: 
     393 
     394    from numpy import inf 
     395    from matplotlib import pyplot as plt 
     396    from sasmodels import weights 
     397 
     398    # reload the user-defined weights 
     399    weights.load_weights() 
     400    x, wx = weights.get_weights('laplace', n=35, width=0.1, nsigmas=3, value=50, 
     401                                limits=[0, inf], relative=True) 
     402 
     403    # plot the weights 
     404    plt.interactive(True) 
     405    plt.plot(x, wx, 'x') 
     406 
     407The *self.nsigmas* and *self.npts* parameters are normally used to control 
     408the accuracy of the distribution integral. The *self._linspace* function 
     409uses them to define the *x* values (along with the *center*, *sigma*, 
     410*lb*, and *ub* which are passed as parameters).  If you repurpose npts or 
     411nsigmas you will need to generate your own *x*.  Be sure to honour the 
     412limits *lb* and *ub*, for example to disallow a negative radius or constrain 
     413the volume fraction to lie between zero and one. 
     414 
     415To activate a user-defined distribution, put it in a file such as *distname.py* 
     416in the *SAS_WEIGHTS_PATH* folder.  This is defined with an environment 
     417variable, defaulting to:: 
     418 
     419    SAS_WEIGHTS_PATH=~/.sasview/weights 
     420 
     421The weights path is loaded on startup.  To update the distribution definition 
     422in a running application you will need to enter the following python commands:: 
     423 
     424    import sasmodels.weights 
     425    sasmodels.weights.load_weights('path/to/distname.py') 
     426 
     427.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
     428 
    336429Note about DLS polydispersity 
    337430^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 
    338431 
    339 Several measures of polydispersity abound in Dynamic Light Scattering (DLS) and  
    340 it should not be assumed that any of the following can be simply equated with  
     432Several measures of polydispersity abound in Dynamic Light Scattering (DLS) and 
     433it should not be assumed that any of the following can be simply equated with 
    341434the polydispersity *PD* parameter used in SasView. 
    342435 
    343 The dimensionless **Polydispersity Index (PI)** is a measure of the width of the  
    344 distribution of autocorrelation function decay rates (*not* the distribution of  
    345 particle sizes itself, though the two are inversely related) and is defined by  
     436The dimensionless **Polydispersity Index (PI)** is a measure of the width of the 
     437distribution of autocorrelation function decay rates (*not* the distribution of 
     438particle sizes itself, though the two are inversely related) and is defined by 
    346439ISO 22412:2017 as 
    347440 
    … …  
    350443    PI = \mu_{2} / \bar \Gamma^2 
    351444 
    352 where $\mu_\text{2}$ is the second cumulant, and $\bar \Gamma^2$ is the  
     445where $\mu_\text{2}$ is the second cumulant, and $\bar \Gamma^2$ is the 
    353446intensity-weighted average value, of the distribution of decay rates. 
    354447 
    … …  
    359452    PI = \sigma^2 / 2\bar \Gamma^2 
    360453 
    361 where $\sigma$ is the standard deviation, allowing a **Relative Polydispersity (RP)**  
     454where $\sigma$ is the standard deviation, allowing a **Relative Polydispersity (RP)** 
    362455to be defined as 
    363456 
    … …  
    366459    RP = \sigma / \bar \Gamma = \sqrt{2 \cdot PI} 
    367460 
    368 PI values smaller than 0.05 indicate a highly monodisperse system. Values  
     461PI values smaller than 0.05 indicate a highly monodisperse system. Values 
    369462greater than 0.7 indicate significant polydispersity. 
    370463 
    371 The **size polydispersity P-parameter** is defined as the relative standard  
    372 deviation coefficient of variation   
     464The **size polydispersity P-parameter** is defined as the relative standard 
     465deviation coefficient of variation 
    373466 
    374467.. math:: 
    … …  
    377470 
    378471where $\nu$ is the variance of the distribution and $\bar R$ is the mean 
    379 value of $R$. Here, the product $P \bar R$ is *equal* to the standard  
     472value of $R$. Here, the product $P \bar R$ is *equal* to the standard 
    380473deviation of the Lognormal distribution. 
    381474 
Note: See TracChangeset for help on using the changeset viewer.