Changes in / [bf6631a:f94d8a2] in sasmodels


Ignore:
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • extra/pylint.rc

    r37a7252 r5ef0633  
    66# Python code to execute, usually for sys.path manipulation such as 
    77# pygtk.require(). 
    8 #init-hook='import os, sys; sys.path.extend([os.getcwd(),os.path.join(os.getcwd(),"extra")])' 
     8#init-hook= 
    99 
    1010# Profiled execution. 

  • sasmodels/__init__.py

    r37a7252 r9890053  
    1 """ 
    2 sasmodels 
    3 ========= 
    4  
    5 **sasmodels** is a package containing models for small angle neutron and xray 
    6 scattering.  Models supported are the one dimensional circular average and 
    7 two dimensional oriented patterns.  As well as the form factor calculations 
    8 for the individual shapes **sasmodels** also provides automatic shape 
    9 polydispersity, angular dispersion and resolution convolution.  SESANS 
    10 patterns can be computed for any model. 
    11  
    12 Models can be written in python or in C.  C models can run on the GPU if 
    13 OpenCL drivers are available.  See :mod:`generate` for details on 
    14 defining new models. 
    15 """ 
    16  
    171__version__ = "0.9" 

  • sasmodels/bumps_model.py

    r37a7252 rec7e360  
    66arguments to set the initial value for each parameter. 
    77 
    8 :class:`Experiment` combines the *Model* function with a data file loaded by 
    9 the sasview data loader.  *Experiment* takes a *cutoff* parameter controlling 
     8:class:`Experiment` combines the *Model* function with a data file loaded by the 
     9sasview data loader.  *Experiment* takes a *cutoff* parameter controlling 
    1010how far the polydispersity integral extends. 
    1111 
    1212""" 
    13  
    14 __all__ = [ 
    15     "Model", "Experiment", 
    16     ] 
    1713 
    1814import warnings 
    … …  
    2420 
    2521# CRUFT: old style bumps wrapper which doesn't separate data and model 
    26 # pylint: disable=invalid-name 
    2722def BumpsModel(data, model, cutoff=1e-5, **kw): 
    28     r""" 
    29     Bind a model to data, along with a polydispersity cutoff. 
    30  
    31     *data* is a :class:`data.Data1D`, :class:`data.Data2D` or 
    32     :class:`data.Sesans` object.  Use :func:`data.empty_data1D` or 
    33     :func:`data.empty_data2D` to define $q, \Delta q$ calculation 
    34     points for displaying the SANS curve when there is no measured data. 
    35  
    36     *model* is a runnable module as returned from :func:`core.load_model`. 
    37  
    38     *cutoff* is the polydispersity weight cutoff. 
    39  
    40     Any additional *key=value* pairs are model dependent parameters. 
    41  
    42     Returns an :class:`Experiment` object. 
    43  
    44     Note that the usual Bumps semantics is not fully supported, since 
    45     assigning *M.name = parameter* on the returned experiment object 
    46     does not set that parameter in the model.  Range setting will still 
    47     work as expected though. 
    48  
    49     .. deprecated:: 0.1 
    50         Use :class:`Experiment` instead. 
    51     """ 
    5223    warnings.warn("Use of BumpsModel is deprecated.  Use bumps_model.Experiment instead.") 
    53  
    54     # Create the model and experiment 
    5524    model = Model(model, **kw) 
    5625    experiment = Experiment(data=data, model=model, cutoff=cutoff) 
    57  
    58     # Copy the model parameters up to the experiment object. 
    59     for k, v in model.parameters().items(): 
    60         setattr(experiment, k, v) 
     26    for k in model._parameter_names: 
     27        setattr(experiment, k, getattr(model, k)) 
    6128    return experiment 
    6229 
    63  
    6430def create_parameters(model_info, **kwargs): 
    65     """ 
    66     Generate Bumps parameters from the model info. 
    67  
    68     *model_info* is returned from :func:`generate.model_info` on the 
    69     model definition module. 
    70  
    71     Any additional *key=value* pairs are initial values for the parameters 
    72     to the models.  Uninitialized parameters will use the model default 
    73     value. 
    74  
    75     Returns a dictionary of *{name: Parameter}* containing the bumps 
    76     parameters for each model parameter, and a dictionary of 
    77     *{name: str}* containing the polydispersity distribution types. 
    78     """ 
    7931    # lazy import; this allows the doc builder and nosetests to run even 
    8032    # when bumps is not on the path. 
    8133    from bumps.names import Parameter 
     34 
     35    partype = model_info['partype'] 
    8236 
    8337    pars = {} 
    … …  
    8640        value = kwargs.pop(name, default) 
    8741        pars[name] = Parameter.default(value, name=name, limits=limits) 
    88     for name in model_info['partype']['pd-2d']: 
     42    for name in partype['pd-2d']: 
    8943        for xpart, xdefault, xlimits in [ 
    90                 ('_pd', 0., limits), 
    91                 ('_pd_n', 35., (0, 1000)), 
    92                 ('_pd_nsigma', 3., (0, 10)), 
    93             ]: 
     44            ('_pd', 0., limits), 
     45            ('_pd_n', 35., (0, 1000)), 
     46            ('_pd_nsigma', 3., (0, 10)), 
     47        ]: 
    9448            xname = name + xpart 
    9549            xvalue = kwargs.pop(xname, xdefault) 
    … …  
    9751 
    9852    pd_types = {} 
    99     for name in model_info['partype']['pd-2d']: 
     53    for name in partype['pd-2d']: 
    10054        xname = name + '_pd_type' 
    10155        xvalue = kwargs.pop(xname, 'gaussian') 
    … …  
    10963 
    11064class Model(object): 
    111     """ 
    112     Bumps wrapper for a SAS model. 
    113  
    114     *model* is a runnable module as returned from :func:`core.load_model`. 
    115  
    116     *cutoff* is the polydispersity weight cutoff. 
    117  
    118     Any additional *key=value* pairs are model dependent parameters. 
    119     """ 
    12065    def __init__(self, model, **kwargs): 
    12166        self._sasmodel = model 
    12267        pars, pd_types = create_parameters(model.info, **kwargs) 
    123         for k, v in pars.items(): 
     68        for k,v in pars.items(): 
    12469            setattr(self, k, v) 
    125         for k, v in pd_types.items(): 
     70        for k,v in pd_types.items(): 
    12671            setattr(self, k, v) 
    12772        self._parameter_names = list(pars.keys()) 
    … …  
    13075    def parameters(self): 
    13176        """ 
    132         Return a dictionary of parameters objects for the parameters, 
    133         excluding polydispersity distribution type. 
     77        Return a dictionary of parameters 
    13478        """ 
    13579        return dict((k, getattr(self, k)) for k in self._parameter_names) 
    13680 
    13781    def state(self): 
    138         """ 
    139         Return a dictionary of current values for all the parameters, 
    140         including polydispersity distribution type. 
    141         """ 
    14282        pars = dict((k, getattr(self, k).value) for k in self._parameter_names) 
    14383        pars.update((k, getattr(self, k)) for k in self._pd_type_names) 
    … …  
    14585 
    14686class Experiment(DataMixin): 
    147     r""" 
    148     Bumps wrapper for a SAS experiment. 
     87    """ 
     88    Return a bumps wrapper for a SAS model. 
    14989 
    150     *data* is a :class:`data.Data1D`, :class:`data.Data2D` or 
    151     :class:`data.Sesans` object.  Use :func:`data.empty_data1D` or 
    152     :func:`data.empty_data2D` to define $q, \Delta q$ calculation 
    153     points for displaying the SANS curve when there is no measured data. 
     90    *data* is the data to be fitted. 
    15491 
    155     *model* is a :class:`Model` object. 
     92    *model* is the SAS model from :func:`core.load_model`. 
    15693 
    15794    *cutoff* is the integration cutoff, which avoids computing the 
    15895    the SAS model where the polydispersity weight is low. 
    15996 
    160     The resulting model can be used directly in a Bumps FitProblem call. 
     97    Model parameters can be initialized with additional keyword 
     98    arguments, or by assigning to model.parameter_name.value. 
     99 
     100    The resulting bumps model can be used directly in a FitProblem call. 
    161101    """ 
    162102    def __init__(self, data, model, cutoff=1e-5): 
    … …  
    169109 
    170110    def update(self): 
    171         """ 
    172         Call when model parameters have changed and theory needs to be 
    173         recalculated. 
    174         """ 
    175111        self._cache = {} 
    176112 
    177113    def numpoints(self): 
    178114        """ 
    179         Return the number of data points 
     115            Return the number of points 
    180116        """ 
    181117        return len(self.Iq) 
    … …  
    188124 
    189125    def theory(self): 
    190         """ 
    191         Return the theory corresponding to the model parameters. 
    192  
    193         This method uses lazy evaluation, and requires model.update() to be 
    194         called when the parameters have changed. 
    195         """ 
    196126        if 'theory' not in self._cache: 
    197127            pars = self.model.state() 
    … …  
    200130 
    201131    def residuals(self): 
    202         """ 
    203         Return theory minus data normalized by uncertainty. 
    204         """ 
    205132        #if np.any(self.err ==0): print("zeros in err") 
    206133        return (self.theory() - self.Iq) / self.dIq 
    207134 
    208135    def nllf(self): 
    209         """ 
    210         Return the negative log likelihood of seeing data given the model 
    211         parameters, up to a normalizing constant which depends on the data 
    212         uncertainty. 
    213         """ 
    214136        delta = self.residuals() 
    215137        #if np.any(np.isnan(R)): print("NaN in residuals") 
    … …  
    227149 
    228150    def simulate_data(self, noise=None): 
    229         """ 
    230         Generate simulated data. 
    231         """ 
    232151        Iq = self.theory() 
    233152        self._set_data(Iq, noise) 
    234153 
    235154    def save(self, basename): 
    236         """ 
    237         Save the model parameters and data into a file. 
    238  
    239         Not Implemented. 
    240         """ 
    241155        pass 
    242156 

  • sasmodels/compare.py

    rd15a908 r608e31e  
    7272# doc string so that we can display it at run time if there is an error. 
    7373# lin 
    74 __doc__ = (__doc__  # pylint: disable=redefined-builtin 
    75            + """ 
     74__doc__ = __doc__ + """ 
    7675Program description 
    7776------------------- 
    7877 
    79 """ 
    80            + USAGE) 
     78""" + USAGE 
    8179 
    8280 
    … …  
    283281            theory = lambda: smearer.get_value() 
    284282        else: 
    285             theory = lambda: model.evalDistribution([data.qx_data[index], 
    286                                                      data.qy_data[index]]) 
     283            theory = lambda: model.evalDistribution([data.qx_data[index], data.qy_data[index]]) 
    287284    elif smearer is not None: 
    288285        theory = lambda: smearer(model.evalDistribution(data.x)) 
    … …  
    419416        try: 
    420417            base_value, base_time = time_calculation(base, pars, Nbase) 
    421             print("%s t=%.1f ms, intensity=%.0f" 
    422                   % (base.engine, base_time, sum(base_value))) 
     418            print("%s t=%.1f ms, intensity=%.0f"%(base.engine, base_time, sum(base_value))) 
    423419        except ImportError: 
    424420            traceback.print_exc() 
    … …  
    430426        try: 
    431427            comp_value, comp_time = time_calculation(comp, pars, Ncomp) 
    432             print("%s t=%.1f ms, intensity=%.0f" 
    433                   % (comp.engine, comp_time, sum(comp_value))) 
     428            print("%s t=%.1f ms, intensity=%.0f"%(comp.engine, comp_time, sum(comp_value))) 
    434429        except ImportError: 
    435430            traceback.print_exc() 
    … …  
    438433    # Compare, but only if computing both forms 
    439434    if Nbase > 0 and Ncomp > 0: 
     435        #print("speedup %.2g"%(comp_time/base_time)) 
     436        #print("max |base/comp|", max(abs(base_value/comp_value)), "%.15g"%max(abs(base_value)), "%.15g"%max(abs(comp_value))) 
     437        #comp *= max(base_value/comp_value) 
    440438        resid = (base_value - comp_value) 
    441439        relerr = resid/comp_value 
    442         _print_stats("|%s-%s|" 
    443                      % (base.engine, comp.engine) + (" "*(3+len(comp.engine))), 
     440        _print_stats("|%s-%s|"%(base.engine, comp.engine) + (" "*(3+len(comp.engine))), 
    444441                     resid) 
    445         _print_stats("|(%s-%s)/%s|" 
    446                      % (base.engine, comp.engine, comp.engine), 
     442        _print_stats("|(%s-%s)/%s|"%(base.engine, comp.engine, comp.engine), 
    447443                     relerr) 
    448444 
    … …  
    595591    invalid = [o[1:] for o in flags 
    596592               if o[1:] not in NAME_OPTIONS 
    597                and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)] 
     593                   and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)] 
    598594    if invalid: 
    599595        print("Invalid options: %s"%(", ".join(invalid))) 
    … …  
    601597 
    602598 
    603     # pylint: disable=bad-whitespace 
    604599    # Interpret the flags 
    605600    opts = { 
    … …  
    656651        elif arg == '-sasview': engines.append(arg[1:]) 
    657652        elif arg == '-edit':    opts['explore'] = True 
    658     # pylint: enable=bad-whitespace 
    659653 
    660654    if len(engines) == 0: 
    … …  
    681675    presets = {} 
    682676    for arg in values: 
    683         k, v = arg.split('=', 1) 
     677        k,v = arg.split('=',1) 
    684678        if k not in pars: 
    685679            # extract base name without polydispersity info 
    686680            s = set(p.split('_pd')[0] for p in pars) 
    687             print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s)))) 
     681            print("%r invalid; parameters are: %s"%(k,", ".join(sorted(s)))) 
    688682            sys.exit(1) 
    689683        presets[k] = float(v) if not k.endswith('type') else v 
    … …  
    703697 
    704698    # Create the computational engines 
    705     data, _ = make_data(opts) 
     699    data, _index = make_data(opts) 
    706700    if n1: 
    707701        base = make_engine(model_definition, data, engines[0], opts['cutoff']) 
    … …  
    713707        comp = None 
    714708 
    715     # pylint: disable=bad-whitespace 
    716709    # Remember it all 
    717710    opts.update({ 
    … …  
    725718        'engines'   : [base, comp], 
    726719    }) 
    727     # pylint: enable=bad-whitespace 
    728720 
    729721    return opts 
    730722 
     723def main(): 
     724    opts = parse_opts() 
     725    if opts['explore']: 
     726        explore(opts) 
     727    else: 
     728        compare(opts) 
     729 
    731730def explore(opts): 
    732     """ 
    733     Explore the model using the Bumps GUI. 
    734     """ 
    735731    import wx 
    736732    from bumps.names import FitProblem 
    … …  
    738734 
    739735    problem = FitProblem(Explore(opts)) 
    740     is_mac = "cocoa" in wx.version() 
     736    isMac = "cocoa" in wx.version() 
    741737    app = wx.App() 
    742738    frame = AppFrame(parent=None, title="explore") 
    743     if not is_mac: frame.Show() 
     739    if not isMac: frame.Show() 
    744740    frame.panel.set_model(model=problem) 
    745741    frame.panel.Layout() 
    746742    frame.panel.aui.Split(0, wx.TOP) 
    747     if is_mac: frame.Show() 
     743    if isMac: frame.Show() 
    748744    app.MainLoop() 
    749745 
    750746class Explore(object): 
    751747    """ 
    752     Bumps wrapper for a SAS model comparison. 
    753  
    754     The resulting object can be used as a Bumps fit problem so that 
    755     parameters can be adjusted in the GUI, with plots updated on the fly. 
     748    Return a bumps wrapper for a SAS model comparison. 
    756749    """ 
    757750    def __init__(self, opts): 
    … …  
    794787        Return cost. 
    795788        """ 
    796         # pylint: disable=no-self-use 
    797789        return 0.  # No nllf 
    798790 
    … …  
    812804 
    813805 
    814 def main(): 
    815     """ 
    816     Main program. 
    817     """ 
    818     opts = parse_opts() 
    819     if opts['explore']: 
    820         explore(opts) 
    821     else: 
    822         compare(opts) 
    823  
    824806if __name__ == "__main__": 
    825807    main() 

  • sasmodels/compare_many.py

    rd15a908 r6458608  
    11#!/usr/bin/env python 
    2 """ 
    3 Program to compare results from many random parameter sets for a given model. 
    4  
    5 The result is a comma separated value (CSV) table that can be redirected 
    6 from standard output into a file and loaded into a spreadsheet. 
    7  
    8 The models are compared for each parameter set and if the difference is 
    9 greater than expected for that precision, the parameter set is labeled 
    10 as bad and written to the output, along with the random seed used to 
    11 generate that parameter value.  This seed can be used with :mod:`compare` 
    12 to reload and display the details of the model. 
    13 """ 
    142from __future__ import print_function 
    153 
    … …  
    2614 
    2715def calc_stats(target, value, index): 
    28     """ 
    29     Calculate statistics between the target value and the computed value. 
    30  
    31     *target* and *value* are the vectors being compared, with the 
    32     difference normalized by *target* to get relative error.  Only 
    33     the elements listed in *index* are used, though index may be 
    34     and empty slice defined by *slice(None, None)*. 
    35  
    36     Returns: 
    37  
    38         *maxrel* the maximum relative difference 
    39  
    40         *rel95* the relative difference with the 5% biggest differences ignored 
    41  
    42         *maxabs* the maximum absolute difference for the 5% biggest differences 
    43  
    44         *maxval* the maximum value for the 5% biggest differences 
    45     """ 
    4616    resid = abs(value-target)[index] 
    4717    relerr = resid/target[index] 
    48     sorted_rel_index = np.argsort(relerr) 
     18    srel = np.argsort(relerr) 
    4919    #p90 = int(len(relerr)*0.90) 
    5020    p95 = int(len(relerr)*0.95) 
    5121    maxrel = np.max(relerr) 
    52     rel95 = relerr[sorted_rel_index[p95]] 
    53     maxabs = np.max(resid[sorted_rel_index[p95:]]) 
    54     maxval = np.max(value[sorted_rel_index[p95:]]) 
    55     return maxrel, rel95, maxabs, maxval 
     22    rel95 = relerr[srel[p95]] 
     23    maxabs = np.max(resid[srel[p95:]]) 
     24    maxval = np.max(value[srel[p95:]]) 
     25    return maxrel,rel95,maxabs,maxval 
    5626 
    5727def print_column_headers(pars, parts): 
    58     """ 
    59     Generate column headers for the differences and for the parameters, 
    60     and print them to standard output. 
    61     """ 
    6228    stats = list('Max rel err|95% rel err|Max abs err above 90% rel|Max value above 90% rel'.split('|')) 
    6329    groups = [''] 
    … …  
    7036    print(','.join('"%s"'%c for c in columns)) 
    7137 
    72 # Target 'good' value for various precision levels. 
    7338PRECISION = { 
    7439    'fast': 1e-3, 
    … …  
    8348def compare_instance(name, data, index, N=1, mono=True, cutoff=1e-5, 
    8449                     base='sasview', comp='double'): 
    85     r""" 
    86     Compare the model under different calculation engines. 
    87  
    88     *name* is the name of the model. 
    89  
    90     *data* is the data object giving $q, \Delta q$ calculation points. 
    91  
    92     *index* is the active set of points. 
    93  
    94     *N* is the number of comparisons to make. 
    95  
    96     *cutoff* is the polydispersity weight cutoff to make the calculation 
    97     a little bit faster. 
    98  
    99     *base* and *comp* are the names of the calculation engines to compare. 
    100     """ 
    101  
    102     is_2d = hasattr(data, 'qx_data') 
     50    is2D = hasattr(data, 'qx_data') 
    10351    model_definition = core.load_model_definition(name) 
    10452    pars = get_demo_pars(model_definition) 
    10553    header = ('\n"Model","%s","Count","%d","Dimension","%s"' 
    106               % (name, N, "2D" if is_2d else "1D")) 
     54              % (name, N, "2D" if is2D else "1D")) 
    10755    if not mono: header += ',"Cutoff",%g'%(cutoff,) 
    10856    print(header) 
    10957 
    110     if is_2d: 
     58    if is2D: 
    11159        info = generate.make_info(model_definition) 
    11260        partype = info['partype'] 
    … …  
    12169    # declarations are not available in python 2.7. 
    12270    def try_model(fn, pars): 
    123         """ 
    124         Return the model evaluated at *pars*.  If there is an exception, 
    125         print it and return NaN of the right shape. 
    126         """ 
    12771        try: 
    12872            result = fn(**pars) 
    … …  
    13882        return result 
    13983    def check_model(pars): 
    140         """ 
    141         Run the two calculators against *pars*, returning statistics 
    142         on the differences.  See :func:`calc_stats` for the list of stats. 
    143         """ 
    14484        base_value = try_model(calc_base, pars) 
    14585        comp_value = try_model(calc_comp, pars) 
    … …  
    168108        good = [True] 
    169109        columns = check_model(pars_i) 
    170         columns += [v for _, v in sorted(pars_i.items())] 
     110        columns += [v for _,v in sorted(pars_i.items())] 
    171111        if first: 
    172112            labels = [" vs. ".join((calc_base.engine, calc_comp.engine))] 
    … …  
    181121 
    182122def print_usage(): 
    183     """ 
    184     Print the command usage string. 
    185     """ 
    186123    print("usage: compare_many.py MODEL COUNT (1dNQ|2dNQ) (CUTOFF|mono) (single|double|quad)") 
    187124 
    188125 
    189126def print_models(): 
    190     """ 
    191     Print the list of available models in columns. 
    192     """ 
    193127    print(columnize(MODELS, indent="  ")) 
    194128 
    195129 
    196130def print_help(): 
    197     """ 
    198     Print usage string, the option description and the list of available models. 
    199     """ 
    200131    print_usage() 
    201132    print("""\ 
    … …  
    227158 
    228159def main(): 
    229     """ 
    230     Main program. 
    231     """ 
    232     if len(sys.argv) not in (6, 7): 
     160    if len(sys.argv) not in (6,7): 
    233161        print_help() 
    234162        sys.exit(1) 
    … …  
    254182 
    255183    data, index = make_data({'qmax':1.0, 'is2d':is2D, 'nq':Nq, 'res':0., 
    256                              'accuracy': 'Low', 'view':'log'}) 
     184                              'accuracy': 'Low', 'view':'log'}) 
    257185    model_list = [model] if model != "all" else MODELS 
    258186    for model in model_list: 

  • sasmodels/convert.py

    rd15a908 r29da213  
    4949    """ 
    5050    return dict((p, (v*1e6 if p.endswith('sld') else v)) 
    51                 for p, v in pars.items()) 
     51                for p,v in pars.items()) 
    5252 
    5353def convert_model(name, pars): 
    … …  
    5555    Convert model from old style parameter names to new style. 
    5656    """ 
    57     _, _ = name, pars # lint 
     57    _,_ = name,pars # lint 
    5858    raise NotImplementedError 
    5959    # need to load all new models in order to determine old=>new 
    … …  
    6767    """ 
    6868    return dict((p, (v*1e-6 if p.endswith('sld') else v)) 
    69                 for p, v in pars.items()) 
     69                for p,v in pars.items()) 
    7070 
    7171def _remove_pd(pars, key, name): 

  • sasmodels/core.py

    rd15a908 rcde11f0f  
    22Core model handling routines. 
    33""" 
    4 __all__ = [ 
    5     "list_models", "load_model_definition", "precompile_dll", 
    6     "load_model", "make_kernel", "call_kernel", "call_ER", "call_VR", 
    7     ] 
     4__all__ = ["list_models", "load_model_definition",  "precompile_dll", 
     5           "load_model", "make_kernel", "call_kernel", "call_ER", "call_VR" ] 
    86 
    97from os.path import basename, dirname, join as joinpath 
    … …  
    6361 
    6462def isstr(s): 
    65     """ 
    66     Return True if *s* is a string-like object. 
    67     """ 
    6863    try: s + '' 
    6964    except: return False 
    … …  
    10499    # source = open(info['name']+'.cl','r').read() 
    105100 
    106     if (platform == "dll" 
     101    if (platform=="dll" 
    107102            or not HAVE_OPENCL 
    108103            or not kernelcl.environment().has_type(dtype)): 
    … …  
    133128    width = pars.get(name+'_pd', 0.0) 
    134129    nsigma = pars.get(name+'_pd_nsigma', 3.0) 
    135     value, weight = weights.get_weights( 
     130    value,weight = weights.get_weights( 
    136131        disperser, npts, width, nsigma, value, limits, relative) 
    137132    return value, weight / np.sum(weight) 
    … …  
    200195                    for name in info['partype']['volume']] 
    201196        value, weight = dispersion_mesh(vol_pars) 
    202         whole, part = VR(*value) 
     197        whole,part = VR(*value) 
    203198        return np.sum(weight*part)/np.sum(weight*whole) 
    204199 

  • sasmodels/data.py

    rd15a908 r013adb7  
    284284                mdata[mdata <= 0] = masked 
    285285            plt.errorbar(data.x/10, scale*mdata, yerr=data.dy, fmt='.') 
    286             all_positive = all_positive and (mdata > 0).all() 
     286            all_positive = all_positive and (mdata>0).all() 
    287287            some_present = some_present or (mdata.count() > 0) 
    288288 
    … …  
    292292            mtheory[~np.isfinite(mtheory)] = masked 
    293293            if view is 'log': 
    294                 mtheory[mtheory <= 0] = masked 
     294                mtheory[mtheory<=0] = masked 
    295295            plt.plot(data.x/10, scale*mtheory, '-', hold=True) 
    296             all_positive = all_positive and (mtheory > 0).all() 
     296            all_positive = all_positive and (mtheory>0).all() 
    297297            some_present = some_present or (mtheory.count() > 0) 
    298298 
    … …  
    396396        plt.title('theory') 
    397397        h = plt.colorbar() 
    398         h.set_label(r'$\log_{10}I(q)$' if view == 'log' 
     398        h.set_label(r'$\log_{10}I(q)$' if view=='log' 
    399399                    else r'$q^4 I(q)$' if view == 'q4' 
    400400                    else '$I(q)$') 
    … …  
    411411        plt.title('residuals') 
    412412        h = plt.colorbar() 
    413         h.set_label(r'$\Delta I(q)$') 
     413        h.set_label('$\Delta I(q)$') 
    414414 
    415415 

  • sasmodels/resolution.py

    rfdc538a r5258859  
    135135    """ 
    136136    #print("apply shapes", theory.shape, weight_matrix.shape) 
    137     Iq = np.dot(theory[None, :], weight_matrix) 
     137    Iq = np.dot(theory[None,:], weight_matrix) 
    138138    #print("result shape",Iq.shape) 
    139139    return Iq.flatten() 
    … …  
    153153    # neither trapezoid nor Simpson's rule improved the accuracy. 
    154154    edges = bin_edges(q_calc) 
    155     edges[edges < 0.0] = 0.0 # clip edges below zero 
    156     G = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :]) 
     155    edges[edges<0.0] = 0.0 # clip edges below zero 
     156    G = erf( (edges[:,None] - q[None,:]) / (sqrt(2.0)*q_width)[None,:] ) 
    157157    weights = G[1:] - G[:-1] 
    158     weights /= np.sum(weights, axis=0)[None, :] 
     158    weights /= np.sum(weights, axis=0)[None,:] 
    159159    return weights 
    160160 
    … …  
    287287    # The current algorithm is a midpoint rectangle rule. 
    288288    q_edges = bin_edges(q_calc) # Note: requires q > 0 
    289     q_edges[q_edges < 0.0] = 0.0 # clip edges below zero 
     289    q_edges[q_edges<0.0] = 0.0 # clip edges below zero 
    290290    weights = np.zeros((len(q), len(q_calc)), 'd') 
    291291 
    … …  
    306306            #print(qi - h, qi + h) 
    307307            #print(in_x + abs_x) 
    308             weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*h) 
     308            weights[i,:] = (in_x + abs_x) * np.diff(q_edges) / (2*h) 
    309309        else: 
    310310            L = n_height 
    311311            for k in range(-L, L+1): 
    312                 weights[i, :] += _q_perp_weights(q_edges, qi+k*h/L, w) 
    313             weights[i, :] /= 2*L + 1 
     312                weights[i,:] += _q_perp_weights(q_edges, qi+k*h/L, w) 
     313            weights[i,:] /= 2*L + 1 
    314314 
    315315    return weights.T 
    … …  
    358358    log-scaled data. 
    359359    """ 
    360     if len(x) < 2 or (np.diff(x) < 0).any(): 
     360    if len(x) < 2 or (np.diff(x)<0).any(): 
    361361        raise ValueError("Expected bins to be an increasing set") 
    362362    edges = np.hstack([ 
    … …  
    373373    """ 
    374374    step = np.diff(q) 
    375     index = step > max_step 
     375    index = step>max_step 
    376376    if np.any(index): 
    377377        inserts = [] 
    378         for q_i, step_i in zip(q[:-1][index], step[index]): 
     378        for q_i,step_i in zip(q[:-1][index],step[index]): 
    379379            n = np.ceil(step_i/max_step) 
    380             inserts.extend(q_i + np.arange(1, n)*(step_i/n)) 
     380            inserts.extend(q_i + np.arange(1,n)*(step_i/n)) 
    381381        # Extend a couple of fringes beyond the end of the data 
    382         inserts.extend(q[-1] + np.arange(1, 8)*max_step) 
    383         q_calc = np.sort(np.hstack((q, inserts))) 
     382        inserts.extend(q[-1] + np.arange(1,8)*max_step) 
     383        q_calc = np.sort(np.hstack((q,inserts))) 
    384384    else: 
    385385        q_calc = q 
    … …  
    399399    if q_min < q[0]: 
    400400        if q_min <= 0: q_min = q_min*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 
    401         n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15 
     401        n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1]>q[0] else 15 
    402402        q_low = np.linspace(q_min, q[0], n_low+1)[:-1] 
    403403    else: 
    404404        q_low = [] 
    405405    if q_max > q[-1]: 
    406         n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1] > q[-2] else 15 
     406        n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1]>q[-2] else 15 
    407407        q_high = np.linspace(q[-1], q_max, n_high+1)[1:] 
    408408    else: 
    … …  
    452452        if q_min < 0: q_min = q[0]*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 
    453453        n_low = log_delta_q * (log(q[0])-log(q_min)) 
    454         q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1] 
     454        q_low  = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1] 
    455455    else: 
    456456        q_low = [] 
    … …  
    489489    that make it slow to evaluate but give it good accuracy. 
    490490    """ 
    491     from scipy.integrate import romberg 
     491    from scipy.integrate import romberg, dblquad 
    492492 
    493493    if any(k not in form.info['defaults'] for k in pars.keys()): 
    … …  
    520520            result[i] = r/(2*h) 
    521521        else: 
    522             w_grid = np.linspace(0, w, 21)[None, :] 
    523             h_grid = np.linspace(-h, h, 23)[:, None] 
     522            w_grid = np.linspace(0, w, 21)[None,:] 
     523            h_grid = np.linspace(-h, h, 23)[:,None] 
    524524            u = sqrt((qi+h_grid)**2 + w_grid**2) 
    525525            Iu = np.interp(u, q_calc, Iq) 
    526526            #print(np.trapz(Iu, w_grid, axis=1)) 
    527             Is = np.trapz(np.trapz(Iu, w_grid, axis=1), h_grid[:, 0]) 
     527            Is = np.trapz(np.trapz(Iu, w_grid, axis=1), h_grid[:,0]) 
    528528            result[i] = Is / (2*h*w) 
    529             # from scipy.integrate import dblquad 
    530             # r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w, 
    531             #                  args=(qi,)) 
    532             # result[i] = r/(w*2*h) 
     529            """ 
     530            r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w, 
     531                             args=(qi,)) 
     532            result[i] = r/(w*2*h) 
     533            """ 
    533534 
    534535    # r should be [float, ...], but it is [array([float]), array([float]),...] 
    … …  
    552553 
    553554    _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq) 
    554     r = [romberg(_fn, max(qi-nsigma*dqi, 1e-10*q[0]), qi+nsigma*dqi, 
    555                  args=(qi, dqi), divmax=100, vec_func=True, tol=0, rtol=1e-8) 
    556          for qi, dqi in zip(q, q_width)] 
     555    r = [romberg(_fn, max(qi-nsigma*dqi,1e-10*q[0]), qi+nsigma*dqi, args=(qi, dqi), 
     556                 divmax=100, vec_func=True, tol=0, rtol=1e-8) 
     557         for qi,dqi in zip(q,q_width)] 
    557558    return np.asarray(r).flatten() 
    558559 
    … …  
    594595        theory = self.Iq(resolution.q_calc) 
    595596        output = resolution.apply(theory) 
    596         answer = [ 
    597             9.0618, 8.6402, 8.1187, 7.1392, 6.1528, 
    598             5.5555, 4.5584, 3.5606, 2.5623, 2.0000, 
    599             ] 
     597        answer = [ 9.0618, 8.6402, 8.1187, 7.1392, 6.1528, 
     598                   5.5555, 4.5584, 3.5606, 2.5623, 2.0000 ] 
    600599        np.testing.assert_allclose(output, answer, atol=1e-4) 
    601600 
    … …  
    609608        theory = 1000*self.Iq(resolution.q_calc**4) 
    610609        output = resolution.apply(theory) 
    611         answer = [ 
    612             8.85785, 8.43012, 7.92687, 6.94566, 6.03660, 
    613             5.40363, 4.40655, 3.40880, 2.41058, 2.00000, 
    614             ] 
     610        answer = [ 8.85785, 8.43012, 7.92687, 6.94566, 6.03660, 
     611                   5.40363, 4.40655, 3.40880, 2.41058, 2.00000 ] 
    615612        np.testing.assert_allclose(output, answer, atol=1e-4) 
    616613 
    … …  
    623620        theory = self.Iq(resolution.q_calc) 
    624621        output = resolution.apply(theory) 
    625         answer = [ 
    626             11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 
    627             ] 
     622        answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ] 
    628623        np.testing.assert_allclose(output, answer, atol=1e-4) 
    629624 
    … …  
    637632        theory = self.Iq(resolution.q_calc) 
    638633        output = resolution.apply(theory) 
    639         answer = [ 
    640             11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 
    641             ] 
     634        answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ] 
    642635        np.testing.assert_allclose(output, answer, atol=1e-4) 
    643636 
    … …  
    659652        theory = 12.0-1000.0*resolution.q_calc 
    660653        output = resolution.apply(theory) 
    661         answer = [ 
    662             10.44785079, 9.84991299, 8.98101708, 
    663             7.99906585, 6.99998311, 6.00001689, 
    664             5.00093415, 4.01898292, 3.15008701, 2.55214921, 
    665             ] 
     654        answer = [ 10.44785079, 9.84991299, 8.98101708, 
     655                  7.99906585, 6.99998311, 6.00001689, 
     656                  5.00093415, 4.01898292, 3.15008701, 2.55214921] 
    666657        np.testing.assert_allclose(output, answer, atol=1e-8) 
    667658 
    … …  
    792783            } 
    793784        form = load_model('ellipsoid', dtype='double') 
    794         q = np.logspace(log10(4e-5), log10(2.5e-2), 68) 
     785        q = np.logspace(log10(4e-5),log10(2.5e-2), 68) 
    795786        width, height = 0.117, 0. 
    796787        resolution = Slit1D(q, width=width, height=height) 
    … …  
    10801071    #h = 0.0277790 
    10811072    resolution = Slit1D(q, w, h) 
    1082     _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w, h)) 
     1073    _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w,h)) 
    10831074 
    10841075def demo(): 

  • sasmodels/sesans.py

    r384d114 r3c56da87  
    1212import numpy as np 
    1313from numpy import pi, exp 
     14 
    1415from scipy.special import jv as besselj 
    1516 
    16 def make_q(q_max, Rmax): 
    17     r""" 
    18     Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$ 
    19     to $q_max$. 
    20     """ 
     17def make_q(q_zmax, Rmax): 
    2118    q_min = dq = 0.1 * 2*pi / Rmax 
    22     return np.arange(q_min, q_max, dq) 
     19    #q_min = 0.00003 
     20    return np.arange(q_min, q_zmax, dq) 
     21 
     22# TODO: dead code; for now the call to the hankel transform happens in BumpsModel 
     23class SesansCalculator: 
     24    def __init__(self, kernel, q_zmax, Rmax, SElength, wavelength, thickness): 
     25        self._set_kernel(kernel, q_zmax, Rmax) 
     26        self.SElength = SElength 
     27        self.wavelength = wavelength 
     28        self.thickness = thickness 
     29 
     30    def _set_kernel(self, kernel, q_zmax, Rmax): 
     31        kernel_input = kernel.make_input([make_q(q_zmax, Rmax)]) 
     32        self.sans_calculator = kernel(kernel_input) 
     33 
     34    def __call__(self, pars, pd_pars, cutoff=1e-5): 
     35        Iq = self.sans_calculator(pars, pd_pars, cutoff) 
     36        P = hankel(self.SElength, self.wavelength, self.thickness, self.q, Iq) 
     37        self.Iq = Iq 
     38        return P 
    2339 
    2440def hankel(SElength, wavelength, thickness, q, Iq): 
    25     r""" 
     41    """ 
    2642    Compute the expected SESANS polarization for a given SANS pattern. 
    2743 
    28     Uses the hankel transform followed by the exponential.  The values for *zz* 
    29     (or spin echo length, or delta), wavelength and sample thickness should 
    30     come from the dataset.  $q$ should be chosen such that the oscillations 
    31     in $I(q)$ are well sampled (e.g., $5 \cdot 2 \pi/d_{\max}$). 
     44    Uses the hankel transform followed by the exponential.  The values 
     45    for zz (or spin echo length, or delta), wavelength and sample thickness 
     46    information should come from the dataset.  *q* should be chosen such 
     47    that the oscillations in *I(q)* are well sampled (e.g., 5*2*pi/d_max). 
    3248 
    33     *SElength* [A] is the set of $z$ points at which to compute the 
    34     Hankel transform 
     49    *SElength* [A] is the set of z points at which to compute the hankel transform 
    3550 
    3651    *wavelength* [m]  is the wavelength of each individual point *zz* 
    … …  
    3853    *thickness* [cm] is the sample thickness. 
    3954 
    40     *q* [A$^{-1}$] is the set of $q$ points at which the model has been 
    41     computed. These should be equally spaced. 
     55    *q* [A^{-1}] is the set of q points at which the model has been computed. 
     56    These should be equally spaced. 
    4257 
    43     *I* [cm$^{-1}$] is the value of the SANS model at *q* 
     58    *I* [cm^{-1}] is the value of the SANS model at *q* 
    4459    """ 
    4560    G = np.zeros(len(SElength), 'd') 
    4661    for i in range(len(SElength)): 
    47         integr = besselj(0, q*SElength[i])*Iq*q 
     62        integr = besselj(0,q*SElength[i])*Iq*q 
    4863        G[i] = np.sum(integr) 
    49  
    50     # [m^-1] step size in q, needed for integration 
    51     dq=(q[1]-q[0])*1e10 
    52  
    53     # integration step, convert q into [m**-1] and 2 pi circle integration 
    54     G *= dq*1e10*2*pi 
    55  
     64    dq=(q[1]-q[0])*1e10   # [m^-1] step size in q, needed for integration 
     65    G *= dq*1e10*2*pi # integr step, conver q into [m**-1] and 2 pi circle integr 
    5666    P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0])) 
    5767 
Note: See TracChangeset for help on using the changeset viewer.