Changeset 3fb3449 in sasmodels


Ignore:
Timestamp:
Apr 6, 2016 3:42:25 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:
398aa94
Parents:
21b116f (diff), c4e7a5f (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

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

Files:
5 added
42 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/kernel_template.c

    rcf85329 rc4e7a5f  
    158158#ifdef IQ_OPEN_LOOPS 
    159159    double ret=0.0, norm=0.0; 
    160   #ifdef VOLUME_PARAMETERS 
    161     double vol=0.0, norm_vol=0.0; 
    162   #endif 
    163160    IQ_OPEN_LOOPS 
    164161    //for (int radius_i=0; radius_i < Nradius; radius_i++) { 
    … …  
    172169      if (!isnan(scattering)) { 
    173170        ret += weight*scattering; 
     171      #ifdef VOLUME_PARAMETERS 
     172        norm += weight * form_volume(VOLUME_PARAMETERS); 
     173      #else 
    174174        norm += weight; 
    175       #ifdef VOLUME_PARAMETERS 
    176         const double vol_weight = VOLUME_WEIGHT_PRODUCT; 
    177         vol += vol_weight*form_volume(VOLUME_PARAMETERS); 
    178         norm_vol += vol_weight; 
    179175      #endif 
    180176      } 
    … …  
    182178    } 
    183179    IQ_CLOSE_LOOPS 
    184   #ifdef VOLUME_PARAMETERS 
    185     if (vol*norm_vol != 0.0) { 
    186       ret *= norm_vol/vol; 
    187     } 
    188   #endif 
    189     result[i] = scale*ret/norm+background; 
     180    // norm can only be zero if volume is zero, so no scattering 
     181    result[i] = (norm > 0. ? scale*ret/norm + background : background); 
    190182#else 
    191183    result[i] = scale*Iq(qi, IQ_PARAMETERS) + background; 
    … …  
    241233    //  const double radius = loops[2*(radius_i)]; 
    242234    //  const double radius_w = loops[2*(radius_i)+1]; 
    243  
    244     const double weight = IQXY_WEIGHT_PRODUCT; 
     235    double weight = IQXY_WEIGHT_PRODUCT; 
    245236    if (weight > cutoff) { 
    246237 
    247238      const double scattering = Iqxy(qxi, qyi, IQXY_PARAMETERS); 
    248239      if (!isnan(scattering)) { // if scattering is bad, exclude it from sum 
    249       //if (scattering >= 0.0) { // scattering cannot be negative 
    250         // TODO: use correct angle for spherical correction 
    251         // Definition of theta and phi are probably reversed relative to the 
    252         // equation which gave rise to this correction, leading to an 
    253         // attenuation of the pattern as theta moves through pi/2.  Either 
    254         // reverse the meanings of phi and theta in the forms, or use phi 
    255         // rather than theta in this correction.  Current code uses cos(theta) 
    256         // so that values match those of sasview. 
    257       #if defined(IQXY_HAS_THETA) // && 0 
    258         const double spherical_correction 
    259           = (Ntheta>1 ? fabs(cos(M_PI_180*theta))*M_PI_2:1.0); 
    260         const double next = spherical_correction * weight * scattering; 
    261       #else 
    262         const double next = weight * scattering; 
    263       #endif 
     240      #if defined(IQXY_HAS_THETA) 
     241        // Force a nominal value for the spherical correction even when 
     242        // theta is +90/-90 so that there are no divide by zero problems. 
     243        // For cos(theta) fixed at 90, we effectively multiply top and bottom 
     244        // by 1e-6, so the effect cancels. 
     245        const double spherical_correction = fmax(fabs(cos(M_PI_180*theta)), 1e-6); 
     246        weight *= spherical_correction; 
     247      #endif 
     248      const double next = weight * scattering; 
    264249      #if USE_KAHAN_SUMMATION 
    265250        const double y = next - accumulated_error; 
    … …  
    270255        ret += next; 
    271256      #endif 
     257      #ifdef VOLUME_PARAMETERS 
     258        norm += weight*form_volume(VOLUME_PARAMETERS); 
     259      #else 
    272260        norm += weight; 
    273       #ifdef VOLUME_PARAMETERS 
    274         const double vol_weight = VOLUME_WEIGHT_PRODUCT; 
    275         vol += vol_weight*form_volume(VOLUME_PARAMETERS); 
    276         norm_vol += vol_weight; 
    277261      #endif 
    278262      } 
    … …  
    280264    } 
    281265    IQXY_CLOSE_LOOPS 
    282   #ifdef VOLUME_PARAMETERS 
    283     if (vol*norm_vol != 0.0) { 
    284       ret *= norm_vol/vol; 
    285     } 
    286   #endif 
    287     result[i] = scale*ret/norm+background; 
     266    // norm can only be zero if volume is zero, so no scattering 
     267    result[i] = (norm>0. ? scale*ret/norm + background : background); 
    288268#else 
    289269    result[i] = scale*Iqxy(qxi, qyi, IQXY_PARAMETERS) + background; 

  • doc/developer/index.rst

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

  • sasmodels/bumps_model.py

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

  • sasmodels/compare.py

    rf247314 r21b116f  
    3838from . import core 
    3939from . import kerneldll 
    40 from . import product 
    4140from .data import plot_theory, empty_data1D, empty_data2D 
    4241from .direct_model import DirectModel 
    … …  
    306305    Format the parameter list for printing. 
    307306    """ 
    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 
    314307    lines = [] 
    315     for p in model_info['parameters']: 
    316         if exclude(p.name): continue 
     308    parameters = model_info['parameters'] 
     309    for p in parameters.user_parameters(pars, is2d): 
    317310        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')) 
     311            value=pars.get(p.id, p.default), 
     312            pd=pars.get(p.id+"_pd", 0.), 
     313            n=int(pars.get(p.id+"_pd_n", 0)), 
     314            nsigma=pars.get(p.id+"_pd_nsgima", 3.), 
     315            type=pars.get(p.id+"_pd_type", 'gaussian')) 
    323316        lines.append(_format_par(p.name, **fields)) 
    324317    return "\n".join(lines) 
    … …  
    454447    """ 
    455448    # initialize the code so time is more accurate 
    456     value = calculator(**suppress_pd(pars)) 
     449    if Nevals > 1: 
     450        value = calculator(**suppress_pd(pars)) 
    457451    toc = tic() 
    458452    for _ in range(max(Nevals, 1)):  # make sure there is at least one eval 
    … …  
    661655    """ 
    662656    # 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" 
     657    pars = {} 
     658    for p in model_info['parameters'].call_parameters: 
     659        parts = [('', p.default)] 
     660        if p.polydisperse: 
     661            parts.append(('_pd', 0.0)) 
     662            parts.append(('_pd_n', 0)) 
     663            parts.append(('_pd_nsigma', 3.0)) 
     664            parts.append(('_pd_type', "gaussian")) 
     665        for ext, val in parts: 
     666            if p.length > 1: 
     667                dict(("%s%d%s"%(p.id,k,ext), val) for k in range(p.length)) 
     668            else: 
     669                pars[p.id+ext] = val 
    672670 
    673671    # Plug in values given in demo 
    … …  
    774772 
    775773    if len(engines) == 0: 
    776         engines.extend(['single', 'sasview']) 
     774        engines.extend(['single', 'double']) 
    777775    elif len(engines) == 1: 
    778         if engines[0][0] != 'sasview': 
    779             engines.append('sasview') 
     776        if engines[0][0] != 'double': 
     777            engines.append('double') 
    780778        else: 
    781779            engines.append('single') 
    … …  
    879877        model_info = opts['def'] 
    880878        pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars']) 
     879        # Initialize parameter ranges, fixing the 2D parameters for 1D data. 
    881880        if not opts['is2d']: 
    882             active = [base + ext 
    883                       for base in model_info['partype']['pd-1d'] 
    884                       for ext in ['', '_pd', '_pd_n', '_pd_nsigma']] 
    885             active.extend(model_info['partype']['fixed-1d']) 
    886             for k in active: 
    887                 v = pars[k] 
    888                 v.range(*parameter_range(k, v.value)) 
     881            for p in model_info['parameters'].user_parameters(is2d=False): 
     882                for ext in ['', '_pd', '_pd_n', '_pd_nsigma']: 
     883                    k = p.name+ext 
     884                    v = pars.get(k, None) 
     885                    if v is not None: 
     886                        v.range(*parameter_range(k, v.value)) 
    889887        else: 
    890888            for k, v in pars.items(): 

  • sasmodels/convert.json

    rec45c4f rd2bb604  
    623623    "SphereSLDModel",  
    624624    { 
     625      "n": "n_shells", 
    625626      "radius_core": "rad_core",  
    626627      "sld_solvent": "sld_solv" 

  • sasmodels/core.py

    r4d76711 rd2bb604  
    170170 
    171171 
    172 def get_weights(model_info, pars, name): 
     172def get_weights(parameter, values): 
    173173    """ 
    174174    Generate the distribution for parameter *name* given the parameter values 
    … …  
    178178    from the *pars* dictionary for parameter value and parameter dispersion. 
    179179    """ 
    180     relative = name in model_info['partype']['pd-rel'] 
    181     limits = model_info['limits'][name] 
    182     disperser = pars.get(name+'_pd_type', 'gaussian') 
    183     value = pars.get(name, model_info['defaults'][name]) 
    184     npts = pars.get(name+'_pd_n', 0) 
    185     width = pars.get(name+'_pd', 0.0) 
    186     nsigma = pars.get(name+'_pd_nsigma', 3.0) 
     180    value = float(values.get(parameter.name, parameter.default)) 
     181    relative = parameter.relative_pd 
     182    limits = parameter.limits 
     183    disperser = values.get(parameter.name+'_pd_type', 'gaussian') 
     184    npts = values.get(parameter.name+'_pd_n', 0) 
     185    width = values.get(parameter.name+'_pd', 0.0) 
     186    nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 
     187    if npts == 0 or width == 0: 
     188        return [value], [] 
    187189    value, weight = weights.get_weights( 
    188190        disperser, npts, width, nsigma, value, limits, relative) 
    … …  
    198200    """ 
    199201    value, weight = zip(*pars) 
     202    weight = [w if w else [1.] for w in weight] 
    200203    value = [v.flatten() for v in meshgrid(*value)] 
    201204    weight = np.vstack([v.flatten() for v in meshgrid(*weight)]) 
    … …  
    216219    *mono* is True if polydispersity should be set to none on all parameters. 
    217220    """ 
    218     fixed_pars = [pars.get(name, kernel.info['defaults'][name]) 
    219                   for name in kernel.fixed_pars] 
     221    parameters = kernel.info['parameters'] 
    220222    if mono: 
    221         pd_pars = [( np.array([pars[name]]), np.array([1.0]) ) 
    222                    for name in kernel.pd_pars] 
    223     else: 
    224         pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars] 
    225     return kernel(fixed_pars, pd_pars, cutoff=cutoff) 
     223        active = lambda name: False 
     224    elif kernel.dim == '1d': 
     225        active = lambda name: name in parameters.pd_1d 
     226    elif kernel.dim == '2d': 
     227        active = lambda name: name in parameters.pd_2d 
     228    else: 
     229        active = lambda name: True 
     230 
     231    vw_pairs = [(get_weights(p, pars) if active(p.name) 
     232                 else ([pars.get(p.name, p.default)], [])) 
     233                for p in parameters.call_parameters] 
     234 
     235    details, weights, values = build_details(kernel, vw_pairs) 
     236    return kernel(details, weights, values, cutoff) 
     237 
     238def build_details(kernel, pairs): 
     239    values, weights = zip(*pairs) 
     240    if max([len(w) for w in weights]) > 1: 
     241        details = generate.poly_details(kernel.info, weights) 
     242    else: 
     243        details = kernel.info['mono_details'] 
     244    weights, values = [np.hstack(v) for v in (weights, values)] 
     245    weights = weights.astype(dtype=kernel.dtype) 
     246    values = values.astype(dtype=kernel.dtype) 
     247    return details, weights, values 
     248 
    226249 
    227250def call_ER_VR(model_info, vol_pars): 
    … …  
    256279        return 1.0 
    257280    else: 
    258         vol_pars = [get_weights(model_info, values, name) 
    259                     for name in model_info['partype']['volume']] 
     281        vol_pars = [get_weights(parameter, values) 
     282                    for parameter in model_info['parameters'].call_parameters 
     283                    if parameter.type == 'volume'] 
    260284        value, weight = dispersion_mesh(vol_pars) 
    261285        individual_radii = ER(*value) 
    262         #print(values[0].shape, weights.shape, fv.shape) 
    263286        return np.sum(weight*individual_radii) / np.sum(weight) 
    264287 
    … …  
    273296        return 1.0 
    274297    else: 
    275         vol_pars = [get_weights(model_info, values, name) 
    276                     for name in model_info['partype']['volume']] 
     298        vol_pars = [get_weights(parameter, values) 
     299                    for parameter in model_info['parameters'].call_parameters 
     300                    if parameter.type == 'volume'] 
    277301        value, weight = dispersion_mesh(vol_pars) 
    278302        whole, part = VR(*value) 

  • sasmodels/direct_model.py

    r4d76711 r68e7f9d  
    6767            self.data_type = 'Iq' 
    6868 
    69         partype = model.info['partype'] 
    70  
    7169        if self.data_type == 'sesans': 
    7270             
    … …  
    8280            q_mono = sesans.make_all_q(data) 
    8381        elif self.data_type == 'Iqxy': 
    84             if not partype['orientation'] and not partype['magnetic']: 
    85                 raise ValueError("not 2D without orientation or magnetic parameters") 
     82            #if not model.info['parameters'].has_2d: 
     83            #    raise ValueError("not 2D without orientation or magnetic parameters") 
    8684            q = np.sqrt(data.qx_data**2 + data.qy_data**2) 
    8785            qmin = getattr(data, 'qmin', 1e-16) 
    … …  
    172170    def _calc_theory(self, pars, cutoff=0.0): 
    173171        if self._kernel is None: 
    174             self._kernel = self._model.make_kernel(self._kernel_inputs)  # pylint: disable=attribute-dedata_type 
     172            self._kernel = self._model.make_kernel(self._kernel_inputs) 
    175173            self._kernel_mono = (self._model.make_kernel(self._kernel_mono_inputs) 
    176174                                 if self._kernel_mono_inputs else None) 

  • sasmodels/generate.py

    r4d76711 rd2bb604  
    2121 
    2222    *VR(p1, p2, ...)* returns the volume ratio for core-shell style forms. 
     23 
     24    #define INVALID(v) (expr)  returns False if v.parameter is invalid 
     25    for some parameter or other (e.g., v.bell_radius < v.radius).  If 
     26    necessary, the expression can call a function. 
    2327 
    2428These functions are defined in a kernel module .py script and an associated 
    … …  
    6569The constructor code will generate the necessary vectors for computing 
    6670them with the desired polydispersity. 
    67  
    68 The available kernel parameters are defined as a list, with each parameter 
    69 defined as a sublist with the following elements: 
    70  
    71     *name* is the name that will be used in the call to the kernel 
    72     function and the name that will be displayed to the user.  Names 
    73     should be lower case, with words separated by underscore.  If 
    74     acronyms are used, the whole acronym should be upper case. 
    75  
    76     *units* should be one of *degrees* for angles, *Ang* for lengths, 
    77     *1e-6/Ang^2* for SLDs. 
    78  
    79     *default value* will be the initial value for  the model when it 
    80     is selected, or when an initial value is not otherwise specified. 
    81  
    82     *limits = [lb, ub]* are the hard limits on the parameter value, used to 
    83     limit the polydispersity density function.  In the fit, the parameter limits 
    84     given to the fit are the limits  on the central value of the parameter. 
    85     If there is polydispersity, it will evaluate parameter values outside 
    86     the fit limits, but not outside the hard limits specified in the model. 
    87     If there are no limits, use +/-inf imported from numpy. 
    88  
    89     *type* indicates how the parameter will be used.  "volume" parameters 
    90     will be used in all functions.  "orientation" parameters will be used 
    91     in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in 
    92     *Imagnetic* only.  If *type* is the empty string, the parameter will 
    93     be used in all of *Iq*, *Iqxy* and *Imagnetic*.  "sld" parameters 
    94     can automatically be promoted to magnetic parameters, each of which 
    95     will have a magnitude and a direction, which may be different from 
    96     other sld parameters. 
    97  
    98     *description* is a short description of the parameter.  This will 
    99     be displayed in the parameter table and used as a tool tip for the 
    100     parameter value in the user interface. 
    101  
    10271The kernel module must set variables defining the kernel meta data: 
    10372 
    … …  
    204173from __future__ import print_function 
    205174 
    206 #TODO: identify model files which have changed since loading and reload them. 
    207175#TODO: determine which functions are useful outside of generate 
    208176#__all__ = ["model_info", "make_doc", "make_source", "convert_type"] 
    209177 
    210 import sys 
    211178from os.path import abspath, dirname, join as joinpath, exists, basename, \ 
    212     splitext 
     179    splitext, getmtime 
    213180import re 
    214181import string 
    215182import warnings 
    216 from collections import namedtuple 
    217183 
    218184import numpy as np 
    219185 
     186from .modelinfo import ModelInfo, Parameter, make_parameter_table, set_demo 
    220187from .custom import load_custom_kernel_module 
    221188 
    222 PARAMETER_FIELDS = ['name', 'units', 'default', 'limits', 'type', 'description'] 
    223 Parameter = namedtuple('Parameter', PARAMETER_FIELDS) 
    224  
    225 C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c') 
     189TEMPLATE_ROOT = dirname(__file__) 
    226190 
    227191F16 = np.dtype('float16') 
    … …  
    232196except TypeError: 
    233197    F128 = None 
    234  
    235 # Scale and background, which are parameters common to every form factor 
    236 COMMON_PARAMETERS = [ 
    237     ["scale", "", 1, [0, np.inf], "", "Source intensity"], 
    238     ["background", "1/cm", 1e-3, [0, np.inf], "", "Source background"], 
    239     ] 
    240198 
    241199# Conversion from units defined in the parameter table for each model 
    … …  
    331289    raise ValueError("%r not found in %s" % (filename, search_path)) 
    332290 
     291 
    333292def model_sources(model_info): 
    334293    """ 
    … …  
    339298    return [_search(search_path, f) for f in model_info['source']] 
    340299 
    341 # Pragmas for enable OpenCL features.  Be sure to protect them so that they 
    342 # still compile even if OpenCL is not present. 
    343 _F16_PRAGMA = """\ 
    344 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 
    345 #  pragma OPENCL EXTENSION cl_khr_fp16: enable 
    346 #endif 
    347 """ 
    348  
    349 _F64_PRAGMA = """\ 
    350 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 
    351 #  pragma OPENCL EXTENSION cl_khr_fp64: enable 
    352 #endif 
    353 """ 
     300def timestamp(model_info): 
     301    """ 
     302    Return a timestamp for the model corresponding to the most recently 
     303    changed file or dependency. 
     304    """ 
     305    source_files = (model_sources(model_info) 
     306                    + model_templates() 
     307                    + [model_info['filename']]) 
     308    newest = max(getmtime(f) for f in source_files) 
     309    return newest 
    354310 
    355311def convert_type(source, dtype): 
    … …  
    362318    if dtype == F16: 
    363319        fbytes = 2 
    364         source = _F16_PRAGMA + _convert_type(source, "half", "f") 
     320        source = _convert_type(source, "half", "f") 
    365321    elif dtype == F32: 
    366322        fbytes = 4 
    … …  
    368324    elif dtype == F64: 
    369325        fbytes = 8 
    370         source = _F64_PRAGMA + source  # Source is already double 
     326        # no need to convert if it is already double 
    371327    elif dtype == F128: 
    372328        fbytes = 16 
    … …  
    411367 
    412368 
    413 LOOP_OPEN = """\ 
    414 for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) { 
    415   const double %(name)s = loops[2*(%(name)s_i%(offset)s)]; 
    416   const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\ 
     369_template_cache = {} 
     370def load_template(filename): 
     371    path = joinpath(TEMPLATE_ROOT, filename) 
     372    mtime = getmtime(path) 
     373    if filename not in _template_cache or mtime > _template_cache[filename][0]: 
     374        with open(path) as fid: 
     375            _template_cache[filename] = (mtime, fid.read(), path) 
     376    return _template_cache[filename][1] 
     377 
     378def model_templates(): 
     379    # TODO: fails DRY; templates are listed in two places. 
     380    # should instead have model_info contain a list of paths 
     381    return [joinpath(TEMPLATE_ROOT, filename) 
     382            for filename in ('kernel_header.c', 'kernel_iq.c')] 
     383 
     384 
     385_FN_TEMPLATE = """\ 
     386double %(name)s(%(pars)s); 
     387double %(name)s(%(pars)s) { 
     388    %(body)s 
     389} 
     390 
     391 
    417392""" 
    418 def build_polydispersity_loops(pd_pars): 
    419     """ 
    420     Build polydispersity loops 
    421  
    422     Returns loop opening and loop closing 
    423     """ 
    424     depth = 4 
    425     offset = "" 
    426     loop_head = [] 
    427     loop_end = [] 
    428     for name in pd_pars: 
    429         subst = {'name': name, 'offset': offset} 
    430         loop_head.append(indent(LOOP_OPEN % subst, depth)) 
    431         loop_end.insert(0, (" "*depth) + "}") 
    432         offset += '+N' + name 
    433         depth += 2 
    434     return "\n".join(loop_head), "\n".join(loop_end) 
    435  
    436 C_KERNEL_TEMPLATE = None 
     393 
     394def _gen_fn(name, pars, body): 
     395    """ 
     396    Generate a function given pars and body. 
     397 
     398    Returns the following string:: 
     399 
     400         double fn(double a, double b, ...); 
     401         double fn(double a, double b, ...) { 
     402             .... 
     403         } 
     404    """ 
     405    par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void' 
     406    return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl} 
     407 
     408def _call_pars(prefix, pars): 
     409    """ 
     410    Return a list of *prefix.parameter* from parameter items. 
     411    """ 
     412    return [p.as_call_reference(prefix) for p in pars] 
     413 
     414_IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)", 
     415                           flags=re.MULTILINE) 
     416def _have_Iqxy(sources): 
     417    """ 
     418    Return true if any file defines Iqxy. 
     419 
     420    Note this is not a C parser, and so can be easily confused by 
     421    non-standard syntax.  Also, it will incorrectly identify the following 
     422    as having Iqxy:: 
     423 
     424        /* 
     425        double Iqxy(qx, qy, ...) { ... fill this in later ... } 
     426        */ 
     427 
     428    If you want to comment out an Iqxy function, use // on the front of the 
     429    line instead. 
     430    """ 
     431    for code in sources: 
     432        if _IQXY_PATTERN.search(code): 
     433            return True 
     434    else: 
     435        return False 
     436 
    437437def make_source(model_info): 
    438438    """ 
    … …  
    453453    # for computing volume even if we allow non-disperse volume parameters. 
    454454 
    455     # Load template 
    456     global C_KERNEL_TEMPLATE 
    457     if C_KERNEL_TEMPLATE is None: 
    458         with open(C_KERNEL_TEMPLATE_PATH) as fid: 
    459             C_KERNEL_TEMPLATE = fid.read() 
    460  
    461     # Load additional sources 
    462     source = [open(f).read() for f in model_sources(model_info)] 
    463  
    464     # Prepare defines 
    465     defines = [] 
    466     partype = model_info['partype'] 
    467     pd_1d = partype['pd-1d'] 
    468     pd_2d = partype['pd-2d'] 
    469     fixed_1d = partype['fixed-1d'] 
    470     fixed_2d = partype['fixed-1d'] 
    471  
    472     iq_parameters = [p.name 
    473                      for p in model_info['parameters'][2:]  # skip scale, background 
    474                      if p.name in set(fixed_1d + pd_1d)] 
    475     iqxy_parameters = [p.name 
    476                        for p in model_info['parameters'][2:]  # skip scale, background 
    477                        if p.name in set(fixed_2d + pd_2d)] 
    478     volume_parameters = [p.name 
    479                          for p in model_info['parameters'] 
    480                          if p.type == 'volume'] 
    481  
    482     # Fill in defintions for volume parameters 
    483     if volume_parameters: 
    484         defines.append(('VOLUME_PARAMETERS', 
    485                         ','.join(volume_parameters))) 
    486         defines.append(('VOLUME_WEIGHT_PRODUCT', 
    487                         '*'.join(p + '_w' for p in volume_parameters))) 
    488  
    489     # Generate form_volume function from body only 
     455    partable = model_info['parameters'] 
     456 
     457    # Identify parameters for Iq, Iqxy, Iq_magnetic and form_volume. 
     458    # Note that scale and volume are not possible types. 
     459 
     460    # Load templates and user code 
     461    kernel_header = load_template('kernel_header.c') 
     462    kernel_code = load_template('kernel_iq.c') 
     463    user_code = [open(f).read() for f in model_sources(model_info)] 
     464 
     465    # Build initial sources 
     466    source = [kernel_header] + user_code 
     467 
     468    # Make parameters for q, qx, qy so that we can use them in declarations 
     469    q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 
     470    # Generate form_volume function, etc. from body only 
    490471    if model_info['form_volume'] is not None: 
    491         if volume_parameters: 
    492             vol_par_decl = ', '.join('double ' + p for p in volume_parameters) 
    493         else: 
    494             vol_par_decl = 'void' 
    495         defines.append(('VOLUME_PARAMETER_DECLARATIONS', 
    496                         vol_par_decl)) 
    497         fn = """\ 
    498 double form_volume(VOLUME_PARAMETER_DECLARATIONS); 
    499 double form_volume(VOLUME_PARAMETER_DECLARATIONS) { 
    500     %(body)s 
    501 } 
    502 """ % {'body':model_info['form_volume']} 
    503         source.append(fn) 
    504  
    505     # Fill in definitions for Iq parameters 
    506     defines.append(('IQ_KERNEL_NAME', model_info['name'] + '_Iq')) 
    507     defines.append(('IQ_PARAMETERS', ', '.join(iq_parameters))) 
    508     if fixed_1d: 
    509         defines.append(('IQ_FIXED_PARAMETER_DECLARATIONS', 
    510                         ', \\\n    '.join('const double %s' % p for p in fixed_1d))) 
    511     if pd_1d: 
    512         defines.append(('IQ_WEIGHT_PRODUCT', 
    513                         '*'.join(p + '_w' for p in pd_1d))) 
    514         defines.append(('IQ_DISPERSION_LENGTH_DECLARATIONS', 
    515                         ', \\\n    '.join('const int N%s' % p for p in pd_1d))) 
    516         defines.append(('IQ_DISPERSION_LENGTH_SUM', 
    517                         '+'.join('N' + p for p in pd_1d))) 
    518         open_loops, close_loops = build_polydispersity_loops(pd_1d) 
    519         defines.append(('IQ_OPEN_LOOPS', 
    520                         open_loops.replace('\n', ' \\\n'))) 
    521         defines.append(('IQ_CLOSE_LOOPS', 
    522                         close_loops.replace('\n', ' \\\n'))) 
     472        pars = partable.form_volume_parameters 
     473        source.append(_gen_fn('form_volume', pars, model_info['form_volume'])) 
    523474    if model_info['Iq'] is not None: 
    524         defines.append(('IQ_PARAMETER_DECLARATIONS', 
    525                         ', '.join('double ' + p for p in iq_parameters))) 
    526         fn = """\ 
    527 double Iq(double q, IQ_PARAMETER_DECLARATIONS); 
    528 double Iq(double q, IQ_PARAMETER_DECLARATIONS) { 
    529     %(body)s 
    530 } 
    531 """ % {'body':model_info['Iq']} 
    532         source.append(fn) 
    533  
    534     # Fill in definitions for Iqxy parameters 
    535     defines.append(('IQXY_KERNEL_NAME', model_info['name'] + '_Iqxy')) 
    536     defines.append(('IQXY_PARAMETERS', ', '.join(iqxy_parameters))) 
    537     if fixed_2d: 
    538         defines.append(('IQXY_FIXED_PARAMETER_DECLARATIONS', 
    539                         ', \\\n    '.join('const double %s' % p for p in fixed_2d))) 
    540     if pd_2d: 
    541         defines.append(('IQXY_WEIGHT_PRODUCT', 
    542                         '*'.join(p + '_w' for p in pd_2d))) 
    543         defines.append(('IQXY_DISPERSION_LENGTH_DECLARATIONS', 
    544                         ', \\\n    '.join('const int N%s' % p for p in pd_2d))) 
    545         defines.append(('IQXY_DISPERSION_LENGTH_SUM', 
    546                         '+'.join('N' + p for p in pd_2d))) 
    547         open_loops, close_loops = build_polydispersity_loops(pd_2d) 
    548         defines.append(('IQXY_OPEN_LOOPS', 
    549                         open_loops.replace('\n', ' \\\n'))) 
    550         defines.append(('IQXY_CLOSE_LOOPS', 
    551                         close_loops.replace('\n', ' \\\n'))) 
     475        pars = [q] + partable.iq_parameters 
     476        source.append(_gen_fn('Iq', pars, model_info['Iq'])) 
    552477    if model_info['Iqxy'] is not None: 
    553         defines.append(('IQXY_PARAMETER_DECLARATIONS', 
    554                         ', '.join('double ' + p for p in iqxy_parameters))) 
    555         fn = """\ 
    556 double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS); 
    557 double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS) { 
    558     %(body)s 
    559 } 
    560 """ % {'body':model_info['Iqxy']} 
    561         source.append(fn) 
    562  
    563     # Need to know if we have a theta parameter for Iqxy; it is not there 
    564     # for the magnetic sphere model, for example, which has a magnetic 
    565     # orientation but no shape orientation. 
    566     if 'theta' in pd_2d: 
    567         defines.append(('IQXY_HAS_THETA', '1')) 
    568  
    569     #for d in defines: print(d) 
    570     defines = '\n'.join('#define %s %s' % (k, v) for k, v in defines) 
    571     sources = '\n\n'.join(source) 
    572     return C_KERNEL_TEMPLATE % { 
    573         'DEFINES': defines, 
    574         'SOURCES': sources, 
    575         } 
    576  
    577 def categorize_parameters(pars): 
    578     """ 
    579     Build parameter categories out of the the parameter definitions. 
    580  
    581     Returns a dictionary of categories. 
    582  
    583     Note: these categories are subject to change, depending on the needs of 
    584     the UI and the needs of the kernel calling function. 
    585  
    586     The categories are as follows: 
    587  
    588     * *volume* list of volume parameter names 
    589     * *orientation* list of orientation parameters 
    590     * *magnetic* list of magnetic parameters 
    591     * *<empty string>* list of parameters that have no type info 
    592  
    593     Each parameter is in one and only one category. 
    594  
    595     The following derived categories are created: 
    596  
    597     * *fixed-1d* list of non-polydisperse parameters for 1D models 
    598     * *pd-1d* list of polydisperse parameters for 1D models 
    599     * *fixed-2d* list of non-polydisperse parameters for 2D models 
    600     * *pd-d2* list of polydisperse parameters for 2D models 
    601     """ 
    602     partype = { 
    603         'volume': [], 'orientation': [], 'magnetic': [], 'sld': [], '': [], 
    604         'fixed-1d': [], 'fixed-2d': [], 'pd-1d': [], 'pd-2d': [], 
    605         'pd-rel': set(), 
    606     } 
    607  
    608     for p in pars: 
    609         if p.type == 'volume': 
    610             partype['pd-1d'].append(p.name) 
    611             partype['pd-2d'].append(p.name) 
    612             partype['pd-rel'].add(p.name) 
    613         elif p.type == 'magnetic': 
    614             partype['fixed-2d'].append(p.name) 
    615         elif p.type == 'orientation': 
    616             partype['pd-2d'].append(p.name) 
    617         elif p.type in ('', 'sld'): 
    618             partype['fixed-1d'].append(p.name) 
    619             partype['fixed-2d'].append(p.name) 
    620         else: 
    621             raise ValueError("unknown parameter type %r" % p.type) 
    622         partype[p.type].append(p.name) 
    623  
    624     return partype 
    625  
    626 def process_parameters(model_info): 
    627     """ 
    628     Process parameter block, precalculating parameter details. 
    629     """ 
    630     # convert parameters into named tuples 
    631     for p in model_info['parameters']: 
    632         if p[4] == '' and (p[0].startswith('sld') or p[0].endswith('sld')): 
    633             p[4] = 'sld' 
    634             # TODO: make sure all models explicitly label their sld parameters 
    635             #raise ValueError("%s.%s needs to be explicitly set to type 'sld'" %(model_info['id'], p[0])) 
    636  
    637     pars = [Parameter(*p) for p in model_info['parameters']] 
    638     # Fill in the derived attributes 
    639     model_info['parameters'] = pars 
    640     partype = categorize_parameters(pars) 
    641     model_info['limits'] = dict((p.name, p.limits) for p in pars) 
    642     model_info['partype'] = partype 
    643     model_info['defaults'] = dict((p.name, p.default) for p in pars) 
    644     if model_info.get('demo', None) is None: 
    645         model_info['demo'] = model_info['defaults'] 
    646     model_info['has_2d'] = partype['orientation'] or partype['magnetic'] 
     478        pars = [qx, qy] + partable.iqxy_parameters 
     479        source.append(_gen_fn('Iqxy', pars, model_info['Iqxy'])) 
     480 
     481    # Define the parameter table 
     482    source.append("#define PARAMETER_TABLE \\") 
     483    source.append("\\\n".join(p.as_definition() 
     484                              for p in partable.kernel_parameters)) 
     485 
     486    # Define the function calls 
     487    if partable.form_volume_parameters: 
     488        refs = _call_pars("_v.", partable.form_volume_parameters) 
     489        call_volume = "#define CALL_VOLUME(_v) form_volume(%s)" % (",".join(refs)) 
     490    else: 
     491        # Model doesn't have volume.  We could make the kernel run a little 
     492        # faster by not using/transferring the volume normalizations, but 
     493        # the ifdef's reduce readability more than is worthwhile. 
     494        call_volume = "#define CALL_VOLUME(v) 1.0" 
     495    source.append(call_volume) 
     496 
     497    refs = ["_q[_i]"] + _call_pars("_v.", partable.iq_parameters) 
     498    call_iq = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(refs)) 
     499    if _have_Iqxy(user_code): 
     500        # Call 2D model 
     501        refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("_v.", partable.iqxy_parameters) 
     502        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iqxy(%s)" % (",".join(refs)) 
     503    else: 
     504        # Call 1D model with sqrt(qx^2 + qy^2) 
     505        warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 
     506        # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 
     507        pars_sqrt = ["sqrt(_q[2*_i]*_q[2*_i]+_q[2*_i+1]*_q[2*_i+1])"] + refs[1:] 
     508        call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt)) 
     509 
     510    # Fill in definitions for numbers of parameters 
     511    source.append("#define MAX_PD %s"%partable.max_pd) 
     512    source.append("#define NPARS %d"%partable.npars) 
     513 
     514    # TODO: allow mixed python/opencl kernels? 
     515 
     516    # define the Iq kernel 
     517    source.append("#define KERNEL_NAME %s_Iq"%model_info['name']) 
     518    source.append(call_iq) 
     519    source.append(kernel_code) 
     520    source.append("#undef CALL_IQ") 
     521    source.append("#undef KERNEL_NAME") 
     522 
     523    # define the Iqxy kernel from the same source with different #defines 
     524    source.append("#define KERNEL_NAME %s_Iqxy"%model_info['name']) 
     525    source.append(call_iqxy) 
     526    source.append(kernel_code) 
     527    source.append("#undef CALL_IQ") 
     528    source.append("#undef KERNEL_NAME") 
     529 
     530    return '\n'.join(source) 
     531 
     532class CoordinationDetails(object): 
     533    def __init__(self, model_info): 
     534        parameters = model_info['parameters'] 
     535        max_pd = parameters.max_pd 
     536        npars = parameters.npars 
     537        par_offset = 4*max_pd 
     538        self.details = np.zeros(par_offset + 3*npars + 4, 'i4') 
     539 
     540        # generate views on different parts of the array 
     541        self._pd_par     = self.details[0*max_pd:1*max_pd] 
     542        self._pd_length  = self.details[1*max_pd:2*max_pd] 
     543        self._pd_offset  = self.details[2*max_pd:3*max_pd] 
     544        self._pd_stride  = self.details[3*max_pd:4*max_pd] 
     545        self._par_offset = self.details[par_offset+0*npars:par_offset+1*npars] 
     546        self._par_coord  = self.details[par_offset+1*npars:par_offset+2*npars] 
     547        self._pd_coord   = self.details[par_offset+2*npars:par_offset+3*npars] 
     548 
     549        # theta_par is fixed 
     550        self.details[-1] = parameters.theta_offset 
     551 
     552    @property 
     553    def ctypes(self): return self.details.ctypes 
     554    @property 
     555    def pd_par(self): return self._pd_par 
     556    @property 
     557    def pd_length(self): return self._pd_length 
     558    @property 
     559    def pd_offset(self): return self._pd_offset 
     560    @property 
     561    def pd_stride(self): return self._pd_stride 
     562    @property 
     563    def pd_coord(self): return self._pd_coord 
     564    @property 
     565    def par_coord(self): return self._par_coord 
     566    @property 
     567    def par_offset(self): return self._par_offset 
     568    @property 
     569    def num_coord(self): return self.details[-2] 
     570    @num_coord.setter 
     571    def num_coord(self, v): self.details[-2] = v 
     572    @property 
     573    def total_pd(self): return self.details[-3] 
     574    @total_pd.setter 
     575    def total_pd(self, v): self.details[-3] = v 
     576    @property 
     577    def num_active(self): return self.details[-4] 
     578    @num_active.setter 
     579    def num_active(self, v): self.details[-4] = v 
     580 
     581    def show(self): 
     582        print("total_pd", self.total_pd) 
     583        print("num_active", self.num_active) 
     584        print("pd_par", self.pd_par) 
     585        print("pd_length", self.pd_length) 
     586        print("pd_offset", self.pd_offset) 
     587        print("pd_stride", self.pd_stride) 
     588        print("par_offsets", self.par_offset) 
     589        print("num_coord", self.num_coord) 
     590        print("par_coord", self.par_coord) 
     591        print("pd_coord", self.pd_coord) 
     592        print("theta par", self.details[-1]) 
     593 
     594def mono_details(model_info): 
     595    details = CoordinationDetails(model_info) 
     596    # The zero defaults for monodisperse systems are mostly fine 
     597    details.par_offset[:] = np.arange(2, len(details.par_offset)+2) 
     598    return details 
     599 
     600def poly_details(model_info, weights): 
     601    #print("weights",weights) 
     602    weights = weights[2:] # Skip scale and background 
     603 
     604    # Decreasing list of polydispersity lengths 
     605    # Note: the reversing view, x[::-1], does not require a copy 
     606    pd_length = np.array([len(w) for w in weights]) 
     607    num_active = np.sum(pd_length>1) 
     608    if num_active > model_info['parameters'].max_pd: 
     609        raise ValueError("Too many polydisperse parameters") 
     610 
     611    pd_offset = np.cumsum(np.hstack((0, pd_length))) 
     612    idx = np.argsort(pd_length)[::-1][:num_active] 
     613    par_length = np.array([max(len(w),1) for w in weights]) 
     614    pd_stride = np.cumprod(np.hstack((1, par_length[idx]))) 
     615    par_offsets = np.cumsum(np.hstack((2, par_length))) 
     616 
     617    details = CoordinationDetails(model_info) 
     618    details.pd_par[:num_active] = idx 
     619    details.pd_length[:num_active] = pd_length[idx] 
     620    details.pd_offset[:num_active] = pd_offset[idx] 
     621    details.pd_stride[:num_active] = pd_stride[:-1] 
     622    details.par_offset[:] = par_offsets[:-1] 
     623    details.total_pd = pd_stride[-1] 
     624    details.num_active = num_active 
     625    # Without constraints coordinated parameters are just the pd parameters 
     626    details.par_coord[:num_active] = idx 
     627    details.pd_coord[:num_active] = 2**np.arange(num_active) 
     628    details.num_coord = num_active 
     629    #details.show() 
     630    return details 
     631 
     632def constrained_poly_details(model_info, weights, constraints): 
     633    # Need to find the independently varying pars and sort them 
     634    # Need to build a coordination list for the dependent variables 
     635    # Need to generate a constraints function which takes values 
     636    # and weights, returning par blocks 
     637    raise NotImplementedError("Can't handle constraints yet") 
     638 
     639 
     640def create_default_functions(model_info): 
     641    """ 
     642    Autogenerate missing functions, such as Iqxy from Iq. 
     643 
     644    This only works for Iqxy when Iq is written in python. :func:`make_source` 
     645    performs a similar role for Iq written in C. 
     646    """ 
     647    if callable(model_info['Iq']) and model_info['Iqxy'] is None: 
     648        partable = model_info['parameters'] 
     649        if partable.has_2d: 
     650            raise ValueError("Iqxy model is missing") 
     651        Iq = model_info['Iq'] 
     652        def Iqxy(qx, qy, **kw): 
     653            return Iq(np.sqrt(qx**2 + qy**2), **kw) 
     654        model_info['Iqxy'] = Iqxy 
    647655 
    648656 
    … …  
    682690      model variants (e.g., the list of cases) or is None if there are no 
    683691      model variants 
    684     * *defaults* is the *{parameter: value}* table built from the parameter 
    685       description table. 
    686     * *limits* is the *{parameter: [min, max]}* table built from the 
    687       parameter description table. 
    688     * *partypes* categorizes the model parameters. See 
     692    * *par_type* categorizes the model parameters. See 
    689693      :func:`categorize_parameters` for details. 
    690694    * *demo* contains the *{parameter: value}* map used in compare (and maybe 
    … …  
    700704      *model_info* blocks for the composition objects.  This allows us to 
    701705      build complete product and mixture models from just the info. 
    702  
    703706    """ 
    704707    # TODO: maybe turn model_info into a class ModelDefinition 
    705     parameters = COMMON_PARAMETERS + kernel_module.parameters 
     708    #print("make parameter table", kernel_module.parameters) 
     709    parameters = make_parameter_table(kernel_module.parameters) 
    706710    filename = abspath(kernel_module.__file__) 
    707711    kernel_id = splitext(basename(filename))[0] 
    … …  
    713717        filename=abspath(kernel_module.__file__), 
    714718        name=name, 
    715         title=kernel_module.title, 
    716         description=kernel_module.description, 
     719        title=getattr(kernel_module, 'title', name+" model"), 
     720        description=getattr(kernel_module, 'description', 'no description'), 
    717721        parameters=parameters, 
    718722        composition=None, 
    … …  
    721725        single=getattr(kernel_module, 'single', True), 
    722726        structure_factor=getattr(kernel_module, 'structure_factor', False), 
     727        profile_axes=getattr(kernel_module, 'profile_axes', ['x','y']), 
    723728        variant_info=getattr(kernel_module, 'invariant_info', None), 
    724729        demo=getattr(kernel_module, 'demo', None), 
    … …  
    726731        tests=getattr(kernel_module, 'tests', []), 
    727732        ) 
    728     process_parameters(model_info) 
     733    set_demo(model_info, getattr(kernel_module, 'demo', None)) 
    729734    # Check for optional functions 
    730     functions = "ER VR form_volume Iq Iqxy shape sesans".split() 
     735    functions = "ER VR form_volume Iq Iqxy profile sesans".split() 
    731736    model_info.update((k, getattr(kernel_module, k, None)) for k in functions) 
     737    create_default_functions(model_info) 
     738    # Precalculate the monodisperse parameters 
     739    # TODO: make this a lazy evaluator 
     740    # make_model_info is called for every model on sasview startup 
     741    model_info['mono_details'] = mono_details(model_info) 
    732742    return model_info 
    733743 
    … …  
    796806    Program which prints the source produced by the model. 
    797807    """ 
     808    import sys 
    798809    if len(sys.argv) <= 1: 
    799810        print("usage: python -m sasmodels.generate modelname") 

  • sasmodels/kernelcl.py

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

  • sasmodels/kerneldll.py

    r4d76711 r1e2a1ba  
    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): 
    … …  
    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

    r4d76711 r1e2a1ba  
    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/list_pars.py

    ra84a0ca r21b116f  
    2525    for name in sorted(MODELS): 
    2626        model_info = load_model_info(name) 
    27         for p in model_info['parameters']: 
     27        for p in model_info['parameters'].kernel_parameters: 
    2828            partable.setdefault(p.name, []) 
    2929            partable[p.name].append(name) 

  • sasmodels/mixture.py

    r72a081d 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

    r4d76711 rd2bb604  
    8989 
    9090        if is_py:  # kernel implemented in python 
     91            continue # TODO: re-enable python tests 
    9192            test_name = "Model: %s, Kernel: python"%model_name 
    9293            test_method_name = "test_%s_python" % model_name 

  • sasmodels/models/cylinder.c

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

  • sasmodels/models/flexible_cylinder.c

    r43b7eea rd2bb604  
    1 double form_volume(double length, double kuhn_length, double radius); 
    2 double Iq(double q, double length, double kuhn_length, double radius, 
    3           double sld, double solvent_sld); 
    4 double Iqxy(double qx, double qy, double length, double kuhn_length, 
    5             double radius, double sld, double solvent_sld); 
    6 double flexible_cylinder_kernel(double q, double length, double kuhn_length, 
    7                                 double radius, double sld, double solvent_sld); 
    8  
    9  
    10 double form_volume(double length, double kuhn_length, double radius) 
     1static double 
     2form_volume(length, kuhn_length, radius) 
    113{ 
    12  
    13       return 0.0; 
     4    return 1.; 
    145} 
    156 
    16 double flexible_cylinder_kernel(double q, 
    17           double length, 
    18           double kuhn_length, 
    19           double radius, 
    20           double sld, 
    21           double solvent_sld) 
     7static double 
     8Iq(double q, 
     9   double length, 
     10   double kuhn_length, 
     11   double radius, 
     12   double sld, 
     13   double solvent_sld) 
    2214{ 
    23  
    24     const double cont = sld-solvent_sld; 
    25     const double qr = q*radius; 
    26     //const double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr); 
    27     const double crossSect = sas_J1c(qr); 
    28     double flex = Sk_WR(q,length,kuhn_length); 
    29     flex *= crossSect*crossSect; 
    30     flex *= M_PI*radius*radius*length; 
    31     flex *= cont*cont; 
    32     flex *= 1.0e-4; 
    33     return flex; 
     15    const double contrast = sld - solvent_sld; 
     16    const double cross_section = sas_J1c(q*radius); 
     17    const double volume = M_PI*radius*radius*length; 
     18    const double flex = Sk_WR(q, length, kuhn_length); 
     19    return 1.0e-4 * volume * square(contrast*cross_section) * flex; 
    3420} 
    35  
    36 double Iq(double q, 
    37           double length, 
    38           double kuhn_length, 
    39           double radius, 
    40           double sld, 
    41           double solvent_sld) 
    42 { 
    43  
    44     double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld); 
    45     return result; 
    46 } 
    47  
    48 double Iqxy(double qx, double qy, 
    49             double length, 
    50             double kuhn_length, 
    51             double radius, 
    52             double sld, 
    53             double solvent_sld) 
    54 { 
    55     double q; 
    56     q = sqrt(qx*qx+qy*qy); 
    57     double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld); 
    58  
    59     return result; 
    60 } 

  • sasmodels/models/gel_fit.c

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

  • sasmodels/models/hardsphere.py

    rec45c4f rd2bb604  
    149149   """ 
    150150 
    151 Iqxy = """ 
    152     // never called since no orientation or magnetic parameters. 
    153     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    154     """ 
    155  
    156151# ER defaults to 0.0 
    157152# VR defaults to 1.0 

  • sasmodels/models/hayter_msa.py

    rec45c4f rd2bb604  
    8787    return 1.0; 
    8888    """ 
    89 Iqxy = """ 
    90     // never called since no orientation or magnetic parameters. 
    91     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    92     """ 
    9389# ER defaults to 0.0 
    9490# VR defaults to 1.0 

  • sasmodels/models/lamellar.py

    rec45c4f rd2bb604  
    8282    """ 
    8383 
    84 Iqxy = """ 
    85     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    86     """ 
    87  
    8884# ER defaults to 0.0 
    8985# VR defaults to 1.0 

  • sasmodels/models/lamellar_hg.py

    rec45c4f rd2bb604  
    101101    """ 
    102102 
    103 Iqxy = """ 
    104     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    105     """ 
    106  
    107103# ER defaults to 0.0 
    108104# VR defaults to 1.0 

  • sasmodels/models/lamellar_hg_stack_caille.py

    rec45c4f rd2bb604  
    120120    """ 
    121121 
    122 Iqxy = """ 
    123     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    124     """ 
    125  
    126122# ER defaults to 0.0 
    127123# VR defaults to 1.0 

  • sasmodels/models/lamellar_stack_caille.py

    rec45c4f rd2bb604  
    104104    """ 
    105105 
    106 Iqxy = """ 
    107     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    108     """ 
    109  
    110106# ER defaults to 0.0 
    111107# VR defaults to 1.0 

  • sasmodels/models/lamellar_stack_paracrystal.py

    rec45c4f rd2bb604  
    132132    """ 
    133133 
    134 Iqxy = """ 
    135     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    136     """ 
    137  
    138134# ER defaults to 0.0 
    139135# VR defaults to 1.0 

  • sasmodels/models/lib/Si.c

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

  • sasmodels/models/lib/polevl.c

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

  • sasmodels/models/lib/sas_J0.c

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

  • sasmodels/models/lib/sas_J1.c

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

  • sasmodels/models/lib/sas_JN.c

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

  • sasmodels/models/lib/sph_j1c.c

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

  • sasmodels/models/lib/sphere_form.c

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

  • sasmodels/models/lib/wrc_cyl.c

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

  • sasmodels/models/onion.c

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

  • sasmodels/models/onion.py

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

  • sasmodels/models/rpa.c

    r13ed84c rd2bb604  
    11double Iq(double q, double case_num, 
    2     double Na, double Phia, double va, double a_sld, double ba, 
    3     double Nb, double Phib, double vb, double b_sld, double bb, 
    4     double Nc, double Phic, double vc, double c_sld, double bc, 
    5     double Nd, double Phid, double vd, double d_sld, double bd, 
     2    double N[], double Phi[], double v[], double L[], double b[], 
    63    double Kab, double Kac, double Kad, 
    74    double Kbc, double Kbd, double Kcd 
    85    ); 
    96 
    10 double Iqxy(double qx, double qy, double case_num, 
    11     double Na, double Phia, double va, double a_sld, double ba, 
    12     double Nb, double Phib, double vb, double b_sld, double bb, 
    13     double Nc, double Phic, double vc, double c_sld, double bc, 
    14     double Nd, double Phid, double vd, double d_sld, double bd, 
    15     double Kab, double Kac, double Kad, 
    16     double Kbc, double Kbd, double Kcd 
    17     ); 
    18  
    19 double form_volume(void); 
    20  
    21 double form_volume(void) 
    22 { 
    23     return 1.0; 
    24 } 
    25  
    267double Iq(double q, double case_num, 
    27     double Na, double Phia, double va, double La, double ba, 
    28     double Nb, double Phib, double vb, double Lb, double bb, 
    29     double Nc, double Phic, double vc, double Lc, double bc, 
    30     double Nd, double Phid, double vd, double Ld, double bd, 
     8    double N[], double Phi[], double v[], double L[], double b[], 
    319    double Kab, double Kac, double Kad, 
    3210    double Kbc, double Kbd, double Kcd 
    … …  
    3614#if 0  // Sasview defaults 
    3715  if (icase <= 1) { 
    38     Na=Nb=1000.0; 
    39     Phia=Phib=0.0000001; 
     16    N[0]=N[1]=1000.0; 
     17    Phi[0]=Phi[1]=0.0000001; 
    4018    Kab=Kac=Kad=Kbc=Kbd=-0.0004; 
    41     La=Lb=1e-12; 
    42     va=vb=100.0; 
    43     ba=bb=5.0; 
     19    L[0]=L[1]=1e-12; 
     20    v[0]=v[1]=100.0; 
     21    b[0]=b[1]=5.0; 
    4422  } else if (icase <= 4) { 
    45     Phia=0.0000001; 
     23    Phi[0]=0.0000001; 
    4624    Kab=Kac=Kad=-0.0004; 
    47     La=1e-12; 
    48     va=100.0; 
    49     ba=5.0; 
     25    L[0]=1e-12; 
     26    v[0]=100.0; 
     27    b[0]=5.0; 
    5028  } 
    5129#else 
    5230  if (icase <= 1) { 
    53     Na=Nb=0.0; 
    54     Phia=Phib=0.0; 
     31    N[0]=N[1]=0.0; 
     32    Phi[0]=Phi[1]=0.0; 
    5533    Kab=Kac=Kad=Kbc=Kbd=0.0; 
    56     La=Lb=Ld; 
    57     va=vb=vd; 
    58     ba=bb=0.0; 
     34    L[0]=L[1]=L[3]; 
     35    v[0]=v[1]=v[3]; 
     36    b[0]=b[1]=0.0; 
    5937  } else if (icase <= 4) { 
    60     Na = 0.0; 
    61     Phia=0.0; 
     38    N[0] = 0.0; 
     39    Phi[0]=0.0; 
    6240    Kab=Kac=Kad=0.0; 
    63     La=Ld; 
    64     va=vd; 
    65     ba=0.0; 
     41    L[0]=L[3]; 
     42    v[0]=v[3]; 
     43    b[0]=0.0; 
    6644  } 
    6745#endif 
    6846 
    69   const double Xa = q*q*ba*ba*Na/6.0; 
    70   const double Xb = q*q*bb*bb*Nb/6.0; 
    71   const double Xc = q*q*bc*bc*Nc/6.0; 
    72   const double Xd = q*q*bd*bd*Nd/6.0; 
     47  const double Xa = q*q*b[0]*b[0]*N[0]/6.0; 
     48  const double Xb = q*q*b[1]*b[1]*N[1]/6.0; 
     49  const double Xc = q*q*b[2]*b[2]*N[2]/6.0; 
     50  const double Xd = q*q*b[3]*b[3]*N[3]/6.0; 
    7351 
    7452  // limit as Xa goes to 0 is 1 
    … …  
    9876#if 0 
    9977  const double S0aa = icase<5 
    100                       ? 1.0 : Na*Phia*va*Paa; 
     78                      ? 1.0 : N[0]*Phi[0]*v[0]*Paa; 
    10179  const double S0bb = icase<2 
    102                       ? 1.0 : Nb*Phib*vb*Pbb; 
    103   const double S0cc = Nc*Phic*vc*Pcc; 
    104   const double S0dd = Nd*Phid*vd*Pdd; 
     80                      ? 1.0 : N[1]*Phi[1]*v[1]*Pbb; 
     81  const double S0cc = N[2]*Phi[2]*v[2]*Pcc; 
     82  const double S0dd = N[3]*Phi[3]*v[3]*Pdd; 
    10583  const double S0ab = icase<8 
    106                       ? 0.0 : sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb; 
     84                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 
    10785  const double S0ac = icase<9 
    108                       ? 0.0 : sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb); 
     86                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 
    10987  const double S0ad = icase<9 
    110                       ? 0.0 : sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc); 
     88                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 
    11189  const double S0bc = (icase!=4 && icase!=7 && icase!= 9) 
    112                       ? 0.0 : sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc; 
     90                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 
    11391  const double S0bd = (icase!=4 && icase!=7 && icase!= 9) 
    114                       ? 0.0 : sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc); 
     92                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 
    11593  const double S0cd = (icase==0 || icase==2 || icase==5) 
    116                       ? 0.0 : sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd; 
     94                      ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 
    11795#else  // sasview equivalent 
    118 //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,Nc,Phic,vc,Pcc); 
    119   double S0aa = Na*Phia*va*Paa; 
    120   double S0bb = Nb*Phib*vb*Pbb; 
    121   double S0cc = Nc*Phic*vc*Pcc; 
    122   double S0dd = Nd*Phid*vd*Pdd; 
    123   double S0ab = sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb; 
    124   double S0ac = sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb); 
    125   double S0ad = sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc); 
    126   double S0bc = sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc; 
    127   double S0bd = sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc); 
    128   double S0cd = sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd; 
     96//printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc); 
     97  double S0aa = N[0]*Phi[0]*v[0]*Paa; 
     98  double S0bb = N[1]*Phi[1]*v[1]*Pbb; 
     99  double S0cc = N[2]*Phi[2]*v[2]*Pcc; 
     100  double S0dd = N[3]*Phi[3]*v[3]*Pdd; 
     101  double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 
     102  double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 
     103  double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 
     104  double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 
     105  double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 
     106  double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 
    129107switch(icase){ 
    130108  case 0: 
    … …  
    311289  // Note: 1e-13 to convert from fm to cm for scattering length 
    312290  const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13; 
    313   const double Lad = icase<5 ? 0.0 : (La/va - Ld/vd)*sqrt_Nav; 
    314   const double Lbd = icase<2 ? 0.0 : (Lb/vb - Ld/vd)*sqrt_Nav; 
    315   const double Lcd = (Lc/vc - Ld/vd)*sqrt_Nav; 
     291  const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav; 
     292  const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav; 
     293  const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav; 
    316294 
    317295  const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33 
    … …  
    321299 
    322300} 
    323  
    324 double Iqxy(double qx, double qy, 
    325     double case_num, 
    326     double Na, double Phia, double va, double a_sld, double ba, 
    327     double Nb, double Phib, double vb, double b_sld, double bb, 
    328     double Nc, double Phic, double vc, double c_sld, double bc, 
    329     double Nd, double Phid, double vd, double d_sld, double bd, 
    330     double Kab, double Kac, double Kad, 
    331     double Kbc, double Kbd, double Kcd 
    332     ) 
    333 { 
    334     double q = sqrt(qx*qx + qy*qy); 
    335     return Iq(q, 
    336         case_num, 
    337         Na, Phia, va, a_sld, ba, 
    338         Nb, Phib, vb, b_sld, bb, 
    339         Nc, Phic, vc, c_sld, bc, 
    340         Nd, Phid, vd, d_sld, bd, 
    341         Kab, Kac, Kad, 
    342         Kbc, Kbd, Kcd); 
    343 } 

  • sasmodels/models/rpa.py

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

  • sasmodels/models/spherical_sld.py

    rec45c4f rd2bb604  
    170170# pylint: disable=bad-whitespace, line-too-long 
    171171#            ["name", "units", default, [lower, upper], "type", "description"], 
    172 parameters = [["n_shells",         "",               1,      [0, 9],         "", "number of shells"], 
     172parameters = [["n",                "",               1,      [0, 9],         "", "number of shells"], 
    173173              ["radius_core",      "Ang",            50.0,   [0, inf],       "", "intern layer thickness"], 
    174174              ["sld_core",         "1e-6/Ang^2",     2.07,   [-inf, inf],    "", "sld function flat"], 
    … …  
    192192 
    193193demo = dict( 
    194     n_shells=4, 
     194    n=4, 
    195195    scale=1.0, 
    196196    solvent_sld=1.0, 

  • sasmodels/models/squarewell.py

    rec45c4f rd2bb604  
    123123""" 
    124124 
    125 Iqxy = """ 
    126     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    127     """ 
    128  
    129125# ER defaults to 0.0 
    130126# VR defaults to 1.0 

  • sasmodels/models/stickyhardsphere.py

    rec45c4f rd2bb604  
    171171""" 
    172172 
    173 Iqxy = """ 
    174     return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    175     """ 
    176  
    177173# ER defaults to 0.0 
    178174# VR defaults to 1.0 

  • sasmodels/product.py

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

  • sasmodels/resolution.py

    r4d76711 rd2bb604  
    502502    from scipy.integrate import romberg 
    503503 
    504     if any(k not in form.info['defaults'] for k in pars.keys()): 
    505         keys = set(form.info['defaults'].keys()) 
    506         extra = set(pars.keys()) - keys 
    507         raise ValueError("bad parameters: [%s] not in [%s]"% 
    508                          (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 
     504    par_set = set([p.name for p in form.info['parameters'].call_parameters]) 
     505    if any(k not in par_set for k in pars.keys()): 
     506        extra = set(pars.keys()) - par_set 
     507        raise ValueError("bad parameters: [%s] not in [%s]" 
     508                         % (", ".join(sorted(extra)), 
     509                            ", ".join(sorted(pars.keys())))) 
    509510 
    510511    if np.isscalar(width): 
    … …  
    556557    from scipy.integrate import romberg 
    557558 
    558     if any(k not in form.info['defaults'] for k in pars.keys()): 
    559         keys = set(form.info['defaults'].keys()) 
    560         extra = set(pars.keys()) - keys 
    561         raise ValueError("bad parameters: [%s] not in [%s]"% 
    562                          (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 
     559    par_set = set([p.name for p in form.info['parameters'].call_parameters]) 
     560    if any(k not in par_set for k in pars.keys()): 
     561        extra = set(pars.keys()) - par_set 
     562        raise ValueError("bad parameters: [%s] not in [%s]" 
     563                         % (", ".join(sorted(extra)), 
     564                            ", ".join(sorted(pars.keys())))) 
    563565 
    564566    _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq) 

  • sasmodels/sasview_model.py

    r4d76711 rd2bb604  
    2626from . import custom 
    2727from . import generate 
     28from . import weights 
    2829 
    2930def load_standard_models(): 
    … …  
    8485        self._model = None 
    8586        model_info = self._model_info 
     87        parameters = model_info['parameters'] 
    8688 
    8789        self.name = model_info['name'] 
    8890        self.description = model_info['description'] 
    8991        self.category = None 
    90         self.multiplicity_info = None 
    91         self.is_multifunc = False 
     92        #self.is_multifunc = False 
     93        for p in parameters.kernel_parameters: 
     94            if p.is_control: 
     95                profile_axes = model_info['profile_axes'] 
     96                self.multiplicity_info = [ 
     97                    p.limits[1], p.name, p.choices, profile_axes[0] 
     98                    ] 
     99                break 
     100        else: 
     101            self.multiplicity_info = [] 
    92102 
    93103        ## interpret the parameters 
    … …  
    96106        self.params = collections.OrderedDict() 
    97107        self.dispersion = dict() 
    98         partype = model_info['partype'] 
    99  
    100         for p in model_info['parameters']: 
     108 
     109        self.orientation_params = [] 
     110        self.magnetic_params = [] 
     111        self.fixed = [] 
     112        for p in parameters.user_parameters(): 
    101113            self.params[p.name] = p.default 
    102114            self.details[p.name] = [p.units] + p.limits 
    103  
    104         for name in partype['pd-2d']: 
    105             self.dispersion[name] = { 
    106                 'width': 0, 
    107                 'npts': 35, 
    108                 'nsigmas': 3, 
    109                 'type': 'gaussian', 
    110             } 
    111  
    112         self.orientation_params = ( 
    113             partype['orientation'] 
    114             + [n + '.width' for n in partype['orientation']] 
    115             + partype['magnetic']) 
    116         self.magnetic_params = partype['magnetic'] 
    117         self.fixed = [n + '.width' for n in partype['pd-2d']] 
     115            if p.polydisperse: 
     116                self.dispersion[p.name] = { 
     117                    'width': 0, 
     118                    'npts': 35, 
     119                    'nsigmas': 3, 
     120                    'type': 'gaussian', 
     121                } 
     122            if p.type == 'orientation': 
     123                self.orientation_params.append(p.name) 
     124                self.orientation_params.append(p.name+".width") 
     125                self.fixed.append(p.name+".width") 
     126            if p.type == 'magnetic': 
     127                self.orientation_params.append(p.name) 
     128                self.magnetic_params.append(p.name) 
     129                self.fixed.append(p.name+".width") 
     130 
    118131        self.non_fittable = [] 
    119132 
    … …  
    234247        """ 
    235248        # TODO: fix test so that parameter order doesn't matter 
    236         ret = ['%s.%s' % (d.lower(), p) 
    237                for d in self._model_info['partype']['pd-2d'] 
    238                for p in ('npts', 'nsigmas', 'width')] 
     249        ret = ['%s.%s' % (p.name.lower(), ext) 
     250               for p in self._model_info['parameters'].user_parameters() 
     251               for ext in ('npts', 'nsigmas', 'width') 
     252               if p.polydisperse] 
    239253        #print(ret) 
    240254        return ret 
    … …  
    309323            # Check whether we have a list of ndarrays [qx,qy] 
    310324            qx, qy = qdist 
    311             partype = self._model_info['partype'] 
    312             if not partype['orientation'] and not partype['magnetic']: 
     325            if not self._model_info['parameters'].has_2d: 
    313326                return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2)) 
    314327            else: 
    … …  
    335348            self._model = core.build_model(self._model_info) 
    336349        q_vectors = [np.asarray(q) for q in args] 
    337         fn = self._model.make_kernel(q_vectors) 
    338         pars = [self.params[v] for v in fn.fixed_pars] 
    339         pd_pars = [self._get_weights(p) for p in fn.pd_pars] 
    340         result = fn(pars, pd_pars, self.cutoff) 
    341         fn.q_input.release() 
    342         fn.release() 
     350        kernel = self._model.make_kernel(q_vectors) 
     351        pairs = [self._get_weights(p) 
     352                 for p in self._model_info['parameters'].call_parameters] 
     353        details, weights, values = core.build_details(kernel, pairs) 
     354        result = kernel(details, weights, values, cutoff=self.cutoff) 
     355        kernel.q_input.release() 
     356        kernel.release() 
    343357        return result 
    344358 
    … …  
    415429        Return dispersion weights for parameter 
    416430        """ 
    417         from . import weights 
    418         relative = self._model_info['partype']['pd-rel'] 
    419         limits = self._model_info['limits'] 
    420         dis = self.dispersion[par] 
    421         value, weight = weights.get_weights( 
    422             dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 
    423             self.params[par], limits[par], par in relative) 
    424         return value, weight / np.sum(weight) 
    425  
     431        if par.polydisperse: 
     432            dis = self.dispersion[par.name] 
     433            value, weight = weights.get_weights( 
     434                dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 
     435                self.params[par.name], par.limits, par.relative_pd) 
     436            return value, weight / np.sum(weight) 
     437        else: 
     438            return [self.params[par.name]], [] 
    426439 
    427440def test_model(): 
Note: See TracChangeset for help on using the changeset viewer.