Changeset 2afc26d in sasmodels


Ignore:
Timestamp:
Apr 8, 2016 12:15:23 PM (9 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
a86f6c0
Parents:
4bfd277 (diff), 416609b (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge remote-tracking branch 'origin/master' into polydisp

Files:
6 added
44 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/elliptical_cylinder.py

    rec45c4f r416609b  
    8181""" 
    8282 
    83 import math 
    84 from numpy import pi, inf 
     83from numpy import pi, inf, sqrt 
    8584 
    8685name = "elliptical_cylinder" 
    … …  
    117116        @param length: Length of the cylinder 
    118117    """ 
    119     radius = math.sqrt(r_minor * r_minor * axis_ratio) 
     118    radius = sqrt(r_minor * r_minor * axis_ratio) 
    120119    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
    121120    return 0.5 * (ddd) ** (1. / 3.) 

  • sasmodels/models/line.py

    rec45c4f r416609b  
    99.. math:: 
    1010 
    11     I(q) = A + B \cdot q 
     11    I(q) = \text{scale} (A + B \cdot q) + \text{background} 
    1212 
    1313.. note:: 
    … …  
    1515 
    1616.. math:: 
    17     I(q) = I(qx) \cdot I(qy) 
     17    I(q) = \text{scale} (I(qx) \cdot I(qy)) + \text{background} 
    1818 
    1919References 
    … …  
    5151    """ 
    5252    inten = intercept + slope*q 
    53     # TODO: In SasView code additional formula for list has been specified. 
    54     # if inten(x) = intercept + slope*x: 
    55     # then if q is a list, Iq=inten(x[0]*math.cos(x[1]))*inten(x[0]*math.sin(x[1])) 
    5653    return inten 
    5754 
    … …  
    6562    :return:     2D-Intensity 
    6663    """ 
    67     # TODO: SasView documention lists 2D intensity as Iq(qx)*Iq(qy) but code says: 
    68     # return self._line(x[1]) 
     64    # TODO: SasView documents 2D intensity as Iq(qx)*Iq(qy), but returns Iq(qy) 
     65    # Note: SasView.run([r, theta]) does return Iq(qx)*Iq(qy) 
    6966    return Iq(qx, *args)*Iq(qy, *args) 
    7067 
    7168Iqxy.vectorized = True  # Iqxy accepts an array of qx, qy values 
    7269 
    73 demo = dict(scale=1.0, background=0, intercept=1.0, slope=1.0) 
    74  
    7570tests = [ 
    76  
    77     [{'intercept':   1.0, 
    78       'slope': 1.0, 
    79      }, 1.0, 2.001], 
    80  
    81     [{'intercept':   1.0, 
    82       'slope': 1.0, 
    83      }, 0.0, 1.001], 
    84  
    85     [{'intercept':   1.0, 
    86       'slope': 1.0, 
    87      }, 0.4, 1.401], 
    88  
    89     [{'intercept':   1.0, 
    90       'slope': 1.0, 
    91      }, 1.3, 2.301], 
    92  
    93     [{'intercept':   1.0, 
    94       'slope': 1.0, 
    95      }, 0.5, 1.501], 
    96  
    97     [{'intercept':   1.0, 
    98       'slope': 1.0, 
    99      }, [0.4, 0.5], [1.401, 1.501]], 
    100  
    101     [{'intercept':   1.0, 
    102       'slope': 1.0, 
    103       'background': 0.0, 
    104      }, (1.3, 1.57), 5.911], 
     71    [{'intercept': 1.0, 'slope': 1.0, }, 1.0, 2.001], 
     72    [{'intercept': 1.0, 'slope': 1.0, }, 0.0, 1.001], 
     73    [{'intercept': 1.0, 'slope': 1.0, }, 0.4, 1.401], 
     74    [{'intercept': 1.0, 'slope': 1.0, }, 1.3, 2.301], 
     75    [{'intercept': 1.0, 'slope': 1.0, }, 0.5, 1.501], 
     76    [{'intercept': 1.0, 'slope': 1.0, }, [0.4, 0.5], [1.401, 1.501]], 
     77    [{'intercept': 1.0, 'slope': 1.0, 'background': 0.0, }, (1.3, 1.57), 5.911], 
    10578] 

  • sasmodels/models/onion.py

    rea05c87 r2afc26d  
    382382    # "A[1]": 0, "A[2]": -1, "A[3]": 1e-4, "A[4]": 1, 
    383383    } 
     384 

  • sasmodels/models/polymer_excl_volume.py

    rec45c4f r416609b  
    9393""" 
    9494 
    95 from math import sqrt 
    96 from numpy import inf, power 
     95from numpy import inf, power, sqrt 
    9796from scipy.special import gammainc, gamma 
    9897 
    … …  
    114113 
    115114 
    116 def Iq(q, 
    117        rg=60.0, 
    118        porod_exp=3.0): 
     115def Iq(q, rg=60.0, porod_exp=3.0): 
    119116    """ 
    120117    :param q:         Input q-value (float or [float, float]) 
    … …  
    144141    :return:     2D-Intensity 
    145142    """ 
    146  
    147143    return Iq(sqrt(qx**2 + qy**2), *args) 
    148144 
    149145Iqxy.vectorized = True  # Iqxy accepts an array of qx, qy values 
    150146 
    151  
    152 demo = dict(scale=1, background=0.0, 
    153             rg=60.0, 
    154             porod_exp=3.0) 
    155147 
    156148tests = [ 

  • sasmodels/models/two_lorentzian.py

    rec45c4f r416609b  
    3434""" 
    3535 
    36 from math import sqrt 
    37 from numpy import inf, power 
     36from numpy import inf, power, sqrt 
    3837 
    3938name = "two_lorentzian" 

  • doc/developer/index.rst

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

  • sasmodels/bumps_model.py

    rea75043 r6d6508e  
    8181    from bumps.names import Parameter 
    8282 
    83     pars = {} 
    84     for p in model_info['parameters']: 
     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 r6d6508e  
    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 
    … …  
    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']: 
    … …  
    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 r8bd7b77  
    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 r6d6508e  
    153153def revert_name(model_info): 
    154154    _read_conversion_table() 
    155     oldname, oldpars = CONVERSION_TABLE.get(model_info['id'], [None, {}]) 
     155    oldname, oldpars = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
    156156    return oldname 
    157157 
    158158def _get_old_pars(model_info): 
    159159    _read_conversion_table() 
    160     oldname, oldpars = CONVERSION_TABLE.get(model_info['id'], [None, {}]) 
     160    oldname, oldpars = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
    161161    return oldpars 
    162162 
    … …  
    165165    Convert model from new style parameter names to old style. 
    166166    """ 
    167     if model_info['composition'] is not None: 
    168         composition_type, parts = model_info['composition'] 
     167    if model_info.composition is not None: 
     168        composition_type, parts = model_info.composition 
    169169        if composition_type == 'product': 
    170170            P, S = parts 
    … …  
    180180 
    181181    # Note: update compare.constrain_pars to match 
    182     name = model_info['id'] 
    183     if name in MODELS_WITHOUT_SCALE or model_info['structure_factor']: 
     182    name = model_info.id 
     183    if name in MODELS_WITHOUT_SCALE or model_info.structure_factor: 
    184184        if oldpars.pop('scale', 1.0) != 1.0: 
    185185            warnings.warn("parameter scale not used in sasview %s"%name) 
    186     if name in MODELS_WITHOUT_BACKGROUND or model_info['structure_factor']: 
     186    if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor: 
    187187        if oldpars.pop('background', 0.0) != 0.0: 
    188188            warnings.warn("parameter background not used in sasview %s"%name) 
    … …  
    206206        elif name == 'rpa': 
    207207            # convert scattering lengths from femtometers to centimeters 
    208             for p in "La", "Lb", "Lc", "Ld": 
     208            for p in "L1", "L2", "L3", "L4": 
    209209                if p in oldpars: oldpars[p] *= 1e-13 
    210210        elif name == 'core_shell_parallelepiped': 
    … …  
    225225    Restrict parameter values to those that will match sasview. 
    226226    """ 
    227     name = model_info['id'] 
     227    name = model_info.id 
    228228    # Note: update convert.revert_model to match 
    229     if name in MODELS_WITHOUT_SCALE or model_info['structure_factor']: 
     229    if name in MODELS_WITHOUT_SCALE or model_info.structure_factor: 
    230230        pars['scale'] = 1 
    231     if name in MODELS_WITHOUT_BACKGROUND or model_info['structure_factor']: 
     231    if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor: 
    232232        pars['background'] = 0 
    233233    # sasview multiplies background by structure factor 

  • sasmodels/core.py

    r4d76711 r6d6508e  
    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 
    1214import numpy as np 
    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 
    28  
    29 try: 
    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] 
    3927 
    4028# TODO: refactor composite model support 
    … …  
    6452    """ 
    6553    try: s + '' 
    66     except: return False 
     54    except Exception: return False 
    6755    return True 
    6856 
    … …  
    8876    parts = model_name.split('*') 
    8977    if len(parts) > 1: 
    90         from . import product 
    91         # Note: currently have circular reference 
    9278        if len(parts) > 2: 
    9379            raise ValueError("use P*S to apply structure factor S to model P") 
    … …  
    9682 
    9783    kernel_module = generate.load_kernel_module(model_name) 
    98     return generate.make_model_info(kernel_module) 
     84    return modelinfo.make_model_info(kernel_module) 
    9985 
    10086 
    … …  
    118104    otherwise it uses the default "ocl". 
    119105    """ 
    120     composition = model_info.get('composition', None) 
     106    composition = model_info.composition 
    121107    if composition is not None: 
    122108        composition_type, parts = composition 
    … …  
    137123    ##  4. rerun "python -m sasmodels.direct_model $MODELNAME" 
    138124    ##  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() 
     125    # open(model_info.name+'.c','w').write(source) 
     126    # source = open(model_info.name+'.cl','r').read() 
    141127    source = generate.make_source(model_info) 
    142128    if dtype is None: 
    143         dtype = 'single' if model_info['single'] else 'double' 
    144     if callable(model_info.get('Iq', None)): 
     129        dtype = 'single' if model_info.single else 'double' 
     130    if callable(model_info.Iq): 
    145131        return kernelpy.PyModel(model_info) 
    146132    if (platform == "dll" 
    … …  
    168154    source = generate.make_source(model_info) 
    169155    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/data.py

    rd6f5da6 ree8f734  
    358358        try: 
    359359            return fn(*args, **kw) 
    360         except KeyboardInterrupt: 
    361             raise 
    362         except: 
     360        except Exception: 
    363361            traceback.print_exc() 
    364362 

  • sasmodels/direct_model.py

    r4d76711 r6d6508e  
    2525import numpy as np 
    2626 
    27 from .core import call_kernel, call_ER_VR 
    2827from . import sesans 
     28from . import weights 
    2929from . import resolution 
    3030from . import resolution2d 
     31from . import details 
     32 
     33 
     34def call_kernel(kernel, pars, cutoff=0, mono=False): 
     35    """ 
     36    Call *kernel* returned from *model.make_kernel* with parameters *pars*. 
     37 
     38    *cutoff* is the limiting value for the product of dispersion weights used 
     39    to perform the multidimensional dispersion calculation more quickly at a 
     40    slight cost to accuracy. The default value of *cutoff=0* integrates over 
     41    the entire dispersion cube.  Using *cutoff=1e-5* can be 50% faster, but 
     42    with an error of about 1%, which is usually less than the measurement 
     43    uncertainty. 
     44 
     45    *mono* is True if polydispersity should be set to none on all parameters. 
     46    """ 
     47    parameters = kernel.info.parameters 
     48    if mono: 
     49        active = lambda name: False 
     50    elif kernel.dim == '1d': 
     51        active = lambda name: name in parameters.pd_1d 
     52    elif kernel.dim == '2d': 
     53        active = lambda name: name in parameters.pd_2d 
     54    else: 
     55        active = lambda name: True 
     56 
     57    vw_pairs = [(get_weights(p, pars) if active(p.name) 
     58                 else ([pars.get(p.name, p.default)], [])) 
     59                for p in parameters.call_parameters] 
     60 
     61    call_details, weights, values = details.build_details(kernel, vw_pairs) 
     62    return kernel(call_details, weights, values, cutoff) 
     63 
     64def get_weights(parameter, values): 
     65    """ 
     66    Generate the distribution for parameter *name* given the parameter values 
     67    in *pars*. 
     68 
     69    Uses "name", "name_pd", "name_pd_type", "name_pd_n", "name_pd_sigma" 
     70    from the *pars* dictionary for parameter value and parameter dispersion. 
     71    """ 
     72    value = float(values.get(parameter.name, parameter.default)) 
     73    relative = parameter.relative_pd 
     74    limits = parameter.limits 
     75    disperser = values.get(parameter.name+'_pd_type', 'gaussian') 
     76    npts = values.get(parameter.name+'_pd_n', 0) 
     77    width = values.get(parameter.name+'_pd', 0.0) 
     78    nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 
     79    if npts == 0 or width == 0: 
     80        return [value], [] 
     81    value, weight = weights.get_weights( 
     82        disperser, npts, width, nsigma, value, limits, relative) 
     83    return value, weight / np.sum(weight) 
    3184 
    3285class DataMixin(object): 
    … …  
    67120            self.data_type = 'Iq' 
    68121 
    69         partype = model.info['partype'] 
    70  
    71122        if self.data_type == 'sesans': 
    72123             
    … …  
    82133            q_mono = sesans.make_all_q(data) 
    83134        elif self.data_type == 'Iqxy': 
    84             if not partype['orientation'] and not partype['magnetic']: 
    85                 raise ValueError("not 2D without orientation or magnetic parameters") 
     135            #if not model.info.parameters.has_2d: 
     136            #    raise ValueError("not 2D without orientation or magnetic parameters") 
    86137            q = np.sqrt(data.qx_data**2 + data.qy_data**2) 
    87138            qmin = getattr(data, 'qmin', 1e-16) 
    … …  
    172223    def _calc_theory(self, pars, cutoff=0.0): 
    173224        if self._kernel is None: 
    174             self._kernel = self._model.make_kernel(self._kernel_inputs)  # pylint: disable=attribute-dedata_type 
     225            self._kernel = self._model.make_kernel(self._kernel_inputs) 
    175226            self._kernel_mono = (self._model.make_kernel(self._kernel_mono_inputs) 
    176227                                 if self._kernel_mono_inputs else None) 
    … …  
    214265        return self._calc_theory(pars, cutoff=self.cutoff) 
    215266 
    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  
    222267    def simulate_data(self, noise=None, **pars): 
    223268        """ 
    … …  
    243288        try: 
    244289            values = [float(v) for v in call.split(',')] 
    245         except: 
     290        except Exception: 
    246291            values = [] 
    247292        if len(values) == 1: 

  • sasmodels/generate.py

    re6408d0 r6d6508e  
    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 
    217161 
    218162import numpy as np 
    219163 
     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') 
     167TEMPLATE_ROOT = dirname(__file__) 
    226168 
    227169F16 = np.dtype('float16') 
    … …  
    232174except TypeError: 
    233175    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     ] 
    240176 
    241177# Conversion from units defined in the parameter table for each model 
    … …  
    331267    raise ValueError("%r not found in %s" % (filename, search_path)) 
    332268 
     269 
    333270def model_sources(model_info): 
    334271    """ 
    335272    Return a list of the sources file paths for the module. 
    336273    """ 
    337     search_path = [dirname(model_info['filename']), 
     274    search_path = [dirname(model_info.filename), 
    338275                   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 """ 
     276    return [_search(search_path, f) for f in model_info.source] 
     277 
     278def timestamp(model_info): 
     279    """ 
     280    Return a timestamp for the model corresponding to the most recently 
     281    changed file or dependency. 
     282    """ 
     283    source_files = (model_sources(model_info) 
     284                    + model_templates() 
     285                    + [model_info.filename]) 
     286    newest = max(getmtime(f) for f in source_files) 
     287    return newest 
    354288 
    355289def convert_type(source, dtype): 
    … …  
    362296    if dtype == F16: 
    363297        fbytes = 2 
    364         source = _F16_PRAGMA + _convert_type(source, "half", "f") 
     298        source = _convert_type(source, "half", "f") 
    365299    elif dtype == F32: 
    366300        fbytes = 4 
    … …  
    368302    elif dtype == F64: 
    369303        fbytes = 8 
    370         source = _F64_PRAGMA + source  # Source is already double 
     304        # no need to convert if it is already double 
    371305    elif dtype == F128: 
    372306        fbytes = 16 
    … …  
    399333    Name of the exported kernel symbol. 
    400334    """ 
    401     return model_info['name'] + "_" + ("Iqxy" if is_2d else "Iq") 
     335    return model_info.name + "_" + ("Iqxy" if is_2d else "Iq") 
    402336 
    403337 
    … …  
    411345 
    412346 
    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];\ 
     347_template_cache = {} 
     348def load_template(filename): 
     349    path = joinpath(TEMPLATE_ROOT, filename) 
     350    mtime = getmtime(path) 
     351    if filename not in _template_cache or mtime > _template_cache[filename][0]: 
     352        with open(path) as fid: 
     353            _template_cache[filename] = (mtime, fid.read(), path) 
     354    return _template_cache[filename][1] 
     355 
     356def model_templates(): 
     357    # TODO: fails DRY; templates are listed in two places. 
     358    # should instead have model_info contain a list of paths 
     359    return [joinpath(TEMPLATE_ROOT, filename) 
     360            for filename in ('kernel_header.c', 'kernel_iq.c')] 
     361 
     362 
     363_FN_TEMPLATE = """\ 
     364double %(name)s(%(pars)s); 
     365double %(name)s(%(pars)s) { 
     366    %(body)s 
     367} 
     368 
     369 
    417370""" 
    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 
     371 
     372def _gen_fn(name, pars, body): 
     373    """ 
     374    Generate a function given pars and body. 
     375 
     376    Returns the following string:: 
     377 
     378         double fn(double a, double b, ...); 
     379         double fn(double a, double b, ...) { 
     380             .... 
     381         } 
     382    """ 
     383    par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void' 
     384    return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl} 
     385 
     386def _call_pars(prefix, pars): 
     387    """ 
     388    Return a list of *prefix.parameter* from parameter items. 
     389    """ 
     390    return [p.as_call_reference(prefix) for p in pars] 
     391 
     392_IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)", 
     393                           flags=re.MULTILINE) 
     394def _have_Iqxy(sources): 
     395    """ 
     396    Return true if any file defines Iqxy. 
     397 
     398    Note this is not a C parser, and so can be easily confused by 
     399    non-standard syntax.  Also, it will incorrectly identify the following 
     400    as having Iqxy:: 
     401 
     402        /* 
     403        double Iqxy(qx, qy, ...) { ... fill this in later ... } 
     404        */ 
     405 
     406    If you want to comment out an Iqxy function, use // on the front of the 
     407    line instead. 
     408    """ 
     409    for code in sources: 
     410        if _IQXY_PATTERN.search(code): 
     411            return True 
     412    else: 
     413        return False 
     414 
    437415def make_source(model_info): 
    438416    """ 
    … …  
    441419    Uses source files found in the given search path. 
    442420    """ 
    443     if callable(model_info['Iq']): 
     421    if callable(model_info.Iq): 
    444422        return None 
    445423 
    … …  
    453431    # for computing volume even if we allow non-disperse volume parameters. 
    454432 
    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  
     433    partable = model_info.parameters 
     434 
     435    # Identify parameters for Iq, Iqxy, Iq_magnetic and form_volume. 
     436    # Note that scale and volume are not possible types. 
     437 
     438    # Load templates and user code 
     439    kernel_header = load_template('kernel_header.c') 
     440    kernel_code = load_template('kernel_iq.c') 
     441    user_code = [open(f).read() for f in model_sources(model_info)] 
     442 
     443    # Build initial sources 
     444    source = [kernel_header] + user_code 
     445 
     446    # Make parameters for q, qx, qy so that we can use them in declarations 
     447    q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 
     448    # Generate form_volume function, etc. from body only 
     449    if model_info.form_volume is not None: 
     450        pars = partable.form_volume_parameters 
     451        source.append(_gen_fn('form_volume', pars, model_info.form_volume)) 
     452    if model_info.Iq is not None: 
     453        pars = [q] + partable.iq_parameters 
     454        source.append(_gen_fn('Iq', pars, model_info.Iq)) 
     455    if model_info.Iqxy is not None: 
     456        pars = [qx, qy] + partable.iqxy_parameters 
     457        source.append(_gen_fn('Iqxy', pars, model_info.Iqxy)) 
     458 
     459    # Define the parameter table 
     460    source.append("#define PARAMETER_TABLE \\") 
     461    source.append("\\\n".join(p.as_definition() 
     462                              for p in partable.kernel_parameters)) 
     463 
     464    # Define the function calls 
     465    if partable.form_volume_parameters: 
     466        refs = _call_pars("_v.", partable.form_volume_parameters) 
     467        call_volume = "#define CALL_VOLUME(_v) form_volume(%s)" % (",".join(refs)) 
     468    else: 
     469        # Model doesn't have volume.  We could make the kernel run a little 
     470        # faster by not using/transferring the volume normalizations, but 
     471        # the ifdef's reduce readability more than is worthwhile. 
     472        call_volume = "#define CALL_VOLUME(v) 1.0" 
     473    source.append(call_volume) 
     474 
     475    refs = ["_q[_i]"] + _call_pars("_v.", partable.iq_parameters) 
     476    call_iq = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(refs)) 
     477    if _have_Iqxy(user_code): 
     478        # Call 2D model 
     479        refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("_v.", partable.iqxy_parameters) 
     480        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iqxy(%s)" % (",".join(refs)) 
     481    else: 
     482        # Call 1D model with sqrt(qx^2 + qy^2) 
     483        warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 
     484        # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 
     485        pars_sqrt = ["sqrt(_q[2*_i]*_q[2*_i]+_q[2*_i+1]*_q[2*_i+1])"] + refs[1:] 
     486        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt)) 
     487 
     488    # Fill in definitions for numbers of parameters 
     489    source.append("#define MAX_PD %s"%partable.max_pd) 
     490    source.append("#define NPARS %d"%partable.npars) 
     491 
     492    # TODO: allow mixed python/opencl kernels? 
     493 
     494    # define the Iq kernel 
     495    source.append("#define KERNEL_NAME %s_Iq"%model_info.name) 
     496    source.append(call_iq) 
     497    source.append(kernel_code) 
     498    source.append("#undef CALL_IQ") 
     499    source.append("#undef KERNEL_NAME") 
     500 
     501    # define the Iqxy kernel from the same source with different #defines 
     502    source.append("#define KERNEL_NAME %s_Iqxy"%model_info.name) 
     503    source.append(call_iqxy) 
     504    source.append(kernel_code) 
     505    source.append("#undef CALL_IQ") 
     506    source.append("#undef KERNEL_NAME") 
     507 
     508    return '\n'.join(source) 
    648509 
    649510def load_kernel_module(model_name): 
    … …  
    657518 
    658519 
    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     * *variant_info* contains the information required to select between 
    682       model variants (e.g., the list of cases) or is None if there are no 
    683       model variants 
    684     * *defaults* is the *{parameter: value}* table built from the parameter 
    685       description table. 
    686     * *limits* is the *{parameter: [min, max]}* table built from the 
    687       parameter description table. 
    688     * *partypes* categorizes the model parameters. See 
    689       :func:`categorize_parameters` for details. 
    690     * *demo* contains the *{parameter: value}* map used in compare (and maybe 
    691       for the demo plot, if plots aren't set up to use the default values). 
    692       If *demo* is not given in the file, then the default values will be used. 
    693     * *tests* is a set of tests that must pass 
    694     * *source* is the list of library files to include in the C model build 
    695     * *Iq*, *Iqxy*, *form_volume*, *ER*, *VR* and *sesans* are python functions 
    696       implementing the kernel for the module, or None if they are not 
    697       defined in python 
    698     * *composition* is None if the model is independent, otherwise it is a 
    699       tuple with composition type ('product' or 'mixture') and a list of 
    700       *model_info* blocks for the composition objects.  This allows us to 
    701       build complete product and mixture models from just the info. 
    702  
    703     """ 
    704     # TODO: maybe turn model_info into a class ModelDefinition 
    705     parameters = COMMON_PARAMETERS + kernel_module.parameters 
    706     filename = abspath(kernel_module.__file__) 
    707     kernel_id = splitext(basename(filename))[0] 
    708     name = getattr(kernel_module, 'name', None) 
    709     if name is None: 
    710         name = " ".join(w.capitalize() for w in kernel_id.split('_')) 
    711     model_info = dict( 
    712         id=kernel_id,  # string used to load the kernel 
    713         filename=abspath(kernel_module.__file__), 
    714         name=name, 
    715         title=kernel_module.title, 
    716         description=kernel_module.description, 
    717         parameters=parameters, 
    718         composition=None, 
    719         docs=kernel_module.__doc__, 
    720         category=getattr(kernel_module, 'category', None), 
    721         single=getattr(kernel_module, 'single', True), 
    722         structure_factor=getattr(kernel_module, 'structure_factor', False), 
    723         variant_info=getattr(kernel_module, 'invariant_info', None), 
    724         demo=getattr(kernel_module, 'demo', None), 
    725         source=getattr(kernel_module, 'source', []), 
    726         tests=getattr(kernel_module, 'tests', []), 
    727         ) 
    728     process_parameters(model_info) 
    729     # Check for optional functions 
    730     functions = "ER VR form_volume Iq Iqxy shape sesans".split() 
    731     model_info.update((k, getattr(kernel_module, k, None)) for k in functions) 
    732     return model_info 
    733520 
    734521section_marker = re.compile(r'\A(?P<first>[%s])(?P=first)*\Z' 
    … …  
    771558    Iq_units = "The returned value is scaled to units of |cm^-1| |sr^-1|, absolute scale." 
    772559    Sq_units = "The returned value is a dimensionless structure factor, $S(q)$." 
    773     docs = convert_section_titles_to_boldface(model_info['docs']) 
    774     subst = dict(id=model_info['id'].replace('_', '-'), 
    775                  name=model_info['name'], 
    776                  title=model_info['title'], 
    777                  parameters=make_partable(model_info['parameters']), 
    778                  returns=Sq_units if model_info['structure_factor'] else Iq_units, 
     560    docs = convert_section_titles_to_boldface(model_info.docs) 
     561    subst = dict(id=model_info.id.replace('_', '-'), 
     562                 name=model_info.name, 
     563                 title=model_info.title, 
     564                 parameters=make_partable(model_info.parameters), 
     565                 returns=Sq_units if model_info.structure_factor else Iq_units, 
    779566                 docs=docs) 
    780567    return DOC_HEADER % subst 
    … …  
    785572    Show how long it takes to process a model. 
    786573    """ 
     574    import datetime 
     575    from .modelinfo import make_model_info 
    787576    from .models import cylinder 
    788     import datetime 
     577 
    789578    tic = datetime.datetime.now() 
    790579    make_source(make_model_info(cylinder)) 
    … …  
    796585    Program which prints the source produced by the model. 
    797586    """ 
     587    import sys 
     588    from .modelinfo import make_model_info 
     589 
    798590    if len(sys.argv) <= 1: 
    799591        print("usage: python -m sasmodels.generate modelname") 

  • sasmodels/kernelcl.py

    r4d76711 r6d6508e  
    4848harmless, albeit annoying. 
    4949""" 
     50from __future__ import print_function 
    5051import os 
    5152import warnings 
    … …  
    5455 
    5556try: 
     57    raise NotImplementedError("OpenCL not yet implemented for new kernel template") 
    5658    import pyopencl as cl 
    5759    # Ask OpenCL for the default context so that we know that one exists 
    … …  
    7375# of polydisperse parameters. 
    7476MAX_LOOPS = 2048 
     77 
     78 
     79# Pragmas for enable OpenCL features.  Be sure to protect them so that they 
     80# still compile even if OpenCL is not present. 
     81_F16_PRAGMA = """\ 
     82#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 
     83#  pragma OPENCL EXTENSION cl_khr_fp16: enable 
     84#endif 
     85""" 
     86 
     87_F64_PRAGMA = """\ 
     88#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 
     89#  pragma OPENCL EXTENSION cl_khr_fp64: enable 
     90#endif 
     91""" 
    7592 
    7693 
    … …  
    142159        raise RuntimeError("%s not supported for devices"%dtype) 
    143160 
    144     source = generate.convert_type(source, dtype) 
     161    source_list = [generate.convert_type(source, dtype)] 
     162 
     163    if dtype == generate.F16: 
     164        source_list.insert(0, _F16_PRAGMA) 
     165    elif dtype == generate.F64: 
     166        source_list.insert(0, _F64_PRAGMA) 
     167 
    145168    # Note: USE_SINCOS makes the intel cpu slower under opencl 
    146169    if context.devices[0].type == cl.device_type.GPU: 
    147         source = "#define USE_SINCOS\n" + source 
     170        source_list.insert(0, "#define USE_SINCOS\n") 
    148171    options = (get_fast_inaccurate_build_options(context.devices[0]) 
    149172               if fast else []) 
     173    source = "\n".join(source_list) 
    150174    program = cl.Program(context, source).build(options=options) 
    151175    return program 
    … …  
    220244        key = "%s-%s-%s"%(name, dtype, fast) 
    221245        if key not in self.compiled: 
    222             #print("compiling",name) 
     246            print("compiling",name) 
    223247            dtype = np.dtype(dtype) 
    224             program = compile_model(self.get_context(dtype), source, dtype, fast) 
     248            program = compile_model(self.get_context(dtype), 
     249                                    str(source), dtype, fast) 
    225250            self.compiled[key] = program 
    226251        return self.compiled[key] 
    … …  
    317342        if self.program is None: 
    318343            compiler = environment().compile_program 
    319             self.program = compiler(self.info['name'], self.source, self.dtype, 
    320                                     self.fast) 
     344            self.program = compiler(self.info.name, self.source, 
     345                                    self.dtype, self.fast) 
    321346        is_2d = len(q_vectors) == 2 
    322347        kernel_name = generate.kernel_name(self.info, is_2d) 
    … …  
    329354        """ 
    330355        if self.program is not None: 
    331             environment().release_program(self.info['name']) 
     356            environment().release_program(self.info.name) 
    332357            self.program = None 
    333358 
    … …  
    365390        # at this point, so instead using 32, which is good on the set of 
    366391        # architectures tested so far. 
    367         self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 
     392        if self.is_2d: 
     393            # Note: 17 rather than 15 because results is 2 elements 
     394            # longer than input. 
     395            width = ((self.nq+17)//16)*16 
     396            self.q = np.empty((width, 2), dtype=dtype) 
     397            self.q[:self.nq, 0] = q_vectors[0] 
     398            self.q[:self.nq, 1] = q_vectors[1] 
     399        else: 
     400            # Note: 33 rather than 31 because results is 2 elements 
     401            # longer than input. 
     402            width = ((self.nq+33)//32)*32 
     403            self.q = np.empty(width, dtype=dtype) 
     404            self.q[:self.nq] = q_vectors[0] 
     405        self.global_size = [self.q.shape[0]] 
    368406        context = env.get_context(self.dtype) 
    369         self.global_size = [self.q_vectors[0].size] 
    370407        #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         ] 
     408        self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     409                             hostbuf=self.q) 
    375410 
    376411    def release(self): 
    … …  
    378413        Free the memory. 
    379414        """ 
    380         for b in self.q_buffers: 
    381             b.release() 
    382         self.q_buffers = [] 
     415        if self.q is not None: 
     416            self.q.release() 
     417            self.q = None 
    383418 
    384419    def __del__(self): 
    … …  
    406441    """ 
    407442    def __init__(self, kernel, model_info, q_vectors, dtype): 
     443        max_pd = model_info.max_pd 
     444        npars = len(model_info.parameters)-2 
    408445        q_input = GpuInput(q_vectors, dtype) 
     446        self.dtype = dtype 
     447        self.dim = '2d' if q_input.is_2d else '1d' 
    409448        self.kernel = kernel 
    410449        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] 
     450        self.pd_stop_index = 4*max_pd-1 
     451        # plus three for the normalization values 
     452        self.result = np.empty(q_input.nq+3, q_input.dtype) 
    415453 
    416454        # Inputs and outputs for each kernel call 
    … …  
    418456        env = environment() 
    419457        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, 
     458 
     459        # details is int32 data, padded to an 8 integer boundary 
     460        size = ((max_pd*5 + npars*3 + 2 + 7)//8)*8 
     461        self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
    423462                               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): 
     463        self.q_input = q_input # allocated by GpuInput above 
     464 
     465        self._need_release = [ self.result_b, self.q_input ] 
     466 
     467    def __call__(self, call_details, weights, values, cutoff): 
    429468        real = (np.float32 if self.q_input.dtype == generate.F32 
    430469                else np.float64 if self.q_input.dtype == generate.F64 
    431470                else np.float16 if self.q_input.dtype == generate.F16 
    432471                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 
     472        assert call_details.dtype == np.int32 
     473        assert weights.dtype == real and values.dtype == real 
     474 
     475        context = self.queue.context 
     476        details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     477                              hostbuf=call_details) 
     478        weights_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     479                              hostbuf=weights) 
     480        values_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     481                             hostbuf=values) 
     482 
     483        start, stop = 0, self.details[self.pd_stop_index] 
     484        args = [ 
     485            np.uint32(self.q_input.nq), np.uint32(start), np.uint32(stop), 
     486            self.details_b, self.weights_b, self.values_b, 
     487            self.q_input.q_b, self.result_b, real(cutoff), 
     488        ] 
    460489        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 
     490        cl.enqueue_copy(self.queue, self.result, self.result_b) 
     491        [v.release() for v in details_b, weights_b, values_b] 
     492 
     493        return self.result[:self.nq] 
    464494 
    465495    def release(self): 

  • sasmodels/kerneldll.py

    r4d76711 r6d6508e  
    5050import tempfile 
    5151import ctypes as ct 
    52 from ctypes import c_void_p, c_int, c_longdouble, c_double, c_float 
    53 import _ctypes 
     52from ctypes import c_void_p, c_int32, c_longdouble, c_double, c_float 
    5453 
    5554import numpy as np 
    5655 
    5756from . import generate 
     57from . import details 
    5858from .kernelpy import PyInput, PyModel 
    5959from .exception import annotate_exception 
    … …  
    8181        COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 
    8282        if "SAS_OPENMP" in os.environ: 
    83             COMPILE = COMPILE + " -fopenmp" 
     83            COMPILE += " -fopenmp" 
    8484else: 
    8585    COMPILE = "cc -shared -fPIC -fopenmp -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 
    … …  
    9090 
    9191 
    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  
     92def dll_name(model_info, dtype): 
     93    """ 
     94    Name of the dll containing the model.  This is the base file name without 
     95    any path or extension, with a form such as 'sas_sphere32'. 
     96    """ 
     97    bits = 8*dtype.itemsize 
     98    return "sas_%s%d"%(model_info.id, bits) 
     99 
     100def dll_path(model_info, dtype): 
     101    """ 
     102    Complete path to the dll for the model.  Note that the dll may not 
     103    exist yet if it hasn't been compiled. 
     104    """ 
     105    return os.path.join(DLL_PATH, dll_name(model_info, dtype)+".so") 
    106106 
    107107def make_dll(source, model_info, dtype="double"): 
    … …  
    124124    models are allowed as DLLs. 
    125125    """ 
    126     if callable(model_info.get('Iq', None)): 
     126    if callable(model_info.Iq): 
    127127        return PyModel(model_info) 
    128128     
    … …  
    133133        dtype = generate.F64  # Force 64-bit dll 
    134134 
    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   
    142135    source = generate.convert_type(source, dtype) 
    143     source_files = generate.model_sources(model_info) + [model_info['filename']] 
     136    newest = generate.timestamp(model_info) 
    144137    dll = dll_path(model_info, dtype) 
    145     newest = max(os.path.getmtime(f) for f in source_files) 
    146138    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) 
     139        basename = dll_name(model_info, dtype) + "_" 
     140        fid, filename = tempfile.mkstemp(suffix=".c", prefix=basename) 
    149141        os.fdopen(fid, "w").write(source) 
    150142        command = COMPILE%{"source":filename, "output":dll} 
    … …  
    171163    return DllModel(filename, model_info, dtype=dtype) 
    172164 
    173  
    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  
    177165class DllModel(object): 
    178166    """ 
    … …  
    197185 
    198186    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  
    204187        #print("dll", self.dllpath) 
    205188        try: 
    … …  
    212195              else c_double if self.dtype == generate.F64 
    213196              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 [] 
     197 
     198        # int, int, int, int*, double*, double*, double*, double*, double*, double 
     199        argtypes = [c_int32]*3 + [c_void_p]*5 + [fp] 
    216200        self.Iq = self.dll[generate.kernel_name(self.info, False)] 
    217         self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d 
    218  
    219201        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() 
     202        self.Iq.argtypes = argtypes 
     203        self.Iqxy.argtypes = argtypes 
    223204 
    224205    def __getstate__(self): 
    … …  
    272253    """ 
    273254    def __init__(self, kernel, model_info, q_input): 
     255        self.kernel = kernel 
    274256        self.info = model_info 
    275257        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): 
     258        self.dtype = q_input.dtype 
     259        self.dim = '2d' if q_input.is_2d else '1d' 
     260        self.result = np.empty(q_input.nq+3, q_input.dtype) 
     261 
     262    def __call__(self, call_details, weights, values, cutoff): 
    286263        real = (np.float32 if self.q_input.dtype == generate.F32 
    287264                else np.float64 if self.q_input.dtype == generate.F64 
    288265                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) 
     266        assert isinstance(call_details, details.CallDetails) 
     267        assert weights.dtype == real and values.dtype == real 
     268 
     269        start, stop = 0, call_details.total_pd 
     270        #print("in kerneldll") 
     271        #print("weights", weights) 
     272        #print("values", values) 
     273        args = [ 
     274            self.q_input.nq, # nq 
     275            start, # pd_start 
     276            stop, # pd_stop pd_stride[MAX_PD] 
     277            call_details.ctypes.data, # problem 
     278            weights.ctypes.data,  # weights 
     279            values.ctypes.data,  #pars 
     280            self.q_input.q.ctypes.data, #q 
     281            self.result.ctypes.data,   # results 
     282            real(cutoff), # cutoff 
     283            ] 
    303284        self.kernel(*args) 
    304  
    305         return self.res 
     285        return self.result[:-3] 
    306286 
    307287    def release(self): 

  • sasmodels/kernelpy.py

    r4d76711 r4bfd277  
    1010from numpy import pi, cos 
    1111 
     12from . import details 
    1213from .generate import F64 
    1314 
    … …  
    1718    """ 
    1819    def __init__(self, model_info): 
     20        # Make sure Iq and Iqxy are available and vectorized 
     21        _create_default_functions(model_info) 
    1922        self.info = model_info 
    2023 
    2124    def make_kernel(self, q_vectors): 
    2225        q_input = PyInput(q_vectors, dtype=F64) 
    23         kernel = self.info['Iqxy'] if q_input.is_2d else self.info['Iq'] 
     26        kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq 
    2427        return PyKernel(kernel, self.info, q_input) 
    2528 
    … …  
    5356        self.dtype = dtype 
    5457        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] 
     58        if self.is_2d: 
     59            self.q = np.empty((self.nq, 2), dtype=dtype) 
     60            self.q[:, 0] = q_vectors[0] 
     61            self.q[:, 1] = q_vectors[1] 
     62        else: 
     63            self.q = np.empty(self.nq, dtype=dtype) 
     64            self.q[:self.nq] = q_vectors[0] 
    5765 
    5866    def release(self): 
    … …  
    6068        Free resources associated with the model inputs. 
    6169        """ 
    62         self.q_vectors = [] 
     70        self.q = None 
    6371 
    6472class PyKernel(object): 
    … …  
    8290    """ 
    8391    def __init__(self, kernel, model_info, q_input): 
     92        self.dtype = np.dtype('d') 
    8493        self.info = model_info 
    8594        self.q_input = q_input 
    8695        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)]) 
     96        self.kernel = kernel 
     97        self.dim = '2d' if q_input.is_2d else '1d' 
     98 
     99        partable = model_info.parameters 
     100        kernel_parameters = (partable.iqxy_parameters if q_input.is_2d 
     101                             else partable.iq_parameters) 
     102        volume_parameters = partable.form_volume_parameters 
     103 
     104        # Create an array to hold the parameter values.  There will be a 
     105        # single array whose values are updated as the calculator goes 
     106        # through the loop.  Arguments to the kernel and volume functions 
     107        # will use views into this vector, relying on the fact that a 
     108        # an array of no dimensions acts like a scalar. 
     109        parameter_vector = np.empty(len(partable.call_parameters)-2, 'd') 
     110 
     111        # Create views into the array to hold the arguments 
     112        offset = 0 
     113        kernel_args, volume_args = [], [] 
     114        for p in partable.kernel_parameters: 
     115            if p.length == 1: 
     116                # Scalar values are length 1 vectors with no dimensions. 
     117                v = parameter_vector[offset:offset+1].reshape(()) 
    98118            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 
     119                # Vector values are simple views. 
     120                v = parameter_vector[offset:offset+p.length] 
     121            offset += p.length 
     122            if p in kernel_parameters: 
     123                kernel_args.append(v) 
     124            if p in volume_parameters: 
     125                volume_args.append(v) 
     126 
     127        # Hold on to the parameter vector so we can use it to call kernel later. 
     128        # This may also be required to preserve the views into the vector. 
     129        self._parameter_vector = parameter_vector 
     130 
     131        # Generate a closure which calls the kernel with the views into the 
     132        # parameter array. 
     133        if q_input.is_2d: 
     134            form = model_info.Iqxy 
     135            qx, qy = q_input.q[:,0], q_input.q[:,1] 
     136            self._form = lambda: form(qx, qy, *kernel_args) 
    106137        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  
     138            form = model_info.Iq 
     139            q = q_input.q 
     140            self._form = lambda: form(q, *kernel_args) 
     141 
     142        # Generate a closure which calls the form_volume if it exists. 
     143        form_volume = model_info.form_volume 
     144        self._volume = ((lambda: form_volume(*volume_args)) if form_volume 
     145                        else (lambda: 1.0)) 
     146 
     147    def __call__(self, call_details, weights, values, cutoff): 
     148        assert isinstance(call_details, details.CallDetails) 
     149        res = _loops(self._parameter_vector, self._form, self._volume, 
     150                     self.q_input.nq, call_details, weights, values, cutoff) 
    141151        return res 
    142152 
    … …  
    147157        self.q_input = None 
    148158 
    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  
     159def _loops(parameters,  # type: np.ndarray 
     160           form,        # type: Callable[[], np.ndarray] 
     161           form_volume, # type: Callable[[], float] 
     162           nq,          # type: int 
     163           call_details,# type: details.CallDetails 
     164           weights,     # type: np.ndarray 
     165           values,      # type: np.ndarray 
     166           cutoff,      # type: float 
     167           ):           # type: (...) -> None 
    179168    ################################################################ 
    180169    #                                                              # 
    … …  
    186175    #                                                              # 
    187176    ################################################################ 
    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] 
     177    parameters[:] = values[call_details.par_offset] 
     178    scale, background = values[0], values[1] 
     179    if call_details.num_active == 0: 
     180        norm = float(form_volume()) 
     181        if norm > 0.0: 
     182            return (scale/norm)*form() + background 
     183        else: 
     184            return np.ones(nq, 'd')*background 
     185 
     186    partial_weight = np.NaN 
     187    spherical_correction = 1.0 
     188    pd_stride = call_details.pd_stride[:call_details.num_active] 
     189    pd_length = call_details.pd_length[:call_details.num_active] 
     190    pd_offset = call_details.pd_offset[:call_details.num_active] 
     191    pd_index = np.empty_like(pd_offset) 
     192    offset = np.empty_like(call_details.par_offset) 
     193    theta = call_details.theta_par 
     194    fast_length = pd_length[0] 
     195    pd_index[0] = fast_length 
     196    total = np.zeros(nq, 'd') 
     197    norm = 0.0 
     198    for loop_index in range(call_details.total_pd): 
     199        # update polydispersity parameter values 
     200        if pd_index[0] == fast_length: 
     201            pd_index[:] = (loop_index/pd_stride)%pd_length 
     202            partial_weight = np.prod(weights[pd_offset+pd_index][1:]) 
     203            for k in range(call_details.num_coord): 
     204                par = call_details.par_coord[k] 
     205                coord = call_details.pd_coord[k] 
     206                this_offset = call_details.par_offset[par] 
     207                block_size = 1 
     208                for bit in range(len(pd_offset)): 
     209                    if coord&1: 
     210                        this_offset += block_size * pd_index[bit] 
     211                        block_size *= pd_length[bit] 
     212                    coord >>= 1 
     213                    if coord == 0: break 
     214                offset[par] = this_offset 
     215                parameters[par] = values[this_offset] 
     216                if par == theta and not (call_details.par_coord[k]&1): 
     217                    spherical_correction = max(abs(cos(pi/180 * parameters[theta])), 1e-6) 
     218        for k in range(call_details.num_coord): 
     219            if call_details.pd_coord[k]&1: 
     220                #par = call_details.par_coord[k] 
     221                parameters[par] = values[offset[par]] 
     222                #print "par",par,offset[par],parameters[par+2] 
     223                offset[par] += 1 
     224                if par == theta: 
     225                    spherical_correction = max(abs(cos(pi/180 * parameters[theta])), 1e-6) 
     226 
     227        weight = partial_weight * weights[pd_offset[0] + pd_index[0]] 
     228        pd_index[0] += 1 
     229        if weight > cutoff: 
     230            # Call the scattering function 
     231            # Assume that NaNs are only generated if the parameters are bad; 
     232            # exclude all q for that NaN.  Even better would be to have an 
     233            # INVALID expression like the C models, but that is too expensive. 
     234            I = form() 
     235            if np.isnan(I).any(): continue 
     236 
     237            # update value and norm 
     238            weight *= spherical_correction 
     239            total += weight * I 
     240            norm += weight * form_volume() 
     241 
     242    if norm > 0.0: 
     243        return (scale/norm)*total + background 
    206244    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 
     245        return np.ones(nq, 'd')*background 
     246 
     247 
     248def _create_default_functions(model_info): 
     249    """ 
     250    Autogenerate missing functions, such as Iqxy from Iq. 
     251 
     252    This only works for Iqxy when Iq is written in python. :func:`make_source` 
     253    performs a similar role for Iq written in C.  This also vectorizes 
     254    any functions that are not already marked as vectorized. 
     255    """ 
     256    _create_vector_Iq(model_info) 
     257    _create_vector_Iqxy(model_info)  # call create_vector_Iq() first 
     258 
     259 
     260def _create_vector_Iq(model_info): 
     261    """ 
     262    Define Iq as a vector function if it exists. 
     263    """ 
     264    Iq = model_info.Iq 
     265    if callable(Iq) and not getattr(Iq, 'vectorized', False): 
     266        #print("vectorizing Iq") 
     267        def vector_Iq(q, *args): 
     268            """ 
     269            Vectorized 1D kernel. 
     270            """ 
     271            return np.array([Iq(qi, *args) for qi in q]) 
     272        vector_Iq.vectorized = True 
     273        model_info.Iq = vector_Iq 
     274 
     275def _create_vector_Iqxy(model_info): 
     276    """ 
     277    Define Iqxy as a vector function if it exists, or default it from Iq(). 
     278    """ 
     279    Iq, Iqxy = model_info.Iq, model_info.Iqxy 
     280    if callable(Iqxy): 
     281        if not getattr(Iqxy, 'vectorized', False): 
     282            #print("vectorizing Iqxy") 
     283            def vector_Iqxy(qx, qy, *args): 
     284                """ 
     285                Vectorized 2D kernel. 
     286                """ 
     287                return np.array([Iqxy(qxi, qyi, *args) for qxi, qyi in zip(qx, qy)]) 
     288            vector_Iqxy.vectorized = True 
     289            model_info.Iqxy = vector_Iqxy 
     290    elif callable(Iq): 
     291        #print("defaulting Iqxy") 
     292        # Iq is vectorized because create_vector_Iq was already called. 
     293        def default_Iqxy(qx, qy, *args): 
     294            """ 
     295            Default 2D kernel. 
     296            """ 
     297            return Iq(np.sqrt(qx**2 + qy**2), *args) 
     298        default_Iqxy.vectorized = True 
     299        model_info.Iqxy = default_Iqxy 
     300 

  • 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 r6d6508e  
    1414import numpy as np 
    1515 
    16 from .generate import process_parameters, COMMON_PARAMETERS, Parameter 
    17  
    18 SCALE=0 
    19 BACKGROUND=1 
    20 EFFECT_RADIUS=2 
    21 VOLFRACTION=3 
     16from .modelinfo import Parameter, ParameterTable, ModelInfo 
    2217 
    2318def make_mixture_info(parts): 
    … …  
    3429 
    3530    # Build new parameter list 
    36     pars = COMMON_PARAMETERS + [] 
     31    pars = [] 
    3732    for k, part in enumerate(parts): 
    3833        # Parameter prefix per model, A_, B_, ... 
     34        # Note that prefix must also be applied to id and length_control 
     35        # to support vector parameters 
    3936        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_, ... 
     37        pars.append(Parameter(prefix+'scale')) 
     38        for p in part['parameters'].kernel_pars: 
     39            p = copy(p) 
     40            p.name = prefix+p.name 
     41            p.id = prefix+p.id 
     42            if p.length_control is not None: 
     43                p.length_control = prefix+p.length_control 
    5044            pars.append(p) 
     45    partable = ParameterTable(pars) 
    5146 
    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'] = [] 
     47    model_info = ModelInfo() 
     48    model_info.id = '+'.join(part['id']) 
     49    model_info.name = ' + '.join(part['name']) 
     50    model_info.filename = None 
     51    model_info.title = 'Mixture model with ' + model_info.name 
     52    model_info.description = model_info.title 
     53    model_info.docs = model_info.title 
     54    model_info.category = "custom" 
     55    model_info.parameters = partable 
     56    #model_info.single = any(part['single'] for part in parts) 
     57    model_info.structure_factor = False 
     58    model_info.variant_info = None 
     59    #model_info.tests = [] 
     60    #model_info.source = [] 
    6661    # Iq, Iqxy, form_volume, ER, VR and sesans 
    6762    # 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 
     63    model_info.composition = ('mixture', parts) 
    7164 
    7265 
    … …  
    110103        if dim == '2d': 
    111104            for p in kernels: 
    112                 partype = p.info['partype'] 
     105                partype = p.info.partype 
    113106                accumulate(partype['fixed-2d'], partype['pd-2d'], partype['volume']) 
    114107        else: 
    115108            for p in kernels: 
    116                 partype = p.info['partype'] 
     109                partype = p.info.partype 
    117110                accumulate(partype['fixed-1d'], partype['pd-1d'], partype['volume']) 
    118111 

  • sasmodels/model_test.py

    r4d76711 r4bfd277  
    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 
    … …  
    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 
     60def call_ER(model_info, pars): 
     61    """ 
     62    Call the model ER function using *values*. *model_info* is either 
     63    *model.info* if you have a loaded model, or *kernel.info* if you 
     64    have a model kernel prepared for evaluation. 
     65    """ 
     66    if model_info.ER is None: 
     67        return 1.0 
     68    else: 
     69        value, weight = _vol_pars(model_info, pars) 
     70        individual_radii = model_info.ER(*value) 
     71        return np.sum(weight*individual_radii) / np.sum(weight) 
     72 
     73def call_VR(model_info, pars): 
     74    """ 
     75    Call the model VR function using *pars*. 
     76    *info* is either *model.info* if you have a loaded model, or *kernel.info* 
     77    if you have a model kernel prepared for evaluation. 
     78    """ 
     79    if model_info.VR is None: 
     80        return 1.0 
     81    else: 
     82        value, weight = _vol_pars(model_info, pars) 
     83        whole, part = model_info.VR(*value) 
     84        return np.sum(weight*part)/np.sum(weight*whole) 
     85 
     86def _vol_pars(model_info, pars): 
     87    vol_pars = [get_weights(p, pars) 
     88                for p in model_info.parameters.call_parameters 
     89                if p.type == 'volume'] 
     90    value, weight = dispersion_mesh(model_info, vol_pars) 
     91    return value, weight 
     92 
    5793 
    5894def make_suite(loaders, models): 
    … …  
    86122        # don't try to call cl kernel since it will not be 
    87123        # available in some environmentes. 
    88         is_py = callable(model_info['Iq']) 
     124        is_py = callable(model_info.Iq) 
    89125 
    90126        if is_py:  # kernel implemented in python 
    … …  
    140176            self.dtype = dtype 
    141177 
    142             setattr(self, test_method_name, self._runTest) 
     178            setattr(self, test_method_name, self.run_all) 
    143179            unittest.TestCase.__init__(self, test_method_name) 
    144180 
    145         def _runTest(self): 
     181        def run_all(self): 
    146182            smoke_tests = [ 
    147183                [{}, 0.1, None], 
    … …  
    151187                ] 
    152188 
    153             tests = self.info['tests'] 
     189            tests = self.info.tests 
    154190            try: 
    155191                model = build_model(self.info, dtype=self.dtype, 
    156192                                    platform=self.platform) 
    157193                for test in smoke_tests + tests: 
    158                     self._run_one_test(model, test) 
     194                    self.run_one(model, test) 
    159195 
    160196                if not tests and self.platform == "dll": 
    … …  
    170206                raise 
    171207 
    172         def _run_one_test(self, model, test): 
     208        def run_one(self, model, test): 
    173209            pars, x, y = test 
     210            pars = expand_pars(self.info.parameters, pars) 
    174211 
    175212            if not isinstance(y, list): 
    … …  
    194231                actual = call_kernel(kernel, pars) 
    195232 
    196             self.assertGreater(len(actual), 0) 
     233            self.assertTrue(len(actual) > 0) 
    197234            self.assertEqual(len(y), len(actual)) 
    198235 
    … …  
    276313    tests = make_suite(['opencl', 'dll'], ['all']) 
    277314    for test_i in tests: 
    278         yield test_i._runTest 
     315        yield test_i.run_all 
    279316 
    280317 

  • sasmodels/models/cylinder.c

    r26141cb re9b1663d  
    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/flexible_cylinder.c

    re6408d0 r4937980  
    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 r4937980  
    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/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 rea05c87  
    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 r6d6508e  
    1313import numpy as np 
    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, Parameter, ModelInfo 
    1717 
    18 SCALE=0 
    19 BACKGROUND=1 
    20 RADIUS_EFFECTIVE=2 
    21 VOLFRACTION=3 
     18# TODO: make estimates available to constraints 
     19#ESTIMATED_PARAMETERS = [ 
     20#    ["est_effect_radius", "A", 0.0, [0, np.inf], "", "Estimated effective radius"], 
     21#    ["est_volume_ratio", "", 1.0, [0, np.inf], "", "Estimated volume ratio"], 
     22#] 
    2223 
    2324# TODO: core_shell_sphere model has suppressed the volume ratio calculation 
    … …  
    2728    Create info block for product model. 
    2829    """ 
    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'] = [] 
     30    p_id, p_name, p_partable = p_info.id, p_info.name, p_info.parameters 
     31    s_id, s_name, s_partable = s_info.id, s_info.name, s_info.parameters 
     32    p_set = set(p.id for p in p_partable) 
     33    s_set = set(p.id for p in s_partable) 
     34 
     35    if p_set & s_set: 
     36        # there is some overlap between the parameter names; tag the 
     37        # overlapping S parameters with name_S 
     38        s_pars = [(suffix_parameter(par, "_S") if par.id in p_set else par) 
     39                  for par in s_partable.kernel_parameters] 
     40        pars = p_partable.kernel_parameters + s_pars 
     41    else: 
     42        pars= p_partable.kernel_parameters + s_partable.kernel_parameters 
     43 
     44    model_info = ModelInfo() 
     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 %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 = ParameterTable(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 = [] 
    5858    # Iq, Iqxy, form_volume, ER, VR and sesans 
    59     model_info['composition'] = ('product', [p_info, s_info]) 
    60     process_parameters(model_info) 
     59    model_info.composition = ('product', [p_info, s_info]) 
    6160    return model_info 
    6261 
    … …  
    8887class ProductKernel(object): 
    8988    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] 
    13489        self.info = model_info 
    13590        self.p_kernel = p_kernel 
    13691        self.s_kernel = s_kernel 
    137         self.par_map = par_map 
    13892 
    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  
     93    def __call__(self, details, weights, values, cutoff): 
    15094        effect_radius, vol_ratio = call_ER_VR(self.p_kernel.info, vol_pars) 
    15195 
    … …  
    157101        s_pd[0] = [effect_radius], [1.0] 
    158102 
    159         p_res = self.p_kernel(p_fixed, p_pd, cutoff) 
    160         s_res = self.s_kernel(s_fixed, s_pd, cutoff) 
     103        p_res = self.p_kernel(p_details, p_weights, p_values, cutoff) 
     104        s_res = self.s_kernel(s_details, s_weights, s_values, cutoff) 
    161105        #print s_fixed, s_pd, p_fixed, p_pd 
    162106 
    … …  
    167111        self.q_kernel.release() 
    168112 
     113def call_ER_VR(model_info, vol_pars): 
     114    """ 
     115    Return effect radius and volume ratio for the model. 
     116    """ 
     117    value, weight = dispersion_mesh(vol_pars) 
     118 
     119    individual_radii = model_info.ER(*value) if model_info.ER else 1.0 
     120    whole, part = model_info.VR(*value) if model_info.VR else (1.0, 1.0) 
     121 
     122    effect_radius = np.sum(weight*individual_radii) / np.sum(weight) 
     123    volume_ratio = np.sum(weight*part)/np.sum(weight*whole) 
     124    return effect_radius, volume_ratio 

  • sasmodels/resolution.py

    r4d76711 r6d6508e  
    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 
    … …  
    502502    from scipy.integrate import romberg 
    503503 
    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)))) 
     504    par_set = set([p.name for p in form.info.parameters.call_parameters]) 
     505    if any(k not in par_set for k in pars.keys()): 
     506        extra = set(pars.keys()) - par_set 
     507        raise ValueError("bad parameters: [%s] not in [%s]" 
     508                         % (", ".join(sorted(extra)), 
     509                            ", ".join(sorted(pars.keys())))) 
    509510 
    510511    if np.isscalar(width): 
    … …  
    556557    from scipy.integrate import romberg 
    557558 
    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)))) 
     559    par_set = set([p.name for p in form.info.parameters.call_parameters]) 
     560    if any(k not in par_set for k in pars.keys()): 
     561        extra = set(pars.keys()) - par_set 
     562        raise ValueError("bad parameters: [%s] not in [%s]" 
     563                         % (", ".join(sorted(extra)), 
     564                            ", ".join(sorted(pars.keys())))) 
    563565 
    564566    _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq) 
    … …  
    692694 
    693695    def _eval_sphere(self, pars, resolution): 
    694         from sasmodels import core 
     696        from sasmodels import direct_model 
    695697        kernel = self.model.make_kernel([resolution.q_calc]) 
    696         theory = core.call_kernel(kernel, pars) 
     698        theory = direct_model.call_kernel(kernel, pars) 
    697699        result = resolution.apply(theory) 
    698700        kernel.release() 
    … …  
    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 

  • sasmodels/sasview_model.py

    r4d76711 r4bfd277  
    2626from . import custom 
    2727from . import generate 
     28from . import weights 
     29from . import details 
     30from . import modelinfo 
    2831 
    2932def load_standard_models(): 
    … …  
    3841        try: 
    3942            models.append(_make_standard_model(name)) 
    40         except: 
     43        except Exception: 
    4144            logging.error(traceback.format_exc()) 
    4245    return models 
    … …  
    4851    """ 
    4952    kernel_module = custom.load_custom_kernel_module(path) 
    50     model_info = generate.make_model_info(kernel_module) 
     53    model_info = modelinfo.make_model_info(kernel_module) 
    5154    return _make_model_from_info(model_info) 
    5255 
    … …  
    6164    """ 
    6265    kernel_module = generate.load_kernel_module(name) 
    63     model_info = generate.make_model_info(kernel_module) 
     66    model_info = modelinfo.make_model_info(kernel_module) 
    6467    return _make_model_from_info(model_info) 
    6568 
    … …  
    7275        SasviewModel.__init__(self) 
    7376    attrs = dict(__init__=__init__, _model_info=model_info) 
    74     ConstructedModel = type(model_info['name'], (SasviewModel,), attrs) 
     77    ConstructedModel = type(model_info.name, (SasviewModel,), attrs) 
    7578    return ConstructedModel 
    7679 
    … …  
    8083    Sasview wrapper for opencl/ctypes model. 
    8184    """ 
    82     _model_info = {} 
     85    _model_info = None # type: modelinfo.ModelInfo 
    8386    def __init__(self): 
    8487        self._model = None 
    8588        model_info = self._model_info 
    86  
    87         self.name = model_info['name'] 
    88         self.description = model_info['description'] 
     89        parameters = model_info.parameters 
     90 
     91        self.name = model_info.name 
     92        self.description = model_info.description 
    8993        self.category = None 
    90         self.multiplicity_info = None 
    91         self.is_multifunc = False 
     94        #self.is_multifunc = False 
     95        for p in parameters.kernel_parameters: 
     96            if p.is_control: 
     97                profile_axes = model_info.profile_axes 
     98                self.multiplicity_info = [ 
     99                    p.limits[1], p.name, p.choices, profile_axes[0] 
     100                    ] 
     101                break 
     102        else: 
     103            self.multiplicity_info = [] 
    92104 
    93105        ## interpret the parameters 
    … …  
    96108        self.params = collections.OrderedDict() 
    97109        self.dispersion = dict() 
    98         partype = model_info['partype'] 
    99  
    100         for p in model_info['parameters']: 
     110 
     111        self.orientation_params = [] 
     112        self.magnetic_params = [] 
     113        self.fixed = [] 
     114        for p in parameters.user_parameters(): 
    101115            self.params[p.name] = p.default 
    102116            self.details[p.name] = [p.units] + p.limits 
    103  
    104         for name in partype['pd-2d']: 
    105             self.dispersion[name] = { 
    106                 'width': 0, 
    107                 'npts': 35, 
    108                 'nsigmas': 3, 
    109                 'type': 'gaussian', 
    110             } 
    111  
    112         self.orientation_params = ( 
    113             partype['orientation'] 
    114             + [n + '.width' for n in partype['orientation']] 
    115             + partype['magnetic']) 
    116         self.magnetic_params = partype['magnetic'] 
    117         self.fixed = [n + '.width' for n in partype['pd-2d']] 
     117            if p.polydisperse: 
     118                self.dispersion[p.name] = { 
     119                    'width': 0, 
     120                    'npts': 35, 
     121                    'nsigmas': 3, 
     122                    'type': 'gaussian', 
     123                } 
     124            if p.type == 'orientation': 
     125                self.orientation_params.append(p.name) 
     126                self.orientation_params.append(p.name+".width") 
     127                self.fixed.append(p.name+".width") 
     128            if p.type == 'magnetic': 
     129                self.orientation_params.append(p.name) 
     130                self.magnetic_params.append(p.name) 
     131                self.fixed.append(p.name+".width") 
     132 
    118133        self.non_fittable = [] 
    119134 
    120135        ## independent parameter name and unit [string] 
    121         self.input_name = model_info.get("input_name", "Q") 
    122         self.input_unit = model_info.get("input_unit", "A^{-1}") 
    123         self.output_name = model_info.get("output_name", "Intensity") 
    124         self.output_unit = model_info.get("output_unit", "cm^{-1}") 
     136        self.input_name = "Q", #model_info.get("input_name", "Q") 
     137        self.input_unit = "A^{-1}" #model_info.get("input_unit", "A^{-1}") 
     138        self.output_name = "Intensity" #model_info.get("output_name", "Intensity") 
     139        self.output_unit = "cm^{-1}" #model_info.get("output_unit", "cm^{-1}") 
    125140 
    126141        ## _persistency_dict is used by sas.perspectives.fitting.basepage 
    … …  
    234249        """ 
    235250        # TODO: fix test so that parameter order doesn't matter 
    236         ret = ['%s.%s' % (d.lower(), p) 
    237                for d in self._model_info['partype']['pd-2d'] 
    238                for p in ('npts', 'nsigmas', 'width')] 
     251        ret = ['%s.%s' % (p.name.lower(), ext) 
     252               for p in self._model_info.parameters.user_parameters() 
     253               for ext in ('npts', 'nsigmas', 'width') 
     254               if p.polydisperse] 
    239255        #print(ret) 
    240256        return ret 
    … …  
    309325            # Check whether we have a list of ndarrays [qx,qy] 
    310326            qx, qy = qdist 
    311             partype = self._model_info['partype'] 
    312             if not partype['orientation'] and not partype['magnetic']: 
     327            if not self._model_info.parameters.has_2d: 
    313328                return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2)) 
    314329            else: 
    … …  
    335350            self._model = core.build_model(self._model_info) 
    336351        q_vectors = [np.asarray(q) for q in args] 
    337         fn = self._model.make_kernel(q_vectors) 
    338         pars = [self.params[v] for v in fn.fixed_pars] 
    339         pd_pars = [self._get_weights(p) for p in fn.pd_pars] 
    340         result = fn(pars, pd_pars, self.cutoff) 
    341         fn.q_input.release() 
    342         fn.release() 
     352        kernel = self._model.make_kernel(q_vectors) 
     353        pairs = [self._get_weights(p) 
     354                 for p in self._model_info.parameters.call_parameters] 
     355        call_details, weight, value = details.build_details(kernel, pairs) 
     356        result = kernel(call_details, weight, value, cutoff=self.cutoff) 
     357        kernel.q_input.release() 
     358        kernel.release() 
    343359        return result 
    344360 
    … …  
    349365        :return: the value of the effective radius 
    350366        """ 
    351         ER = self._model_info.get('ER', None) 
    352         if ER is None: 
     367        if self._model_info.ER is None: 
    353368            return 1.0 
    354369        else: 
    355             values, weights = self._dispersion_mesh() 
    356             fv = ER(*values) 
     370            value, weight = self._dispersion_mesh() 
     371            fv = self._model_info.ER(*value) 
    357372            #print(values[0].shape, weights.shape, fv.shape) 
    358             return np.sum(weights * fv) / np.sum(weights) 
     373            return np.sum(weight * fv) / np.sum(weight) 
    359374 
    360375    def calculate_VR(self): 
    … …  
    364379        :return: the value of the volf ratio 
    365380        """ 
    366         VR = self._model_info.get('VR', None) 
    367         if VR is None: 
     381        if self._model_info.VR is None: 
    368382            return 1.0 
    369383        else: 
    370             values, weights = self._dispersion_mesh() 
    371             whole, part = VR(*values) 
    372             return np.sum(weights * part) / np.sum(weights * whole) 
     384            value, weight = self._dispersion_mesh() 
     385            whole, part = self._model_info.VR(*value) 
     386            return np.sum(weight * part) / np.sum(weight * whole) 
    373387 
    374388    def set_dispersion(self, parameter, dispersion): 
    … …  
    408422        parameter set in the vector. 
    409423        """ 
    410         pars = self._model_info['partype']['volume'] 
    411         return core.dispersion_mesh([self._get_weights(p) for p in pars]) 
     424        pars = [self._get_weights(p) 
     425                for p in self._model_info.parameters.call_parameters 
     426                if p.type == 'volume'] 
     427        return details.dispersion_mesh(self._model_info, pars) 
    412428 
    413429    def _get_weights(self, par): 
    … …  
    415431        Return dispersion weights for parameter 
    416432        """ 
    417         from . import weights 
    418         relative = self._model_info['partype']['pd-rel'] 
    419         limits = self._model_info['limits'] 
    420         dis = self.dispersion[par] 
    421         value, weight = weights.get_weights( 
    422             dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 
    423             self.params[par], limits[par], par in relative) 
    424         return value, weight / np.sum(weight) 
    425  
     433        if par.polydisperse: 
     434            dis = self.dispersion[par.name] 
     435            value, weight = weights.get_weights( 
     436                dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 
     437                self.params[par.name], par.limits, par.relative_pd) 
     438            return value, weight / np.sum(weight) 
     439        else: 
     440            return [self.params[par.name]], [] 
    426441 
    427442def test_model(): 
Note: See TracChangeset for help on using the changeset viewer.