Changes in / [e9b1663d:34edbb8] in sasmodels


Ignore:
Files:
6 deleted
16 edited

Legend:

Unmodified
Added
Removed

  • doc/developer/index.rst

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

  • sasmodels/bumps_model.py

    r9217ef8 r9217ef8  
    8181    from bumps.names import Parameter 
    8282 
    83     pars = {} # => floating point parameters 
     83    pars = {} 
    8484    for p in model_info['parameters']: 
    8585        value = kwargs.pop(p.name, p.default) 
    8686        pars[p.name] = Parameter.default(value, name=p.name, limits=p.limits) 
    87     for name in model_info['par_type']['pd']: 
     87    for name in model_info['partype']['pd-2d']: 
    8888        for xpart, xdefault, xlimits in [ 
    8989                ('_pd', 0., pars[name].limits), 
    … …  
    9595            pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits) 
    9696 
    97     pd_types = {}  # => distribution names 
    98     for name in model_info['par_type']['pd']: 
     97    pd_types = {} 
     98    for name in model_info['partype']['pd-2d']: 
    9999        xname = name + '_pd_type' 
    100100        xvalue = kwargs.pop(xname, 'gaussian') 

  • sasmodels/compare.py

    r69aa451 r1671636  
    309309        exclude = lambda n: False 
    310310    else: 
    311         par1d = model_info['par_type']['1d'] 
     311        partype = model_info['partype'] 
     312        par1d = set(partype['fixed-1d']+partype['pd-1d']) 
    312313        exclude = lambda n: n not in par1d 
    313314    lines = [] 
    … …  
    453454    """ 
    454455    # initialize the code so time is more accurate 
    455     if Nevals > 1: 
    456         value = calculator(**suppress_pd(pars)) 
     456    value = calculator(**suppress_pd(pars)) 
    457457    toc = tic() 
    458458    for _ in range(max(Nevals, 1)):  # make sure there is at least one eval 
    … …  
    877877        pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars']) 
    878878        if not opts['is2d']: 
    879             for p in model_info['parameters'].type['1d']: 
    880                 for ext in ['', '_pd', '_pd_n', '_pd_nsigma']: 
    881                     k = p.name+ext 
    882                     v = pars.get(k, None) 
    883                     if v is not None: 
    884                         v.range(*parameter_range(k, v.value)) 
     879            active = [base + ext 
     880                      for base in model_info['partype']['pd-1d'] 
     881                      for ext in ['', '_pd', '_pd_n', '_pd_nsigma']] 
     882            active.extend(model_info['partype']['fixed-1d']) 
     883            for k in active: 
     884                v = pars[k] 
     885                v.range(*parameter_range(k, v.value)) 
    885886        else: 
    886887            for k, v in pars.items(): 

  • sasmodels/core.py

    r69aa451 re79f0a5  
    2626__all__ = [ 
    2727    "list_models", "load_model_info", "precompile_dll", 
    28     "build_model", "call_kernel", "call_ER_VR", 
     28    "build_model", "make_kernel", "call_kernel", "call_ER_VR", 
    2929] 
    3030 
    … …  
    167167 
    168168 
    169 def get_weights(parameter, values): 
     169def make_kernel(model, q_vectors): 
     170    """ 
     171    Return a computation kernel from the model definition and the q input. 
     172    """ 
     173    return model(q_vectors) 
     174 
     175def get_weights(model_info, pars, name): 
    170176    """ 
    171177    Generate the distribution for parameter *name* given the parameter values 
    … …  
    175181    from the *pars* dictionary for parameter value and parameter dispersion. 
    176182    """ 
    177     value = values.get(parameter.name, parameter.default) 
    178     if parameter.type not in ('volume', 'orientation'): 
    179         return np.array([value]), np.array([1.0]) 
    180     relative = parameter.type == 'volume' 
    181     limits = parameter.limits 
    182     disperser = values.get(parameter.name+'_pd_type', 'gaussian') 
    183     npts = values.get(parameter.name+'_pd_n', 0) 
    184     width = values.get(parameter.name+'_pd', 0.0) 
    185     nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 
     183    relative = name in model_info['partype']['pd-rel'] 
     184    limits = model_info['limits'][name] 
     185    disperser = pars.get(name+'_pd_type', 'gaussian') 
     186    value = pars.get(name, model_info['defaults'][name]) 
     187    npts = pars.get(name+'_pd_n', 0) 
     188    width = pars.get(name+'_pd', 0.0) 
     189    nsigma = pars.get(name+'_pd_nsigma', 3.0) 
    186190    value, weight = weights.get_weights( 
    187191        disperser, npts, width, nsigma, value, limits, relative) 
    … …  
    204208def call_kernel(kernel, pars, cutoff=0, mono=False): 
    205209    """ 
    206     Call *kernel* returned from *model.make_kernel* with parameters *pars*. 
     210    Call *kernel* returned from :func:`make_kernel` with parameters *pars*. 
    207211 
    208212    *cutoff* is the limiting value for the product of dispersion weights used 
    … …  
    212216    with an error of about 1%, which is usually less than the measurement 
    213217    uncertainty. 
    214  
    215     *mono* is True if polydispersity should be set to none on all parameters. 
    216     """ 
     218    """ 
     219    fixed_pars = [pars.get(name, kernel.info['defaults'][name]) 
     220                  for name in kernel.fixed_pars] 
    217221    if mono: 
    218         active = lambda name: False 
    219     elif kernel.dim == '1d': 
    220         pars_1d = set(p.name for p in kernel.info['parameters'].type['1d']) 
    221         active = lambda name: name in pars_1d 
    222     elif kernel.dim == '2d': 
    223         pars_2d = set(p.name for p in kernel.info['parameters'].type['2d']) 
    224         active = lambda name: name in pars_1d 
    225     else: 
    226         active = lambda name: True 
    227  
    228     vw_pairs = [(get_weights(p, pars) if active(p.name) else ([p.default], [1])) 
    229                 for p in kernel.info['parameters']] 
    230     values, weights = zip(*vw_pairs) 
    231  
    232     if max([len(w) for w in weights]) > 1: 
    233         details = generate.poly_details(kernel.info, weights) 
    234     else: 
    235         details = kernel.info['mono_details'] 
    236  
    237     weights, values = [np.hstack(v) for v in (weights, values)] 
    238  
    239     weights = weights.astype(dtype=kernel.dtype) 
    240     values = values.astype(dtype=kernel.dtype) 
    241     return kernel(details, weights, values, cutoff) 
     222        pd_pars = [( np.array([pars[name]]), np.array([1.0]) ) 
     223                   for name in kernel.pd_pars] 
     224    else: 
     225        pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars] 
     226    return kernel(fixed_pars, pd_pars, cutoff=cutoff) 
    242227 
    243228def call_ER_VR(model_info, vol_pars): 
    … …  
    262247 
    263248 
    264 def call_ER(model_info, values): 
    265     """ 
    266     Call the model ER function using *values*. *model_info* is either 
    267     *model.info* if you have a loaded model, or *kernel.info* if you 
    268     have a model kernel prepared for evaluation. 
    269     """ 
    270     ER = model_info.get('ER', None) 
     249def call_ER(info, pars): 
     250    """ 
     251    Call the model ER function using *pars*. 
     252    *info* is either *model.info* if you have a loaded model, or *kernel.info* 
     253    if you have a model kernel prepared for evaluation. 
     254    """ 
     255    ER = info.get('ER', None) 
    271256    if ER is None: 
    272257        return 1.0 
    273258    else: 
    274         vol_pars = [get_weights(parameter, values) 
    275                     for parameter in model_info['parameters'] 
    276                     if parameter.type == 'volume'] 
     259        vol_pars = [get_weights(info, pars, name) 
     260                    for name in info['partype']['volume']] 
    277261        value, weight = dispersion_mesh(vol_pars) 
    278262        individual_radii = ER(*value) 
    … …  
    280264        return np.sum(weight*individual_radii) / np.sum(weight) 
    281265 
    282 def call_VR(model_info, values): 
     266def call_VR(info, pars): 
    283267    """ 
    284268    Call the model VR function using *pars*. 
    … …  
    286270    if you have a model kernel prepared for evaluation. 
    287271    """ 
    288     VR = model_info.get('VR', None) 
     272    VR = info.get('VR', None) 
    289273    if VR is None: 
    290274        return 1.0 
    291275    else: 
    292         vol_pars = [get_weights(parameter, values) 
    293                     for parameter in model_info['parameters'] 
    294                     if parameter.type == 'volume'] 
     276        vol_pars = [get_weights(info, pars, name) 
     277                    for name in info['partype']['volume']] 
    295278        value, weight = dispersion_mesh(vol_pars) 
    296279        whole, part = VR(*value) 

  • sasmodels/direct_model.py

    r48fbd50 r02e70ff  
    2525import numpy as np 
    2626 
     27from .core import make_kernel 
    2728from .core import call_kernel, call_ER_VR 
    2829from . import sesans 
    … …  
    6566            self.data_type = 'Iq' 
    6667 
     68        partype = model.info['partype'] 
     69 
    6770        if self.data_type == 'sesans': 
    6871             
    … …  
    7881            q_mono = sesans.make_all_q(data) 
    7982        elif self.data_type == 'Iqxy': 
    80             partype = model.info['par_type'] 
    8183            if not partype['orientation'] and not partype['magnetic']: 
    8284                raise ValueError("not 2D without orientation or magnetic parameters") 
    … …  
    152154    def _calc_theory(self, pars, cutoff=0.0): 
    153155        if self._kernel is None: 
    154             self._kernel = self._model.make_kernel(self._kernel_inputs) 
    155             self._kernel_mono = ( 
    156                 self._model.make_kernel(self._kernel_mono_inputs) 
    157                 if self._kernel_mono_inputs else None 
    158             ) 
     156            self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-dedata_type 
     157            self._kernel_mono = make_kernel(self._model, self._kernel_mono_inputs) if self._kernel_mono_inputs else None 
    159158 
    160159        Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) 
    161         Iq_mono = (call_kernel(self._kernel_mono, pars, mono=True) 
    162                    if self._kernel_mono_inputs else None) 
     160        Iq_mono = call_kernel(self._kernel_mono, pars, mono=True) if self._kernel_mono_inputs else None 
    163161        if self.data_type == 'sesans': 
    164162            result = sesans.transform(self._data, 

  • sasmodels/generate.py

    r69aa451 r9ef9dd9  
    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. 
    2723 
    2824These functions are defined in a kernel module .py script and an associated 
    … …  
    6965The constructor code will generate the necessary vectors for computing 
    7066them with the desired polydispersity. 
     67 
     68The available kernel parameters are defined as a list, with each parameter 
     69defined 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 
    71102The kernel module must set variables defining the kernel meta data: 
    72103 
    … …  
    181212from __future__ import print_function 
    182213 
    183 #TODO: determine which functions are useful outside of generate 
    184 #__all__ = ["model_info", "make_doc", "make_source", "convert_type"] 
     214# TODO: identify model files which have changed since loading and reload them. 
    185215 
    186216import sys 
    187217from os.path import abspath, dirname, join as joinpath, exists, basename, \ 
    188     splitext, getmtime 
     218    splitext 
    189219import re 
    190220import string 
    191221import warnings 
     222from collections import namedtuple 
    192223 
    193224import numpy as np 
    194225 
    195 from .modelinfo import ModelInfo, Parameter, make_parameter_table 
    196  
    197 # TODO: identify model files which have changed since loading and reload them. 
    198  
    199 TEMPLATE_ROOT = dirname(__file__) 
    200  
    201 MAX_PD = 4 
     226PARAMETER_FIELDS = ['name', 'units', 'default', 'limits', 'type', 'description'] 
     227Parameter = namedtuple('Parameter', PARAMETER_FIELDS) 
     228 
     229#TODO: determine which functions are useful outside of generate 
     230#__all__ = ["model_info", "make_doc", "make_source", "convert_type"] 
     231 
     232C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c') 
    202233 
    203234F16 = np.dtype('float16') 
    … …  
    208239except TypeError: 
    209240    F128 = None 
     241 
     242# Scale and background, which are parameters common to every form factor 
     243COMMON_PARAMETERS = [ 
     244    ["scale", "", 1, [0, np.inf], "", "Source intensity"], 
     245    ["background", "1/cm", 1e-3, [0, np.inf], "", "Source background"], 
     246    ] 
    210247 
    211248# Conversion from units defined in the parameter table for each model 
    … …  
    301338    raise ValueError("%r not found in %s" % (filename, search_path)) 
    302339 
    303  
    304340def model_sources(model_info): 
    305341    """ 
    … …  
    310346    return [_search(search_path, f) for f in model_info['source']] 
    311347 
    312 def timestamp(model_info): 
    313     """ 
    314     Return a timestamp for the model corresponding to the most recently 
    315     changed file or dependency. 
    316     """ 
    317     source_files = (model_sources(model_info) 
    318                     + model_templates() 
    319                     + [model_info['filename']]) 
    320     newest = max(getmtime(f) for f in source_files) 
    321     return newest 
     348# Pragmas for enable OpenCL features.  Be sure to protect them so that they 
     349# still compile even if OpenCL is not present. 
     350_F16_PRAGMA = """\ 
     351#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 
     352#  pragma OPENCL EXTENSION cl_khr_fp16: enable 
     353#endif 
     354""" 
     355 
     356_F64_PRAGMA = """\ 
     357#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 
     358#  pragma OPENCL EXTENSION cl_khr_fp64: enable 
     359#endif 
     360""" 
    322361 
    323362def convert_type(source, dtype): 
    … …  
    330369    if dtype == F16: 
    331370        fbytes = 2 
    332         source = _convert_type(source, "float", "f") 
     371        source = _F16_PRAGMA + _convert_type(source, "half", "f") 
    333372    elif dtype == F32: 
    334373        fbytes = 4 
    … …  
    336375    elif dtype == F64: 
    337376        fbytes = 8 
    338         # no need to convert if it is already double 
     377        source = _F64_PRAGMA + source  # Source is already double 
    339378    elif dtype == F128: 
    340379        fbytes = 16 
    … …  
    379418 
    380419 
    381 _template_cache = {} 
    382 def load_template(filename): 
    383     path = joinpath(TEMPLATE_ROOT, filename) 
    384     mtime = getmtime(path) 
    385     if filename not in _template_cache or mtime > _template_cache[filename][0]: 
    386         with open(path) as fid: 
    387             _template_cache[filename] = (mtime, fid.read(), path) 
    388     return _template_cache[filename][1] 
    389  
    390 def model_templates(): 
    391     # TODO: fails DRY; templates are listed in two places. 
    392     # should instead have model_info contain a list of paths 
    393     return [joinpath(TEMPLATE_ROOT, filename) 
    394             for filename in ('kernel_header.c', 'kernel_iq.c')] 
    395  
    396  
    397 _FN_TEMPLATE = """\ 
    398 double %(name)s(%(pars)s); 
    399 double %(name)s(%(pars)s) { 
    400     %(body)s 
    401 } 
    402  
    403  
     420LOOP_OPEN = """\ 
     421for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) { 
     422  const double %(name)s = loops[2*(%(name)s_i%(offset)s)]; 
     423  const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\ 
    404424""" 
    405  
    406 def _gen_fn(name, pars, body): 
    407     """ 
    408     Generate a function given pars and body. 
    409  
    410     Returns the following string:: 
    411  
    412          double fn(double a, double b, ...); 
    413          double fn(double a, double b, ...) { 
    414              .... 
    415          } 
    416     """ 
    417     par_decl = ', '.join(p.as_argument() for p in pars) if pars else 'void' 
    418     return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl} 
    419  
    420 def _call_pars(prefix, pars): 
    421     """ 
    422     Return a list of *prefix.parameter* from parameter items. 
    423     """ 
    424     return [p.as_call_reference(prefix) for p in pars] 
    425  
    426 _IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)", 
    427                            flags=re.MULTILINE) 
    428 def _have_Iqxy(sources): 
    429     """ 
    430     Return true if any file defines Iqxy. 
    431  
    432     Note this is not a C parser, and so can be easily confused by 
    433     non-standard syntax.  Also, it will incorrectly identify the following 
    434     as having Iqxy:: 
    435  
    436         /* 
    437         double Iqxy(qx, qy, ...) { ... fill this in later ... } 
    438         */ 
    439  
    440     If you want to comment out an Iqxy function, use // on the front of the 
    441     line instead. 
    442     """ 
    443     for code in sources: 
    444         if _IQXY_PATTERN.search(code): 
    445             return True 
    446     else: 
    447         return False 
    448  
     425def build_polydispersity_loops(pd_pars): 
     426    """ 
     427    Build polydispersity loops 
     428 
     429    Returns loop opening and loop closing 
     430    """ 
     431    depth = 4 
     432    offset = "" 
     433    loop_head = [] 
     434    loop_end = [] 
     435    for name in pd_pars: 
     436        subst = {'name': name, 'offset': offset} 
     437        loop_head.append(indent(LOOP_OPEN % subst, depth)) 
     438        loop_end.insert(0, (" "*depth) + "}") 
     439        offset += '+N' + name 
     440        depth += 2 
     441    return "\n".join(loop_head), "\n".join(loop_end) 
     442 
     443C_KERNEL_TEMPLATE = None 
    449444def make_source(model_info): 
    450445    """ 
    … …  
    465460    # for computing volume even if we allow non-disperse volume parameters. 
    466461 
    467     partable = model_info['parameters'] 
    468  
    469     # Identify parameters for Iq, Iqxy, Iq_magnetic and form_volume. 
    470     # Note that scale and volume are not possible types. 
    471  
    472     # Load templates and user code 
    473     kernel_header = load_template('kernel_header.c') 
    474     kernel_code = load_template('kernel_iq.c') 
    475     user_code = [open(f).read() for f in model_sources(model_info)] 
    476  
    477     # Build initial sources 
    478     source = [kernel_header] + user_code 
    479  
    480     vol_parameters = partable.kernel_pars('volume') 
    481     iq_parameters = partable.kernel_pars('1d') 
    482     iqxy_parameters = partable.kernel_pars('2d') 
    483  
    484     # Make parameters for q, qx, qy so that we can use them in declarations 
    485     q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 
    486     # Generate form_volume function, etc. from body only 
     462    # Load template 
     463    global C_KERNEL_TEMPLATE 
     464    if C_KERNEL_TEMPLATE is None: 
     465        with open(C_KERNEL_TEMPLATE_PATH) as fid: 
     466            C_KERNEL_TEMPLATE = fid.read() 
     467 
     468    # Load additional sources 
     469    source = [open(f).read() for f in model_sources(model_info)] 
     470 
     471    # Prepare defines 
     472    defines = [] 
     473    partype = model_info['partype'] 
     474    pd_1d = partype['pd-1d'] 
     475    pd_2d = partype['pd-2d'] 
     476    fixed_1d = partype['fixed-1d'] 
     477    fixed_2d = partype['fixed-1d'] 
     478 
     479    iq_parameters = [p.name 
     480                     for p in model_info['parameters'][2:]  # skip scale, background 
     481                     if p.name in set(fixed_1d + pd_1d)] 
     482    iqxy_parameters = [p.name 
     483                       for p in model_info['parameters'][2:]  # skip scale, background 
     484                       if p.name in set(fixed_2d + pd_2d)] 
     485    volume_parameters = [p.name 
     486                         for p in model_info['parameters'] 
     487                         if p.type == 'volume'] 
     488 
     489    # Fill in defintions for volume parameters 
     490    if volume_parameters: 
     491        defines.append(('VOLUME_PARAMETERS', 
     492                        ','.join(volume_parameters))) 
     493        defines.append(('VOLUME_WEIGHT_PRODUCT', 
     494                        '*'.join(p + '_w' for p in volume_parameters))) 
     495 
     496    # Generate form_volume function from body only 
    487497    if model_info['form_volume'] is not None: 
    488         pars = vol_parameters 
    489         source.append(_gen_fn('form_volume', pars, model_info['form_volume'])) 
     498        if volume_parameters: 
     499            vol_par_decl = ', '.join('double ' + p for p in volume_parameters) 
     500        else: 
     501            vol_par_decl = 'void' 
     502        defines.append(('VOLUME_PARAMETER_DECLARATIONS', 
     503                        vol_par_decl)) 
     504        fn = """\ 
     505double form_volume(VOLUME_PARAMETER_DECLARATIONS); 
     506double form_volume(VOLUME_PARAMETER_DECLARATIONS) { 
     507    %(body)s 
     508} 
     509""" % {'body':model_info['form_volume']} 
     510        source.append(fn) 
     511 
     512    # Fill in definitions for Iq parameters 
     513    defines.append(('IQ_KERNEL_NAME', model_info['name'] + '_Iq')) 
     514    defines.append(('IQ_PARAMETERS', ', '.join(iq_parameters))) 
     515    if fixed_1d: 
     516        defines.append(('IQ_FIXED_PARAMETER_DECLARATIONS', 
     517                        ', \\\n    '.join('const double %s' % p for p in fixed_1d))) 
     518    if pd_1d: 
     519        defines.append(('IQ_WEIGHT_PRODUCT', 
     520                        '*'.join(p + '_w' for p in pd_1d))) 
     521        defines.append(('IQ_DISPERSION_LENGTH_DECLARATIONS', 
     522                        ', \\\n    '.join('const int N%s' % p for p in pd_1d))) 
     523        defines.append(('IQ_DISPERSION_LENGTH_SUM', 
     524                        '+'.join('N' + p for p in pd_1d))) 
     525        open_loops, close_loops = build_polydispersity_loops(pd_1d) 
     526        defines.append(('IQ_OPEN_LOOPS', 
     527                        open_loops.replace('\n', ' \\\n'))) 
     528        defines.append(('IQ_CLOSE_LOOPS', 
     529                        close_loops.replace('\n', ' \\\n'))) 
    490530    if model_info['Iq'] is not None: 
    491         pars = [q] + iq_parameters 
    492         source.append(_gen_fn('Iq', pars, model_info['Iq'])) 
     531        defines.append(('IQ_PARAMETER_DECLARATIONS', 
     532                        ', '.join('double ' + p for p in iq_parameters))) 
     533        fn = """\ 
     534double Iq(double q, IQ_PARAMETER_DECLARATIONS); 
     535double Iq(double q, IQ_PARAMETER_DECLARATIONS) { 
     536    %(body)s 
     537} 
     538""" % {'body':model_info['Iq']} 
     539        source.append(fn) 
     540 
     541    # Fill in definitions for Iqxy parameters 
     542    defines.append(('IQXY_KERNEL_NAME', model_info['name'] + '_Iqxy')) 
     543    defines.append(('IQXY_PARAMETERS', ', '.join(iqxy_parameters))) 
     544    if fixed_2d: 
     545        defines.append(('IQXY_FIXED_PARAMETER_DECLARATIONS', 
     546                        ', \\\n    '.join('const double %s' % p for p in fixed_2d))) 
     547    if pd_2d: 
     548        defines.append(('IQXY_WEIGHT_PRODUCT', 
     549                        '*'.join(p + '_w' for p in pd_2d))) 
     550        defines.append(('IQXY_DISPERSION_LENGTH_DECLARATIONS', 
     551                        ', \\\n    '.join('const int N%s' % p for p in pd_2d))) 
     552        defines.append(('IQXY_DISPERSION_LENGTH_SUM', 
     553                        '+'.join('N' + p for p in pd_2d))) 
     554        open_loops, close_loops = build_polydispersity_loops(pd_2d) 
     555        defines.append(('IQXY_OPEN_LOOPS', 
     556                        open_loops.replace('\n', ' \\\n'))) 
     557        defines.append(('IQXY_CLOSE_LOOPS', 
     558                        close_loops.replace('\n', ' \\\n'))) 
    493559    if model_info['Iqxy'] is not None: 
    494         pars = [qx, qy] + iqxy_parameters 
    495         source.append(_gen_fn('Iqxy', pars, model_info['Iqxy'])) 
    496  
    497     # Define the parameter table 
    498     source.append("#define PARAMETER_TABLE \\") 
    499     source.append("\\\n".join(p.as_definition() 
    500                                   for p in model_info['parameters'][2:])) 
    501  
    502     # Define the function calls 
    503     if vol_parameters: 
    504         refs = _call_pars("v.", vol_parameters) 
    505         call_volume = "#define CALL_VOLUME(v) form_volume(%s)" % (",".join(refs)) 
    506     else: 
    507         # Model doesn't have volume.  We could make the kernel run a little 
    508         # faster by not using/transferring the volume normalizations, but 
    509         # the ifdef's reduce readability more than is worthwhile. 
    510         call_volume = "#define CALL_VOLUME(v) 0.0" 
    511     source.append(call_volume) 
    512  
    513     refs = ["q[i]"] + _call_pars("v.", iq_parameters) 
    514     call_iq = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(refs)) 
    515     if _have_Iqxy(user_code): 
    516         # Call 2D model 
    517         refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("v.", iqxy_parameters) 
    518         call_iqxy = "#define CALL_IQ(q,i,v) Iqxy(%s)" % (",".join(refs)) 
    519     else: 
    520         # Call 1D model with sqrt(qx^2 + qy^2) 
    521         warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 
    522         # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 
    523         pars_sqrt = ["sqrt(q[2*i]*q[2*i]+q[2*i+1]*q[2*i+1])"] + refs[1:] 
    524         call_iqxy = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(pars_sqrt)) 
    525  
    526     # Fill in definitions for numbers of parameters 
    527     source.append("#define MAX_PD %s"%model_info['max_pd']) 
    528     source.append("#define NPARS %d"%(len(partable.kernel_pars()))) 
    529  
    530     # TODO: allow mixed python/opencl kernels? 
    531  
    532     # define the Iq kernel 
    533     source.append("#define KERNEL_NAME %s_Iq"%model_info['name']) 
    534     source.append(call_iq) 
    535     source.append(kernel_code) 
    536     source.append("#undef CALL_IQ") 
    537     source.append("#undef KERNEL_NAME") 
    538  
    539     # define the Iqxy kernel from the same source with different #defines 
    540     source.append("#define KERNEL_NAME %s_Iqxy"%model_info['name']) 
    541     source.append(call_iqxy) 
    542     source.append(kernel_code) 
    543     source.append("#undef CALL_IQ") 
    544     source.append("#undef KERNEL_NAME") 
    545  
    546     return '\n'.join(source) 
     560        defines.append(('IQXY_PARAMETER_DECLARATIONS', 
     561                        ', '.join('double ' + p for p in iqxy_parameters))) 
     562        fn = """\ 
     563double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS); 
     564double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS) { 
     565    %(body)s 
     566} 
     567""" % {'body':model_info['Iqxy']} 
     568        source.append(fn) 
     569 
     570    # Need to know if we have a theta parameter for Iqxy; it is not there 
     571    # for the magnetic sphere model, for example, which has a magnetic 
     572    # orientation but no shape orientation. 
     573    if 'theta' in pd_2d: 
     574        defines.append(('IQXY_HAS_THETA', '1')) 
     575 
     576    #for d in defines: print(d) 
     577    defines = '\n'.join('#define %s %s' % (k, v) for k, v in defines) 
     578    sources = '\n\n'.join(source) 
     579    return C_KERNEL_TEMPLATE % { 
     580        'DEFINES': defines, 
     581        'SOURCES': sources, 
     582        } 
    547583 
    548584def categorize_parameters(pars): 
    549585    """ 
    550     Categorize the parameters by use: 
    551  
    552     * *pd* list of polydisperse parameters in order; gui should test whether 
    553       they are in *2d* or *magnetic* as appropriate for the data 
    554     * *1d* set of parameters that are used to compute 1D patterns 
    555     * *2d* set of parameters that are used to compute 2D patterns (which 
    556       includes all 1D parameters) 
    557     * *magnetic* set of parameters that are used to compute magnetic 
    558       patterns (which includes all 1D and 2D parameters) 
    559     * *pd_relative* is the set of parameters with relative distribution 
    560       width (e.g., radius +/- 10%) rather than absolute distribution 
    561       width (e.g., theta +/- 6 degrees). 
    562     * *theta_par* is the index of the polar angle polydispersion parameter 
    563       or -1 if no such parameter exists 
    564     """ 
    565     par_set = {} 
     586    Build parameter categories out of the the parameter definitions. 
     587 
     588    Returns a dictionary of categories. 
     589 
     590    Note: these categories are subject to change, depending on the needs of 
     591    the UI and the needs of the kernel calling function. 
     592 
     593    The categories are as follows: 
     594 
     595    * *volume* list of volume parameter names 
     596    * *orientation* list of orientation parameters 
     597    * *magnetic* list of magnetic parameters 
     598    * *<empty string>* list of parameters that have no type info 
     599 
     600    Each parameter is in one and only one category. 
     601 
     602    The following derived categories are created: 
     603 
     604    * *fixed-1d* list of non-polydisperse parameters for 1D models 
     605    * *pd-1d* list of polydisperse parameters for 1D models 
     606    * *fixed-2d* list of non-polydisperse parameters for 2D models 
     607    * *pd-d2* list of polydisperse parameters for 2D models 
     608    """ 
     609    partype = { 
     610        'volume': [], 'orientation': [], 'magnetic': [], 'sld': [], '': [], 
     611        'fixed-1d': [], 'fixed-2d': [], 'pd-1d': [], 'pd-2d': [], 
     612        'pd-rel': set(), 
     613    } 
     614 
     615    for p in pars: 
     616        if p.type == 'volume': 
     617            partype['pd-1d'].append(p.name) 
     618            partype['pd-2d'].append(p.name) 
     619            partype['pd-rel'].add(p.name) 
     620        elif p.type == 'magnetic': 
     621            partype['fixed-2d'].append(p.name) 
     622        elif p.type == 'orientation': 
     623            partype['pd-2d'].append(p.name) 
     624        elif p.type in ('', 'sld'): 
     625            partype['fixed-1d'].append(p.name) 
     626            partype['fixed-2d'].append(p.name) 
     627        else: 
     628            raise ValueError("unknown parameter type %r" % p.type) 
     629        partype[p.type].append(p.name) 
     630 
     631    return partype 
    566632 
    567633def process_parameters(model_info): 
    … …  
    569635    Process parameter block, precalculating parameter details. 
    570636    """ 
    571     partable = model_info['parameters'] 
     637    # convert parameters into named tuples 
     638    for p in model_info['parameters']: 
     639        if p[4] == '' and (p[0].startswith('sld') or p[0].endswith('sld')): 
     640            p[4] = 'sld' 
     641            # TODO: make sure all models explicitly label their sld parameters 
     642            #raise ValueError("%s.%s needs to be explicitly set to type 'sld'" %(model_info['id'], p[0])) 
     643 
     644    pars = [Parameter(*p) for p in model_info['parameters']] 
     645    # Fill in the derived attributes 
     646    model_info['parameters'] = pars 
     647    partype = categorize_parameters(pars) 
     648    model_info['limits'] = dict((p.name, p.limits) for p in pars) 
     649    model_info['partype'] = partype 
     650    model_info['defaults'] = dict((p.name, p.default) for p in pars) 
    572651    if model_info.get('demo', None) is None: 
    573         model_info['demo'] = partable.defaults 
    574  
    575     # Don't use more polydisperse parameters than are available in the model 
    576     # Note: we can do polydispersity on arbitrary parameters, so it is not 
    577     # clear that this is a good idea; it does however make the poly_details 
    578     # code easier to write, so we will leave it in for now. 
    579     model_info['max_pd'] = min(partable.num_pd, MAX_PD) 
    580  
    581 def mono_details(model_info): 
    582     # TODO: move max_pd into ParameterTable? 
    583     max_pd = model_info['max_pd'] 
    584     pars = model_info['parameters'].kernel_pars() 
    585     npars = len(pars) 
    586     par_offset = 5*max_pd 
    587     constants_offset = par_offset + 3*npars 
    588  
    589     details = np.zeros(constants_offset + 2, 'int32') 
    590     details[0*max_pd:1*max_pd] = range(max_pd)       # pd_par: arbitrary order; use first 
    591     details[1*max_pd:2*max_pd] = [1]*max_pd          # pd_length: only one element 
    592     details[2*max_pd:3*max_pd] = range(max_pd)       # pd_offset: consecutive 1.0 weights 
    593     details[3*max_pd:4*max_pd] = [1]*max_pd          # pd_stride: vectors of length 1 
    594     details[4*max_pd:5*max_pd] = [0]*max_pd          # pd_isvol: doens't matter if no norm 
    595     details[par_offset+0*npars:par_offset+1*npars] = range(2, npars+2) # par_offset: skip scale and background 
    596     details[par_offset+1*npars:par_offset+2*npars] = [0]*npars         # no coordination 
    597     #details[p+npars] = 1 # par_coord[0] is coordinated with the first par? 
    598     details[par_offset+2*npars:par_offset+3*npars] = 0 # fast coord with 0 
    599     details[constants_offset]   = 1     # fast_coord_count: one fast index 
    600     details[constants_offset+1] = -1    # theta_par: None 
    601     return details 
    602  
    603 def poly_details(model_info, weights): 
    604     weights = weights[2:] 
    605  
    606     # TODO: move max_pd into ParameterTable? 
    607     max_pd = model_info['max_pd'] 
    608     pars = model_info['parameters'].kernel_pars 
    609     npars = len(pars) 
    610     par_offset = 5*max_pd 
    611     constants_offset = par_offset + 3*npars 
    612  
    613     # Decreasing list of polydispersity lengths 
    614     # Note: the reversing view, x[::-1], does not require a copy 
    615     pd_length = np.array([len(w) for w in weights]) 
    616     print (pd_length) 
    617     print (weights) 
    618     pd_offset = np.cumsum(np.hstack((0, pd_length))) 
    619     pd_isvol = np.array([p.type=='volume' for p in pars]) 
    620     idx = np.argsort(pd_length)[::-1][:max_pd] 
    621     print (idx) 
    622     pd_stride = np.cumprod(np.hstack((1, pd_length[idx][:-1]))) 
    623     par_offsets = np.cumsum(np.hstack((2, pd_length)))[:-1] 
    624  
    625     theta_par = -1 
    626     if 'theta_par' in model_info: 
    627         theta_par = model_info['theta_par'] 
    628         if theta_par >= 0 and pd_length[theta_par] <= 1: 
    629             theta_par = -1 
    630  
    631     details = np.empty(constants_offset + 2, 'int32') 
    632     details[0*max_pd:1*max_pd] = idx             # pd_par 
    633     details[1*max_pd:2*max_pd] = pd_length[idx] 
    634     details[2*max_pd:3*max_pd] = pd_offset[idx] 
    635     details[3*max_pd:4*max_pd] = pd_stride 
    636     details[4*max_pd:5*max_pd] = pd_isvol[idx] 
    637     details[par_offset+0*npars:par_offset+1*npars] = par_offsets 
    638     details[par_offset+1*npars:par_offset+2*npars] = 0  # no coordination for most 
    639     details[par_offset+2*npars:par_offset+3*npars] = 0  # no fast coord with 0 
    640     coord_offset = par_offset+1*npars 
    641     for k,parameter_num in enumerate(idx): 
    642         details[coord_offset+parameter_num] = 2**k 
    643     details[constants_offset] = 1   # fast_coord_count: one fast index 
    644     details[constants_offset+1] = theta_par 
    645     print ("details",details) 
    646     return details 
    647  
    648 def constrained_poly_details(model_info, weights, constraints): 
    649     # Need to find the independently varying pars and sort them 
    650     # Need to build a coordination list for the dependent variables 
    651     # Need to generate a constraints function which takes values 
    652     # and weights, returning par blocks 
    653     raise NotImplementedError("Can't handle constraints yet") 
    654  
    655  
    656 def create_default_functions(model_info): 
    657     """ 
    658     Autogenerate missing functions, such as Iqxy from Iq. 
    659  
    660     This only works for Iqxy when Iq is written in python. :func:`make_source` 
    661     performs a similar role for Iq written in C. 
    662     """ 
    663     if model_info['Iq'] is not None and model_info['Iqxy'] is None: 
    664         partable = model_info['parameters'] 
    665         if partable.type['1d'] != partable.type['2d']: 
    666             raise ValueError("Iqxy model is missing") 
    667         Iq = model_info['Iq'] 
    668         def Iqxy(qx, qy, **kw): 
    669             return Iq(np.sqrt(qx**2 + qy**2), **kw) 
    670         model_info['Iqxy'] = Iqxy 
    671  
     652        model_info['demo'] = model_info['defaults'] 
     653    model_info['has_2d'] = partype['orientation'] or partype['magnetic'] 
    672654 
    673655def make_model_info(kernel_module): 
    … …  
    696678      model variants (e.g., the list of cases) or is None if there are no 
    697679      model variants 
    698     * *par_type* categorizes the model parameters. See 
     680    * *defaults* is the *{parameter: value}* table built from the parameter 
     681      description table. 
     682    * *limits* is the *{parameter: [min, max]}* table built from the 
     683      parameter description table. 
     684    * *partypes* categorizes the model parameters. See 
    699685      :func:`categorize_parameters` for details. 
    700686    * *demo* contains the *{parameter: value}* map used in compare (and maybe 
    … …  
    713699      *model_info* blocks for the composition objects.  This allows us to 
    714700      build complete product and mixture models from just the info. 
    715     * *max_pd* is the max polydispersity dimension.  This is constant and 
    716       should not be reset.  You may be able to change it when the program 
    717       starts by setting *sasmodels.generate.MAX_PD*. 
    718701 
    719702    """ 
    720703    # TODO: maybe turn model_info into a class ModelDefinition 
    721     #print("make parameter table", kernel_module.parameters) 
    722     parameters = make_parameter_table(kernel_module.parameters) 
     704    parameters = COMMON_PARAMETERS + kernel_module.parameters 
    723705    filename = abspath(kernel_module.__file__) 
    724706    kernel_id = splitext(basename(filename))[0] 
    … …  
    749731    functions = "ER VR form_volume Iq Iqxy shape sesans".split() 
    750732    model_info.update((k, getattr(kernel_module, k, None)) for k in functions) 
    751     create_default_functions(model_info) 
    752     # Precalculate the monodisperse parameters 
    753     # TODO: make this a lazy evaluator 
    754     # make_model_info is called for every model on sasview startup 
    755     model_info['mono_details'] = mono_details(model_info) 
    756733    return model_info 
    757734 

  • sasmodels/kernelcl.py

    rfec69dd r8e0d974  
    7575 
    7676 
    77 # Pragmas for enable OpenCL features.  Be sure to protect them so that they 
    78 # still compile even if OpenCL is not present. 
    79 _F16_PRAGMA = """\ 
    80 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 
    81 #  pragma OPENCL EXTENSION cl_khr_fp16: enable 
    82 #endif 
    83 """ 
    84  
    85 _F64_PRAGMA = """\ 
    86 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 
    87 #  pragma OPENCL EXTENSION cl_khr_fp64: enable 
    88 #endif 
    89 """ 
    90  
    91  
    9277ENV = None 
    9378def environment(): 
    … …  
    157142        raise RuntimeError("%s not supported for devices"%dtype) 
    158143 
    159     source_list = [generate.convert_type(source, dtype)] 
    160  
    161     if dtype == generate.F16: 
    162         source_list.insert(0, _F16_PRAGMA) 
    163     elif dtype == generate.F64: 
    164         source_list.insert(0, _F64_PRAGMA) 
    165  
     144    source = generate.convert_type(source, dtype) 
    166145    # Note: USE_SINCOS makes the intel cpu slower under opencl 
    167146    if context.devices[0].type == cl.device_type.GPU: 
    168         source_list.insert(0, "#define USE_SINCOS\n") 
     147        source = "#define USE_SINCOS\n" + source 
    169148    options = (get_fast_inaccurate_build_options(context.devices[0]) 
    170149               if fast else []) 
    171     source = "\n".join(source) 
    172150    program = cl.Program(context, source).build(options=options) 
    173151    return program 
    … …  
    242220        key = "%s-%s-%s"%(name, dtype, fast) 
    243221        if key not in self.compiled: 
    244             print("compiling",name) 
     222            #print("compiling",name) 
    245223            dtype = np.dtype(dtype) 
    246             program = compile_model(self.get_context(dtype), 
    247                                     str(source), dtype, fast) 
     224            program = compile_model(self.get_context(dtype), source, dtype, fast) 
    248225            self.compiled[key] = program 
    249226        return self.compiled[key] 
    … …  
    337314        self.program = None 
    338315 
    339     def make_kernel(self, q_vectors): 
     316    def __call__(self, q_vectors): 
    340317        if self.program is None: 
    341318            compiler = environment().compile_program 
    342             self.program = compiler(self.info['name'], self.source, 
    343                                     self.dtype, self.fast) 
     319            self.program = compiler(self.info['name'], self.source, self.dtype, 
     320                                    self.fast) 
    344321        is_2d = len(q_vectors) == 2 
    345322        kernel_name = generate.kernel_name(self.info, is_2d) 
    … …  
    388365        # at this point, so instead using 32, which is good on the set of 
    389366        # architectures tested so far. 
    390         if self.is_2d: 
    391             # Note: 17 rather than 15 because results is 2 elements 
    392             # longer than input. 
    393             width = ((self.nq+17)//16)*16 
    394             self.q = np.empty((width, 2), dtype=dtype) 
    395             self.q[:self.nq, 0] = q_vectors[0] 
    396             self.q[:self.nq, 1] = q_vectors[1] 
    397         else: 
    398             # Note: 33 rather than 31 because results is 2 elements 
    399             # longer than input. 
    400             width = ((self.nq+33)//32)*32 
    401             self.q = np.empty(width, dtype=dtype) 
    402             self.q[:self.nq] = q_vectors[0] 
    403         self.global_size = [self.q.shape[0]] 
     367        self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 
    404368        context = env.get_context(self.dtype) 
     369        self.global_size = [self.q_vectors[0].size] 
    405370        #print("creating inputs of size", self.global_size) 
    406         self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    407                              hostbuf=self.q) 
     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        ] 
    408375 
    409376    def release(self): 
    … …  
    411378        Free the memory. 
    412379        """ 
    413         if self.q is not None: 
    414             self.q.release() 
    415             self.q = None 
     380        for b in self.q_buffers: 
     381            b.release() 
     382        self.q_buffers = [] 
    416383 
    417384    def __del__(self): 
    … …  
    439406    """ 
    440407    def __init__(self, kernel, model_info, q_vectors, dtype): 
    441         max_pd = self.info['max_pd'] 
    442         npars = len(model_info['parameters'])-2 
    443408        q_input = GpuInput(q_vectors, dtype) 
    444         self.dtype = dtype 
    445         self.dim = '2d' if q_input.is_2d else '1d' 
    446409        self.kernel = kernel 
    447410        self.info = model_info 
    448         self.pd_stop_index = 4*max_pd-1 
    449         # plus three for the normalization values 
    450         self.result = np.empty(q_input.nq+3, q_input.dtype) 
     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] 
    451415 
    452416        # Inputs and outputs for each kernel call 
    … …  
    454418        env = environment() 
    455419        self.queue = env.get_queue(dtype) 
    456  
    457         # details is int32 data, padded to an 8 integer boundary 
    458         size = ((max_pd*5 + npars*3 + 2 + 7)//8)*8 
    459         self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
     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, 
    460423                               q_input.global_size[0] * q_input.dtype.itemsize) 
    461         self.q_input = q_input # allocated by GpuInput above 
    462  
    463         self._need_release = [ self.result_b, self.q_input ] 
    464  
    465     def __call__(self, details, weights, values, cutoff): 
     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): 
    466429        real = (np.float32 if self.q_input.dtype == generate.F32 
    467430                else np.float64 if self.q_input.dtype == generate.F64 
    468431                else np.float16 if self.q_input.dtype == generate.F16 
    469432                else np.float32)  # will never get here, so use np.float32 
    470         assert details.dtype == np.int32 
    471         assert weights.dtype == real and values.dtype == real 
    472  
    473         context = self.queue.context 
    474         details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    475                               hostbuf=details) 
    476         weights_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    477                               hostbuf=weights) 
    478         values_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    479                              hostbuf=values) 
    480  
    481         start, stop = 0, self.details[self.pd_stop_index] 
    482         args = [ 
    483             np.uint32(self.q_input.nq), np.uint32(start), np.uint32(stop), 
    484             self.details_b, self.weights_b, self.values_b, 
    485             self.q_input.q_b, self.result_b, real(cutoff), 
    486         ] 
     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 
    487460        self.kernel(self.queue, self.q_input.global_size, None, *args) 
    488         cl.enqueue_copy(self.queue, self.result, self.result_b) 
    489         [v.release() for v in details_b, weights_b, values_b] 
    490  
    491         return self.result[:self.nq] 
     461        cl.enqueue_copy(self.queue, self.res, res_bi) 
     462 
     463        return self.res 
    492464 
    493465    def release(self): 

  • sasmodels/kerneldll.py

    r69aa451 r6ad0e87  
    5050import tempfile 
    5151import ctypes as ct 
    52 from ctypes import c_void_p, c_int32, c_longdouble, c_double, c_float 
     52from ctypes import c_void_p, c_int, c_longdouble, c_double, c_float 
     53import _ctypes 
    5354 
    5455import numpy as np 
    … …  
    8081        COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 
    8182        if "SAS_OPENMP" in os.environ: 
    82             COMPILE += " -fopenmp" 
     83            COMPILE = COMPILE + " -fopenmp" 
    8384else: 
    8485    COMPILE = "cc -shared -fPIC -fopenmp -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 
    … …  
    140141  
    141142    source = generate.convert_type(source, dtype) 
    142     newest = generate.timestamp(model_info) 
     143    source_files = generate.model_sources(model_info) + [model_info['filename']] 
    143144    dll = dll_path(model_info, dtype) 
     145    newest = max(os.path.getmtime(f) for f in source_files) 
    144146    if not os.path.exists(dll) or os.path.getmtime(dll) < newest: 
    145147        # Replace with a proper temp file 
    … …  
    169171    return DllModel(filename, model_info, dtype=dtype) 
    170172 
     173 
     174IQ_ARGS = [c_void_p, c_void_p, c_int] 
     175IQXY_ARGS = [c_void_p, c_void_p, c_void_p, c_int] 
     176 
    171177class DllModel(object): 
    172178    """ 
    … …  
    191197 
    192198    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 
    193204        #print("dll", self.dllpath) 
    194205        try: 
    … …  
    201212              else c_double if self.dtype == generate.F64 
    202213              else c_longdouble) 
    203  
    204         # int, int, int, int*, double*, double*, double*, double*, double*, double 
    205         argtypes = [c_int32]*3 + [c_void_p]*5 + [fp] 
     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 [] 
    206216        self.Iq = self.dll[generate.kernel_name(self.info, False)] 
     217        self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d 
     218 
    207219        self.Iqxy = self.dll[generate.kernel_name(self.info, True)] 
    208         self.Iq.argtypes = argtypes 
    209         self.Iqxy.argtypes = argtypes 
     220        self.Iqxy.argtypes = IQXY_ARGS + pd_args_2d + [fp]*Nfixed2d 
     221         
     222        self.release() 
    210223 
    211224    def __getstate__(self): 
    … …  
    216229        self.dll = None 
    217230 
    218     def make_kernel(self, q_vectors): 
     231    def __call__(self, q_vectors): 
    219232        q_input = PyInput(q_vectors, self.dtype) 
    220233        if self.dll is None: self._load_dll() 
    … …  
    259272    """ 
    260273    def __init__(self, kernel, model_info, q_input): 
    261         self.kernel = kernel 
    262274        self.info = model_info 
    263275        self.q_input = q_input 
    264         self.dtype = q_input.dtype 
    265         self.dim = '2d' if q_input.is_2d else '1d' 
    266         self.result = np.empty(q_input.nq+3, q_input.dtype) 
    267  
    268     def __call__(self, details, weights, values, cutoff): 
     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): 
    269286        real = (np.float32 if self.q_input.dtype == generate.F32 
    270287                else np.float64 if self.q_input.dtype == generate.F64 
    271288                else np.float128) 
    272         assert details.dtype == np.int32 
    273         assert weights.dtype == real and values.dtype == real 
    274  
    275         max_pd = self.info['max_pd'] 
    276         start, stop = 0, details[4*max_pd-1] 
    277         print(details) 
    278         args = [ 
    279             self.q_input.nq, # nq 
    280             start, # pd_start 
    281             stop, # pd_stop pd_stride[MAX_PD] 
    282             details.ctypes.data, # problem 
    283             weights.ctypes.data,  # weights 
    284             values.ctypes.data,  #pars 
    285             self.q_input.q.ctypes.data, #q 
    286             self.result.ctypes.data,   # results 
    287             real(cutoff), # cutoff 
    288             ] 
     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) 
    289303        self.kernel(*args) 
    290         return self.result[:-3] 
     304 
     305        return self.res 
    291306 
    292307    def release(self): 

  • sasmodels/kernelpy.py

    r48fbd50 ra84a0ca  
    5353        self.dtype = dtype 
    5454        self.is_2d = (len(q_vectors) == 2) 
    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] 
     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] 
    6257 
    6358    def release(self): 
    … …  
    6560        Free resources associated with the model inputs. 
    6661        """ 
    67         self.q = None 
     62        self.q_vectors = [] 
    6863 
    6964class PyKernel(object): 

  • sasmodels/mixture.py

    r69aa451 r72a081d  
    1414import numpy as np 
    1515 
    16 from .modelinfo import Parameter, ParameterTable 
     16from .generate import process_parameters, COMMON_PARAMETERS, Parameter 
    1717 
    1818SCALE=0 
    … …  
    3434 
    3535    # Build new parameter list 
    36     pars = [] 
     36    pars = COMMON_PARAMETERS + [] 
    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 
    4139        prefix = chr(ord('A')+k) + '_' 
    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 
     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_, ... 
    4950            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'] = partable 
     60    model_info['parameters'] = pars 
    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 
    6971 
    7072 

  • sasmodels/model_test.py

    r48fbd50 ra84a0ca  
    5151 
    5252from .core import list_models, load_model_info, build_model, HAVE_OPENCL 
    53 from .core import call_kernel, call_ER, call_VR 
     53from .core import make_kernel, 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 = model.make_kernel(q_vectors) 
     189                kernel = make_kernel(model, q_vectors) 
    190190                actual = call_kernel(kernel, pars) 
    191191            else: 
    192192                q_vectors = [np.array(x)] 
    193                 kernel = model.make_kernel(q_vectors) 
     193                kernel = make_kernel(model, q_vectors) 
    194194                actual = call_kernel(kernel, pars) 
    195195 

  • sasmodels/models/cylinder.c

    r03cac08 r26141cb  
    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) 
    75 
    86double form_volume(double radius, double length) 
    … …  
    1715    double length) 
    1816{ 
     17    // TODO: return NaN if radius<0 or length<0? 
    1918    // precompute qr and qh to save time in the loop 
    2019    const double qr = q*radius; 
    … …  
    4847    double phi) 
    4948{ 
     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

    r03cac08 r30b4ddf  
    1 static double Iq(double q, 
     1double form_volume(void); 
     2 
     3double Iq(double q, 
     4          double guinier_scale, 
     5          double lorentzian_scale, 
     6          double gyration_radius, 
     7          double fractal_exp, 
     8          double cor_length); 
     9 
     10double 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 
     17static double _gel_fit_kernel(double q, 
    218          double guinier_scale, 
    319          double lorentzian_scale, 
    … …  
    824    // Lorentzian Term 
    925    ////////////////////////double a(x[i]*x[i]*zeta*zeta); 
    10     double lorentzian_term = square(q*cor_length); 
     26    double lorentzian_term = q*q*cor_length*cor_length; 
    1127    lorentzian_term = 1.0 + ((fractal_exp + 1.0)/3.0)*lorentzian_term; 
    1228    lorentzian_term = pow(lorentzian_term, (fractal_exp/2.0) ); 
    … …  
    1430    // Exponential Term 
    1531    ////////////////////////double d(x[i]*x[i]*rg*rg); 
    16     double exp_term = square(q*gyration_radius); 
     32    double exp_term = q*q*gyration_radius*gyration_radius; 
    1733    exp_term = exp(-1.0*(exp_term/3.0)); 
    1834 
    … …  
    2137    return result; 
    2238} 
     39double form_volume(void){ 
     40    // Unused, so free to return garbage. 
     41    return NAN; 
     42} 
     43 
     44double 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. 
     60double 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/rpa.py

    r69aa451 raad336c  
    8686#   ["name", "units", default, [lower, upper], "type","description"], 
    8787parameters = [ 
    88     ["case_num", "", 1, CASES, "", "Component organization"], 
     88    ["case_num", CASES, 0, [0, 10], "", "Component organization"], 
    8989 
    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"], 
     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"], 
    95113 
    96114    ["Kab", "", -0.0004, [-inf, inf], "", "Interaction parameter"], 

  • sasmodels/product.py

    r8e0d974 r8e0d974  
    9999        # a parameter map. 
    100100        par_map = {} 
    101         p_info = p_kernel.info['par_type'] 
    102         s_info = s_kernel.info['par_type'] 
     101        p_info = p_kernel.info['partype'] 
     102        s_info = s_kernel.info['partype'] 
    103103        vol_pars = set(p_info['volume']) 
    104104        if dim == '2d': 

  • sasmodels/resolution.py

    r48fbd50 r8b935d1  
    477477    """ 
    478478    from sasmodels import core 
    479     kernel = form.make_kernel([q]) 
     479    kernel = core.make_kernel(form, [q]) 
    480480    theory = core.call_kernel(kernel, pars) 
    481481    kernel.release() 
    … …  
    502502    from scipy.integrate import romberg 
    503503 
    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 
     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 
    507507        raise ValueError("bad parameters: [%s] not in [%s]"% 
    508508                         (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 
    … …  
    556556    from scipy.integrate import romberg 
    557557 
    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 
     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 
    561561        raise ValueError("bad parameters: [%s] not in [%s]"% 
    562                          (", ".join(sorted(extra)), 
    563                           ", ".join(sorted(pars.keys())))) 
     562                         (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 
    564563 
    565564    _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq) 
    … …  
    694693    def _eval_sphere(self, pars, resolution): 
    695694        from sasmodels import core 
    696         kernel = self.model.make_kernel([resolution.q_calc]) 
     695        kernel = core.make_kernel(self.model, [resolution.q_calc]) 
    697696        theory = core.call_kernel(kernel, pars) 
    698697        result = resolution.apply(theory) 
    … …  
    10631062    model = core.build_model(model_info) 
    10641063 
    1065     kernel = model.make_kernel([resolution.q_calc]) 
     1064    kernel = core.make_kernel(model, [resolution.q_calc]) 
    10661065    theory = core.call_kernel(kernel, pars) 
    10671066    Iq = resolution.apply(theory) 
Note: See TracChangeset for help on using the changeset viewer.