Changeset 5d4777d in sasmodels


Ignore:
Timestamp:
Sep 2, 2014 1:24:38 AM (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:
f4cf580
Parents:
ff7119b
Message:

reorganize, check and update models

Files:
10 added
1 deleted
14 edited

Legend:

Unmodified
Added
Removed

  • compare-new.py

    rff7119b r5d4777d  
    77 
    88from sasmodels.bumps_model import BumpsModel, plot_data, tic 
    9 from sasmodels.gen import opencl_model, dll_model 
     9from sasmodels import gpu, dll 
     10from sasmodels.convert import revert_model 
     11 
    1012 
    1113def sasview_model(modelname, **pars): 
    … …  
    1315    Load a sasview model given the model name. 
    1416    """ 
    15     modelname = modelname+"Model" 
     17    modelname = modelname 
    1618    sans = __import__('sans.models.'+modelname) 
    1719    ModelClass = getattr(getattr(sans.models,modelname,None),modelname,None) 
    … …  
    3638    sasmodels = __import__('sasmodels.models.'+modelname) 
    3739    module = getattr(sasmodels.models, modelname, None) 
    38     kernel = opencl_model(module, dtype=dtype) 
     40    kernel = gpu.load_model(module, dtype=dtype) 
    3941    return kernel 
    4042 
    … …  
    4244    sasmodels = __import__('sasmodels.models.'+modelname) 
    4345    module = getattr(sasmodels.models, modelname, None) 
    44     kernel = dll_model(module, dtype=dtype) 
     46    kernel = dll.load_model(module, dtype=dtype) 
    4547    return kernel 
    4648 
    47  
    48 def compare(Ncpu, cpuname, cpupars, Ngpu, gpuname, gpupars): 
    49  
     49def randomize(p, v): 
     50    """ 
     51    Randomizing pars using sasview names. 
     52    Guessing at angles and slds. 
     53    """ 
     54    if any(p.endswith(s) for s in ('_pd_n','_pd_nsigmas','_pd_type')): 
     55        return v 
     56    if any(s in p for s in ('theta','phi','psi')): 
     57        return np.random.randint(0,90) 
     58    elif any(s in p for s in ('sld','Sld','SLD')): 
     59        return np.random.rand()*1e-5 
     60    else: 
     61        return np.random.rand() 
     62 
     63def parlist(pars): 
     64    return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items())) 
     65 
     66def suppress_pd(pars): 
     67    """ 
     68    Suppress theta_pd for now until the normalization is resolved. 
     69 
     70    May also suppress complete polydispersity of the model to test 
     71    models more quickly. 
     72    """ 
     73    for p in pars: 
     74        if p.endswith("_pd"): pars[p] = 0 
     75 
     76def compare(gpuname, gpupars, Ncpu, Ngpu, opts): 
    5077    #from sasmodels.core import load_data 
    5178    #data = load_data('December/DEC07098.DAT') 
    52     from sasmodels.core import empty_data1D 
    53     data = empty_data1D(np.logspace(-4, -1, 128)) 
    54     #from sasmodels.core import empty_2D, set_beam_stop 
    55     #data = empty_data2D(np.linspace(-0.05, 0.05, 128)) 
    56     #set_beam_stop(data, 0.004) 
     79    qmax = 1.0 if '-highq' in opts else (0.2 if '-midq' in opts else 0.05) 
     80    if "-1d" in opts: 
     81        from sasmodels.bumps_model import empty_data1D 
     82        qmax = np.log10(qmax) 
     83        data = empty_data1D(np.logspace(qmax-3, qmax, 128)) 
     84    else: 
     85        from sasmodels.bumps_model import empty_data2D, set_beam_stop 
     86        data = empty_data2D(np.linspace(-qmax, qmax, 128)) 
     87        set_beam_stop(data, 0.004) 
    5788    is2D = hasattr(data, 'qx_data') 
    58  
     89    dtype = 'double' if '-double' in opts else 'single' 
     90    cutoff_opts = [s for s in opts if s.startswith('-cutoff')] 
     91    cutoff = float(cutoff_opts[0].split('=')[1]) if cutoff_opts else 1e-5 
     92 
     93 
     94    if '-random' in opts: 
     95        gpupars = dict((p,randomize(p,v)) for p,v in gpupars.items()) 
     96    if '-pars' in opts: print "pars",parlist(gpupars) 
     97    if '-mono' in opts: suppress_pd(gpupars) 
     98 
     99    cpuname, cpupars = revert_model(gpuname, gpupars) 
     100 
     101    try: 
     102        gpumodel = load_opencl(gpuname, dtype=dtype) 
     103    except Exception,exc: 
     104        print exc 
     105        print "... trying again with single precision" 
     106        gpumodel = load_opencl(gpuname, dtype='single') 
     107    model = BumpsModel(data, gpumodel, cutoff=cutoff, **gpupars) 
    59108    if Ngpu > 0: 
    60         gpumodel = load_opencl(gpuname, dtype='single') 
    61         model = BumpsModel(data, gpumodel, **gpupars) 
    62109        toc = tic() 
    63110        for i in range(Ngpu): 
    … …  
    66113            gpu = model.theory() 
    67114        gpu_time = toc()*1000./Ngpu 
    68         print "ocl t=%.1f ms, intensity=%.0f"%(gpu_time, sum(gpu[~np.isnan(gpu)])) 
     115        print "ocl t=%.1f ms, intensity=%f"%(gpu_time, sum(gpu[~np.isnan(gpu)])) 
    69116        #print max(gpu), min(gpu) 
    70117 
    71     if 0 and Ncpu > 0: # Hack to compare ctypes vs. opencl 
    72         dllmodel = load_dll(gpuname, dtype='double') 
    73         model = BumpsModel(data, dllmodel, **gpupars) 
     118    comp = None 
     119    if Ncpu > 0 and "-dll" in opts: 
     120        dllmodel = load_dll(gpuname, dtype=dtype) 
     121        model = BumpsModel(data, dllmodel, cutoff=cutoff, **gpupars) 
    74122        toc = tic() 
    75123        for i in range(Ncpu): 
    … …  
    77125            cpu = model.theory() 
    78126        cpu_time = toc()*1000./Ncpu 
    79         print "dll t=%.1f ms"%cpu_time 
     127        comp = "dll" 
    80128 
    81129    elif 0: # Hack to check new vs old for GpuCylinder 
    82130        from Models.code_cylinder_f import GpuCylinder as oldgpu 
    83131        from sasmodel import SasModel 
    84         oldmodel = SasModel(data, oldgpu, dtype='single', **cpupars) 
     132        oldmodel = SasModel(data, oldgpu, dtype=dtype, **cpupars) 
    85133        toc = tic() 
    86134        for i in range(Ngpu): 
    … …  
    88136            cpu = oldmodel.theory() 
    89137        cpu_time = toc()*1000./Ngpu 
    90         print "old t=%.1f ms"%cpu_time 
     138        comp = "old" 
    91139 
    92140    elif Ncpu > 0: 
    … …  
    100148                cpu = cpumodel.evalDistribution(data.x) 
    101149        cpu_time = toc()*1000./Ncpu 
    102         print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[model.index])) 
    103  
     150        comp = 'sasview' 
     151        #print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[model.index])) 
     152 
     153    if comp: 
     154        print "%s t=%.1f ms, intensity=%f"%(comp, cpu_time, sum(cpu[model.index])) 
    104155    if Ngpu > 0 and Ncpu > 0: 
    105         print "gpu/cpu", max(abs(gpu/cpu)), "%.15g"%max(abs(gpu)), "%.15g"%max(abs(cpu)) 
     156        #print "speedup %.2g"%(cpu_time/gpu_time) 
     157        #print "max |gpu/cpu|", max(abs(gpu/cpu)), "%.15g"%max(abs(gpu)), "%.15g"%max(abs(cpu)) 
    106158        #cpu *= max(gpu/cpu) 
    107         abserr = (gpu - cpu) 
    108         relerr = (gpu - cpu)/cpu 
    109         print "max(|ocl-omp|)", max(abs(abserr[model.index])) 
    110         print "max(|(ocl-omp)/ocl|)", max(abs(relerr[model.index])) 
    111     #return 
     159        resid, relerr = np.zeros_like(gpu), np.zeros_like(gpu) 
     160        resid[model.index] = (gpu - cpu)[model.index] 
     161        relerr[model.index] = resid[model.index]/cpu[model.index] 
     162        print "max(|ocl-%s|)"%comp, max(abs(resid[model.index])) 
     163        print "max(|(ocl-%s)/ocl|)"%comp, max(abs(relerr[model.index])) 
     164 
     165    if '-noplot' in opts: return 
    112166 
    113167    import matplotlib.pyplot as plt 
    … …  
    115169        if Ngpu > 0: plt.subplot(131) 
    116170        plot_data(data, cpu, scale='log') 
    117         plt.title("omp t=%.1f ms"%cpu_time) 
     171        plt.title("%s t=%.1f ms"%(comp,cpu_time)) 
    118172    if Ngpu > 0: 
    119173        if Ncpu > 0: plt.subplot(132) 
    … …  
    122176    if Ncpu > 0 and Ngpu > 0: 
    123177        plt.subplot(133) 
    124         plot_data(data, 1e8*relerr, scale='linear') 
    125         plt.title("max rel err = %.3g"%max(abs(relerr))) 
     178        err = resid if '-abs' in opts else relerr 
     179        plot_data(data, err, scale='linear') 
     180        plt.title("max rel err = %.3g"%max(abs(err[model.index]))) 
    126181        if is2D: plt.colorbar() 
    127182    plt.show() 
    128183 
    129 def rename(pars, **names): 
    130     newpars = pars.copy() 
    131     for new,old in names.items(): 
    132         for variant in ("", "_pd", "_pd_n", "_pd_nsigma"): 
    133             if old+variant in newpars: 
    134                 newpars[new+variant] = pars[old+variant] 
    135                 del newpars[old+variant] 
    136     return newpars 
    137  
    138 def rescale_sld(pars, names): 
    139     newpars = pars.copy() 
    140     for p in names: 
    141         newpars[p] *= 1e6 
    142     return newpars 
    143  
    144184# =========================================================================== 
    145185# 
     186USAGE=""" 
     187usage: compare.py model [Nopencl] [Nsasview] [options...] [key=val] 
     188 
     189Compare the speed and value for a model between the SasView original and the 
     190OpenCL rewrite. 
     191 
     192model is the name of the model to compare (see below). 
     193Nopencl is the number of times to run the OpenCL model (default=5) 
     194Nsasview is the number of times to run the Sasview model (default=1) 
     195 
     196Options (* for default): 
     197 
     198    -plot*/-noplot plots or suppress the plot of the model 
     199    -single/-double uses double precision for comparison 
     200    -lowq/-midq/-highq use q values up to 0.05, 0.2 or 1.0 
     201    -2d/-1d uses 1d or 2d random data 
     202    -preset/-random randomizes the parameters 
     203    -poly/-mono force monodisperse/polydisperse 
     204    -sasview/-dll whether cpu is tested using sasview or dll 
     205    -cutoff=1e-5/value cutoff for including a point in polydispersity 
     206    -nopars/-pars prints the parameter set or not 
     207    -rel/-abs plot relative or absolute error 
     208 
     209Key=value pairs allow you to set specific values to any of the model 
     210parameters. 
     211 
     212Available models: 
     213 
     214    %s 
     215""" 
     216 
     217def main(): 
     218    opts = [arg for arg in sys.argv[1:] if arg.startswith('-')] 
     219    args = [arg for arg in sys.argv[1:] if not arg.startswith('-')] 
     220    models = "\n    ".join("%-7s: %s"%(k,v.__name__.replace('_',' ')) 
     221                           for k,v in sorted(MODELS.items())) 
     222    if len(args) == 0: 
     223        print(USAGE%models) 
     224        sys.exit(1) 
     225    if args[0] not in MODELS: 
     226        print "Model %r not available. Use one of:\n    %s"%(args[0],models) 
     227        sys.exit(1) 
     228 
     229    name, pars = MODELS[args[0]]() 
     230    Nopencl = int(args[1]) if len(args) > 1 else 5 
     231    Nsasview = int(args[2]) if len(args) > 3 else 1 
     232 
     233    # Fill in default polydispersity parameters 
     234    pds = set(p.split('_pd')[0] for p in pars if p.endswith('_pd')) 
     235    for p in pds: 
     236        if p+"_pd_nsigma" not in pars: pars[p+"_pd_nsigma"] = 3 
     237        if p+"_pd_type" not in pars: pars[p+"_pd_type"] = "gaussian" 
     238 
     239    # Fill in parameters given on the command line 
     240    for arg in args[3:]: 
     241        k,v = arg.split('=') 
     242        if k not in pars: 
     243            # extract base name without distribution 
     244            # style may be either a.d or a_pd_d 
     245            s = set((p.split('_pd')[0]).split('.')[0] for p in pars) 
     246            print "%r invalid; parameters are: %s"%(k,", ".join(sorted(s))) 
     247            sys.exit(1) 
     248        pars[k] = float(v) if not v.endswith('type') else v 
     249 
     250    compare(name, pars, Nsasview, Nopencl, opts) 
     251 
     252# =========================================================================== 
     253# 
     254 
    146255MODELS = {} 
    147256def model(name): 
    … …  
    151260    return gather_function 
    152261 
    153 USAGE=""" 
    154 usage: compare model [Nopencl] [Nsasview] 
    155  
    156 Compare the speed and value for a model between the SasView original and the 
    157 OpenCL rewrite. 
    158  
    159 * Nopencl is the number of times to run the OpenCL model (default=5) 
    160  
    161 * Nsasview is the number of times to run the Sasview model (default=1) 
    162  
    163 * model is the name of the model to compare: 
    164  
    165     %s 
    166 """ 
    167  
    168 def main(): 
    169     if len(sys.argv) == 1: 
    170         models = "\n    ".join("%-7s: %s"%(k,v.__name__.replace('_',' ')) 
    171                                for k,v in sorted(MODELS.items())) 
    172         print(USAGE%models) 
    173         sys.exit(1) 
    174  
    175     cpuname, cpupars, gpuname, gpupars = MODELS[sys.argv[1]]() 
    176     Nopencl = int(sys.argv[2]) if len(sys.argv) > 2 else 5 
    177     Nsasview = int(sys.argv[3]) if len(sys.argv) > 3 else 1 
    178  
    179     compare(Nsasview, cpuname, cpupars, Nopencl, gpuname, gpupars) 
    180262 
    181263@model('cyl') 
    182264def cylinder(): 
    183     cpupars = dict( 
    184         scale=.003, background=.1, 
    185         sldCyl=.291e-6, sldSolv=5.77e-6, 
    186         radius=264.1, length=66.96, 
    187         cyl_theta=85, cyl_phi=0, 
    188         radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 
    189         length_pd=0.1,length_pd_n=1, length_pd_nsigma=3, 
    190         cyl_theta_pd=45, cyl_theta_pd_n=50, cyl_theta_pd_nsigma=3, 
    191         cyl_phi_pd=0.1, cyl_phi_pd_n=5, cyl_phi_pd_nsigma=3, 
    192         ) 
    193     cpuname = 'Cylinder' 
    194  
    195     gpupars = rename(cpupars, theta='cyl_theta', phi='cyl_phi', sld='sldCyl', solvent_sld='sldSolv') 
    196     gpupars = rescale_sld(gpupars, ['sld', 'solvent_sld']) 
    197     gpuname = 'cylinder' 
    198     return cpuname, cpupars, gpuname, gpupars 
     265    pars = dict( 
     266        scale=1, background=0, 
     267        sld=6e-6, solvent_sld=1e-6, 
     268        #radius=5, length=20, 
     269        radius=260, length=290, 
     270        theta=30, phi=15, 
     271        radius_pd=.2, radius_pd_n=1, 
     272        length_pd=.2,length_pd_n=1, 
     273        theta_pd=15, theta_pd_n=25, 
     274        phi_pd=15, phi_pd_n=1, 
     275        ) 
     276    return 'cylinder', pars 
     277 
     278@model('capcyl') 
     279def capped_cylinder(): 
     280    pars = dict( 
     281        scale=1, background=0, 
     282        sld=6e-6, solvent_sld=1e-6, 
     283        radius=260, cap_radius=80000, length=290, 
     284        theta=30, phi=15, 
     285        radius_pd=.2, radius_pd_n=1, 
     286        cap_radius_pd=.2, cap_radius_pd_n=1, 
     287        length_pd=.2, length_pd_n=1, 
     288        theta_pd=15, theta_pd_n=25, 
     289        phi_pd=15, phi_pd_n=1, 
     290        ) 
     291    return 'capped_cylinder', pars 
     292 
     293 
     294@model('cscyl') 
     295def core_shell_cylinder(): 
     296    pars = dict( 
     297        scale=1, background=0, 
     298        core_sld=6e-6, shell_sld=8e-6, solvent_sld=1e-6, 
     299        radius=325, thickness=25, length=34.2709, 
     300        theta=90, phi=0, 
     301        radius_pd=0.1, radius_pd_n=10, 
     302        length_pd=0.1, length_pd_n=5, 
     303        thickness_pd=0.1, thickness_pd_n=5, 
     304        theta_pd=15.8, theta_pd_n=5, 
     305        phi_pd=0.0008748, phi_pd_n=5, 
     306        ) 
     307    return 'core_shell_cylinder', pars 
    199308 
    200309 
    201310@model('ell') 
    202 def ellipse(): 
    203     pars = dict( 
    204         scale=.027, background=4.9, 
    205         sldEll=.297e-6, sldSolv=5.773e-6, 
    206         radius_a=60, radius_b=180, 
    207         axis_theta=0, axis_phi=90, 
    208         radius_a_pd=0.1, radius_a_pd_n=10, radius_a_pd_nsigma=3, 
    209         radius_b_pd=0.1, radius_b_pd_n=10, radius_b_pd_nsigma=3, 
    210         axis_theta_pd=0.1, axis_theta_pd_n=6, axis_theta_pd_nsigma=3, 
    211         axis_phi_pd=0.1, axis_phi_pd_n=6, axis_phi_pd_nsigma=3, 
    212         ) 
    213  
    214     from Models.code_ellipse import GpuEllipse as gpumodel 
    215     model = sasview_model('Ellipsoid', **pars) 
    216  
    217     pars = rename(pars, theta='axis_theta', phi='axis_phi', sld='sldEll', solvent_sld='sldSolv') 
    218     pars = rescale_sld(pars, ['sld', 'solvent_sld']) 
    219     return model, gpumodel, pars 
    220  
    221  
    222 @model('cscyl') 
    223 def core_shell_cylinder(N=1): 
    224     pars = dict( 
    225         scale= 1.77881e-06, background=223.827, 
    226         core_sld=1e-6, shell_sld=.291e-6, solvent_sld=7.105e-6, 
    227         radius=325, thickness=25, length=34.2709, 
    228         axis_theta=90, axis_phi=0, 
    229         radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 
    230         length_pd=0.1, length_pd_n=10, length_pd_nsigma=3, 
    231         thickness_pd=0.1, thickness_pd_n=5, thickness_pd_nsigma=3, 
    232         axis_theta_pd=15.8, axis_theta_pd_n=20, axis_theta_pd_nsigma=5, 
    233         axis_phi_pd=0.0008748, axis_phi_pd_n=5, axis_phi_pd_nsigma=3, 
    234         ) 
    235  
    236     model = sasview_model('CoreShellCylinder', **pars) 
    237     from Models.code_coreshellcyl_f import GpuCoreShellCylinder as gpumodel 
    238  
    239     pars = rename(pars, theta='axis_theta', phi='axis_phi') 
    240     pars = rescale_sld(pars, ['core_sld', 'shell_sld', 'solvent_sld']) 
    241     return model, gpumodel, pars 
     311def ellipsoid(): 
     312    pars = dict( 
     313        scale=1, background=0, 
     314        sld=6e-6, solvent_sld=1e-6, 
     315        rpolar=50, requatorial=30, 
     316        theta=0, phi=0, 
     317        rpolar_pd=0.3, rpolar_pd_n=10, 
     318        requatorial_pd=0, requatorial_pd_n=10, 
     319        theta_pd=0, theta_pd_n=45, 
     320        phi_pd=0, phi_pd_n=45, 
     321        ) 
     322    return 'ellipsoid', pars 
    242323 
    243324 
    244325@model('ell3') 
    245 def triaxial_ellipse(N=1): 
    246     pars = dict( 
    247         scale=0.08, background=5, 
    248         sldEll=7.105e-6, sldSolv=.291e-6, 
    249         axis_theta=0, axis_phi=0, axis_psi=0, 
    250         semi_axisA=15, semi_axisB=20, semi_axisC=500, 
    251         axis_theta_pd=20, axis_theta_pd_n=10, axis_theta_pd_nsigma=3, 
    252         axis_phi_pd=.1, axis_phi_pd_n=10, axis_phi_pd_nsigma=3, 
    253         axis_psi_pd=30, axis_psi_pd_n=5, axis_psi_pd_nsigma=3, 
    254         semi_axisA_pd=.1, semi_axisA_pd_n=5, semi_axisA_pd_nsigma=3, 
    255         semi_axisB_pd=.1, semi_axisB_pd_n=5, semi_axisB_pd_nsigma=3, 
    256         semi_axisC_pd=.1, semi_axisC_pd_n=5, semi_axisC_pd_nsigma=3, 
    257         ) 
    258  
    259     model = sasview_model('TriaxialEllipsoid', **pars) 
    260     from Models.code_triaxialellipse import GpuTriEllipse as gpumodel 
    261  
    262     pars = rename(pars, 
    263                   theta='axis_theta', phi='axis_phi', psi='axis_psi', 
    264                   sld='sldEll', solvent_sld='sldSolv', 
    265                   radius_a='semi_axisA', radius_b='semi_axisB', 
    266                   radius_c='semi_axisC', 
    267                   ) 
    268     pars = rescale_sld(pars, ['sld', 'solvent_sld']) 
    269     return model, gpumodel, pars 
     326def triaxial_ellipsoid(): 
     327    pars = dict( 
     328        scale=1, background=0, 
     329        sld=6e-6, solvent_sld=1e-6, 
     330        theta=0, phi=0, psi=0, 
     331        req_minor=25, req_major=36, rpolar=50, 
     332        theta_pd=0, theta_pd_n=5, 
     333        phi_pd=0, phi_pd_n=5, 
     334        psi_pd=0, psi_pd_n=5, 
     335        req_minor_pd=0, req_minor_pd_n=5, 
     336        req_major_pd=0, req_major_pd_n=5, 
     337        rpolar_pd=.3, rpolar_pd_n=25, 
     338        ) 
     339    return 'triaxial_ellipsoid', pars 
     340 
     341@model('sph') 
     342def sphere(): 
     343    pars = dict( 
     344        scale=1, background=0, 
     345        sld=6e-6, solvent_sld=1e-6, 
     346        radius=120, 
     347        radius_pd=.3, radius_pd_n=5, 
     348        ) 
     349    return 'sphere', pars 
    270350 
    271351@model('lam') 
    272 def lamellar(N=1): 
    273     pars = dict( 
    274         scale=0.08, background=0.003, 
    275         sld_bi=5.38e-6,sld_sol=7.105e-6, 
    276         bi_thick=19.2946, 
    277         bi_thick_pd= 0.37765, bi_thick_pd_n=40, bi_thick_pd_nsigma=3, 
    278         ) 
    279  
    280     model = sasview_model('Lamellar', **pars) 
    281     from Models.code_lamellar import GpuLamellar as gpumodel 
    282  
    283     pars = rename(pars, sld='sld_bi', solvent_sld='sld_sol', thickness='bi_thick') 
    284     pars = rescale_sld(pars, ['sld', 'solvent_sld']) 
    285     return model, gpumodel, pars 
    286  
    287 @model('capcyl') 
    288 def capped_cylinder(N=1): 
    289     pars = dict( 
    290         scale=.08, background=0, 
    291         sld_capcyl=1e-6, sld_solv=6.3e-6, 
    292         rad_cyl=20, rad_cap=40, len_cyl=400, 
    293         theta=0, phi=0, 
    294         rad_cyl_pd=.1, rad_cyl_pd_n=10, rad_cyl_pd_nsigma=3, 
    295         rad_cap_pd=.1, rad_cap_pd_n=10, rad_cap_pd_nsigma=3, 
    296         len_cyl_pd=.1, len_cyl_pd_n=3, len_cyl_pd_nsigma=3, 
    297         theta_pd=.1, theta_pd_n=3, theta_pd_nsigma=3, 
    298         phi_pd=.1, phi_pd_n=3, phi_pd_nsigma=3, 
    299         ) 
    300  
    301  
    302     model = sasview_model('CappedCylinder', **pars) 
    303     from Models.code_capcyl import GpuCapCylinder as gpumodel 
    304  
    305     pars = rename(pars, 
    306                   sld='sld_capcyl', solvent_sld='sld_solv', 
    307                   length='len_cyl', radius='rad_cyl', 
    308                   cap_radius='rad_cap') 
    309     pars = rescale_sld(pars, ['sld', 'solvent_sld']) 
    310     return model, gpumodel, pars 
    311  
     352def lamellar(): 
     353    pars = dict( 
     354        scale=1, background=0, 
     355        sld=6e-6, solvent_sld=1e-6, 
     356        thickness=40, 
     357        thickness_pd= 0.3, thickness_pd_n=40, 
     358        ) 
     359    return 'lamellar', pars 
    312360 
    313361if __name__ == "__main__": 

  • sasmodels/bumps_model.py

    rff7119b r5d4777d  
    140140    import matplotlib.pyplot as plt 
    141141    if hasattr(data, 'qx_data'): 
    142         img = masked_array(iq, data.mask) 
     142        iq = iq[:] 
     143        valid = np.isfinite(iq) 
    143144        if scale == 'log': 
    144             img[(img <= 0) | ~np.isfinite(img)] = masked 
    145             img = np.log10(img) 
     145            valid[valid] = (iq[valid] > 0) 
     146            iq[valid] = np.log10(iq[valid]) 
     147        iq[~valid|data.mask] = 0 
     148        #plottable = iq 
     149        plottable = masked_array(iq, ~valid|data.mask) 
    146150        xmin, xmax = min(data.qx_data), max(data.qx_data) 
    147151        ymin, ymax = min(data.qy_data), max(data.qy_data) 
    148         plt.imshow(img.reshape(128,128), 
     152        plt.imshow(plottable.reshape(128,128), 
    149153                   interpolation='nearest', aspect=1, origin='upper', 
    150154                   extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax) 
    … …  
    154158            plt.plot(data.x[idx], iq[idx]) 
    155159        else: 
    156             idx = np.isfinite(iq) & (iq>0) 
     160            idx = np.isfinite(iq) 
     161            idx[idx] = (iq[idx]>0) 
    157162            plt.loglog(data.x[idx], iq[idx]) 
    158163 
    … …  
    228233    def __init__(self, data, model, cutoff=1e-5, **kw): 
    229234        from bumps.names import Parameter 
     235        partype = model.info['partype'] 
    230236 
    231237        # interpret data 
    … …  
    236242            self.diq = data.err_data[self.index] 
    237243            self._theory = np.zeros_like(data.data) 
    238             q_vectors = [data.qx_data, data.qy_data] 
     244            if not partype['orientation'] and not partype['magnetic']: 
     245                q_vectors = [np.sqrt(data.qx_data**2+data.qy_data**2)] 
     246            else: 
     247                q_vectors = [data.qx_data, data.qy_data] 
    239248        else: 
    240249            self.index = (data.x>=data.qmin) & (data.x<=data.qmax) & ~np.isnan(data.y) 
    … …  
    258267            setattr(self, name, Parameter.default(value, name=name, limits=limits)) 
    259268            pars.append(name) 
    260         for name in model.info['partype']['pd-2d']: 
     269        for name in partype['pd-2d']: 
    261270            for xpart,xdefault,xlimits in [ 
    262271                    ('_pd', 0, limits), 

  • sasmodels/dll.py

    rff7119b r5d4777d  
    22C types wrapper for sasview models. 
    33""" 
    4  
     4import sys 
     5import os 
    56import ctypes as ct 
    67from ctypes import c_void_p, c_int, c_double 
    … …  
    1112 
    1213from .gen import F32, F64 
     14# Compiler platform details 
     15if sys.platform == 'darwin': 
     16    COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp" 
     17elif os.name == 'nt': 
     18    COMPILE = "gcc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm" 
     19else: 
     20    COMPILE = "cc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm" 
     21DLL_PATH = "/tmp" 
     22 
     23 
     24def dll_path(info): 
     25    """ 
     26    Path to the compiled model defined by *info*. 
     27    """ 
     28    from os.path import join as joinpath, split as splitpath, splitext 
     29    basename = splitext(splitpath(info['filename'])[1])[0] 
     30    return joinpath(DLL_PATH, basename+'.so') 
     31 
     32 
     33def load_model(kernel_module, dtype=None): 
     34    """ 
     35    Load the compiled model defined by *kernel_module*. 
     36 
     37    Recompile if any files are newer than the model file. 
     38 
     39    *dtype* is ignored.  Compiled files are always double. 
     40 
     41    The DLL is not loaded until the kernel is called so models an 
     42    be defined without using too many resources. 
     43    """ 
     44    import tempfile 
     45 
     46    source, info = gen.make(kernel_module) 
     47    source_files = gen.sources(info) + [info['filename']] 
     48    newest = max(os.path.getmtime(f) for f in source_files) 
     49    dllpath = dll_path(info) 
     50    if not os.path.exists(dllpath) or os.path.getmtime(dllpath)<newest: 
     51        # Replace with a proper temp file 
     52        fid, filename = tempfile.mkstemp(suffix=".c",prefix="sas_"+info['name']) 
     53        os.fdopen(fid,"w").write(source) 
     54        status = os.system(COMPILE%(filename, dllpath)) 
     55        if status != 0: 
     56            print "compile failed.  File is in %r"%filename 
     57        else: 
     58            ## uncomment the following to keep the generated c file 
     59            #os.unlink(filename); print "saving compiled file in %r"%filename 
     60            pass 
     61    return DllModel(dllpath, info) 
     62 
    1363 
    1464IQ_ARGS = [c_void_p, c_void_p, c_int, c_void_p, c_double] 
    … …  
    124174    """ 
    125175    def __init__(self, kernel, info, input): 
     176        self.info = info 
    126177        self.input = input 
    127178        self.kernel = kernel 

  • sasmodels/gen.py

    rff7119b r5d4777d  
    2828 
    2929*Iq*, *Iqxy*, *Imagnetic* and *form_volume* should be stylized C-99 
    30 functions written for OpenCL.  Floating point values should be 
    31 declared as *real*.  Depending on how the function is called, a macro 
    32 will replace *real* with *float* or *double*.  Unfortunately, MacOSX 
    33 is picky about floating point constants, which should be defined with 
    34 value + 'f' if they are of type *float* or just a bare value if they 
    35 are of type *double*.  The solution is a macro *REAL(value)* which 
    36 adds the 'f' if compiling for single precision floating point.  This 
    37 does make the code ugly, and may someday be replaced by a clever 
    38 regular expression which does the same job.  OpenCL has a *sincos* 
    39 function which can improve performance when both the *sin* and *cos* 
    40 values are needed for a particular argument.  Since this function 
    41 does not exist in C-99, all use of *sincos* should be replaced by the 
    42 macro *SINCOS(value,sn,cn)* where *sn* and *cn* are previously declared 
    43 *real* values.  *value* may be an expression.  When compiled for systems 
    44 without OpenCL, *SINCOS* will be replaced by *sin* and *cos* calls.  All 
    45 functions need prototype declarations even if the are defined before they 
    46 are used -- another present from MacOSX.  OpenCL does not support 
    47 *#include* preprocessor directives; instead the includes must be listed 
    48 in the kernel metadata, with functions defined before they are used. 
    49 The included files should be listed using relative path to the kernel 
    50 source file, or if using one of the standard models, relative to the 
    51 sasmodels source files. 
     30functions written for OpenCL.  All functions need prototype declarations 
     31even if the are defined before they are used.  OpenCL does not support 
     32*#include* preprocessor directives, so instead the list of includes needs 
     33to be given as part of the metadata in the kernel module definition. 
     34The included files should be listed using a path relative to the kernel 
     35module, or if using "lib/file.c" if it is one of the standard includes 
     36provided with the sasmodels source.  The includes need to be listed in 
     37order so that functions are defined before they are used. 
     38 
     39Floating point values should be declared as *real*.  Depending on how the 
     40function is called, a macro will replace *real* with *float* or *double*. 
     41Unfortunately, MacOSX is picky about floating point constants, which 
     42should be defined with value + 'f' if they are of type *float* or just 
     43a bare value if they are of type *double*.  The solution is a macro 
     44*REAL(value)* which adds the 'f' if compiling for single precision 
     45floating point.  [Note: we could write a clever regular expression 
     46which automatically detects real-valued constants.  If we wanted to 
     47make our code more C-like, we could define variables with double but 
     48replace them with float before compiling for single precision.] 
     49 
     50OpenCL has a *sincos* function which can improve performance when both 
     51the *sin* and *cos* values are needed for a particular argument.  Since 
     52this function does not exist in C99, all use of *sincos* should be 
     53replaced by the macro *SINCOS(value,sn,cn)* where *sn* and *cn* are 
     54previously declared *real* values.  *value* may be an expression.  When 
     55compiled for systems without OpenCL, *SINCOS* will be replaced by 
     56*sin* and *cos* calls. 
     57 
     58If the input parameters are invalid, the scattering calculator should 
     59return a negative number. Particularly with polydispersity, there are 
     60some sets of shape parameters which lead to nonsensical forms, such 
     61as a capped cylinder where the cap radius is smaller than the 
     62cylinder radius.  The polydispersity calculation will ignore these points, 
     63effectively chopping the parameter weight distributions at the boundary 
     64of the infeasible region.  The resulting scattering will be set to 
     65background.  This will work correctly even when polydispersity is off. 
    5266 
    5367*ER* and *VR* are python functions which operate on parameter vectors. 
    … …  
    112126    *VR* is a python function defining the volume ratio.  If it is not 
    113127    present, the volume ratio is 1. 
     128 
     129    *form_volume*, *Iq*, *Iqxy*, *Imagnetic* are strings containing the 
     130    C source code for the body of the volume, Iq, and Iqxy functions 
     131    respectively.  These can also be defined in the last source file. 
    114132 
    115133An *info* dictionary is constructed from the kernel meta data and 
    … …  
    156174__all__ = ["make, doc", "sources"] 
    157175 
     176import sys 
     177import os 
    158178import os.path 
    159179 
    … …  
    226246#endif 
    227247 
    228 // Standard mathematical constants, prefixed with M_: 
    229 //   E, LOG2E, LOG10E, LN2, LN10, PI, PI_2, PI_4, 1_PI, 2_PI, 
    230 //   2_SQRTPI, SQRT2, SQRT1_2 
     248// Standard mathematical constants: 
     249//   M_E, M_LOG2E, M_LOG10E, M_LN2, M_LN10, M_PI, M_PI_2=pi/2, M_PI_4=pi/4, 
     250//   M_1_PI=1/pi, M_2_PI=2/pi, M_2_SQRTPI=2/sqrt(pi), SQRT2, SQRT1_2=sqrt(1/2) 
    231251// OpenCL defines M_constant_F for float constants, and nothing if double 
    232252// is not enabled on the card, which is why these constants may be missing 
    … …  
    256276    'qinit': "const real qi = q[i];", 
    257277    'qcall': "qi", 
     278    'qwork': ["q"], 
    258279    } 
    259280 
    … …  
    263284    'qinit': "const real qxi = qx[i];\n    const real qyi = qy[i];", 
    264285    'qcall': "qxi, qyi", 
     286    'qwork': ["qx", "qy"], 
    265287    } 
    266288 
    … …  
    318340for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) { 
    319341  const real %(name)s = loops[2*(%(name)s_i%(offset)s)]; 
    320   const real %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];""" 
     342  const real %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\ 
     343""" 
    321344 
    322345# Polydispersity loop body. 
    … …  
    327350const real weight = %(weight_product)s; 
    328351if (weight > cutoff) { 
    329   ret += weight*%(fn)s(%(qcall)s, %(pcall)s); 
    330   norm += weight; 
    331   %(volume_norm)s 
    332 }""" 
     352  const real I = %(fn)s(%(qcall)s, %(pcall)s); 
     353  if (I>=REAL(0.0)) { // scattering cannot be negative 
     354    ret += weight*I; 
     355    norm += weight; 
     356    %(volume_norm)s 
     357  } 
     358  //else { printf("exclude qx,qy,I:%%g,%%g,%%g\\n",%(qcall)s,I); } 
     359} 
     360//else { printf("exclude weight:%%g\\n",weight); }\ 
     361""" 
     362 
     363# Use this when integrating over orientation 
     364SPHERICAL_CORRECTION="""\ 
     365// Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta 
     366real spherical_correction = (Ntheta>1 ? fabs(cos(M_PI_180*phi)) : REAL(1.0));\ 
     367""" 
    333368 
    334369# Volume normalization. 
    … …  
    337372# a normalized weight. 
    338373VOLUME_NORM="""const real vol_weight = %(weight)s; 
    339   vol += vol_weight*form_volume(%(pars)s); 
    340   norm_vol += vol_weight;""" 
     374    vol += vol_weight*form_volume(%(pars)s); 
     375    norm_vol += vol_weight;\ 
     376""" 
     377 
     378# functions defined as strings in the .py module 
     379WORK_FUNCTION="""\ 
     380real %(name)s(%(pars)s); 
     381real %(name)s(%(pars)s) 
     382{ 
     383%(body)s 
     384}\ 
     385""" 
    341386 
    342387# Documentation header for the module, giving the model name, its short 
    … …  
    393438    vol_pars = info['partype']['volume'] 
    394439    q_pars = KERNEL_2D if is_2D else KERNEL_1D 
     440    fn = q_pars['fn'] 
    395441 
    396442    # Build polydispersity loops 
    … …  
    426472    fq_pars = [p[0] for p in info['parameters'][len(COMMON_PARAMETERS):] 
    427473               if p[0] in set(fixed_pars+pd_pars)] 
     474    if False and "phi" in pd_pars: 
     475        spherical_correction = [indent(SPHERICAL_CORRECTION, depth)] 
     476        weights = [p+"_w" for p in pd_pars]+['spherical_correction'] 
     477    else: 
     478        spherical_correction = [] 
     479        weights = [p+"_w" for p in pd_pars] 
    428480    subst = { 
    429         'weight_product': "*".join(p+"_w" for p in pd_pars), 
     481        'weight_product': "*".join(weights), 
    430482        'volume_norm': volume_norm, 
    431         'fn': q_pars['fn'], 
     483        'fn': fn, 
    432484        'qcall': q_pars['qcall'], 
    433485        'pcall': ", ".join(fq_pars), # skip scale and background 
    434486        } 
    435487    loop_body = [indent(LOOP_BODY%subst, depth)] 
    436     loops = "\n".join(loop_head+loop_body+loop_end) 
     488    loops = "\n".join(loop_head+spherical_correction+loop_body+loop_end) 
    437489 
    438490    # declarations for non-pd followed by pd pars 
    … …  
    465517        } 
    466518    kernel = KERNEL_TEMPLATE%subst 
     519 
     520    # If the working function is defined in the kernel metadata as a 
     521    # string, translate the string to an actual function definition 
     522    # and put it before the kernel. 
     523    if info[fn]: 
     524        subst = { 
     525            'name': fn, 
     526            'pars': ", ".join("real "+p for p in q_pars['qwork']+fq_pars), 
     527            'body': info[fn], 
     528            } 
     529        kernel = "\n".join((WORK_FUNCTION%subst, kernel)) 
    467530    return kernel 
    468531 
    … …  
    514577    search_path = [ dirname(info['filename']), 
    515578                    abspath(joinpath(dirname(__file__),'models')) ] 
    516     return [_search(search_path) for f in info['source']] 
     579    return [_search(search_path, f) for f in info['source']] 
    517580 
    518581def make_model(info): 
    … …  
    521584    found in the given search path. 
    522585    """ 
     586    source = [open(f).read() for f in sources(info)] 
     587    # If the form volume is defined as a string, then wrap it in a 
     588    # function definition and place it after the external sources but 
     589    # before the kernel functions.  If the kernel functions are strings, 
     590    # they will be translated in the make_kernel call. 
     591    if info['form_volume']: 
     592        subst = { 
     593            'name': "form_volume", 
     594            'pars': ", ".join("real "+p for p in info['partype']['volume']), 
     595            'body': info['form_volume'], 
     596            } 
     597        source.append(WORK_FUNCTION%subst) 
    523598    kernel_Iq = make_kernel(info, is_2D=False) 
    524599    kernel_Iqxy = make_kernel(info, is_2D=True) 
    525     source = [open(f).read() for f in sources(info)] 
    526600    kernel = "\n\n".join([KERNEL_HEADER]+source+[kernel_Iq, kernel_Iqxy]) 
    527601    return kernel 
    … …  
    573647        name = kernel_module.name, 
    574648        title = kernel_module.title, 
    575         source = kernel_module.source, 
    576649        description = kernel_module.description, 
    577650        parameters = COMMON_PARAMETERS + kernel_module.parameters, 
    578         ER = getattr(kernel_module, 'ER', None), 
    579         VR = getattr(kernel_module, 'VR', None), 
     651        source = getattr(kernel_module, 'source', []), 
    580652        ) 
     653    # Fill in attributes which default to None 
     654    info.update((k,getattr(kernel_module, k, None)) 
     655                for k in ('ER', 'VR', 'form_volume', 'Iq', 'Iqxy')) 
     656    # Fill in the derived attributes 
    581657    info['limits'] = dict((p[0],p[3]) for p in info['parameters']) 
    582658    info['partype'] = categorize_parameters(info['parameters']) 
    … …  
    596672    return DOC_HEADER%subst 
    597673 
    598 # Compiler platform details 
    599 if sys.platform == 'darwin': 
    600     COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp" 
    601 elif os.name == 'nt': 
    602     COMPILE = "gcc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm" 
    603 else: 
    604     COMPILE = "cc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm" 
    605 DLL_PATH = "/tmp" 
    606  
    607  
    608 def dll_path(info): 
    609     """ 
    610     Path to the compiled model defined by *info*. 
    611     """ 
    612     from os.path import join as joinpath, split as splitpath, splitext 
    613     basename = splitext(splitpath(info['filename'])[1])[0] 
    614     return joinpath(DLL_PATH, basename+'.so') 
    615  
    616  
    617 def dll_model(kernel_module, dtype=None): 
    618     """ 
    619     Load the compiled model defined by *kernel_module*. 
    620  
    621     Recompile if any files are newer than the model file. 
    622  
    623     *dtype* is ignored.  Compiled files are always double. 
    624  
    625     The DLL is not loaded until the kernel is called so models an 
    626     be defined without using too many resources. 
    627     """ 
    628     import tempfile 
    629     from sasmodels import dll 
    630  
    631     source, info = make(kernel_module) 
    632     source_files = sources(info) + [info['filename']] 
    633     newest = max(os.path.getmtime(f) for f in source_files) 
    634     dllpath = dll_path(info) 
    635     if not os.path.exists(dllpath) or os.path.getmtime(dllpath)<newest: 
    636         # Replace with a proper temp file 
    637         srcfile = tempfile.mkstemp(suffix=".c",prefix="sas_"+info['name']) 
    638         open(srcfile, 'w').write(source) 
    639         os.system(COMPILE%(srcfile, dllpath)) 
    640         ## comment the following to keep the generated c file 
    641         os.unlink(srcfile) 
    642     return dll.DllModel(dllpath, info) 
    643  
    644  
    645 def opencl_model(kernel_module, dtype="single"): 
    646     """ 
    647     Load the OpenCL model defined by *kernel_module*. 
    648  
    649     Access to the OpenCL device is delayed until the kernel is called 
    650     so models can be defined without using too many resources. 
    651     """ 
    652     from sasmodels import gpu 
    653  
    654     source, info = make(kernel_module) 
    655     ## for debugging, save source to a .cl file, edit it, and reload as model 
    656     #open(modelname+'.cl','w').write(source) 
    657     #source = open(modelname+'.cl','r').read() 
    658     return gpu.GpuModel(source, info, dtype) 
    659674 
    660675 

  • sasmodels/gpu.py

    rff7119b r5d4777d  
    2525 
    2626""" 
    27 import warnings 
    28  
    2927import numpy as np 
    3028import pyopencl as cl 
    … …  
    3230 
    3331from . import gen 
    34  
    35 from .gen import F32, F64 
    3632 
    3733F32_DEFS = """\ 
    … …  
    5248# larger than necessary given that cost grows as npts^k where k is the number 
    5349# of polydisperse parameters. 
    54 MAX_LOOPS = 1024 
     50MAX_LOOPS = 2048 
     51 
     52def load_model(kernel_module, dtype="single"): 
     53    """ 
     54    Load the OpenCL model defined by *kernel_module*. 
     55 
     56    Access to the OpenCL device is delayed until the kernel is called 
     57    so models can be defined without using too many resources. 
     58    """ 
     59    source, info = gen.make(kernel_module) 
     60    ## for debugging, save source to a .cl file, edit it, and reload as model 
     61    open(info['name']+'.cl','w').write(source) 
     62    #source = open(info['name']+'.cl','r').read() 
     63    return GpuModel(source, info, dtype) 
    5564 
    5665ENV = None 
    … …  
    103112    """ 
    104113    dtype = np.dtype(dtype) 
    105     if dtype==F64 and not all(has_double(d) for d in context.devices): 
     114    if dtype==gen.F64 and not all(has_double(d) for d in context.devices): 
    106115        raise RuntimeError("Double precision not supported for devices") 
    107116 
    108     header = F64_DEFS if dtype == F64 else F32_DEFS 
     117    header = F64_DEFS if dtype == gen.F64 else F32_DEFS 
    109118    # Note: USE_SINCOS makes the intel cpu slower under opencl 
    110119    if context.devices[0].type == cl.device_type.GPU: 
    … …  
    158167    is an optional extension which may not be available on all devices. 
    159168    """ 
    160     def __init__(self, source, info, dtype=F32): 
     169    def __init__(self, source, info, dtype=gen.F32): 
    161170        self.info = info 
    162171        self.source = source 
    … …  
    221230    buffer will be released when the data object is freed. 
    222231    """ 
    223     def __init__(self, q_vectors, dtype=F32): 
     232    def __init__(self, q_vectors, dtype=gen.F32): 
    224233        env = environment() 
    225234        self.nq = q_vectors[0].size 
    … …  
    273282        env = environment() 
    274283        self.loops_b = [cl.Buffer(env.context, mf.READ_WRITE, 
    275                                   MAX_LOOPS*input.dtype.itemsize) 
     284                                  2*MAX_LOOPS*input.dtype.itemsize) 
    276285                        for _ in env.queues] 
    277286        self.res_b = [cl.Buffer(env.context, mf.READ_WRITE, 
    … …  
    281290 
    282291    def __call__(self, pars, pd_pars, cutoff=1e-5): 
    283         real = np.float32 if self.input.dtype == F32 else np.float64 
     292        real = np.float32 if self.input.dtype == gen.F32 else np.float64 
    284293        fixed = [real(p) for p in pars] 
    285294        cutoff = real(cutoff) 
    286295        loops = np.hstack(pd_pars) 
    287296        loops = np.ascontiguousarray(loops.T, self.input.dtype).flatten() 
    288         loops_N = [np.uint32(len(p[0])) for p in pd_pars] 
     297        Nloops = [np.uint32(len(p[0])) for p in pd_pars] 
     298        #print "loops",Nloops, loops 
    289299 
    290300        #import sys; print >>sys.stderr,"opencl eval",pars 
    291301        #print "opencl eval",pars 
    292         if len(loops) > MAX_LOOPS: 
     302        if len(loops) > 2*MAX_LOOPS: 
    293303            raise ValueError("too many polydispersity points") 
    294304        device_num = 0 
    … …  
    300310        #ctx = environment().context 
    301311        #loops_bi = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=loops) 
    302         args = self.input.q_buffers + [res_bi,loops_bi,loops_l,cutoff] + fixed + loops_N 
     312        args = self.input.q_buffers + [res_bi,loops_bi,loops_l,cutoff] + fixed + Nloops 
    303313        self.kernel(queuei, self.input.global_size, None, *args) 
    304314        cl.enqueue_copy(queuei, self.res, res_bi) 

  • sasmodels/models/README.rst

    rff7119b r5d4777d  
    77model. 
    88 
    9 cylindermodel_onefile.py is cylinder.c+cylinder.py merged into one file. 
    10 This doesn't actually run yet since sasmodels.gen has not been updated 
    11 to support it.  It exists as a proposal.  Note that the function declaration 
    12 has been removed since there is enough information in the parameter 
    13 definitions to generate it automatically.  Note also that "source" which 
    14 used to be all the source has been renamed "includes". 
    15  
    16 One-file models could coexist with the py+c file models by checking for the 
    17 existence of c_blah and creating the appropriate function wrappers.  These 
    18 would be appended after any include files.  You shouldn't mix the two forms 
    19 within a single model since form_volume needs to be defined before 
    20 Iq/Iqxy but after the libraries. 
     9lamellar.py is an example of a single file model with embedded C code. 
    2110 
    2211Note: may want to rename form_volume to calc_volume and Iq/Iqxy to 
    23 calc_Iq/calc_Iqxy in both the c+py and the one file forms so that the 
    24 names are more predictable.  Similarly ER/VR go to calc_ER/calc_VR. 
     12calc_Iq/calc_Iqxy. Similarly ER/VR go to calc_ER/calc_VR. 
    2513 
    2614Note: It is possible to translate python code automatically to opencl, using 
    27 something like numba, clyther, shedskin or pypy. 
     15something like numba, clyther, shedskin or pypy, so maybe the kernel functions 
     16could be implemented without any C syntax. 
    2817 
    2918Magnetism hasn't been implemented yet.  We may want a separate Imagnetic 
    … …  
    3524 
    3625Need to write code to generate the polydispersity loops in python for 
    37 kernels that are only implemented in python. 
     26kernels that are only implemented in python.  Also, need to implement 
     27an example kernel directly in python. 

  • sasmodels/models/cylinder.c

    rff7119b r5d4777d  
    11real form_volume(real radius, real length); 
    22real Iq(real q, real sld, real solvent_sld, real radius, real length); 
    3 real Iqxy(real qx, real qy, real sld, real solvent_sld, real radius, real length, real theta, real phi); 
     3real Iqxy(real qx, real qy, real sld, real solvent_sld, 
     4    real radius, real length, real theta, real phi); 
     5 
     6// twovd = 2 * volume * delta_rho 
     7// besarg = q * R * sin(alpha) 
     8// siarg = q * L/2 * cos(alpha) 
     9real _cyl(real twovd, real besarg, real siarg); 
     10real _cyl(real twovd, real besarg, real siarg) 
     11{ 
     12    const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); 
     13    const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg); 
     14    return twovd*si*bj; 
     15} 
    416 
    517real form_volume(real radius, real length) 
    … …  
    1426    real length) 
    1527{ 
    16     const real halflength = REAL(0.5)*length; 
    17     real summ = REAL(0.0); 
     28    const real qr = q*radius; 
     29    const real qh = q*REAL(0.5)*length; 
     30    const real twovd = REAL(2.0)*(sld-solvent_sld)*form_volume(radius, length); 
     31    real total = REAL(0.0); 
    1832    // real lower=0, upper=M_PI_2; 
    1933    for (int i=0; i<76 ;i++) { 
    2034        // translate a point in [-1,1] to a point in [lower,upper] 
    21         //const real zi = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    22         const real zi = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
    23         summ += Gauss76Wt[i] * CylKernel(q, radius, halflength, zi); 
     35        //const real alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
     36        const real alpha = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
     37        real sn, cn; 
     38        SINCOS(alpha, sn, cn); 
     39        const real fq = _cyl(twovd, qr*sn, qh*cn); 
     40        total += Gauss76Wt[i] * fq * fq * sn; 
    2441    } 
    2542    // translate dx in [-1,1] to dx in [lower,upper] 
    26     //const real form = (upper-lower)/2.0*summ; 
    27     const real form = summ * M_PI_4; 
    28  
    29     // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 
    30     // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 
    31     // The additional volume factor is for polydisperse volume normalization. 
    32     const real s = (sld - solvent_sld) * form_volume(radius, length); 
    33     return REAL(1.0e-4) * form * s * s; 
     43    //const real form = (upper-lower)/2.0*total; 
     44    return REAL(1.0e-4) * total * M_PI_4; 
    3445} 
    3546 
    … …  
    4354    real phi) 
    4455{ 
     56    // TODO: check that radius<0 and length<0 give zero scattering. 
     57    // This should be the case since the polydispersity weight vector should 
     58    // be zero length, and this function never called. 
    4559    real sn, cn; // slots to hold sincos function output 
    4660 
    … …  
    5468    const real alpha = acos(cos_val); 
    5569 
    56     // The following is CylKernel() / sin(alpha), but we are doing it in place 
    57     // to avoid sin(alpha)/sin(alpha) for alpha = 0.  It is also a teensy bit 
    58     // faster since we don't mulitply and divide sin(alpha). 
     70    const real twovd = REAL(2.0)*(sld-solvent_sld)*form_volume(radius, length); 
    5971    SINCOS(alpha, sn, cn); 
    60     const real besarg = q*radius*sn; 
    61     const real siarg = REAL(0.5)*q*length*cn; 
    62     // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1 
    63     const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); 
    64     const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg); 
    65     const real form = REAL(4.0)*bj*bj*si*si; 
    66  
    67     // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 
    68     // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 
    69     // The additional volume factor is for polydisperse volume normalization. 
    70     const real s = (sld - solvent_sld) * form_volume(radius, length); 
    71     return REAL(1.0e-4) * form * s * s; // * correction; 
     72    const real fq = _cyl(twovd, q*radius*sn, q*REAL(0.5)*length*cn); 
     73    return REAL(1.0e-4) * fq * fq; // * correction; 
    7274} 

  • sasmodels/models/cylinder.py

    ra7684e5 r5d4777d  
     1# cylinder model 
    12# Note: model title and parameter table are inserted automatically 
    23r""" 
    … …  
    3132To provide easy access to the orientation of the cylinder, we define the 
    3233axis of the cylinder using two angles $\theta$ and $\phi$. Those angles 
    33 are defined in Figure :num:`figure #CylinderModel-orientation`. 
     34are defined in Figure :num:`figure #cylinder-orientation`. 
    3435 
    35 .. _CylinderModel-orientation: 
     36.. _cylinder-orientation: 
    3637 
    3738.. figure:: img/image061.JPG   (should be img/cylinder-1.jpg, or img/cylinder-orientation.jpg) 
    … …  
    6465Validation of our code was done by comparing the output of the 1D model 
    6566to the output of the software provided by the NIST (Kline, 2006). 
    66 Figure :num:`figure #CylinderModel-compare` shows a comparison of 
     67Figure :num:`figure #cylinder-compare` shows a comparison of 
    6768the 1D output of our model and the output of the NIST software. 
    6869 
    69 .. _CylinderModel-compare: 
     70.. _cylinder-compare: 
    7071 
    7172.. figure:: img/image065.JPG 
    … …  
    9192the intensity for fully oriented cylinders, we can compare the result of 
    9293averaging our 2D output using a uniform distribution $p(\theta, \phi) = 1.0$. 
    93 Figure :num:`figure #CylinderModel-crosscheck` shows the result of 
     94Figure :num:`figure #cylinder-crosscheck` shows the result of 
    9495such a cross-check. 
    9596 
    96 .. _CylinderModel-crosscheck: 
     97.. _cylinder-crosscheck: 
    9798 
    9899.. figure:: img/image066.JPG 
    … …  
    111112title = "Right circular cylinder with uniform scattering length density." 
    112113description = """ 
    113      f(q)= 2*(sldCyl - sldSolv)*V*sin(qLcos(alpha/2)) 
     114     P(q)= 2*(sld - solvent_sld)*V*sin(qLcos(alpha/2)) 
    114115            /[qLcos(alpha/2)]*J1(qRsin(alpha/2))/[qRsin(alpha)] 
    115116 
    … …  
    119120            L: Length of the cylinder 
    120121            J1: The bessel function 
    121             alpha: angle betweenthe axis of the 
     122            alpha: angle between the axis of the 
    122123            cylinder and the q-vector for 1D 
    123124            :the ouput is P(q)=scale/V*integral 
    124125            from pi/2 to zero of... 
    125             f(q)^(2)*sin(alpha)*dalpha+ bkg 
    126     """ 
     126            f(q)^(2)*sin(alpha)*dalpha + background 
     127""" 
    127128 
    128129parameters = [ 
    … …  
    143144    ] 
    144145 
    145 source = [ "lib/J1.c", "lib/gauss76.c", "lib/cylkernel.c", "cylinder.c"] 
     146source = [ "lib/J1.c", "lib/gauss76.c", "cylinder.c" ] 
    146147 
    147148def ER(radius, length): 

  • sasmodels/models/cylinder_clone.c

    rff7119b r5d4777d  
    22real Iq(real q, real sld, real solvent_sld, real radius, real length); 
    33real Iqxy(real qx, real qy, real sld, real solvent_sld, real radius, real length, real theta, real phi); 
     4 
     5 
     6// twovd = 2 * volume * delta_rho 
     7// besarg = q * R * sin(alpha) 
     8// siarg = q * L/2 * cos(alpha) 
     9real _cyl(real twovd, real besarg, real siarg, real alpha); 
     10real _cyl(real twovd, real besarg, real siarg, real alpha) 
     11{ 
     12    const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); 
     13    const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg); 
     14    return twovd*si*bj; 
     15} 
    416 
    517real form_volume(real radius, real length) 
    … …  
    719    return M_PI*radius*radius*length; 
    820} 
    9  
    1021real Iq(real q, 
    1122    real sldCyl, 
    … …  
    1425    real length) 
    1526{ 
    16     const real h = REAL(0.5)*length; 
    17     real summ = REAL(0.0); 
     27    const real qr = q*radius; 
     28    const real qh = q*REAL(0.5)*length; 
     29    const real twovd = REAL(2.0)*(sldCyl-sldSolv)*form_volume(radius, length); 
     30    real total = REAL(0.0); 
     31    // real lower=0, upper=M_PI_2; 
    1832    for (int i=0; i<76 ;i++) { 
    19         //const real zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 
    20         const real zi = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
    21         summ += Gauss76Wt[i] * CylKernel(q, radius, h, zi); 
     33        // translate a point in [-1,1] to a point in [lower,upper] 
     34        //const real alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
     35        const real alpha = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
     36        real sn, cn; 
     37        SINCOS(alpha, sn, cn); 
     38        const real fq = _cyl(twovd, qr*sn, qh*cn, alpha); 
     39        total += Gauss76Wt[i] * fq * fq * sn; 
    2240    } 
    23     //const real form = (uplim-lolim)/2.0*summ; 
    24     const real form = summ * M_PI_4; 
    25  
    26     // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 
    27     // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 
    28     // The additional volume factor is for polydisperse volume normalization. 
    29     const real s = (sldCyl - sldSolv) * form_volume(radius, length); 
    30     return REAL(1.0e8) * form * s * s; 
     41    // translate dx in [-1,1] to dx in [lower,upper] 
     42    //const real form = (upper-lower)/2.0*total; 
     43    return REAL(1.0e8) * total * M_PI_4; 
    3144} 
    3245 
    … …  
    5063    const real alpha = acos(cos_val); 
    5164 
    52     // The following is CylKernel() / sin(alpha), but we are doing it in place 
    53     // to avoid sin(alpha)/sin(alpha) for alpha = 0.  It is also a teensy bit 
    54     // faster since we don't mulitply and divide sin(alpha). 
     65    const real qr = q*radius; 
     66    const real qh = q*REAL(0.5)*length; 
     67    const real twovd = REAL(2.0)*(sldCyl-sldSolv)*form_volume(radius, length); 
    5568    SINCOS(alpha, sn, cn); 
    56     const real besarg = q*radius*sn; 
    57     const real siarg = REAL(0.5)*q*length*cn; 
    58     // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1 
    59     const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); 
    60     const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg); 
    61     const real form = REAL(4.0)*bj*bj*si*si; 
    62  
    63     // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 
    64     // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl 
    65     // The additional volume factor is for polydisperse volume normalization. 
    66     const real s = (sldCyl - sldSolv) * form_volume(radius, length); 
    67     return REAL(1.0e8) * form * s * s * spherical_integration; 
     69    const real fq = _cyl(twovd, qr*sn, qh*cn, alpha); 
     70    return REAL(1.0e8) * fq * fq * spherical_integration; 
    6871} 

  • sasmodels/models/cylinder_clone.py

    ra7684e5 r5d4777d  
    150150      "Out of plane angle" ], 
    151151    ] 
    152 source = [ "lib/J1.c", "lib/gauss76.c", "lib/cylkernel.c", "cylinder_clone.c"] 
     152source = [ "lib/J1.c", "lib/gauss76.c", "cylinder_clone.c" ] 
    153153 
    154154def ER(radius, length): 

  • sasmodels/models/cylinder_onefile.py

    ra7684e5 r5d4777d  
    154154source = [ "lib/J1.c", "lib/gauss76.c", "lib/cylkernel.c" ] 
    155155 
    156 c_form_volume = """ 
     156form_volume = """ 
    157157    return M_PI*radius*radius*length; 
    158 """ 
    159  
    160 c_Iq = """ 
     158    """ 
     159 
     160Iq = """ 
    161161    const real h = REAL(0.5)*length; 
    162162    real summ = REAL(0.0); 
    … …  
    174174    const real s = (sld - solvent_sld) * form_volume(radius, length); 
    175175    return REAL(1.0e-4) * form * s * s; 
    176 """ 
    177  
    178 c_Iqxy = """ 
     176    """ 
     177 
     178Iqxy = """ 
    179179    real sn, cn; // slots to hold sincos function output 
    180180 
    … …  
    204204    const real s = (sld - solvent_sld) * form_volume(radius, length); 
    205205    return REAL(1.0e-4) * form * s * s; // * correction; 
    206 """ 
     206    """ 
    207207 
    208208def ER(radius, length): 

  • sasmodels/models/ellipsoid.c

    rce27e21 r5d4777d  
    1 /* PARAMETERS 
     1real form_volume(real rpolar, real requatorial); 
     2real Iq(real q, real sld, real solvent_sld, real rpolar, real requatorial); 
     3real Iqxy(real qx, real qy, real sld, real solvent_sld, 
     4    real rpolar, real requatorial, real theta, real phi); 
     5 
     6real _ellipsoid_kernel(real q, real rpolar, real requatorial, real cos_alpha); 
     7real _ellipsoid_kernel(real q, real rpolar, real requatorial, real cos_alpha) 
    28{ 
    3 name: "ellipsoid", 
    4 title: "Ellipsoid with uniform scattering length density", 
    5 include: [ "lib/gauss76.c" ], 
    6 parameters: [ 
    7    // [ "name", "units", default, [lower, upper], "type", "description" ], 
    8    [ "sld", "1e-6/Ang^2", 4, [-Infinity,Infinity], "", 
    9      "Cylinder scattering length density" ], 
    10    [ "solvent_sld", "1e-6/Ang^2", 1, [-Infinity,Infinity], "", 
    11      "Solvent scattering length density" ], 
    12    [ "a", "Ang",  20, [0, Infinity], "volume", 
    13      "Cylinder radius" ], 
    14    [ "b", "Ang",  20, [0, Infinity], "volume", 
    15      "Cylinder length" ], 
    16    [ "theta", "degrees", 60, [-Infinity, Infinity], "orientation", 
    17      "In plane angle" ], 
    18    [ "phi", "degrees", 60, [-Infinity, Infinity], "orientation", 
    19      "Out of plane angle" ], 
    20 ], 
    21 } 
    22 PARAMETERS END 
    23  
    24 DOCUMENTATION 
    25 .. _EllipseModel: 
    26  
    27 DOCUMENTATION END 
    28 */ 
    29  
    30 real form_volume(real a, real b); 
    31 real Iq(real qx, real qy, real sld, real solvent_sld, real a, real b); 
    32 real Iqxy(real qx, real qy, real sld, real solvent_sld, real a, real b, real theta, real phi); 
    33  
    34 real form_volume(real a, real b) 
    35 { 
    36     return REAL(1.333333333333333)*M_PI_2*a*b*b; 
     9    real sn, cn; 
     10    real ratio = rpolar/requatorial; 
     11    const real u = q*requatorial*sqrt(REAL(1.0) 
     12                   + cos_alpha*cos_alpha*(ratio*ratio - REAL(1.0))); 
     13    SINCOS(u, sn, cn); 
     14    const real f = ( u==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(sn-u*cn)/(u*u*u) ); 
     15    return f*f; 
    3716} 
    3817 
    39 real ellipsoid_kernel(double q, double b, double a, double dum) 
     18real form_volume(real rpolar, real requatorial) 
    4019{ 
    41     real sn, cn; 
    42     const real nu = a/b; 
    43     const real arg = q * b * sqrt(REAL(1.0)+(dum*dum*(nu*nu--REAL(1.0)))); 
    44     SINCOS(arg, sn, cn); 
    45     const real f = (arg==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(sn-arg*cn)/(arg*arg*arg); 
    46     return f*f; 
     20    return REAL(1.333333333333333)*M_PI*rpolar*requatorial*requatorial; 
    4721} 
    4822 
    … …  
    5024    real sld, 
    5125    real solvent_sld, 
    52     real a, 
    53     real b) 
     26    real rpolar, 
     27    real requatorial) 
    5428{ 
    55     real summ = REAL(0.0); 
     29    //const real lower = REAL(0.0); 
     30    //const real upper = REAL(1.0); 
     31    real total = REAL(0.0); 
    5632    for (int i=0;i<76;i++) { 
    57         //const real zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 
    58         zi = ( Gauss76Z[i] + REAL(1.0))/REAL(2.0); 
    59         summ += Gauss76Wt[i] * ellipsoid_kernel(q, b, a, zi); 
     33        //const real cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
     34        const real cos_alpha = REAL(0.5)*(Gauss76Z[i] + REAL(1.0)); 
     35        total += Gauss76Wt[i] * _ellipsoid_kernel(q, rpolar, requatorial, cos_alpha); 
    6036    } 
    61     //const real form = (uplim-lolim)/2.0*summ; 
    62     const real form = REAL(0.5)*summ 
    63     const real s = (sld - sld_solvent) * form_volume(a, b); 
     37    //const real form = (upper-lower)/2*total; 
     38    const real form = REAL(0.5)*total; 
     39    const real s = (sld - solvent_sld) * form_volume(rpolar, requatorial); 
    6440    return REAL(1.0e-4) * form * s * s; 
    6541} 
    … …  
    6844    real sld, 
    6945    real solvent_sld, 
    70     real a, 
    71     real b, 
     46    real rpolar, 
     47    real requatorial, 
    7248    real theta, 
    7349    real phi) 
    … …  
    7753    const real q = sqrt(qx*qx + qy*qy); 
    7854    SINCOS(theta*M_PI_180, sn, cn); 
    79     const real cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
    80     const real form = ellipsoid_kernel(q, b, a, cos_val); 
    81     const real s = (sld - solvent_sld) * form_volume(a, b); 
     55    const real cos_alpha = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 
     56    const real form = _ellipsoid_kernel(q, rpolar, requatorial, cos_alpha); 
     57    const real s = (sld - solvent_sld) * form_volume(rpolar, requatorial); 
    8258 
    8359    return REAL(1.0e-4) * form * s * s; 

  • sasmodels/sasview_model.py

    rff7119b r5d4777d  
    66 
    77try: 
    8     import pyopencl 
    9     from .gen import opencl_model as load_model 
    10 except ImportError: 
     8    from .gpu import load_model 
     9except ImportError,exc: 
     10    warnings.warn(str(exc)) 
    1111    warnings.warn("OpenCL not available --- using ctypes instead") 
    12     from .gen import dll_model as load_model 
     12    from .dll import load_model 
    1313 
    1414 
    … …  
    246246            # Check whether we have a list of ndarrays [qx,qy] 
    247247            qx, qy = qdist 
    248             return self.calculate_Iq(qx, qy) 
     248            partype = self._model.info['partype'] 
     249            if not partype['orientation'] and not partype['magnetic']: 
     250                return self.calculate_Iq(np.sqrt(qx**2+qy**2)) 
     251            else: 
     252                return self.calculate_Iq(qx, qy) 
    249253 
    250254        elif isinstance(qdist, np.ndarray): 

  • sasmodels/weights.py

    rff7119b r5d4777d  
    4242        """ 
    4343        sigma = self.width * center if relative else self.width 
    44         if sigma == 0: 
    45             return np.array([center], 'd'), np.array([1.], 'd') 
     44        if sigma == 0 or self.npts < 2: 
     45            if lb <= center <= ub: 
     46                return np.array([center], 'd'), np.array([1.], 'd') 
     47            else: 
     48                return np.array([], 'd'), np.array([], 'd') 
    4649        return self._weights(center, sigma, lb, ub) 
    4750 
Note: See TracChangeset for help on using the changeset viewer.