Changes in / [ec45c4f:ea05c87] in sasmodels


Ignore:
Files:
7 added
27 edited

Legend:

Unmodified
Added
Removed

  • doc/developer/index.rst

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

  • sasmodels/bumps_model.py

    rea75043 rea75043  
    8181    from bumps.names import Parameter 
    8282 
    83     pars = {} 
     83    pars = {}     # floating point parameters 
     84    pd_types = {} # distribution names 
    8485    for p in model_info['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            pd_types[p.name+'_pd_type'] = kwargs.pop(name, 'gaussian') 
    10298 
    10399    if kwargs:  # args not corresponding to parameters 

  • sasmodels/compare.py

    rf247314 rf247314  
    306306    Format the parameter list for printing. 
    307307    """ 
    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 
    314308    lines = [] 
    315     for p in model_info['parameters']: 
    316         if exclude(p.name): continue 
     309    parameters = model_info['parameters'] 
     310    for p in parameters.user_parameters(pars, is2d): 
    317311        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')) 
     312            value=pars.get(p.id, p.default), 
     313            pd=pars.get(p.id+"_pd", 0.), 
     314            n=int(pars.get(p.id+"_pd_n", 0)), 
     315            nsigma=pars.get(p.id+"_pd_nsgima", 3.), 
     316            type=pars.get(p.id+"_pd_type", 'gaussian')) 
    323317        lines.append(_format_par(p.name, **fields)) 
    324318    return "\n".join(lines) 
    … …  
    454448    """ 
    455449    # initialize the code so time is more accurate 
    456     value = calculator(**suppress_pd(pars)) 
     450    if Nevals > 1: 
     451        value = calculator(**suppress_pd(pars)) 
    457452    toc = tic() 
    458453    for _ in range(max(Nevals, 1)):  # make sure there is at least one eval 
    … …  
    661656    """ 
    662657    # 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" 
     658    pars = {} 
     659    for p in model_info['parameters'].call_parameters: 
     660        parts = [('', p.default)] 
     661        if p.polydisperse: 
     662            parts.append(('_pd', 0.0)) 
     663            parts.append(('_pd_n', 0)) 
     664            parts.append(('_pd_nsigma', 3.0)) 
     665            parts.append(('_pd_type', "gaussian")) 
     666        for ext, val in parts: 
     667            if p.length > 1: 
     668                dict(("%s%d%s"%(p.id,k,ext), val) for k in range(p.length)) 
     669            else: 
     670                pars[p.id+ext] = val 
    672671 
    673672    # Plug in values given in demo 
    … …  
    774773 
    775774    if len(engines) == 0: 
    776         engines.extend(['single', 'sasview']) 
     775        engines.extend(['single', 'double']) 
    777776    elif len(engines) == 1: 
    778         if engines[0][0] != 'sasview': 
    779             engines.append('sasview') 
     777        if engines[0][0] != 'double': 
     778            engines.append('double') 
    780779        else: 
    781780            engines.append('single') 
    … …  
    880879        pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars']) 
    881880        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)) 
     881            for p in model_info['parameters'].type['1d']: 
     882                for ext in ['', '_pd', '_pd_n', '_pd_nsigma']: 
     883                    k = p.name+ext 
     884                    v = pars.get(k, None) 
     885                    if v is not None: 
     886                        v.range(*parameter_range(k, v.value)) 
    889887        else: 
    890888            for k, v in pars.items(): 

  • sasmodels/core.py

    rb7172bb r68e7f9d  
    77    ] 
    88 
    9  
    109from os.path import basename, dirname, join as joinpath, splitext 
    1110from glob import glob 
    12 import imp 
    1311 
    1412import numpy as np 
    … …  
    4038            return [np.asarray(v) for v in args] 
    4139 
    42  
    4340def list_models(): 
    4441    """ 
    … …  
    6360    """ 
    6461    return build_model(load_model_info(model_name), **kw) 
    65  
    66 def load_model_info_from_path(path): 
    67     # Pull off the last .ext if it exists; there may be others 
    68     name = basename(splitext(path)[0]) 
    69  
    70     # Not cleaning name since don't need to be able to reload this 
    71     # model later 
    72     # Should probably turf the model from sys.modules after we are done... 
    73  
    74     # Placing the model in the 'sasmodels.custom' name space, even 
    75     # though it doesn't actually exist.  imp.load_source doesn't seem 
    76     # to care. 
    77     kernel_module = imp.load_source('sasmodels.custom.'+name, path) 
    78  
    79     # Now that we have the module, we can load it as usual 
    80     model_info = generate.make_model_info(kernel_module) 
    81     return model_info 
    8262 
    8363def load_model_info(model_name): 
    … …  
    10282        return product.make_product_info(P_info, Q_info) 
    10383 
    104     #import sys; print "\n".join(sys.path) 
    105     __import__('sasmodels.models.'+model_name) 
    106     kernel_module = getattr(models, model_name, None) 
     84    kernel_module = generate.load_kernel_module(model_name) 
    10785    return generate.make_model_info(kernel_module) 
    10886 
    … …  
    179157 
    180158 
    181 def make_kernel(model, q_vectors): 
    182     """ 
    183     Return a computation kernel from the model definition and the q input. 
    184     """ 
    185     return model(q_vectors) 
    186  
    187 def get_weights(model_info, pars, name): 
     159def get_weights(parameter, values): 
    188160    """ 
    189161    Generate the distribution for parameter *name* given the parameter values 
    … …  
    193165    from the *pars* dictionary for parameter value and parameter dispersion. 
    194166    """ 
    195     relative = name in model_info['partype']['pd-rel'] 
    196     limits = model_info['limits'][name] 
    197     disperser = pars.get(name+'_pd_type', 'gaussian') 
    198     value = pars.get(name, model_info['defaults'][name]) 
    199     npts = pars.get(name+'_pd_n', 0) 
    200     width = pars.get(name+'_pd', 0.0) 
    201     nsigma = pars.get(name+'_pd_nsigma', 3.0) 
     167    value = values.get(parameter.name, parameter.default) 
     168    relative = parameter.relative_pd 
     169    limits = parameter.limits 
     170    disperser = values.get(parameter.name+'_pd_type', 'gaussian') 
     171    npts = values.get(parameter.name+'_pd_n', 0) 
     172    width = values.get(parameter.name+'_pd', 0.0) 
     173    nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 
     174    if npts == 0 or width == 0: 
     175        return [value], [] 
    202176    value, weight = weights.get_weights( 
    203177        disperser, npts, width, nsigma, value, limits, relative) 
    … …  
    220194def call_kernel(kernel, pars, cutoff=0, mono=False): 
    221195    """ 
    222     Call *kernel* returned from :func:`make_kernel` with parameters *pars*. 
     196    Call *kernel* returned from *model.make_kernel* with parameters *pars*. 
    223197 
    224198    *cutoff* is the limiting value for the product of dispersion weights used 
    … …  
    228202    with an error of about 1%, which is usually less than the measurement 
    229203    uncertainty. 
    230     """ 
    231     fixed_pars = [pars.get(name, kernel.info['defaults'][name]) 
    232                   for name in kernel.fixed_pars] 
     204 
     205    *mono* is True if polydispersity should be set to none on all parameters. 
     206    """ 
     207    parameters = kernel.info['parameters'] 
    233208    if mono: 
    234         pd_pars = [( np.array([pars[name]]), np.array([1.0]) ) 
    235                    for name in kernel.pd_pars] 
    236     else: 
    237         pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars] 
    238     return kernel(fixed_pars, pd_pars, cutoff=cutoff) 
     209        active = lambda name: False 
     210    elif kernel.dim == '1d': 
     211        active = lambda name: name in parameters.pd_1d 
     212    elif kernel.dim == '2d': 
     213        active = lambda name: name in parameters.pd_2d 
     214    else: 
     215        active = lambda name: True 
     216 
     217    vw_pairs = [(get_weights(p, pars) if active(p.name) 
     218                 else ([pars.get(p.name, p.default)], [])) 
     219                for p in parameters.call_parameters] 
     220 
     221    details, weights, values = build_details(kernel, vw_pairs) 
     222    return kernel(details, weights, values, cutoff) 
     223 
     224def build_details(kernel, pairs): 
     225    values, weights = zip(*pairs) 
     226    if max([len(w) for w in weights]) > 1: 
     227        details = generate.poly_details(kernel.info, weights) 
     228    else: 
     229        details = kernel.info['mono_details'] 
     230    weights, values = [np.hstack(v) for v in (weights, values)] 
     231    weights = weights.astype(dtype=kernel.dtype) 
     232    values = values.astype(dtype=kernel.dtype) 
     233    return details, weights, values 
     234 
    239235 
    240236def call_ER_VR(model_info, vol_pars): 
    … …  
    259255 
    260256 
    261 def call_ER(info, pars): 
    262     """ 
    263     Call the model ER function using *pars*. 
    264     *info* is either *model.info* if you have a loaded model, or *kernel.info* 
    265     if you have a model kernel prepared for evaluation. 
    266     """ 
    267     ER = info.get('ER', None) 
     257def call_ER(model_info, values): 
     258    """ 
     259    Call the model ER function using *values*. *model_info* is either 
     260    *model.info* if you have a loaded model, or *kernel.info* if you 
     261    have a model kernel prepared for evaluation. 
     262    """ 
     263    ER = model_info.get('ER', None) 
    268264    if ER is None: 
    269265        return 1.0 
    270266    else: 
    271         vol_pars = [get_weights(info, pars, name) 
    272                     for name in info['partype']['volume']] 
     267        vol_pars = [get_weights(parameter, values) 
     268                    for parameter in model_info['parameters'] 
     269                    if parameter.type == 'volume'] 
    273270        value, weight = dispersion_mesh(vol_pars) 
    274271        individual_radii = ER(*value) 
    … …  
    276273        return np.sum(weight*individual_radii) / np.sum(weight) 
    277274 
    278 def call_VR(info, pars): 
     275def call_VR(model_info, values): 
    279276    """ 
    280277    Call the model VR function using *pars*. 
    … …  
    282279    if you have a model kernel prepared for evaluation. 
    283280    """ 
    284     VR = info.get('VR', None) 
     281    VR = model_info.get('VR', None) 
    285282    if VR is None: 
    286283        return 1.0 
    287284    else: 
    288         vol_pars = [get_weights(info, pars, name) 
    289                     for name in info['partype']['volume']] 
     285        vol_pars = [get_weights(parameter, values) 
     286                    for parameter in model_info['parameters'] 
     287                    if parameter.type == 'volume'] 
    290288        value, weight = dispersion_mesh(vol_pars) 
    291289        whole, part = VR(*value) 

  • sasmodels/direct_model.py

    rea75043 r68e7f9d  
    2525import numpy as np 
    2626 
    27 from .core import make_kernel 
    2827from .core import call_kernel, call_ER_VR 
    2928from . import sesans 
    … …  
    6867            self.data_type = 'Iq' 
    6968 
    70         partype = model.info['partype'] 
    71  
    7269        if self.data_type == 'sesans': 
    7370             
    … …  
    8380            q_mono = sesans.make_all_q(data) 
    8481        elif self.data_type == 'Iqxy': 
    85             if not partype['orientation'] and not partype['magnetic']: 
    86                 raise ValueError("not 2D without orientation or magnetic parameters") 
     82            #if not model.info['parameters'].has_2d: 
     83            #    raise ValueError("not 2D without orientation or magnetic parameters") 
    8784            q = np.sqrt(data.qx_data**2 + data.qy_data**2) 
    8885            qmin = getattr(data, 'qmin', 1e-16) 
    … …  
    173170    def _calc_theory(self, pars, cutoff=0.0): 
    174171        if self._kernel is None: 
    175             self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-dedata_type 
    176             self._kernel_mono = (make_kernel(self._model, self._kernel_mono_inputs) 
     172            self._kernel = self._model.make_kernel(self._kernel_inputs) 
     173            self._kernel_mono = (self._model.make_kernel(self._kernel_mono_inputs) 
    177174                                 if self._kernel_mono_inputs else None) 
    178175 

  • sasmodels/generate.py

    rf247314 rf247314  
    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 
    … …  
    204173from __future__ import print_function 
    205174 
    206 # TODO: identify model files which have changed since loading and reload them. 
    207  
    208 import sys 
     175#TODO: determine which functions are useful outside of generate 
     176#__all__ = ["model_info", "make_doc", "make_source", "convert_type"] 
     177 
    209178from os.path import abspath, dirname, join as joinpath, exists, basename, \ 
    210     splitext 
     179    splitext, getmtime 
    211180import re 
    212181import string 
    213182import warnings 
    214 from collections import namedtuple 
    215183 
    216184import numpy as np 
    217185 
    218 PARAMETER_FIELDS = ['name', 'units', 'default', 'limits', 'type', 'description'] 
    219 Parameter = namedtuple('Parameter', PARAMETER_FIELDS) 
    220  
    221 #TODO: determine which functions are useful outside of generate 
    222 #__all__ = ["model_info", "make_doc", "make_source", "convert_type"] 
    223  
    224 C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c') 
     186from .modelinfo import ModelInfo, Parameter, make_parameter_table, set_demo 
     187from .custom import load_custom_kernel_module 
     188 
     189# TODO: identify model files which have changed since loading and reload them. 
     190 
     191TEMPLATE_ROOT = dirname(__file__) 
    225192 
    226193F16 = np.dtype('float16') 
    … …  
    231198except TypeError: 
    232199    F128 = None 
    233  
    234 # Scale and background, which are parameters common to every form factor 
    235 COMMON_PARAMETERS = [ 
    236     ["scale", "", 1, [0, np.inf], "", "Source intensity"], 
    237     ["background", "1/cm", 1e-3, [0, np.inf], "", "Source background"], 
    238     ] 
    239200 
    240201# Conversion from units defined in the parameter table for each model 
    … …  
    330291    raise ValueError("%r not found in %s" % (filename, search_path)) 
    331292 
     293 
    332294def model_sources(model_info): 
    333295    """ 
    … …  
    338300    return [_search(search_path, f) for f in model_info['source']] 
    339301 
    340 # Pragmas for enable OpenCL features.  Be sure to protect them so that they 
    341 # still compile even if OpenCL is not present. 
    342 _F16_PRAGMA = """\ 
    343 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 
    344 #  pragma OPENCL EXTENSION cl_khr_fp16: enable 
    345 #endif 
    346 """ 
    347  
    348 _F64_PRAGMA = """\ 
    349 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 
    350 #  pragma OPENCL EXTENSION cl_khr_fp64: enable 
    351 #endif 
    352 """ 
     302def timestamp(model_info): 
     303    """ 
     304    Return a timestamp for the model corresponding to the most recently 
     305    changed file or dependency. 
     306    """ 
     307    source_files = (model_sources(model_info) 
     308                    + model_templates() 
     309                    + [model_info['filename']]) 
     310    newest = max(getmtime(f) for f in source_files) 
     311    return newest 
    353312 
    354313def convert_type(source, dtype): 
    … …  
    361320    if dtype == F16: 
    362321        fbytes = 2 
    363         source = _F16_PRAGMA + _convert_type(source, "half", "f") 
     322        source = _convert_type(source, "float", "f") 
    364323    elif dtype == F32: 
    365324        fbytes = 4 
    … …  
    367326    elif dtype == F64: 
    368327        fbytes = 8 
    369         source = _F64_PRAGMA + source  # Source is already double 
     328        # no need to convert if it is already double 
    370329    elif dtype == F128: 
    371330        fbytes = 16 
    … …  
    410369 
    411370 
    412 LOOP_OPEN = """\ 
    413 for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) { 
    414   const double %(name)s = loops[2*(%(name)s_i%(offset)s)]; 
    415   const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\ 
     371_template_cache = {} 
     372def load_template(filename): 
     373    path = joinpath(TEMPLATE_ROOT, filename) 
     374    mtime = getmtime(path) 
     375    if filename not in _template_cache or mtime > _template_cache[filename][0]: 
     376        with open(path) as fid: 
     377            _template_cache[filename] = (mtime, fid.read(), path) 
     378    return _template_cache[filename][1] 
     379 
     380def model_templates(): 
     381    # TODO: fails DRY; templates are listed in two places. 
     382    # should instead have model_info contain a list of paths 
     383    return [joinpath(TEMPLATE_ROOT, filename) 
     384            for filename in ('kernel_header.c', 'kernel_iq.c')] 
     385 
     386 
     387_FN_TEMPLATE = """\ 
     388double %(name)s(%(pars)s); 
     389double %(name)s(%(pars)s) { 
     390    %(body)s 
     391} 
     392 
     393 
    416394""" 
    417 def build_polydispersity_loops(pd_pars): 
    418     """ 
    419     Build polydispersity loops 
    420  
    421     Returns loop opening and loop closing 
    422     """ 
    423     depth = 4 
    424     offset = "" 
    425     loop_head = [] 
    426     loop_end = [] 
    427     for name in pd_pars: 
    428         subst = {'name': name, 'offset': offset} 
    429         loop_head.append(indent(LOOP_OPEN % subst, depth)) 
    430         loop_end.insert(0, (" "*depth) + "}") 
    431         offset += '+N' + name 
    432         depth += 2 
    433     return "\n".join(loop_head), "\n".join(loop_end) 
    434  
    435 C_KERNEL_TEMPLATE = None 
     395 
     396def _gen_fn(name, pars, body): 
     397    """ 
     398    Generate a function given pars and body. 
     399 
     400    Returns the following string:: 
     401 
     402         double fn(double a, double b, ...); 
     403         double fn(double a, double b, ...) { 
     404             .... 
     405         } 
     406    """ 
     407    par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void' 
     408    return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl} 
     409 
     410def _call_pars(prefix, pars): 
     411    """ 
     412    Return a list of *prefix.parameter* from parameter items. 
     413    """ 
     414    return [p.as_call_reference(prefix) for p in pars] 
     415 
     416_IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)", 
     417                           flags=re.MULTILINE) 
     418def _have_Iqxy(sources): 
     419    """ 
     420    Return true if any file defines Iqxy. 
     421 
     422    Note this is not a C parser, and so can be easily confused by 
     423    non-standard syntax.  Also, it will incorrectly identify the following 
     424    as having Iqxy:: 
     425 
     426        /* 
     427        double Iqxy(qx, qy, ...) { ... fill this in later ... } 
     428        */ 
     429 
     430    If you want to comment out an Iqxy function, use // on the front of the 
     431    line instead. 
     432    """ 
     433    for code in sources: 
     434        if _IQXY_PATTERN.search(code): 
     435            return True 
     436    else: 
     437        return False 
     438 
    436439def make_source(model_info): 
    437440    """ 
    … …  
    452455    # for computing volume even if we allow non-disperse volume parameters. 
    453456 
    454     # Load template 
    455     global C_KERNEL_TEMPLATE 
    456     if C_KERNEL_TEMPLATE is None: 
    457         with open(C_KERNEL_TEMPLATE_PATH) as fid: 
    458             C_KERNEL_TEMPLATE = fid.read() 
    459  
    460     # Load additional sources 
    461     source = [open(f).read() for f in model_sources(model_info)] 
    462  
    463     # Prepare defines 
    464     defines = [] 
    465     partype = model_info['partype'] 
    466     pd_1d = partype['pd-1d'] 
    467     pd_2d = partype['pd-2d'] 
    468     fixed_1d = partype['fixed-1d'] 
    469     fixed_2d = partype['fixed-1d'] 
    470  
    471     iq_parameters = [p.name 
    472                      for p in model_info['parameters'][2:]  # skip scale, background 
    473                      if p.name in set(fixed_1d + pd_1d)] 
    474     iqxy_parameters = [p.name 
    475                        for p in model_info['parameters'][2:]  # skip scale, background 
    476                        if p.name in set(fixed_2d + pd_2d)] 
    477     volume_parameters = [p.name 
    478                          for p in model_info['parameters'] 
    479                          if p.type == 'volume'] 
    480  
    481     # Fill in defintions for volume parameters 
    482     if volume_parameters: 
    483         defines.append(('VOLUME_PARAMETERS', 
    484                         ','.join(volume_parameters))) 
    485         defines.append(('VOLUME_WEIGHT_PRODUCT', 
    486                         '*'.join(p + '_w' for p in volume_parameters))) 
    487  
    488     # Generate form_volume function from body only 
     457    partable = model_info['parameters'] 
     458 
     459    # Identify parameters for Iq, Iqxy, Iq_magnetic and form_volume. 
     460    # Note that scale and volume are not possible types. 
     461 
     462    # Load templates and user code 
     463    kernel_header = load_template('kernel_header.c') 
     464    kernel_code = load_template('kernel_iq.c') 
     465    user_code = [open(f).read() for f in model_sources(model_info)] 
     466 
     467    # Build initial sources 
     468    source = [kernel_header] + user_code 
     469 
     470    # Make parameters for q, qx, qy so that we can use them in declarations 
     471    q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 
     472    # Generate form_volume function, etc. from body only 
    489473    if model_info['form_volume'] is not None: 
    490         if volume_parameters: 
    491             vol_par_decl = ', '.join('double ' + p for p in volume_parameters) 
    492         else: 
    493             vol_par_decl = 'void' 
    494         defines.append(('VOLUME_PARAMETER_DECLARATIONS', 
    495                         vol_par_decl)) 
    496         fn = """\ 
    497 double form_volume(VOLUME_PARAMETER_DECLARATIONS); 
    498 double form_volume(VOLUME_PARAMETER_DECLARATIONS) { 
    499     %(body)s 
    500 } 
    501 """ % {'body':model_info['form_volume']} 
    502         source.append(fn) 
    503  
    504     # Fill in definitions for Iq parameters 
    505     defines.append(('IQ_KERNEL_NAME', model_info['name'] + '_Iq')) 
    506     defines.append(('IQ_PARAMETERS', ', '.join(iq_parameters))) 
    507     if fixed_1d: 
    508         defines.append(('IQ_FIXED_PARAMETER_DECLARATIONS', 
    509                         ', \\\n    '.join('const double %s' % p for p in fixed_1d))) 
    510     if pd_1d: 
    511         defines.append(('IQ_WEIGHT_PRODUCT', 
    512                         '*'.join(p + '_w' for p in pd_1d))) 
    513         defines.append(('IQ_DISPERSION_LENGTH_DECLARATIONS', 
    514                         ', \\\n    '.join('const int N%s' % p for p in pd_1d))) 
    515         defines.append(('IQ_DISPERSION_LENGTH_SUM', 
    516                         '+'.join('N' + p for p in pd_1d))) 
    517         open_loops, close_loops = build_polydispersity_loops(pd_1d) 
    518         defines.append(('IQ_OPEN_LOOPS', 
    519                         open_loops.replace('\n', ' \\\n'))) 
    520         defines.append(('IQ_CLOSE_LOOPS', 
    521                         close_loops.replace('\n', ' \\\n'))) 
     474        pars = partable.form_volume_parameters 
     475        source.append(_gen_fn('form_volume', pars, model_info['form_volume'])) 
    522476    if model_info['Iq'] is not None: 
    523         defines.append(('IQ_PARAMETER_DECLARATIONS', 
    524                         ', '.join('double ' + p for p in iq_parameters))) 
    525         fn = """\ 
    526 double Iq(double q, IQ_PARAMETER_DECLARATIONS); 
    527 double Iq(double q, IQ_PARAMETER_DECLARATIONS) { 
    528     %(body)s 
    529 } 
    530 """ % {'body':model_info['Iq']} 
    531         source.append(fn) 
    532  
    533     # Fill in definitions for Iqxy parameters 
    534     defines.append(('IQXY_KERNEL_NAME', model_info['name'] + '_Iqxy')) 
    535     defines.append(('IQXY_PARAMETERS', ', '.join(iqxy_parameters))) 
    536     if fixed_2d: 
    537         defines.append(('IQXY_FIXED_PARAMETER_DECLARATIONS', 
    538                         ', \\\n    '.join('const double %s' % p for p in fixed_2d))) 
    539     if pd_2d: 
    540         defines.append(('IQXY_WEIGHT_PRODUCT', 
    541                         '*'.join(p + '_w' for p in pd_2d))) 
    542         defines.append(('IQXY_DISPERSION_LENGTH_DECLARATIONS', 
    543                         ', \\\n    '.join('const int N%s' % p for p in pd_2d))) 
    544         defines.append(('IQXY_DISPERSION_LENGTH_SUM', 
    545                         '+'.join('N' + p for p in pd_2d))) 
    546         open_loops, close_loops = build_polydispersity_loops(pd_2d) 
    547         defines.append(('IQXY_OPEN_LOOPS', 
    548                         open_loops.replace('\n', ' \\\n'))) 
    549         defines.append(('IQXY_CLOSE_LOOPS', 
    550                         close_loops.replace('\n', ' \\\n'))) 
     477        pars = [q] + partable.iq_parameters 
     478        source.append(_gen_fn('Iq', pars, model_info['Iq'])) 
    551479    if model_info['Iqxy'] is not None: 
    552         defines.append(('IQXY_PARAMETER_DECLARATIONS', 
    553                         ', '.join('double ' + p for p in iqxy_parameters))) 
    554         fn = """\ 
    555 double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS); 
    556 double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS) { 
    557     %(body)s 
    558 } 
    559 """ % {'body':model_info['Iqxy']} 
    560         source.append(fn) 
    561  
    562     # Need to know if we have a theta parameter for Iqxy; it is not there 
    563     # for the magnetic sphere model, for example, which has a magnetic 
    564     # orientation but no shape orientation. 
    565     if 'theta' in pd_2d: 
    566         defines.append(('IQXY_HAS_THETA', '1')) 
    567  
    568     #for d in defines: print(d) 
    569     defines = '\n'.join('#define %s %s' % (k, v) for k, v in defines) 
    570     sources = '\n\n'.join(source) 
    571     return C_KERNEL_TEMPLATE % { 
    572         'DEFINES': defines, 
    573         'SOURCES': sources, 
    574         } 
    575  
    576 def categorize_parameters(pars): 
    577     """ 
    578     Build parameter categories out of the the parameter definitions. 
    579  
    580     Returns a dictionary of categories. 
    581  
    582     Note: these categories are subject to change, depending on the needs of 
    583     the UI and the needs of the kernel calling function. 
    584  
    585     The categories are as follows: 
    586  
    587     * *volume* list of volume parameter names 
    588     * *orientation* list of orientation parameters 
    589     * *magnetic* list of magnetic parameters 
    590     * *<empty string>* list of parameters that have no type info 
    591  
    592     Each parameter is in one and only one category. 
    593  
    594     The following derived categories are created: 
    595  
    596     * *fixed-1d* list of non-polydisperse parameters for 1D models 
    597     * *pd-1d* list of polydisperse parameters for 1D models 
    598     * *fixed-2d* list of non-polydisperse parameters for 2D models 
    599     * *pd-d2* list of polydisperse parameters for 2D models 
    600     """ 
    601     partype = { 
    602         'volume': [], 'orientation': [], 'magnetic': [], 'sld': [], '': [], 
    603         'fixed-1d': [], 'fixed-2d': [], 'pd-1d': [], 'pd-2d': [], 
    604         'pd-rel': set(), 
    605     } 
    606  
    607     for p in pars: 
    608         if p.type == 'volume': 
    609             partype['pd-1d'].append(p.name) 
    610             partype['pd-2d'].append(p.name) 
    611             partype['pd-rel'].add(p.name) 
    612         elif p.type == 'magnetic': 
    613             partype['fixed-2d'].append(p.name) 
    614         elif p.type == 'orientation': 
    615             partype['pd-2d'].append(p.name) 
    616         elif p.type in ('', 'sld'): 
    617             partype['fixed-1d'].append(p.name) 
    618             partype['fixed-2d'].append(p.name) 
    619         else: 
    620             raise ValueError("unknown parameter type %r" % p.type) 
    621         partype[p.type].append(p.name) 
    622  
    623     return partype 
    624  
    625 def process_parameters(model_info): 
    626     """ 
    627     Process parameter block, precalculating parameter details. 
    628     """ 
    629     # convert parameters into named tuples 
    630     for p in model_info['parameters']: 
    631         if p[4] == '' and (p[0].startswith('sld') or p[0].endswith('sld')): 
    632             p[4] = 'sld' 
    633             # TODO: make sure all models explicitly label their sld parameters 
    634             #raise ValueError("%s.%s needs to be explicitly set to type 'sld'" %(model_info['id'], p[0])) 
    635  
    636     pars = [Parameter(*p) for p in model_info['parameters']] 
    637     # Fill in the derived attributes 
    638     model_info['parameters'] = pars 
    639     partype = categorize_parameters(pars) 
    640     model_info['limits'] = dict((p.name, p.limits) for p in pars) 
    641     model_info['partype'] = partype 
    642     model_info['defaults'] = dict((p.name, p.default) for p in pars) 
    643     if model_info.get('demo', None) is None: 
    644         model_info['demo'] = model_info['defaults'] 
    645     model_info['has_2d'] = partype['orientation'] or partype['magnetic'] 
     480        pars = [qx, qy] + partable.iqxy_parameters 
     481        source.append(_gen_fn('Iqxy', pars, model_info['Iqxy'])) 
     482 
     483    # Define the parameter table 
     484    source.append("#define PARAMETER_TABLE \\") 
     485    source.append("\\\n".join(p.as_definition() 
     486                              for p in partable.kernel_parameters)) 
     487 
     488    # Define the function calls 
     489    if partable.form_volume_parameters: 
     490        refs = _call_pars("v.", partable.form_volume_parameters) 
     491        call_volume = "#define CALL_VOLUME(v) form_volume(%s)" % (",".join(refs)) 
     492    else: 
     493        # Model doesn't have volume.  We could make the kernel run a little 
     494        # faster by not using/transferring the volume normalizations, but 
     495        # the ifdef's reduce readability more than is worthwhile. 
     496        call_volume = "#define CALL_VOLUME(v) 1.0" 
     497    source.append(call_volume) 
     498 
     499    refs = ["q[i]"] + _call_pars("v.", partable.iq_parameters) 
     500    call_iq = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(refs)) 
     501    if _have_Iqxy(user_code): 
     502        # Call 2D model 
     503        refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("v.", partable.iqxy_parameters) 
     504        call_iqxy = "#define CALL_IQ(q,i,v) Iqxy(%s)" % (",".join(refs)) 
     505    else: 
     506        # Call 1D model with sqrt(qx^2 + qy^2) 
     507        warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 
     508        # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 
     509        pars_sqrt = ["sqrt(q[2*i]*q[2*i]+q[2*i+1]*q[2*i+1])"] + refs[1:] 
     510        call_iqxy = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(pars_sqrt)) 
     511 
     512    # Fill in definitions for numbers of parameters 
     513    source.append("#define MAX_PD %s"%partable.max_pd) 
     514    source.append("#define NPARS %d"%partable.npars) 
     515 
     516    # TODO: allow mixed python/opencl kernels? 
     517 
     518    # define the Iq kernel 
     519    source.append("#define KERNEL_NAME %s_Iq"%model_info['name']) 
     520    source.append(call_iq) 
     521    source.append(kernel_code) 
     522    source.append("#undef CALL_IQ") 
     523    source.append("#undef KERNEL_NAME") 
     524 
     525    # define the Iqxy kernel from the same source with different #defines 
     526    source.append("#define KERNEL_NAME %s_Iqxy"%model_info['name']) 
     527    source.append(call_iqxy) 
     528    source.append(kernel_code) 
     529    source.append("#undef CALL_IQ") 
     530    source.append("#undef KERNEL_NAME") 
     531 
     532    return '\n'.join(source) 
     533 
     534class CoordinationDetails(object): 
     535    def __init__(self, model_info): 
     536        parameters = model_info['parameters'] 
     537        max_pd = parameters.max_pd 
     538        npars = parameters.npars 
     539        par_offset = 4*max_pd 
     540        self.details = np.zeros(par_offset + 3*npars + 4, 'i4') 
     541 
     542        # generate views on different parts of the array 
     543        self._pd_par     = self.details[0*max_pd:1*max_pd] 
     544        self._pd_length  = self.details[1*max_pd:2*max_pd] 
     545        self._pd_offset  = self.details[2*max_pd:3*max_pd] 
     546        self._pd_stride  = self.details[3*max_pd:4*max_pd] 
     547        self._par_offset = self.details[par_offset+0*npars:par_offset+1*npars] 
     548        self._par_coord  = self.details[par_offset+1*npars:par_offset+2*npars] 
     549        self._pd_coord   = self.details[par_offset+2*npars:par_offset+3*npars] 
     550 
     551        # theta_par is fixed 
     552        self.details[-1] = parameters.theta_offset 
     553 
     554    @property 
     555    def ctypes(self): return self.details.ctypes 
     556    @property 
     557    def pd_par(self): return self._pd_par 
     558    @property 
     559    def pd_length(self): return self._pd_length 
     560    @property 
     561    def pd_offset(self): return self._pd_offset 
     562    @property 
     563    def pd_stride(self): return self._pd_stride 
     564    @property 
     565    def pd_coord(self): return self._pd_coord 
     566    @property 
     567    def par_coord(self): return self._par_coord 
     568    @property 
     569    def par_offset(self): return self._par_offset 
     570    @property 
     571    def num_coord(self): return self.details[-2] 
     572    @num_coord.setter 
     573    def num_coord(self, v): self.details[-2] = v 
     574    @property 
     575    def total_pd(self): return self.details[-3] 
     576    @total_pd.setter 
     577    def total_pd(self, v): self.details[-3] = v 
     578    @property 
     579    def num_active(self): return self.details[-4] 
     580    @num_active.setter 
     581    def num_active(self, v): self.details[-4] = v 
     582 
     583    def show(self): 
     584        print("total_pd", self.total_pd) 
     585        print("num_active", self.num_active) 
     586        print("pd_par", self.pd_par) 
     587        print("pd_length", self.pd_length) 
     588        print("pd_offset", self.pd_offset) 
     589        print("pd_stride", self.pd_stride) 
     590        print("par_offsets", self.par_offset) 
     591        print("num_coord", self.num_coord) 
     592        print("par_coord", self.par_coord) 
     593        print("pd_coord", self.pd_coord) 
     594        print("theta par", self.details[-1]) 
     595 
     596def mono_details(model_info): 
     597    details = CoordinationDetails(model_info) 
     598    # The zero defaults for monodisperse systems are mostly fine 
     599    details.par_offset[:] = np.arange(2, len(details.par_offset)+2) 
     600    return details 
     601 
     602def poly_details(model_info, weights): 
     603    #print("weights",weights) 
     604    weights = weights[2:] # Skip scale and background 
     605 
     606    # Decreasing list of polydispersity lengths 
     607    # Note: the reversing view, x[::-1], does not require a copy 
     608    pd_length = np.array([len(w) for w in weights]) 
     609    num_active = np.sum(pd_length>1) 
     610    if num_active > model_info['parameters'].max_pd: 
     611        raise ValueError("Too many polydisperse parameters") 
     612 
     613    pd_offset = np.cumsum(np.hstack((0, pd_length))) 
     614    idx = np.argsort(pd_length)[::-1][:num_active] 
     615    par_length = np.array([max(len(w),1) for w in weights]) 
     616    pd_stride = np.cumprod(np.hstack((1, par_length[idx]))) 
     617    par_offsets = np.cumsum(np.hstack((2, par_length))) 
     618 
     619    details = CoordinationDetails(model_info) 
     620    details.pd_par[:num_active] = idx 
     621    details.pd_length[:num_active] = pd_length[idx] 
     622    details.pd_offset[:num_active] = pd_offset[idx] 
     623    details.pd_stride[:num_active] = pd_stride[:-1] 
     624    details.par_offset[:] = par_offsets[:-1] 
     625    details.total_pd = pd_stride[-1] 
     626    details.num_active = num_active 
     627    # Without constraints coordinated parameters are just the pd parameters 
     628    details.par_coord[:num_active] = idx 
     629    details.pd_coord[:num_active] = 2**np.arange(num_active) 
     630    details.num_coord = num_active 
     631    #details.show() 
     632    return details 
     633 
     634def constrained_poly_details(model_info, weights, constraints): 
     635    # Need to find the independently varying pars and sort them 
     636    # Need to build a coordination list for the dependent variables 
     637    # Need to generate a constraints function which takes values 
     638    # and weights, returning par blocks 
     639    raise NotImplementedError("Can't handle constraints yet") 
     640 
     641 
     642def create_default_functions(model_info): 
     643    """ 
     644    Autogenerate missing functions, such as Iqxy from Iq. 
     645 
     646    This only works for Iqxy when Iq is written in python. :func:`make_source` 
     647    performs a similar role for Iq written in C. 
     648    """ 
     649    if callable(model_info['Iq']) and model_info['Iqxy'] is None: 
     650        partable = model_info['parameters'] 
     651        if partable.has_2d: 
     652            raise ValueError("Iqxy model is missing") 
     653        Iq = model_info['Iq'] 
     654        def Iqxy(qx, qy, **kw): 
     655            return Iq(np.sqrt(qx**2 + qy**2), **kw) 
     656        model_info['Iqxy'] = Iqxy 
     657 
     658 
     659def load_kernel_module(model_name): 
     660    if model_name.endswith('.py'): 
     661        kernel_module = load_custom_kernel_module(model_name) 
     662    else: 
     663        from sasmodels import models 
     664        __import__('sasmodels.models.'+model_name) 
     665        kernel_module = getattr(models, model_name, None) 
     666    return kernel_module 
     667 
    646668 
    647669def make_model_info(kernel_module): 
    … …  
    670692      model variants (e.g., the list of cases) or is None if there are no 
    671693      model variants 
    672     * *defaults* is the *{parameter: value}* table built from the parameter 
    673       description table. 
    674     * *limits* is the *{parameter: [min, max]}* table built from the 
    675       parameter description table. 
    676     * *partypes* categorizes the model parameters. See 
     694    * *par_type* categorizes the model parameters. See 
    677695      :func:`categorize_parameters` for details. 
    678696    * *demo* contains the *{parameter: value}* map used in compare (and maybe 
    … …  
    688706      *model_info* blocks for the composition objects.  This allows us to 
    689707      build complete product and mixture models from just the info. 
    690  
    691708    """ 
    692709    # TODO: maybe turn model_info into a class ModelDefinition 
    693     parameters = COMMON_PARAMETERS + kernel_module.parameters 
     710    #print("make parameter table", kernel_module.parameters) 
     711    parameters = make_parameter_table(kernel_module.parameters) 
    694712    filename = abspath(kernel_module.__file__) 
    695713    kernel_id = splitext(basename(filename))[0] 
    … …  
    701719        filename=abspath(kernel_module.__file__), 
    702720        name=name, 
    703         title=kernel_module.title, 
    704         description=kernel_module.description, 
     721        title=getattr(kernel_module, 'title', name+" model"), 
     722        description=getattr(kernel_module, 'description', 'no description'), 
    705723        parameters=parameters, 
    706724        composition=None, 
    … …  
    709727        single=getattr(kernel_module, 'single', True), 
    710728        structure_factor=getattr(kernel_module, 'structure_factor', False), 
     729        profile_axes=getattr(kernel_module, 'profile_axes', ['x','y']), 
    711730        variant_info=getattr(kernel_module, 'invariant_info', None), 
    712731        demo=getattr(kernel_module, 'demo', None), 
    … …  
    714733        tests=getattr(kernel_module, 'tests', []), 
    715734        ) 
    716     process_parameters(model_info) 
     735    set_demo(model_info, getattr(kernel_module, 'demo', None)) 
    717736    # Check for optional functions 
    718     functions = "ER VR form_volume Iq Iqxy shape sesans".split() 
     737    functions = "ER VR form_volume Iq Iqxy profile sesans".split() 
    719738    model_info.update((k, getattr(kernel_module, k, None)) for k in functions) 
     739    create_default_functions(model_info) 
     740    # Precalculate the monodisperse parameters 
     741    # TODO: make this a lazy evaluator 
     742    # make_model_info is called for every model on sasview startup 
     743    model_info['mono_details'] = mono_details(model_info) 
    720744    return model_info 
    721745 
    … …  
    769793 
    770794 
    771  
    772795def demo_time(): 
    773796    """ 
    … …  
    785808    Program which prints the source produced by the model. 
    786809    """ 
     810    import sys 
    787811    if len(sys.argv) <= 1: 
    788812        print("usage: python -m sasmodels.generate modelname") 
    789813    else: 
    790814        name = sys.argv[1] 
    791         import sasmodels.models 
    792         __import__('sasmodels.models.' + name) 
    793         model = getattr(sasmodels.models, name) 
    794         model_info = make_model_info(model) 
     815        kernel_module = load_kernel_module(name) 
     816        model_info = make_model_info(kernel_module) 
    795817        source = make_source(model_info) 
    796818        print(source) 

  • sasmodels/kernelcl.py

    r8e0d974 rba32cdd  
    4848harmless, albeit annoying. 
    4949""" 
     50from __future__ import print_function 
    5051import os 
    5152import warnings 
    … …  
    7374# of polydisperse parameters. 
    7475MAX_LOOPS = 2048 
     76 
     77 
     78# Pragmas for enable OpenCL features.  Be sure to protect them so that they 
     79# still compile even if OpenCL is not present. 
     80_F16_PRAGMA = """\ 
     81#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 
     82#  pragma OPENCL EXTENSION cl_khr_fp16: enable 
     83#endif 
     84""" 
     85 
     86_F64_PRAGMA = """\ 
     87#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 
     88#  pragma OPENCL EXTENSION cl_khr_fp64: enable 
     89#endif 
     90""" 
    7591 
    7692 
    … …  
    142158        raise RuntimeError("%s not supported for devices"%dtype) 
    143159 
    144     source = generate.convert_type(source, dtype) 
     160    source_list = [generate.convert_type(source, dtype)] 
     161 
     162    if dtype == generate.F16: 
     163        source_list.insert(0, _F16_PRAGMA) 
     164    elif dtype == generate.F64: 
     165        source_list.insert(0, _F64_PRAGMA) 
     166 
    145167    # Note: USE_SINCOS makes the intel cpu slower under opencl 
    146168    if context.devices[0].type == cl.device_type.GPU: 
    147         source = "#define USE_SINCOS\n" + source 
     169        source_list.insert(0, "#define USE_SINCOS\n") 
    148170    options = (get_fast_inaccurate_build_options(context.devices[0]) 
    149171               if fast else []) 
     172    source = "\n".join(source_list) 
    150173    program = cl.Program(context, source).build(options=options) 
    151174    return program 
    … …  
    220243        key = "%s-%s-%s"%(name, dtype, fast) 
    221244        if key not in self.compiled: 
    222             #print("compiling",name) 
     245            print("compiling",name) 
    223246            dtype = np.dtype(dtype) 
    224             program = compile_model(self.get_context(dtype), source, dtype, fast) 
     247            program = compile_model(self.get_context(dtype), 
     248                                    str(source), dtype, fast) 
    225249            self.compiled[key] = program 
    226250        return self.compiled[key] 
    … …  
    314338        self.program = None 
    315339 
    316     def __call__(self, q_vectors): 
     340    def make_kernel(self, q_vectors): 
    317341        if self.program is None: 
    318342            compiler = environment().compile_program 
    319             self.program = compiler(self.info['name'], self.source, self.dtype, 
    320                                     self.fast) 
     343            self.program = compiler(self.info['name'], self.source, 
     344                                    self.dtype, self.fast) 
    321345        is_2d = len(q_vectors) == 2 
    322346        kernel_name = generate.kernel_name(self.info, is_2d) 
    … …  
    365389        # at this point, so instead using 32, which is good on the set of 
    366390        # architectures tested so far. 
    367         self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 
     391        if self.is_2d: 
     392            # Note: 17 rather than 15 because results is 2 elements 
     393            # longer than input. 
     394            width = ((self.nq+17)//16)*16 
     395            self.q = np.empty((width, 2), dtype=dtype) 
     396            self.q[:self.nq, 0] = q_vectors[0] 
     397            self.q[:self.nq, 1] = q_vectors[1] 
     398        else: 
     399            # Note: 33 rather than 31 because results is 2 elements 
     400            # longer than input. 
     401            width = ((self.nq+33)//32)*32 
     402            self.q = np.empty(width, dtype=dtype) 
     403            self.q[:self.nq] = q_vectors[0] 
     404        self.global_size = [self.q.shape[0]] 
    368405        context = env.get_context(self.dtype) 
    369         self.global_size = [self.q_vectors[0].size] 
    370406        #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         ] 
     407        self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     408                             hostbuf=self.q) 
    375409 
    376410    def release(self): 
    … …  
    378412        Free the memory. 
    379413        """ 
    380         for b in self.q_buffers: 
    381             b.release() 
    382         self.q_buffers = [] 
     414        if self.q is not None: 
     415            self.q.release() 
     416            self.q = None 
    383417 
    384418    def __del__(self): 
    … …  
    406440    """ 
    407441    def __init__(self, kernel, model_info, q_vectors, dtype): 
     442        max_pd = model_info['max_pd'] 
     443        npars = len(model_info['parameters'])-2 
    408444        q_input = GpuInput(q_vectors, dtype) 
     445        self.dtype = dtype 
     446        self.dim = '2d' if q_input.is_2d else '1d' 
    409447        self.kernel = kernel 
    410448        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] 
     449        self.pd_stop_index = 4*max_pd-1 
     450        # plus three for the normalization values 
     451        self.result = np.empty(q_input.nq+3, q_input.dtype) 
    415452 
    416453        # Inputs and outputs for each kernel call 
    … …  
    418455        env = environment() 
    419456        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, 
     457 
     458        # details is int32 data, padded to an 8 integer boundary 
     459        size = ((max_pd*5 + npars*3 + 2 + 7)//8)*8 
     460        self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
    423461                               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): 
     462        self.q_input = q_input # allocated by GpuInput above 
     463 
     464        self._need_release = [ self.result_b, self.q_input ] 
     465 
     466    def __call__(self, details, weights, values, cutoff): 
    429467        real = (np.float32 if self.q_input.dtype == generate.F32 
    430468                else np.float64 if self.q_input.dtype == generate.F64 
    431469                else np.float16 if self.q_input.dtype == generate.F16 
    432470                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 
     471        assert details.dtype == np.int32 
     472        assert weights.dtype == real and values.dtype == real 
     473 
     474        context = self.queue.context 
     475        details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     476                              hostbuf=details) 
     477        weights_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     478                              hostbuf=weights) 
     479        values_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     480                             hostbuf=values) 
     481 
     482        start, stop = 0, self.details[self.pd_stop_index] 
     483        args = [ 
     484            np.uint32(self.q_input.nq), np.uint32(start), np.uint32(stop), 
     485            self.details_b, self.weights_b, self.values_b, 
     486            self.q_input.q_b, self.result_b, real(cutoff), 
     487        ] 
    460488        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 
     489        cl.enqueue_copy(self.queue, self.result, self.result_b) 
     490        [v.release() for v in details_b, weights_b, values_b] 
     491 
     492        return self.result[:self.nq] 
    464493 
    465494    def release(self): 

  • sasmodels/kerneldll.py

    r6ad0e87 rd19962c  
    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 
    … …  
    8180        COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 
    8281        if "SAS_OPENMP" in os.environ: 
    83             COMPILE = COMPILE + " -fopenmp" 
     82            COMPILE += " -fopenmp" 
    8483else: 
    8584    COMPILE = "cc -shared -fPIC -fopenmp -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 
    … …  
    9089 
    9190 
    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  
     91def dll_name(model_info, dtype): 
     92    """ 
     93    Name of the dll containing the model.  This is the base file name without 
     94    any path or extension, with a form such as 'sas_sphere32'. 
     95    """ 
     96    bits = 8*dtype.itemsize 
     97    return "sas_%s%d"%(model_info['id'], bits) 
     98 
     99def dll_path(model_info, dtype): 
     100    """ 
     101    Complete path to the dll for the model.  Note that the dll may not 
     102    exist yet if it hasn't been compiled. 
     103    """ 
     104    return os.path.join(DLL_PATH, dll_name(model_info, dtype)+".so") 
    106105 
    107106def make_dll(source, model_info, dtype="double"): 
    … …  
    133132        dtype = generate.F64  # Force 64-bit dll 
    134133 
    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   
    142134    source = generate.convert_type(source, dtype) 
    143     source_files = generate.model_sources(model_info) + [model_info['filename']] 
     135    newest = generate.timestamp(model_info) 
    144136    dll = dll_path(model_info, dtype) 
    145     newest = max(os.path.getmtime(f) for f in source_files) 
    146137    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) 
     138        basename = dll_name(model_info, dtype) + "_" 
     139        fid, filename = tempfile.mkstemp(suffix=".c", prefix=basename) 
    149140        os.fdopen(fid, "w").write(source) 
    150141        command = COMPILE%{"source":filename, "output":dll} 
    … …  
    171162    return DllModel(filename, model_info, dtype=dtype) 
    172163 
    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  
    177164class DllModel(object): 
    178165    """ 
    … …  
    197184 
    198185    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  
    204186        #print("dll", self.dllpath) 
    205187        try: 
    … …  
    212194              else c_double if self.dtype == generate.F64 
    213195              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 [] 
     196 
     197        # int, int, int, int*, double*, double*, double*, double*, double*, double 
     198        argtypes = [c_int32]*3 + [c_void_p]*5 + [fp] 
    216199        self.Iq = self.dll[generate.kernel_name(self.info, False)] 
    217         self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d 
    218  
    219200        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() 
     201        self.Iq.argtypes = argtypes 
     202        self.Iqxy.argtypes = argtypes 
    223203 
    224204    def __getstate__(self): 
    … …  
    229209        self.dll = None 
    230210 
    231     def __call__(self, q_vectors): 
     211    def make_kernel(self, q_vectors): 
    232212        q_input = PyInput(q_vectors, self.dtype) 
    233213        if self.dll is None: self._load_dll() 
    … …  
    272252    """ 
    273253    def __init__(self, kernel, model_info, q_input): 
     254        self.kernel = kernel 
    274255        self.info = model_info 
    275256        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): 
     257        self.dtype = q_input.dtype 
     258        self.dim = '2d' if q_input.is_2d else '1d' 
     259        self.result = np.empty(q_input.nq+3, q_input.dtype) 
     260 
     261    def __call__(self, details, weights, values, cutoff): 
    286262        real = (np.float32 if self.q_input.dtype == generate.F32 
    287263                else np.float64 if self.q_input.dtype == generate.F64 
    288264                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) 
     265        assert isinstance(details, generate.CoordinationDetails) 
     266        assert weights.dtype == real and values.dtype == real 
     267 
     268        start, stop = 0, details.total_pd 
     269        #print("in kerneldll") 
     270        #print("weights", weights) 
     271        #print("values", values) 
     272        args = [ 
     273            self.q_input.nq, # nq 
     274            start, # pd_start 
     275            stop, # pd_stop pd_stride[MAX_PD] 
     276            details.ctypes.data, # problem 
     277            weights.ctypes.data,  # weights 
     278            values.ctypes.data,  #pars 
     279            self.q_input.q.ctypes.data, #q 
     280            self.result.ctypes.data,   # results 
     281            real(cutoff), # cutoff 
     282            ] 
    303283        self.kernel(*args) 
    304  
    305         return self.res 
     284        return self.result[:-3] 
    306285 
    307286    def release(self): 

  • sasmodels/kernelpy.py

    ra84a0ca r48fbd50  
    5353        self.dtype = dtype 
    5454        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] 
     55        if self.is_2d: 
     56            self.q = np.empty((self.nq, 2), dtype=dtype) 
     57            self.q[:, 0] = q_vectors[0] 
     58            self.q[:, 1] = q_vectors[1] 
     59        else: 
     60            self.q = np.empty(self.nq, dtype=dtype) 
     61            self.q[:self.nq] = q_vectors[0] 
    5762 
    5863    def release(self): 
    … …  
    6065        Free resources associated with the model inputs. 
    6166        """ 
    62         self.q_vectors = [] 
     67        self.q = None 
    6368 
    6469class PyKernel(object): 

  • sasmodels/mixture.py

    r72a081d r69aa451  
    1414import numpy as np 
    1515 
    16 from .generate import process_parameters, COMMON_PARAMETERS, Parameter 
     16from .modelinfo import Parameter, ParameterTable 
    1717 
    1818SCALE=0 
    … …  
    3434 
    3535    # Build new parameter list 
    36     pars = COMMON_PARAMETERS + [] 
     36    pars = [] 
    3737    for k, part in enumerate(parts): 
    3838        # Parameter prefix per model, A_, B_, ... 
     39        # Note that prefix must also be applied to id and length_control 
     40        # to support vector parameters 
    3941        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_, ... 
     42        pars.append(Parameter(prefix+'scale')) 
     43        for p in part['parameters'].kernel_pars: 
     44            p = copy(p) 
     45            p.name = prefix+p.name 
     46            p.id = prefix+p.id 
     47            if p.length_control is not None: 
     48                p.length_control = prefix+p.length_control 
    5049            pars.append(p) 
     50    partable = ParameterTable(pars) 
    5151 
    5252    model_info = {} 
    … …  
    5858    model_info['docs'] = model_info['title'] 
    5959    model_info['category'] = "custom" 
    60     model_info['parameters'] = pars 
     60    model_info['parameters'] = partable 
    6161    #model_info['single'] = any(part['single'] for part in parts) 
    6262    model_info['structure_factor'] = False 
    … …  
    6767    # Remember the component info blocks so we can build the model 
    6868    model_info['composition'] = ('mixture', parts) 
    69     process_parameters(model_info) 
    70     return model_info 
    7169 
    7270 

  • sasmodels/model_test.py

    ra84a0ca r48fbd50  
    5151 
    5252from .core import list_models, load_model_info, build_model, HAVE_OPENCL 
    53 from .core import make_kernel, call_kernel, call_ER, call_VR 
     53from .core import call_kernel, call_ER, call_VR 
    5454from .exception import annotate_exception 
    5555 
    … …  
    187187                Qx, Qy = zip(*x) 
    188188                q_vectors = [np.array(Qx), np.array(Qy)] 
    189                 kernel = make_kernel(model, q_vectors) 
     189                kernel = model.make_kernel(q_vectors) 
    190190                actual = call_kernel(kernel, pars) 
    191191            else: 
    192192                q_vectors = [np.array(x)] 
    193                 kernel = make_kernel(model, q_vectors) 
     193                kernel = model.make_kernel(q_vectors) 
    194194                actual = call_kernel(kernel, pars) 
    195195 

  • sasmodels/models/cylinder.c

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

  • sasmodels/models/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/lib/Si.c

    re7678b2 rba32cdd  
     1// integral of sin(x)/x Taylor series approximated to w/i 0.1% 
    12double Si(double x); 
    23 
    3 // integral of sin(x)/x Taylor series approximated to w/i 0.1% 
    44double Si(double x) 
    55{ 

  • sasmodels/models/lib/polevl.c

    r0b05c24 rba32cdd  
    5151*/ 
    5252 
    53 double polevl( double x, constant double *coef, int N ) { 
     53double polevl( double x, constant double *coef, int N ); 
     54double p1evl( double x, constant double *coef, int N ); 
     55 
     56 
     57double polevl( double x, constant double *coef, int N ) 
     58{ 
    5459 
    5560    int i = 0; 
    … …  
    7075 */ 
    7176 
    72 double p1evl( double x, constant double *coef, int N ) { 
     77double p1evl( double x, constant double *coef, int N ) 
     78{ 
    7379    int i=0; 
    7480    double ans = x+coef[i]; 

  • sasmodels/models/lib/sas_J0.c

    ra776125 rba32cdd  
    4949Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier 
    5050*/ 
     51 
     52double sas_J0(double x); 
    5153 
    5254/* Note: all coefficients satisfy the relative error criterion 
    … …  
    177179 }; 
    178180 
    179 double sas_J0(double x) { 
     181double sas_J0(double x) 
     182{ 
    180183 
    181184//Cephes single precission 

  • sasmodels/models/lib/sas_J1.c

    r19b9005 rba32cdd  
    3939Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier 
    4040*/ 
     41double sas_J1(double x); 
     42double sas_J1c(double x); 
    4143 
    4244constant double RPJ1[8] = { 
    … …  
    135137    }; 
    136138 
    137 double sas_J1(double x) { 
     139double sas_J1(double x) 
     140{ 
    138141 
    139142//Cephes double pression function 

  • sasmodels/models/lib/sas_JN.c

    ra776125 rba32cdd  
    4747Copyright 1984, 1987, 2000 by Stephen L. Moshier 
    4848*/ 
    49  
    5049double sas_JN( int n, double x ); 
    5150 
    52 double sas_JN( int n, double x ) { 
     51double sas_JN( int n, double x ) 
     52{ 
    5353 
    5454    const double MACHEP = 1.11022302462515654042E-16; 

  • sasmodels/models/lib/sph_j1c.c

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

  • sasmodels/models/lib/sphere_form.c

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

  • sasmodels/models/lib/wrc_cyl.c

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

  • sasmodels/models/onion.c

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

  • sasmodels/models/onion.py

    rec45c4f rec45c4f  
    293293 
    294294#             ["name", "units", default, [lower, upper], "type","description"], 
    295 parameters = [["core_sld", "1e-6/Ang^2", 1.0, [-inf, inf], "", 
     295parameters = [["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "", 
    296296               "Core scattering length density"], 
    297297              ["core_radius", "Ang", 200., [0, inf], "volume", 
    298298               "Radius of the core"], 
    299               ["solvent_sld", "1e-6/Ang^2", 6.4, [-inf, inf], "", 
     299              ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "", 
    300300               "Solvent scattering length density"], 
    301301              ["n", "", 1, [0, 10], "volume", 
    302302               "number of shells"], 
    303               ["in_sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "", 
     303              ["sld_in[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "", 
    304304               "scattering length density at the inner radius of shell k"], 
    305               ["out_sld[n]", "1e-6/Ang^2", 2.0, [-inf, inf], "", 
     305              ["sld_out[n]", "1e-6/Ang^2", 2.0, [-inf, inf], "", 
    306306               "scattering length density at the outer radius of shell k"], 
    307307              ["thickness[n]", "Ang", 40., [0, inf], "volume", 
    … …  
    311311              ] 
    312312 
    313 #source = ["lib/sph_j1c.c", "onion.c"] 
    314  
    315 def Iq(q, *args, **kw): 
    316     return q 
    317  
    318 def Iqxy(qx, *args, **kw): 
    319     return qx 
    320  
    321  
    322 def shape(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A): 
     313source = ["lib/sph_j1c.c", "onion.c"] 
     314 
     315#def Iq(q, *args, **kw): 
     316#    return q 
     317 
     318profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 
     319def profile(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A): 
    323320    """ 
    324321    SLD profile 
    … …  
    374371 
    375372demo = { 
    376     "solvent_sld": 2.2, 
    377     "core_sld": 1.0, 
     373    "sld_solvent": 2.2, 
     374    "sld_core": 1.0, 
    378375    "core_radius": 100, 
    379376    "n": 4, 
    380     "in_sld": [0.5, 1.5, 0.9, 2.0], 
    381     "out_sld": [nan, 0.9, 1.2, 1.6], 
     377    "sld_in": [0.5, 1.5, 0.9, 2.0], 
     378    "sld_out": [nan, 0.9, 1.2, 1.6], 
    382379    "thickness": [50, 75, 150, 75], 
    383380    "A": [0, -1, 1e-4, 1], 

  • sasmodels/models/rpa.py

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

  • sasmodels/product.py

    rf247314 rf247314  
    1414 
    1515from .core import call_ER_VR 
    16 from .generate import process_parameters 
    1716 
    1817SCALE=0 
    … …  
    5857    # Iq, Iqxy, form_volume, ER, VR and sesans 
    5958    model_info['composition'] = ('product', [p_info, s_info]) 
    60     process_parameters(model_info) 
    6159    return model_info 
    6260 
    … …  
    9391        # a parameter map. 
    9492        par_map = {} 
    95         p_info = p_kernel.info['partype'] 
    96         s_info = s_kernel.info['partype'] 
     93        p_info = p_kernel.info['par_type'] 
     94        s_info = s_kernel.info['par_type'] 
    9795        vol_pars = set(p_info['volume']) 
    9896        if dim == '2d': 

  • sasmodels/resolution.py

    r486fcf6 r486fcf6  
    477477    """ 
    478478    from sasmodels import core 
    479     kernel = core.make_kernel(form, [q]) 
     479    kernel = form.make_kernel([q]) 
    480480    theory = core.call_kernel(kernel, pars) 
    481481    kernel.release() 
    … …  
    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 
     504    par_set = set([p.name for p in form.info['parameters']]) 
     505    if any(k not in par_set for k in pars.keys()): 
     506        extra = set(pars.keys()) - par_set 
    507507        raise ValueError("bad parameters: [%s] not in [%s]"% 
    508508                         (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 
    … …  
    556556    from scipy.integrate import romberg 
    557557 
    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 
     558    par_set = set([p.name for p in form.info['parameters']]) 
     559    if any(k not in par_set for k in pars.keys()): 
     560        extra = set(pars.keys()) - par_set 
    561561        raise ValueError("bad parameters: [%s] not in [%s]"% 
    562                          (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 
     562                         (", ".join(sorted(extra)), 
     563                          ", ".join(sorted(pars.keys())))) 
    563564 
    564565    _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq) 
    … …  
    693694    def _eval_sphere(self, pars, resolution): 
    694695        from sasmodels import core 
    695         kernel = core.make_kernel(self.model, [resolution.q_calc]) 
     696        kernel = self.model.make_kernel([resolution.q_calc]) 
    696697        theory = core.call_kernel(kernel, pars) 
    697698        result = resolution.apply(theory) 
    … …  
    10621063    model = core.build_model(model_info) 
    10631064 
    1064     kernel = core.make_kernel(model, [resolution.q_calc]) 
     1065    kernel = model.make_kernel([resolution.q_calc]) 
    10651066    theory = core.call_kernel(kernel, pars) 
    10661067    Iq = resolution.apply(theory) 

  • sasmodels/sasview_model.py

    rf247314 rf247314  
    2121 
    2222from . import core 
    23 from . import generate 
     23from . import weights 
    2424 
    2525def standard_models(): 
    … …  
    5252    """ 
    5353    def __init__(self): 
    54         self._kernel = None 
     54        self._model = None 
    5555        model_info = self._model_info 
     56        parameters = model_info['parameters'] 
    5657 
    5758        self.name = model_info['name'] 
    5859        self.description = model_info['description'] 
    5960        self.category = None 
    60         self.multiplicity_info = None 
    61         self.is_multifunc = False 
     61        #self.is_multifunc = False 
     62        for p in parameters.kernel_parameters: 
     63            if p.is_control: 
     64                profile_axes = model_info['profile_axes'] 
     65                self.multiplicity_info = [ 
     66                    p.limits[1], p.name, p.choices, profile_axes[0] 
     67                    ] 
     68                break 
     69        else: 
     70            self.multiplicity_info = [] 
    6271 
    6372        ## interpret the parameters 
    … …  
    6675        self.params = collections.OrderedDict() 
    6776        self.dispersion = dict() 
    68         partype = model_info['partype'] 
    69  
    70         for p in model_info['parameters']: 
     77 
     78        self.orientation_params = [] 
     79        self.magnetic_params = [] 
     80        self.fixed = [] 
     81        for p in parameters.user_parameters(): 
    7182            self.params[p.name] = p.default 
    7283            self.details[p.name] = [p.units] + p.limits 
    73  
    74         for name in partype['pd-2d']: 
    75             self.dispersion[name] = { 
    76                 'width': 0, 
    77                 'npts': 35, 
    78                 'nsigmas': 3, 
    79                 'type': 'gaussian', 
    80             } 
    81  
    82         self.orientation_params = ( 
    83             partype['orientation'] 
    84             + [n + '.width' for n in partype['orientation']] 
    85             + partype['magnetic']) 
    86         self.magnetic_params = partype['magnetic'] 
    87         self.fixed = [n + '.width' for n in partype['pd-2d']] 
     84            if p.polydisperse: 
     85                self.dispersion[p.name] = { 
     86                    'width': 0, 
     87                    'npts': 35, 
     88                    'nsigmas': 3, 
     89                    'type': 'gaussian', 
     90                } 
     91            if p.type == 'orientation': 
     92                self.orientation_params.append(p.name) 
     93                self.orientation_params.append(p.name+".width") 
     94                self.fixed.append(p.name+".width") 
     95            if p.type == 'magnetic': 
     96                self.orientation_params.append(p.name) 
     97                self.magnetic_params.append(p.name) 
     98                self.fixed.append(p.name+".width") 
     99 
    88100        self.non_fittable = [] 
    89101 
    … …  
    104116    def __get_state__(self): 
    105117        state = self.__dict__.copy() 
    106         model_id = self._model_info['id'] 
    107118        state.pop('_kernel') 
    108119        # May need to reload model info on set state since it has pointers 
    … …  
    113124    def __set_state__(self, state): 
    114125        self.__dict__ = state 
    115         self._kernel = None 
     126        self._model = None 
    116127 
    117128    def __str__(self): 
    … …  
    202213    def getDispParamList(self): 
    203214        """ 
    204         Return a list of all available parameters for the model 
     215        Return a list of polydispersity parameters for the model 
    205216        """ 
    206217        # TODO: fix test so that parameter order doesn't matter 
    207         ret = ['%s.%s' % (d.lower(), p) 
    208                for d in self._model_info['partype']['pd-2d'] 
    209                for p in ('npts', 'nsigmas', 'width')] 
     218        ret = ['%s.%s' % (p.name.lower(), ext) 
     219               for p in self._model_info['parameters'].user_parameters() 
     220               for ext in ('npts', 'nsigmas', 'width') 
     221               if p.polydisperse] 
    210222        #print(ret) 
    211223        return ret 
    … …  
    280292            # Check whether we have a list of ndarrays [qx,qy] 
    281293            qx, qy = qdist 
    282             partype = self._model_info['partype'] 
    283             if not partype['orientation'] and not partype['magnetic']: 
     294            if not self._model_info['parameters'].has_2d: 
    284295                return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2)) 
    285296            else: 
    … …  
    303314        to the card for each evaluation. 
    304315        """ 
    305         if self._kernel is None: 
    306             self._kernel = core.build_model(self._model_info) 
     316        if self._model is None: 
     317            self._model = core.build_model(self._model_info, platform='dll') 
    307318        q_vectors = [np.asarray(q) for q in args] 
    308         fn = self._kernel(q_vectors) 
    309         pars = [self.params[v] for v in fn.fixed_pars] 
    310         pd_pars = [self._get_weights(p) for p in fn.pd_pars] 
    311         result = fn(pars, pd_pars, self.cutoff) 
    312         fn.q_input.release() 
    313         fn.release() 
     319        kernel = self._model.make_kernel(q_vectors) 
     320        pairs = [self._get_weights(p) 
     321                 for p in self._model_info['parameters'].call_parameters] 
     322        details, weights, values = core.build_details(kernel, pairs) 
     323        return kernel(details, weights, values, cutoff=self.cutoff) 
     324        kernel.q_input.release() 
     325        kernel.release() 
    314326        return result 
    315327 
    … …  
    384396    def _get_weights(self, par): 
    385397        """ 
    386             Return dispersion weights 
    387             :param par parameter name 
    388         """ 
    389         from . import weights 
    390  
    391         relative = self._model_info['partype']['pd-rel'] 
    392         limits = self._model_info['limits'] 
    393         dis = self.dispersion[par] 
    394         value, weight = weights.get_weights( 
    395             dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 
    396             self.params[par], limits[par], par in relative) 
    397         return value, weight / np.sum(weight) 
    398  
     398        Return dispersion weights for parameter 
     399        """ 
     400        if par.polydisperse: 
     401            dis = self.dispersion[par.name] 
     402            value, weight = weights.get_weights( 
     403                dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 
     404                self.params[par.name], par.limits, par.relative_pd) 
     405            return value, weight / np.sum(weight) 
     406        else: 
     407            return [self.params[par.name]], [] 
     408 
     409def test_model(): 
     410    Cylinder = make_class('cylinder') 
     411    cylinder = Cylinder() 
     412    return cylinder.evalDistribution([0.1,0.1]) 
     413 
     414if __name__ == "__main__": 
     415    print("cylinder(0.1,0.1)=%g"%test_model()) 
Note: See TracChangeset for help on using the changeset viewer.