Changeset 8696a87 in sasmodels


Ignore:
Timestamp:
Feb 1, 2016 4:38:33 AM (9 years ago)
Author:
piotr
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
f12357f, 577912b
Parents:
504abee (diff), 5c962df (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 remote-tracking branch 'origin/master'

Files:
5 added
19 edited

Legend:

Unmodified
Added
Removed

  • extra/pylint.rc

    rd4666ca r823e620  
    1010 
    1111# Profiled execution. 
    12 profile=no 
     12#profile=no 
    1313 
    1414# Add files or directories to the blacklist. They should be base names, not 
    … …  
    9292# Add a comment according to your evaluation note. This is used by the global 
    9393# evaluation report (RP0004). 
    94 comment=no 
     94#comment=no 
    9595 
    9696# Template used to display messages. This is a python new-style format string 
    … …  
    101101 
    102102# Required attributes for module, separated by a comma 
    103 required-attributes= 
     103#required-attributes= 
    104104 
    105105# List of builtins function names that should not be used, separated by a comma 
    … …  
    108108# Good variable names which should always be accepted, separated by a comma 
    109109good-names=_, 
    110     input, 
    111     h,i,j,k,n,w,x,y,z, 
    112     q,qx,qy,qz, 
    113     dt,dx,dy,dz,id, 
    114     Iq,dIq,Qx,Qy,Qz,Iqxy, 
    115     p, 
     110    input, id, 
     111    h, i, j, k, n, p, v, w, x, y, z, 
     112    q, qx, qy, qz, Q, Qx, Qy, Qz, 
     113    dt, dx, dy, dz, dq, 
     114    Iq, dIq, Iqxy, Iq_calc, 
    116115    ER, call_ER, VR, call_VR, 
    117116    Rmax, SElength, 
    … …  
    281280# (useful for modules/projects where namespaces are manipulated during runtime 
    282281# and thus existing member attributes cannot be deduced by static analysis 
    283 ignored-modules=numpy,np 
     282ignored-modules=numpy,np,numpy.random, 
     283    bumps,sas, 
    284284 
    285285# List of classes names for which member attributes should not be checked 
    … …  
    289289# When zope mode is activated, add a predefined set of Zope acquired attributes 
    290290# to generated-members. 
    291 zope=no 
     291#zope=no 
    292292 
    293293# List of members which are set dynamically and missed by pylint inference 
    … …  
    319319# List of interface methods to ignore, separated by a comma. This is used for 
    320320# instance to not check methods defines in Zope's Interface base class. 
    321 ignore-iface-methods=isImplementedBy,deferred,extends,names,namesAndDescriptions,queryDescriptionFor,getBases,getDescriptionFor,getDoc,getName,getTaggedValue,getTaggedValueTags,isEqualOrExtendedBy,setTaggedValue,isImplementedByInstancesOf,adaptWith,is_implemented_by 
     321#ignore-iface-methods=isImplementedBy,deferred,extends,names,namesAndDescriptions,queryDescriptionFor,getBases,getDescriptionFor,getDoc,getName,getTaggedValue,getTaggedValueTags,isEqualOrExtendedBy,setTaggedValue,isImplementedByInstancesOf,adaptWith,is_implemented_by 
    322322 
    323323# List of method names used to declare (i.e. assign) instance attributes. 

  • extra/pylint_numpy.py

    r63b32bb r823e620  
    88    #print("processing",module.name) 
    99    if module.name.startswith('numpy'): 
    10         if module.name == 'numpy':  import numpy 
    11         elif module.name == 'numpy.random': import numpy.random 
     10        if module.name == 'numpy': 
     11            import numpy 
     12        elif module.name == 'numpy.random': 
     13            import numpy.random 
     14            from numpy.random import seed, get_state, set_state 
    1215 

  • extra/pylint_pyopencl.py

    r3c56da87 r823e620  
    99    if module.name == 'pyopencl': 
    1010        import pyopencl 
     11        import pyopencl as cl 

  • extra/run-pylint.py

    rf01e53d r823e620  
    66 
    77def main(): 
     8    extra = abspath(dirname(__file__)) 
     9    root = abspath(joinpath(extra, '..')) 
     10 
    811    envpath = os.environ.get('PYTHONPATH',None) 
    912    path = [envpath] if envpath else [] 
    10     path.append(abspath(dirname(__file__))) # so we can find the plugins 
     13    path.append(extra) 
     14 
     15    #bumps = abspath(joinpath(root, "..", "bumps")) 
     16    #periodictable = abspath(joinpath(root, "..", "periodictable")) 
     17    #sasview = abspath(joinpath(root, "..", "sasview", "src")) 
     18    #path.extend((bumps, periodictable, sasview)) 
     19 
    1120    os.environ['PYTHONPATH'] = ':'.join(path) 
    12     root = abspath(joinpath(dirname(__file__), '..')) 
     21 
     22    # Run the lint command 
     23    cmd = "pylint --rcfile extra/pylint.rc -f parseable sasmodels" 
    1324    os.chdir(root) 
    14     cmd = "pylint --rcfile extra/pylint.rc -f parseable sasmodels" 
    1525    status = os.system(cmd) 
    1626    sys.exit(status) 

  • sasmodels/compare.py

    r190fc2b rd07c883  
    604604def columnize(L, indent="", width=79): 
    605605    """ 
    606     Format a list of strings into columns for printing. 
     606    Format a list of strings into columns. 
     607 
     608    Returns a string with carriage returns ready for printing. 
    607609    """ 
    608610    column_width = max(len(w) for w in L) + 1 

  • sasmodels/core.py

    r190fc2b reafc9fa  
    3838    """ 
    3939    Load a model definition given the model name. 
     40 
     41    This returns a handle to the module defining the model.  This can be 
     42    used with functions in generate to build the docs or extract model info. 
    4043    """ 
    4144    __import__('sasmodels.models.'+model_name) 
    … …  
    115118    Return a computation kernel from the model definition and the q input. 
    116119    """ 
    117     model_input = model.make_input(q_vectors) 
    118     return model(model_input) 
     120    return model(q_vectors) 
    119121 
    120122def get_weights(info, pars, name): 

  • sasmodels/data.py

    r69ec80f r5c962df  
    8787 
    8888class Data1D(object): 
     89    """ 
     90    1D data object. 
     91 
     92    Note that this definition matches the attributes from sasview, with 
     93    some generic 1D data vectors and some SAS specific definitions.  Some 
     94    refactoring to allow consistent naming conventions between 1D, 2D and 
     95    SESANS data would be helpful. 
     96 
     97    **Attributes** 
     98 
     99    *x*, *dx*: $q$ vector and gaussian resolution 
     100 
     101    *y*, *dy*: $I(q)$ vector and measurement uncertainty 
     102 
     103    *mask*: values to include in plotting/analysis 
     104 
     105    *dxl*: slit widths for slit smeared data, with *dx* ignored 
     106 
     107    *qmin*, *qmax*: range of $q$ values in *x* 
     108 
     109    *filename*: label for the data line 
     110 
     111    *_xaxis*, *_xunit*: label and units for the *x* axis 
     112 
     113    *_yaxis*, *_yunit*: label and units for the *y* axis 
     114    """ 
    89115    def __init__(self, x=None, y=None, dx=None, dy=None): 
    90116        self.x, self.y, self.dx, self.dy = x, y, dx, dy 
    … …  
    93119        self.qmin = x.min() if x is not None else np.NaN 
    94120        self.qmax = x.max() if x is not None else np.NaN 
    95         self.mask = np.isnan(y) if y is not None else None 
     121        # TODO: why is 1D mask False and 2D mask True? 
     122        self.mask = (np.isnan(y) if y is not None 
     123                     else np.zeros_like(x, 'b') if x is not None 
     124                     else None) 
    96125        self._xaxis, self._xunit = "x", "" 
    97126        self._yaxis, self._yunit = "y", "" 
    … …  
    114143 
    115144class Data2D(object): 
     145    """ 
     146    2D data object. 
     147 
     148    Note that this definition matches the attributes from sasview. Some 
     149    refactoring to allow consistent naming conventions between 1D, 2D and 
     150    SESANS data would be helpful. 
     151 
     152    **Attributes** 
     153 
     154    *qx_data*, *dqx_data*: $q_x$ matrix and gaussian resolution 
     155 
     156    *qy_data*, *dqy_data*: $q_y$ matrix and gaussian resolution 
     157 
     158    *data*, *err_data*: $I(q)$ matrix and measurement uncertainty 
     159 
     160    *mask*: values to exclude from plotting/analysis 
     161 
     162    *qmin*, *qmax*: range of $q$ values in *x* 
     163 
     164    *filename*: label for the data line 
     165 
     166    *_xaxis*, *_xunit*: label and units for the *x* axis 
     167 
     168    *_yaxis*, *_yunit*: label and units for the *y* axis 
     169 
     170    *_zaxis*, *_zunit*: label and units for the *y* axis 
     171 
     172    *Q_unit*, *I_unit*: units for Q and intensity 
     173 
     174    *x_bins*, *y_bins*: grid steps in *x* and *y* directions 
     175    """ 
    116176    def __init__(self, x=None, y=None, z=None, dx=None, dy=None, dz=None): 
    117177        self.qx_data, self.dqx_data = x, dx 
    118178        self.qy_data, self.dqy_data = y, dy 
    119179        self.data, self.err_data = z, dz 
    120         self.mask = ~np.isnan(z) if z is not None else None 
     180        self.mask = (~np.isnan(z) if z is not None 
     181                     else np.ones_like(x) if x is not None 
     182                     else None) 
    121183        self.q_data = np.sqrt(x**2 + y**2) 
    122184        self.qmin = 1e-16 
    … …  
    126188        self.Q_unit = "1/A" 
    127189        self.I_unit = "1/cm" 
    128         self.xaxis("Q_x", "A^{-1}") 
    129         self.yaxis("Q_y", "A^{-1}") 
    130         self.zaxis("Intensity", r"\text{cm}^{-1}") 
     190        self.xaxis("Q_x", "1/A") 
     191        self.yaxis("Q_y", "1/A") 
     192        self.zaxis("Intensity", "1/cm") 
    131193        self._xaxis, self._xunit = "x", "" 
    132194        self._yaxis, self._yunit = "y", "" 
    … …  
    157219 
    158220class Vector(object): 
     221    """ 
     222    3-space vector of *x*, *y*, *z* 
     223    """ 
    159224    def __init__(self, x=None, y=None, z=None): 
    160225        self.x, self.y, self.z = x, y, z 
    … …  
    235300    """ 
    236301    Plot data loaded by the sasview loader. 
     302 
     303    *data* is a sasview data object, either 1D, 2D or SESANS. 
     304 
     305    *view* is log or linear. 
     306 
     307    *limits* sets the intensity limits on the plot; if None then the limits 
     308    are inferred from the data. 
    237309    """ 
    238310    # Note: kind of weird using the plot result functions to plot just the 
    … …  
    249321def plot_theory(data, theory, resid=None, view='log', 
    250322                use_data=True, limits=None): 
     323    """ 
     324    Plot theory calculation. 
     325 
     326    *data* is needed to define the graph properties such as labels and 
     327    units, and to define the data mask. 
     328 
     329    *theory* is a matrix of the same shape as the data. 
     330 
     331    *view* is log or linear 
     332 
     333    *use_data* is True if the data should be plotted as well as the theory. 
     334 
     335    *limits* sets the intensity limits on the plot; if None then the limits 
     336    are inferred from the data. 
     337    """ 
    251338    if hasattr(data, 'lam'): 
    252339        _plot_result_sesans(data, theory, resid, use_data=True, limits=limits) 
    … …  
    258345 
    259346def protect(fn): 
     347    """ 
     348    Decorator to wrap calls in an exception trapper which prints the 
     349    exception and continues.  Keyboard interrupts are ignored. 
     350    """ 
    260351    def wrapper(*args, **kw): 
     352        """ 
     353        Trap and print errors from function. 
     354        """ 
    261355        try: 
    262356            return fn(*args, **kw) 
     357        except KeyboardInterrupt: 
     358            raise 
    263359        except: 
    264360            traceback.print_exc() 
    … …  
    332428@protect 
    333429def _plot_result_sesans(data, theory, resid, use_data, limits=None): 
     430    """ 
     431    Plot SESANS results. 
     432    """ 
    334433    import matplotlib.pyplot as plt 
    335434    use_data = use_data and data.y is not None 
    … …  
    459558 
    460559def demo(): 
     560    """ 
     561    Load and plot a SAS dataset. 
     562    """ 
    461563    data = load_data('DEC07086.DAT') 
    462564    set_beam_stop(data, 0.004) 

  • sasmodels/direct_model.py

    r9404dd3 reafc9fa  
    1 import warnings 
     1""" 
     2Class interface to the model calculator. 
     3 
     4Calling a model is somewhat non-trivial since the functions called depend 
     5on the data type.  For 1D data the *Iq* kernel needs to be called, for 
     62D data the *Iqxy* kernel needs to be called, and for SESANS data the 
     7*Iq* kernel needs to be called followed by a Hankel transform.  Before 
     8the kernel is called an appropriate *q* calculation vector needs to be 
     9constructed.  This is not the simple *q* vector where you have measured 
     10the data since the resolution calculation will require values beyond the 
     11range of the measured data.  After the calculation the resolution calculator 
     12must be called to return the predicted value for each measured data point. 
     13 
     14:class:`DirectModel` is a callable object that takes *parameter=value* 
     15keyword arguments and returns the appropriate theory values for the data. 
     16 
     17:class:`DataMixin` does the real work of interpreting the data and calling 
     18the model calculator.  This is used by :class:`DirectModel`, which uses 
     19direct parameter values and by :class:`bumps_model.Experiment` which wraps 
     20the parameter values in boxes so that the user can set fitting ranges, etc. 
     21on the individual parameters and send the model to the Bumps optimizers. 
     22""" 
     23from __future__ import print_function 
    224 
    325import numpy as np 
    … …  
    1941    currently used by :class:`sasview_model.SasviewModel` since this will 
    2042    require a number of changes to SasView before we can do it. 
     43 
     44    :meth:`_interpret_data` initializes the data structures necessary 
     45    to manage the calculations.  This sets attributes in the child class 
     46    such as *data_type* and *resolution*. 
     47 
     48    :meth:`_calc_theory` evaluates the model at the given control values. 
     49 
     50    :meth:`_set_data` sets the intensity data in the data object, 
     51    possibly with random noise added.  This is useful for simulating a 
     52    dataset with the results from :meth:`_calc_theory`. 
    2153    """ 
    2254    def _interpret_data(self, data, model): 
     55        # pylint: disable=attribute-defined-outside-init 
     56 
    2357        self._data = data 
    2458        self._model = model 
    … …  
    3670        if self.data_type == 'sesans': 
    3771            q = sesans.make_q(data.sample.zacceptance, data.Rmax) 
    38             self.index = slice(None, None) 
     72            index = slice(None, None) 
     73            res = None 
    3974            if data.y is not None: 
    40                 self.Iq = data.y 
    41                 self.dIq = data.dy 
     75                Iq, dIq = data.y, data.dy 
     76            else: 
     77                Iq, dIq = None, None 
    4278            #self._theory = np.zeros_like(q) 
    4379            q_vectors = [q] 
    … …  
    4985            qmax = getattr(data, 'qmax', np.inf) 
    5086            accuracy = getattr(data, 'accuracy', 'Low') 
    51             self.index = ~data.mask & (q >= qmin) & (q <= qmax) 
     87            index = ~data.mask & (q >= qmin) & (q <= qmax) 
    5288            if data.data is not None: 
    53                 self.index &= ~np.isnan(data.data) 
    54                 self.Iq = data.data[self.index] 
    55                 self.dIq = data.err_data[self.index] 
    56             self.resolution = resolution2d.Pinhole2D(data=data, index=self.index, 
    57                                                      nsigma=3.0, accuracy=accuracy) 
     89                index &= ~np.isnan(data.data) 
     90                Iq = data.data[index] 
     91                dIq = data.err_data[index] 
     92            else: 
     93                Iq, dIq = None, None 
     94            res = resolution2d.Pinhole2D(data=data, index=index, 
     95                                         nsigma=3.0, accuracy=accuracy) 
    5896            #self._theory = np.zeros_like(self.Iq) 
    59             q_vectors = self.resolution.q_calc 
     97            q_vectors = res.q_calc 
    6098        elif self.data_type == 'Iq': 
    61             self.index = (data.x >= data.qmin) & (data.x <= data.qmax) 
     99            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
    62100            if data.y is not None: 
    63                 self.index &= ~np.isnan(data.y) 
    64                 self.Iq = data.y[self.index] 
    65                 self.dIq = data.dy[self.index] 
     101                index &= ~np.isnan(data.y) 
     102                Iq = data.y[index] 
     103                dIq = data.dy[index] 
     104            else: 
     105                Iq, dIq = None, None 
    66106            if getattr(data, 'dx', None) is not None: 
    67                 q, dq = data.x[self.index], data.dx[self.index] 
    68                 if (dq>0).any(): 
    69                     self.resolution = resolution.Pinhole1D(q, dq) 
     107                q, dq = data.x[index], data.dx[index] 
     108                if (dq > 0).any(): 
     109                    res = resolution.Pinhole1D(q, dq) 
    70110                else: 
    71                     self.resolution = resolution.Perfect1D(q) 
    72             elif (getattr(data, 'dxl', None) is not None and 
    73                           getattr(data, 'dxw', None) is not None): 
    74                 self.resolution = resolution.Slit1D(data.x[self.index], 
    75                                                     width=data.dxh[self.index], 
    76                                                     height=data.dxw[self.index]) 
    77             else: 
    78                 self.resolution = resolution.Perfect1D(data.x[self.index]) 
     111                    res = resolution.Perfect1D(q) 
     112            elif (getattr(data, 'dxl', None) is not None 
     113                  and getattr(data, 'dxw', None) is not None): 
     114                res = resolution.Slit1D(data.x[index], 
     115                                        width=data.dxh[index], 
     116                                        height=data.dxw[index]) 
     117            else: 
     118                res = resolution.Perfect1D(data.x[index]) 
    79119 
    80120            #self._theory = np.zeros_like(self.Iq) 
    81             q_vectors = [self.resolution.q_calc] 
     121            q_vectors = [res.q_calc] 
    82122        else: 
    83123            raise ValueError("Unknown data type") # never gets here 
    … …  
    87127        self._kernel_inputs = [v for v in q_vectors] 
    88128        self._kernel = None 
     129        self.Iq, self.dIq, self.index = Iq, dIq, index 
     130        self.resolution = res 
    89131 
    90132    def _set_data(self, Iq, noise=None): 
     133        # pylint: disable=attribute-defined-outside-init 
    91134        if noise is not None: 
    92135            self.dIq = Iq*noise*0.01 
    … …  
    106149    def _calc_theory(self, pars, cutoff=0.0): 
    107150        if self._kernel is None: 
    108             q_input = self._model.make_input(self._kernel_inputs) 
    109             self._kernel = self._model(q_input) 
     151            self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-defined-outside-init 
    110152 
    111153        Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) 
    … …  
    120162 
    121163class DirectModel(DataMixin): 
     164    """ 
     165    Create a calculator object for a model. 
     166 
     167    *data* is 1D SAS, 2D SAS or SESANS data 
     168 
     169    *model* is a model calculator return from :func:`generate.load_model` 
     170 
     171    *cutoff* is the polydispersity weight cutoff. 
     172    """ 
    122173    def __init__(self, data, model, cutoff=1e-5): 
    123174        self.model = model 
    124175        self.cutoff = cutoff 
     176        # Note: _interpret_data defines the model attributes 
    125177        self._interpret_data(data, model) 
    126         self.kernel = make_kernel(self.model, self._kernel_inputs) 
     178 
    127179    def __call__(self, **pars): 
    128180        return self._calc_theory(pars, cutoff=self.cutoff) 
     181 
    129182    def ER(self, **pars): 
     183        """ 
     184        Compute the equivalent radius for the model. 
     185 
     186        Return 0. if not defined. 
     187        """ 
    130188        return call_ER(self.model.info, pars) 
     189 
    131190    def VR(self, **pars): 
     191        """ 
     192        Compute the equivalent volume for the model, including polydispersity 
     193        effects. 
     194 
     195        Return 1. if not defined. 
     196        """ 
    132197        return call_VR(self.model.info, pars) 
     198 
    133199    def simulate_data(self, noise=None, **pars): 
     200        """ 
     201        Generate simulated data for the model. 
     202        """ 
    134203        Iq = self.__call__(**pars) 
    135204        self._set_data(Iq, noise=noise) 
    136205 
    137 def demo(): 
     206def main(): 
     207    """ 
     208    Program to evaluate a particular model at a set of q values. 
     209    """ 
    138210    import sys 
    139211    from .data import empty_data1D, empty_data2D 
    … …  
    144216    model_name = sys.argv[1] 
    145217    call = sys.argv[2].upper() 
    146     if call not in ("ER","VR"): 
     218    if call not in ("ER", "VR"): 
    147219        try: 
    148220            values = [float(v) for v in call.split(',')] 
    … …  
    153225            data = empty_data1D([q]) 
    154226        elif len(values) == 2: 
    155             qx,qy = values 
    156             data = empty_data2D([qx],[qy]) 
     227            qx, qy = values 
     228            data = empty_data2D([qx], [qy]) 
    157229        else: 
    158230            print("use q or qx,qy or ER or VR") 
    … …  
    164236    model = load_model(model_definition, dtype='single') 
    165237    calculator = DirectModel(data, model) 
    166     pars = dict((k,float(v)) 
     238    pars = dict((k, float(v)) 
    167239                for pair in sys.argv[3:] 
    168                 for k,v in [pair.split('=')]) 
     240                for k, v in [pair.split('=')]) 
    169241    if call == "ER": 
    170242        print(calculator.ER(**pars)) 
    … …  
    176248 
    177249if __name__ == "__main__": 
    178     demo() 
     250    main() 

  • sasmodels/exception.py

    reaca9eb r823e620  
    88except NameError: 
    99    class WindowsError(Exception): 
     10        """ 
     11        Fake WindowsException when not on Windows. 
     12        """ 
    1013        pass 
    1114 
    … …  
    3134    # TODO: try to incorporate current stack trace in the raised exception 
    3235    if isinstance(exc, WindowsError): 
    33         raise OSError(str(exc)+" "+msg) 
     36        raise OSError(str(exc) + " " + msg) 
    3437 
    3538    args = exc.args 
    … …  
    3841    else: 
    3942        try: 
    40             arg0 = " ".join((args[0],msg)) 
     43            arg0 = " ".join((args[0], msg)) 
    4144            exc.args = tuple([arg0] + list(args[1:])) 
    4245        except: 
    43             exc.args = (" ".join((str(exc),msg)),) 
     46            exc.args = (" ".join((str(exc), msg)),) 

  • sasmodels/generate.py

    r190fc2b reafc9fa  
    77    particular dimensions averaged over all orientations. 
    88 
    9     *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx,qy for a form 
     9    *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy for a form 
    1010    with particular dimensions for a single orientation. 
    1111 
    … …  
    4747the *sin* and *cos* values are needed for a particular argument.  Since 
    4848this function does not exist in C99, all use of *sincos* should be 
    49 replaced by the macro *SINCOS(value,sn,cn)* where *sn* and *cn* are 
     49replaced by the macro *SINCOS(value, sn, cn)* where *sn* and *cn* are 
    5050previously declared *double* variables.  When compiled for systems without 
    5151OpenCL, *SINCOS* will be replaced by *sin* and *cos* calls.   If *value* is 
    … …  
    194194returns a list of files required by the model. 
    195195""" 
     196from __future__ import print_function 
    196197 
    197198# TODO: identify model files which have changed since loading and reload them. 
    … …  
    312313    raise ValueError("%r not found in %s" % (filename, search_path)) 
    313314 
    314 def sources(info): 
     315def model_sources(info): 
    315316    """ 
    316317    Return a list of the sources file paths for the module. 
    … …  
    372373 
    373374 
    374 def kernel_name(info, is_2D): 
     375def kernel_name(info, is_2d): 
    375376    """ 
    376377    Name of the exported kernel symbol. 
    377378    """ 
    378     return info['name'] + "_" + ("Iqxy" if is_2D else "Iq") 
     379    return info['name'] + "_" + ("Iqxy" if is_2d else "Iq") 
    379380 
    380381 
    … …  
    419420 
    420421 
    421 def build_polydispersity_loops(pd_pars): 
    422     """ 
    423     Build polydispersity loops 
    424  
    425     Returns loop opening and loop closing 
    426     """ 
    427     LOOP_OPEN = """\ 
     422LOOP_OPEN = """\ 
    428423for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) { 
    429424  const double %(name)s = loops[2*(%(name)s_i%(offset)s)]; 
    430425  const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\ 
    431426""" 
     427def build_polydispersity_loops(pd_pars): 
     428    """ 
     429    Build polydispersity loops 
     430 
     431    Returns loop opening and loop closing 
     432    """ 
    432433    depth = 4 
    433434    offset = "" 
    … …  
    464465 
    465466    # Load additional sources 
    466     source = [open(f).read() for f in sources(info)] 
     467    source = [open(f).read() for f in model_sources(info)] 
    467468 
    468469    # Prepare defines 
    … …  
    572573 
    573574    #for d in defines: print(d) 
    574     DEFINES = '\n'.join('#define %s %s' % (k, v) for k, v in defines) 
    575     SOURCES = '\n\n'.join(source) 
     575    defines = '\n'.join('#define %s %s' % (k, v) for k, v in defines) 
     576    sources = '\n\n'.join(source) 
    576577    return C_KERNEL_TEMPLATE % { 
    577         'DEFINES':DEFINES, 
    578         'SOURCES':SOURCES, 
     578        'DEFINES': defines, 
     579        'SOURCES': sources, 
    579580        } 
    580581 
    … …  
    590591    demo_parameters = getattr(kernel_module, 'demo', None) 
    591592    if demo_parameters is None: 
    592         demo_parameters = dict((p[0],p[2]) for p in parameters) 
     593        demo_parameters = dict((p[0], p[2]) for p in parameters) 
    593594    filename = abspath(kernel_module.__file__) 
    594595    kernel_id = splitext(basename(filename))[0] 
    … …  
    597598        name = " ".join(w.capitalize() for w in kernel_id.split('_')) 
    598599    info = dict( 
    599         id = kernel_id,  # string used to load the kernel 
     600        id=kernel_id,  # string used to load the kernel 
    600601        filename=abspath(kernel_module.__file__), 
    601602        name=name, 
    … …  
    643644        elif section_marker.match(line): 
    644645            if len(line) >= len(prior): 
    645                 yield "".join( ("**",prior,"**") ) 
     646                yield "".join(("**", prior, "**")) 
    646647                prior = None 
    647648            else: 

  • sasmodels/kernelcl.py

    rcde11f0f reafc9fa  
    11""" 
    2 GPU support through OpenCL 
     2GPU driver for C kernels 
    33 
    44There should be a single GPU environment running on the system.  This 
    … …  
    152152 
    153153 
    154 def make_result(self, size): 
    155     self.res = np.empty(size, dtype=self.dtype) 
    156     self.res_b = cl.Buffer(self.program.context, mf.READ_WRITE, self.res.nbytes) 
    157     return self.res, self.res_b 
    158  
    159  
    160154# for now, this returns one device in the context 
    161155# TODO: create a context that contains all devices on all platforms 
    … …  
    183177 
    184178    def has_type(self, dtype): 
     179        """ 
     180        Return True if all devices support a given type. 
     181        """ 
    185182        dtype = generate.F32 if dtype == 'fast' else np.dtype(dtype) 
    186183        return all(has_type(d, dtype) for d in self.context.devices) 
    187184 
    188185    def _create_some_context(self): 
     186        """ 
     187        Protected call to cl.create_some_context without interactivity.  Use 
     188        this if PYOPENCL_CTX is set in the environment.  Sets the *context* 
     189        attribute. 
     190        """ 
    189191        try: 
    190192            self.context = cl.create_some_context(interactive=False) 
    … …  
    195197 
    196198    def compile_program(self, name, source, dtype, fast=False): 
     199        """ 
     200        Compile the program for the device in the given context. 
     201        """ 
    197202        key = "%s-%s-%s"%(name, dtype, fast) 
    198203        if key not in self.compiled: 
    … …  
    204209 
    205210    def release_program(self, name): 
     211        """ 
     212        Free memory associated with the program on the device. 
     213        """ 
    206214        if name in self.compiled: 
    207215            self.compiled[name].release() 
    … …  
    209217 
    210218def _get_default_context(): 
     219    """ 
     220    Get an OpenCL context, preferring GPU over CPU. 
     221    """ 
    211222    default = None 
    212223    for platform in cl.get_platforms(): 
    … …  
    241252        self.info = info 
    242253        self.source = source 
    243         self.dtype = generate.F32 if dtype=='fast' else np.dtype(dtype) 
     254        self.dtype = generate.F32 if dtype == 'fast' else np.dtype(dtype) 
    244255        self.fast = (dtype == 'fast') 
    245256        self.program = None # delay program creation 
    246257 
    247258    def __getstate__(self): 
    248         state = self.__dict__.copy() 
    249         state['program'] = None 
    250         return state 
     259        return self.info, self.source, self.dtype, self.fast 
    251260 
    252261    def __setstate__(self, state): 
    253         self.__dict__ = state.copy() 
    254  
    255     def __call__(self, q_input): 
    256         if self.dtype != q_input.dtype: 
    257             raise TypeError("data is %s kernel is %s" 
    258                             % (q_input.dtype, self.dtype)) 
     262        self.info, self.source, self.dtype, self.fast = state 
     263        self.program = None 
     264 
     265    def __call__(self, q_vectors): 
    259266        if self.program is None: 
    260267            compiler = environment().compile_program 
    261268            self.program = compiler(self.info['name'], self.source, self.dtype, 
    262269                                    self.fast) 
    263         kernel_name = generate.kernel_name(self.info, q_input.is_2D) 
     270        is_2d = len(q_vectors) == 2 
     271        kernel_name = generate.kernel_name(self.info, is_2d) 
    264272        kernel = getattr(self.program, kernel_name) 
    265         return GpuKernel(kernel, self.info, q_input) 
     273        return GpuKernel(kernel, self.info, q_vectors, self.dtype) 
    266274 
    267275    def release(self): 
     276        """ 
     277        Free the resources associated with the model. 
     278        """ 
    268279        if self.program is not None: 
    269280            environment().release_program(self.info['name']) 
    270281            self.program = None 
    271282 
    272     def make_input(self, q_vectors): 
    273         """ 
    274         Make q input vectors available to the model. 
    275  
    276         Note that each model needs its own q vector even if the case of 
    277         mixture models because some models may be OpenCL, some may be 
    278         ctypes and some may be pure python. 
    279         """ 
    280         return GpuInput(q_vectors, dtype=self.dtype) 
     283    def __del__(self): 
     284        self.release() 
    281285 
    282286# TODO: check that we don't need a destructor for buffers which go out of scope 
    … …  
    304308        self.nq = q_vectors[0].size 
    305309        self.dtype = np.dtype(dtype) 
    306         self.is_2D = (len(q_vectors) == 2) 
     310        self.is_2d = (len(q_vectors) == 2) 
    307311        # TODO: stretch input based on get_warp() 
    308312        # not doing it now since warp depends on kernel, which is not known 
    … …  
    317321 
    318322    def release(self): 
     323        """ 
     324        Free the memory. 
     325        """ 
    319326        for b in self.q_buffers: 
    320327            b.release() 
    321328        self.q_buffers = [] 
    322329 
     330    def __del__(self): 
     331        self.release() 
     332 
    323333class GpuKernel(object): 
    324334    """ 
    325335    Callable SAS kernel. 
    326336 
    327     *kernel* is the GpuKernel object to call. 
     337    *kernel* is the GpuKernel object to call 
    328338 
    329339    *info* is the module information 
    330340 
    331     *q_input* is the DllInput q vectors at which the kernel should be 
    332     evaluated. 
     341    *q_vectors* is the q vectors at which the kernel should be evaluated 
     342 
     343    *dtype* is the kernel precision 
    333344 
    334345    The resulting call method takes the *pars*, a list of values for 
    … …  
    340351    Call :meth:`release` when done with the kernel instance. 
    341352    """ 
    342     def __init__(self, kernel, info, q_input): 
    343         self.q_input = q_input 
     353    def __init__(self, kernel, info, q_vectors, dtype): 
     354        q_input = GpuInput(q_vectors, dtype) 
    344355        self.kernel = kernel 
    345356        self.info = info 
    346357        self.res = np.empty(q_input.nq, q_input.dtype) 
    347         dim = '2d' if q_input.is_2D else '1d' 
     358        dim = '2d' if q_input.is_2d else '1d' 
    348359        self.fixed_pars = info['partype']['fixed-' + dim] 
    349360        self.pd_pars = info['partype']['pd-' + dim] 
    … …  
    358369                                q_input.global_size[0] * q_input.dtype.itemsize) 
    359370                      for _ in env.queues] 
    360  
     371        self.q_input = q_input 
    361372 
    362373    def __call__(self, fixed_pars, pd_pars, cutoff=1e-5): 
    … …  
    399410 
    400411    def release(self): 
     412        """ 
     413        Release resources associated with the kernel. 
     414        """ 
    401415        for b in self.loops_b: 
    402416            b.release() 
    … …  
    405419            b.release() 
    406420        self.res_b = [] 
     421        self.q_input.release() 
    407422 
    408423    def __del__(self): 

  • sasmodels/kerneldll.py

    r74667d3 reafc9fa  
    11r""" 
    2 C types wrapper for sasview models. 
     2DLL driver for C kernels 
    33 
    44The global attribute *ALLOW_SINGLE_PRECISION_DLLS* should be set to *True* if 
    … …  
    4444directory for this application, then OpenMP should be supported. 
    4545""" 
     46from __future__ import print_function 
    4647 
    4748import sys 
    … …  
    122123    models are allowed as DLLs. 
    123124    """ 
    124     if callable(info.get('Iq',None)): 
     125    if callable(info.get('Iq', None)): 
    125126        return PyModel(info) 
    126127 
    … …  
    139140 
    140141    source = generate.convert_type(source, dtype) 
    141     source_files = generate.sources(info) + [info['filename']] 
    142     dll= dll_path(info, dtype) 
     142    source_files = generate.model_sources(info) + [info['filename']] 
     143    dll = dll_path(info, dtype) 
    143144    newest = max(os.path.getmtime(f) for f in source_files) 
    144     if not os.path.exists(dll) or os.path.getmtime(dll)<newest: 
     145    if not os.path.exists(dll) or os.path.getmtime(dll) < newest: 
    145146        # Replace with a proper temp file 
    146         fid, filename = tempfile.mkstemp(suffix=".c",prefix=tempfile_prefix) 
    147         os.fdopen(fid,"w").write(source) 
     147        fid, filename = tempfile.mkstemp(suffix=".c", prefix=tempfile_prefix) 
     148        os.fdopen(fid, "w").write(source) 
    148149        command = COMPILE%{"source":filename, "output":dll} 
    149150        print("Compile command: "+command) 
    … …  
    160161def load_dll(source, info, dtype="double"): 
    161162    """ 
    162     Create and load a dll corresponding to the source,info pair returned 
     163    Create and load a dll corresponding to the source, info pair returned 
    163164    from :func:`sasmodels.generate.make` compiled for the target precision. 
    164165 
    … …  
    199200        Npd2d = len(self.info['partype']['pd-2d']) 
    200201 
    201         #print("dll",self.dllpath) 
     202        #print("dll", self.dllpath) 
    202203        try: 
    203204            self.dll = ct.CDLL(self.dllpath) 
    … …  
    210211              else c_longdouble) 
    211212        pd_args_1d = [c_void_p, fp] + [c_int]*Npd1d if Npd1d else [] 
    212         pd_args_2d= [c_void_p, fp] + [c_int]*Npd2d if Npd2d else [] 
     213        pd_args_2d = [c_void_p, fp] + [c_int]*Npd2d if Npd2d else [] 
    213214        self.Iq = self.dll[generate.kernel_name(self.info, False)] 
    214215        self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d 
    … …  
    218219 
    219220    def __getstate__(self): 
    220         return {'info': self.info, 'dllpath': self.dllpath, 'dll': None} 
     221        return self.info, self.dllpath 
    221222 
    222223    def __setstate__(self, state): 
    223         self.__dict__ = state 
    224  
    225     def __call__(self, q_input): 
    226         if self.dtype != q_input.dtype: 
    227             raise TypeError("data is %s kernel is %s" % (q_input.dtype, self.dtype)) 
     224        self.info, self.dllpath = state 
     225        self.dll = None 
     226 
     227    def __call__(self, q_vectors): 
     228        q_input = PyInput(q_vectors, self.dtype) 
    228229        if self.dll is None: self._load_dll() 
    229         kernel = self.Iqxy if q_input.is_2D else self.Iq 
     230        kernel = self.Iqxy if q_input.is_2d else self.Iq 
    230231        return DllKernel(kernel, self.info, q_input) 
    231232 
    232     # pylint: disable=no-self-use 
    233     def make_input(self, q_vectors): 
    234         """ 
    235         Make q input vectors available to the model. 
    236  
    237         Note that each model needs its own q vector even if the case of 
    238         mixture models because some models may be OpenCL, some may be 
    239         ctypes and some may be pure python. 
    240         """ 
    241         return PyInput(q_vectors, dtype=self.dtype) 
    242  
    243233    def release(self): 
     234        """ 
     235        Release any resources associated with the model. 
     236        """ 
    244237        pass # TODO: should release the dll 
    245238 
    … …  
    257250 
    258251    The resulting call method takes the *pars*, a list of values for 
    259     the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight) 
     252    the fixed parameters to the kernel, and *pd_pars*, a list of (value, weight) 
    260253    vectors for the polydisperse parameters.  *cutoff* determines the 
    261254    integration limits: any points with combined weight less than *cutoff* 
    … …  
    269262        self.kernel = kernel 
    270263        self.res = np.empty(q_input.nq, q_input.dtype) 
    271         dim = '2d' if q_input.is_2D else '1d' 
     264        dim = '2d' if q_input.is_2d else '1d' 
    272265        self.fixed_pars = info['partype']['fixed-'+dim] 
    273266        self.pd_pars = info['partype']['pd-'+dim] 
    … …  
    299292 
    300293    def release(self): 
     294        """ 
     295        Release any resources associated with the kernel. 
     296        """ 
    301297        pass 

  • sasmodels/kernelpy.py

    r4c2c535 reafc9fa  
     1""" 
     2Python driver for python kernels 
     3 
     4Calls the kernel with a vector of $q$ values for a single parameter set. 
     5Polydispersity is supported by looping over different parameter sets and 
     6summing the results.  The interface to :class:`PyModel` matches those for 
     7:class:`kernelcl.GpuModel` and :class:`kerneldll.DllModel`. 
     8""" 
    19import numpy as np 
    210from numpy import pi, cos 
    … …  
    513 
    614class PyModel(object): 
     15    """ 
     16    Wrapper for pure python models. 
     17    """ 
    718    def __init__(self, info): 
    819        self.info = info 
    920 
    10     def __call__(self, q_input): 
    11         kernel = self.info['Iqxy'] if q_input.is_2D else self.info['Iq'] 
     21    def __call__(self, q_vectors): 
     22        q_input = PyInput(q_vectors, dtype=F64) 
     23        kernel = self.info['Iqxy'] if q_input.is_2d else self.info['Iq'] 
    1224        return PyKernel(kernel, self.info, q_input) 
    1325 
    14     # pylint: disable=no-self-use 
    15     def make_input(self, q_vectors): 
    16         return PyInput(q_vectors, dtype=F64) 
    17  
    1826    def release(self): 
     27        """ 
     28        Free resources associated with the model. 
     29        """ 
    1930        pass 
    2031 
    … …  
    4152        self.nq = q_vectors[0].size 
    4253        self.dtype = dtype 
    43         self.is_2D = (len(q_vectors) == 2) 
     54        self.is_2d = (len(q_vectors) == 2) 
    4455        self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors] 
    4556        self.q_pointers = [q.ctypes.data for q in self.q_vectors] 
    4657 
    4758    def release(self): 
     59        """ 
     60        Free resources associated with the model inputs. 
     61        """ 
    4862        self.q_vectors = [] 
    4963 
    … …  
    7185        self.q_input = q_input 
    7286        self.res = np.empty(q_input.nq, q_input.dtype) 
    73         dim = '2d' if q_input.is_2D else '1d' 
     87        dim = '2d' if q_input.is_2d else '1d' 
    7488        # Loop over q unless user promises that the kernel is vectorized by 
    7589        # taggining it with vectorized=True 
    … …  
    7791            if dim == '2d': 
    7892                def vector_kernel(qx, qy, *args): 
     93                    """ 
     94                    Vectorized 2D kernel. 
     95                    """ 
    7996                    return np.array([kernel(qxi, qyi, *args) 
    8097                                     for qxi, qyi in zip(qx, qy)]) 
    8198            else: 
    8299                def vector_kernel(q, *args): 
     100                    """ 
     101                    Vectorized 1D kernel. 
     102                    """ 
    83103                    return np.array([kernel(qi, *args) 
    84104                                     for qi in q]) 
    … …  
    122142 
    123143    def release(self): 
     144        """ 
     145        Free resources associated with the kernel. 
     146        """ 
    124147        self.q_input = None 
    125148 

  • sasmodels/model_test.py

    r9404dd3 r5c962df  
    99    if model1 is 'all', then all except the remaining models will be tested 
    1010 
    11 Each model is tested using the default parameters at q=0.1, (qx,qy)=(0.1,0.1), 
     11Each model is tested using the default parameters at q=0.1, (qx, qy)=(0.1, 0.1), 
    1212and the ER and VR are computed.  The return values at these points are not 
    1313considered.  The test is only to verify that the models run to completion, 
    … …  
    3030        [ {parameters}, (qx, qy), I(qx, Iqy)], 
    3131        [ {parameters}, [(qx1, qy1), (qx2, qy2), ...], 
    32                         [I(qx1,qy1), I(qx2,qy2), ...]], 
     32                        [I(qx1, qy1), I(qx2, qy2), ...]], 
    3333 
    3434        [ {parameters}, 'ER', ER(pars) ], 
    … …  
    4343Precision defaults to 5 digits (relative). 
    4444""" 
     45from __future__ import print_function 
    4546 
    4647import sys 
    … …  
    5556 
    5657def make_suite(loaders, models): 
     58    """ 
     59    Construct the pyunit test suite. 
     60 
     61    *loaders* is the list of kernel drivers to use, which is one of 
     62    *["dll", "opencl"]*, *["dll"]* or *["opencl"]*.  For python models, 
     63    the python driver is always used. 
     64 
     65    *models* is the list of models to test, or *["all"]* to test all models. 
     66    """ 
    5767 
    5868    ModelTestCase = _hide_model_case_from_nosetests() 
    … …  
    7585        # don't try to call cl kernel since it will not be 
    7686        # available in some environmentes. 
    77         is_py = callable(getattr(model_definition,'Iq', None)) 
     87        is_py = callable(getattr(model_definition, 'Iq', None)) 
    7888 
    7989        if is_py:  # kernel implemented in python 
    … …  
    111121def _hide_model_case_from_nosetests(): 
    112122    class ModelTestCase(unittest.TestCase): 
     123        """ 
     124        Test suit for a particular model with a particular kernel driver. 
     125 
     126        The test suite runs a simple smoke test to make sure the model 
     127        functions, then runs the list of tests at the bottom of the model 
     128        description file. 
     129        """ 
    113130        def __init__(self, test_name, definition, test_method_name, 
    114131                     platform, dtype): 
    … …  
    123140        def _runTest(self): 
    124141            smoke_tests = [ 
    125                 [{},0.1,None], 
    126                 [{},(0.1,0.1),None], 
    127                 [{},'ER',None], 
    128                 [{},'VR',None], 
     142                [{}, 0.1, None], 
     143                [{}, (0.1, 0.1), None], 
     144                [{}, 'ER', None], 
     145                [{}, 'VR', None], 
    129146                ] 
    130147 
    … …  
    163180                actual = [call_VR(model.info, pars)] 
    164181            elif isinstance(x[0], tuple): 
    165                 Qx,Qy = zip(*x) 
     182                Qx, Qy = zip(*x) 
    166183                q_vectors = [np.array(Qx), np.array(Qy)] 
    167184                kernel = make_kernel(model, q_vectors) 
    … …  
    179196                    # smoke test --- make sure it runs and produces a value 
    180197                    self.assertTrue(np.isfinite(actual_yi), 
    181                         'invalid f(%s): %s' % (xi, actual_yi)) 
     198                                    'invalid f(%s): %s' % (xi, actual_yi)) 
    182199                else: 
    183200                    err = abs(yi - actual_yi) 
    184201                    nrm = abs(yi) 
    185202                    self.assertLess(err * 10**5, nrm, 
    186                         'f(%s); expected:%s; actual:%s' % (xi, yi, actual_yi)) 
     203                                    'f(%s); expected:%s; actual:%s' 
     204                                    % (xi, yi, actual_yi)) 
    187205 
    188206    return ModelTestCase 
    … …  
    237255    Run "nosetests sasmodels" on the command line to invoke it. 
    238256    """ 
    239     tests = make_suite(['opencl','dll'],['all']) 
     257    tests = make_suite(['opencl', 'dll'], ['all']) 
    240258    for test_i in tests: 
    241259        yield test_i._runTest 

  • sasmodels/resolution.py

    r190fc2b r5925e90  
    1414           "pinhole_extend_q", "slit_extend_q", "bin_edges", 
    1515           "interpolate", "linear_extrapolation", "geometric_extrapolation", 
    16            ] 
     16          ] 
    1717 
    1818MINIMUM_RESOLUTION = 1e-8 
    … …  
    8080        # default to Perfect1D if the pinhole geometry is not defined. 
    8181        self.q, self.q_width = q, q_width 
    82         self.q_calc = pinhole_extend_q(q, q_width, nsigma=nsigma) \ 
    83             if q_calc is None else np.sort(q_calc) 
    84         self.weight_matrix = pinhole_resolution(self.q_calc, 
    85                 self.q, np.maximum(q_width, MINIMUM_RESOLUTION)) 
     82        self.q_calc = (pinhole_extend_q(q, q_width, nsigma=nsigma) 
     83                       if q_calc is None else np.sort(q_calc)) 
     84        self.weight_matrix = pinhole_resolution( 
     85            self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION)) 
    8686 
    8787    def apply(self, theory): 
    … …  
    470470 
    471471def eval_form(q, form, pars): 
     472    """ 
     473    Return the SAS model evaluated at *q*. 
     474 
     475    *form* is the SAS model returned from :fun:`core.load_model`. 
     476 
     477    *pars* are the parameter values to use when evaluating. 
     478    """ 
    472479    from sasmodels import core 
    473480    kernel = core.make_kernel(form, [q]) 
    … …  
    478485 
    479486def gaussian(q, q0, dq): 
     487    """ 
     488    Return the Gaussian resolution function. 
     489 
     490    *q0* is the center, *dq* is the width and *q* are the points to evaluate. 
     491    """ 
    480492    from numpy import exp, pi 
    481493    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq) 
    … …  
    559571 
    560572class ResolutionTest(unittest.TestCase): 
     573    """ 
     574    Test the resolution calculations. 
     575    """ 
    561576 
    562577    def setUp(self): 
    … …  
    668683 
    669684class IgorComparisonTest(unittest.TestCase): 
     685    """ 
     686    Test resolution calculations against those returned by Igor. 
     687    """ 
    670688 
    671689    def setUp(self): 
    … …  
    675693        self.model = core.load_model(sphere, dtype='double') 
    676694 
    677     def Iq_sphere(self, pars, resolution): 
     695    def _eval_sphere(self, pars, resolution): 
    678696        from sasmodels import core 
    679697        kernel = core.make_kernel(self.model, [resolution.q_calc]) 
    … …  
    683701        return result 
    684702 
    685     def compare(self, q, output, answer, tolerance): 
     703    def _compare(self, q, output, answer, tolerance): 
    686704        #err = (output - answer)/answer 
    687705        #idx = abs(err) >= tolerance 
    … …  
    700718        q, width, answer, _ = data 
    701719        resolution = Perfect1D(q) 
    702         output = self.Iq_sphere(pars, resolution) 
    703         self.compare(q, output, answer, 1e-6) 
     720        output = self._eval_sphere(pars, resolution) 
     721        self._compare(q, output, answer, 1e-6) 
    704722 
    705723    def test_pinhole(self): 
    … …  
    713731        q, q_width, answer = data 
    714732        resolution = Pinhole1D(q, q_width) 
    715         output = self.Iq_sphere(pars, resolution) 
     733        output = self._eval_sphere(pars, resolution) 
    716734        # TODO: relative error should be lower 
    717         self.compare(q, output, answer, 3e-4) 
     735        self._compare(q, output, answer, 3e-4) 
    718736 
    719737    def test_pinhole_romberg(self): 
    … …  
    736754        q_calc, tol = None, 0.01 
    737755        resolution = Pinhole1D(q, q_width, q_calc=q_calc) 
    738         output = self.Iq_sphere(pars, resolution) 
     756        output = self._eval_sphere(pars, resolution) 
    739757        if 0: # debug plot 
    740758            import matplotlib.pyplot as plt 
    741759            resolution = Perfect1D(q) 
    742             source = self.Iq_sphere(pars, resolution) 
     760            source = self._eval_sphere(pars, resolution) 
    743761            plt.loglog(q, source, '.') 
    744762            plt.loglog(q, answer, '-', hold=True) 
    745763            plt.loglog(q, output, '-', hold=True) 
    746764            plt.show() 
    747         self.compare(q, output, answer, tol) 
     765        self._compare(q, output, answer, tol) 
    748766 
    749767    def test_slit(self): 
    … …  
    757775        q, delta_qv, _, answer = data 
    758776        resolution = Slit1D(q, width=delta_qv, height=0) 
    759         output = self.Iq_sphere(pars, resolution) 
     777        output = self._eval_sphere(pars, resolution) 
    760778        # TODO: eliminate Igor test since it is too inaccurate to be useful. 
    761779        # This means we can eliminate the test data as well, and instead 
    762780        # use a generated q vector. 
    763         self.compare(q, output, answer, 0.5) 
     781        self._compare(q, output, answer, 0.5) 
    764782 
    765783    def test_slit_romberg(self): 
    … …  
    777795                               delta_qv[0], 0.) 
    778796        resolution = Slit1D(q, width=delta_qv, height=0, q_calc=q_calc) 
    779         output = self.Iq_sphere(pars, resolution) 
     797        output = self._eval_sphere(pars, resolution) 
    780798        # TODO: relative error should be lower 
    781         self.compare(q, output, answer, 0.025) 
     799        self._compare(q, output, answer, 0.025) 
    782800 
    783801    def test_ellipsoid(self): 
    … …  
    799817        # TODO: 10% is too much error; use better algorithm 
    800818        #print(np.max(abs(answer-output)/answer)) 
    801         self.compare(q, output, answer, 0.1) 
     819        self._compare(q, output, answer, 0.1) 
    802820 
    803821    #TODO: can sas q spacing be too sparse for the resolution calculation? 
    … …  
    813831        q, q_width, answer = data[:, ::20] # Take every nth point 
    814832        resolution = Pinhole1D(q, q_width) 
    815         output = self.Iq_sphere(pars, resolution) 
    816         self.compare(q, output, answer, 1e-6) 
     833        output = self._eval_sphere(pars, resolution) 
     834        self._compare(q, output, answer, 1e-6) 
    817835 
    818836 
    … …  
    10671085 
    10681086def demo_pinhole_1d(): 
     1087    """ 
     1088    Show example of pinhole smearing. 
     1089    """ 
    10691090    q = np.logspace(-4, np.log10(0.2), 400) 
    10701091    q_width = 0.1*q 
    … …  
    10731094 
    10741095def demo_slit_1d(): 
     1096    """ 
     1097    Show example of slit smearing. 
     1098    """ 
    10751099    q = np.logspace(-4, np.log10(0.2), 100) 
    10761100    w = h = 0. 
    … …  
    10831107 
    10841108def demo(): 
     1109    """ 
     1110    Run the resolution demos. 
     1111    """ 
    10851112    import matplotlib.pyplot as plt 
    10861113    plt.subplot(121) 

  • sasmodels/resolution2d.py

    r841753c r823e620  
    143143        if self.coords == 'polar': 
    144144            q_r = sqrt(qx**2 + qy**2) 
    145             qx_res = ((dqx*cos(dphi) + q_r) * cos(-q_phi) + 
    146                           dqy*sin(dphi) * sin(-q_phi)) 
    147             qy_res = (-(dqx*cos(dphi) + q_r) * sin(-q_phi) + 
    148                           dqy*sin(dphi) * cos(-q_phi)) 
     145            qx_res = ((dqx*cos(dphi) + q_r) * cos(-q_phi) 
     146                      + dqy*sin(dphi) * sin(-q_phi)) 
     147            qy_res = (-(dqx*cos(dphi) + q_r) * sin(-q_phi) 
     148                      + dqy*sin(dphi) * cos(-q_phi)) 
    149149        else: 
    150             qx_res = qx +  dqx*cos(dphi) 
    151             qy_res = qy +  dqy*sin(dphi) 
     150            qx_res = qx + dqx*cos(dphi) 
     151            qy_res = qy + dqy*sin(dphi) 
    152152 
    153153 
    … …  
    162162            theory = np.reshape(theory, (nbins, nq)) 
    163163            ## Averaging with Gaussian weighting: normalization included. 
    164             value =np.average(theory, axis=0, weights=self.q_calc_weights) 
     164            value = np.average(theory, axis=0, weights=self.q_calc_weights) 
    165165            ## Return the smeared values in the range of self.index 
    166166            return value 

  • sasmodels/sasview_model.py

    r9404dd3 reafc9fa  
    283283        """ 
    284284        q_vectors = [np.asarray(q) for q in args] 
    285         fn = self._model(self._model.make_input(q_vectors)) 
     285        fn = self._model(q_vectors) 
    286286        pars = [self.params[v] for v in fn.fixed_pars] 
    287287        pd_pars = [self._get_weights(p) for p in fn.pd_pars] 

  • sasmodels/weights.py

    r3c56da87 r5c962df  
    2222 
    2323    def get_pars(self): 
     24        """ 
     25        Return the parameters to the disperser as a dictionary. 
     26        """ 
    2427        pars = {'type': self.type} 
    2528        pars.update(self.__dict__) 
    … …  
    2831    # pylint: disable=no-self-use 
    2932    def set_weights(self, values, weights): 
     33        """ 
     34        Set the weights on the disperser if it is :class:`ArrayDispersion`. 
     35        """ 
    3036        raise RuntimeError("set_weights is only available for ArrayDispersion") 
    3137 
    … …  
    3642        *center* is the center of the distribution 
    3743 
    38         *lb*,*ub* are the min and max allowed values 
     44        *lb*, *ub* are the min and max allowed values 
    3945 
    4046        *relative* is True if the distribution width is proportional to the 
    … …  
    6369 
    6470class GaussianDispersion(Dispersion): 
     71    r""" 
     72    Gaussian dispersion, with 1-\ $\sigma$ width. 
     73 
     74    .. math:: 
     75 
     76        w = \exp\left(-\tfrac12 (x - c)^2/\sigma^2\right) 
     77    """ 
    6578    type = "gaussian" 
    6679    default = dict(npts=35, width=0, nsigmas=3) 
    … …  
    7285 
    7386class RectangleDispersion(Dispersion): 
     87    r""" 
     88    Uniform dispersion, with width $\sqrt{3}\sigma$. 
     89 
     90    .. math:: 
     91 
     92        w = 1 
     93    """ 
    7494    type = "rectangle" 
    7595    default = dict(npts=35, width=0, nsigmas=1.70325) 
    … …  
    81101 
    82102class LogNormalDispersion(Dispersion): 
     103    r""" 
     104    log Gaussian dispersion, with 1-\ $\sigma$ width. 
     105 
     106    .. math:: 
     107 
     108        w = \frac{\exp\left(-\tfrac12 (\ln x - c)^2/\sigma^2\right)}{x\sigma} 
     109    """ 
    83110    type = "lognormal" 
    84111    default = dict(npts=80, width=0, nsigmas=8) 
    85112    def _weights(self, center, sigma, lb, ub): 
    86         x = self._linspace(center, sigma, max(lb,1e-8), max(ub,1e-8)) 
    87         px = np.exp(-0.5*(np.log(x)-center)**2)/sigma**2/(x*sigma) 
     113        x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8)) 
     114        px = np.exp(-0.5*(np.log(x)-center)**2/sigma**2)/(x*sigma) 
    88115        return x, px 
    89116 
    90117 
    91118class SchulzDispersion(Dispersion): 
     119    r""" 
     120    Schultz dispersion, with 1-\ $\sigma$ width. 
     121 
     122    .. math:: 
     123 
     124        w = \frac{z^z\,R^{z-1}}{e^{Rz}\,c \Gamma(z)} 
     125 
     126    where $c$ is the center of the distribution, $R = x/c$ and $z=(c/\sigma)^2$. 
     127 
     128    This is evaluated using logarithms as 
     129 
     130    .. math:: 
     131 
     132        w = \exp\left(z \ln z + (z-1)\ln R - Rz - \ln c - \ln \Gamma(z) \right) 
     133    """ 
    92134    type = "schulz" 
    93135    default = dict(npts=80, width=0, nsigmas=8) 
    94136    def _weights(self, center, sigma, lb, ub): 
    95         x = self._linspace(center, sigma, max(lb,1e-8), max(ub,1e-8)) 
    96         R= x/center 
     137        x = self._linspace(center, sigma, max(lb, 1e-8), max(ub, 1e-8)) 
     138        R = x/center 
    97139        z = (center/sigma)**2 
    98140        arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z) 
    … …  
    102144 
    103145class ArrayDispersion(Dispersion): 
     146    r""" 
     147    Empirical dispersion curve. 
     148 
     149    Use :meth:`set_weights` to set $w = f(x)$. 
     150    """ 
    104151    type = "array" 
    105152    default = dict(npts=35, width=0, nsigmas=1) 
    … …  
    110157 
    111158    def set_weights(self, values, weights): 
     159        """ 
     160        Set the weights for the given x values. 
     161        """ 
    112162        self.values = np.ascontiguousarray(values, 'd') 
    113163        self.weights = np.ascontiguousarray(weights, 'd') 
    … …  
    117167        # TODO: interpolate into the array dispersion using npts 
    118168        x = center + self.values*sigma 
    119         idx = (x>=lb)&(x<=ub) 
     169        idx = (x >= lb) & (x <= ub) 
    120170        x = x[idx] 
    121171        px = self.weights[idx] 
    … …  
    124174 
    125175# dispersion name -> disperser lookup table. 
    126 models = dict((d.type,d) for d in ( 
     176models = dict((d.type, d) for d in ( 
    127177    GaussianDispersion, RectangleDispersion, 
    128178    ArrayDispersion, SchulzDispersion, LogNormalDispersion 
    … …  
    149199    of the parameter, and false if it is an absolute width. 
    150200 
    151     Returns *(value,weight)*, where *value* and *weight* are vectors. 
     201    Returns *(value, weight)*, where *value* and *weight* are vectors. 
    152202    """ 
    153203    cls = models[disperser] 
    154204    obj = cls(n, width, nsigmas) 
    155     v,w = obj.get_weights(value, limits[0], limits[1], relative) 
    156     return v,w 
     205    v, w = obj.get_weights(value, limits[0], limits[1], relative) 
     206    return v, w 
    157207 
    158208# Hack to allow sasview dispersion objects to interoperate with sasmodels 
    159 dispersers = dict((v.__name__,k) for k,v in models.items()) 
     209dispersers = dict((v.__name__, k) for k, v in models.items()) 
    160210dispersers['DispersionModel'] = RectangleDispersion.type 
    161211 

  • sasmodels/models/lib/wrc_cyl.c

    rf94d8a2 r504abee  
    297297} 
    298298 
    299  
     299static double 
     300Sdebye_kernel(double arg) 
     301{ 
     302    // ORIGINAL 
     303    double result = 2.0*(exp(-arg) + arg -1.0)/(pow((arg),2)); 
     304 
     305    // CONVERSION 1 from http://herbie.uwplse.org/ 
     306    // 
     307    // exhibits discontinuity - needs more investigation 
     308    //double a1 = 1.0/6.0; 
     309    //double a2 = 1.0/72.0; 
     310    //double a3 = 1.0/24.0; 
     311    //double result = pow((1.0 - a1*arg - (a2+a3)*arg*arg), 2); 
     312 
     313    return result; 
     314} 
    300315static double 
    301316Sdebye(double q, double L, double b) 
    302317{ 
    303     return 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2)); 
     318    double arg = u_WR(q,L,b); 
     319    return Sdebye_kernel(arg); 
    304320} 
    305321 
    … …  
    308324Sdebye1(double q, double L, double b) 
    309325{ 
    310     return 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) ); 
     326    double arg = u1(q,L,b); 
     327    return Sdebye_kernel(arg); 
     328 
    311329} 
    312330 
Note: See TracChangeset for help on using the changeset viewer.