source: sasmodels/compare.py @ a217f7d

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

allow test of dll using single precision

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