Changes in / [81ec7c8:793beb3] in sasmodels


Ignore:
Files:
7 added
46 edited

Legend:

Unmodified
Added
Removed

  • doc/developer/index.rst

    r56fc97a rb85be2d  
    33*************************** 
    44 
     5.. toctree:: 
     6   :numbered: 4 
     7   :maxdepth: 4 
    58 
     9   calculator.rst 

  • sasmodels/alignment.py

    r3c56da87 r7ae2b7f  
    1818to decide whether it is really required. 
    1919""" 
    20 import numpy as np 
     20import numpy as np  # type: ignore 
    2121 
    2222def align_empty(shape, dtype, alignment=128): 

  • sasmodels/bumps_model.py

    rea75043 r7ae2b7f  
    1414import warnings 
    1515 
    16 import numpy as np 
     16import numpy as np  # type: ignore 
    1717 
    1818from .data import plot_theory 
    … …  
    7979    # lazy import; this allows the doc builder and nosetests to run even 
    8080    # when bumps is not on the path. 
    81     from bumps.names import Parameter 
    82  
    83     pars = {} 
    84     for p in model_info['parameters']: 
     81    from bumps.names import Parameter  # type: ignore 
     82 
     83    pars = {}     # floating point parameters 
     84    pd_types = {} # distribution names 
     85    for p in model_info.parameters.call_parameters: 
    8586        value = kwargs.pop(p.name, p.default) 
    8687        pars[p.name] = Parameter.default(value, name=p.name, limits=p.limits) 
    87     for name in model_info['partype']['pd-2d']: 
    88         for xpart, xdefault, xlimits in [ 
    89                 ('_pd', 0., pars[name].limits), 
    90                 ('_pd_n', 35., (0, 1000)), 
    91                 ('_pd_nsigma', 3., (0, 10)), 
    92             ]: 
    93             xname = name + xpart 
    94             xvalue = kwargs.pop(xname, xdefault) 
    95             pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits) 
    96  
    97     pd_types = {} 
    98     for name in model_info['partype']['pd-2d']: 
    99         xname = name + '_pd_type' 
    100         xvalue = kwargs.pop(xname, 'gaussian') 
    101         pd_types[xname] = xvalue 
     88        if p.polydisperse: 
     89            for part, default, limits in [ 
     90                    ('_pd', 0., pars[p.name].limits), 
     91                    ('_pd_n', 35., (0, 1000)), 
     92                    ('_pd_nsigma', 3., (0, 10)), 
     93                ]: 
     94                name = p.name + part 
     95                value = kwargs.pop(name, default) 
     96                pars[name] = Parameter.default(value, name=name, limits=limits) 
     97            name = p.name + '_pd_type' 
     98            pd_types[name] = kwargs.pop(name, 'gaussian') 
    10299 
    103100    if kwargs:  # args not corresponding to parameters 

  • sasmodels/compare.py

    rf247314 r7ae2b7f  
    3434import traceback 
    3535 
    36 import numpy as np 
     36import numpy as np  # type: ignore 
    3737 
    3838from . import core 
    3939from . import kerneldll 
    40 from . import product 
    4140from .data import plot_theory, empty_data1D, empty_data2D 
    4241from .direct_model import DirectModel 
    … …  
    210209    Choose a parameter range based on parameter name and initial value. 
    211210    """ 
     211    # process the polydispersity options 
    212212    if p.endswith('_pd_n'): 
    213213        return [0, 100] 
    … …  
    222222        else: 
    223223            return [-180, 180] 
     224    elif p.endswith('_pd'): 
     225        return [0, 1] 
    224226    elif 'sld' in p: 
    225227        return [-0.5, 10] 
    226     elif p.endswith('_pd'): 
    227         return [0, 1] 
    228228    elif p == 'background': 
    229229        return [0, 10] 
    230230    elif p == 'scale': 
    231231        return [0, 1e3] 
    232     elif p == 'case_num': 
    233         # RPA hack 
    234         return [0, 10] 
    235232    elif v < 0: 
    236         # Kxy parameters in rpa model can be negative 
    237233        return [2*v, -2*v] 
    238234    else: 
    … …  
    240236 
    241237 
    242 def _randomize_one(p, v): 
     238def _randomize_one(model_info, p, v): 
    243239    """ 
    244240    Randomize a single parameter. 
    … …  
    246242    if any(p.endswith(s) for s in ('_pd_n', '_pd_nsigma', '_pd_type')): 
    247243        return v 
     244 
     245    # Find the parameter definition 
     246    for par in model_info.parameters.call_parameters: 
     247        if par.name == p: 
     248            break 
    248249    else: 
    249         return np.random.uniform(*parameter_range(p, v)) 
    250  
    251  
    252 def randomize_pars(pars, seed=None): 
     250        raise ValueError("unknown parameter %r"%p) 
     251 
     252    if len(par.limits) > 2:  # choice list 
     253        return np.random.randint(len(par.limits)) 
     254 
     255    limits = par.limits 
     256    if not np.isfinite(limits).all(): 
     257        low, high = parameter_range(p, v) 
     258        limits = (max(limits[0], low), min(limits[1], high)) 
     259 
     260    return np.random.uniform(*limits) 
     261 
     262 
     263def randomize_pars(model_info, pars, seed=None): 
    253264    """ 
    254265    Generate random values for all of the parameters. 
    … …  
    261272    with push_seed(seed): 
    262273        # Note: the sort guarantees order `of calls to random number generator 
    263         pars = dict((p, _randomize_one(p, v)) 
     274        pars = dict((p, _randomize_one(model_info, p, v)) 
    264275                    for p, v in sorted(pars.items())) 
    265276    return pars 
    … …  
    273284    cylinder radius in this case). 
    274285    """ 
    275     name = model_info['id'] 
     286    name = model_info.id 
    276287    # if it is a product model, then just look at the form factor since 
    277288    # none of the structure factors need any constraints. 
    … …  
    294305        # Make sure phi sums to 1.0 
    295306        if pars['case_num'] < 2: 
    296             pars['Phia'] = 0. 
    297             pars['Phib'] = 0. 
     307            pars['Phi1'] = 0. 
     308            pars['Phi2'] = 0. 
    298309        elif pars['case_num'] < 5: 
    299             pars['Phia'] = 0. 
    300         total = sum(pars['Phi'+c] for c in 'abcd') 
    301         for c in 'abcd': 
     310            pars['Phi1'] = 0. 
     311        total = sum(pars['Phi'+c] for c in '1234') 
     312        for c in '1234': 
    302313            pars['Phi'+c] /= total 
    303314 
    … …  
    306317    Format the parameter list for printing. 
    307318    """ 
    308     if is2d: 
    309         exclude = lambda n: False 
    310     else: 
    311         partype = model_info['partype'] 
    312         par1d = set(partype['fixed-1d']+partype['pd-1d']) 
    313         exclude = lambda n: n not in par1d 
    314319    lines = [] 
    315     for p in model_info['parameters']: 
    316         if exclude(p.name): continue 
     320    parameters = model_info.parameters 
     321    for p in parameters.user_parameters(pars, is2d): 
    317322        fields = dict( 
    318             value=pars.get(p.name, p.default), 
    319             pd=pars.get(p.name+"_pd", 0.), 
    320             n=int(pars.get(p.name+"_pd_n", 0)), 
    321             nsigma=pars.get(p.name+"_pd_nsgima", 3.), 
    322             type=pars.get(p.name+"_pd_type", 'gaussian')) 
     323            value=pars.get(p.id, p.default), 
     324            pd=pars.get(p.id+"_pd", 0.), 
     325            n=int(pars.get(p.id+"_pd_n", 0)), 
     326            nsigma=pars.get(p.id+"_pd_nsgima", 3.), 
     327            type=pars.get(p.id+"_pd_type", 'gaussian')) 
    323328        lines.append(_format_par(p.name, **fields)) 
    324329    return "\n".join(lines) 
    … …  
    363368 
    364369    # grab the sasview model, or create it if it is a product model 
    365     if model_info['composition']: 
    366         composition_type, parts = model_info['composition'] 
     370    if model_info.composition: 
     371        composition_type, parts = model_info.composition 
    367372        if composition_type == 'product': 
    368373            from sas.models.MultiplicationModel import MultiplicationModel 
    … …  
    454459    """ 
    455460    # initialize the code so time is more accurate 
    456     value = calculator(**suppress_pd(pars)) 
     461    if Nevals > 1: 
     462        value = calculator(**suppress_pd(pars)) 
    457463    toc = tic() 
    458464    for _ in range(max(Nevals, 1)):  # make sure there is at least one eval 
    … …  
    661667    """ 
    662668    # Get the default values for the parameters 
    663     pars = dict((p.name, p.default) for p in model_info['parameters']) 
    664  
    665     # Fill in default values for the polydispersity parameters 
    666     for p in model_info['parameters']: 
    667         if p.type in ('volume', 'orientation'): 
    668             pars[p.name+'_pd'] = 0.0 
    669             pars[p.name+'_pd_n'] = 0 
    670             pars[p.name+'_pd_nsigma'] = 3.0 
    671             pars[p.name+'_pd_type'] = "gaussian" 
     669    pars = {} 
     670    for p in model_info.parameters.call_parameters: 
     671        parts = [('', p.default)] 
     672        if p.polydisperse: 
     673            parts.append(('_pd', 0.0)) 
     674            parts.append(('_pd_n', 0)) 
     675            parts.append(('_pd_nsigma', 3.0)) 
     676            parts.append(('_pd_type', "gaussian")) 
     677        for ext, val in parts: 
     678            if p.length > 1: 
     679                dict(("%s%d%s"%(p.id,k,ext), val) for k in range(p.length)) 
     680            else: 
     681                pars[p.id+ext] = val 
    672682 
    673683    # Plug in values given in demo 
    674684    if use_demo: 
    675         pars.update(model_info['demo']) 
     685        pars.update(model_info.demo) 
    676686    return pars 
    677687 
    … …  
    700710    try: 
    701711        model_info = core.load_model_info(name) 
    702     except ImportError, exc: 
     712    except ImportError as exc: 
    703713        print(str(exc)) 
    704714        print("Could not find model; use one of:\n    " + models) 
    … …  
    774784 
    775785    if len(engines) == 0: 
    776         engines.extend(['single', 'sasview']) 
     786        engines.extend(['single', 'double']) 
    777787    elif len(engines) == 1: 
    778         if engines[0][0] != 'sasview': 
    779             engines.append('sasview') 
     788        if engines[0][0] != 'double': 
     789            engines.append('double') 
    780790        else: 
    781791            engines.append('single') 
    … …  
    807817    #pars.update(set_pars)  # set value before random to control range 
    808818    if opts['seed'] > -1: 
    809         pars = randomize_pars(pars, seed=opts['seed']) 
     819        pars = randomize_pars(model_info, pars, seed=opts['seed']) 
    810820        print("Randomize using -random=%i"%opts['seed']) 
    811821    if opts['mono']: 
    … …  
    850860    Explore the model using the Bumps GUI. 
    851861    """ 
    852     import wx 
    853     from bumps.names import FitProblem 
    854     from bumps.gui.app_frame import AppFrame 
     862    import wx  # type: ignore 
     863    from bumps.names import FitProblem  # type: ignore 
     864    from bumps.gui.app_frame import AppFrame  # type: ignore 
    855865 
    856866    problem = FitProblem(Explore(opts)) 
    … …  
    873883    """ 
    874884    def __init__(self, opts): 
    875         from bumps.cli import config_matplotlib 
     885        from bumps.cli import config_matplotlib  # type: ignore 
    876886        from . import bumps_model 
    877887        config_matplotlib() 
    … …  
    879889        model_info = opts['def'] 
    880890        pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars']) 
     891        # Initialize parameter ranges, fixing the 2D parameters for 1D data. 
    881892        if not opts['is2d']: 
    882             active = [base + ext 
    883                       for base in model_info['partype']['pd-1d'] 
    884                       for ext in ['', '_pd', '_pd_n', '_pd_nsigma']] 
    885             active.extend(model_info['partype']['fixed-1d']) 
    886             for k in active: 
    887                 v = pars[k] 
    888                 v.range(*parameter_range(k, v.value)) 
     893            for p in model_info.parameters.user_parameters(is2d=False): 
     894                for ext in ['', '_pd', '_pd_n', '_pd_nsigma']: 
     895                    k = p.name+ext 
     896                    v = pars.get(k, None) 
     897                    if v is not None: 
     898                        v.range(*parameter_range(k, v.value)) 
    889899        else: 
    890900            for k, v in pars.items(): 

  • sasmodels/compare_many.py

    rce346b6 r7ae2b7f  
    1717import traceback 
    1818 
    19 import numpy as np 
     19import numpy as np  # type: ignore 
    2020 
    2121from . import core 
    22 from . import generate 
    2322from .compare import (MODELS, randomize_pars, suppress_pd, make_data, 
    2423                      make_engine, get_pars, columnize, 
    … …  
    109108 
    110109    if is_2d: 
    111         if not model_info['has_2d']: 
     110        if not model_info['parameters'].has_2d: 
    112111            print(',"1-D only"') 
    113112            return 
    … …  
    125124        try: 
    126125            result = fn(**pars) 
    127         except KeyboardInterrupt: 
    128             raise 
    129         except: 
     126        except Exception: 
    130127            traceback.print_exc() 
    131128            print("when comparing %s for %d"%(name, seed)) 
    … …  
    246243        base = sys.argv[5] 
    247244        comp = sys.argv[6] if len(sys.argv) > 6 else "sasview" 
    248     except: 
     245    except Exception: 
    249246        traceback.print_exc() 
    250247        print_usage() 

  • sasmodels/convert.json

    rec45c4f r6e7ff6d  
    602602    "RPAModel",  
    603603    { 
     604      "N1": "Na", "N2": "Nb", "N3": "Nc", "N4": "Nd", 
     605      "L1": "La", "L2": "Lb", "L3": "Lc", "L4": "Ld", 
     606      "v1": "va", "v2": "vb", "v3": "vc", "v4": "vd", 
     607      "b1": "ba", "b2": "bb", "b3": "bc", "b4": "bd", 
     608      "Phi1": "Phia", "Phi2": "Phib", "Phi3": "Phic", "Phi4": "Phid", 
    604609      "case_num": "lcase_n" 
    605610    } 
    … …  
    620625    } 
    621626  ],  
     627  "_spherepy": [ 
     628    "SphereModel", 
     629    { 
     630      "sld": "sldSph", 
     631      "radius": "radius", 
     632      "sld_solvent": "sldSolv" 
     633    } 
     634  ], 
    622635  "spherical_sld": [ 
    623636    "SphereSLDModel",  
    624637    { 
     638      "n": "n_shells", 
    625639      "radius_core": "rad_core",  
    626640      "sld_solvent": "sld_solv" 

  • sasmodels/convert.py

    rf247314 r7ae2b7f  
    6262    # but only if we haven't already byteified it 
    6363    if isinstance(data, dict) and not ignore_dicts: 
    64         return { 
    65             _byteify(key, ignore_dicts=True): _byteify(value, ignore_dicts=True) 
    66             for key, value in data.iteritems() 
    67         } 
     64        return dict((_byteify(key, ignore_dicts=True), 
     65                     _byteify(value, ignore_dicts=True)) 
     66                    for key, value in data.items()) 
    6867    # if it's anything else, return it in its original form 
    6968    return data 
    … …  
    153152def revert_name(model_info): 
    154153    _read_conversion_table() 
    155     oldname, oldpars = CONVERSION_TABLE.get(model_info['id'], [None, {}]) 
     154    oldname, oldpars = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
    156155    return oldname 
    157156 
    158157def _get_old_pars(model_info): 
    159158    _read_conversion_table() 
    160     oldname, oldpars = CONVERSION_TABLE.get(model_info['id'], [None, {}]) 
     159    oldname, oldpars = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
    161160    return oldpars 
    162161 
    … …  
    165164    Convert model from new style parameter names to old style. 
    166165    """ 
    167     if model_info['composition'] is not None: 
    168         composition_type, parts = model_info['composition'] 
     166    if model_info.composition is not None: 
     167        composition_type, parts = model_info.composition 
    169168        if composition_type == 'product': 
    170169            P, S = parts 
    … …  
    180179 
    181180    # Note: update compare.constrain_pars to match 
    182     name = model_info['id'] 
    183     if name in MODELS_WITHOUT_SCALE or model_info['structure_factor']: 
     181    name = model_info.id 
     182    if name in MODELS_WITHOUT_SCALE or model_info.structure_factor: 
    184183        if oldpars.pop('scale', 1.0) != 1.0: 
    185184            warnings.warn("parameter scale not used in sasview %s"%name) 
    186     if name in MODELS_WITHOUT_BACKGROUND or model_info['structure_factor']: 
     185    if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor: 
    187186        if oldpars.pop('background', 0.0) != 0.0: 
    188187            warnings.warn("parameter background not used in sasview %s"%name) 
    … …  
    206205        elif name == 'rpa': 
    207206            # convert scattering lengths from femtometers to centimeters 
    208             for p in "La", "Lb", "Lc", "Ld": 
     207            for p in "L1", "L2", "L3", "L4": 
    209208                if p in oldpars: oldpars[p] *= 1e-13 
    210209        elif name == 'core_shell_parallelepiped': 
    … …  
    225224    Restrict parameter values to those that will match sasview. 
    226225    """ 
    227     name = model_info['id'] 
     226    name = model_info.id 
    228227    # Note: update convert.revert_model to match 
    229     if name in MODELS_WITHOUT_SCALE or model_info['structure_factor']: 
     228    if name in MODELS_WITHOUT_SCALE or model_info.structure_factor: 
    230229        pars['scale'] = 1 
    231     if name in MODELS_WITHOUT_BACKGROUND or model_info['structure_factor']: 
     230    if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor: 
    232231        pars['background'] = 0 
    233232    # sasview multiplies background by structure factor 

  • sasmodels/core.py

    r4d76711 rfa5fd8d  
    22Core model handling routines. 
    33""" 
     4from __future__ import print_function 
     5 
    46__all__ = [ 
    5     "list_models", "load_model_info", "precompile_dll", 
    6     "build_model", "make_kernel", "call_kernel", "call_ER_VR", 
     7    "list_models", "load_model", "load_model_info", 
     8    "build_model", "precompile_dll", 
    79    ] 
    810 
    9 from os.path import basename, dirname, join as joinpath, splitext 
     11from os.path import basename, dirname, join as joinpath 
    1012from glob import glob 
    1113 
    12 import numpy as np 
     14import numpy as np # type: ignore 
    1315 
    14 from . import models 
    15 from . import weights 
    1616from . import generate 
    17 # TODO: remove circular references between product and core 
    18 # product uses call_ER/call_VR, core uses make_product_info/ProductModel 
    19 #from . import product 
     17from . import modelinfo 
     18from . import product 
    2019from . import mixture 
    2120from . import kernelpy 
    … …  
    2423    from . import kernelcl 
    2524    HAVE_OPENCL = True 
    26 except: 
     25except Exception: 
    2726    HAVE_OPENCL = False 
    2827 
    2928try: 
    30     np.meshgrid([]) 
    31     meshgrid = np.meshgrid 
    32 except ValueError: 
    33     # CRUFT: np.meshgrid requires multiple vectors 
    34     def meshgrid(*args): 
    35         if len(args) > 1: 
    36             return np.meshgrid(*args) 
    37         else: 
    38             return [np.asarray(v) for v in args] 
     29    from typing import List, Union, Optional, Any 
     30    DType = Union[None, str, np.dtype] 
     31    from .kernel import KernelModel 
     32except ImportError: 
     33    pass 
     34 
    3935 
    4036# TODO: refactor composite model support 
    … …  
    5147 
    5248def list_models(): 
     49    # type: () -> List[str] 
    5350    """ 
    5451    Return the list of available models on the model path. 
    … …  
    6057 
    6158def isstr(s): 
     59    # type: (Any) -> bool 
    6260    """ 
    6361    Return True if *s* is a string-like object. 
    6462    """ 
    6563    try: s + '' 
    66     except: return False 
     64    except Exception: return False 
    6765    return True 
    6866 
    69 def load_model(model_name, **kw): 
     67def load_model(model_name, dtype=None, platform='ocl'): 
     68    # type: (str, DType, str) -> KernelModel 
    7069    """ 
    7170    Load model info and build model. 
     71 
     72    *model_name* is the name of the model as used by :func:`load_model_info`. 
     73    Additional keyword arguments are passed directly to :func:`build_model`. 
    7274    """ 
    73     return build_model(load_model_info(model_name), **kw) 
     75    return build_model(load_model_info(model_name), 
     76                       dtype=dtype, platform=platform) 
    7477 
    7578 
    7679def load_model_info(model_name): 
     80    # type: (str) -> modelinfo.ModelInfo 
    7781    """ 
    7882    Load a model definition given the model name. 
    … …  
    8892    parts = model_name.split('*') 
    8993    if len(parts) > 1: 
    90         from . import product 
    91         # Note: currently have circular reference 
    9294        if len(parts) > 2: 
    9395            raise ValueError("use P*S to apply structure factor S to model P") 
    … …  
    9698 
    9799    kernel_module = generate.load_kernel_module(model_name) 
    98     return generate.make_model_info(kernel_module) 
     100    return modelinfo.make_model_info(kernel_module) 
    99101 
    100102 
    101103def build_model(model_info, dtype=None, platform="ocl"): 
     104    # type: (modelinfo.ModelInfo, np.dtype, str) -> KernelModel 
    102105    """ 
    103106    Prepare the model for the default execution platform. 
    … …  
    118121    otherwise it uses the default "ocl". 
    119122    """ 
    120     composition = model_info.get('composition', None) 
     123    composition = model_info.composition 
    121124    if composition is not None: 
    122125        composition_type, parts = composition 
    … …  
    131134            raise ValueError('unknown mixture type %s'%composition_type) 
    132135 
     136    # If it is a python model, return it immediately 
     137    if callable(model_info.Iq): 
     138        return kernelpy.PyModel(model_info) 
     139 
    133140    ## for debugging: 
    134141    ##  1. uncomment open().write so that the source will be saved next time 
    … …  
    137144    ##  4. rerun "python -m sasmodels.direct_model $MODELNAME" 
    138145    ##  5. uncomment open().read() so that source will be regenerated from model 
    139     # open(model_info['name']+'.c','w').write(source) 
    140     # source = open(model_info['name']+'.cl','r').read() 
     146    # open(model_info.name+'.c','w').write(source) 
     147    # source = open(model_info.name+'.cl','r').read() 
    141148    source = generate.make_source(model_info) 
    142149    if dtype is None: 
    143         dtype = 'single' if model_info['single'] else 'double' 
    144     if callable(model_info.get('Iq', None)): 
    145         return kernelpy.PyModel(model_info) 
     150        dtype = generate.F32 if model_info.single else generate.F64 
    146151    if (platform == "dll" 
    147152            or not HAVE_OPENCL 
    … …  
    152157 
    153158def precompile_dll(model_name, dtype="double"): 
     159    # type: (str, DType) -> Optional[str] 
    154160    """ 
    155161    Precompile the dll for a model. 
    … …  
    168174    source = generate.make_source(model_info) 
    169175    return kerneldll.make_dll(source, model_info, dtype=dtype) if source else None 
    170  
    171  
    172 def get_weights(model_info, pars, name): 
    173     """ 
    174     Generate the distribution for parameter *name* given the parameter values 
    175     in *pars*. 
    176  
    177     Uses "name", "name_pd", "name_pd_type", "name_pd_n", "name_pd_sigma" 
    178     from the *pars* dictionary for parameter value and parameter dispersion. 
    179     """ 
    180     relative = name in model_info['partype']['pd-rel'] 
    181     limits = model_info['limits'][name] 
    182     disperser = pars.get(name+'_pd_type', 'gaussian') 
    183     value = pars.get(name, model_info['defaults'][name]) 
    184     npts = pars.get(name+'_pd_n', 0) 
    185     width = pars.get(name+'_pd', 0.0) 
    186     nsigma = pars.get(name+'_pd_nsigma', 3.0) 
    187     value, weight = weights.get_weights( 
    188         disperser, npts, width, nsigma, value, limits, relative) 
    189     return value, weight / np.sum(weight) 
    190  
    191 def dispersion_mesh(pars): 
    192     """ 
    193     Create a mesh grid of dispersion parameters and weights. 
    194  
    195     Returns [p1,p2,...],w where pj is a vector of values for parameter j 
    196     and w is a vector containing the products for weights for each 
    197     parameter set in the vector. 
    198     """ 
    199     value, weight = zip(*pars) 
    200     value = [v.flatten() for v in meshgrid(*value)] 
    201     weight = np.vstack([v.flatten() for v in meshgrid(*weight)]) 
    202     weight = np.prod(weight, axis=0) 
    203     return value, weight 
    204  
    205 def call_kernel(kernel, pars, cutoff=0, mono=False): 
    206     """ 
    207     Call *kernel* returned from *model.make_kernel* with parameters *pars*. 
    208  
    209     *cutoff* is the limiting value for the product of dispersion weights used 
    210     to perform the multidimensional dispersion calculation more quickly at a 
    211     slight cost to accuracy. The default value of *cutoff=0* integrates over 
    212     the entire dispersion cube.  Using *cutoff=1e-5* can be 50% faster, but 
    213     with an error of about 1%, which is usually less than the measurement 
    214     uncertainty. 
    215  
    216     *mono* is True if polydispersity should be set to none on all parameters. 
    217     """ 
    218     fixed_pars = [pars.get(name, kernel.info['defaults'][name]) 
    219                   for name in kernel.fixed_pars] 
    220     if mono: 
    221         pd_pars = [( np.array([pars[name]]), np.array([1.0]) ) 
    222                    for name in kernel.pd_pars] 
    223     else: 
    224         pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars] 
    225     return kernel(fixed_pars, pd_pars, cutoff=cutoff) 
    226  
    227 def call_ER_VR(model_info, vol_pars): 
    228     """ 
    229     Return effect radius and volume ratio for the model. 
    230  
    231     *info* is either *kernel.info* for *kernel=make_kernel(model,q)* 
    232     or *model.info*. 
    233  
    234     *pars* are the parameters as expected by :func:`call_kernel`. 
    235     """ 
    236     ER = model_info.get('ER', None) 
    237     VR = model_info.get('VR', None) 
    238     value, weight = dispersion_mesh(vol_pars) 
    239  
    240     individual_radii = ER(*value) if ER else 1.0 
    241     whole, part = VR(*value) if VR else (1.0, 1.0) 
    242  
    243     effect_radius = np.sum(weight*individual_radii) / np.sum(weight) 
    244     volume_ratio = np.sum(weight*part)/np.sum(weight*whole) 
    245     return effect_radius, volume_ratio 
    246  
    247  
    248 def call_ER(model_info, values): 
    249     """ 
    250     Call the model ER function using *values*. *model_info* is either 
    251     *model.info* if you have a loaded model, or *kernel.info* if you 
    252     have a model kernel prepared for evaluation. 
    253     """ 
    254     ER = model_info.get('ER', None) 
    255     if ER is None: 
    256         return 1.0 
    257     else: 
    258         vol_pars = [get_weights(model_info, values, name) 
    259                     for name in model_info['partype']['volume']] 
    260         value, weight = dispersion_mesh(vol_pars) 
    261         individual_radii = ER(*value) 
    262         #print(values[0].shape, weights.shape, fv.shape) 
    263         return np.sum(weight*individual_radii) / np.sum(weight) 
    264  
    265 def call_VR(model_info, values): 
    266     """ 
    267     Call the model VR function using *pars*. 
    268     *info* is either *model.info* if you have a loaded model, or *kernel.info* 
    269     if you have a model kernel prepared for evaluation. 
    270     """ 
    271     VR = model_info.get('VR', None) 
    272     if VR is None: 
    273         return 1.0 
    274     else: 
    275         vol_pars = [get_weights(model_info, values, name) 
    276                     for name in model_info['partype']['volume']] 
    277         value, weight = dispersion_mesh(vol_pars) 
    278         whole, part = VR(*value) 
    279         return np.sum(weight*part)/np.sum(weight*whole) 
    280  
    281 # TODO: remove call_ER, call_VR 
    282  

  • sasmodels/custom/__init__.py

    rcd0a808 r7ae2b7f  
    1414try: 
    1515    # Python 3.5 and up 
    16     from importlib.util import spec_from_file_location, module_from_spec 
     16    from importlib.util import spec_from_file_location, module_from_spec  # type: ignore 
    1717    def load_module_from_path(fullname, path): 
    1818        spec = spec_from_file_location(fullname, path) 

  • sasmodels/data.py

    rd6f5da6 r7ae2b7f  
    3535import traceback 
    3636 
    37 import numpy as np 
     37import numpy as np  # type: ignore 
    3838 
    3939def load_data(filename): 
    … …  
    4141    Load data using a sasview loader. 
    4242    """ 
    43     from sas.sascalc.dataloader.loader import Loader 
     43    from sas.sascalc.dataloader.loader import Loader  # type: ignore 
    4444    loader = Loader() 
    4545    data = loader.load(filename) 
    … …  
    5353    Add a beam stop of the given *radius*.  If *outer*, make an annulus. 
    5454    """ 
    55     from sas.dataloader.manipulations import Ringcut 
     55    from sas.dataloader.manipulations import Ringcut  # type: ignore 
    5656    if hasattr(data, 'qx_data'): 
    5757        data.mask = Ringcut(0, radius)(data) 
    … …  
    6868    Select half of the data, either "right" or "left". 
    6969    """ 
    70     from sas.dataloader.manipulations import Boxcut 
     70    from sas.dataloader.manipulations import Boxcut  # type: ignore 
    7171    if half == 'right': 
    7272        data.mask += \ 
    … …  
    8181    Chop the top off the data, above *cutoff*. 
    8282    """ 
    83     from sas.dataloader.manipulations import Boxcut 
     83    from sas.dataloader.manipulations import Boxcut  # type: ignore 
    8484    data.mask += \ 
    8585        Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=cutoff)(data) 
    … …  
    358358        try: 
    359359            return fn(*args, **kw) 
    360         except KeyboardInterrupt: 
    361             raise 
    362         except: 
     360        except Exception: 
    363361            traceback.print_exc() 
    364362 
    … …  
    372370    Plot the data and residuals for 1D data. 
    373371    """ 
    374     import matplotlib.pyplot as plt 
    375     from numpy.ma import masked_array, masked 
     372    import matplotlib.pyplot as plt  # type: ignore 
     373    from numpy.ma import masked_array, masked  # type: ignore 
    376374 
    377375    use_data = use_data and data.y is not None 
    … …  
    449447    Plot SESANS results. 
    450448    """ 
    451     import matplotlib.pyplot as plt 
     449    import matplotlib.pyplot as plt  # type: ignore 
    452450    use_data = use_data and data.y is not None 
    453451    use_theory = theory is not None 
    … …  
    492490    Plot the data and residuals for 2D data. 
    493491    """ 
    494     import matplotlib.pyplot as plt 
     492    import matplotlib.pyplot as plt  # type: ignore 
    495493    use_data = use_data and data.data is not None 
    496494    use_theory = theory is not None 
    … …  
    554552    *scale* can be 'log' for log scale data, or 'linear'. 
    555553    """ 
    556     import matplotlib.pyplot as plt 
    557     from numpy.ma import masked_array 
     554    import matplotlib.pyplot as plt  # type: ignore 
     555    from numpy.ma import masked_array  # type: ignore 
    558556 
    559557    image = np.zeros_like(data.qx_data) 
    … …  
    595593    set_beam_stop(data, 0.004) 
    596594    plot_data(data) 
    597     import matplotlib.pyplot as plt; plt.show() 
     595    import matplotlib.pyplot as plt  # type: ignore 
     596    plt.show() 
    598597 
    599598 

  • sasmodels/direct_model.py

    r4d76711 r7ae2b7f  
    2323from __future__ import print_function 
    2424 
    25 import numpy as np 
    26  
    27 from .core import call_kernel, call_ER_VR 
    28 from . import sesans 
     25import numpy as np  # type: ignore 
     26 
     27# TODO: fix sesans module 
     28from . import sesans  # type: ignore 
     29from . import weights 
    2930from . import resolution 
    3031from . import resolution2d 
     32from . import details 
     33 
     34 
     35def call_kernel(kernel, pars, cutoff=0, mono=False): 
     36    """ 
     37    Call *kernel* returned from *model.make_kernel* with parameters *pars*. 
     38 
     39    *cutoff* is the limiting value for the product of dispersion weights used 
     40    to perform the multidimensional dispersion calculation more quickly at a 
     41    slight cost to accuracy. The default value of *cutoff=0* integrates over 
     42    the entire dispersion cube.  Using *cutoff=1e-5* can be 50% faster, but 
     43    with an error of about 1%, which is usually less than the measurement 
     44    uncertainty. 
     45 
     46    *mono* is True if polydispersity should be set to none on all parameters. 
     47    """ 
     48    parameters = kernel.info.parameters 
     49    if mono: 
     50        active = lambda name: False 
     51    elif kernel.dim == '1d': 
     52        active = lambda name: name in parameters.pd_1d 
     53    elif kernel.dim == '2d': 
     54        active = lambda name: name in parameters.pd_2d 
     55    else: 
     56        active = lambda name: True 
     57 
     58    vw_pairs = [(get_weights(p, pars) if active(p.name) 
     59                 else ([pars.get(p.name, p.default)], [])) 
     60                for p in parameters.call_parameters] 
     61 
     62    call_details, weights, values = details.build_details(kernel, vw_pairs) 
     63    return kernel(call_details, weights, values, cutoff) 
     64 
     65def get_weights(parameter, values): 
     66    """ 
     67    Generate the distribution for parameter *name* given the parameter values 
     68    in *pars*. 
     69 
     70    Uses "name", "name_pd", "name_pd_type", "name_pd_n", "name_pd_sigma" 
     71    from the *pars* dictionary for parameter value and parameter dispersion. 
     72    """ 
     73    value = float(values.get(parameter.name, parameter.default)) 
     74    relative = parameter.relative_pd 
     75    limits = parameter.limits 
     76    disperser = values.get(parameter.name+'_pd_type', 'gaussian') 
     77    npts = values.get(parameter.name+'_pd_n', 0) 
     78    width = values.get(parameter.name+'_pd', 0.0) 
     79    nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 
     80    if npts == 0 or width == 0: 
     81        return [value], [] 
     82    value, weight = weights.get_weights( 
     83        disperser, npts, width, nsigma, value, limits, relative) 
     84    return value, weight / np.sum(weight) 
    3185 
    3286class DataMixin(object): 
    … …  
    67121            self.data_type = 'Iq' 
    68122 
    69         partype = model.info['partype'] 
    70  
    71123        if self.data_type == 'sesans': 
    72124             
    … …  
    82134            q_mono = sesans.make_all_q(data) 
    83135        elif self.data_type == 'Iqxy': 
    84             if not partype['orientation'] and not partype['magnetic']: 
    85                 raise ValueError("not 2D without orientation or magnetic parameters") 
     136            #if not model.info.parameters.has_2d: 
     137            #    raise ValueError("not 2D without orientation or magnetic parameters") 
    86138            q = np.sqrt(data.qx_data**2 + data.qy_data**2) 
    87139            qmin = getattr(data, 'qmin', 1e-16) 
    … …  
    172224    def _calc_theory(self, pars, cutoff=0.0): 
    173225        if self._kernel is None: 
    174             self._kernel = self._model.make_kernel(self._kernel_inputs)  # pylint: disable=attribute-dedata_type 
     226            self._kernel = self._model.make_kernel(self._kernel_inputs) 
    175227            self._kernel_mono = (self._model.make_kernel(self._kernel_mono_inputs) 
    176228                                 if self._kernel_mono_inputs else None) 
    … …  
    214266        return self._calc_theory(pars, cutoff=self.cutoff) 
    215267 
    216     def ER_VR(self, **pars): 
    217         """ 
    218         Compute the equivalent radius and volume ratio for the model. 
    219         """ 
    220         return call_ER_VR(self.model.info, pars) 
    221  
    222268    def simulate_data(self, noise=None, **pars): 
    223269        """ 
    … …  
    243289        try: 
    244290            values = [float(v) for v in call.split(',')] 
    245         except: 
     291        except Exception: 
    246292            values = [] 
    247293        if len(values) == 1: 

  • sasmodels/generate.py

    r81ec7c8 r81ec7c8  
    2121 
    2222    *VR(p1, p2, ...)* returns the volume ratio for core-shell style forms. 
     23 
     24    #define INVALID(v) (expr)  returns False if v.parameter is invalid 
     25    for some parameter or other (e.g., v.bell_radius < v.radius).  If 
     26    necessary, the expression can call a function. 
    2327 
    2428These functions are defined in a kernel module .py script and an associated 
    … …  
    6569The constructor code will generate the necessary vectors for computing 
    6670them with the desired polydispersity. 
    67  
    68 The available kernel parameters are defined as a list, with each parameter 
    69 defined as a sublist with the following elements: 
    70  
    71     *name* is the name that will be used in the call to the kernel 
    72     function and the name that will be displayed to the user.  Names 
    73     should be lower case, with words separated by underscore.  If 
    74     acronyms are used, the whole acronym should be upper case. 
    75  
    76     *units* should be one of *degrees* for angles, *Ang* for lengths, 
    77     *1e-6/Ang^2* for SLDs. 
    78  
    79     *default value* will be the initial value for  the model when it 
    80     is selected, or when an initial value is not otherwise specified. 
    81  
    82     *limits = [lb, ub]* are the hard limits on the parameter value, used to 
    83     limit the polydispersity density function.  In the fit, the parameter limits 
    84     given to the fit are the limits  on the central value of the parameter. 
    85     If there is polydispersity, it will evaluate parameter values outside 
    86     the fit limits, but not outside the hard limits specified in the model. 
    87     If there are no limits, use +/-inf imported from numpy. 
    88  
    89     *type* indicates how the parameter will be used.  "volume" parameters 
    90     will be used in all functions.  "orientation" parameters will be used 
    91     in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in 
    92     *Imagnetic* only.  If *type* is the empty string, the parameter will 
    93     be used in all of *Iq*, *Iqxy* and *Imagnetic*.  "sld" parameters 
    94     can automatically be promoted to magnetic parameters, each of which 
    95     will have a magnitude and a direction, which may be different from 
    96     other sld parameters. 
    97  
    98     *description* is a short description of the parameter.  This will 
    99     be displayed in the parameter table and used as a tool tip for the 
    100     parameter value in the user interface. 
    101  
    10271The kernel module must set variables defining the kernel meta data: 
    10372 
    … …  
    156125An *model_info* dictionary is constructed from the kernel meta data and 
    157126returned to the caller. 
    158  
    159 The model evaluator, function call sequence consists of q inputs and the return vector, 
    160 followed by the loop value/weight vector, followed by the values for 
    161 the non-polydisperse parameters, followed by the lengths of the 
    162 polydispersity loops.  To construct the call for 1D models, the 
    163 categories *fixed-1d* and *pd-1d* list the names of the parameters 
    164 of the non-polydisperse and the polydisperse parameters respectively. 
    165 Similarly, *fixed-2d* and *pd-2d* provide parameter names for 2D models. 
    166 The *pd-rel* category is a set of those parameters which give 
    167 polydispersitiy as a portion of the value (so a 10% length dispersity 
    168 would use a polydispersity value of 0.1) rather than absolute 
    169 dispersity such as an angle plus or minus 15 degrees. 
    170  
    171 The *volume* category lists the volume parameters in order for calls 
    172 to volume within the kernel (used for volume normalization) and for 
    173 calls to ER and VR for effective radius and volume ratio respectively. 
    174  
    175 The *orientation* and *magnetic* categories list the orientation and 
    176 magnetic parameters.  These are used by the sasview interface.  The 
    177 blank category is for parameters such as scale which don't have any 
    178 other marking. 
    179127 
    180128The doc string at the start of the kernel module will be used to 
    … …  
    204152from __future__ import print_function 
    205153 
    206 #TODO: identify model files which have changed since loading and reload them. 
    207154#TODO: determine which functions are useful outside of generate 
    208155#__all__ = ["model_info", "make_doc", "make_source", "convert_type"] 
    209156 
    210 import sys 
    211 from os.path import abspath, dirname, join as joinpath, exists, basename, \ 
    212     splitext 
     157from os.path import abspath, dirname, join as joinpath, exists, getmtime 
    213158import re 
    214159import string 
    215160import warnings 
    216 from collections import namedtuple 
    217  
    218 import numpy as np 
    219  
     161 
     162import numpy as np  # type: ignore 
     163 
     164from .modelinfo import Parameter 
    220165from .custom import load_custom_kernel_module 
    221166 
    222 PARAMETER_FIELDS = ['name', 'units', 'default', 'limits', 'type', 'description'] 
    223 Parameter = namedtuple('Parameter', PARAMETER_FIELDS) 
    224  
    225 C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c') 
     167try: 
     168    from typing import Tuple, Sequence, Iterator 
     169    from .modelinfo import ModelInfo 
     170except ImportError: 
     171    pass 
     172 
     173TEMPLATE_ROOT = dirname(__file__) 
    226174 
    227175F16 = np.dtype('float16') 
    … …  
    232180except TypeError: 
    233181    F128 = None 
    234  
    235 # Scale and background, which are parameters common to every form factor 
    236 COMMON_PARAMETERS = [ 
    237     ["scale", "", 1, [0, np.inf], "", "Source intensity"], 
    238     ["background", "1/cm", 1e-3, [0, np.inf], "", "Source background"], 
    239     ] 
    240182 
    241183# Conversion from units defined in the parameter table for each model 
    … …  
    284226 
    285227def format_units(units): 
     228    # type: (str) -> str 
    286229    """ 
    287230    Convert units into ReStructured Text format. 
    … …  
    290233 
    291234def make_partable(pars): 
     235    # type: (List[Parameter]) -> str 
    292236    """ 
    293237    Generate the parameter table to include in the sphinx documentation. 
    … …  
    320264 
    321265def _search(search_path, filename): 
     266    # type: (List[str], str) -> str 
    322267    """ 
    323268    Find *filename* in *search_path*. 
    … …  
    331276    raise ValueError("%r not found in %s" % (filename, search_path)) 
    332277 
     278 
    333279def model_sources(model_info): 
     280    # type: (ModelInfo) -> List[str] 
    334281    """ 
    335282    Return a list of the sources file paths for the module. 
    336283    """ 
    337     search_path = [dirname(model_info['filename']), 
     284    search_path = [dirname(model_info.filename), 
    338285                   abspath(joinpath(dirname(__file__), 'models'))] 
    339     return [_search(search_path, f) for f in model_info['source']] 
    340  
    341 # Pragmas for enable OpenCL features.  Be sure to protect them so that they 
    342 # still compile even if OpenCL is not present. 
    343 _F16_PRAGMA = """\ 
    344 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 
    345 #  pragma OPENCL EXTENSION cl_khr_fp16: enable 
    346 #endif 
    347 """ 
    348  
    349 _F64_PRAGMA = """\ 
    350 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 
    351 #  pragma OPENCL EXTENSION cl_khr_fp64: enable 
    352 #endif 
    353 """ 
     286    return [_search(search_path, f) for f in model_info.source] 
     287 
     288def timestamp(model_info): 
     289    # type: (ModelInfo) -> int 
     290    """ 
     291    Return a timestamp for the model corresponding to the most recently 
     292    changed file or dependency. 
     293    """ 
     294    source_files = (model_sources(model_info) 
     295                    + model_templates() 
     296                    + [model_info.filename]) 
     297    newest = max(getmtime(f) for f in source_files) 
     298    return newest 
    354299 
    355300def convert_type(source, dtype): 
     301    # type: (str, np.dtype) -> str 
    356302    """ 
    357303    Convert code from double precision to the desired type. 
    … …  
    362308    if dtype == F16: 
    363309        fbytes = 2 
    364         source = _F16_PRAGMA + _convert_type(source, "half", "f") 
     310        source = _convert_type(source, "half", "f") 
    365311    elif dtype == F32: 
    366312        fbytes = 4 
    … …  
    368314    elif dtype == F64: 
    369315        fbytes = 8 
    370         source = _F64_PRAGMA + source  # Source is already double 
     316        # no need to convert if it is already double 
    371317    elif dtype == F128: 
    372318        fbytes = 16 
    … …  
    378324 
    379325def _convert_type(source, type_name, constant_flag): 
     326    # type: (str, str, str) -> str 
    380327    """ 
    381328    Replace 'double' with *type_name* in *source*, tagging floating point 
    … …  
    396343 
    397344def kernel_name(model_info, is_2d): 
     345    # type: (ModelInfo, bool) -> str 
    398346    """ 
    399347    Name of the exported kernel symbol. 
    400348    """ 
    401     return model_info['name'] + "_" + ("Iqxy" if is_2d else "Iq") 
     349    return model_info.name + "_" + ("Iqxy" if is_2d else "Iq") 
    402350 
    403351 
    404352def indent(s, depth): 
     353    # type: (str, int) -> str 
    405354    """ 
    406355    Indent a string of text with *depth* additional spaces on each line. 
    … …  
    411360 
    412361 
    413 LOOP_OPEN = """\ 
    414 for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) { 
    415   const double %(name)s = loops[2*(%(name)s_i%(offset)s)]; 
    416   const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\ 
     362_template_cache = {}  # type: Dict[str, Tuple[int, str, str]] 
     363def load_template(filename): 
     364    # type: (str) -> str 
     365    path = joinpath(TEMPLATE_ROOT, filename) 
     366    mtime = getmtime(path) 
     367    if filename not in _template_cache or mtime > _template_cache[filename][0]: 
     368        with open(path) as fid: 
     369            _template_cache[filename] = (mtime, fid.read(), path) 
     370    return _template_cache[filename][1] 
     371 
     372def model_templates(): 
     373    # type: () -> List[str] 
     374    # TODO: fails DRY; templates are listed in two places. 
     375    # should instead have model_info contain a list of paths 
     376    return [joinpath(TEMPLATE_ROOT, filename) 
     377            for filename in ('kernel_header.c', 'kernel_iq.c')] 
     378 
     379 
     380_FN_TEMPLATE = """\ 
     381double %(name)s(%(pars)s); 
     382double %(name)s(%(pars)s) { 
     383    %(body)s 
     384} 
     385 
     386 
    417387""" 
    418 def build_polydispersity_loops(pd_pars): 
    419     """ 
    420     Build polydispersity loops 
    421  
    422     Returns loop opening and loop closing 
    423     """ 
    424     depth = 4 
    425     offset = "" 
    426     loop_head = [] 
    427     loop_end = [] 
    428     for name in pd_pars: 
    429         subst = {'name': name, 'offset': offset} 
    430         loop_head.append(indent(LOOP_OPEN % subst, depth)) 
    431         loop_end.insert(0, (" "*depth) + "}") 
    432         offset += '+N' + name 
    433         depth += 2 
    434     return "\n".join(loop_head), "\n".join(loop_end) 
    435  
    436 C_KERNEL_TEMPLATE = None 
     388 
     389def _gen_fn(name, pars, body): 
     390    # type: (str, List[Parameter], str) -> str 
     391    """ 
     392    Generate a function given pars and body. 
     393 
     394    Returns the following string:: 
     395 
     396         double fn(double a, double b, ...); 
     397         double fn(double a, double b, ...) { 
     398             .... 
     399         } 
     400    """ 
     401    par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void' 
     402    return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl} 
     403 
     404def _call_pars(prefix, pars): 
     405    # type: (str, List[Parameter]) -> List[str] 
     406    """ 
     407    Return a list of *prefix.parameter* from parameter items. 
     408    """ 
     409    return [p.as_call_reference(prefix) for p in pars] 
     410 
     411_IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)", 
     412                           flags=re.MULTILINE) 
     413def _have_Iqxy(sources): 
     414    # type: (List[str]) -> bool 
     415    """ 
     416    Return true if any file defines Iqxy. 
     417 
     418    Note this is not a C parser, and so can be easily confused by 
     419    non-standard syntax.  Also, it will incorrectly identify the following 
     420    as having Iqxy:: 
     421 
     422        /* 
     423        double Iqxy(qx, qy, ...) { ... fill this in later ... } 
     424        */ 
     425 
     426    If you want to comment out an Iqxy function, use // on the front of the 
     427    line instead. 
     428    """ 
     429    for code in sources: 
     430        if _IQXY_PATTERN.search(code): 
     431            return True 
     432    else: 
     433        return False 
     434 
    437435def make_source(model_info): 
     436    # type: (ModelInfo) -> str 
    438437    """ 
    439438    Generate the OpenCL/ctypes kernel from the module info. 
    440439 
    441     Uses source files found in the given search path. 
    442     """ 
    443     if callable(model_info['Iq']): 
    444         return None 
     440    Uses source files found in the given search path.  Returns None if this 
     441    is a pure python model, with no C source components. 
     442    """ 
     443    if callable(model_info.Iq): 
     444        raise ValueError("can't compile python model") 
    445445 
    446446    # TODO: need something other than volume to indicate dispersion parameters 
    … …  
    453453    # for computing volume even if we allow non-disperse volume parameters. 
    454454 
    455     # Load template 
    456     global C_KERNEL_TEMPLATE 
    457     if C_KERNEL_TEMPLATE is None: 
    458         with open(C_KERNEL_TEMPLATE_PATH) as fid: 
    459             C_KERNEL_TEMPLATE = fid.read() 
    460  
    461     # Load additional sources 
    462     source = [open(f).read() for f in model_sources(model_info)] 
    463  
    464     # Prepare defines 
    465     defines = [] 
    466     partype = model_info['partype'] 
    467     pd_1d = partype['pd-1d'] 
    468     pd_2d = partype['pd-2d'] 
    469     fixed_1d = partype['fixed-1d'] 
    470     fixed_2d = partype['fixed-1d'] 
    471  
    472     iq_parameters = [p.name 
    473                      for p in model_info['parameters'][2:]  # skip scale, background 
    474                      if p.name in set(fixed_1d + pd_1d)] 
    475     iqxy_parameters = [p.name 
    476                        for p in model_info['parameters'][2:]  # skip scale, background 
    477                        if p.name in set(fixed_2d + pd_2d)] 
    478     volume_parameters = [p.name 
    479                          for p in model_info['parameters'] 
    480                          if p.type == 'volume'] 
    481  
    482     # Fill in defintions for volume parameters 
    483     if volume_parameters: 
    484         defines.append(('VOLUME_PARAMETERS', 
    485                         ','.join(volume_parameters))) 
    486         defines.append(('VOLUME_WEIGHT_PRODUCT', 
    487                         '*'.join(p + '_w' for p in volume_parameters))) 
    488  
    489     # Generate form_volume function from body only 
    490     if model_info['form_volume'] is not None: 
    491         if volume_parameters: 
    492             vol_par_decl = ', '.join('double ' + p for p in volume_parameters) 
    493         else: 
    494             vol_par_decl = 'void' 
    495         defines.append(('VOLUME_PARAMETER_DECLARATIONS', 
    496                         vol_par_decl)) 
    497         fn = """\ 
    498 double form_volume(VOLUME_PARAMETER_DECLARATIONS); 
    499 double form_volume(VOLUME_PARAMETER_DECLARATIONS) { 
    500     %(body)s 
    501 } 
    502 """ % {'body':model_info['form_volume']} 
    503         source.append(fn) 
    504  
    505     # Fill in definitions for Iq parameters 
    506     defines.append(('IQ_KERNEL_NAME', model_info['name'] + '_Iq')) 
    507     defines.append(('IQ_PARAMETERS', ', '.join(iq_parameters))) 
    508     if fixed_1d: 
    509         defines.append(('IQ_FIXED_PARAMETER_DECLARATIONS', 
    510                         ', \\\n    '.join('const double %s' % p for p in fixed_1d))) 
    511     if pd_1d: 
    512         defines.append(('IQ_WEIGHT_PRODUCT', 
    513                         '*'.join(p + '_w' for p in pd_1d))) 
    514         defines.append(('IQ_DISPERSION_LENGTH_DECLARATIONS', 
    515                         ', \\\n    '.join('const int N%s' % p for p in pd_1d))) 
    516         defines.append(('IQ_DISPERSION_LENGTH_SUM', 
    517                         '+'.join('N' + p for p in pd_1d))) 
    518         open_loops, close_loops = build_polydispersity_loops(pd_1d) 
    519         defines.append(('IQ_OPEN_LOOPS', 
    520                         open_loops.replace('\n', ' \\\n'))) 
    521         defines.append(('IQ_CLOSE_LOOPS', 
    522                         close_loops.replace('\n', ' \\\n'))) 
    523     if model_info['Iq'] is not None: 
    524         defines.append(('IQ_PARAMETER_DECLARATIONS', 
    525                         ', '.join('double ' + p for p in iq_parameters))) 
    526         fn = """\ 
    527 double Iq(double q, IQ_PARAMETER_DECLARATIONS); 
    528 double Iq(double q, IQ_PARAMETER_DECLARATIONS) { 
    529     %(body)s 
    530 } 
    531 """ % {'body':model_info['Iq']} 
    532         source.append(fn) 
    533  
    534     # Fill in definitions for Iqxy parameters 
    535     defines.append(('IQXY_KERNEL_NAME', model_info['name'] + '_Iqxy')) 
    536     defines.append(('IQXY_PARAMETERS', ', '.join(iqxy_parameters))) 
    537     if fixed_2d: 
    538         defines.append(('IQXY_FIXED_PARAMETER_DECLARATIONS', 
    539                         ', \\\n    '.join('const double %s' % p for p in fixed_2d))) 
    540     if pd_2d: 
    541         defines.append(('IQXY_WEIGHT_PRODUCT', 
    542                         '*'.join(p + '_w' for p in pd_2d))) 
    543         defines.append(('IQXY_DISPERSION_LENGTH_DECLARATIONS', 
    544                         ', \\\n    '.join('const int N%s' % p for p in pd_2d))) 
    545         defines.append(('IQXY_DISPERSION_LENGTH_SUM', 
    546                         '+'.join('N' + p for p in pd_2d))) 
    547         open_loops, close_loops = build_polydispersity_loops(pd_2d) 
    548         defines.append(('IQXY_OPEN_LOOPS', 
    549                         open_loops.replace('\n', ' \\\n'))) 
    550         defines.append(('IQXY_CLOSE_LOOPS', 
    551                         close_loops.replace('\n', ' \\\n'))) 
    552     if model_info['Iqxy'] is not None: 
    553         defines.append(('IQXY_PARAMETER_DECLARATIONS', 
    554                         ', '.join('double ' + p for p in iqxy_parameters))) 
    555         fn = """\ 
    556 double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS); 
    557 double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS) { 
    558     %(body)s 
    559 } 
    560 """ % {'body':model_info['Iqxy']} 
    561         source.append(fn) 
    562  
    563     # Need to know if we have a theta parameter for Iqxy; it is not there 
    564     # for the magnetic sphere model, for example, which has a magnetic 
    565     # orientation but no shape orientation. 
    566     if 'theta' in pd_2d: 
    567         defines.append(('IQXY_HAS_THETA', '1')) 
    568  
    569     #for d in defines: print(d) 
    570     defines = '\n'.join('#define %s %s' % (k, v) for k, v in defines) 
    571     sources = '\n\n'.join(source) 
    572     return C_KERNEL_TEMPLATE % { 
    573         'DEFINES': defines, 
    574         'SOURCES': sources, 
    575         } 
    576  
    577 def categorize_parameters(pars): 
    578     """ 
    579     Build parameter categories out of the the parameter definitions. 
    580  
    581     Returns a dictionary of categories. 
    582  
    583     Note: these categories are subject to change, depending on the needs of 
    584     the UI and the needs of the kernel calling function. 
    585  
    586     The categories are as follows: 
    587  
    588     * *volume* list of volume parameter names 
    589     * *orientation* list of orientation parameters 
    590     * *magnetic* list of magnetic parameters 
    591     * *<empty string>* list of parameters that have no type info 
    592  
    593     Each parameter is in one and only one category. 
    594  
    595     The following derived categories are created: 
    596  
    597     * *fixed-1d* list of non-polydisperse parameters for 1D models 
    598     * *pd-1d* list of polydisperse parameters for 1D models 
    599     * *fixed-2d* list of non-polydisperse parameters for 2D models 
    600     * *pd-d2* list of polydisperse parameters for 2D models 
    601     """ 
    602     partype = { 
    603         'volume': [], 'orientation': [], 'magnetic': [], 'sld': [], '': [], 
    604         'fixed-1d': [], 'fixed-2d': [], 'pd-1d': [], 'pd-2d': [], 
    605         'pd-rel': set(), 
    606     } 
    607  
    608     for p in pars: 
    609         if p.type == 'volume': 
    610             partype['pd-1d'].append(p.name) 
    611             partype['pd-2d'].append(p.name) 
    612             partype['pd-rel'].add(p.name) 
    613         elif p.type == 'magnetic': 
    614             partype['fixed-2d'].append(p.name) 
    615         elif p.type == 'orientation': 
    616             partype['pd-2d'].append(p.name) 
    617         elif p.type in ('', 'sld'): 
    618             partype['fixed-1d'].append(p.name) 
    619             partype['fixed-2d'].append(p.name) 
    620         else: 
    621             raise ValueError("unknown parameter type %r" % p.type) 
    622         partype[p.type].append(p.name) 
    623  
    624     return partype 
    625  
    626 def process_parameters(model_info): 
    627     """ 
    628     Process parameter block, precalculating parameter details. 
    629     """ 
    630     # convert parameters into named tuples 
    631     for p in model_info['parameters']: 
    632         if p[4] == '' and (p[0].startswith('sld') or p[0].endswith('sld')): 
    633             p[4] = 'sld' 
    634             # TODO: make sure all models explicitly label their sld parameters 
    635             #raise ValueError("%s.%s needs to be explicitly set to type 'sld'" %(model_info['id'], p[0])) 
    636  
    637     pars = [Parameter(*p) for p in model_info['parameters']] 
    638     # Fill in the derived attributes 
    639     model_info['parameters'] = pars 
    640     partype = categorize_parameters(pars) 
    641     model_info['limits'] = dict((p.name, p.limits) for p in pars) 
    642     model_info['partype'] = partype 
    643     model_info['defaults'] = dict((p.name, p.default) for p in pars) 
    644     if model_info.get('demo', None) is None: 
    645         model_info['demo'] = model_info['defaults'] 
    646     model_info['has_2d'] = partype['orientation'] or partype['magnetic'] 
    647  
     455    partable = model_info.parameters 
     456 
     457    # Identify parameters for Iq, Iqxy, Iq_magnetic and form_volume. 
     458    # Note that scale and volume are not possible types. 
     459 
     460    # Load templates and user code 
     461    kernel_header = load_template('kernel_header.c') 
     462    kernel_code = load_template('kernel_iq.c') 
     463    user_code = [open(f).read() for f in model_sources(model_info)] 
     464 
     465    # Build initial sources 
     466    source = [kernel_header] + user_code 
     467 
     468    # Make parameters for q, qx, qy so that we can use them in declarations 
     469    q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 
     470    # Generate form_volume function, etc. from body only 
     471    if isinstance(model_info.form_volume, str): 
     472        pars = partable.form_volume_parameters 
     473        source.append(_gen_fn('form_volume', pars, model_info.form_volume)) 
     474    if isinstance(model_info.Iq, str): 
     475        pars = [q] + partable.iq_parameters 
     476        source.append(_gen_fn('Iq', pars, model_info.Iq)) 
     477    if isinstance(model_info.Iqxy, str): 
     478        pars = [qx, qy] + partable.iqxy_parameters 
     479        source.append(_gen_fn('Iqxy', pars, model_info.Iqxy)) 
     480 
     481    # Define the parameter table 
     482    source.append("#define PARAMETER_TABLE \\") 
     483    source.append("\\\n".join(p.as_definition() 
     484                              for p in partable.kernel_parameters)) 
     485 
     486    # Define the function calls 
     487    if partable.form_volume_parameters: 
     488        refs = _call_pars("_v.", partable.form_volume_parameters) 
     489        call_volume = "#define CALL_VOLUME(_v) form_volume(%s)" % (",".join(refs)) 
     490    else: 
     491        # Model doesn't have volume.  We could make the kernel run a little 
     492        # faster by not using/transferring the volume normalizations, but 
     493        # the ifdef's reduce readability more than is worthwhile. 
     494        call_volume = "#define CALL_VOLUME(v) 1.0" 
     495    source.append(call_volume) 
     496 
     497    refs = ["_q[_i]"] + _call_pars("_v.", partable.iq_parameters) 
     498    call_iq = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(refs)) 
     499    if _have_Iqxy(user_code): 
     500        # Call 2D model 
     501        refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("_v.", partable.iqxy_parameters) 
     502        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iqxy(%s)" % (",".join(refs)) 
     503    else: 
     504        # Call 1D model with sqrt(qx^2 + qy^2) 
     505        warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 
     506        # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 
     507        pars_sqrt = ["sqrt(_q[2*_i]*_q[2*_i]+_q[2*_i+1]*_q[2*_i+1])"] + refs[1:] 
     508        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt)) 
     509 
     510    # Fill in definitions for numbers of parameters 
     511    source.append("#define MAX_PD %s"%partable.max_pd) 
     512    source.append("#define NPARS %d"%partable.npars) 
     513 
     514    # TODO: allow mixed python/opencl kernels? 
     515 
     516    # define the Iq kernel 
     517    source.append("#define KERNEL_NAME %s_Iq"%model_info.name) 
     518    source.append(call_iq) 
     519    source.append(kernel_code) 
     520    source.append("#undef CALL_IQ") 
     521    source.append("#undef KERNEL_NAME") 
     522 
     523    # define the Iqxy kernel from the same source with different #defines 
     524    source.append("#define KERNEL_NAME %s_Iqxy"%model_info.name) 
     525    source.append(call_iqxy) 
     526    source.append(kernel_code) 
     527    source.append("#undef CALL_IQ") 
     528    source.append("#undef KERNEL_NAME") 
     529 
     530    return '\n'.join(source) 
    648531 
    649532def load_kernel_module(model_name): 
     533    # type: (str) -> module 
    650534    if model_name.endswith('.py'): 
    651535        kernel_module = load_custom_kernel_module(model_name) 
    … …  
    657541 
    658542 
    659 def make_model_info(kernel_module): 
    660     """ 
    661     Interpret the model definition file, categorizing the parameters. 
    662  
    663     The module can be loaded with a normal python import statement if you 
    664     know which module you need, or with __import__('sasmodels.model.'+name) 
    665     if the name is in a string. 
    666  
    667     The *model_info* structure contains the following fields: 
    668  
    669     * *id* is the id of the kernel 
    670     * *name* is the display name of the kernel 
    671     * *filename* is the full path to the module defining the file (if any) 
    672     * *title* is a short description of the kernel 
    673     * *description* is a long description of the kernel (this doesn't seem 
    674       very useful since the Help button on the model page brings you directly 
    675       to the documentation page) 
    676     * *docs* is the docstring from the module.  Use :func:`make_doc` to 
    677     * *category* specifies the model location in the docs 
    678     * *parameters* is the model parameter table 
    679     * *single* is True if the model allows single precision 
    680     * *structure_factor* is True if the model is useable in a product 
    681     * *defaults* is the *{parameter: value}* table built from the parameter 
    682       description table. 
    683     * *limits* is the *{parameter: [min, max]}* table built from the 
    684       parameter description table. 
    685     * *partypes* categorizes the model parameters. See 
    686       :func:`categorize_parameters` for details. 
    687     * *demo* contains the *{parameter: value}* map used in compare (and maybe 
    688       for the demo plot, if plots aren't set up to use the default values). 
    689       If *demo* is not given in the file, then the default values will be used. 
    690     * *tests* is a set of tests that must pass 
    691     * *source* is the list of library files to include in the C model build 
    692     * *Iq*, *Iqxy*, *form_volume*, *ER*, *VR* and *sesans* are python functions 
    693       implementing the kernel for the module, or None if they are not 
    694       defined in python 
    695     * *composition* is None if the model is independent, otherwise it is a 
    696       tuple with composition type ('product' or 'mixture') and a list of 
    697       *model_info* blocks for the composition objects.  This allows us to 
    698       build complete product and mixture models from just the info. 
    699  
    700     """ 
    701     parameters = COMMON_PARAMETERS + kernel_module.parameters 
    702     filename = abspath(kernel_module.__file__) 
    703     kernel_id = splitext(basename(filename))[0] 
    704     name = getattr(kernel_module, 'name', None) 
    705     if name is None: 
    706         name = " ".join(w.capitalize() for w in kernel_id.split('_')) 
    707     model_info = dict( 
    708         id=kernel_id,  # string used to load the kernel 
    709         filename=abspath(kernel_module.__file__), 
    710         name=name, 
    711         title=kernel_module.title, 
    712         description=kernel_module.description, 
    713         parameters=parameters, 
    714         composition=None, 
    715         docs=kernel_module.__doc__, 
    716         category=getattr(kernel_module, 'category', None), 
    717         single=getattr(kernel_module, 'single', True), 
    718         structure_factor=getattr(kernel_module, 'structure_factor', False), 
    719         control=getattr(kernel_module, 'control', None), 
    720         demo=getattr(kernel_module, 'demo', None), 
    721         source=getattr(kernel_module, 'source', []), 
    722         tests=getattr(kernel_module, 'tests', []), 
    723         ) 
    724     process_parameters(model_info) 
    725     # Check for optional functions 
    726     functions = "ER VR form_volume Iq Iqxy shape sesans".split() 
    727     model_info.update((k, getattr(kernel_module, k, None)) for k in functions) 
    728     return model_info 
    729543 
    730544section_marker = re.compile(r'\A(?P<first>[%s])(?P=first)*\Z' 
    731545                            %re.escape(string.punctuation)) 
    732546def _convert_section_titles_to_boldface(lines): 
     547    # type: (Sequence[str]) -> Iterator[str] 
    733548    """ 
    734549    Do the actual work of identifying and converting section headings. 
    … …  
    752567 
    753568def convert_section_titles_to_boldface(s): 
     569    # type: (str) -> str 
    754570    """ 
    755571    Use explicit bold-face rather than section headings so that the table of 
    … …  
    762578 
    763579def make_doc(model_info): 
     580    # type: (ModelInfo) -> str 
    764581    """ 
    765582    Return the documentation for the model. 
    … …  
    767584    Iq_units = "The returned value is scaled to units of |cm^-1| |sr^-1|, absolute scale." 
    768585    Sq_units = "The returned value is a dimensionless structure factor, $S(q)$." 
    769     docs = convert_section_titles_to_boldface(model_info['docs']) 
    770     subst = dict(id=model_info['id'].replace('_', '-'), 
    771                  name=model_info['name'], 
    772                  title=model_info['title'], 
    773                  parameters=make_partable(model_info['parameters']), 
    774                  returns=Sq_units if model_info['structure_factor'] else Iq_units, 
     586    docs = convert_section_titles_to_boldface(model_info.docs) 
     587    subst = dict(id=model_info.id.replace('_', '-'), 
     588                 name=model_info.name, 
     589                 title=model_info.title, 
     590                 parameters=make_partable(model_info.parameters.kernel_parameters), 
     591                 returns=Sq_units if model_info.structure_factor else Iq_units, 
    775592                 docs=docs) 
    776593    return DOC_HEADER % subst 
    … …  
    778595 
    779596def demo_time(): 
     597    # type: () -> None 
    780598    """ 
    781599    Show how long it takes to process a model. 
    782600    """ 
     601    import datetime 
     602    from .modelinfo import make_model_info 
    783603    from .models import cylinder 
    784     import datetime 
     604 
    785605    tic = datetime.datetime.now() 
    786606    make_source(make_model_info(cylinder)) 
    … …  
    789609 
    790610def main(): 
     611    # type: () -> None 
    791612    """ 
    792613    Program which prints the source produced by the model. 
    793614    """ 
     615    import sys 
     616    from .modelinfo import make_model_info 
     617 
    794618    if len(sys.argv) <= 1: 
    795619        print("usage: python -m sasmodels.generate modelname") 

  • sasmodels/kernelcl.py

    r4d76711 r7ae2b7f  
    4848harmless, albeit annoying. 
    4949""" 
     50from __future__ import print_function 
    5051import os 
    5152import warnings 
    5253 
    53 import numpy as np 
     54import numpy as np  # type: ignore 
    5455 
    5556try: 
    56     import pyopencl as cl 
     57    raise NotImplementedError("OpenCL not yet implemented for new kernel template") 
     58    import pyopencl as cl  # type: ignore 
    5759    # Ask OpenCL for the default context so that we know that one exists 
    5860    cl.create_some_context(interactive=False) 
    … …  
    6163    raise RuntimeError("OpenCL not available") 
    6264 
    63 from pyopencl import mem_flags as mf 
    64 from pyopencl.characterize import get_fast_inaccurate_build_options 
     65from pyopencl import mem_flags as mf  # type: ignore 
     66from pyopencl.characterize import get_fast_inaccurate_build_options  # type: ignore 
    6567 
    6668from . import generate 
     69from .kernel import KernelModel, Kernel 
    6770 
    6871# The max loops number is limited by the amount of local memory available 
    … …  
    7376# of polydisperse parameters. 
    7477MAX_LOOPS = 2048 
     78 
     79 
     80# Pragmas for enable OpenCL features.  Be sure to protect them so that they 
     81# still compile even if OpenCL is not present. 
     82_F16_PRAGMA = """\ 
     83#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 
     84#  pragma OPENCL EXTENSION cl_khr_fp16: enable 
     85#endif 
     86""" 
     87 
     88_F64_PRAGMA = """\ 
     89#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 
     90#  pragma OPENCL EXTENSION cl_khr_fp64: enable 
     91#endif 
     92""" 
    7593 
    7694 
    … …  
    142160        raise RuntimeError("%s not supported for devices"%dtype) 
    143161 
    144     source = generate.convert_type(source, dtype) 
     162    source_list = [generate.convert_type(source, dtype)] 
     163 
     164    if dtype == generate.F16: 
     165        source_list.insert(0, _F16_PRAGMA) 
     166    elif dtype == generate.F64: 
     167        source_list.insert(0, _F64_PRAGMA) 
     168 
    145169    # Note: USE_SINCOS makes the intel cpu slower under opencl 
    146170    if context.devices[0].type == cl.device_type.GPU: 
    147         source = "#define USE_SINCOS\n" + source 
     171        source_list.insert(0, "#define USE_SINCOS\n") 
    148172    options = (get_fast_inaccurate_build_options(context.devices[0]) 
    149173               if fast else []) 
     174    source = "\n".join(source_list) 
    150175    program = cl.Program(context, source).build(options=options) 
    151176    return program 
    … …  
    220245        key = "%s-%s-%s"%(name, dtype, fast) 
    221246        if key not in self.compiled: 
    222             #print("compiling",name) 
     247            print("compiling",name) 
    223248            dtype = np.dtype(dtype) 
    224             program = compile_model(self.get_context(dtype), source, dtype, fast) 
     249            program = compile_model(self.get_context(dtype), 
     250                                    str(source), dtype, fast) 
    225251            self.compiled[key] = program 
    226252        return self.compiled[key] 
    … …  
    285311 
    286312 
    287 class GpuModel(object): 
     313class GpuModel(KernelModel): 
    288314    """ 
    289315    GPU wrapper for a single model. 
    … …  
    317343        if self.program is None: 
    318344            compiler = environment().compile_program 
    319             self.program = compiler(self.info['name'], self.source, self.dtype, 
    320                                     self.fast) 
     345            self.program = compiler(self.info.name, self.source, 
     346                                    self.dtype, self.fast) 
    321347        is_2d = len(q_vectors) == 2 
    322348        kernel_name = generate.kernel_name(self.info, is_2d) 
    … …  
    329355        """ 
    330356        if self.program is not None: 
    331             environment().release_program(self.info['name']) 
     357            environment().release_program(self.info.name) 
    332358            self.program = None 
    333359 
    … …  
    365391        # at this point, so instead using 32, which is good on the set of 
    366392        # architectures tested so far. 
    367         self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 
     393        if self.is_2d: 
     394            # Note: 17 rather than 15 because results is 2 elements 
     395            # longer than input. 
     396            width = ((self.nq+17)//16)*16 
     397            self.q = np.empty((width, 2), dtype=dtype) 
     398            self.q[:self.nq, 0] = q_vectors[0] 
     399            self.q[:self.nq, 1] = q_vectors[1] 
     400        else: 
     401            # Note: 33 rather than 31 because results is 2 elements 
     402            # longer than input. 
     403            width = ((self.nq+33)//32)*32 
     404            self.q = np.empty(width, dtype=dtype) 
     405            self.q[:self.nq] = q_vectors[0] 
     406        self.global_size = [self.q.shape[0]] 
    368407        context = env.get_context(self.dtype) 
    369         self.global_size = [self.q_vectors[0].size] 
    370408        #print("creating inputs of size", self.global_size) 
    371         self.q_buffers = [ 
    372             cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q) 
    373             for q in self.q_vectors 
    374         ] 
     409        self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     410                             hostbuf=self.q) 
    375411 
    376412    def release(self): 
    … …  
    378414        Free the memory. 
    379415        """ 
    380         for b in self.q_buffers: 
    381             b.release() 
    382         self.q_buffers = [] 
     416        if self.q is not None: 
     417            self.q.release() 
     418            self.q = None 
    383419 
    384420    def __del__(self): 
    385421        self.release() 
    386422 
    387 class GpuKernel(object): 
     423class GpuKernel(Kernel): 
    388424    """ 
    389425    Callable SAS kernel. 
    … …  
    406442    """ 
    407443    def __init__(self, kernel, model_info, q_vectors, dtype): 
     444        max_pd = model_info.max_pd 
     445        npars = len(model_info.parameters)-2 
    408446        q_input = GpuInput(q_vectors, dtype) 
     447        self.dtype = dtype 
     448        self.dim = '2d' if q_input.is_2d else '1d' 
    409449        self.kernel = kernel 
    410450        self.info = model_info 
    411         self.res = np.empty(q_input.nq, q_input.dtype) 
    412         dim = '2d' if q_input.is_2d else '1d' 
    413         self.fixed_pars = model_info['partype']['fixed-' + dim] 
    414         self.pd_pars = model_info['partype']['pd-' + dim] 
     451        self.pd_stop_index = 4*max_pd-1 
     452        # plus three for the normalization values 
     453        self.result = np.empty(q_input.nq+3, q_input.dtype) 
    415454 
    416455        # Inputs and outputs for each kernel call 
    … …  
    418457        env = environment() 
    419458        self.queue = env.get_queue(dtype) 
    420         self.loops_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
    421                                  2 * MAX_LOOPS * q_input.dtype.itemsize) 
    422         self.res_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
     459 
     460        # details is int32 data, padded to an 8 integer boundary 
     461        size = ((max_pd*5 + npars*3 + 2 + 7)//8)*8 
     462        self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
    423463                               q_input.global_size[0] * q_input.dtype.itemsize) 
    424         self.q_input = q_input 
    425  
    426         self._need_release = [self.loops_b, self.res_b, self.q_input] 
    427  
    428     def __call__(self, fixed_pars, pd_pars, cutoff): 
     464        self.q_input = q_input # allocated by GpuInput above 
     465 
     466        self._need_release = [ self.result_b, self.q_input ] 
     467 
     468    def __call__(self, call_details, weights, values, cutoff): 
    429469        real = (np.float32 if self.q_input.dtype == generate.F32 
    430470                else np.float64 if self.q_input.dtype == generate.F64 
    431471                else np.float16 if self.q_input.dtype == generate.F16 
    432472                else np.float32)  # will never get here, so use np.float32 
    433  
    434         #print "pars", fixed_pars, pd_pars 
    435         res_bi = self.res_b 
    436         nq = np.uint32(self.q_input.nq) 
    437         if pd_pars: 
    438             cutoff = real(cutoff) 
    439             loops_N = [np.uint32(len(p[0])) for p in pd_pars] 
    440             loops = np.hstack(pd_pars) \ 
    441                 if pd_pars else np.empty(0, dtype=self.q_input.dtype) 
    442             loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten() 
    443             #print("loops",Nloops, loops) 
    444  
    445             #import sys; print("opencl eval",pars) 
    446             #print("opencl eval",pars) 
    447             if len(loops) > 2 * MAX_LOOPS: 
    448                 raise ValueError("too many polydispersity points") 
    449  
    450             loops_bi = self.loops_b 
    451             cl.enqueue_copy(self.queue, loops_bi, loops) 
    452             loops_l = cl.LocalMemory(len(loops.data)) 
    453             #ctx = environment().context 
    454             #loops_bi = cl.Buffer(ctx, mf.READ_ONLY|mf.COPY_HOST_PTR, hostbuf=loops) 
    455             dispersed = [loops_bi, loops_l, cutoff] + loops_N 
    456         else: 
    457             dispersed = [] 
    458         fixed = [real(p) for p in fixed_pars] 
    459         args = self.q_input.q_buffers + [res_bi, nq] + dispersed + fixed 
     473        assert call_details.dtype == np.int32 
     474        assert weights.dtype == real and values.dtype == real 
     475 
     476        context = self.queue.context 
     477        details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     478                              hostbuf=call_details) 
     479        weights_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     480                              hostbuf=weights) 
     481        values_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     482                             hostbuf=values) 
     483 
     484        start, stop = 0, self.details[self.pd_stop_index] 
     485        args = [ 
     486            np.uint32(self.q_input.nq), np.uint32(start), np.uint32(stop), 
     487            self.details_b, self.weights_b, self.values_b, 
     488            self.q_input.q_b, self.result_b, real(cutoff), 
     489        ] 
    460490        self.kernel(self.queue, self.q_input.global_size, None, *args) 
    461         cl.enqueue_copy(self.queue, self.res, res_bi) 
    462  
    463         return self.res 
     491        cl.enqueue_copy(self.queue, self.result, self.result_b) 
     492        [v.release() for v in (details_b, weights_b, values_b)] 
     493 
     494        return self.result[:self.nq] 
    464495 
    465496    def release(self): 

  • sasmodels/kerneldll.py

    r4d76711 r7ae2b7f  
    4949import os 
    5050import tempfile 
    51 import ctypes as ct 
    52 from ctypes import c_void_p, c_int, c_longdouble, c_double, c_float 
    53 import _ctypes 
    54  
    55 import numpy as np 
     51import ctypes as ct  # type: ignore 
     52from ctypes import c_void_p, c_int32, c_longdouble, c_double, c_float  # type: ignore 
     53 
     54import numpy as np  # type: ignore 
    5655 
    5756from . import generate 
    58 from .kernelpy import PyInput, PyModel 
     57from . import details 
     58from .kernel import KernelModel, Kernel 
     59from .kernelpy import PyInput 
    5960from .exception import annotate_exception 
     61from .generate import F16, F32, F64 
     62 
     63try: 
     64    from typing import Tuple, Callable, Any 
     65    from .modelinfo import ModelInfo 
     66    from .details import CallDetails 
     67except ImportError: 
     68    pass 
    6069 
    6170# Compiler platform details 
    … …  
    8190        COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 
    8291        if "SAS_OPENMP" in os.environ: 
    83             COMPILE = COMPILE + " -fopenmp" 
     92            COMPILE += " -fopenmp" 
    8493else: 
    8594    COMPILE = "cc -shared -fPIC -fopenmp -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 
    … …  
    9099 
    91100 
    92 def dll_path(model_info, dtype="double"): 
    93     """ 
    94     Path to the compiled model defined by *model_info*. 
    95     """ 
    96     from os.path import join as joinpath, split as splitpath, splitext 
    97     basename = splitext(splitpath(model_info['filename'])[1])[0] 
    98     if np.dtype(dtype) == generate.F32: 
    99         basename += "32" 
    100     elif np.dtype(dtype) == generate.F64: 
    101         basename += "64" 
    102     else: 
    103         basename += "128" 
    104     return joinpath(DLL_PATH, basename+'.so') 
    105  
    106  
    107 def make_dll(source, model_info, dtype="double"): 
    108     """ 
    109     Load the compiled model defined by *kernel_module*. 
    110  
    111     Recompile if any files are newer than the model file. 
     101def dll_name(model_info, dtype): 
     102    # type: (ModelInfo, np.dtype) ->  str 
     103    """ 
     104    Name of the dll containing the model.  This is the base file name without 
     105    any path or extension, with a form such as 'sas_sphere32'. 
     106    """ 
     107    bits = 8*dtype.itemsize 
     108    return "sas_%s%d"%(model_info.id, bits) 
     109 
     110 
     111def dll_path(model_info, dtype): 
     112    # type: (ModelInfo, np.dtype) -> str 
     113    """ 
     114    Complete path to the dll for the model.  Note that the dll may not 
     115    exist yet if it hasn't been compiled. 
     116    """ 
     117    return os.path.join(DLL_PATH, dll_name(model_info, dtype)+".so") 
     118 
     119 
     120def make_dll(source, model_info, dtype=F64): 
     121    # type: (str, ModelInfo, np.dtype) -> str 
     122    """ 
     123    Returns the path to the compiled model defined by *kernel_module*. 
     124 
     125    If the model has not been compiled, or if the source file(s) are newer 
     126    than the dll, then *make_dll* will compile the model before returning. 
     127    This routine does not load the resulting dll. 
    112128 
    113129    *dtype* is a numpy floating point precision specifier indicating whether 
    114     the model should be single or double precision.  The default is double 
    115     precision. 
    116  
    117     The DLL is not loaded until the kernel is called so models can 
    118     be defined without using too many resources. 
     130    the model should be single, double or long double precision.  The default 
     131    is double precision, *np.dtype('d')*. 
     132 
     133    Set *sasmodels.ALLOW_SINGLE_PRECISION_DLLS* to False if single precision 
     134    models are not allowed as DLLs. 
    119135 
    120136    Set *sasmodels.kerneldll.DLL_PATH* to the compiled dll output path. 
    121137    The default is the system temporary directory. 
    122  
    123     Set *sasmodels.ALLOW_SINGLE_PRECISION_DLLS* to True if single precision 
    124     models are allowed as DLLs. 
    125     """ 
    126     if callable(model_info.get('Iq', None)): 
    127         return PyModel(model_info) 
    128      
    129     dtype = np.dtype(dtype) 
    130     if dtype == generate.F16: 
     138    """ 
     139    if dtype == F16: 
    131140        raise ValueError("16 bit floats not supported") 
    132     if dtype == generate.F32 and not ALLOW_SINGLE_PRECISION_DLLS: 
    133         dtype = generate.F64  # Force 64-bit dll 
    134  
    135     if dtype == generate.F32: # 32-bit dll 
    136         tempfile_prefix = 'sas_' + model_info['name'] + '32_' 
    137     elif dtype == generate.F64: 
    138         tempfile_prefix = 'sas_' + model_info['name'] + '64_' 
    139     else: 
    140         tempfile_prefix = 'sas_' + model_info['name'] + '128_' 
    141   
    142     source = generate.convert_type(source, dtype) 
    143     source_files = generate.model_sources(model_info) + [model_info['filename']] 
     141    if dtype == F32 and not ALLOW_SINGLE_PRECISION_DLLS: 
     142        dtype = F64  # Force 64-bit dll 
     143    # Note: dtype may be F128 for long double precision 
     144 
     145    newest = generate.timestamp(model_info) 
    144146    dll = dll_path(model_info, dtype) 
    145     newest = max(os.path.getmtime(f) for f in source_files) 
    146147    if not os.path.exists(dll) or os.path.getmtime(dll) < newest: 
    147         # Replace with a proper temp file 
    148         fid, filename = tempfile.mkstemp(suffix=".c", prefix=tempfile_prefix) 
     148        basename = dll_name(model_info, dtype) + "_" 
     149        fid, filename = tempfile.mkstemp(suffix=".c", prefix=basename) 
     150        source = generate.convert_type(source, dtype) 
    149151        os.fdopen(fid, "w").write(source) 
    150152        command = COMPILE%{"source":filename, "output":dll} 
    … …  
    160162 
    161163 
    162 def load_dll(source, model_info, dtype="double"): 
     164def load_dll(source, model_info, dtype=F64): 
     165    # type: (str, ModelInfo, np.dtype) -> "DllModel" 
    163166    """ 
    164167    Create and load a dll corresponding to the source, info pair returned 
    … …  
    172175 
    173176 
    174 IQ_ARGS = [c_void_p, c_void_p, c_int] 
    175 IQXY_ARGS = [c_void_p, c_void_p, c_void_p, c_int] 
    176  
    177 class DllModel(object): 
     177class DllModel(KernelModel): 
    178178    """ 
    179179    ctypes wrapper for a single model. 
    … …  
    191191     
    192192    def __init__(self, dllpath, model_info, dtype=generate.F32): 
     193        # type: (str, ModelInfo, np.dtype) -> None 
    193194        self.info = model_info 
    194195        self.dllpath = dllpath 
    195         self.dll = None 
     196        self._dll = None  # type: ct.CDLL 
    196197        self.dtype = np.dtype(dtype) 
    197198 
    198199    def _load_dll(self): 
    199         Nfixed1d = len(self.info['partype']['fixed-1d']) 
    200         Nfixed2d = len(self.info['partype']['fixed-2d']) 
    201         Npd1d = len(self.info['partype']['pd-1d']) 
    202         Npd2d = len(self.info['partype']['pd-2d']) 
    203  
     200        # type: () -> None 
    204201        #print("dll", self.dllpath) 
    205202        try: 
    206             self.dll = ct.CDLL(self.dllpath) 
     203            self._dll = ct.CDLL(self.dllpath) 
    207204        except: 
    208205            annotate_exception("while loading "+self.dllpath) 
    … …  
    212209              else c_double if self.dtype == generate.F64 
    213210              else c_longdouble) 
    214         pd_args_1d = [c_void_p, fp] + [c_int]*Npd1d if Npd1d else [] 
    215         pd_args_2d = [c_void_p, fp] + [c_int]*Npd2d if Npd2d else [] 
    216         self.Iq = self.dll[generate.kernel_name(self.info, False)] 
    217         self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d 
    218  
    219         self.Iqxy = self.dll[generate.kernel_name(self.info, True)] 
    220         self.Iqxy.argtypes = IQXY_ARGS + pd_args_2d + [fp]*Nfixed2d 
    221          
    222         self.release() 
     211 
     212        # int, int, int, int*, double*, double*, double*, double*, double*, double 
     213        argtypes = [c_int32]*3 + [c_void_p]*5 + [fp] 
     214        self._Iq = self._dll[generate.kernel_name(self.info, is_2d=False)] 
     215        self._Iqxy = self._dll[generate.kernel_name(self.info, is_2d=True)] 
     216        self._Iq.argtypes = argtypes 
     217        self._Iqxy.argtypes = argtypes 
    223218 
    224219    def __getstate__(self): 
     220        # type: () -> Tuple[ModelInfo, str] 
    225221        return self.info, self.dllpath 
    226222 
    227223    def __setstate__(self, state): 
     224        # type: (Tuple[ModelInfo, str]) -> None 
    228225        self.info, self.dllpath = state 
    229         self.dll = None 
     226        self._dll = None 
    230227 
    231228    def make_kernel(self, q_vectors): 
     229        # type: (List[np.ndarray]) -> DllKernel 
    232230        q_input = PyInput(q_vectors, self.dtype) 
    233         if self.dll is None: self._load_dll() 
    234         kernel = self.Iqxy if q_input.is_2d else self.Iq 
     231        # Note: pickle not supported for DllKernel 
     232        if self._dll is None: 
     233            self._load_dll() 
     234        kernel = self._Iqxy if q_input.is_2d else self._Iq 
    235235        return DllKernel(kernel, self.info, q_input) 
    236236 
    237237    def release(self): 
     238        # type: () -> None 
    238239        """ 
    239240        Release any resources associated with the model. 
    … …  
    244245            libHandle = dll._handle 
    245246            #libHandle = ct.c_void_p(dll._handle) 
    246             del dll, self.dll 
    247             self.dll = None 
     247            del dll, self._dll 
     248            self._dll = None 
    248249            #_ctypes.FreeLibrary(libHandle) 
    249250            ct.windll.kernel32.FreeLibrary(libHandle) 
    … …  
    252253 
    253254 
    254 class DllKernel(object): 
     255class DllKernel(Kernel): 
    255256    """ 
    256257    Callable SAS kernel. 
    … …  
    272273    """ 
    273274    def __init__(self, kernel, model_info, q_input): 
     275        # type: (Callable[[], np.ndarray], ModelInfo, PyInput) -> None 
     276        self.kernel = kernel 
    274277        self.info = model_info 
    275278        self.q_input = q_input 
    276         self.kernel = kernel 
    277         self.res = np.empty(q_input.nq, q_input.dtype) 
    278         dim = '2d' if q_input.is_2d else '1d' 
    279         self.fixed_pars = model_info['partype']['fixed-' + dim] 
    280         self.pd_pars = model_info['partype']['pd-' + dim] 
    281  
    282         # In dll kernel, but not in opencl kernel 
    283         self.p_res = self.res.ctypes.data 
    284  
    285     def __call__(self, fixed_pars, pd_pars, cutoff): 
     279        self.dtype = q_input.dtype 
     280        self.dim = '2d' if q_input.is_2d else '1d' 
     281        self.result = np.empty(q_input.nq+3, q_input.dtype) 
     282 
     283    def __call__(self, call_details, weights, values, cutoff): 
     284        # type: (CallDetails, np.ndarray, np.ndarray, float) -> np.ndarray 
    286285        real = (np.float32 if self.q_input.dtype == generate.F32 
    287286                else np.float64 if self.q_input.dtype == generate.F64 
    288287                else np.float128) 
    289  
    290         nq = c_int(self.q_input.nq) 
    291         if pd_pars: 
    292             cutoff = real(cutoff) 
    293             loops_N = [np.uint32(len(p[0])) for p in pd_pars] 
    294             loops = np.hstack(pd_pars) 
    295             loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten() 
    296             p_loops = loops.ctypes.data 
    297             dispersed = [p_loops, cutoff] + loops_N 
    298         else: 
    299             dispersed = [] 
    300         fixed = [real(p) for p in fixed_pars] 
    301         args = self.q_input.q_pointers + [self.p_res, nq] + dispersed + fixed 
    302         #print(pars) 
    303         self.kernel(*args) 
    304  
    305         return self.res 
     288        assert isinstance(call_details, details.CallDetails) 
     289        assert weights.dtype == real and values.dtype == real 
     290 
     291        start, stop = 0, call_details.total_pd 
     292        #print("in kerneldll") 
     293        #print("weights", weights) 
     294        #print("values", values) 
     295        args = [ 
     296            self.q_input.nq, # nq 
     297            start, # pd_start 
     298            stop, # pd_stop pd_stride[MAX_PD] 
     299            call_details.ctypes.data, # problem 
     300            weights.ctypes.data,  # weights 
     301            values.ctypes.data,  #pars 
     302            self.q_input.q.ctypes.data, #q 
     303            self.result.ctypes.data,   # results 
     304            real(cutoff), # cutoff 
     305            ] 
     306        self.kernel(*args) # type: ignore 
     307        return self.result[:-3] 
    306308 
    307309    def release(self): 
     310        # type: () -> None 
    308311        """ 
    309312        Release any resources associated with the kernel. 
    310313        """ 
    311         pass 
     314        self.q_input.release() 

  • sasmodels/kernelpy.py

    r4d76711 r7ae2b7f  
    77:class:`kernelcl.GpuModel` and :class:`kerneldll.DllModel`. 
    88""" 
    9 import numpy as np 
    10 from numpy import pi, cos 
    11  
     9import numpy as np  # type: ignore 
     10from numpy import pi, cos  #type: ignore 
     11 
     12from . import details 
    1213from .generate import F64 
    13  
    14 class PyModel(object): 
     14from .kernel import KernelModel, Kernel 
     15 
     16try: 
     17    from typing import Union, Callable 
     18except: 
     19    pass 
     20else: 
     21    DType = Union[None, str, np.dtype] 
     22 
     23class PyModel(KernelModel): 
    1524    """ 
    1625    Wrapper for pure python models. 
    1726    """ 
    1827    def __init__(self, model_info): 
     28        # Make sure Iq and Iqxy are available and vectorized 
     29        _create_default_functions(model_info) 
    1930        self.info = model_info 
    2031 
    2132    def make_kernel(self, q_vectors): 
    2233        q_input = PyInput(q_vectors, dtype=F64) 
    23         kernel = self.info['Iqxy'] if q_input.is_2d else self.info['Iq'] 
     34        kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq 
    2435        return PyKernel(kernel, self.info, q_input) 
    2536 
    … …  
    5364        self.dtype = dtype 
    5465        self.is_2d = (len(q_vectors) == 2) 
    55         self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors] 
    56         self.q_pointers = [q.ctypes.data for q in self.q_vectors] 
     66        if self.is_2d: 
     67            self.q = np.empty((self.nq, 2), dtype=dtype) 
     68            self.q[:, 0] = q_vectors[0] 
     69            self.q[:, 1] = q_vectors[1] 
     70        else: 
     71            self.q = np.empty(self.nq, dtype=dtype) 
     72            self.q[:self.nq] = q_vectors[0] 
    5773 
    5874    def release(self): 
    … …  
    6076        Free resources associated with the model inputs. 
    6177        """ 
    62         self.q_vectors = [] 
    63  
    64 class PyKernel(object): 
     78        self.q = None 
     79 
     80class PyKernel(Kernel): 
    6581    """ 
    6682    Callable SAS kernel. 
    … …  
    8298    """ 
    8399    def __init__(self, kernel, model_info, q_input): 
     100        self.dtype = np.dtype('d') 
    84101        self.info = model_info 
    85102        self.q_input = q_input 
    86103        self.res = np.empty(q_input.nq, q_input.dtype) 
    87         dim = '2d' if q_input.is_2d else '1d' 
    88         # Loop over q unless user promises that the kernel is vectorized by 
    89         # taggining it with vectorized=True 
    90         if not getattr(kernel, 'vectorized', False): 
    91             if dim == '2d': 
    92                 def vector_kernel(qx, qy, *args): 
    93                     """ 
    94                     Vectorized 2D kernel. 
    95                     """ 
    96                     return np.array([kernel(qxi, qyi, *args) 
    97                                      for qxi, qyi in zip(qx, qy)]) 
     104        self.kernel = kernel 
     105        self.dim = '2d' if q_input.is_2d else '1d' 
     106 
     107        partable = model_info.parameters 
     108        kernel_parameters = (partable.iqxy_parameters if q_input.is_2d 
     109                             else partable.iq_parameters) 
     110        volume_parameters = partable.form_volume_parameters 
     111 
     112        # Create an array to hold the parameter values.  There will be a 
     113        # single array whose values are updated as the calculator goes 
     114        # through the loop.  Arguments to the kernel and volume functions 
     115        # will use views into this vector, relying on the fact that a 
     116        # an array of no dimensions acts like a scalar. 
     117        parameter_vector = np.empty(len(partable.call_parameters)-2, 'd') 
     118 
     119        # Create views into the array to hold the arguments 
     120        offset = 0 
     121        kernel_args, volume_args = [], [] 
     122        for p in partable.kernel_parameters: 
     123            if p.length == 1: 
     124                # Scalar values are length 1 vectors with no dimensions. 
     125                v = parameter_vector[offset:offset+1].reshape(()) 
    98126            else: 
    99                 def vector_kernel(q, *args): 
    100                     """ 
    101                     Vectorized 1D kernel. 
    102                     """ 
    103                     return np.array([kernel(qi, *args) 
    104                                      for qi in q]) 
    105             self.kernel = vector_kernel 
     127                # Vector values are simple views. 
     128                v = parameter_vector[offset:offset+p.length] 
     129            offset += p.length 
     130            if p in kernel_parameters: 
     131                kernel_args.append(v) 
     132            if p in volume_parameters: 
     133                volume_args.append(v) 
     134 
     135        # Hold on to the parameter vector so we can use it to call kernel later. 
     136        # This may also be required to preserve the views into the vector. 
     137        self._parameter_vector = parameter_vector 
     138 
     139        # Generate a closure which calls the kernel with the views into the 
     140        # parameter array. 
     141        if q_input.is_2d: 
     142            form = model_info.Iqxy 
     143            qx, qy = q_input.q[:,0], q_input.q[:,1] 
     144            self._form = lambda: form(qx, qy, *kernel_args) 
    106145        else: 
    107             self.kernel = kernel 
    108         fixed_pars = model_info['partype']['fixed-' + dim] 
    109         pd_pars = model_info['partype']['pd-' + dim] 
    110         vol_pars = model_info['partype']['volume'] 
    111  
    112         # First two fixed pars are scale and background 
    113         pars = [p.name for p in model_info['parameters'][2:]] 
    114         offset = len(self.q_input.q_vectors) 
    115         self.args = self.q_input.q_vectors + [None] * len(pars) 
    116         self.fixed_index = np.array([pars.index(p) + offset 
    117                                      for p in fixed_pars[2:]]) 
    118         self.pd_index = np.array([pars.index(p) + offset 
    119                                   for p in pd_pars]) 
    120         self.vol_index = np.array([pars.index(p) + offset 
    121                                    for p in vol_pars]) 
    122         try: self.theta_index = pars.index('theta') + offset 
    123         except ValueError: self.theta_index = -1 
    124  
    125         # Caller needs fixed_pars and pd_pars 
    126         self.fixed_pars = fixed_pars 
    127         self.pd_pars = pd_pars 
    128  
    129     def __call__(self, fixed, pd, cutoff=1e-5): 
    130         #print("fixed",fixed) 
    131         #print("pd", pd) 
    132         args = self.args[:]  # grab a copy of the args 
    133         form, form_volume = self.kernel, self.info['form_volume'] 
    134         # First two fixed 
    135         scale, background = fixed[:2] 
    136         for index, value in zip(self.fixed_index, fixed[2:]): 
    137             args[index] = float(value) 
    138         res = _loops(form, form_volume, cutoff, scale, background, args, 
    139                      pd, self.pd_index, self.vol_index, self.theta_index) 
    140  
     146            form = model_info.Iq 
     147            q = q_input.q 
     148            self._form = lambda: form(q, *kernel_args) 
     149 
     150        # Generate a closure which calls the form_volume if it exists. 
     151        form_volume = model_info.form_volume 
     152        self._volume = ((lambda: form_volume(*volume_args)) if form_volume 
     153                        else (lambda: 1.0)) 
     154 
     155    def __call__(self, call_details, weights, values, cutoff): 
     156        assert isinstance(call_details, details.CallDetails) 
     157        res = _loops(self._parameter_vector, self._form, self._volume, 
     158                     self.q_input.nq, call_details, weights, values, cutoff) 
    141159        return res 
    142160 
    … …  
    145163        Free resources associated with the kernel. 
    146164        """ 
     165        self.q_input.release() 
    147166        self.q_input = None 
    148167 
    149 def _loops(form, form_volume, cutoff, scale, background, 
    150            args, pd, pd_index, vol_index, theta_index): 
    151     """ 
    152     *form* is the name of the form function, which should be vectorized over 
    153     q, but otherwise have an interface like the opencl kernels, with the 
    154     q parameters first followed by the individual parameters in the order 
    155     given in model.parameters (see :mod:`sasmodels.generate`). 
    156  
    157     *form_volume* calculates the volume of the shape.  *vol_index* gives 
    158     the list of volume parameters 
    159  
    160     *cutoff* ignores the corners of the dispersion hypercube 
    161  
    162     *scale*, *background* multiplies the resulting form and adds an offset 
    163  
    164     *args* is the prepopulated set of arguments to the form function, starting 
    165     with the q vectors, and including slots for all the parameters.  The 
    166     values for the parameters will be substituted with values from the 
    167     dispersion functions. 
    168  
    169     *pd* is the list of dispersion parameters 
    170  
    171     *pd_index* are the indices of the dispersion parameters in *args* 
    172  
    173     *vol_index* are the indices of the volume parameters in *args* 
    174  
    175     *theta_index* is the index of the theta parameter for the sasview 
    176     spherical correction, or -1 if there is no angular dispersion 
    177     """ 
    178  
     168def _loops(parameters, form, form_volume, nq, call_details, 
     169           weights, values, cutoff): 
     170    # type: (np.ndarray, Callable[[], np.ndarray], Callable[[], float], int, details.CallDetails, np.ndarray, np.ndarray, float) -> None 
    179171    ################################################################ 
    180172    #                                                              # 
    … …  
    186178    #                                                              # 
    187179    ################################################################ 
    188  
    189     weight = np.empty(len(pd), 'd') 
    190     if weight.size > 0: 
    191         # weight vector, to be populated by polydispersity loops 
    192  
    193         # identify which pd parameters are volume parameters 
    194         vol_weight_index = np.array([(index in vol_index) for index in pd_index]) 
    195  
    196         # Sort parameters in decreasing order of pd length 
    197         Npd = np.array([len(pdi[0]) for pdi in pd], 'i') 
    198         order = np.argsort(Npd)[::-1] 
    199         stride = np.cumprod(Npd[order]) 
    200         pd = [pd[index] for index in order] 
    201         pd_index = pd_index[order] 
    202         vol_weight_index = vol_weight_index[order] 
    203  
    204         fast_value = pd[0][0] 
    205         fast_weight = pd[0][1] 
     180    parameters[:] = values[call_details.par_offset] 
     181    scale, background = values[0], values[1] 
     182    if call_details.num_active == 0: 
     183        norm = float(form_volume()) 
     184        if norm > 0.0: 
     185            return (scale/norm)*form() + background 
     186        else: 
     187            return np.ones(nq, 'd')*background 
     188 
     189    partial_weight = np.NaN 
     190    spherical_correction = 1.0 
     191    pd_stride = call_details.pd_stride[:call_details.num_active] 
     192    pd_length = call_details.pd_length[:call_details.num_active] 
     193    pd_offset = call_details.pd_offset[:call_details.num_active] 
     194    pd_index = np.empty_like(pd_offset) 
     195    offset = np.empty_like(call_details.par_offset) 
     196    theta = call_details.theta_par 
     197    fast_length = pd_length[0] 
     198    pd_index[0] = fast_length 
     199    total = np.zeros(nq, 'd') 
     200    norm = 0.0 
     201    for loop_index in range(call_details.total_pd): 
     202        # update polydispersity parameter values 
     203        if pd_index[0] == fast_length: 
     204            pd_index[:] = (loop_index/pd_stride)%pd_length 
     205            partial_weight = np.prod(weights[pd_offset+pd_index][1:]) 
     206            for k in range(call_details.num_coord): 
     207                par = call_details.par_coord[k] 
     208                coord = call_details.pd_coord[k] 
     209                this_offset = call_details.par_offset[par] 
     210                block_size = 1 
     211                for bit in range(len(pd_offset)): 
     212                    if coord&1: 
     213                        this_offset += block_size * pd_index[bit] 
     214                        block_size *= pd_length[bit] 
     215                    coord >>= 1 
     216                    if coord == 0: break 
     217                offset[par] = this_offset 
     218                parameters[par] = values[this_offset] 
     219                if par == theta and not (call_details.par_coord[k]&1): 
     220                    spherical_correction = max(abs(cos(pi/180 * parameters[theta])), 1e-6) 
     221        for k in range(call_details.num_coord): 
     222            if call_details.pd_coord[k]&1: 
     223                #par = call_details.par_coord[k] 
     224                parameters[par] = values[offset[par]] 
     225                #print "par",par,offset[par],parameters[par+2] 
     226                offset[par] += 1 
     227                if par == theta: 
     228                    spherical_correction = max(abs(cos(pi/180 * parameters[theta])), 1e-6) 
     229 
     230        weight = partial_weight * weights[pd_offset[0] + pd_index[0]] 
     231        pd_index[0] += 1 
     232        if weight > cutoff: 
     233            # Call the scattering function 
     234            # Assume that NaNs are only generated if the parameters are bad; 
     235            # exclude all q for that NaN.  Even better would be to have an 
     236            # INVALID expression like the C models, but that is too expensive. 
     237            I = form() 
     238            if np.isnan(I).any(): continue 
     239 
     240            # update value and norm 
     241            weight *= spherical_correction 
     242            total += weight * I 
     243            norm += weight * form_volume() 
     244 
     245    if norm > 0.0: 
     246        return (scale/norm)*total + background 
    206247    else: 
    207         stride = np.array([1]) 
    208         vol_weight_index = slice(None, None) 
    209         # keep lint happy 
    210         fast_value = [None] 
    211         fast_weight = [None] 
    212  
    213     ret = np.zeros_like(args[0]) 
    214     norm = np.zeros_like(ret) 
    215     vol = np.zeros_like(ret) 
    216     vol_norm = np.zeros_like(ret) 
    217     for k in range(stride[-1]): 
    218         # update polydispersity parameter values 
    219         fast_index = k % stride[0] 
    220         if fast_index == 0:  # bottom loop complete ... check all other loops 
    221             if weight.size > 0: 
    222                 for i, index, in enumerate(k % stride): 
    223                     args[pd_index[i]] = pd[i][0][index] 
    224                     weight[i] = pd[i][1][index] 
    225         else: 
    226             args[pd_index[0]] = fast_value[fast_index] 
    227             weight[0] = fast_weight[fast_index] 
    228  
    229         # Computes the weight, and if it is not sufficient then ignore this 
    230         # parameter set. 
    231         # Note: could precompute w1*...*wn so we only need to multiply by w0 
    232         w = np.prod(weight) 
    233         if w > cutoff: 
    234             # Note: can precompute spherical correction if theta_index is not 
    235             # the fast index. Correction factor for spherical integration 
    236             #spherical_correction = abs(cos(pi*args[phi_index])) if phi_index>=0 else 1.0 
    237             spherical_correction = (abs(cos(pi * args[theta_index])) * pi / 2 
    238                                     if theta_index >= 0 else 1.0) 
    239             #spherical_correction = 1.0 
    240  
    241             # Call the scattering function and adds it to the total. 
    242             I = form(*args) 
    243             if np.isnan(I).any(): continue 
    244             ret += w * I * spherical_correction 
    245             norm += w 
    246  
    247             # Volume normalization. 
    248             # If there are "volume" polydispersity parameters, then these 
    249             # will be used to call the form_volume function from the user 
    250             # supplied kernel, and accumulate a normalized weight. 
    251             # Note: can precompute volume norm if fast index is not a volume 
    252             if form_volume: 
    253                 vol_args = [args[index] for index in vol_index] 
    254                 vol_weight = np.prod(weight[vol_weight_index]) 
    255                 vol += vol_weight * form_volume(*vol_args) 
    256                 vol_norm += vol_weight 
    257  
    258     positive = (vol * vol_norm != 0.0) 
    259     ret[positive] *= vol_norm[positive] / vol[positive] 
    260     result = scale * ret / norm + background 
    261     return result 
     248        return np.ones(nq, 'd')*background 
     249 
     250 
     251def _create_default_functions(model_info): 
     252    """ 
     253    Autogenerate missing functions, such as Iqxy from Iq. 
     254 
     255    This only works for Iqxy when Iq is written in python. :func:`make_source` 
     256    performs a similar role for Iq written in C.  This also vectorizes 
     257    any functions that are not already marked as vectorized. 
     258    """ 
     259    _create_vector_Iq(model_info) 
     260    _create_vector_Iqxy(model_info)  # call create_vector_Iq() first 
     261 
     262 
     263def _create_vector_Iq(model_info): 
     264    """ 
     265    Define Iq as a vector function if it exists. 
     266    """ 
     267    Iq = model_info.Iq 
     268    if callable(Iq) and not getattr(Iq, 'vectorized', False): 
     269        #print("vectorizing Iq") 
     270        def vector_Iq(q, *args): 
     271            """ 
     272            Vectorized 1D kernel. 
     273            """ 
     274            return np.array([Iq(qi, *args) for qi in q]) 
     275        vector_Iq.vectorized = True 
     276        model_info.Iq = vector_Iq 
     277 
     278def _create_vector_Iqxy(model_info): 
     279    """ 
     280    Define Iqxy as a vector function if it exists, or default it from Iq(). 
     281    """ 
     282    Iq, Iqxy = model_info.Iq, model_info.Iqxy 
     283    if callable(Iqxy): 
     284        if not getattr(Iqxy, 'vectorized', False): 
     285            #print("vectorizing Iqxy") 
     286            def vector_Iqxy(qx, qy, *args): 
     287                """ 
     288                Vectorized 2D kernel. 
     289                """ 
     290                return np.array([Iqxy(qxi, qyi, *args) for qxi, qyi in zip(qx, qy)]) 
     291            vector_Iqxy.vectorized = True 
     292            model_info.Iqxy = vector_Iqxy 
     293    elif callable(Iq): 
     294        #print("defaulting Iqxy") 
     295        # Iq is vectorized because create_vector_Iq was already called. 
     296        def default_Iqxy(qx, qy, *args): 
     297            """ 
     298            Default 2D kernel. 
     299            """ 
     300            return Iq(np.sqrt(qx**2 + qy**2), *args) 
     301        default_Iqxy.vectorized = True 
     302        model_info.Iqxy = default_Iqxy 
     303 

  • sasmodels/list_pars.py

    ra84a0ca r6d6508e  
    2525    for name in sorted(MODELS): 
    2626        model_info = load_model_info(name) 
    27         for p in model_info['parameters']: 
     27        for p in model_info.parameters.kernel_parameters: 
    2828            partable.setdefault(p.name, []) 
    2929            partable[p.name].append(name) 

  • sasmodels/mixture.py

    r72a081d r7ae2b7f  
    1212""" 
    1313from copy import copy 
    14 import numpy as np 
     14import numpy as np  # type: ignore 
    1515 
    16 from .generate import process_parameters, COMMON_PARAMETERS, Parameter 
     16from .modelinfo import Parameter, ParameterTable, ModelInfo 
     17from .kernel import KernelModel, Kernel 
    1718 
    18 SCALE=0 
    19 BACKGROUND=1 
    20 EFFECT_RADIUS=2 
    21 VOLFRACTION=3 
     19try: 
     20    from typing import List 
     21    from .details import CallDetails 
     22except ImportError: 
     23    pass 
    2224 
    2325def make_mixture_info(parts): 
     26    # type: (List[ModelInfo]) -> ModelInfo 
    2427    """ 
    2528    Create info block for product model. 
    … …  
    2730    flatten = [] 
    2831    for part in parts: 
    29         if part['composition'] and part['composition'][0] == 'mixture': 
    30             flatten.extend(part['compostion'][1]) 
     32        if part.composition and part.composition[0] == 'mixture': 
     33            flatten.extend(part.composition[1]) 
    3134        else: 
    3235            flatten.append(part) 
    … …  
    3437 
    3538    # Build new parameter list 
    36     pars = COMMON_PARAMETERS + [] 
     39    combined_pars = [] 
    3740    for k, part in enumerate(parts): 
    3841        # Parameter prefix per model, A_, B_, ... 
     42        # Note that prefix must also be applied to id and length_control 
     43        # to support vector parameters 
    3944        prefix = chr(ord('A')+k) + '_' 
    40         for p in part['parameters']: 
    41             # No background on the individual mixture elements 
    42             if p.name == 'background': 
    43                 continue 
    44             # TODO: promote Parameter to a full class 
    45             # this code knows too much about the implementation! 
    46             p = list(p) 
    47             if p[0] == 'scale':  # allow model subtraction 
    48                 p[3] = [-np.inf, np.inf]  # limits 
    49             p[0] = prefix+p[0]   # relabel parameters with A_, ... 
    50             pars.append(p) 
     45        combined_pars.append(Parameter(prefix+'scale')) 
     46        for p in part.parameters.kernel_parameters: 
     47            p = copy(p) 
     48            p.name = prefix + p.name 
     49            p.id = prefix + p.id 
     50            if p.length_control is not None: 
     51                p.length_control = prefix + p.length_control 
     52            combined_pars.append(p) 
     53    parameters = ParameterTable(combined_pars) 
    5154 
    52     model_info = {} 
    53     model_info['id'] = '+'.join(part['id']) 
    54     model_info['name'] = ' + '.join(part['name']) 
    55     model_info['filename'] = None 
    56     model_info['title'] = 'Mixture model with ' + model_info['name'] 
    57     model_info['description'] = model_info['title'] 
    58     model_info['docs'] = model_info['title'] 
    59     model_info['category'] = "custom" 
    60     model_info['parameters'] = pars 
    61     #model_info['single'] = any(part['single'] for part in parts) 
    62     model_info['structure_factor'] = False 
    63     model_info['variant_info'] = None 
    64     #model_info['tests'] = [] 
    65     #model_info['source'] = [] 
     55    model_info = ModelInfo() 
     56    model_info.id = '+'.join(part.id for part in parts) 
     57    model_info.name = ' + '.join(part.name for part in parts) 
     58    model_info.filename = None 
     59    model_info.title = 'Mixture model with ' + model_info.name 
     60    model_info.description = model_info.title 
     61    model_info.docs = model_info.title 
     62    model_info.category = "custom" 
     63    model_info.parameters = parameters 
     64    #model_info.single = any(part['single'] for part in parts) 
     65    model_info.structure_factor = False 
     66    model_info.variant_info = None 
     67    #model_info.tests = [] 
     68    #model_info.source = [] 
    6669    # Iq, Iqxy, form_volume, ER, VR and sesans 
    6770    # Remember the component info blocks so we can build the model 
    68     model_info['composition'] = ('mixture', parts) 
    69     process_parameters(model_info) 
    70     return model_info 
     71    model_info.composition = ('mixture', parts) 
    7172 
    7273 
    73 class MixtureModel(object): 
     74class MixtureModel(KernelModel): 
    7475    def __init__(self, model_info, parts): 
     76        # type: (ModelInfo, List[KernelModel]) -> None 
    7577        self.info = model_info 
    7678        self.parts = parts 
    7779 
    7880    def __call__(self, q_vectors): 
     81        # type: (List[np.ndarray]) -> MixtureKernel 
    7982        # Note: may be sending the q_vectors to the n times even though they 
    8083        # are only needed once.  It would mess up modularity quite a bit to 
    … …  
    8386        # in opencl; or both in opencl, but one in single precision and the 
    8487        # other in double precision). 
    85         kernels = [part(q_vectors) for part in self.parts] 
     88        kernels = [part.make_kernel(q_vectors) for part in self.parts] 
    8689        return MixtureKernel(self.info, kernels) 
    8790 
    8891    def release(self): 
     92        # type: () -> None 
    8993        """ 
    9094        Free resources associated with the model. 
    … …  
    9498 
    9599 
    96 class MixtureKernel(object): 
     100class MixtureKernel(Kernel): 
    97101    def __init__(self, model_info, kernels): 
    98         dim = '2d' if kernels[0].q_input.is_2d else '1d' 
     102        # type: (ModelInfo, List[Kernel]) -> None 
     103        self.dim = kernels[0].dim 
     104        self.info =  model_info 
     105        self.kernels = kernels 
    99106 
    100         # fixed offsets starts at 2 for scale and background 
    101         fixed_pars, pd_pars = [], [] 
    102         offsets = [[2, 0]] 
    103         #vol_index = [] 
    104         def accumulate(fixed, pd, volume): 
    105             # subtract 1 from fixed since we are removing background 
    106             fixed_offset, pd_offset = offsets[-1] 
    107             #vol_index.extend(k+pd_offset for k,v in pd if v in volume) 
    108             offsets.append([fixed_offset + len(fixed) - 1, pd_offset + len(pd)]) 
    109             pd_pars.append(pd) 
    110         if dim == '2d': 
    111             for p in kernels: 
    112                 partype = p.info['partype'] 
    113                 accumulate(partype['fixed-2d'], partype['pd-2d'], partype['volume']) 
    114         else: 
    115             for p in kernels: 
    116                 partype = p.info['partype'] 
    117                 accumulate(partype['fixed-1d'], partype['pd-1d'], partype['volume']) 
    118  
    119         #self.vol_index = vol_index 
    120         self.offsets = offsets 
    121         self.fixed_pars = fixed_pars 
    122         self.pd_pars = pd_pars 
    123         self.info = model_info 
    124         self.kernels = kernels 
    125         self.results = None 
    126  
    127     def __call__(self, fixed_pars, pd_pars, cutoff=1e-5): 
    128         scale, background = fixed_pars[0:2] 
     107    def __call__(self, call_details, value, weight, cutoff): 
     108        # type: (CallDetails, np.ndarray, np.ndarry, float) -> np.ndarray 
     109        scale, background = value[0:2] 
    129110        total = 0.0 
    130         self.results = []  # remember the parts for plotting later 
    131         for k in range(len(self.offsets)-1): 
    132             start_fixed, start_pd = self.offsets[k] 
    133             end_fixed, end_pd = self.offsets[k+1] 
    134             part_fixed = [fixed_pars[start_fixed], 0.0] + fixed_pars[start_fixed+1:end_fixed] 
    135             part_pd = [pd_pars[start_pd], 0.0] + pd_pars[start_pd+1:end_pd] 
    136             part_result = self.kernels[k](part_fixed, part_pd) 
     111        # remember the parts for plotting later 
     112        self.results = [] 
     113        for kernel, kernel_details in zip(self.kernels, call_details.parts): 
     114            part_result = kernel(kernel_details, value, weight, cutoff) 
    137115            total += part_result 
    138             self.results.append(scale*sum+background) 
     116            self.results.append(part_result) 
    139117 
    140118        return scale*total + background 
    141119 
    142120    def release(self): 
    143         self.p_kernel.release() 
    144         self.q_kernel.release() 
     121        # type: () -> None 
     122        for k in self.kernels: 
     123            k.release() 
    145124 

  • sasmodels/model_test.py

    r4d76711 r7ae2b7f  
    4343Precision defaults to 5 digits (relative). 
    4444""" 
     45#TODO: rename to tests so that tab completion works better for models directory 
     46 
    4547from __future__ import print_function 
    4648 
    … …  
    4850import unittest 
    4951 
    50 import numpy as np 
     52import numpy as np  # type: ignore 
    5153 
    5254from .core import list_models, load_model_info, build_model, HAVE_OPENCL 
    53 from .core import call_kernel, call_ER, call_VR 
     55from .details import dispersion_mesh 
     56from .direct_model import call_kernel, get_weights 
    5457from .exception import annotate_exception 
    55  
    56 #TODO: rename to tests so that tab completion works better for models directory 
     58from .modelinfo import expand_pars 
     59 
     60try: 
     61    from typing import List, Iterator, Callable 
     62except ImportError: 
     63    pass 
     64else: 
     65    from .modelinfo import ParameterTable, ParameterSet, TestCondition, ModelInfo 
     66    from .kernelpy import PyModel, PyInput, PyKernel, DType 
     67 
     68def call_ER(model_info, pars): 
     69    # type: (ModelInfo, ParameterSet) -> float 
     70    """ 
     71    Call the model ER function using *values*. 
     72 
     73    *model_info* is either *model.info* if you have a loaded model, 
     74    or *kernel.info* if you have a model kernel prepared for evaluation. 
     75    """ 
     76    if model_info.ER is None: 
     77        return 1.0 
     78    else: 
     79        value, weight = _vol_pars(model_info, pars) 
     80        individual_radii = model_info.ER(*value) 
     81        return np.sum(weight*individual_radii) / np.sum(weight) 
     82 
     83def call_VR(model_info, pars): 
     84    # type: (ModelInfo, ParameterSet) -> float 
     85    """ 
     86    Call the model VR function using *pars*. 
     87 
     88    *model_info* is either *model.info* if you have a loaded model, 
     89    or *kernel.info* if you have a model kernel prepared for evaluation. 
     90    """ 
     91    if model_info.VR is None: 
     92        return 1.0 
     93    else: 
     94        value, weight = _vol_pars(model_info, pars) 
     95        whole, part = model_info.VR(*value) 
     96        return np.sum(weight*part)/np.sum(weight*whole) 
     97 
     98def _vol_pars(model_info, pars): 
     99    # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray] 
     100    vol_pars = [get_weights(p, pars) 
     101                for p in model_info.parameters.call_parameters 
     102                if p.type == 'volume'] 
     103    value, weight = dispersion_mesh(model_info, vol_pars) 
     104    return value, weight 
     105 
    57106 
    58107def make_suite(loaders, models): 
     108    # type: (List[str], List[str]) -> unittest.TestSuite 
    59109    """ 
    60110    Construct the pyunit test suite. 
    … …  
    66116    *models* is the list of models to test, or *["all"]* to test all models. 
    67117    """ 
    68  
    69     ModelTestCase = _hide_model_case_from_nosetests() 
     118    ModelTestCase = _hide_model_case_from_nose() 
    70119    suite = unittest.TestSuite() 
    71120 
    … …  
    86135        # don't try to call cl kernel since it will not be 
    87136        # available in some environmentes. 
    88         is_py = callable(model_info['Iq']) 
     137        is_py = callable(model_info.Iq) 
    89138 
    90139        if is_py:  # kernel implemented in python 
    … …  
    124173 
    125174 
    126 def _hide_model_case_from_nosetests(): 
     175def _hide_model_case_from_nose(): 
     176    # type: () -> type 
    127177    class ModelTestCase(unittest.TestCase): 
    128178        """ 
    … …  
    135185        def __init__(self, test_name, model_info, test_method_name, 
    136186                     platform, dtype): 
     187            # type: (str, ModelInfo, str, str, DType) -> None 
    137188            self.test_name = test_name 
    138189            self.info = model_info 
    … …  
    140191            self.dtype = dtype 
    141192 
    142             setattr(self, test_method_name, self._runTest) 
     193            setattr(self, test_method_name, self.run_all) 
    143194            unittest.TestCase.__init__(self, test_method_name) 
    144195 
    145         def _runTest(self): 
     196        def run_all(self): 
     197            # type: () -> None 
    146198            smoke_tests = [ 
    147                 [{}, 0.1, None], 
    148                 [{}, (0.1, 0.1), None], 
    149                 [{}, 'ER', None], 
    150                 [{}, 'VR', None], 
     199                ({}, 0.1, None), 
     200                ({}, (0.1, 0.1), None), 
     201                ({}, 'ER', None), 
     202                ({}, 'VR', None), 
    151203                ] 
    152204 
    153             tests = self.info['tests'] 
     205            tests = self.info.tests 
    154206            try: 
    155207                model = build_model(self.info, dtype=self.dtype, 
    156208                                    platform=self.platform) 
    157209                for test in smoke_tests + tests: 
    158                     self._run_one_test(model, test) 
     210                    self.run_one(model, test) 
    159211 
    160212                if not tests and self.platform == "dll": 
    … …  
    170222                raise 
    171223 
    172         def _run_one_test(self, model, test): 
    173             pars, x, y = test 
     224        def run_one(self, model, test): 
     225            # type: (PyModel, TestCondition) -> None 
     226            user_pars, x, y = test 
     227            pars = expand_pars(self.info.parameters, user_pars) 
    174228 
    175229            if not isinstance(y, list): 
    … …  
    185239                actual = [call_VR(model.info, pars)] 
    186240            elif isinstance(x[0], tuple): 
    187                 Qx, Qy = zip(*x) 
    188                 q_vectors = [np.array(Qx), np.array(Qy)] 
     241                qx, qy = zip(*x) 
     242                q_vectors = [np.array(qx), np.array(qy)] 
    189243                kernel = model.make_kernel(q_vectors) 
    190244                actual = call_kernel(kernel, pars) 
    … …  
    194248                actual = call_kernel(kernel, pars) 
    195249 
    196             self.assertGreater(len(actual), 0) 
     250            self.assertTrue(len(actual) > 0) 
    197251            self.assertEqual(len(y), len(actual)) 
    198252 
    … …  
    210264 
    211265def is_near(target, actual, digits=5): 
     266    # type: (float, float, int) -> bool 
    212267    """ 
    213268    Returns true if *actual* is within *digits* significant digits of *target*. 
    … …  
    218273 
    219274def main(): 
     275    # type: () -> int 
    220276    """ 
    221277    Run tests given is sys.argv. 
    … …  
    251307  python -m sasmodels.model_test [-v] [opencl|dll] model1 model2 ... 
    252308 
    253 If -v is included on the command line, then use verboe output. 
     309If -v is included on the command line, then use verbose output. 
    254310 
    255311If neither opencl nor dll is specified, then models will be tested with 
    256 both opencl and dll; the compute target is ignored for pure python models. 
     312both OpenCL and dll; the compute target is ignored for pure python models. 
    257313 
    258314If model1 is 'all', then all except the remaining models will be tested. 
    … …  
    269325 
    270326def model_tests(): 
     327    # type: () -> Iterator[Callable[[], None]] 
    271328    """ 
    272329    Test runner visible to nosetests. 
    … …  
    276333    tests = make_suite(['opencl', 'dll'], ['all']) 
    277334    for test_i in tests: 
    278         yield test_i._runTest 
     335        yield test_i.run_all 
    279336 
    280337 

  • sasmodels/models/cylinder.c

    r26141cb r03cac08  
    33double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    44    double radius, double length, double theta, double phi); 
     5 
     6#define INVALID(v) (v.radius<0 || v.length<0) 
    57 
    68double form_volume(double radius, double length) 
    … …  
    1517    double length) 
    1618{ 
    17     // TODO: return NaN if radius<0 or length<0? 
    1819    // precompute qr and qh to save time in the loop 
    1920    const double qr = q*radius; 
    … …  
    4748    double phi) 
    4849{ 
    49     // TODO: return NaN if radius<0 or length<0? 
    5050    double sn, cn; // slots to hold sincos function output 
    5151 

  • sasmodels/models/cylinder.py

    rf247314 r7ae2b7f  
    8282""" 
    8383 
    84 import numpy as np 
    85 from numpy import pi, inf 
     84import numpy as np  # type: ignore 
     85from numpy import pi, inf  # type: ignore 
    8686 
    8787name = "cylinder" 

  • sasmodels/models/flexible_cylinder.c

    re6408d0 re6408d0  
    1 double form_volume(double length, double kuhn_length, double radius); 
    2 double Iq(double q, double length, double kuhn_length, double radius, 
    3           double sld, double solvent_sld); 
    4 double Iqxy(double qx, double qy, double length, double kuhn_length, 
    5             double radius, double sld, double solvent_sld); 
    6 double flexible_cylinder_kernel(double q, double length, double kuhn_length, 
    7                                 double radius, double sld, double solvent_sld); 
    8  
    9  
    10 double form_volume(double length, double kuhn_length, double radius) 
     1static double 
     2form_volume(length, kuhn_length, radius) 
    113{ 
    124    return 1.0; 
    135} 
    146 
    15 double flexible_cylinder_kernel(double q, 
    16           double length, 
    17           double kuhn_length, 
    18           double radius, 
    19           double sld, 
    20           double solvent_sld) 
     7static double 
     8Iq(double q, 
     9   double length, 
     10   double kuhn_length, 
     11   double radius, 
     12   double sld, 
     13   double solvent_sld) 
    2114{ 
    22  
    23     const double cont = sld-solvent_sld; 
    24     const double qr = q*radius; 
    25     //const double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr); 
    26     const double crossSect = sas_J1c(qr); 
    27     double flex = Sk_WR(q,length,kuhn_length); 
    28     flex *= crossSect*crossSect; 
    29     flex *= M_PI*radius*radius*length; 
    30     flex *= cont*cont; 
    31     flex *= 1.0e-4; 
    32     return flex; 
     15    const double contrast = sld - solvent_sld; 
     16    const double cross_section = sas_J1c(q*radius); 
     17    const double volume = M_PI*radius*radius*length; 
     18    const double flex = Sk_WR(q, length, kuhn_length); 
     19    return 1.0e-4 * volume * square(contrast*cross_section) * flex; 
    3320} 
    34  
    35 double Iq(double q, 
    36           double length, 
    37           double kuhn_length, 
    38           double radius, 
    39           double sld, 
    40           double solvent_sld) 
    41 { 
    42  
    43     double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld); 
    44     return result; 
    45 } 
    46  
    47 double Iqxy(double qx, double qy, 
    48             double length, 
    49             double kuhn_length, 
    50             double radius, 
    51             double sld, 
    52             double solvent_sld) 
    53 { 
    54     double q; 
    55     q = sqrt(qx*qx+qy*qy); 
    56     double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld); 
    57  
    58     return result; 
    59 } 

  • sasmodels/models/gel_fit.c

    r30b4ddf r03cac08  
    1 double form_volume(void); 
    2  
    3 double Iq(double q, 
    4           double guinier_scale, 
    5           double lorentzian_scale, 
    6           double gyration_radius, 
    7           double fractal_exp, 
    8           double cor_length); 
    9  
    10 double Iqxy(double qx, double qy, 
    11           double guinier_scale, 
    12           double lorentzian_scale, 
    13           double gyration_radius, 
    14           double fractal_exp, 
    15           double cor_length); 
    16  
    17 static double _gel_fit_kernel(double q, 
     1static double Iq(double q, 
    182          double guinier_scale, 
    193          double lorentzian_scale, 
    … …  
    248    // Lorentzian Term 
    259    ////////////////////////double a(x[i]*x[i]*zeta*zeta); 
    26     double lorentzian_term = q*q*cor_length*cor_length; 
     10    double lorentzian_term = square(q*cor_length); 
    2711    lorentzian_term = 1.0 + ((fractal_exp + 1.0)/3.0)*lorentzian_term; 
    2812    lorentzian_term = pow(lorentzian_term, (fractal_exp/2.0) ); 
    … …  
    3014    // Exponential Term 
    3115    ////////////////////////double d(x[i]*x[i]*rg*rg); 
    32     double exp_term = q*q*gyration_radius*gyration_radius; 
     16    double exp_term = square(q*gyration_radius); 
    3317    exp_term = exp(-1.0*(exp_term/3.0)); 
    3418 
    … …  
    3721    return result; 
    3822} 
    39 double form_volume(void){ 
    40     // Unused, so free to return garbage. 
    41     return NAN; 
    42 } 
    43  
    44 double Iq(double q, 
    45           double guinier_scale, 
    46           double lorentzian_scale, 
    47           double gyration_radius, 
    48           double fractal_exp, 
    49           double cor_length) 
    50 { 
    51     return _gel_fit_kernel(q, 
    52                           guinier_scale, 
    53                           lorentzian_scale, 
    54                           gyration_radius, 
    55                           fractal_exp, 
    56                           cor_length); 
    57 } 
    58  
    59 // Iqxy is never called since no orientation or magnetic parameters. 
    60 double Iqxy(double qx, double qy, 
    61           double guinier_scale, 
    62           double lorentzian_scale, 
    63           double gyration_radius, 
    64           double fractal_exp, 
    65           double cor_length) 
    66 { 
    67     double q = sqrt(qx*qx + qy*qy); 
    68     return _gel_fit_kernel(q, 
    69                           guinier_scale, 
    70                           lorentzian_scale, 
    71                           gyration_radius, 
    72                           fractal_exp, 
    73                           cor_length); 
    74 } 
    75  

  • sasmodels/models/hardsphere.py

    rec45c4f rd2bb604  
    149149   """ 
    150150 
    151 Iqxy = """ 
    152     // never called since no orientation or magnetic parameters. 
    153     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    154     """ 
    155  
    156151# ER defaults to 0.0 
    157152# VR defaults to 1.0 

  • sasmodels/models/hayter_msa.py

    rec45c4f rd2bb604  
    8787    return 1.0; 
    8888    """ 
    89 Iqxy = """ 
    90     // never called since no orientation or magnetic parameters. 
    91     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    92     """ 
    9389# ER defaults to 0.0 
    9490# VR defaults to 1.0 

  • sasmodels/models/lamellar.py

    rec45c4f rd2bb604  
    8282    """ 
    8383 
    84 Iqxy = """ 
    85     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    86     """ 
    87  
    8884# ER defaults to 0.0 
    8985# VR defaults to 1.0 

  • sasmodels/models/lamellar_hg.py

    rec45c4f rd2bb604  
    101101    """ 
    102102 
    103 Iqxy = """ 
    104     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    105     """ 
    106  
    107103# ER defaults to 0.0 
    108104# VR defaults to 1.0 

  • sasmodels/models/lamellar_hg_stack_caille.py

    rec45c4f rd2bb604  
    120120    """ 
    121121 
    122 Iqxy = """ 
    123     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    124     """ 
    125  
    126122# ER defaults to 0.0 
    127123# VR defaults to 1.0 

  • sasmodels/models/lamellar_stack_caille.py

    rec45c4f rd2bb604  
    104104    """ 
    105105 
    106 Iqxy = """ 
    107     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    108     """ 
    109  
    110106# ER defaults to 0.0 
    111107# VR defaults to 1.0 

  • sasmodels/models/lamellar_stack_paracrystal.py

    rec45c4f rd2bb604  
    132132    """ 
    133133 
    134 Iqxy = """ 
    135     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    136     """ 
    137  
    138134# ER defaults to 0.0 
    139135# VR defaults to 1.0 

  • sasmodels/models/lib/sas_JN.c

    re6408d0 re6408d0  
    4848*/ 
    4949 
    50 static double 
    51 sas_JN( int n, double x ) { 
     50double sas_JN( int n, double x ); 
     51 
     52double sas_JN( int n, double x ) { 
    5253 
    5354    const double MACHEP = 1.11022302462515654042E-16; 

  • sasmodels/models/lib/sph_j1c.c

    re6f1410 rba32cdd  
    77* using double precision that are the source. 
    88*/ 
     9double sph_j1c(double q); 
    910 
    1011// The choice of the number of terms in the series and the cutoff value for 
    … …  
    4344#endif 
    4445 
    45 double sph_j1c(double q); 
    4646double sph_j1c(double q) 
    4747{ 

  • sasmodels/models/lib/sphere_form.c

    rad90df9 rba32cdd  
    1 inline double 
    2 sphere_volume(double radius) 
     1double sphere_volume(double radius); 
     2double sphere_form(double q, double radius, double sld, double solvent_sld); 
     3 
     4double sphere_volume(double radius) 
    35{ 
    46    return M_4PI_3*cube(radius); 
    57} 
    68 
    7 inline double 
    8 sphere_form(double q, double radius, double sld, double solvent_sld) 
     9double sphere_form(double q, double radius, double sld, double solvent_sld) 
    910{ 
    1011    const double fq = sphere_volume(radius) * sph_j1c(q*radius); 

  • sasmodels/models/lib/wrc_cyl.c

    re7678b2 rba32cdd  
    22    Functions for WRC implementation of flexible cylinders 
    33*/ 
     4double Sk_WR(double q, double L, double b); 
     5 
     6 
    47static double 
    58AlphaSquare(double x) 
    … …  
    363366} 
    364367 
    365 double Sk_WR(double q, double L, double b); 
    366368double Sk_WR(double q, double L, double b) 
    367369{ 

  • sasmodels/models/onion.c

    rfdb1487 rce896fd  
    44    double thickness, double A) 
    55{ 
    6   const double vol = 4.0/3.0 * M_PI * r * r * r; 
     6  const double vol = M_4PI_3 * cube(r); 
    77  const double qr = q * r; 
    88  const double alpha = A * r/thickness; 
    … …  
    1919    double sinqr, cosqr; 
    2020    SINCOS(qr, sinqr, cosqr); 
    21     fun = -3.0( 
     21    fun = -3.0*( 
    2222            ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq 
    2323                - (alpha*sinqr/qr - cosqr) 
    … …  
    3232f_linear(double q, double r, double sld, double slope) 
    3333{ 
    34   const double vol = 4.0/3.0 * M_PI * r * r * r; 
     34  const double vol = M_4PI_3 * cube(r); 
    3535  const double qr = q * r; 
    3636  const double bes = sph_j1c(qr); 
    … …  
    5252{ 
    5353  const double bes = sph_j1c(q * r); 
    54   const double vol = 4.0/3.0 * M_PI * r * r * r; 
     54  const double vol = M_4PI_3 * cube(r); 
    5555  return sld * vol * bes; 
    5656} 
    … …  
    6464    r += thickness[i]; 
    6565  } 
    66   return 4.0/3.0 * M_PI * r * r * r; 
     66  return M_4PI_3*cube(r); 
    6767} 
    6868 
    6969static double 
    70 Iq(double q, double core_sld, double core_radius, double solvent_sld, 
    71     double n, double in_sld[], double out_sld[], double thickness[], 
     70Iq(double q, double sld_core, double core_radius, double sld_solvent, 
     71    double n, double sld_in[], double sld_out[], double thickness[], 
    7272    double A[]) 
    7373{ 
    7474  int i; 
    75   r = core_radius; 
    76   f = f_constant(q, r, core_sld); 
     75  double r = core_radius; 
     76  double f = f_constant(q, r, sld_core); 
    7777  for (i=0; i<n; i++){ 
    7878    const double r0 = r; 
    … …  
    9292    } 
    9393  } 
    94   f -= f_constant(q, r, solvent_sld); 
    95   f2 = f * f * 1.0e-4; 
     94  f -= f_constant(q, r, sld_solvent); 
     95  const double f2 = f * f * 1.0e-4; 
    9696 
    9797  return f2; 

  • sasmodels/models/onion.py

    r416609b rfa5fd8d  
    293293 
    294294#             ["name", "units", default, [lower, upper], "type","description"], 
    295 parameters = [["core_sld", "1e-6/Ang^2", 1.0, [-inf, inf], "", 
     295parameters = [["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "", 
    296296               "Core scattering length density"], 
    297297              ["core_radius", "Ang", 200., [0, inf], "volume", 
    298298               "Radius of the core"], 
    299               ["solvent_sld", "1e-6/Ang^2", 6.4, [-inf, inf], "", 
     299              ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "", 
    300300               "Solvent scattering length density"], 
    301               ["n", "", 1, [0, 10], "volume", 
     301              ["n_shells", "", 1, [0, 10], "volume", 
    302302               "number of shells"], 
    303               ["in_sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "", 
     303              ["sld_in[n_shells]", "1e-6/Ang^2", 1.7, [-inf, inf], "", 
    304304               "scattering length density at the inner radius of shell k"], 
    305               ["out_sld[n]", "1e-6/Ang^2", 2.0, [-inf, inf], "", 
     305              ["sld_out[n_shells]", "1e-6/Ang^2", 2.0, [-inf, inf], "", 
    306306               "scattering length density at the outer radius of shell k"], 
    307               ["thickness[n]", "Ang", 40., [0, inf], "volume", 
     307              ["thickness[n_shells]", "Ang", 40., [0, inf], "volume", 
    308308               "Thickness of shell k"], 
    309               ["A[n]", "", 1.0, [-inf, inf], "", 
     309              ["A[n_shells]", "", 1.0, [-inf, inf], "", 
    310310               "Decay rate of shell k"], 
    311311              ] 
    312312 
    313 #source = ["lib/sph_j1c.c", "onion.c"] 
    314  
    315 def Iq(q, *args, **kw): 
    316     return q 
    317  
    318 def Iqxy(qx, *args, **kw): 
    319     return qx 
    320  
    321  
    322 def shape(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A): 
     313source = ["lib/sph_j1c.c", "onion.c"] 
     314 
     315#def Iq(q, *args, **kw): 
     316#    return q 
     317 
     318profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 
     319def profile(core_sld, core_radius, solvent_sld, n_shells, 
     320            in_sld, out_sld, thickness, A): 
    323321    """ 
    324     SLD profile 
     322    Returns shape profile with x=radius, y=SLD. 
    325323    """ 
    326324 
    327     total_radius = 1.25*(sum(thickness[:n]) + core_radius + 1) 
     325    total_radius = 1.25*(sum(thickness[:n_shells]) + core_radius + 1) 
    328326    dr = total_radius/400  # 400 points for a smooth plot 
    329327 
    … …  
    338336 
    339337    # add in the shells 
    340     for k in range(n): 
     338    for k in range(n_shells): 
    341339        # Left side of each shells 
    342340        r0 = r[-1] 
    … …  
    365363    beta.append(solvent_sld) 
    366364 
    367     return np.asarray(r), np.asarray(beta) 
     365    return np.asarray(r), np.asarray(beta)*1e-6 
    368366 
    369367def ER(core_radius, n, thickness): 
    … …  
    374372 
    375373demo = { 
    376     "solvent_sld": 2.2, 
    377     "core_sld": 1.0, 
     374    "sld_solvent": 2.2, 
     375    "sld_core": 1.0, 
    378376    "core_radius": 100, 
    379377    "n": 4, 
    380     "in_sld": [0.5, 1.5, 0.9, 2.0], 
    381     "out_sld": [nan, 0.9, 1.2, 1.6], 
     378    "sld_in": [0.5, 1.5, 0.9, 2.0], 
     379    "sld_out": [nan, 0.9, 1.2, 1.6], 
    382380    "thickness": [50, 75, 150, 75], 
    383381    "A": [0, -1, 1e-4, 1], 

  • sasmodels/models/rpa.c

    r13ed84c rd2bb604  
    11double Iq(double q, double case_num, 
    2     double Na, double Phia, double va, double a_sld, double ba, 
    3     double Nb, double Phib, double vb, double b_sld, double bb, 
    4     double Nc, double Phic, double vc, double c_sld, double bc, 
    5     double Nd, double Phid, double vd, double d_sld, double bd, 
     2    double N[], double Phi[], double v[], double L[], double b[], 
    63    double Kab, double Kac, double Kad, 
    74    double Kbc, double Kbd, double Kcd 
    85    ); 
    96 
    10 double Iqxy(double qx, double qy, double case_num, 
    11     double Na, double Phia, double va, double a_sld, double ba, 
    12     double Nb, double Phib, double vb, double b_sld, double bb, 
    13     double Nc, double Phic, double vc, double c_sld, double bc, 
    14     double Nd, double Phid, double vd, double d_sld, double bd, 
    15     double Kab, double Kac, double Kad, 
    16     double Kbc, double Kbd, double Kcd 
    17     ); 
    18  
    19 double form_volume(void); 
    20  
    21 double form_volume(void) 
    22 { 
    23     return 1.0; 
    24 } 
    25  
    267double Iq(double q, double case_num, 
    27     double Na, double Phia, double va, double La, double ba, 
    28     double Nb, double Phib, double vb, double Lb, double bb, 
    29     double Nc, double Phic, double vc, double Lc, double bc, 
    30     double Nd, double Phid, double vd, double Ld, double bd, 
     8    double N[], double Phi[], double v[], double L[], double b[], 
    319    double Kab, double Kac, double Kad, 
    3210    double Kbc, double Kbd, double Kcd 
    … …  
    3614#if 0  // Sasview defaults 
    3715  if (icase <= 1) { 
    38     Na=Nb=1000.0; 
    39     Phia=Phib=0.0000001; 
     16    N[0]=N[1]=1000.0; 
     17    Phi[0]=Phi[1]=0.0000001; 
    4018    Kab=Kac=Kad=Kbc=Kbd=-0.0004; 
    41     La=Lb=1e-12; 
    42     va=vb=100.0; 
    43     ba=bb=5.0; 
     19    L[0]=L[1]=1e-12; 
     20    v[0]=v[1]=100.0; 
     21    b[0]=b[1]=5.0; 
    4422  } else if (icase <= 4) { 
    45     Phia=0.0000001; 
     23    Phi[0]=0.0000001; 
    4624    Kab=Kac=Kad=-0.0004; 
    47     La=1e-12; 
    48     va=100.0; 
    49     ba=5.0; 
     25    L[0]=1e-12; 
     26    v[0]=100.0; 
     27    b[0]=5.0; 
    5028  } 
    5129#else 
    5230  if (icase <= 1) { 
    53     Na=Nb=0.0; 
    54     Phia=Phib=0.0; 
     31    N[0]=N[1]=0.0; 
     32    Phi[0]=Phi[1]=0.0; 
    5533    Kab=Kac=Kad=Kbc=Kbd=0.0; 
    56     La=Lb=Ld; 
    57     va=vb=vd; 
    58     ba=bb=0.0; 
     34    L[0]=L[1]=L[3]; 
     35    v[0]=v[1]=v[3]; 
     36    b[0]=b[1]=0.0; 
    5937  } else if (icase <= 4) { 
    60     Na = 0.0; 
    61     Phia=0.0; 
     38    N[0] = 0.0; 
     39    Phi[0]=0.0; 
    6240    Kab=Kac=Kad=0.0; 
    63     La=Ld; 
    64     va=vd; 
    65     ba=0.0; 
     41    L[0]=L[3]; 
     42    v[0]=v[3]; 
     43    b[0]=0.0; 
    6644  } 
    6745#endif 
    6846 
    69   const double Xa = q*q*ba*ba*Na/6.0; 
    70   const double Xb = q*q*bb*bb*Nb/6.0; 
    71   const double Xc = q*q*bc*bc*Nc/6.0; 
    72   const double Xd = q*q*bd*bd*Nd/6.0; 
     47  const double Xa = q*q*b[0]*b[0]*N[0]/6.0; 
     48  const double Xb = q*q*b[1]*b[1]*N[1]/6.0; 
     49  const double Xc = q*q*b[2]*b[2]*N[2]/6.0; 
     50  const double Xd = q*q*b[3]*b[3]*N[3]/6.0; 
    7351 
    7452  // limit as Xa goes to 0 is 1 
    … …  
    9876#if 0 
    9977  const double S0aa = icase<5 
    100                       ? 1.0 : Na*Phia*va*Paa; 
     78                      ? 1.0 : N[0]*Phi[0]*v[0]*Paa; 
    10179  const double S0bb = icase<2 
    102                       ? 1.0 : Nb*Phib*vb*Pbb; 
    103   const double S0cc = Nc*Phic*vc*Pcc; 
    104   const double S0dd = Nd*Phid*vd*Pdd; 
     80                      ? 1.0 : N[1]*Phi[1]*v[1]*Pbb; 
     81  const double S0cc = N[2]*Phi[2]*v[2]*Pcc; 
     82  const double S0dd = N[3]*Phi[3]*v[3]*Pdd; 
    10583  const double S0ab = icase<8 
    106                       ? 0.0 : sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb; 
     84                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 
    10785  const double S0ac = icase<9 
    108                       ? 0.0 : sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb); 
     86                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 
    10987  const double S0ad = icase<9 
    110                       ? 0.0 : sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc); 
     88                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 
    11189  const double S0bc = (icase!=4 && icase!=7 && icase!= 9) 
    112                       ? 0.0 : sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc; 
     90                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 
    11391  const double S0bd = (icase!=4 && icase!=7 && icase!= 9) 
    114                       ? 0.0 : sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc); 
     92                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 
    11593  const double S0cd = (icase==0 || icase==2 || icase==5) 
    116                       ? 0.0 : sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd; 
     94                      ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 
    11795#else  // sasview equivalent 
    118 //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,Nc,Phic,vc,Pcc); 
    119   double S0aa = Na*Phia*va*Paa; 
    120   double S0bb = Nb*Phib*vb*Pbb; 
    121   double S0cc = Nc*Phic*vc*Pcc; 
    122   double S0dd = Nd*Phid*vd*Pdd; 
    123   double S0ab = sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb; 
    124   double S0ac = sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb); 
    125   double S0ad = sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc); 
    126   double S0bc = sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc; 
    127   double S0bd = sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc); 
    128   double S0cd = sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd; 
     96//printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc); 
     97  double S0aa = N[0]*Phi[0]*v[0]*Paa; 
     98  double S0bb = N[1]*Phi[1]*v[1]*Pbb; 
     99  double S0cc = N[2]*Phi[2]*v[2]*Pcc; 
     100  double S0dd = N[3]*Phi[3]*v[3]*Pdd; 
     101  double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 
     102  double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 
     103  double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 
     104  double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 
     105  double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 
     106  double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 
    129107switch(icase){ 
    130108  case 0: 
    … …  
    311289  // Note: 1e-13 to convert from fm to cm for scattering length 
    312290  const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13; 
    313   const double Lad = icase<5 ? 0.0 : (La/va - Ld/vd)*sqrt_Nav; 
    314   const double Lbd = icase<2 ? 0.0 : (Lb/vb - Ld/vd)*sqrt_Nav; 
    315   const double Lcd = (Lc/vc - Ld/vd)*sqrt_Nav; 
     291  const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav; 
     292  const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav; 
     293  const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav; 
    316294 
    317295  const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33 
    … …  
    321299 
    322300} 
    323  
    324 double Iqxy(double qx, double qy, 
    325     double case_num, 
    326     double Na, double Phia, double va, double a_sld, double ba, 
    327     double Nb, double Phib, double vb, double b_sld, double bb, 
    328     double Nc, double Phic, double vc, double c_sld, double bc, 
    329     double Nd, double Phid, double vd, double d_sld, double bd, 
    330     double Kab, double Kac, double Kad, 
    331     double Kbc, double Kbd, double Kcd 
    332     ) 
    333 { 
    334     double q = sqrt(qx*qx + qy*qy); 
    335     return Iq(q, 
    336         case_num, 
    337         Na, Phia, va, a_sld, ba, 
    338         Nb, Phib, vb, b_sld, bb, 
    339         Nc, Phic, vc, c_sld, bc, 
    340         Nd, Phid, vd, d_sld, bd, 
    341         Kab, Kac, Kad, 
    342         Kbc, Kbd, Kcd); 
    343 } 

  • sasmodels/models/rpa.py

    rec45c4f rec45c4f  
    8686#   ["name", "units", default, [lower, upper], "type","description"], 
    8787parameters = [ 
    88     ["case_num", CASES, 0, [0, 10], "", "Component organization"], 
     88    ["case_num", "", 1, CASES, "", "Component organization"], 
    8989 
    90     ["Na", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
    91     ["Phia", "", 0.25, [0, 1], "", "volume fraction"], 
    92     ["va", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
    93     ["La", "fm", 10.0, [-inf, inf], "", "scattering length"], 
    94     ["ba", "Ang", 5.0, [0, inf], "", "segment length"], 
    95  
    96     ["Nb", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
    97     ["Phib", "", 0.25, [0, 1], "", "volume fraction"], 
    98     ["vb", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
    99     ["Lb", "fm", 10.0, [-inf, inf], "", "scattering length"], 
    100     ["bb", "Ang", 5.0, [0, inf], "", "segment length"], 
    101  
    102     ["Nc", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
    103     ["Phic", "", 0.25, [0, 1], "", "volume fraction"], 
    104     ["vc", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
    105     ["Lc", "fm", 10.0, [-inf, inf], "", "scattering length"], 
    106     ["bc", "Ang", 5.0, [0, inf], "", "segment length"], 
    107  
    108     ["Nd", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
    109     ["Phid", "", 0.25, [0, 1], "", "volume fraction"], 
    110     ["vd", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
    111     ["Ld", "fm", 10.0, [-inf, inf], "", "scattering length"], 
    112     ["bd", "Ang", 5.0, [0, inf], "", "segment length"], 
     90    ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
     91    ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"], 
     92    ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
     93    ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"], 
     94    ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"], 
    11395 
    11496    ["Kab", "", -0.0004, [-inf, inf], "", "Interaction parameter"], 

  • sasmodels/models/spherical_sld.py

    rec45c4f rd2bb604  
    170170# pylint: disable=bad-whitespace, line-too-long 
    171171#            ["name", "units", default, [lower, upper], "type", "description"], 
    172 parameters = [["n_shells",         "",               1,      [0, 9],         "", "number of shells"], 
     172parameters = [["n",                "",               1,      [0, 9],         "", "number of shells"], 
    173173              ["radius_core",      "Ang",            50.0,   [0, inf],       "", "intern layer thickness"], 
    174174              ["sld_core",         "1e-6/Ang^2",     2.07,   [-inf, inf],    "", "sld function flat"], 
    … …  
    192192 
    193193demo = dict( 
    194     n_shells=4, 
     194    n=4, 
    195195    scale=1.0, 
    196196    solvent_sld=1.0, 

  • sasmodels/models/squarewell.py

    rec45c4f rd2bb604  
    123123""" 
    124124 
    125 Iqxy = """ 
    126     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    127     """ 
    128  
    129125# ER defaults to 0.0 
    130126# VR defaults to 1.0 

  • sasmodels/models/stickyhardsphere.py

    rec45c4f rd2bb604  
    171171""" 
    172172 
    173 Iqxy = """ 
    174     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    175     """ 
    176  
    177173# ER defaults to 0.0 
    178174# VR defaults to 1.0 

  • sasmodels/product.py

    rf247314 r7ae2b7f  
    1111*ProductModel(P, S)*. 
    1212""" 
    13 import numpy as np 
     13import numpy as np  # type: ignore 
    1414 
    15 from .core import call_ER_VR 
    16 from .generate import process_parameters 
     15from .details import dispersion_mesh 
     16from .modelinfo import suffix_parameter, ParameterTable, ModelInfo 
     17from .kernel import KernelModel, Kernel 
    1718 
    18 SCALE=0 
    19 BACKGROUND=1 
    20 RADIUS_EFFECTIVE=2 
    21 VOLFRACTION=3 
     19try: 
     20    from typing import Tuple 
     21    from .modelinfo import ParameterSet 
     22    from .details import CallDetails 
     23except ImportError: 
     24    pass 
     25 
     26# TODO: make estimates available to constraints 
     27#ESTIMATED_PARAMETERS = [ 
     28#    ["est_effect_radius", "A", 0.0, [0, np.inf], "", "Estimated effective radius"], 
     29#    ["est_volume_ratio", "", 1.0, [0, np.inf], "", "Estimated volume ratio"], 
     30#] 
    2231 
    2332# TODO: core_shell_sphere model has suppressed the volume ratio calculation 
    2433# revert it after making VR and ER available at run time as constraints. 
    2534def make_product_info(p_info, s_info): 
     35    # type: (ModelInfo, ModelInfo) -> ModelInfo 
    2636    """ 
    2737    Create info block for product model. 
    2838    """ 
    29     p_id, p_name, p_pars = p_info['id'], p_info['name'], p_info['parameters'] 
    30     s_id, s_name, s_pars = s_info['id'], s_info['name'], s_info['parameters'] 
    31     # We require models to start with scale and background 
    32     assert s_pars[SCALE].name == 'scale' 
    33     assert s_pars[BACKGROUND].name == 'background' 
    34     # We require structure factors to start with effect radius and volfraction 
    35     assert s_pars[RADIUS_EFFECTIVE].name == 'radius_effective' 
    36     assert s_pars[VOLFRACTION].name == 'volfraction' 
    37     # Combine the parameter sets.  We are skipping the first three 
    38     # parameters of S since scale, background are defined in P and 
    39     # effect_radius is set from P.ER(). 
    40     pars = p_pars + s_pars[3:] 
    41     # check for duplicates; can't use assertion since they may never be checked 
    42     if len(set(p.name for p in pars)) != len(pars): 
    43         raise ValueError("Duplicate parameters in %s and %s"%(p_id)) 
    44     model_info = {} 
    45     model_info['id'] = '*'.join((p_id, s_id)) 
    46     model_info['name'] = ' X '.join((p_name, s_name)) 
    47     model_info['filename'] = None 
    48     model_info['title'] = 'Product of %s and structure factor %s'%(p_name, s_name) 
    49     model_info['description'] = model_info['title'] 
    50     model_info['docs'] = model_info['title'] 
    51     model_info['category'] = "custom" 
    52     model_info['parameters'] = pars 
    53     #model_info['single'] = p_info['single'] and s_info['single'] 
    54     model_info['structure_factor'] = False 
    55     model_info['variant_info'] = None 
    56     #model_info['tests'] = [] 
    57     #model_info['source'] = [] 
     39    p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters 
     40    s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters 
     41    p_set = set(p.id for p in p_pars.call_parameters) 
     42    s_set = set(p.id for p in s_pars.call_parameters) 
     43 
     44    if p_set & s_set: 
     45        # there is some overlap between the parameter names; tag the 
     46        # overlapping S parameters with name_S 
     47        s_list = [(suffix_parameter(par, "_S") if par.id in p_set else par) 
     48                  for par in s_pars.kernel_parameters] 
     49        combined_pars = p_pars.kernel_parameters + s_list 
     50    else: 
     51        combined_pars = p_pars.kernel_parameters + s_pars.kernel_parameters 
     52    parameters = ParameterTable(combined_pars) 
     53 
     54    model_info = ModelInfo() 
     55    model_info.id = '*'.join((p_id, s_id)) 
     56    model_info.name = ' X '.join((p_name, s_name)) 
     57    model_info.filename = None 
     58    model_info.title = 'Product of %s and %s'%(p_name, s_name) 
     59    model_info.description = model_info.title 
     60    model_info.docs = model_info.title 
     61    model_info.category = "custom" 
     62    model_info.parameters = parameters 
     63    #model_info.single = p_info.single and s_info.single 
     64    model_info.structure_factor = False 
     65    model_info.variant_info = None 
     66    #model_info.tests = [] 
     67    #model_info.source = [] 
    5868    # Iq, Iqxy, form_volume, ER, VR and sesans 
    59     model_info['composition'] = ('product', [p_info, s_info]) 
    60     process_parameters(model_info) 
     69    model_info.composition = ('product', [p_info, s_info]) 
    6170    return model_info 
    6271 
    63 class ProductModel(object): 
     72class ProductModel(KernelModel): 
    6473    def __init__(self, model_info, P, S): 
     74        # type: (ModelInfo, KernelModel, KernelModel) -> None 
    6575        self.info = model_info 
    6676        self.P = P 
    … …  
    6878 
    6979    def __call__(self, q_vectors): 
     80        # type: (List[np.ndarray]) -> Kernel 
    7081        # Note: may be sending the q_vectors to the GPU twice even though they 
    7182        # are only needed once.  It would mess up modularity quite a bit to 
    … …  
    7485        # in opencl; or both in opencl, but one in single precision and the 
    7586        # other in double precision). 
    76         p_kernel = self.P(q_vectors) 
    77         s_kernel = self.S(q_vectors) 
     87        p_kernel = self.P.make_kernel(q_vectors) 
     88        s_kernel = self.S.make_kernel(q_vectors) 
    7889        return ProductKernel(self.info, p_kernel, s_kernel) 
    7990 
    8091    def release(self): 
     92        # type: (None) -> None 
    8193        """ 
    8294        Free resources associated with the model. 
    … …  
    8698 
    8799 
    88 class ProductKernel(object): 
     100class ProductKernel(Kernel): 
    89101    def __init__(self, model_info, p_kernel, s_kernel): 
    90         dim = '2d' if p_kernel.q_input.is_2d else '1d' 
    91  
    92         # Need to know if we want 2D and magnetic parameters when constructing 
    93         # a parameter map. 
    94         par_map = {} 
    95         p_info = p_kernel.info['partype'] 
    96         s_info = s_kernel.info['partype'] 
    97         vol_pars = set(p_info['volume']) 
    98         if dim == '2d': 
    99             num_p_fixed = len(p_info['fixed-2d']) 
    100             num_p_pd = len(p_info['pd-2d']) 
    101             num_s_fixed = len(s_info['fixed-2d']) 
    102             num_s_pd = len(s_info['pd-2d']) - 1 # exclude effect_radius 
    103             # volume parameters are amongst the pd pars for P, not S 
    104             vol_par_idx = [k for k,v in enumerate(p_info['pd-2d']) 
    105                            if v in vol_pars] 
    106         else: 
    107             num_p_fixed = len(p_info['fixed-1d']) 
    108             num_p_pd = len(p_info['pd-1d']) 
    109             num_s_fixed = len(s_info['fixed-1d']) 
    110             num_s_pd = len(s_info['pd-1d']) - 1  # exclude effect_radius 
    111             # volume parameters are amongst the pd pars for P, not S 
    112             vol_par_idx = [k for k,v in enumerate(p_info['pd-1d']) 
    113                            if v in vol_pars] 
    114  
    115         start = 0 
    116         par_map['p_fixed'] = np.arange(start, start+num_p_fixed) 
    117         # User doesn't set scale, background or effect_radius for S in P*S, 
    118         # so borrow values from end of p_fixed.  This makes volfraction the 
    119         # first S parameter. 
    120         start += num_p_fixed 
    121         par_map['s_fixed'] = np.hstack(([start,start], 
    122                                         np.arange(start, start+num_s_fixed-2))) 
    123         par_map['volfraction'] = num_p_fixed 
    124         start += num_s_fixed-2 
    125         # vol pars offset from the start of pd pars 
    126         par_map['vol_pars'] = [start+k for k in vol_par_idx] 
    127         par_map['p_pd'] = np.arange(start, start+num_p_pd) 
    128         start += num_p_pd-1 
    129         par_map['s_pd'] = np.hstack((start, 
    130                                      np.arange(start, start+num_s_pd-1))) 
    131  
    132         self.fixed_pars = model_info['partype']['fixed-' + dim] 
    133         self.pd_pars = model_info['partype']['pd-' + dim] 
     102        # type: (ModelInfo, Kernel, Kernel) -> None 
    134103        self.info = model_info 
    135104        self.p_kernel = p_kernel 
    136105        self.s_kernel = s_kernel 
    137         self.par_map = par_map 
    138106 
    139     def __call__(self, fixed_pars, pd_pars, cutoff=1e-5): 
    140         pars = fixed_pars + pd_pars 
    141         scale = pars[SCALE] 
    142         background = pars[BACKGROUND] 
    143         s_volfraction = pars[self.par_map['volfraction']] 
    144         p_fixed = [pars[k] for k in self.par_map['p_fixed']] 
    145         s_fixed = [pars[k] for k in self.par_map['s_fixed']] 
    146         p_pd = [pars[k] for k in self.par_map['p_pd']] 
    147         s_pd = [pars[k] for k in self.par_map['s_pd']] 
    148         vol_pars = [pars[k] for k in self.par_map['vol_pars']] 
    149  
     107    def __call__(self, details, weights, values, cutoff): 
     108        # type: (CallDetails, np.ndarray, np.ndarray, float) -> np.ndarray 
    150109        effect_radius, vol_ratio = call_ER_VR(self.p_kernel.info, vol_pars) 
    151110 
    152         p_fixed[SCALE] = s_volfraction 
    153         p_fixed[BACKGROUND] = 0.0 
    154         s_fixed[SCALE] = scale 
    155         s_fixed[BACKGROUND] = 0.0 
    156         s_fixed[2] = s_volfraction/vol_ratio 
    157         s_pd[0] = [effect_radius], [1.0] 
    158  
    159         p_res = self.p_kernel(p_fixed, p_pd, cutoff) 
    160         s_res = self.s_kernel(s_fixed, s_pd, cutoff) 
     111        p_details, s_details = details.parts 
     112        p_res = self.p_kernel(p_details, weights, values, cutoff) 
     113        s_res = self.s_kernel(s_details, weights, values, cutoff) 
    161114        #print s_fixed, s_pd, p_fixed, p_pd 
    162115 
    163         return p_res*s_res + background 
     116        return values[0]*(p_res*s_res) + values[1] 
    164117 
    165118    def release(self): 
     119        # type: () -> None 
    166120        self.p_kernel.release() 
    167         self.q_kernel.release() 
     121        self.s_kernel.release() 
    168122 
     123def call_ER_VR(model_info, pars): 
     124    """ 
     125    Return effect radius and volume ratio for the model. 
     126    """ 
     127    if model_info.ER is None and model_info.VR is None: 
     128        return 1.0, 1.0 
     129 
     130    value, weight = _vol_pars(model_info, pars) 
     131 
     132    if model_info.ER is not None: 
     133        individual_radii = model_info.ER(*value) 
     134        effect_radius = np.sum(weight*individual_radii) / np.sum(weight) 
     135    else: 
     136        effect_radius = 1.0 
     137 
     138    if model_info.VR is not None: 
     139        whole, part = model_info.VR(*value) 
     140        volume_ratio = np.sum(weight*part)/np.sum(weight*whole) 
     141    else: 
     142        volume_ratio = 1.0 
     143 
     144    return effect_radius, volume_ratio 
     145 
     146def _vol_pars(model_info, pars): 
     147    # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray] 
     148    vol_pars = [get_weights(p, pars) 
     149                for p in model_info.parameters.call_parameters 
     150                if p.type == 'volume'] 
     151    value, weight = dispersion_mesh(model_info, vol_pars) 
     152    return value, weight 
     153 

  • sasmodels/resolution.py

    r4d76711 r7ae2b7f  
    66from __future__ import division 
    77 
    8 from scipy.special import erf 
    9 from numpy import sqrt, log, log10 
    10 import numpy as np 
     8from scipy.special import erf  # type: ignore 
     9from numpy import sqrt, log, log10, exp, pi  # type: ignore 
     10import numpy as np  # type: ignore 
    1111 
    1212__all__ = ["Resolution", "Perfect1D", "Pinhole1D", "Slit1D", 
    … …  
    3535    smeared theory I(q). 
    3636    """ 
    37     q = None 
    38     q_calc = None 
     37    q = None  # type: np.ndarray 
     38    q_calc = None  # type: np.ndarray 
    3939    def apply(self, theory): 
    4040        """ 
    … …  
    476476    *pars* are the parameter values to use when evaluating. 
    477477    """ 
    478     from sasmodels import core 
     478    from sasmodels import direct_model 
    479479    kernel = form.make_kernel([q]) 
    480     theory = core.call_kernel(kernel, pars) 
     480    theory = direct_model.call_kernel(kernel, pars) 
    481481    kernel.release() 
    482482    return theory 
    … …  
    489489    *q0* is the center, *dq* is the width and *q* are the points to evaluate. 
    490490    """ 
    491     from numpy import exp, pi 
    492491    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq) 
    493492 
    … …  
    500499    that make it slow to evaluate but give it good accuracy. 
    501500    """ 
    502     from scipy.integrate import romberg 
    503  
    504     if any(k not in form.info['defaults'] for k in pars.keys()): 
    505         keys = set(form.info['defaults'].keys()) 
    506         extra = set(pars.keys()) - keys 
    507         raise ValueError("bad parameters: [%s] not in [%s]"% 
    508                          (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 
     501    from scipy.integrate import romberg  # type: ignore 
     502 
     503    par_set = set([p.name for p in form.info.parameters.call_parameters]) 
     504    if any(k not in par_set for k in pars.keys()): 
     505        extra = set(pars.keys()) - par_set 
     506        raise ValueError("bad parameters: [%s] not in [%s]" 
     507                         % (", ".join(sorted(extra)), 
     508                            ", ".join(sorted(pars.keys())))) 
    509509 
    510510    if np.isscalar(width): 
    … …  
    554554    that make it slow to evaluate but give it good accuracy. 
    555555    """ 
    556     from scipy.integrate import romberg 
    557  
    558     if any(k not in form.info['defaults'] for k in pars.keys()): 
    559         keys = set(form.info['defaults'].keys()) 
    560         extra = set(pars.keys()) - keys 
    561         raise ValueError("bad parameters: [%s] not in [%s]"% 
    562                          (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 
     556    from scipy.integrate import romberg  # type: ignore 
     557 
     558    par_set = set([p.name for p in form.info.parameters.call_parameters]) 
     559    if any(k not in par_set for k in pars.keys()): 
     560        extra = set(pars.keys()) - par_set 
     561        raise ValueError("bad parameters: [%s] not in [%s]" 
     562                         % (", ".join(sorted(extra)), 
     563                            ", ".join(sorted(pars.keys())))) 
    563564 
    564565    _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq) 
    … …  
    692693 
    693694    def _eval_sphere(self, pars, resolution): 
    694         from sasmodels import core 
     695        from sasmodels import direct_model 
    695696        kernel = self.model.make_kernel([resolution.q_calc]) 
    696         theory = core.call_kernel(kernel, pars) 
     697        theory = direct_model.call_kernel(kernel, pars) 
    697698        result = resolution.apply(theory) 
    698699        kernel.release() 
    … …  
    750751        #tol = 0.001 
    751752        ## The default 3 sigma and no extra points gets 1% 
    752         q_calc, tol = None, 0.01 
     753        q_calc = None  # type: np.ndarray 
     754        tol = 0.01 
    753755        resolution = Pinhole1D(q, q_width, q_calc=q_calc) 
    754756        output = self._eval_sphere(pars, resolution) 
    755757        if 0: # debug plot 
    756             import matplotlib.pyplot as plt 
     758            import matplotlib.pyplot as plt  # type: ignore 
    757759            resolution = Perfect1D(q) 
    758760            source = self._eval_sphere(pars, resolution) 
    … …  
    10261028    """ 
    10271029    import sys 
    1028     import xmlrunner 
     1030    import xmlrunner  # type: ignore 
    10291031 
    10301032    suite = unittest.TestSuite() 
    … …  
    10431045    import sys 
    10441046    from sasmodels import core 
     1047    from sasmodels import direct_model 
    10451048    name = sys.argv[1] if len(sys.argv) > 1 else 'cylinder' 
    10461049 
    … …  
    10631066 
    10641067    kernel = model.make_kernel([resolution.q_calc]) 
    1065     theory = core.call_kernel(kernel, pars) 
     1068    theory = direct_model.call_kernel(kernel, pars) 
    10661069    Iq = resolution.apply(theory) 
    10671070 
    … …  
    10731076        Iq_romb = romberg_pinhole_1d(resolution.q, dq, model, pars) 
    10741077 
    1075     import matplotlib.pyplot as plt 
     1078    import matplotlib.pyplot as plt  # type: ignore 
    10761079    plt.loglog(resolution.q_calc, theory, label='unsmeared') 
    10771080    plt.loglog(resolution.q, Iq, label='smeared', hold=True) 
    … …  
    11081111    Run the resolution demos. 
    11091112    """ 
    1110     import matplotlib.pyplot as plt 
     1113    import matplotlib.pyplot as plt  # type: ignore 
    11111114    plt.subplot(121) 
    11121115    demo_pinhole_1d() 

  • sasmodels/resolution2d.py

    rd6f5da6 r7ae2b7f  
    77from __future__ import division 
    88 
    9 import numpy as np 
    10 from numpy import pi, cos, sin, sqrt 
     9import numpy as np  # type: ignore 
     10from numpy import pi, cos, sin, sqrt  # type: ignore 
    1111 
    1212from . import resolution 

  • sasmodels/sasview_model.py

    r81ec7c8 r81ec7c8  
    2121import logging 
    2222 
    23 import numpy as np 
     23import numpy as np  # type: ignore 
    2424 
    2525from . import core 
    2626from . import custom 
    2727from . import generate 
     28from . import weights 
     29from . import details 
     30from . import modelinfo 
    2831 
    2932try: 
    3033    from typing import Dict, Mapping, Any, Sequence, Tuple, NamedTuple, List, Optional 
     34    from .modelinfo import ModelInfo, Parameter 
    3135    from .kernel import KernelModel 
    3236    MultiplicityInfoType = NamedTuple( 
    … …  
    4448) 
    4549 
     50# TODO: figure out how to say that the return type is a subclass 
    4651def load_standard_models(): 
     52    # type: () -> List[type] 
    4753    """ 
    4854    Load and return the list of predefined models. 
    … …  
    5561        try: 
    5662            models.append(_make_standard_model(name)) 
    57         except: 
     63        except Exception: 
    5864            logging.error(traceback.format_exc()) 
    5965    return models 
    … …  
    6167 
    6268def load_custom_model(path): 
     69    # type: (str) -> type 
    6370    """ 
    6471    Load a custom model given the model path. 
    6572    """ 
    6673    kernel_module = custom.load_custom_kernel_module(path) 
    67     model_info = generate.make_model_info(kernel_module) 
     74    model_info = modelinfo.make_model_info(kernel_module) 
    6875    return _make_model_from_info(model_info) 
    6976 
    7077 
    7178def _make_standard_model(name): 
     79    # type: (str) -> type 
    7280    """ 
    7381    Load the sasview model defined by *name*. 
    … …  
    7886    """ 
    7987    kernel_module = generate.load_kernel_module(name) 
    80     model_info = generate.make_model_info(kernel_module) 
     88    model_info = modelinfo.make_model_info(kernel_module) 
    8189    return _make_model_from_info(model_info) 
    8290 
    8391 
    8492def _make_model_from_info(model_info): 
     93    # type: (ModelInfo) -> type 
    8594    """ 
    8695    Convert *model_info* into a SasView model wrapper. 
    8796    """ 
    88     model_info['variant_info'] = None  # temporary hack for older sasview 
    89     def __init__(self, multiplicity=1): 
     97    def __init__(self, multiplicity=None): 
    9098        SasviewModel.__init__(self, multiplicity=multiplicity) 
    9199    attrs = _generate_model_attributes(model_info) 
    92100    attrs['__init__'] = __init__ 
    93     ConstructedModel = type(model_info['id'], (SasviewModel,), attrs) 
     101    ConstructedModel = type(model_info.name, (SasviewModel,), attrs) 
    94102    return ConstructedModel 
    95103 
    … …  
    106114    attrs = {}  # type: Dict[str, Any] 
    107115    attrs['_model_info'] = model_info 
    108     attrs['name'] = model_info['name'] 
    109     attrs['id'] = model_info['id'] 
    110     attrs['description'] = model_info['description'] 
    111     attrs['category'] = model_info['category'] 
     116    attrs['name'] = model_info.name 
     117    attrs['id'] = model_info.id 
     118    attrs['description'] = model_info.description 
     119    attrs['category'] = None 
    112120 
    113121    # TODO: allow model to override axis labels input/output name/unit 
     122 
     123    parameters = model_info.parameters 
    114124 
    115125    #self.is_multifunc = False 
    116126    non_fittable = []  # type: List[str] 
    117     variants = MultiplicityInfo(0, "", [], "") 
    118     attrs['is_structure_factor'] = model_info['structure_factor'] 
    119     attrs['is_form_factor'] = model_info['ER'] is not None 
     127    xlabel = model_info.profile_axes[0] if model_info.profile is not None else "" 
     128    variants = MultiplicityInfo(0, "", [], xlabel) 
     129    for p in parameters.kernel_parameters: 
     130        if p.name == model_info.control: 
     131            non_fittable.append(p.name) 
     132            variants = MultiplicityInfo( 
     133                len(p.choices), p.name, p.choices, xlabel 
     134            ) 
     135            break 
     136        elif p.is_control: 
     137            non_fittable.append(p.name) 
     138            variants = MultiplicityInfo( 
     139                int(p.limits[1]), p.name, p.choices, xlabel 
     140            ) 
     141            break 
     142 
     143    attrs['is_structure_factor'] = model_info.structure_factor 
     144    attrs['is_form_factor'] = model_info.ER is not None 
    120145    attrs['is_multiplicity_model'] = variants[0] > 1 
    121146    attrs['multiplicity_info'] = variants 
    122147 
    123     partype = model_info['partype'] 
    124     orientation_params = ( 
    125             partype['orientation'] 
    126             + [n + '.width' for n in partype['orientation']] 
    127             + partype['magnetic']) 
    128     magnetic_params = partype['magnetic'] 
    129     fixed = [n + '.width' for n in partype['pd-2d']] 
    130  
     148    orientation_params = [] 
     149    magnetic_params = [] 
     150    fixed = [] 
     151    for p in parameters.user_parameters(): 
     152        if p.type == 'orientation': 
     153            orientation_params.append(p.name) 
     154            orientation_params.append(p.name+".width") 
     155            fixed.append(p.name+".width") 
     156        if p.type == 'magnetic': 
     157            orientation_params.append(p.name) 
     158            magnetic_params.append(p.name) 
     159            fixed.append(p.name+".width") 
    131160    attrs['orientation_params'] = tuple(orientation_params) 
    132161    attrs['magnetic_params'] = tuple(magnetic_params) 
    … …  
    196225    multiplicity = None     # type: Optional[int] 
    197226 
    198     def __init__(self, multiplicity): 
     227    def __init__(self, multiplicity=None): 
    199228        # type: () -> None 
    200229        print("initializing", self.name) 
    … …  
    206235        self._persistency_dict = {} 
    207236 
     237        # TODO: _persistency_dict to persistency_dict throughout sasview 
     238        # TODO: refactor multiplicity to encompass variants 
     239        # TODO: dispersion should be a class 
     240        # TODO: refactor multiplicity info 
     241        # TODO: separate profile view from multiplicity 
     242        # The button label, x and y axis labels and scale need to be under 
     243        # the control of the model, not the fit page.  Maximum flexibility, 
     244        # the fit page would supply the canvas and the profile could plot 
     245        # how it wants, but this assumes matplotlib.  Next level is that 
     246        # we provide some sort of data description including title, labels 
     247        # and lines to plot. 
     248 
     249        # Get the list of hidden parameters given the mulitplicity 
     250        # Don't include multiplicity in the list of parameters 
    208251        self.multiplicity = multiplicity 
     252        if multiplicity is not None: 
     253            hidden = self._model_info.get_hidden_parameters(multiplicity) 
     254            hidden |= set([self.multiplicity_info.control]) 
     255        else: 
     256            hidden = set() 
    209257 
    210258        self.params = collections.OrderedDict() 
    211259        self.dispersion = {} 
    212260        self.details = {} 
    213  
    214         for p in self._model_info['parameters']: 
     261        for p in self._model_info.parameters.user_parameters(): 
     262            if p.name in hidden: 
     263                continue 
    215264            self.params[p.name] = p.default 
    216             self.details[p.name] = [p.units] + p.limits 
    217  
    218         for name in self._model_info['partype']['pd-2d']: 
    219             self.dispersion[name] = { 
    220                 'width': 0, 
    221                 'npts': 35, 
    222                 'nsigmas': 3, 
    223                 'type': 'gaussian', 
    224             } 
     265            self.details[p.id] = [p.units, p.limits[0], p.limits[1]] 
     266            if p.polydisperse: 
     267                self.details[p.id+".width"] = [ 
     268                    "", 0.0, 1.0 if p.relative_pd else np.inf 
     269                ] 
     270                self.dispersion[p.name] = { 
     271                    'width': 0, 
     272                    'npts': 35, 
     273                    'nsigmas': 3, 
     274                    'type': 'gaussian', 
     275                } 
    225276 
    226277    def __get_state__(self): 
     278        # type: () -> Dict[str, Any] 
    227279        state = self.__dict__.copy() 
    228280        state.pop('_model') 
    … …  
    233285 
    234286    def __set_state__(self, state): 
     287        # type: (Dict[str, Any]) -> None 
    235288        self.__dict__ = state 
    236289        self._model = None 
    237290 
    238291    def __str__(self): 
     292        # type: () -> str 
    239293        """ 
    240294        :return: string representation 
    … …  
    243297 
    244298    def is_fittable(self, par_name): 
     299        # type: (str) -> bool 
    245300        """ 
    246301        Check if a given parameter is fittable or not 
    … …  
    253308 
    254309 
    255     # pylint: disable=no-self-use 
    256310    def getProfile(self): 
     311        # type: () -> (np.ndarray, np.ndarray) 
    257312        """ 
    258313        Get SLD profile 
    … …  
    261316                beta is a list of the corresponding SLD values 
    262317        """ 
    263         return None, None 
     318        args = [] 
     319        for p in self._model_info.parameters.kernel_parameters: 
     320            if p.id == self.multiplicity_info.control: 
     321                args.append(self.multiplicity) 
     322            elif p.length == 1: 
     323                args.append(self.params.get(p.id, np.NaN)) 
     324            else: 
     325                args.append([self.params.get(p.id+str(k), np.NaN) 
     326                             for k in range(1,p.length+1)]) 
     327        return self._model_info.profile(*args) 
    264328 
    265329    def setParam(self, name, value): 
     330        # type: (str, float) -> None 
    266331        """ 
    267332        Set the value of a model parameter 
    … …  
    290355 
    291356    def getParam(self, name): 
     357        # type: (str) -> float 
    292358        """ 
    293359        Set the value of a model parameter 
    … …  
    313379 
    314380    def getParamList(self): 
     381        # type: () - > Sequence[str] 
    315382        """ 
    316383        Return a list of all available parameters for the model 
    … …  
    322389 
    323390    def getDispParamList(self): 
     391        # type: () - > Sequence[str] 
    324392        """ 
    325393        Return a list of polydispersity parameters for the model 
    326394        """ 
    327395        # TODO: fix test so that parameter order doesn't matter 
    328         ret = ['%s.%s' % (d.lower(), p) 
    329                for d in self._model_info['partype']['pd-2d'] 
    330                for p in ('npts', 'nsigmas', 'width')] 
     396        ret = ['%s.%s' % (p.name.lower(), ext) 
     397               for p in self._model_info.parameters.user_parameters() 
     398               for ext in ('npts', 'nsigmas', 'width') 
     399               if p.polydisperse] 
    331400        #print(ret) 
    332401        return ret 
    333402 
    334403    def clone(self): 
     404        # type: () - > "SasviewModel" 
    335405        """ Return a identical copy of self """ 
    336406        return deepcopy(self) 
    337407 
    338408    def run(self, x=0.0): 
     409        # type: (Union[float, (float, float), List[float]]) -> float 
    339410        """ 
    340411        Evaluate the model 
    … …  
    356427 
    357428    def runXY(self, x=0.0): 
     429        # type: (Union[float, (float, float), List[float]]) -> float 
    358430        """ 
    359431        Evaluate the model in cartesian coordinates 
    … …  
    371443 
    372444    def evalDistribution(self, qdist): 
     445        # type: (Union[np.ndarray, Tuple[np.ndarray, np.ndarray], List[np.ndarray]) -> np.ndarray 
    373446        r""" 
    374447        Evaluate a distribution of q-values. 
    … …  
    401474            # Check whether we have a list of ndarrays [qx,qy] 
    402475            qx, qy = qdist 
    403             partype = self._model_info['partype'] 
    404             if not partype['orientation'] and not partype['magnetic']: 
     476            if not self._model_info.parameters.has_2d: 
    405477                return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2)) 
    406478            else: 
    … …  
    415487                            % type(qdist)) 
    416488 
    417     def calculate_Iq(self, *args): 
     489    def calculate_Iq(self, qx, qy=None): 
     490        # type: (Sequence[float], Optional[Sequence[float]]) -> np.ndarray 
    418491        """ 
    419492        Calculate Iq for one set of q with the current parameters. 
    … …  
    426499        if self._model is None: 
    427500            self._model = core.build_model(self._model_info) 
    428         q_vectors = [np.asarray(q) for q in args] 
    429         fn = self._model.make_kernel(q_vectors) 
    430         pars = [self.params[v] for v in fn.fixed_pars] 
    431         pd_pars = [self._get_weights(p) for p in fn.pd_pars] 
    432         result = fn(pars, pd_pars, self.cutoff) 
    433         fn.q_input.release() 
    434         fn.release() 
     501        if qy is not None: 
     502            q_vectors = [np.asarray(qx), np.asarray(qy)] 
     503        else: 
     504            q_vectors = [np.asarray(qx)] 
     505        kernel = self._model.make_kernel(q_vectors) 
     506        pairs = [self._get_weights(p) 
     507                 for p in self._model_info.parameters.call_parameters] 
     508        call_details, weight, value = details.build_details(kernel, pairs) 
     509        result = kernel(call_details, weight, value, cutoff=self.cutoff) 
     510        kernel.release() 
    435511        return result 
    436512 
    437513    def calculate_ER(self): 
     514        # type: () -> float 
    438515        """ 
    439516        Calculate the effective radius for P(q)*S(q) 
    … …  
    441518        :return: the value of the effective radius 
    442519        """ 
    443         ER = self._model_info.get('ER', None) 
    444         if ER is None: 
     520        if self._model_info.ER is None: 
    445521            return 1.0 
    446522        else: 
    447             values, weights = self._dispersion_mesh() 
    448             fv = ER(*values) 
     523            value, weight = self._dispersion_mesh() 
     524            fv = self._model_info.ER(*value) 
    449525            #print(values[0].shape, weights.shape, fv.shape) 
    450             return np.sum(weights * fv) / np.sum(weights) 
     526            return np.sum(weight * fv) / np.sum(weight) 
    451527 
    452528    def calculate_VR(self): 
     529        # type: () -> float 
    453530        """ 
    454531        Calculate the volf ratio for P(q)*S(q) 
    … …  
    456533        :return: the value of the volf ratio 
    457534        """ 
    458         VR = self._model_info.get('VR', None) 
    459         if VR is None: 
     535        if self._model_info.VR is None: 
    460536            return 1.0 
    461537        else: 
    462             values, weights = self._dispersion_mesh() 
    463             whole, part = VR(*values) 
    464             return np.sum(weights * part) / np.sum(weights * whole) 
     538            value, weight = self._dispersion_mesh() 
     539            whole, part = self._model_info.VR(*value) 
     540            return np.sum(weight * part) / np.sum(weight * whole) 
    465541 
    466542    def set_dispersion(self, parameter, dispersion): 
     543        # type: (str, weights.Dispersion) -> Dict[str, Any] 
    467544        """ 
    468545        Set the dispersion object for a model parameter 
    … …  
    493570 
    494571    def _dispersion_mesh(self): 
     572        # type: () -> List[Tuple[np.ndarray, np.ndarray]] 
    495573        """ 
    496574        Create a mesh grid of dispersion parameters and weights. 
    … …  
    500578        parameter set in the vector. 
    501579        """ 
    502         pars = self._model_info['partype']['volume'] 
    503         return core.dispersion_mesh([self._get_weights(p) for p in pars]) 
     580        pars = [self._get_weights(p) 
     581                for p in self._model_info.parameters.call_parameters 
     582                if p.type == 'volume'] 
     583        return details.dispersion_mesh(self._model_info, pars) 
    504584 
    505585    def _get_weights(self, par): 
     586        # type: (Parameter) -> Tuple[np.ndarray, np.ndarray] 
    506587        """ 
    507588        Return dispersion weights for parameter 
    508589        """ 
    509         from . import weights 
    510         relative = self._model_info['partype']['pd-rel'] 
    511         limits = self._model_info['limits'] 
    512         dis = self.dispersion[par] 
    513         value, weight = weights.get_weights( 
    514             dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 
    515             self.params[par], limits[par], par in relative) 
    516         return value, weight / np.sum(weight) 
    517  
     590        if par.name not in self.params: 
     591            if par.name == self.multiplicity_info.control: 
     592                return [self.multiplicity], [] 
     593            else: 
     594                return [np.NaN], [] 
     595        elif par.polydisperse: 
     596            dis = self.dispersion[par.name] 
     597            value, weight = weights.get_weights( 
     598                dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 
     599                self.params[par.name], par.limits, par.relative_pd) 
     600            return value, weight / np.sum(weight) 
     601        else: 
     602            return [self.params[par.name]], [] 
    518603 
    519604def test_model(): 
     605    # type: () -> float 
    520606    """ 
    521607    Test that a sasview model (cylinder) can be run. 
    … …  
    525611    return cylinder.evalDistribution([0.1,0.1]) 
    526612 
     613def test_rpa(): 
     614    # type: () -> float 
     615    """ 
     616    Test that a sasview model (cylinder) can be run. 
     617    """ 
     618    RPA = _make_standard_model('rpa') 
     619    rpa = RPA(3) 
     620    return rpa.evalDistribution([0.1,0.1]) 
     621 
    527622 
    528623def test_model_list(): 
     624    # type: () -> None 
    529625    """ 
    530626    Make sure that all models build as sasview models. 

  • sasmodels/sesans.py

    r2684d45 r7ae2b7f  
    1212from __future__ import division 
    1313 
    14 import numpy as np 
    15 from numpy import pi, exp 
     14import numpy as np  # type: ignore 
     15from numpy import pi, exp  # type: ignore 
    1616from scipy.special import jv as besselj 
    1717#import direct_model.DataMixin as model 

  • sasmodels/weights.py

    ra936688 ra936688  
    33""" 
    44# TODO: include dispersion docs with the disperser models 
    5 from math import sqrt 
    6 import numpy as np 
    7 from scipy.special import gammaln 
     5from math import sqrt  # type: ignore 
     6import numpy as np  # type: ignore 
     7from scipy.special import gammaln  # type: ignore 
    88 
    99class Dispersion(object): 
Note: See TracChangeset for help on using the changeset viewer.