Changes in / [34edbb8:e9b1663d] in sasmodels


Ignore:
Files:
6 added
16 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

    r9217ef8 r9217ef8  
    8181    from bumps.names import Parameter 
    8282 
    83     pars = {} 
     83    pars = {} # => floating point parameters 
    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['partype']['pd-2d']: 
     87    for name in model_info['par_type']['pd']: 
    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 = {} 
    98     for name in model_info['partype']['pd-2d']: 
     97    pd_types = {}  # => distribution names 
     98    for name in model_info['par_type']['pd']: 
    9999        xname = name + '_pd_type' 
    100100        xvalue = kwargs.pop(xname, 'gaussian') 

  • sasmodels/compare.py

    r1671636 r69aa451  
    309309        exclude = lambda n: False 
    310310    else: 
    311         partype = model_info['partype'] 
    312         par1d = set(partype['fixed-1d']+partype['pd-1d']) 
     311        par1d = model_info['par_type']['1d'] 
    313312        exclude = lambda n: n not in par1d 
    314313    lines = [] 
    … …  
    454453    """ 
    455454    # initialize the code so time is more accurate 
    456     value = calculator(**suppress_pd(pars)) 
     455    if Nevals > 1: 
     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             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)) 
     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)) 
    886885        else: 
    887886            for k, v in pars.items(): 

  • sasmodels/core.py

    re79f0a5 r69aa451  
    2626__all__ = [ 
    2727    "list_models", "load_model_info", "precompile_dll", 
    28     "build_model", "make_kernel", "call_kernel", "call_ER_VR", 
     28    "build_model", "call_kernel", "call_ER_VR", 
    2929] 
    3030 
    … …  
    167167 
    168168 
    169 def 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  
    175 def get_weights(model_info, pars, name): 
     169def get_weights(parameter, values): 
    176170    """ 
    177171    Generate the distribution for parameter *name* given the parameter values 
    … …  
    181175    from the *pars* dictionary for parameter value and parameter dispersion. 
    182176    """ 
    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) 
     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) 
    190186    value, weight = weights.get_weights( 
    191187        disperser, npts, width, nsigma, value, limits, relative) 
    … …  
    208204def call_kernel(kernel, pars, cutoff=0, mono=False): 
    209205    """ 
    210     Call *kernel* returned from :func:`make_kernel` with parameters *pars*. 
     206    Call *kernel* returned from *model.make_kernel* with parameters *pars*. 
    211207 
    212208    *cutoff* is the limiting value for the product of dispersion weights used 
    … …  
    216212    with an error of about 1%, which is usually less than the measurement 
    217213    uncertainty. 
    218     """ 
    219     fixed_pars = [pars.get(name, kernel.info['defaults'][name]) 
    220                   for name in kernel.fixed_pars] 
     214 
     215    *mono* is True if polydispersity should be set to none on all parameters. 
     216    """ 
    221217    if mono: 
    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) 
     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) 
    227242 
    228243def call_ER_VR(model_info, vol_pars): 
    … …  
    247262 
    248263 
    249 def 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) 
     264def 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) 
    256271    if ER is None: 
    257272        return 1.0 
    258273    else: 
    259         vol_pars = [get_weights(info, pars, name) 
    260                     for name in info['partype']['volume']] 
     274        vol_pars = [get_weights(parameter, values) 
     275                    for parameter in model_info['parameters'] 
     276                    if parameter.type == 'volume'] 
    261277        value, weight = dispersion_mesh(vol_pars) 
    262278        individual_radii = ER(*value) 
    … …  
    264280        return np.sum(weight*individual_radii) / np.sum(weight) 
    265281 
    266 def call_VR(info, pars): 
     282def call_VR(model_info, values): 
    267283    """ 
    268284    Call the model VR function using *pars*. 
    … …  
    270286    if you have a model kernel prepared for evaluation. 
    271287    """ 
    272     VR = info.get('VR', None) 
     288    VR = model_info.get('VR', None) 
    273289    if VR is None: 
    274290        return 1.0 
    275291    else: 
    276         vol_pars = [get_weights(info, pars, name) 
    277                     for name in info['partype']['volume']] 
     292        vol_pars = [get_weights(parameter, values) 
     293                    for parameter in model_info['parameters'] 
     294                    if parameter.type == 'volume'] 
    278295        value, weight = dispersion_mesh(vol_pars) 
    279296        whole, part = VR(*value) 

  • sasmodels/direct_model.py

    r02e70ff r48fbd50  
    2525import numpy as np 
    2626 
    27 from .core import make_kernel 
    2827from .core import call_kernel, call_ER_VR 
    2928from . import sesans 
    … …  
    6665            self.data_type = 'Iq' 
    6766 
    68         partype = model.info['partype'] 
    69  
    7067        if self.data_type == 'sesans': 
    7168             
    … …  
    8178            q_mono = sesans.make_all_q(data) 
    8279        elif self.data_type == 'Iqxy': 
     80            partype = model.info['par_type'] 
    8381            if not partype['orientation'] and not partype['magnetic']: 
    8482                raise ValueError("not 2D without orientation or magnetic parameters") 
    … …  
    154152    def _calc_theory(self, pars, cutoff=0.0): 
    155153        if self._kernel is None: 
    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 
     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            ) 
    158159 
    159160        Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) 
    160         Iq_mono = call_kernel(self._kernel_mono, pars, mono=True) if self._kernel_mono_inputs else None 
     161        Iq_mono = (call_kernel(self._kernel_mono, pars, mono=True) 
     162                   if self._kernel_mono_inputs else None) 
    161163        if self.data_type == 'sesans': 
    162164            result = sesans.transform(self._data, 

  • sasmodels/generate.py

    r9ef9dd9 r69aa451  
    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 
    … …  
    212181from __future__ import print_function 
    213182 
    214 # TODO: identify model files which have changed since loading and reload them. 
     183#TODO: determine which functions are useful outside of generate 
     184#__all__ = ["model_info", "make_doc", "make_source", "convert_type"] 
    215185 
    216186import sys 
    217187from os.path import abspath, dirname, join as joinpath, exists, basename, \ 
    218     splitext 
     188    splitext, getmtime 
    219189import re 
    220190import string 
    221191import warnings 
    222 from collections import namedtuple 
    223192 
    224193import numpy as np 
    225194 
    226 PARAMETER_FIELDS = ['name', 'units', 'default', 'limits', 'type', 'description'] 
    227 Parameter = 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  
    232 C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c') 
     195from .modelinfo import ModelInfo, Parameter, make_parameter_table 
     196 
     197# TODO: identify model files which have changed since loading and reload them. 
     198 
     199TEMPLATE_ROOT = dirname(__file__) 
     200 
     201MAX_PD = 4 
    233202 
    234203F16 = np.dtype('float16') 
    … …  
    239208except TypeError: 
    240209    F128 = None 
    241  
    242 # Scale and background, which are parameters common to every form factor 
    243 COMMON_PARAMETERS = [ 
    244     ["scale", "", 1, [0, np.inf], "", "Source intensity"], 
    245     ["background", "1/cm", 1e-3, [0, np.inf], "", "Source background"], 
    246     ] 
    247210 
    248211# Conversion from units defined in the parameter table for each model 
    … …  
    338301    raise ValueError("%r not found in %s" % (filename, search_path)) 
    339302 
     303 
    340304def model_sources(model_info): 
    341305    """ 
    … …  
    346310    return [_search(search_path, f) for f in model_info['source']] 
    347311 
    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 """ 
     312def 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 
    361322 
    362323def convert_type(source, dtype): 
    … …  
    369330    if dtype == F16: 
    370331        fbytes = 2 
    371         source = _F16_PRAGMA + _convert_type(source, "half", "f") 
     332        source = _convert_type(source, "float", "f") 
    372333    elif dtype == F32: 
    373334        fbytes = 4 
    … …  
    375336    elif dtype == F64: 
    376337        fbytes = 8 
    377         source = _F64_PRAGMA + source  # Source is already double 
     338        # no need to convert if it is already double 
    378339    elif dtype == F128: 
    379340        fbytes = 16 
    … …  
    418379 
    419380 
    420 LOOP_OPEN = """\ 
    421 for (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];\ 
     381_template_cache = {} 
     382def 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 
     390def 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 = """\ 
     398double %(name)s(%(pars)s); 
     399double %(name)s(%(pars)s) { 
     400    %(body)s 
     401} 
     402 
     403 
    424404""" 
    425 def 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  
    443 C_KERNEL_TEMPLATE = None 
     405 
     406def _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 
     420def _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) 
     428def _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 
    444449def make_source(model_info): 
    445450    """ 
    … …  
    460465    # for computing volume even if we allow non-disperse volume parameters. 
    461466 
    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 
     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 
    497487    if model_info['form_volume'] is not None: 
    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 = """\ 
    505 double form_volume(VOLUME_PARAMETER_DECLARATIONS); 
    506 double 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'))) 
     488        pars = vol_parameters 
     489        source.append(_gen_fn('form_volume', pars, model_info['form_volume'])) 
    530490    if model_info['Iq'] is not None: 
    531         defines.append(('IQ_PARAMETER_DECLARATIONS', 
    532                         ', '.join('double ' + p for p in iq_parameters))) 
    533         fn = """\ 
    534 double Iq(double q, IQ_PARAMETER_DECLARATIONS); 
    535 double 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'))) 
     491        pars = [q] + iq_parameters 
     492        source.append(_gen_fn('Iq', pars, model_info['Iq'])) 
    559493    if model_info['Iqxy'] is not None: 
    560         defines.append(('IQXY_PARAMETER_DECLARATIONS', 
    561                         ', '.join('double ' + p for p in iqxy_parameters))) 
    562         fn = """\ 
    563 double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS); 
    564 double 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         } 
     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) 
    583547 
    584548def categorize_parameters(pars): 
    585549    """ 
    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 
     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 = {} 
    632566 
    633567def process_parameters(model_info): 
    … …  
    635569    Process parameter block, precalculating parameter details. 
    636570    """ 
    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) 
     571    partable = model_info['parameters'] 
    651572    if model_info.get('demo', None) is None: 
    652         model_info['demo'] = model_info['defaults'] 
    653     model_info['has_2d'] = partype['orientation'] or partype['magnetic'] 
     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 
     581def 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 
     603def 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 
     648def 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 
     656def 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 
    654672 
    655673def make_model_info(kernel_module): 
    … …  
    678696      model variants (e.g., the list of cases) or is None if there are no 
    679697      model variants 
    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 
     698    * *par_type* categorizes the model parameters. See 
    685699      :func:`categorize_parameters` for details. 
    686700    * *demo* contains the *{parameter: value}* map used in compare (and maybe 
    … …  
    699713      *model_info* blocks for the composition objects.  This allows us to 
    700714      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*. 
    701718 
    702719    """ 
    703720    # TODO: maybe turn model_info into a class ModelDefinition 
    704     parameters = COMMON_PARAMETERS + kernel_module.parameters 
     721    #print("make parameter table", kernel_module.parameters) 
     722    parameters = make_parameter_table(kernel_module.parameters) 
    705723    filename = abspath(kernel_module.__file__) 
    706724    kernel_id = splitext(basename(filename))[0] 
    … …  
    731749    functions = "ER VR form_volume Iq Iqxy shape sesans".split() 
    732750    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) 
    733756    return model_info 
    734757 

  • sasmodels/kernelcl.py

    r8e0d974 rfec69dd  
    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 
    7792ENV = None 
    7893def environment(): 
    … …  
    142157        raise RuntimeError("%s not supported for devices"%dtype) 
    143158 
    144     source = generate.convert_type(source, dtype) 
     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 
    145166    # Note: USE_SINCOS makes the intel cpu slower under opencl 
    146167    if context.devices[0].type == cl.device_type.GPU: 
    147         source = "#define USE_SINCOS\n" + source 
     168        source_list.insert(0, "#define USE_SINCOS\n") 
    148169    options = (get_fast_inaccurate_build_options(context.devices[0]) 
    149170               if fast else []) 
     171    source = "\n".join(source) 
    150172    program = cl.Program(context, source).build(options=options) 
    151173    return program 
    … …  
    220242        key = "%s-%s-%s"%(name, dtype, fast) 
    221243        if key not in self.compiled: 
    222             #print("compiling",name) 
     244            print("compiling",name) 
    223245            dtype = np.dtype(dtype) 
    224             program = compile_model(self.get_context(dtype), source, dtype, fast) 
     246            program = compile_model(self.get_context(dtype), 
     247                                    str(source), dtype, fast) 
    225248            self.compiled[key] = program 
    226249        return self.compiled[key] 
    … …  
    314337        self.program = None 
    315338 
    316     def __call__(self, q_vectors): 
     339    def make_kernel(self, q_vectors): 
    317340        if self.program is None: 
    318341            compiler = environment().compile_program 
    319             self.program = compiler(self.info['name'], self.source, self.dtype, 
    320                                     self.fast) 
     342            self.program = compiler(self.info['name'], self.source, 
     343                                    self.dtype, self.fast) 
    321344        is_2d = len(q_vectors) == 2 
    322345        kernel_name = generate.kernel_name(self.info, is_2d) 
    … …  
    365388        # at this point, so instead using 32, which is good on the set of 
    366389        # architectures tested so far. 
    367         self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 
     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]] 
    368404        context = env.get_context(self.dtype) 
    369         self.global_size = [self.q_vectors[0].size] 
    370405        #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         ] 
     406        self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     407                             hostbuf=self.q) 
    375408 
    376409    def release(self): 
    … …  
    378411        Free the memory. 
    379412        """ 
    380         for b in self.q_buffers: 
    381             b.release() 
    382         self.q_buffers = [] 
     413        if self.q is not None: 
     414            self.q.release() 
     415            self.q = None 
    383416 
    384417    def __del__(self): 
    … …  
    406439    """ 
    407440    def __init__(self, kernel, model_info, q_vectors, dtype): 
     441        max_pd = self.info['max_pd'] 
     442        npars = len(model_info['parameters'])-2 
    408443        q_input = GpuInput(q_vectors, dtype) 
     444        self.dtype = dtype 
     445        self.dim = '2d' if q_input.is_2d else '1d' 
    409446        self.kernel = kernel 
    410447        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] 
     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) 
    415451 
    416452        # Inputs and outputs for each kernel call 
    … …  
    418454        env = environment() 
    419455        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, 
     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, 
    423460                               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): 
     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): 
    429466        real = (np.float32 if self.q_input.dtype == generate.F32 
    430467                else np.float64 if self.q_input.dtype == generate.F64 
    431468                else np.float16 if self.q_input.dtype == generate.F16 
    432469                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 
     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        ] 
    460487        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 
     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] 
    464492 
    465493    def release(self): 

  • sasmodels/kerneldll.py

    r6ad0e87 r69aa451  
    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" 
    … …  
    141140  
    142141    source = generate.convert_type(source, dtype) 
    143     source_files = generate.model_sources(model_info) + [model_info['filename']] 
     142    newest = generate.timestamp(model_info) 
    144143    dll = dll_path(model_info, dtype) 
    145     newest = max(os.path.getmtime(f) for f in source_files) 
    146144    if not os.path.exists(dll) or os.path.getmtime(dll) < newest: 
    147145        # Replace with a proper temp file 
    … …  
    171169    return DllModel(filename, model_info, dtype=dtype) 
    172170 
    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  
    177171class DllModel(object): 
    178172    """ 
    … …  
    197191 
    198192    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  
    204193        #print("dll", self.dllpath) 
    205194        try: 
    … …  
    212201              else c_double if self.dtype == generate.F64 
    213202              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 [] 
     203 
     204        # int, int, int, int*, double*, double*, double*, double*, double*, double 
     205        argtypes = [c_int32]*3 + [c_void_p]*5 + [fp] 
    216206        self.Iq = self.dll[generate.kernel_name(self.info, False)] 
    217         self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d 
    218  
    219207        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() 
     208        self.Iq.argtypes = argtypes 
     209        self.Iqxy.argtypes = argtypes 
    223210 
    224211    def __getstate__(self): 
    … …  
    229216        self.dll = None 
    230217 
    231     def __call__(self, q_vectors): 
     218    def make_kernel(self, q_vectors): 
    232219        q_input = PyInput(q_vectors, self.dtype) 
    233220        if self.dll is None: self._load_dll() 
    … …  
    272259    """ 
    273260    def __init__(self, kernel, model_info, q_input): 
     261        self.kernel = kernel 
    274262        self.info = model_info 
    275263        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): 
     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): 
    286269        real = (np.float32 if self.q_input.dtype == generate.F32 
    287270                else np.float64 if self.q_input.dtype == generate.F64 
    288271                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) 
     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            ] 
    303289        self.kernel(*args) 
    304  
    305         return self.res 
     290        return self.result[:-3] 
    306291 
    307292    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/rpa.py

    raad336c r69aa451  
    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

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

  • sasmodels/resolution.py

    r8b935d1 r48fbd50  
    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) 
Note: See TracChangeset for help on using the changeset viewer.