source: sasmodels/compare-new.py @ 13d86bc

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 13d86bc was 13d86bc, checked in by Paul Kienzle <pkienzle@…>, 10 years ago

1D model comparison with sasview

  • Property mode set to 100755
File size: 10.5 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3
4import sys
5
6import numpy as np
7
8from sasmodels.core import BumpsModel, plot_data, tic, opencl_model, dll_model
9
10def sasview_model(modelname, **pars):
11    """
12    Load a sasview model given the model name.
13    """
14    modelname = modelname+"Model"
15    sans = __import__('sans.models.'+modelname)
16    ModelClass = getattr(getattr(sans.models,modelname,None),modelname,None)
17    if ModelClass is None:
18        raise ValueError("could not find model %r in sans.models"%modelname)
19    model = ModelClass()
20
21    for k,v in pars.items():
22        if k.endswith("_pd"):
23            model.dispersion[k[:-3]]['width'] = v
24        elif k.endswith("_pd_n"):
25            model.dispersion[k[:-5]]['npts'] = v
26        elif k.endswith("_pd_nsigma"):
27            model.dispersion[k[:-10]]['nsigmas'] = v
28        elif k.endswith("_pd_type"):
29            model.dispersion[k[:-8]]['type'] = v
30        else:
31            model.setParam(k, v)
32    return model
33
34def load_opencl(modelname, dtype='single'):
35    sasmodels = __import__('sasmodels.models.'+modelname)
36    module = getattr(sasmodels.models, modelname, None)
37    kernel = opencl_model(module, dtype=dtype)
38    return kernel
39
40def load_dll(modelname, dtype='single'):
41    sasmodels = __import__('sasmodels.models.'+modelname)
42    module = getattr(sasmodels.models, modelname, None)
43    kernel = dll_model(module, dtype=dtype)
44    return kernel
45
46
47def compare(Ncpu, cpuname, cpupars, Ngpu, gpuname, gpupars):
48
49    #from sasmodels.core import load_data
50    #data = load_data('December/DEC07098.DAT')
51    from sasmodels.core import empty_data1D
52    data = empty_data1D(np.logspace(-4, -1, 128))
53    #from sasmodels.core import empty_2D, set_beam_stop
54    #data = empty_data2D(np.linspace(-0.05, 0.05, 128))
55    #set_beam_stop(data, 0.004)
56    is2D = hasattr(data, 'qx_data')
57
58    if Ngpu > 0:
59        gpumodel = load_opencl(gpuname, dtype='single')
60        model = BumpsModel(data, gpumodel, **gpupars)
61        toc = tic()
62        for i in range(Ngpu):
63            #pars['scale'] = np.random.rand()
64            model.update()
65            gpu = model.theory()
66        gpu_time = toc()*1000./Ngpu
67        print "ocl t=%.1f ms, intensity=%.0f"%(gpu_time, sum(gpu[~np.isnan(gpu)]))
68        #print max(gpu), min(gpu)
69
70    if 0 and Ncpu > 0: # Hack to compare ctypes vs. opencl
71        dllmodel = load_dll(gpuname, dtype='double')
72        model = BumpsModel(data, dllmodel, **gpupars)
73        toc = tic()
74        for i in range(Ncpu):
75            model.update()
76            cpu = model.theory()
77        cpu_time = toc()*1000./Ncpu
78        print "dll t=%.1f ms"%cpu_time
79
80    elif 0: # Hack to check new vs old for GpuCylinder
81        from Models.code_cylinder_f import GpuCylinder as oldgpu
82        from sasmodel import SasModel
83        oldmodel = SasModel(data, oldgpu, dtype='single', **cpupars)
84        toc = tic()
85        for i in range(Ngpu):
86            oldmodel.update()
87            cpu = oldmodel.theory()
88        cpu_time = toc()*1000./Ngpu
89        print "old t=%.1f ms"%cpu_time
90
91    elif Ncpu > 0:
92        cpumodel = sasview_model(cpuname, **cpupars)
93        toc = tic()
94        if is2D:
95            for i in range(Ncpu):
96                cpu = cpumodel.evalDistribution([data.qx_data, data.qy_data])
97        else:
98            for i in range(Ncpu):
99                cpu = cpumodel.evalDistribution(data.x)
100        cpu_time = toc()*1000./Ncpu
101        print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[model.index]))
102
103    if Ngpu > 0 and Ncpu > 0:
104        print "gpu/cpu", max(abs(gpu/cpu)), "%.15g"%max(abs(gpu)), "%.15g"%max(abs(cpu))
105        #cpu *= max(gpu/cpu)
106        abserr = (gpu - cpu)
107        relerr = (gpu - cpu)/cpu
108        print "max(|ocl-omp|)", max(abs(abserr[model.index]))
109        print "max(|(ocl-omp)/ocl|)", max(abs(relerr[model.index]))
110    #return
111
112    import matplotlib.pyplot as plt
113    if Ncpu > 0:
114        if Ngpu > 0: plt.subplot(131)
115        plot_data(data, cpu, scale='log')
116        plt.title("omp t=%.1f ms"%cpu_time)
117    if Ngpu > 0:
118        if Ncpu > 0: plt.subplot(132)
119        plot_data(data, gpu, scale='log')
120        plt.title("ocl t=%.1f ms"%gpu_time)
121    if Ncpu > 0 and Ngpu > 0:
122        plt.subplot(133)
123        plot_data(data, 1e8*relerr, scale='linear')
124        plt.title("max rel err = %.3g"%max(abs(relerr)))
125        if is2D: plt.colorbar()
126    plt.show()
127
128def rename(pars, **names):
129    newpars = pars.copy()
130    for new,old in names.items():
131        for variant in ("", "_pd", "_pd_n", "_pd_nsigma"):
132            if old+variant in newpars:
133                newpars[new+variant] = pars[old+variant]
134                del newpars[old+variant]
135    return newpars
136
137def rescale_sld(pars, names):
138    newpars = pars.copy()
139    for p in names:
140        newpars[p] *= 1e6
141    return newpars
142
143# ===========================================================================
144#
145MODELS = {}
146def model(name):
147    def gather_function(fn):
148        MODELS[name] = fn
149        return fn
150    return gather_function
151
152USAGE="""
153usage: compare model [Nopencl] [Nsasview]
154
155Compare the speed and value for a model between the SasView original and the
156OpenCL rewrite.
157
158* Nopencl is the number of times to run the OpenCL model (default=5)
159
160* Nsasview is the number of times to run the Sasview model (default=1)
161
162* model is the name of the model to compare:
163
164    %s
165"""
166
167def main():
168    if len(sys.argv) == 1:
169        models = "\n    ".join("%-7s: %s"%(k,v.__name__.replace('_',' '))
170                               for k,v in sorted(MODELS.items()))
171        print(USAGE%models)
172        sys.exit(1)
173
174    cpuname, cpupars, gpuname, gpupars = MODELS[sys.argv[1]]()
175    Nopencl = int(sys.argv[2]) if len(sys.argv) > 2 else 5
176    Nsasview = int(sys.argv[3]) if len(sys.argv) > 3 else 1
177
178    compare(Nsasview, cpuname, cpupars, Nopencl, gpuname, gpupars)
179
180@model('cyl')
181def cylinder():
182    cpupars = dict(
183        scale=.003, background=.1,
184        sldCyl=.291e-6, sldSolv=5.77e-6,
185        radius=264.1, length=66.96,
186        cyl_theta=85, cyl_phi=0,
187        radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,
188        length_pd=0.1,length_pd_n=1, length_pd_nsigma=3,
189        cyl_theta_pd=45, cyl_theta_pd_n=50, cyl_theta_pd_nsigma=3,
190        cyl_phi_pd=0.1, cyl_phi_pd_n=5, cyl_phi_pd_nsigma=3,
191        )
192    cpuname = 'Cylinder'
193
194    gpupars = rename(cpupars, theta='cyl_theta', phi='cyl_phi', sld='sldCyl', solvent_sld='sldSolv')
195    gpupars = rescale_sld(gpupars, ['sld', 'solvent_sld'])
196    gpuname = 'cylinder'
197    return cpuname, cpupars, gpuname, gpupars
198
199
200@model('ell')
201def ellipse():
202    pars = dict(
203        scale=.027, background=4.9,
204        sldEll=.297e-6, sldSolv=5.773e-6,
205        radius_a=60, radius_b=180,
206        axis_theta=0, axis_phi=90,
207        radius_a_pd=0.1, radius_a_pd_n=10, radius_a_pd_nsigma=3,
208        radius_b_pd=0.1, radius_b_pd_n=10, radius_b_pd_nsigma=3,
209        axis_theta_pd=0.1, axis_theta_pd_n=6, axis_theta_pd_nsigma=3,
210        axis_phi_pd=0.1, axis_phi_pd_n=6, axis_phi_pd_nsigma=3,
211        )
212
213    from Models.code_ellipse import GpuEllipse as gpumodel
214    model = sasview_model('Ellipsoid', **pars)
215
216    pars = rename(pars, theta='axis_theta', phi='axis_phi', sld='sldEll', solvent_sld='sldSolv')
217    pars = rescale_sld(pars, ['sld', 'solvent_sld'])
218    return model, gpumodel, pars
219
220
221@model('cscyl')
222def core_shell_cylinder(N=1):
223    pars = dict(
224        scale= 1.77881e-06, background=223.827,
225        core_sld=1e-6, shell_sld=.291e-6, solvent_sld=7.105e-6,
226        radius=325, thickness=25, length=34.2709,
227        axis_theta=90, axis_phi=0,
228        radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3,
229        length_pd=0.1, length_pd_n=10, length_pd_nsigma=3,
230        thickness_pd=0.1, thickness_pd_n=5, thickness_pd_nsigma=3,
231        axis_theta_pd=15.8, axis_theta_pd_n=20, axis_theta_pd_nsigma=5,
232        axis_phi_pd=0.0008748, axis_phi_pd_n=5, axis_phi_pd_nsigma=3,
233        )
234
235    model = sasview_model('CoreShellCylinder', **pars)
236    from Models.code_coreshellcyl_f import GpuCoreShellCylinder as gpumodel
237
238    pars = rename(pars, theta='axis_theta', phi='axis_phi')
239    pars = rescale_sld(pars, ['core_sld', 'shell_sld', 'solvent_sld'])
240    return model, gpumodel, pars
241
242
243@model('ell3')
244def triaxial_ellipse(N=1):
245    pars = dict(
246        scale=0.08, background=5,
247        sldEll=7.105e-6, sldSolv=.291e-6,
248        axis_theta=0, axis_phi=0, axis_psi=0,
249        semi_axisA=15, semi_axisB=20, semi_axisC=500,
250        axis_theta_pd=20, axis_theta_pd_n=10, axis_theta_pd_nsigma=3,
251        axis_phi_pd=.1, axis_phi_pd_n=10, axis_phi_pd_nsigma=3,
252        axis_psi_pd=30, axis_psi_pd_n=5, axis_psi_pd_nsigma=3,
253        semi_axisA_pd=.1, semi_axisA_pd_n=5, semi_axisA_pd_nsigma=3,
254        semi_axisB_pd=.1, semi_axisB_pd_n=5, semi_axisB_pd_nsigma=3,
255        semi_axisC_pd=.1, semi_axisC_pd_n=5, semi_axisC_pd_nsigma=3,
256        )
257
258    model = sasview_model('TriaxialEllipsoid', **pars)
259    from Models.code_triaxialellipse import GpuTriEllipse as gpumodel
260
261    pars = rename(pars,
262                  theta='axis_theta', phi='axis_phi', psi='axis_psi',
263                  sld='sldEll', solvent_sld='sldSolv',
264                  radius_a='semi_axisA', radius_b='semi_axisB',
265                  radius_c='semi_axisC',
266                  )
267    pars = rescale_sld(pars, ['sld', 'solvent_sld'])
268    return model, gpumodel, pars
269
270@model('lam')
271def lamellar(N=1):
272    pars = dict(
273        scale=0.08, background=0.003,
274        sld_bi=5.38e-6,sld_sol=7.105e-6,
275        bi_thick=19.2946,
276        bi_thick_pd= 0.37765, bi_thick_pd_n=40, bi_thick_pd_nsigma=3,
277        )
278
279    model = sasview_model('Lamellar', **pars)
280    from Models.code_lamellar import GpuLamellar as gpumodel
281
282    pars = rename(pars, sld='sld_bi', solvent_sld='sld_sol', thickness='bi_thick')
283    pars = rescale_sld(pars, ['sld', 'solvent_sld'])
284    return model, gpumodel, pars
285
286@model('capcyl')
287def capped_cylinder(N=1):
288    pars = dict(
289        scale=.08, background=0,
290        sld_capcyl=1e-6, sld_solv=6.3e-6,
291        rad_cyl=20, rad_cap=40, len_cyl=400,
292        theta=0, phi=0,
293        rad_cyl_pd=.1, rad_cyl_pd_n=10, rad_cyl_pd_nsigma=3,
294        rad_cap_pd=.1, rad_cap_pd_n=10, rad_cap_pd_nsigma=3,
295        len_cyl_pd=.1, len_cyl_pd_n=3, len_cyl_pd_nsigma=3,
296        theta_pd=.1, theta_pd_n=3, theta_pd_nsigma=3,
297        phi_pd=.1, phi_pd_n=3, phi_pd_nsigma=3,
298        )
299
300
301    model = sasview_model('CappedCylinder', **pars)
302    from Models.code_capcyl import GpuCapCylinder as gpumodel
303
304    pars = rename(pars,
305                  sld='sld_capcyl', solvent_sld='sld_solv',
306                  length='len_cyl', radius='rad_cyl',
307                  cap_radius='rad_cap')
308    pars = rescale_sld(pars, ['sld', 'solvent_sld'])
309    return model, gpumodel, pars
310
311
312if __name__ == "__main__":
313    main()
Note: See TracBrowser for help on using the repository browser.