source: sasmodels/compare-new.py @ 14de349

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

add autogenerated polydispersity loops

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