source: sasmodels/compare.py @ 31819c5

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

add relerr histogram to compare

  • Property mode set to 100755
File size: 13.4 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3
4import sys
5import math
6
7import numpy as np
8
9from sasmodels.bumps_model import BumpsModel, plot_data, tic
10from sasmodels import gpu, dll
11from sasmodels.convert import revert_model
12
13
14def sasview_model(modelname, **pars):
15    """
16    Load a sasview model given the model name.
17    """
18    # convert model parameters from sasmodel form to sasview form
19    #print "old",sorted(pars.items())
20    modelname, pars = revert_model(modelname, pars)
21    #print "new",sorted(pars.items())
22    sans = __import__('sans.models.'+modelname)
23    ModelClass = getattr(getattr(sans.models,modelname,None),modelname,None)
24    if ModelClass is None:
25        raise ValueError("could not find model %r in sans.models"%modelname)
26    model = ModelClass()
27
28    for k,v in pars.items():
29        if k.endswith("_pd"):
30            model.dispersion[k[:-3]]['width'] = v
31        elif k.endswith("_pd_n"):
32            model.dispersion[k[:-5]]['npts'] = v
33        elif k.endswith("_pd_nsigma"):
34            model.dispersion[k[:-10]]['nsigmas'] = v
35        elif k.endswith("_pd_type"):
36            model.dispersion[k[:-8]]['type'] = v
37        else:
38            model.setParam(k, v)
39    return model
40
41def load_opencl(modelname, dtype='single'):
42    sasmodels = __import__('sasmodels.models.'+modelname)
43    module = getattr(sasmodels.models, modelname, None)
44    kernel = gpu.load_model(module, dtype=dtype)
45    return kernel
46
47def load_ctypes(modelname, dtype='single'):
48    sasmodels = __import__('sasmodels.models.'+modelname)
49    module = getattr(sasmodels.models, modelname, None)
50    kernel = dll.load_model(module, dtype=dtype)
51    return kernel
52
53def randomize(p, v):
54    """
55    Randomizing parameter.
56
57    Guess the parameter type from name.
58    """
59    if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')):
60        return v
61    elif any(s in p for s in ('theta','phi','psi')):
62        # orientation in [-180,180], orientation pd in [0,45]
63        if p.endswith('_pd'):
64            return 45*np.random.rand()
65        else:
66            return 360*np.random.rand() - 180
67    elif 'sld' in p:
68        # sld in in [-0.5,10]
69        return 10.5*np.random.rand() - 0.5
70    elif p.endswith('_pd'):
71        # length pd in [0,1]
72        return np.random.rand()
73    else:
74        # length, scale, background in [0,200]
75        return 200*np.random.rand()
76
77def parlist(pars):
78    return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items()))
79
80def suppress_pd(pars):
81    """
82    Suppress theta_pd for now until the normalization is resolved.
83
84    May also suppress complete polydispersity of the model to test
85    models more quickly.
86    """
87    for p in pars:
88        if p.endswith("_pd"): pars[p] = 0
89
90def compare(name, pars, Ncpu, Ngpu, opts, set_pars):
91    # Sort out data
92    qmax = 1.0 if '-highq' in opts else (0.2 if '-midq' in opts else 0.05)
93    if "-1d" in opts:
94        from sasmodels.bumps_model import empty_data1D
95        qmax = math.log10(qmax)
96        data = empty_data1D(np.logspace(qmax-3, qmax, 128))
97        index = slice(None, None)
98    else:
99        from sasmodels.bumps_model import empty_data2D, set_beam_stop
100        data = empty_data2D(np.linspace(-qmax, qmax, 128))
101        set_beam_stop(data, 0.004)
102        index = ~data.mask
103    is2D = hasattr(data, 'qx_data')
104
105    # modelling accuracy is determined by dtype and cutoff
106    dtype = 'double' if '-double' in opts else 'single'
107    cutoff_opts = [s for s in opts if s.startswith('-cutoff')]
108    cutoff = float(cutoff_opts[0].split('=')[1]) if cutoff_opts else 1e-5
109
110    # randomize parameters
111    random_opts = [s for s in opts if s.startswith('-random')]
112    if random_opts:
113        if '=' in random_opts[0]:
114            seed = int(random_opts[0].split('=')[1])
115        else:
116            seed = int(np.random.rand()*10000)
117        np.random.seed(seed)
118        print "Randomize using -random=%i"%seed
119        # Note: the sort guarantees order of calls to random number generator
120        pars = dict((p,randomize(p,v)) for p,v in sorted(pars.items()))
121        # The capped cylinder model has a constraint on its parameters
122        if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']:
123            pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius']
124    pars.update(set_pars)
125
126    # parameter selection
127    if '-mono' in opts:
128        suppress_pd(pars)
129    if '-pars' in opts:
130        print "pars",parlist(pars)
131
132    # OpenCl calculation
133    if Ngpu > 0:
134        try:
135            model = load_opencl(name, dtype=dtype)
136        except Exception,exc:
137            print exc
138            print "... trying again with single precision"
139            model = load_opencl(name, dtype='single')
140        problem = BumpsModel(data, model, cutoff=cutoff, **pars)
141        toc = tic()
142        for _ in range(Ngpu):
143            #pars['scale'] = np.random.rand()
144            problem.update()
145            gpu = problem.theory()
146        gpu_time = toc()*1000./Ngpu
147        print "opencl t=%.1f ms, intensity=%.0f"%(gpu_time, sum(gpu[index]))
148        #print max(gpu), min(gpu)
149
150    # ctypes/sasview calculation
151    if Ncpu > 0 and "-ctypes" in opts:
152        model = load_ctypes(name, dtype=dtype)
153        problem = BumpsModel(data, model, cutoff=cutoff, **pars)
154        toc = tic()
155        for _ in range(Ncpu):
156            problem.update()
157            cpu = problem.theory()
158        cpu_time = toc()*1000./Ncpu
159        comp = "ctypes"
160        print "ctypes t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index]))
161    elif Ncpu > 0:
162        model = sasview_model(name, **pars)
163        toc = tic()
164        for _ in range(Ncpu):
165            if is2D:
166                cpu = model.evalDistribution([data.qx_data, data.qy_data])
167            else:
168                cpu = model.evalDistribution(data.x)
169        cpu_time = toc()*1000./Ncpu
170        comp = "sasview"
171        print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index]))
172
173    # Compare, but only if computing both forms
174    if Ngpu > 0 and Ncpu > 0:
175        #print "speedup %.2g"%(cpu_time/gpu_time)
176        #print "max |gpu/cpu|", max(abs(gpu/cpu)), "%.15g"%max(abs(gpu)), "%.15g"%max(abs(cpu))
177        #cpu *= max(gpu/cpu)
178        resid, relerr = np.zeros_like(gpu), np.zeros_like(gpu)
179        resid[index] = (gpu - cpu)[index]
180        relerr[index] = resid[index]/cpu[index]
181        #bad = (relerr>1e-4)
182        #print relerr[bad],cpu[bad],gpu[bad],data.qx_data[bad],data.qy_data[bad]
183        print "max(|ocl-%s|)"%comp, max(abs(resid[index]))
184        print "max(|(ocl-%s)/%s|)"%(comp,comp), max(abs(relerr[index]))
185        p98 = int(len(relerr[index])*0.98)
186        print "98%% (|(ocl-%s)/%s|) <"%(comp,comp), np.sort(abs(relerr[index]))[p98]
187
188
189    # Plot if requested
190    if '-noplot' in opts: return
191    import matplotlib.pyplot as plt
192    if Ncpu > 0:
193        if Ngpu > 0: plt.subplot(131)
194        plot_data(data, cpu, scale='log')
195        plt.title("%s t=%.1f ms"%(comp,cpu_time))
196    if Ngpu > 0:
197        if Ncpu > 0: plt.subplot(132)
198        plot_data(data, gpu, scale='log')
199        plt.title("opencl t=%.1f ms"%gpu_time)
200    if Ncpu > 0 and Ngpu > 0:
201        plt.subplot(133)
202        err = resid if '-abs' in opts else relerr
203        errstr = "abs err" if '-abs' in opts else "rel err"
204        #err,errstr = gpu/cpu,"ratio"
205        plot_data(data, err, scale='linear')
206        plt.title("max %s = %.3g"%(errstr, max(abs(err[index]))))
207    if is2D: plt.colorbar()
208
209    if Ncpu > 0 and Ngpu > 0 and '-hist' in opts:
210        plt.figure()
211        v = relerr[index]
212        v[v==0] = 0.5*np.min(np.abs(v[v!=0]))
213        plt.hist(np.log10(np.abs(v)), normed=1, bins=50);
214        plt.xlabel('log10(err), err = | F(q) single - F(q) double| / | F(q) double |');
215        plt.ylabel('P(err)')
216        plt.title('Comparison of single and double precision models for %s'%name)
217
218    plt.show()
219
220# ===========================================================================
221#
222USAGE="""
223usage: compare.py model [Nopencl] [Nsasview] [options...] [key=val]
224
225Compare the speed and value for a model between the SasView original and the
226OpenCL rewrite.
227
228model is the name of the model to compare (see below).
229Nopencl is the number of times to run the OpenCL model (default=5)
230Nsasview is the number of times to run the Sasview model (default=1)
231
232Options (* for default):
233
234    -plot*/-noplot plots or suppress the plot of the model
235    -single*/-double uses double precision for comparison
236    -lowq*/-midq/-highq use q values up to 0.05, 0.2 or 1.0
237    -1d/-2d* uses 1d or 2d random data
238    -preset*/-random[=seed] preset or random parameters
239    -mono/-poly* force monodisperse/polydisperse
240    -ctypes/-sasview* whether cpu is tested using sasview or ctypes
241    -cutoff=1e-5*/value cutoff for including a point in polydispersity
242    -pars/-nopars* prints the parameter set or not
243    -abs/-rel* plot relative or absolute error
244    -hist/-nohist* plot histogram of relative error
245
246Key=value pairs allow you to set specific values to any of the model
247parameters.
248
249Available models:
250
251    %s
252"""
253
254VALID_OPTIONS = [
255    'plot','noplot',
256    'single','double',
257    'lowq','midq','highq',
258    '2d','1d',
259    'preset','random',
260    'poly','mono',
261    'sasview','ctypes',
262    'nopars','pars',
263    'rel','abs',
264    'hist','nohist',
265    ]
266
267def main():
268    opts = [arg for arg in sys.argv[1:] if arg.startswith('-')]
269    args = [arg for arg in sys.argv[1:] if not arg.startswith('-')]
270    models = "\n    ".join("%-7s: %s"%(k,v.__name__.replace('_',' '))
271                           for k,v in sorted(MODELS.items()))
272    if len(args) == 0:
273        print(USAGE%models)
274        sys.exit(1)
275    if args[0] not in MODELS:
276        print "Model %r not available. Use one of:\n    %s"%(args[0],models)
277        sys.exit(1)
278
279    valid_opts = set(VALID_OPTIONS)
280    invalid = [o[1:] for o in opts
281               if o[1:] not in valid_opts
282                  and not o.startswith('-cutoff=')
283                  and not o.startswith('-random=')]
284    if invalid:
285        print "Invalid options: %s"%(", ".join(invalid))
286        sys.exit(1)
287
288    name, pars = MODELS[args[0]]()
289    Nopencl = int(args[1]) if len(args) > 1 else 5
290    Nsasview = int(args[2]) if len(args) > 2 else 1
291
292    # Fill in default polydispersity parameters
293    pds = set(p.split('_pd')[0] for p in pars if p.endswith('_pd'))
294    for p in pds:
295        if p+"_pd_nsigma" not in pars: pars[p+"_pd_nsigma"] = 3
296        if p+"_pd_type" not in pars: pars[p+"_pd_type"] = "gaussian"
297
298    # Fill in parameters given on the command line
299    set_pars = {}
300    for arg in args[3:]:
301        k,v = arg.split('=')
302        if k not in pars:
303            # extract base name without distribution
304            s = set(p.split('_pd')[0] for p in pars)
305            print "%r invalid; parameters are: %s"%(k,", ".join(sorted(s)))
306            sys.exit(1)
307        set_pars[k] = float(v) if not v.endswith('type') else v
308
309    compare(name, pars, Nsasview, Nopencl, opts, set_pars)
310
311# ===========================================================================
312#
313
314MODELS = {}
315def model(name):
316    def gather_function(fn):
317        MODELS[name] = fn
318        return fn
319    return gather_function
320
321
322@model('cyl')
323def cylinder():
324    pars = dict(
325        scale=1, background=0,
326        sld=6, solvent_sld=1,
327        #radius=5, length=20,
328        radius=260, length=290,
329        theta=30, phi=0,
330        radius_pd=.2, radius_pd_n=1,
331        length_pd=.2,length_pd_n=1,
332        theta_pd=15, theta_pd_n=45,
333        phi_pd=15, phi_pd_n=1,
334        )
335    return 'cylinder', pars
336
337@model('capcyl')
338def capped_cylinder():
339    pars = dict(
340        scale=1, background=0,
341        sld=6, solvent_sld=1,
342        radius=260, cap_radius=290, length=290,
343        theta=30, phi=15,
344        radius_pd=.2, radius_pd_n=1,
345        cap_radius_pd=.2, cap_radius_pd_n=1,
346        length_pd=.2, length_pd_n=1,
347        theta_pd=15, theta_pd_n=45,
348        phi_pd=15, phi_pd_n=1,
349        )
350    return 'capped_cylinder', pars
351
352
353@model('cscyl')
354def core_shell_cylinder():
355    pars = dict(
356        scale=1, background=0,
357        core_sld=6, shell_sld=8, solvent_sld=1,
358        radius=45, thickness=25, length=340,
359        theta=30, phi=15,
360        radius_pd=.2, radius_pd_n=1,
361        length_pd=.2, length_pd_n=1,
362        thickness_pd=.2, thickness_pd_n=1,
363        theta_pd=15, theta_pd_n=45,
364        phi_pd=15, phi_pd_n=1,
365        )
366    return 'core_shell_cylinder', pars
367
368
369@model('ell')
370def ellipsoid():
371    pars = dict(
372        scale=1, background=0,
373        sld=6, solvent_sld=1,
374        rpolar=50, requatorial=30,
375        theta=30, phi=15,
376        rpolar_pd=.2, rpolar_pd_n=1,
377        requatorial_pd=.2, requatorial_pd_n=1,
378        theta_pd=15, theta_pd_n=45,
379        phi_pd=15, phi_pd_n=1,
380        )
381    return 'ellipsoid', pars
382
383
384@model('ell3')
385def triaxial_ellipsoid():
386    pars = dict(
387        scale=1, background=0,
388        sld=6, solvent_sld=1,
389        theta=30, phi=15, psi=5,
390        req_minor=25, req_major=36, rpolar=50,
391        req_minor_pd=0, req_minor_pd_n=1,
392        req_major_pd=0, req_major_pd_n=1,
393        rpolar_pd=.2, rpolar_pd_n=1,
394        theta_pd=15, theta_pd_n=45,
395        phi_pd=15, phi_pd_n=1,
396        psi_pd=15, psi_pd_n=1,
397        )
398    return 'triaxial_ellipsoid', pars
399
400@model('sph')
401def sphere():
402    pars = dict(
403        scale=1, background=0,
404        sld=6, solvent_sld=1,
405        radius=120,
406        radius_pd=.2, radius_pd_n=45,
407        )
408    return 'sphere', pars
409
410@model('lam')
411def lamellar():
412    pars = dict(
413        scale=1, background=0,
414        sld=6, solvent_sld=1,
415        thickness=40,
416        thickness_pd= 0.2, thickness_pd_n=40,
417        )
418    return 'lamellar', pars
419
420if __name__ == "__main__":
421    main()
Note: See TracBrowser for help on using the repository browser.