Changeset 5c028e3 in sasmodels


Ignore:
Timestamp:
Mar 27, 2016 6:58:51 PM (9 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
ce896fd
Parents:
d19962c (diff), b274bec (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into polydisp

Files:
7 added
25 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/core_multi_shell.py

    r34edbb8 rb274bec  
    2626 
    2727For information about polarised and magnetic scattering, see  
    28 the :ref:`mag_help` documentation. 
     28the :doc:`magnetic help <../sasgui/perspectives/fitting/mag_help>` documentation. 
    2929 
    3030Our model uses the form factor calculations implemented in a c-library provided 

  • 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 rd19962c  
    8181    from bumps.names import Parameter 
    8282 
    83     pars = {} 
     83    pars = {}     # floating point parameters 
     84    pd_types = {} # distribution names 
    8485    for p in model_info['parameters']: 
    8586        value = kwargs.pop(p.name, p.default) 
    8687        pars[p.name] = Parameter.default(value, name=p.name, limits=p.limits) 
    87     for name in model_info['partype']['pd-2d']: 
    88         for xpart, xdefault, xlimits in [ 
    89                 ('_pd', 0., pars[name].limits), 
    90                 ('_pd_n', 35., (0, 1000)), 
    91                 ('_pd_nsigma', 3., (0, 10)), 
    92             ]: 
    93             xname = name + xpart 
    94             xvalue = kwargs.pop(xname, xdefault) 
    95             pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits) 
    96  
    97     pd_types = {} 
    98     for name in model_info['partype']['pd-2d']: 
    99         xname = name + '_pd_type' 
    100         xvalue = kwargs.pop(xname, 'gaussian') 
    101         pd_types[xname] = xvalue 
     88        if p.polydisperse: 
     89            for part, default, limits in [ 
     90                    ('_pd', 0., pars[p.name].limits), 
     91                    ('_pd_n', 35., (0, 1000)), 
     92                    ('_pd_nsigma', 3., (0, 10)), 
     93                ]: 
     94                name = p.name + part 
     95                value = kwargs.pop(name, default) 
     96                pars[name] = Parameter.default(value, name=name, limits=limits) 
     97            pd_types[p.name+'_pd_type'] = kwargs.pop(name, 'gaussian') 
    10298 
    10399    if kwargs:  # args not corresponding to parameters 

  • sasmodels/compare.py

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

  • sasmodels/core.py

    re79f0a5 rd19962c  
    2424    HAVE_OPENCL = False 
    2525 
     26try: 
     27    # Python 3.5 and up 
     28    from importlib.util import spec_from_file_location, module_from_spec 
     29    def load_module(fullname, path): 
     30        spec = spec_from_file_location(fullname, path) 
     31        module = module_from_spec(spec) 
     32        spec.loader.exec_module(module) 
     33        return module 
     34except ImportError: 
     35    # CRUFT: python 2 
     36    import imp 
     37    def load_module(fullname, path): 
     38        module = imp.load_source(fullname, path) 
     39        return module 
     40 
     41 
     42 
    2643__all__ = [ 
    2744    "list_models", "load_model_info", "precompile_dll", 
    28     "build_model", "make_kernel", "call_kernel", "call_ER_VR", 
     45    "build_model", "call_kernel", "call_ER_VR", 
    2946] 
    3047 
    … …  
    5168    """ 
    5269    return build_model(load_model_info(model_name), **kw) 
    53  
    54 def load_model_info_from_path(path): 
    55     # Pull off the last .ext if it exists; there may be others 
    56     name = basename(splitext(path)[0]) 
    57  
    58     # Not cleaning name since don't need to be able to reload this 
    59     # model later 
    60     # Should probably turf the model from sys.modules after we are done... 
    61  
    62     # Placing the model in the 'sasmodels.custom' name space, even 
    63     # though it doesn't actually exist.  imp.load_source doesn't seem 
    64     # to care. 
    65     kernel_module = imp.load_source('sasmodels.custom.'+name, path) 
    66  
    67     # Now that we have the module, we can load it as usual 
    68     model_info = generate.make_model_info(kernel_module) 
    69     return model_info 
    7070 
    7171def load_model_info(model_name): 
    … …  
    9090        return product.make_product_info(P_info, Q_info) 
    9191 
     92    return make_model_by_name(model_name) 
     93 
     94 
     95def make_model_by_name(model_name): 
     96    if model_name.endswith('.py'): 
     97        path = model_name 
     98        # Pull off the last .ext if it exists; there may be others 
     99        name = basename(splitext(path)[0]) 
     100        # Placing the model in the 'sasmodels.custom' name space. 
     101        from sasmodels import custom 
     102        kernel_module = load_module('sasmodels.custom.'+name, path) 
     103    else: 
     104        from sasmodels import models 
     105        __import__('sasmodels.models.'+model_name) 
     106        kernel_module = getattr(models, model_name, None) 
    92107    #import sys; print "\n".join(sys.path) 
    93     __import__('sasmodels.models.'+model_name) 
    94     kernel_module = getattr(models, model_name, None) 
    95108    return generate.make_model_info(kernel_module) 
    96109 
    … …  
    167180 
    168181 
    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): 
     182def get_weights(parameter, values): 
    176183    """ 
    177184    Generate the distribution for parameter *name* given the parameter values 
    … …  
    181188    from the *pars* dictionary for parameter value and parameter dispersion. 
    182189    """ 
    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) 
     190    value = values.get(parameter.name, parameter.default) 
     191    relative = parameter.relative_pd 
     192    limits = parameter.limits 
     193    disperser = values.get(parameter.name+'_pd_type', 'gaussian') 
     194    npts = values.get(parameter.name+'_pd_n', 0) 
     195    width = values.get(parameter.name+'_pd', 0.0) 
     196    nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 
     197    if npts == 0 or width == 0: 
     198        return [value], [] 
    190199    value, weight = weights.get_weights( 
    191200        disperser, npts, width, nsigma, value, limits, relative) 
    … …  
    208217def call_kernel(kernel, pars, cutoff=0, mono=False): 
    209218    """ 
    210     Call *kernel* returned from :func:`make_kernel` with parameters *pars*. 
     219    Call *kernel* returned from *model.make_kernel* with parameters *pars*. 
    211220 
    212221    *cutoff* is the limiting value for the product of dispersion weights used 
    … …  
    216225    with an error of about 1%, which is usually less than the measurement 
    217226    uncertainty. 
    218     """ 
    219     fixed_pars = [pars.get(name, kernel.info['defaults'][name]) 
    220                   for name in kernel.fixed_pars] 
     227 
     228    *mono* is True if polydispersity should be set to none on all parameters. 
     229    """ 
     230    parameters = kernel.info['parameters'] 
    221231    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) 
     232        active = lambda name: False 
     233    elif kernel.dim == '1d': 
     234        active = lambda name: name in parameters.pd_1d 
     235    elif kernel.dim == '2d': 
     236        active = lambda name: name in parameters.pd_2d 
     237    else: 
     238        active = lambda name: True 
     239 
     240    vw_pairs = [(get_weights(p, pars) if active(p.name) 
     241                 else ([pars.get(p.name, p.default)], [])) 
     242                for p in parameters.call_parameters] 
     243    values, weights = zip(*vw_pairs) 
     244 
     245    if max([len(w) for w in weights]) > 1: 
     246        details = generate.poly_details(kernel.info, weights) 
     247    else: 
     248        details = kernel.info['mono_details'] 
     249 
     250    weights, values = [np.hstack(v) for v in (weights, values)] 
     251 
     252    weights = weights.astype(dtype=kernel.dtype) 
     253    values = values.astype(dtype=kernel.dtype) 
     254    return kernel(details, weights, values, cutoff) 
    227255 
    228256def call_ER_VR(model_info, vol_pars): 
    … …  
    247275 
    248276 
    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) 
     277def call_ER(model_info, values): 
     278    """ 
     279    Call the model ER function using *values*. *model_info* is either 
     280    *model.info* if you have a loaded model, or *kernel.info* if you 
     281    have a model kernel prepared for evaluation. 
     282    """ 
     283    ER = model_info.get('ER', None) 
    256284    if ER is None: 
    257285        return 1.0 
    258286    else: 
    259         vol_pars = [get_weights(info, pars, name) 
    260                     for name in info['partype']['volume']] 
     287        vol_pars = [get_weights(parameter, values) 
     288                    for parameter in model_info['parameters'] 
     289                    if parameter.type == 'volume'] 
    261290        value, weight = dispersion_mesh(vol_pars) 
    262291        individual_radii = ER(*value) 
    … …  
    264293        return np.sum(weight*individual_radii) / np.sum(weight) 
    265294 
    266 def call_VR(info, pars): 
     295def call_VR(model_info, values): 
    267296    """ 
    268297    Call the model VR function using *pars*. 
    … …  
    270299    if you have a model kernel prepared for evaluation. 
    271300    """ 
    272     VR = info.get('VR', None) 
     301    VR = model_info.get('VR', None) 
    273302    if VR is None: 
    274303        return 1.0 
    275304    else: 
    276         vol_pars = [get_weights(info, pars, name) 
    277                     for name in info['partype']['volume']] 
     305        vol_pars = [get_weights(parameter, values) 
     306                    for parameter in model_info['parameters'] 
     307                    if parameter.type == 'volume'] 
    278308        value, weight = dispersion_mesh(vol_pars) 
    279309        whole, part = VR(*value) 

  • sasmodels/direct_model.py

    r02e70ff r60eab2a  
    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': 
    83             if not partype['orientation'] and not partype['magnetic']: 
    84                 raise ValueError("not 2D without orientation or magnetic parameters") 
     80            #if not model.info['parameters'].has_2d: 
     81            #    raise ValueError("not 2D without orientation or magnetic parameters") 
    8582            q = np.sqrt(data.qx_data**2 + data.qy_data**2) 
    8683            qmin = getattr(data, 'qmin', 1e-16) 
    … …  
    154151    def _calc_theory(self, pars, cutoff=0.0): 
    155152        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 
     153            self._kernel = self._model.make_kernel(self._kernel_inputs) 
     154            self._kernel_mono = ( 
     155                self._model.make_kernel(self._kernel_mono_inputs) 
     156                if self._kernel_mono_inputs else None 
     157            ) 
    158158 
    159159        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 
     160        Iq_mono = (call_kernel(self._kernel_mono, pars, mono=True) 
     161                   if self._kernel_mono_inputs else None) 
    161162        if self.data_type == 'sesans': 
    162163            result = sesans.transform(self._data, 

  • sasmodels/generate.py

    r9ef9dd9 rd19962c  
    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. 
    215  
    216 import sys 
     183#TODO: determine which functions are useful outside of generate 
     184#__all__ = ["model_info", "make_doc", "make_source", "convert_type"] 
     185 
    217186from os.path import abspath, dirname, join as joinpath, exists, basename, \ 
    218     splitext 
     187    splitext, getmtime 
    219188import re 
    220189import string 
    221190import warnings 
    222 from collections import namedtuple 
    223191 
    224192import numpy as np 
    225193 
    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') 
     194from .modelinfo import ModelInfo, Parameter, make_parameter_table 
     195 
     196# TODO: identify model files which have changed since loading and reload them. 
     197 
     198TEMPLATE_ROOT = dirname(__file__) 
    233199 
    234200F16 = np.dtype('float16') 
    … …  
    239205except TypeError: 
    240206    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     ] 
    247207 
    248208# Conversion from units defined in the parameter table for each model 
    … …  
    338298    raise ValueError("%r not found in %s" % (filename, search_path)) 
    339299 
     300 
    340301def model_sources(model_info): 
    341302    """ 
    … …  
    346307    return [_search(search_path, f) for f in model_info['source']] 
    347308 
    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 """ 
     309def timestamp(model_info): 
     310    """ 
     311    Return a timestamp for the model corresponding to the most recently 
     312    changed file or dependency. 
     313    """ 
     314    source_files = (model_sources(model_info) 
     315                    + model_templates() 
     316                    + [model_info['filename']]) 
     317    newest = max(getmtime(f) for f in source_files) 
     318    return newest 
    361319 
    362320def convert_type(source, dtype): 
    … …  
    369327    if dtype == F16: 
    370328        fbytes = 2 
    371         source = _F16_PRAGMA + _convert_type(source, "half", "f") 
     329        source = _convert_type(source, "float", "f") 
    372330    elif dtype == F32: 
    373331        fbytes = 4 
    … …  
    375333    elif dtype == F64: 
    376334        fbytes = 8 
    377         source = _F64_PRAGMA + source  # Source is already double 
     335        # no need to convert if it is already double 
    378336    elif dtype == F128: 
    379337        fbytes = 16 
    … …  
    418376 
    419377 
    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];\ 
     378_template_cache = {} 
     379def load_template(filename): 
     380    path = joinpath(TEMPLATE_ROOT, filename) 
     381    mtime = getmtime(path) 
     382    if filename not in _template_cache or mtime > _template_cache[filename][0]: 
     383        with open(path) as fid: 
     384            _template_cache[filename] = (mtime, fid.read(), path) 
     385    return _template_cache[filename][1] 
     386 
     387def model_templates(): 
     388    # TODO: fails DRY; templates are listed in two places. 
     389    # should instead have model_info contain a list of paths 
     390    return [joinpath(TEMPLATE_ROOT, filename) 
     391            for filename in ('kernel_header.c', 'kernel_iq.c')] 
     392 
     393 
     394_FN_TEMPLATE = """\ 
     395double %(name)s(%(pars)s); 
     396double %(name)s(%(pars)s) { 
     397    %(body)s 
     398} 
     399 
     400 
    424401""" 
    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 
     402 
     403def _gen_fn(name, pars, body): 
     404    """ 
     405    Generate a function given pars and body. 
     406 
     407    Returns the following string:: 
     408 
     409         double fn(double a, double b, ...); 
     410         double fn(double a, double b, ...) { 
     411             .... 
     412         } 
     413    """ 
     414    par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void' 
     415    return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl} 
     416 
     417def _call_pars(prefix, pars): 
     418    """ 
     419    Return a list of *prefix.parameter* from parameter items. 
     420    """ 
     421    return [p.as_call_reference(prefix) for p in pars] 
     422 
     423_IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)", 
     424                           flags=re.MULTILINE) 
     425def _have_Iqxy(sources): 
     426    """ 
     427    Return true if any file defines Iqxy. 
     428 
     429    Note this is not a C parser, and so can be easily confused by 
     430    non-standard syntax.  Also, it will incorrectly identify the following 
     431    as having Iqxy:: 
     432 
     433        /* 
     434        double Iqxy(qx, qy, ...) { ... fill this in later ... } 
     435        */ 
     436 
     437    If you want to comment out an Iqxy function, use // on the front of the 
     438    line instead. 
     439    """ 
     440    for code in sources: 
     441        if _IQXY_PATTERN.search(code): 
     442            return True 
     443    else: 
     444        return False 
     445 
    444446def make_source(model_info): 
    445447    """ 
    … …  
    460462    # for computing volume even if we allow non-disperse volume parameters. 
    461463 
    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 
     464    partable = model_info['parameters'] 
     465 
     466    # Identify parameters for Iq, Iqxy, Iq_magnetic and form_volume. 
     467    # Note that scale and volume are not possible types. 
     468 
     469    # Load templates and user code 
     470    kernel_header = load_template('kernel_header.c') 
     471    kernel_code = load_template('kernel_iq.c') 
     472    user_code = [open(f).read() for f in model_sources(model_info)] 
     473 
     474    # Build initial sources 
     475    source = [kernel_header] + user_code 
     476 
     477    # Make parameters for q, qx, qy so that we can use them in declarations 
     478    q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 
     479    # Generate form_volume function, etc. from body only 
    497480    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'))) 
     481        pars = partable.form_volume_parameters 
     482        source.append(_gen_fn('form_volume', pars, model_info['form_volume'])) 
    530483    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'))) 
     484        pars = [q] + partable.iq_parameters 
     485        source.append(_gen_fn('Iq', pars, model_info['Iq'])) 
    559486    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         } 
    583  
    584 def categorize_parameters(pars): 
    585     """ 
    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 
     487        pars = [qx, qy] + partable.iqxy_parameters 
     488        source.append(_gen_fn('Iqxy', pars, model_info['Iqxy'])) 
     489 
     490    # Define the parameter table 
     491    source.append("#define PARAMETER_TABLE \\") 
     492    source.append("\\\n".join(p.as_definition() 
     493                              for p in partable.kernel_parameters)) 
     494 
     495    # Define the function calls 
     496    if partable.form_volume_parameters: 
     497        refs = _call_pars("v.", partable.form_volume_parameters) 
     498        call_volume = "#define CALL_VOLUME(v) form_volume(%s)" % (",".join(refs)) 
     499    else: 
     500        # Model doesn't have volume.  We could make the kernel run a little 
     501        # faster by not using/transferring the volume normalizations, but 
     502        # the ifdef's reduce readability more than is worthwhile. 
     503        call_volume = "#define CALL_VOLUME(v) 1.0" 
     504    source.append(call_volume) 
     505 
     506    refs = ["q[i]"] + _call_pars("v.", partable.iq_parameters) 
     507    call_iq = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(refs)) 
     508    if _have_Iqxy(user_code): 
     509        # Call 2D model 
     510        refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("v.", partable.iqxy_parameters) 
     511        call_iqxy = "#define CALL_IQ(q,i,v) Iqxy(%s)" % (",".join(refs)) 
     512    else: 
     513        # Call 1D model with sqrt(qx^2 + qy^2) 
     514        warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 
     515        # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 
     516        pars_sqrt = ["sqrt(q[2*i]*q[2*i]+q[2*i+1]*q[2*i+1])"] + refs[1:] 
     517        call_iqxy = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(pars_sqrt)) 
     518 
     519    # Fill in definitions for numbers of parameters 
     520    source.append("#define MAX_PD %s"%partable.max_pd) 
     521    source.append("#define NPARS %d"%partable.npars) 
     522 
     523    # TODO: allow mixed python/opencl kernels? 
     524 
     525    # define the Iq kernel 
     526    source.append("#define KERNEL_NAME %s_Iq"%model_info['name']) 
     527    source.append(call_iq) 
     528    source.append(kernel_code) 
     529    source.append("#undef CALL_IQ") 
     530    source.append("#undef KERNEL_NAME") 
     531 
     532    # define the Iqxy kernel from the same source with different #defines 
     533    source.append("#define KERNEL_NAME %s_Iqxy"%model_info['name']) 
     534    source.append(call_iqxy) 
     535    source.append(kernel_code) 
     536    source.append("#undef CALL_IQ") 
     537    source.append("#undef KERNEL_NAME") 
     538 
     539    return '\n'.join(source) 
    632540 
    633541def process_parameters(model_info): 
    … …  
    635543    Process parameter block, precalculating parameter details. 
    636544    """ 
    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) 
     545    partable = model_info['parameters'] 
    651546    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'] 
     547        model_info['demo'] = partable.defaults 
     548 
     549class CoordinationDetails(object): 
     550    def __init__(self, model_info): 
     551        parameters = model_info['parameters'] 
     552        max_pd = parameters.max_pd 
     553        npars = parameters.npars 
     554        par_offset = 4*max_pd 
     555        self.details = np.zeros(par_offset + 3*npars + 4, 'i4') 
     556 
     557        # generate views on different parts of the array 
     558        self._pd_par     = self.details[0*max_pd:1*max_pd] 
     559        self._pd_length  = self.details[1*max_pd:2*max_pd] 
     560        self._pd_offset  = self.details[2*max_pd:3*max_pd] 
     561        self._pd_stride  = self.details[3*max_pd:4*max_pd] 
     562        self._par_offset = self.details[par_offset+0*npars:par_offset+1*npars] 
     563        self._par_coord  = self.details[par_offset+1*npars:par_offset+2*npars] 
     564        self._pd_coord   = self.details[par_offset+2*npars:par_offset+3*npars] 
     565 
     566        # theta_par is fixed 
     567        self.details[-1] = parameters.theta_offset 
     568 
     569    @property 
     570    def ctypes(self): return self.details.ctypes 
     571    @property 
     572    def pd_par(self): return self._pd_par 
     573    @property 
     574    def pd_length(self): return self._pd_length 
     575    @property 
     576    def pd_offset(self): return self._pd_offset 
     577    @property 
     578    def pd_stride(self): return self._pd_stride 
     579    @property 
     580    def pd_coord(self): return self._pd_coord 
     581    @property 
     582    def par_coord(self): return self._par_coord 
     583    @property 
     584    def par_offset(self): return self._par_offset 
     585    @property 
     586    def num_coord(self): return self.details[-2] 
     587    @num_coord.setter 
     588    def num_coord(self, v): self.details[-2] = v 
     589    @property 
     590    def total_pd(self): return self.details[-3] 
     591    @total_pd.setter 
     592    def total_pd(self, v): self.details[-3] = v 
     593    @property 
     594    def num_active(self): return self.details[-4] 
     595    @num_active.setter 
     596    def num_active(self, v): self.details[-4] = v 
     597 
     598    def show(self): 
     599        print("total_pd", self.total_pd) 
     600        print("num_active", self.num_active) 
     601        print("pd_par", self.pd_par) 
     602        print("pd_length", self.pd_length) 
     603        print("pd_offset", self.pd_offset) 
     604        print("pd_stride", self.pd_stride) 
     605        print("par_offsets", self.par_offset) 
     606        print("num_coord", self.num_coord) 
     607        print("par_coord", self.par_coord) 
     608        print("pd_coord", self.pd_coord) 
     609        print("theta par", self.details[-1]) 
     610 
     611def mono_details(model_info): 
     612    details = CoordinationDetails(model_info) 
     613    # The zero defaults for monodisperse systems are mostly fine 
     614    details.par_offset[:] = np.arange(2, len(details.par_offset)+2) 
     615    return details 
     616 
     617def poly_details(model_info, weights): 
     618    #print("weights",weights) 
     619    weights = weights[2:] # Skip scale and background 
     620 
     621    # Decreasing list of polydispersity lengths 
     622    # Note: the reversing view, x[::-1], does not require a copy 
     623    pd_length = np.array([len(w) for w in weights]) 
     624    num_active = np.sum(pd_length>1) 
     625    if num_active > model_info['parameters'].max_pd: 
     626        raise ValueError("Too many polydisperse parameters") 
     627 
     628    pd_offset = np.cumsum(np.hstack((0, pd_length))) 
     629    idx = np.argsort(pd_length)[::-1][:num_active] 
     630    par_length = np.array([max(len(w),1) for w in weights]) 
     631    pd_stride = np.cumprod(np.hstack((1, par_length[idx]))) 
     632    par_offsets = np.cumsum(np.hstack((2, par_length))) 
     633 
     634    details = CoordinationDetails(model_info) 
     635    details.pd_par[:num_active] = idx 
     636    details.pd_length[:num_active] = pd_length[idx] 
     637    details.pd_offset[:num_active] = pd_offset[idx] 
     638    details.pd_stride[:num_active] = pd_stride[:-1] 
     639    details.par_offset[:] = par_offsets[:-1] 
     640    details.total_pd = pd_stride[-1] 
     641    details.num_active = num_active 
     642    # Without constraints coordinated parameters are just the pd parameters 
     643    details.par_coord[:num_active] = idx 
     644    details.pd_coord[:num_active] = 2**np.arange(num_active) 
     645    details.num_coord = num_active 
     646    #details.show() 
     647    return details 
     648 
     649def constrained_poly_details(model_info, weights, constraints): 
     650    # Need to find the independently varying pars and sort them 
     651    # Need to build a coordination list for the dependent variables 
     652    # Need to generate a constraints function which takes values 
     653    # and weights, returning par blocks 
     654    raise NotImplementedError("Can't handle constraints yet") 
     655 
     656 
     657def create_default_functions(model_info): 
     658    """ 
     659    Autogenerate missing functions, such as Iqxy from Iq. 
     660 
     661    This only works for Iqxy when Iq is written in python. :func:`make_source` 
     662    performs a similar role for Iq written in C. 
     663    """ 
     664    if callable(model_info['Iq']) and model_info['Iqxy'] is None: 
     665        partable = model_info['parameters'] 
     666        if partable.type['1d'] != partable.type['2d']: 
     667            raise ValueError("Iqxy model is missing") 
     668        Iq = model_info['Iq'] 
     669        def Iqxy(qx, qy, **kw): 
     670            return Iq(np.sqrt(qx**2 + qy**2), **kw) 
     671        model_info['Iqxy'] = Iqxy 
     672 
    654673 
    655674def make_model_info(kernel_module): 
    … …  
    678697      model variants (e.g., the list of cases) or is None if there are no 
    679698      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 
     699    * *par_type* categorizes the model parameters. See 
    685700      :func:`categorize_parameters` for details. 
    686701    * *demo* contains the *{parameter: value}* map used in compare (and maybe 
    … …  
    699714      *model_info* blocks for the composition objects.  This allows us to 
    700715      build complete product and mixture models from just the info. 
    701  
    702716    """ 
    703717    # TODO: maybe turn model_info into a class ModelDefinition 
    704     parameters = COMMON_PARAMETERS + kernel_module.parameters 
     718    #print("make parameter table", kernel_module.parameters) 
     719    parameters = make_parameter_table(kernel_module.parameters) 
    705720    filename = abspath(kernel_module.__file__) 
    706721    kernel_id = splitext(basename(filename))[0] 
    … …  
    712727        filename=abspath(kernel_module.__file__), 
    713728        name=name, 
    714         title=kernel_module.title, 
    715         description=kernel_module.description, 
     729        title=getattr(kernel_module, 'title', name+" model"), 
     730        description=getattr(kernel_module, 'description', 'no description'), 
    716731        parameters=parameters, 
    717732        composition=None, 
    … …  
    731746    functions = "ER VR form_volume Iq Iqxy shape sesans".split() 
    732747    model_info.update((k, getattr(kernel_module, k, None)) for k in functions) 
     748    create_default_functions(model_info) 
     749    # Precalculate the monodisperse parameters 
     750    # TODO: make this a lazy evaluator 
     751    # make_model_info is called for every model on sasview startup 
     752    model_info['mono_details'] = mono_details(model_info) 
    733753    return model_info 
    734754 
    … …  
    782802 
    783803 
    784  
    785804def demo_time(): 
    786805    """ 
    … …  
    798817    Program which prints the source produced by the model. 
    799818    """ 
     819    import sys 
     820    from sasmodels.core import make_model_by_name 
    800821    if len(sys.argv) <= 1: 
    801822        print("usage: python -m sasmodels.generate modelname") 
    802823    else: 
    803824        name = sys.argv[1] 
    804         import sasmodels.models 
    805         __import__('sasmodels.models.' + name) 
    806         model = getattr(sasmodels.models, name) 
    807         model_info = make_model_info(model) 
     825        model_info = make_model_by_name(name) 
    808826        source = make_source(model_info) 
    809827        print(source) 

  • sasmodels/kernelcl.py

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

  • sasmodels/kerneldll.py

    r6ad0e87 rd19962c  
    5050import tempfile 
    5151import ctypes as ct 
    52 from ctypes import c_void_p, c_int, c_longdouble, c_double, c_float 
    53 import _ctypes 
     52from ctypes import c_void_p, c_int32, c_longdouble, c_double, c_float 
    5453 
    5554import numpy as np 
    … …  
    8180        COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 
    8281        if "SAS_OPENMP" in os.environ: 
    83             COMPILE = COMPILE + " -fopenmp" 
     82            COMPILE += " -fopenmp" 
    8483else: 
    8584    COMPILE = "cc -shared -fPIC -fopenmp -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 
    … …  
    9089 
    9190 
    92 def dll_path(model_info, dtype="double"): 
    93     """ 
    94     Path to the compiled model defined by *model_info*. 
    95     """ 
    96     from os.path import join as joinpath, split as splitpath, splitext 
    97     basename = splitext(splitpath(model_info['filename'])[1])[0] 
    98     if np.dtype(dtype) == generate.F32: 
    99         basename += "32" 
    100     elif np.dtype(dtype) == generate.F64: 
    101         basename += "64" 
    102     else: 
    103         basename += "128" 
    104     return joinpath(DLL_PATH, basename+'.so') 
    105  
     91def dll_name(model_info, dtype): 
     92    """ 
     93    Name of the dll containing the model.  This is the base file name without 
     94    any path or extension, with a form such as 'sas_sphere32'. 
     95    """ 
     96    bits = 8*dtype.itemsize 
     97    return "sas_%s%d"%(model_info['id'], bits) 
     98 
     99def dll_path(model_info, dtype): 
     100    """ 
     101    Complete path to the dll for the model.  Note that the dll may not 
     102    exist yet if it hasn't been compiled. 
     103    """ 
     104    return os.path.join(DLL_PATH, dll_name(model_info, dtype)+".so") 
    106105 
    107106def make_dll(source, model_info, dtype="double"): 
    … …  
    133132        dtype = generate.F64  # Force 64-bit dll 
    134133 
    135     if dtype == generate.F32: # 32-bit dll 
    136         tempfile_prefix = 'sas_' + model_info['name'] + '32_' 
    137     elif dtype == generate.F64: 
    138         tempfile_prefix = 'sas_' + model_info['name'] + '64_' 
    139     else: 
    140         tempfile_prefix = 'sas_' + model_info['name'] + '128_' 
    141   
    142134    source = generate.convert_type(source, dtype) 
    143     source_files = generate.model_sources(model_info) + [model_info['filename']] 
     135    newest = generate.timestamp(model_info) 
    144136    dll = dll_path(model_info, dtype) 
    145     newest = max(os.path.getmtime(f) for f in source_files) 
    146137    if not os.path.exists(dll) or os.path.getmtime(dll) < newest: 
    147         # Replace with a proper temp file 
    148         fid, filename = tempfile.mkstemp(suffix=".c", prefix=tempfile_prefix) 
     138        basename = dll_name(model_info, dtype) + "_" 
     139        fid, filename = tempfile.mkstemp(suffix=".c", prefix=basename) 
    149140        os.fdopen(fid, "w").write(source) 
    150141        command = COMPILE%{"source":filename, "output":dll} 
    … …  
    171162    return DllModel(filename, model_info, dtype=dtype) 
    172163 
    173  
    174 IQ_ARGS = [c_void_p, c_void_p, c_int] 
    175 IQXY_ARGS = [c_void_p, c_void_p, c_void_p, c_int] 
    176  
    177164class DllModel(object): 
    178165    """ 
    … …  
    197184 
    198185    def _load_dll(self): 
    199         Nfixed1d = len(self.info['partype']['fixed-1d']) 
    200         Nfixed2d = len(self.info['partype']['fixed-2d']) 
    201         Npd1d = len(self.info['partype']['pd-1d']) 
    202         Npd2d = len(self.info['partype']['pd-2d']) 
    203  
    204186        #print("dll", self.dllpath) 
    205187        try: 
    … …  
    212194              else c_double if self.dtype == generate.F64 
    213195              else c_longdouble) 
    214         pd_args_1d = [c_void_p, fp] + [c_int]*Npd1d if Npd1d else [] 
    215         pd_args_2d = [c_void_p, fp] + [c_int]*Npd2d if Npd2d else [] 
     196 
     197        # int, int, int, int*, double*, double*, double*, double*, double*, double 
     198        argtypes = [c_int32]*3 + [c_void_p]*5 + [fp] 
    216199        self.Iq = self.dll[generate.kernel_name(self.info, False)] 
    217         self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d 
    218  
    219200        self.Iqxy = self.dll[generate.kernel_name(self.info, True)] 
    220         self.Iqxy.argtypes = IQXY_ARGS + pd_args_2d + [fp]*Nfixed2d 
    221          
    222         self.release() 
     201        self.Iq.argtypes = argtypes 
     202        self.Iqxy.argtypes = argtypes 
    223203 
    224204    def __getstate__(self): 
    … …  
    229209        self.dll = None 
    230210 
    231     def __call__(self, q_vectors): 
     211    def make_kernel(self, q_vectors): 
    232212        q_input = PyInput(q_vectors, self.dtype) 
    233213        if self.dll is None: self._load_dll() 
    … …  
    272252    """ 
    273253    def __init__(self, kernel, model_info, q_input): 
     254        self.kernel = kernel 
    274255        self.info = model_info 
    275256        self.q_input = q_input 
    276         self.kernel = kernel 
    277         self.res = np.empty(q_input.nq, q_input.dtype) 
    278         dim = '2d' if q_input.is_2d else '1d' 
    279         self.fixed_pars = model_info['partype']['fixed-' + dim] 
    280         self.pd_pars = model_info['partype']['pd-' + dim] 
    281  
    282         # In dll kernel, but not in opencl kernel 
    283         self.p_res = self.res.ctypes.data 
    284  
    285     def __call__(self, fixed_pars, pd_pars, cutoff): 
     257        self.dtype = q_input.dtype 
     258        self.dim = '2d' if q_input.is_2d else '1d' 
     259        self.result = np.empty(q_input.nq+3, q_input.dtype) 
     260 
     261    def __call__(self, details, weights, values, cutoff): 
    286262        real = (np.float32 if self.q_input.dtype == generate.F32 
    287263                else np.float64 if self.q_input.dtype == generate.F64 
    288264                else np.float128) 
    289  
    290         nq = c_int(self.q_input.nq) 
    291         if pd_pars: 
    292             cutoff = real(cutoff) 
    293             loops_N = [np.uint32(len(p[0])) for p in pd_pars] 
    294             loops = np.hstack(pd_pars) 
    295             loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten() 
    296             p_loops = loops.ctypes.data 
    297             dispersed = [p_loops, cutoff] + loops_N 
    298         else: 
    299             dispersed = [] 
    300         fixed = [real(p) for p in fixed_pars] 
    301         args = self.q_input.q_pointers + [self.p_res, nq] + dispersed + fixed 
    302         #print(pars) 
     265        assert isinstance(details, generate.CoordinationDetails) 
     266        assert weights.dtype == real and values.dtype == real 
     267 
     268        start, stop = 0, details.total_pd 
     269        #print("in kerneldll") 
     270        #print("weights", weights) 
     271        #print("values", values) 
     272        args = [ 
     273            self.q_input.nq, # nq 
     274            start, # pd_start 
     275            stop, # pd_stop pd_stride[MAX_PD] 
     276            details.ctypes.data, # problem 
     277            weights.ctypes.data,  # weights 
     278            values.ctypes.data,  #pars 
     279            self.q_input.q.ctypes.data, #q 
     280            self.result.ctypes.data,   # results 
     281            real(cutoff), # cutoff 
     282            ] 
    303283        self.kernel(*args) 
    304  
    305         return self.res 
     284        return self.result[:-3] 
    306285 
    307286    def release(self): 

  • sasmodels/kernelpy.py

    ra84a0ca r48fbd50  
    5353        self.dtype = dtype 
    5454        self.is_2d = (len(q_vectors) == 2) 
    55         self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors] 
    56         self.q_pointers = [q.ctypes.data for q in self.q_vectors] 
     55        if self.is_2d: 
     56            self.q = np.empty((self.nq, 2), dtype=dtype) 
     57            self.q[:, 0] = q_vectors[0] 
     58            self.q[:, 1] = q_vectors[1] 
     59        else: 
     60            self.q = np.empty(self.nq, dtype=dtype) 
     61            self.q[:self.nq] = q_vectors[0] 
    5762 
    5863    def release(self): 
    … …  
    6065        Free resources associated with the model inputs. 
    6166        """ 
    62         self.q_vectors = [] 
     67        self.q = None 
    6368 
    6469class PyKernel(object): 

  • sasmodels/mixture.py

    r72a081d r69aa451  
    1414import numpy as np 
    1515 
    16 from .generate import process_parameters, COMMON_PARAMETERS, Parameter 
     16from .modelinfo import Parameter, ParameterTable 
    1717 
    1818SCALE=0 
    … …  
    3434 
    3535    # Build new parameter list 
    36     pars = COMMON_PARAMETERS + [] 
     36    pars = [] 
    3737    for k, part in enumerate(parts): 
    3838        # Parameter prefix per model, A_, B_, ... 
     39        # Note that prefix must also be applied to id and length_control 
     40        # to support vector parameters 
    3941        prefix = chr(ord('A')+k) + '_' 
    40         for p in part['parameters']: 
    41             # No background on the individual mixture elements 
    42             if p.name == 'background': 
    43                 continue 
    44             # TODO: promote Parameter to a full class 
    45             # this code knows too much about the implementation! 
    46             p = list(p) 
    47             if p[0] == 'scale':  # allow model subtraction 
    48                 p[3] = [-np.inf, np.inf]  # limits 
    49             p[0] = prefix+p[0]   # relabel parameters with A_, ... 
     42        pars.append(Parameter(prefix+'scale')) 
     43        for p in part['parameters'].kernel_pars: 
     44            p = copy(p) 
     45            p.name = prefix+p.name 
     46            p.id = prefix+p.id 
     47            if p.length_control is not None: 
     48                p.length_control = prefix+p.length_control 
    5049            pars.append(p) 
     50    partable = ParameterTable(pars) 
    5151 
    5252    model_info = {} 
    … …  
    5858    model_info['docs'] = model_info['title'] 
    5959    model_info['category'] = "custom" 
    60     model_info['parameters'] = pars 
     60    model_info['parameters'] = partable 
    6161    #model_info['single'] = any(part['single'] for part in parts) 
    6262    model_info['structure_factor'] = False 
    … …  
    6767    # Remember the component info blocks so we can build the model 
    6868    model_info['composition'] = ('mixture', parts) 
    69     process_parameters(model_info) 
    70     return model_info 
    7169 
    7270 

  • sasmodels/model_test.py

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

  • sasmodels/models/cylinder.c

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

  • sasmodels/models/gel_fit.c

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

  • sasmodels/models/lib/Si.c

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

  • sasmodels/models/lib/polevl.c

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

  • sasmodels/models/lib/sas_J0.c

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

  • sasmodels/models/lib/sas_J1.c

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

  • sasmodels/models/lib/sas_JN.c

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

  • sasmodels/models/lib/sph_j1c.c

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

  • sasmodels/models/lib/sphere_form.c

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

  • sasmodels/models/lib/wrc_cyl.c

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

  • sasmodels/models/rpa.py

    raad336c re9b1663d  
    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 re9b1663d  
    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 re9b1663d  
    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.