Changeset ce27e21 in sasmodels


Ignore:
Timestamp:
Aug 24, 2014 5:18:14 PM (10 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:
1780d59
Parents:
14de349
Message:

first pass for sasview wrapper around opencl models

Files:
4 added
8 edited

Legend:

Unmodified
Added
Removed

  • compare-new.py

    r14de349 rce27e21  
    4848            gpu = model.theory() 
    4949        gpu_time = toc()*1000./Ngpu 
    50         print "ocl t=%.1f ms"%gpu_time 
     50        print "ocl t=%.1f ms, intensity=%.0f"%(gpu_time, sum(gpu[~np.isnan(gpu)])) 
    5151        #print max(gpu), min(gpu) 
    5252 
    … …  
    6161        print "dll t=%.1f ms"%cpu_time 
    6262 
    63     elif 1: # Hack to check new vs old for GpuCylinder 
     63    elif 0: # Hack to check new vs old for GpuCylinder 
    6464        from Models.code_cylinder_f import GpuCylinder as oldgpu 
    6565        from sasmodel import SasModel 
    … …  
    7878            cpu = cpumodel.evalDistribution([data.qx_data, data.qy_data]) 
    7979        cpu_time = toc()*1000./Ncpu 
    80         print "sasview t=%.1f ms"%cpu_time 
     80        print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[model.index])) 
    8181 
    8282    if Ngpu > 0 and Ncpu > 0: 

  • sasmodel.py

    ra953943 rce27e21  
    44import datetime 
    55import numpy as np 
    6 import pyopencl as cl 
    7 from bumps.names import Parameter 
    8 from sans.dataloader.loader import Loader 
    9 from sans.dataloader.manipulations import Ringcut, Boxcut 
    106 
    117 
    … …  
    2016 
    2117def load_data(filename): 
     18    from sans.dataloader.loader import Loader 
    2219    loader = Loader() 
    2320    data = loader.load(filename) 
    … …  
    5754 
    5855def set_beam_stop(data, radius, outer=None): 
     56    from sans.dataloader.manipulations import Ringcut 
    5957    if hasattr(data, 'qx_data'): 
    6058        data.mask = Ringcut(0, radius)(data) 
    … …  
    6765 
    6866def set_half(data, half): 
     67    from sans.dataloader.manipulations import Boxcut 
    6968    if half == 'right': 
    7069        data.mask += Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data) 
    … …  
    7372 
    7473def set_top(data, max): 
     74    from sans.dataloader.manipulations import Boxcut 
    7575    data.mask += Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=max)(data) 
    7676 
    … …  
    144144    global GPU_CONTEXT, GPU_QUEUE 
    145145    if GPU_CONTEXT is None: 
     146        import pyopencl as cl 
    146147        GPU_CONTEXT = cl.create_some_context() 
    147148        GPU_QUEUE = cl.CommandQueue(GPU_CONTEXT) 
    … …  
    151152class SasModel(object): 
    152153    def __init__(self, data, model, dtype='float32', **kw): 
     154        from bumps.names import Parameter 
     155 
    153156        self.__dict__['_parameters'] = {} 
    154157        #self.name = data.filename 

  • sasmodels/core.py

    r14de349 rce27e21  
    44import sys, os 
    55import datetime 
     6import warnings 
    67 
    78import numpy as np 
    … …  
    1415    return gen.make(modelpath) 
    1516 
     17 
     18 
    1619def opencl_model(modelname, dtype="single"): 
    1720    from sasmodels import gpu 
    1821 
    19     source, meta, _ = load_model(modelname) 
     22    source, info, _ = load_model(modelname) 
    2023    # for debugging, save source to a .cl file, edit it, and reload as model 
    2124    #open(modelname+'.cl','w').write(source) 
    2225    #source = open(modelname+'.cl','r').read() 
    23     return gpu.GpuModel(source, meta, dtype) 
     26    return gpu.GpuModel(source, info, dtype) 
    2427 
    2528 
    … …  
    3134    COMPILE = "cc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm" 
    3235DLL_PATH = "/tmp" 
    33 def dll_path(meta): 
     36 
     37 
     38def dll_path(info): 
    3439    from os.path import join as joinpath, split as splitpath, splitext 
    35     basename = splitext(splitpath(meta['filename'])[1])[0] 
     40    basename = splitext(splitpath(info['filename'])[1])[0] 
    3641    return joinpath(DLL_PATH, basename+'.so') 
     42 
    3743 
    3844def dll_model(modelname): 
    … …  
    4046    from sasmodels import dll 
    4147 
    42     source, meta, _ = load_model(modelname) 
    43     dllpath = dll_path(meta) 
     48    source, info, _ = load_model(modelname) 
     49    dllpath = dll_path(info) 
    4450    if not os.path.exists(dllpath) \ 
    45             or (os.path.getmtime(dllpath) < os.path.getmtime(meta['filename'])): 
     51            or (os.path.getmtime(dllpath) < os.path.getmtime(info['filename'])): 
    4652        # Replace with a proper temp file 
    4753        srcfile = '/tmp/%s.c'%modelname 
    4854        open(srcfile, 'w').write(source) 
    4955        os.system(COMPILE%(srcfile, dllpath)) 
    50     return dll.DllModel(dllpath, meta) 
     56    return dll.DllModel(dllpath, info) 
     57 
    5158 
    5259TIC = None 
    … …  
    5764    return TIC 
    5865 
     66 
    5967def toc(): 
    6068    return TIC() 
     69 
    6170 
    6271def load_data(filename): 
    … …  
    6877    return data 
    6978 
     79 
    7080def fake_data2D(qx, qy=None): 
    7181    from sans.dataloader.data_info import Data2D, Detector 
    72  
    7382 
    7483    if qy is None: 
    … …  
    121130            data.mask &= (data.x<outer) 
    122131 
     132 
    123133def set_half(data, half): 
    124134    from sans.dataloader.manipulations import Boxcut 
    … …  
    128138        data.mask += Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data) 
    129139 
     140 
    130141def set_top(data, max): 
    131142    from sans.dataloader.manipulations import Boxcut 
    132143    data.mask += Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=max)(data) 
     144 
    133145 
    134146def plot_data(data, iq, vmin=None, vmax=None): 
    … …  
    141153               interpolation='nearest', aspect=1, origin='upper', 
    142154               extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax) 
     155 
    143156 
    144157def plot_result2D(data, theory, view='linear'): 
    … …  
    184197    mresid = masked_array((theory-data.y)/data.dy, mdata.mask) 
    185198 
    186     plt.subplot(1,2,1) 
     199    plt.subplot(121) 
    187200    plt.errorbar(data.x, mdata, yerr=data.dy) 
    188201    plt.plot(data.x, mtheory, '-', hold=True) 
    189202    plt.yscale(view) 
    190     plt.subplot(1, 2, 2) 
     203    plt.subplot(122) 
    191204    plt.plot(data.x, mresid, 'x') 
    192205    #plt.axhline(1, color='black', ls='--',lw=1, hold=True) 
    … …  
    219232 
    220233        # create model 
    221         self.fn = model(input, cutoff=cutoff) 
     234        self.fn = model(input) 
     235        self.cutoff = cutoff 
    222236 
    223237        # define bumps parameters 
    224238        pars = [] 
    225         for p in model.meta['parameters']: 
     239        extras = [] 
     240        for p in model.info['parameters']: 
    226241            name, default, limits, ptype = p[0], p[2], p[3], p[4] 
    227242            value = kw.pop(name, default) 
    228243            setattr(self, name, Parameter.default(value, name=name, limits=limits)) 
    229244            pars.append(name) 
    230             if ptype != "": 
    231                 for xpart,xdefault,xlimits in [ 
    232                         ('_pd', 0, limits), 
    233                         ('_pd_n', 35, (0,1000)), 
    234                         ('_pd_nsigma', 3, (0,10)), 
    235                         ]: 
    236                     xname = name+xpart 
    237                     xvalue = kw.pop(xname, xdefault) 
    238                     setattr(self, xname, Parameter.default(xvalue, name=xname)) 
     245        for name in model.info['partype']['pd-2d']: 
     246            for xpart,xdefault,xlimits in [ 
     247                    ('_pd', 0, limits), 
     248                    ('_pd_n', 35, (0,1000)), 
     249                    ('_pd_nsigma', 3, (0, 10)), 
     250                    ('_pd_type', 'gaussian', None), 
     251                ]: 
     252                xname = name+xpart 
     253                xvalue = kw.pop(xname, xdefault) 
     254                if xlimits is not None: 
     255                    xvalue = Parameter.default(xvalue, name=xname, limits=xlimits) 
    239256                    pars.append(xname) 
     257                setattr(self, xname, xvalue) 
     258        self._parameter_names = pars 
    240259        if kw: 
    241260            raise TypeError("unexpected parameters: %s"%(", ".join(sorted(kw.keys())))) 
    242         self._parameter_names = pars 
    243261        self.update() 
    244262 
    … …  
    254272    def theory(self): 
    255273        if 'theory' not in self._cache: 
    256             pars = dict((k,getattr(self,k).value) for k in self._parameter_names) 
     274            pars = [getattr(self,p).value for p in self.fn.fixed_pars] 
     275            pd_pars = [self._get_weights(p) for p in self.fn.pd_pars] 
    257276            #print pars 
    258             self._theory[self.index] = self.fn.eval(pars) 
    259             #self._theory[:] = self.fn.eval(pars) 
     277            self._theory[self.index] = self.fn(pars, pd_pars, self.cutoff) 
     278            #self._theory[:] = self.fn.eval(pars, pd_pars) 
    260279            self._cache['theory'] = self._theory 
    261280        return self._cache['theory'] 
    … …  
    282301        pass 
    283302 
     303    def _get_weights(self, par): 
     304        from . import weights 
     305 
     306        relative = self.fn.info['partype']['pd-rel'] 
     307        limits = self.fn.info['limits'] 
     308        disperser,value,npts,width,nsigma = [getattr(self, par+ext) 
     309                for ext in ('_pd_type','','_pd_n','_pd','_pd_nsigma')] 
     310        v,w = weights.get_weights( 
     311            disperser, int(npts.value), width.value, nsigma.value, 
     312            value.value, limits[par], par in relative) 
     313        return v,w/w.max() 
     314 
     315 
    284316def demo(): 
    285317    data = load_data('DEC07086.DAT') 
    … …  
    288320    import matplotlib.pyplot as plt; plt.show() 
    289321 
     322 
    290323if __name__ == "__main__": 
    291324    demo() 

  • sasmodels/dll.py

    r14de349 rce27e21  
    1919    ctypes wrapper for a single model. 
    2020 
    21     *source* and *meta* are the model source and interface as returned 
     21    *source* and *info* are the model source and interface as returned 
    2222    from :func:`gen.make`. 
    2323 
    … …  
    2727    is an optional extension which may not be available on all devices. 
    2828    """ 
    29     def __init__(self, dllpath, meta): 
    30         self.meta = meta 
    31         self.dll = ct.CDLL(dllpath) 
    32         self.Iq = self.dll[gen.kernel_name(self.meta, False)] 
    33         self.Iqxy = self.dll[gen.kernel_name(self.meta, True)] 
     29    def __init__(self, dllpath, info): 
     30        self.info = info 
     31        self.dllpath = dllpath 
     32        self.dll = None 
    3433 
     34    def _load_dll(self): 
     35        Nfixed1d = len(self.info['partype']['fixed-1d']) 
     36        Nfixed2d = len(self.info['partype']['fixed-2d']) 
     37        Npd1d = len(self.info['partype']['pd-1d']) 
     38        Npd2d = len(self.info['partype']['pd-2d']) 
    3539 
    36         self.PARS = dict((p[0],p[2]) for p in meta['parameters']) 
    37         self.PD_PARS = [p[0] for p in meta['parameters'] if p[4] != ""] 
     40        self.dll = ct.CDLL(self.dllpath) 
    3841 
    39         # Determine the set of fixed and polydisperse parameters 
    40         Nfixed = len([p[0] for p in meta['parameters'] if p[4] == ""]) 
    41         N1D = len([p for p in meta['parameters'] if p[4]=="volume"]) 
    42         N2D = len([p for p in meta['parameters'] if p[4]!=""]) 
    43         self.Iq.argtypes = IQ_ARGS + [c_double]*Nfixed + [c_int]*N1D 
    44         self.Iqxy.argtypes = IQXY_ARGS + [c_double]*Nfixed + [c_int]*N2D 
     42        self.Iq = self.dll[gen.kernel_name(self.info, False)] 
     43        self.Iq.argtypes = IQ_ARGS + [c_double]*Nfixed1d + [c_int]*Npd1d 
    4544 
    46     def __call__(self, input, cutoff=1e-5): 
     45        self.Iqxy = self.dll[gen.kernel_name(self.info, True)] 
     46        self.Iqxy.argtypes = IQXY_ARGS + [c_double]*Nfixed2d + [c_int]*Npd2d 
     47 
     48    def __getstate__(self): 
     49        return {'info': self.info, 'dllpath': self.dllpath, 'dll': None} 
     50 
     51    def __setstate__(self, state): 
     52        self.__dict__ = state 
     53 
     54    def __call__(self, input): 
     55        if self.dll is None: self._load_dll() 
     56 
    4757        kernel = self.Iqxy if input.is_2D else self.Iq 
    48         return DllKernel(kernel, self.meta, input, cutoff) 
     58        return DllKernel(kernel, self.info, input) 
    4959 
    5060    def make_input(self, q_vectors): 
    … …  
    90100 
    91101class DllKernel(object): 
    92     def __init__(self, kernel, meta, input, cutoff): 
    93         self.cutoff = cutoff 
     102    def __init__(self, kernel, info, input): 
    94103        self.input = input 
    95104        self.kernel = kernel 
    96         self.meta = meta 
     105        self.info = info 
     106        self.res = np.empty(input.nq, input.dtype) 
     107        dim = '2d' if input.is_2D else '1d' 
     108        self.fixed_pars = info['partype']['fixed-'+dim] 
     109        self.pd_pars = info['partype']['pd-'+dim] 
    97110 
    98         self.res = np.empty(input.nq, input.dtype) 
     111        # In dll kernel, but not in opencl kernel 
    99112        self.p_res = self.res.ctypes.data 
    100113 
    101         # Determine the set of fixed and polydisperse parameters 
    102         self.fixed_pars = [p[0] for p in meta['parameters'] if p[4] == ""] 
    103         self.pd_pars = [p for p in meta['parameters'] 
    104                if p[4]=="volume" or (p[4]=="orientation" and input.is_2D)] 
     114    def __call__(self, pars, pd_pars, cutoff): 
     115        real = np.float32 if self.input.dtype == F32 else np.float64 
     116        fixed = [real(p) for p in pars] 
     117        cutoff = real(cutoff) 
     118        loops = np.hstack(pd_pars) 
     119        loops = np.ascontiguousarray(loops.T, self.input.dtype).flatten() 
     120        loops_N = [np.uint32(len(p[0])) for p in pd_pars] 
    105121 
    106     def eval(self, pars): 
    107         fixed, loops, loop_n = \ 
    108             gen.kernel_pars(pars, self.meta, self.input.is_2D, dtype=self.input.dtype) 
    109         real = np.float32 if self.input.dtype == F32 else np.float64 
    110122        nq = c_int(self.input.nq) 
    111         cutoff = real(self.cutoff) 
    112  
    113123        p_loops = loops.ctypes.data 
    114         pars = self.input.q_pointers + [self.p_res, nq, p_loops, cutoff] + fixed + loop_n 
     124        args = self.input.q_pointers + [self.p_res, nq, p_loops, cutoff] + fixed + loops_N 
    115125        #print pars 
    116         self.kernel(*pars) 
     126        self.kernel(*args) 
    117127 
    118128        return self.res 
    … …  
    120130    def release(self): 
    121131        pass 
    122  
    123     def __del__(self): 
    124         self.release() 

  • sasmodels/gen.py

    r14de349 rce27e21  
    99F64 = np.dtype('float64') 
    1010F32 = np.dtype('float32') 
     11 
     12# Scale and background, which are parameters common to every form factor 
     13COMMON_PARAMETERS = [ 
     14    [ "scale", "", 1, [0, np.inf], "", "Source intensity" ], 
     15    [ "background", "1/cm", 0, [0, np.inf], "", "Source background" ], 
     16    ] 
     17 
    1118 
    1219# Conversion from units defined in the parameter table for each model 
    … …  
    2936 
    3037PARTABLE_VALUE_WIDTH = 10 
    31  
    32 # Scale and background, which are parameters common to every form factor 
    33 COMMON_PARAMETERS = [ 
    34     [ "scale", "", 0, [0, np.inf], "", "Source intensity" ], 
    35     [ "background", "1/cm", 0, [0, np.inf], "", "Source background" ], 
    36     ] 
    37  
    3838 
    3939# Header included before every kernel. 
    … …  
    6060#  define kernel 
    6161#  define SINCOS(angle,svar,cvar) do {svar=sin(angle);cvar=cos(angle);} while (0) 
     62#  define powr(a,b) pow(a,b) 
    6263#else 
    6364#  ifdef USE_SINCOS 
    … …  
    9394# respectively, so the template builder will need to do extra work to 
    9495# declare, initialize and pass the q parameters. 
    95 IQ_KERNEL = { 
     96KERNEL_1D = { 
    9697    'fn': "Iq", 
    9798    'q_par_decl': "global const real *q,", 
    … …  
    100101    } 
    101102 
    102 IQXY_KERNEL = { 
     103KERNEL_2D = { 
    103104    'fn': "Iqxy", 
    104105    'q_par_decl': "global const real *qx,\n    global const real *qy,", 
    … …  
    176177# Volume normalization. 
    177178# If there are "volume" polydispersity parameters, then these will be used 
    178 # to call the volume function from the user supplied kernel, and accumulate 
     179# to call the form_volume function from the user supplied kernel, and accumulate 
    179180# a normalized weight. 
    180181VOLUME_NORM="""const real vol_weight = %(weight)s; 
    181   vol += vol_weight*volume(%(pars)s); 
     182  vol += vol_weight*form_volume(%(pars)s); 
    182183  norm_vol += vol_weight;""" 
     184 
    183185 
    184186def indent(s, depth): 
    … …  
    191193 
    192194 
    193 def make_kernel(meta, form): 
     195def kernel_name(info, is_2D): 
     196    return info['name'] + "_" + ("Iqxy" if is_2D else "Iq") 
     197 
     198 
     199def make_kernel(info, is_2D): 
    194200    """ 
    195201    Build a kernel call from metadata supplied by the user. 
    196202 
    197     *meta* is the json object defined in the kernel file. 
     203    *info* is the json object defined in the kernel file. 
    198204 
    199205    *form* is either "Iq" or "Iqxy". 
    … …  
    206212    # If we are building the Iqxy kernel, we need to propagate qx,qy 
    207213    # parameters, otherwise we can 
    208     if form == "Iqxy": 
    209         qpars = IQXY_KERNEL 
    210     else: 
    211         qpars = IQ_KERNEL 
    212  
     214    dim = "2d" if is_2D else "1d" 
     215    fixed_pars = info['partype']['fixed-'+dim] 
     216    pd_pars = info['partype']['pd-'+dim] 
     217    vol_pars = info['partype']['volume'] 
     218    q_pars = KERNEL_2D if is_2D else KERNEL_1D 
     219 
     220    # Build polydispersity loops 
    213221    depth = 4 
    214222    offset = "" 
    215223    loop_head = [] 
    216224    loop_end = [] 
    217     vol_pars = [] 
    218     fixed_pars = [] 
    219     pd_pars = [] 
    220     fn_pars = [] 
    221     for i,p in enumerate(meta['parameters']): 
    222         name = p[0] 
    223         ptype = p[4] 
    224         if ptype == "volume": 
    225             vol_pars.append(name) 
    226         elif ptype == "orientation": 
    227             if form != "Iqxy": continue  # no orientation for 1D kernels 
    228         elif ptype == "magnetic": 
    229             raise NotImplementedError("no magnetic parameters yet") 
    230         if name not in ['scale','background']: fn_pars.append(name) 
    231         if ptype == "": 
    232             fixed_pars.append(name) 
    233             continue 
    234         else: 
    235             pd_pars.append(name) 
     225    for name in pd_pars: 
    236226        subst = { 'name': name, 'offset': offset } 
    237227        loop_head.append(indent(LOOP_OPEN%subst, depth)) 
    … …  
    254244 
    255245    # Define the inner loop function call 
     246    # The parameters to the f(q,p1,p2...) call should occur in the same 
     247    # order as given in the parameter info structure.  This may be different 
     248    # from the parameter order in the call to the kernel since the kernel 
     249    # call places all fixed parameters before all polydisperse parameters. 
     250    fq_pars = [p[0] for p in info['parameters'][len(COMMON_PARAMETERS):] 
     251               if p[0] in set(fixed_pars+pd_pars)] 
    256252    subst = { 
    257253        'weight_product': "*".join(p+"_w" for p in pd_pars), 
    258254        'volume_norm': volume_norm, 
    259         'fn': qpars['fn'], 
    260         'qcall': qpars['qcall'], 
    261         'pcall': ", ".join(fn_pars), 
     255        'fn': q_pars['fn'], 
     256        'qcall': q_pars['qcall'], 
     257        'pcall': ", ".join(fq_pars), # skip scale and background 
    262258        } 
    263259    loop_body = [indent(LOOP_BODY%subst, depth)] 
    … …  
    280276    subst = { 
    281277        # kernel name is, e.g., cylinder_Iq 
    282         'name': "_".join((meta['name'], qpars['fn'])), 
     278        'name': kernel_name(info, is_2D), 
    283279        # to declare, e.g., global real q[], 
    284         'q_par_decl': qpars['q_par_decl'], 
     280        'q_par_decl': q_pars['q_par_decl'], 
    285281        # to declare, e.g., real sld, int Nradius, int Nlength 
    286282        'par_decl': par_decl, 
    … …  
    288284        'pd_length': "+".join('N'+p for p in pd_pars), 
    289285        # the q initializers, e.g., real qi = q[i]; 
    290         'qinit': qpars['qinit'], 
     286        'qinit': q_pars['qinit'], 
    291287        # the actual polydispersity loop 
    292288        'loops': loops, 
    … …  
    295291    return kernel 
    296292 
    297 def make_partable(meta): 
    298     pars = meta['parameters'] 
     293def make_partable(info): 
     294    pars = info['parameters'] 
    299295    column_widths = [ 
    300296        max(len(p[0]) for p in pars), 
    … …  
    320316    return "\n".join(lines) 
    321317 
    322 def make_doc(kernelfile, meta, doc): 
    323     doc = doc%{'parameters': make_partable(meta)} 
     318def make_doc(kernelfile, info, doc): 
     319    doc = doc%{'parameters': make_partable(info)} 
    324320    return doc 
    325321 
    326 def make_model(kernelfile, meta, source): 
    327     kernel_Iq = make_kernel(meta, "Iq") 
    328     kernel_Iqxy = make_kernel(meta, "Iqxy") 
     322def make_model(kernelfile, info, source): 
     323    kernel_Iq = make_kernel(info, is_2D=False) 
     324    kernel_Iqxy = make_kernel(info, is_2D=True) 
    329325    path = os.path.dirname(kernelfile) 
    330     extra = [open("%s/%s"%(path,f)).read() for f in meta['include']] 
     326    extra = [open("%s/%s"%(path,f)).read() for f in info['include']] 
    331327    kernel = "\n\n".join([KERNEL_HEADER]+extra+[source, kernel_Iq, kernel_Iqxy]) 
    332328    return kernel 
    … …  
    339335    if len(parts) != 3: 
    340336        raise ValueError("PARAMETERS block missing from %r"%kernelfile) 
    341     meta_source = parts[1].strip() 
     337    info_source = parts[1].strip() 
    342338    try: 
    343         meta = relaxed_loads(meta_source) 
     339        info = relaxed_loads(info_source) 
    344340    except: 
    345341        print "in json text:" 
    346342        print "\n".join("%2d: %s"%(i+1,s) 
    347                         for i,s in enumerate(meta_source.split('\n'))) 
     343                        for i,s in enumerate(info_source.split('\n'))) 
    348344        raise 
    349345        #raise ValueError("PARAMETERS block could not be parsed from %r"%kernelfile) 
    350     meta['parameters'] = COMMON_PARAMETERS + meta['parameters'] 
    351     meta['filename'] = kernelfile 
    352346 
    353347    # select documentation out of the source file 
    354348    parts = source.split("DOCUMENTATION") 
    355349    if len(parts) == 3: 
    356         doc = make_doc(kernelfile, meta, parts[1].strip()) 
     350        doc = make_doc(kernelfile, info, parts[1].strip()) 
    357351    elif len(parts) == 1: 
    358352        raise ValueError("DOCUMENTATION block is missing from %r"%kernelfile) 
    … …  
    360354        raise ValueError("DOCUMENTATION block incorrect from %r"%kernelfile) 
    361355 
    362     return source, meta, doc 
     356    return source, info, doc 
     357 
     358def categorize_parameters(pars): 
     359    """ 
     360    Build parameter categories out of the the parameter definitions. 
     361 
     362    Returns a dictionary of categories. 
     363 
     364    The function call sequence consists of q inputs and the return vector, 
     365    followed by the loop value/weight vector, followed by the values for 
     366    the non-polydisperse parameters, followed by the lengths of the 
     367    polydispersity loops.  To construct the call for 1D models, the 
     368    categories *fixed-1d* and *pd-1d* list the names of the parameters 
     369    of the non-polydisperse and the polydisperse parameters respectively. 
     370    Similarly, *fixed-2d* and *pd-2d* provide parameter names for 2D models. 
     371    The *pd-rel* category is a set of those parameters which give 
     372    polydispersitiy as a portion of the value (so a 10% length dispersity 
     373    would use a polydispersity value of 0.1) rather than absolute 
     374    dispersity such as an angle plus or minus 15 degrees. 
     375 
     376    The *volume* category lists the volume parameters in order for calls 
     377    to volume within the kernel (used for volume normalization) and for 
     378    calls to ER and VR for effective radius and volume ratio respectively. 
     379 
     380    The *orientation* and *magnetic* categories list the orientation and 
     381    magnetic parameters.  These are used by the sasview interface.  The 
     382    blank category is for parameters such as scale which don't have any 
     383    other marking. 
     384    """ 
     385    partype = { 
     386        'volume': [], 'orientation': [], 'magnetic': [], '': [], 
     387        'fixed-1d': [], 'fixed-2d': [], 'pd-1d': [], 'pd-2d': [], 
     388        'pd-rel': set(), 
     389    } 
     390 
     391    for p in pars: 
     392        name,ptype = p[0],p[4] 
     393        if ptype == 'volume': 
     394            partype['pd-1d'].append(name) 
     395            partype['pd-2d'].append(name) 
     396            partype['pd-rel'].add(name) 
     397        elif ptype == 'magnetic': 
     398            partype['fixed-2d'].append(name) 
     399        elif ptype == 'orientation': 
     400            partype['pd-2d'].append(name) 
     401        elif ptype == '': 
     402            partype['fixed-1d'].append(name) 
     403            partype['fixed-2d'].append(name) 
     404        else: 
     405            raise ValueError("unknown parameter type %r"%ptype) 
     406        partype[ptype].append(name) 
     407 
     408    return partype 
    363409 
    364410def make(kernelfile): 
    … …  
    370416    """ 
    371417    #print kernelfile 
    372     source, meta, doc = parse_file(kernelfile) 
    373     doc = make_doc(kernelfile, meta, doc) 
    374     model = make_model(kernelfile, meta, source) 
    375     return model, meta, doc 
    376  
    377  
    378  
    379 # Convert from python float to C float or double, depending on dtype 
    380 FLOAT_CONVERTER = { 
    381     F32: np.float32, 
    382     F64: np.float64, 
    383     } 
    384  
    385 def kernel_name(meta, is_2D): 
    386     return meta['name'] + "_" + ("Iqxy" if is_2D else "Iq") 
    387  
    388  
    389 def kernel_pars(pars, par_info, is_2D, dtype=F32): 
    390     """ 
    391     Convert parameter dictionary into arguments for the kernel. 
    392  
    393     *pars* is a dictionary of *{ name: value }*, with *name_pd* for the 
    394     polydispersity width, *name_pd_n* for the number of pd steps, and 
    395     *name_pd_nsigma* for the polydispersity limits. 
    396  
    397     *par_info* is the parameter info structure from the kernel metadata. 
    398  
    399     *is_2D* is True if the dataset represents 2D data, with the corresponding 
    400     orientation parameters. 
    401  
    402     *dtype* is F32 or F64, the numpy single and double precision floating 
    403     point dtypes.  These should not be the strings. 
    404     """ 
    405     from .weights import GaussianDispersion 
    406     real = np.float32 if dtype == F32 else np.float64 
    407     fixed = [] 
    408     parts = [] 
    409     for p in par_info['parameters']: 
    410         name, ptype = p[0],p[4] 
    411         value = pars[name] 
    412         if ptype == "": 
    413             fixed.append(real(value)) 
    414         elif is_2D or ptype != "orientation": 
    415             limits, width = p[3], pars[name+'_pd'] 
    416             n, nsigma = pars[name+'_pd_n'], pars[name+'_pd_nsigma'] 
    417             relative = (ptype != "orientation") 
    418             dist = GaussianDispersion(int(n), width, nsigma) 
    419             # Make sure that weights are normalized to peaks at 1 so that 
    420             # the tolerance term can be used properly on truncated distributions 
    421             v,w = dist.get_weights(value, limits[0], limits[1], relative) 
    422             parts.append((v, w/w.max())) 
    423     loops = np.hstack(parts) 
    424     loops = np.ascontiguousarray(loops.T, dtype).flatten() 
    425     loopsN = [np.uint32(len(p[0])) for p in parts] 
    426     return fixed, loops, loopsN 
     418    source, info, doc = parse_file(kernelfile) 
     419    info['filename'] = kernelfile 
     420    info['parameters'] = COMMON_PARAMETERS + info['parameters'] 
     421    info['partype'] = categorize_parameters(info['parameters']) 
     422    info['limits'] = dict((p[0],p[3]) for p in info['parameters']) 
     423    doc = make_doc(kernelfile, info, doc) 
     424    model = make_model(kernelfile, info, source) 
     425    return model, info, doc 
    427426 
    428427 
    … …  
    437436def demo(): 
    438437    from os.path import join as joinpath, dirname 
    439     c, meta, doc = make_model(joinpath(dirname(__file__), "models", "cylinder.c")) 
     438    c, info, doc = make_model(joinpath(dirname(__file__), "models", "cylinder.c")) 
    440439    #print doc 
    441440    #print c 

  • sasmodels/gpu.py

    r14de349 rce27e21  
    4545#define real double 
    4646""" 
     47 
     48# The max loops number is limited by the amount of local memory available 
     49# on the device.  You don't want to make this value too big because it will 
     50# waste resources, nor too small because it may interfere with users trying 
     51# to do their polydispersity calculations.  A value of 1024 should be much 
     52# larger than necessary given that cost grows as npts^k where k is the number 
     53# of polydisperse parameters. 
     54MAX_LOOPS = 1024 
    4755 
    4856ENV = None 
    … …  
    96104    dtype = np.dtype(dtype) 
    97105    if dtype==F64 and not all(has_double(d) for d in context.devices): 
    98         warnings.warn(RuntimeWarning("Double precision not support; using single precision instead")) 
    99         dtype = F32 
     106        raise RuntimeError("Double precision not supported for devices") 
    100107 
    101108    header = F64_DEFS if dtype == F64 else F32_DEFS 
    … …  
    104111        header += "#define USE_SINCOS\n" 
    105112    program  = cl.Program(context, header+source).build() 
    106     return program, dtype 
     113    return program 
    107114 
    108115 
    … …  
    125132        self.boundary = max(d.min_data_type_align_size 
    126133                            for d in self.context.devices) 
     134        self.has_double = all(has_double(d) for d in self.context.devices) 
     135        self.compiled = {} 
     136 
     137    def compile_program(self, name, source, dtype): 
     138        if name not in self.compiled: 
     139            #print "compiling",name 
     140            self.compiled[name] = compile_model(self.context, source, dtype) 
     141        return self.compiled[name] 
     142 
     143    def release_program(self, name): 
     144        if name in self.compiled: 
     145            self.compiled[name].release() 
     146            del self.compiled[name] 
    127147 
    128148class GpuModel(object): 
    … …  
    130150    GPU wrapper for a single model. 
    131151 
    132     *source* and *meta* are the model source and interface as returned 
     152    *source* and *info* are the model source and interface as returned 
    133153    from :func:`gen.make`. 
    134154 
    … …  
    138158    is an optional extension which may not be available on all devices. 
    139159    """ 
    140     def __init__(self, source, meta, dtype=F32): 
    141         context = environment().context 
    142         self.meta = meta 
    143         self.program, self.dtype = compile_model(context, source, dtype) 
    144  
    145         #TODO: need better interface 
    146         self.PARS = dict((p[0],p[2]) for p in meta['parameters']) 
    147         self.PD_PARS = [p[0] for p in meta['parameters'] if p[4] != ""] 
    148  
    149     def __call__(self, input, cutoff=1e-5): 
     160    def __init__(self, source, info, dtype=F32): 
     161        self.info = info 
     162        self.source = source 
     163        self.dtype = dtype 
     164        self.program = None # delay program creation 
     165 
     166    def __getstate__(self): 
     167        state = self.__dict__.copy() 
     168        state['program'] = None 
     169        return state 
     170 
     171    def __setstate__(self, state): 
     172        self.__dict__ = state.copy() 
     173 
     174    def __call__(self, input): 
    150175        if self.dtype != input.dtype: 
    151176            raise TypeError("data and kernel have different types") 
    152         kernel_name = gen.kernel_name(self.meta, input.is_2D) 
     177        if self.program is None: 
     178            self.program = environment().compile_program(self.info['name'],self.source, self.dtype) 
     179        kernel_name = gen.kernel_name(self.info, input.is_2D) 
    153180        kernel = getattr(self.program, kernel_name) 
    154         return GpuKernel(kernel, self.meta, input, cutoff) 
     181        return GpuKernel(kernel, self.info, input) 
     182 
     183    def release(self): 
     184        if self.program is not None: 
     185            environment().release_program(self.info['name']) 
     186            self.program = None 
    155187 
    156188    def make_input(self, q_vectors): 
    … …  
    210242 
    211243class GpuKernel(object): 
    212     def __init__(self, kernel, meta, input, cutoff): 
    213         env = environment() 
    214  
    215         self.cutoff = cutoff 
     244    def __init__(self, kernel, info, input): 
    216245        self.input = input 
    217246        self.kernel = kernel 
    218         self.meta = meta 
     247        self.info = info 
     248        self.res = np.empty(input.nq, input.dtype) 
     249        dim = '2d' if input.is_2D else '1d' 
     250        self.fixed_pars = info['partype']['fixed-'+dim] 
     251        self.pd_pars = info['partype']['pd-'+dim] 
    219252 
    220253        # Inputs and outputs for each kernel call 
     254        # Note: res may be shorter than res_b if global_size != nq 
     255        env = environment() 
    221256        self.loops_b = [cl.Buffer(env.context, mf.READ_WRITE, 
    222                                   1024*input.dtype.itemsize) 
     257                                  MAX_LOOPS*input.dtype.itemsize) 
    223258                        for _ in env.queues] 
    224259        self.res_b = [cl.Buffer(env.context, mf.READ_WRITE, 
    … …  
    226261                      for _ in env.queues] 
    227262 
    228         # Note: may be shorter than res_b if global_size != nq 
    229         self.res = np.empty(input.nq, input.dtype) 
    230  
    231         # Determine the set of fixed and polydisperse parameters 
    232         self.fixed_pars = [p[0] for p in meta['parameters'] if p[4] == ""] 
    233         self.pd_pars = [p for p in meta['parameters'] 
    234                if p[4]=="volume" or (p[4]=="orientation" and input.is_2D)] 
    235  
    236     def eval(self, pars): 
     263 
     264    def __call__(self, pars, pd_pars, cutoff=1e-5): 
     265        real = np.float32 if self.input.dtype == F32 else np.float64 
     266        fixed = [real(p) for p in pars] 
     267        cutoff = real(cutoff) 
     268        loops = np.hstack(pd_pars) 
     269        loops = np.ascontiguousarray(loops.T, self.input.dtype).flatten() 
     270        loops_N = [np.uint32(len(p[0])) for p in pd_pars] 
     271 
     272        #print "opencl eval",pars 
     273        if len(loops) > MAX_LOOPS: 
     274            raise ValueError("too many polydispersity points") 
    237275        device_num = 0 
    238276        res_bi = self.res_b[device_num] 
    239277        queuei = environment().queues[device_num] 
    240         fixed, loops, loop_n = \ 
    241             gen.kernel_pars(pars, self.meta, self.input.is_2D, dtype=self.input.dtype) 
     278        loops_bi = self.loops_b[device_num] 
    242279        loops_l = cl.LocalMemory(len(loops.data)) 
    243         real = np.float32 if self.input.dtype == F32 else np.float64 
    244         cutoff = real(self.cutoff) 
    245  
    246         loops_bi = self.loops_b[device_num] 
    247280        cl.enqueue_copy(queuei, loops_bi, loops) 
    248281        #ctx = environment().context 
    249282        #loops_bi = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=loops) 
    250         pars = self.input.q_buffers + [res_bi,loops_bi,loops_l,cutoff] + fixed + loop_n 
    251         self.kernel(queuei, self.input.global_size, None, *pars) 
     283        args = self.input.q_buffers + [res_bi,loops_bi,loops_l,cutoff] + fixed + loops_N 
     284        self.kernel(queuei, self.input.global_size, None, *args) 
    252285        cl.enqueue_copy(queuei, self.res, res_bi) 
    253286 

  • sasmodels/models/cylinder.c

    r14de349 rce27e21  
    140140DOCUMENTATION END 
    141141*/ 
    142 real volume(real radius, real length); 
     142real form_volume(real radius, real length); 
    143143real Iq(real q, real sld, real solvent_sld, real radius, real length); 
    144144real Iqxy(real qx, real qy, real sld, real solvent_sld, real radius, real length, real theta, real phi); 
    145145 
    146146 
    147 real volume(real radius, real length) 
     147real form_volume(real radius, real length) 
    148148{ 
    149149    return M_PI*radius*radius*length; 
    150150} 
    151  
    152151 
    153152real Iq(real q, 
    … …  
    170169    // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 
    171170    // The additional volume factor is for polydisperse volume normalization. 
    172     const real s = (sld - solvent_sld) * volume(radius, length); 
     171    const real s = (sld - solvent_sld) * form_volume(radius, length); 
    173172    return REAL(1.0e-4) * form * s * s; 
    174173} 
    … …  
    203202    const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); 
    204203    const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg); 
    205     const real form = REAL(4.0)*bj*bj*si*si; 
     204    const real form = bj*bj*si*si; 
    206205 
    207206    // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 
    208207    // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 
    209208    // The additional volume factor is for polydisperse volume normalization. 
    210     const real s = (sld - solvent_sld) * volume(radius, length); 
     209    const real s = (sld - solvent_sld) * form_volume(radius, length); 
    211210    return REAL(1.0e-4) * form * s * s; // * correction; 
    212211} 

  • sasmodels/weights.py

    r14de349 rce27e21  
     1# TODO: include dispersion docs with the disperser models 
     2from math import sqrt 
    13import numpy as np 
     4from scipy.special import gammaln 
    25 
    3 class GaussianDispersion(object): 
    4     def __init__(self, npts=35, width=0, nsigmas=3): #number want, percent deviation, #standard deviations from mean 
    5         self.type = 'gaussian' 
    6         self.npts = npts 
    7         self.width = width 
    8         self.nsigmas = nsigmas 
     6class Dispersion(object): 
     7    """ 
     8    Base dispersion object. 
     9 
     10    Subclasses should define *_weights(center, sigma, lb, ub)* 
     11    which returns the x points and their corresponding weights. 
     12    """ 
     13    type = "base disperser" 
     14    default = dict(npts=35, width=0, nsigmas=3) 
     15    def __init__(self, npts=None, width=None, nsigmas=None): 
     16        self.npts = self.default['npts'] if npts is None else npts 
     17        self.width = self.default['width'] if width is None else width 
     18        self.nsigmas = self.default['nsigmas'] if nsigmas is None else nsigmas 
    919 
    1020    def get_pars(self): 
    11         return self.__dict__ 
     21        pars = {'type': self.type} 
     22        pars.update(self.__dict__) 
     23        return pars 
    1224 
    13     def get_weights(self, center, min, max, relative): 
    14         """ *center* is the center of the distribution 
    15         *min*,*max* are the min, max allowed values 
    16         *relative* is True if the width is relative to the center instead of absolute 
    17         For polydispersity use relative.  For orientation parameters use absolute.""" 
     25    def set_weights(self, values, weights): 
     26        raise RuntimeError("set_weights is only available for ArrayDispersion") 
     27 
     28    def get_weights(self, center, lb, ub, relative): 
     29        """ 
     30        Return the weights for the distribution. 
     31 
     32        *center* is the center of the distribution 
     33 
     34        *lb*,*ub* are the min and max allowed values 
     35 
     36        *relative* is True if the distribution width is proportional to the 
     37        center value instead of absolute.  For polydispersity use relative. 
     38        For orientation parameters use absolute. 
     39        """ 
    1840        npts, width, nsigmas = self.npts, self.width, self.nsigmas 
    19         sigma = width * center if relative else width 
     41        sigma = self.width * center if relative else self.width 
    2042        if sigma == 0: 
    21             return np.array([center],'d'), np.array([1.], 'd') 
    22         x = center + np.linspace(-nsigmas * sigma, +nsigmas * sigma, npts) 
    23         x = x[(x >= min) & (x <= max)] 
     43            return np.array([center], 'd'), np.array([1.], 'd') 
     44        return self._weights(center, sigma, lb, ub) 
     45 
     46    def _weights(self, center, sigma, lb, ub): 
     47        """actual work of computing the weights""" 
     48        raise NotImplementedError 
     49 
     50    def _linspace(self, center, sigma, lb, ub): 
     51        """helper function to provide linear spaced weight points within range""" 
     52        npts, nsigmas = self.npts, self.nsigmas 
     53        x = center + np.linspace(-nsigmas*sigma, +nsigmas*sigma, npts) 
     54        x = x[(x >= lb) & (x <= ub)] 
     55        return x 
     56 
     57 
     58class GaussianDispersion(Dispersion): 
     59    type = "gaussian" 
     60    default = dict(npts=35, width=0, nsigmas=3) 
     61    def _weights(self, center, sigma, lb, ub): 
     62        x = self._linspace(center, sigma, lb, ub) 
    2463        px = np.exp((x-center)**2 / (-2.0 * sigma * sigma)) 
    2564        return x, px 
    2665 
     66 
     67class RectangleDispersion(Dispersion): 
     68    type = "rectangle" 
     69    default = dict(npts=35, width=0, nsigmas=1.70325) 
     70    def _weights(self, center, sigma, lb, ub): 
     71        x = self._linspace(center, sigma*sqrt(3.0), lb, ub) 
     72        px = np.ones_like(x) 
     73        return x, px 
     74 
     75 
     76class LogNormalDispersion(Dispersion): 
     77    type = "lognormal" 
     78    default = dict(npts=80, width=0, nsigmas=8) 
     79    def _weights(self, center, sigma, lb, ub): 
     80        x = self._linspace(center, sigma, max(lb,1e-8), max(ub,1e-8)) 
     81        px = np.exp(-0.5*(np.log(x)-center)**2)/sigma**2/(x*sigma) 
     82        return x, px 
     83 
     84 
     85class SchulzDispersion(Dispersion): 
     86    type = "schulz" 
     87    default = dict(npts=80, width=0, nsigmas=8) 
     88    def _weights(self, center, sigma, lb, ub): 
     89        x = self._linspace(center, sigma, max(lb,1e-8), max(ub,1e-8)) 
     90        R= x/center 
     91        z = (center/sigma)**2 
     92        arg = z*np.log(z) + (z-1)*np.log(R) - R*z - np.log(center) - gammaln(z) 
     93        px = np.exp(arg) 
     94        return x, px 
     95 
     96 
     97class ArrayDispersion(Dispersion): 
     98    type = "array" 
     99    default = dict(npts=35, width=0, nsigmas=1) 
     100    def __init__(self, npts=None, width=None, nsigmas=None): 
     101        Dispersion.__init__(self, npts, width, nsigmas) 
     102        self.values = np.array([0.], 'd') 
     103        self.weights = np.array([1.], 'd') 
     104 
     105    def set_weights(self, values, weights): 
     106        self.values = np.ascontiguousarray(values, 'd') 
     107        self.weights = np.ascontiguousarray(weights, 'd') 
     108        self.npts = len(values) 
     109 
     110    def _weights(self, center, sigma, lb, ub): 
     111        # TODO: interpolate into the array dispersion using npts 
     112        x = center + self.values*sigma 
     113        idx = (x>=lb)&(x<=ub) 
     114        x = x[idx] 
     115        px = self.weights[idx] 
     116        return x, px 
     117 
     118 
     119models = dict((d.type,d) for d in ( 
     120    GaussianDispersion, RectangleDispersion, 
     121    ArrayDispersion, SchulzDispersion, LogNormalDispersion 
     122)) 
     123 
     124def get_weights(disperser, n, width, nsigma, value, limits, relative): 
     125    cls = models[disperser] 
     126    obj = cls(n, width, nsigma) 
     127    v,w = obj.get_weights(value, limits[0], limits[1], relative) 
     128    return v,w 
Note: See TracChangeset for help on using the changeset viewer.