source: sasmodels/sasmodels/compare.py @ 9a66e65

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

refactor sasmodels to sasview parameter conversion

  • Property mode set to 100755
File size: 17.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
8import datetime
9import traceback
10
11import numpy as np
12
13ROOT = dirname(__file__)
14sys.path.insert(0, ROOT)  # Make sure sasmodels is first on the path
15
16
17from . import core
18from . import kerneldll
19from . import generate
20from .data import plot_theory, empty_data1D, empty_data2D
21from .direct_model import DirectModel
22from .convert import revert_model, constrain_new_to_old
23kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True
24
25# List of available models
26MODELS = [basename(f)[:-3]
27          for f in sorted(glob.glob(joinpath(ROOT,"models","[a-zA-Z]*.py")))]
28
29# CRUFT python 2.6
30if not hasattr(datetime.timedelta, 'total_seconds'):
31    def delay(dt):
32        """Return number date-time delta as number seconds"""
33        return dt.days * 86400 + dt.seconds + 1e-6 * dt.microseconds
34else:
35    def delay(dt):
36        """Return number date-time delta as number seconds"""
37        return dt.total_seconds()
38
39
40def tic():
41    """
42    Timer function.
43
44    Use "toc=tic()" to start the clock and "toc()" to measure
45    a time interval.
46    """
47    then = datetime.datetime.now()
48    return lambda: delay(datetime.datetime.now() - then)
49
50
51def set_beam_stop(data, radius, outer=None):
52    """
53    Add a beam stop of the given *radius*.  If *outer*, make an annulus.
54
55    Note: this function does not use the sasview package
56    """
57    if hasattr(data, 'qx_data'):
58        q = np.sqrt(data.qx_data**2 + data.qy_data**2)
59        data.mask = (q < radius)
60        if outer is not None:
61            data.mask |= (q >= outer)
62    else:
63        data.mask = (data.x < radius)
64        if outer is not None:
65            data.mask |= (data.x >= outer)
66
67
68def sasview_model(model_definition, **pars):
69    """
70    Load a sasview model given the model name.
71    """
72    # convert model parameters from sasmodel form to sasview form
73    #print("old",sorted(pars.items()))
74    modelname, pars = revert_model(model_definition, pars)
75    #print("new",sorted(pars.items()))
76    sas = __import__('sas.models.'+modelname)
77    ModelClass = getattr(getattr(sas.models,modelname,None),modelname,None)
78    if ModelClass is None:
79        raise ValueError("could not find model %r in sas.models"%modelname)
80    model = ModelClass()
81
82    for k,v in pars.items():
83        parts = k.split('.')  # polydispersity components
84        if len(parts) == 2:
85            model.dispersion[parts[0]][parts[1]] = v
86        else:
87            model.setParam(k, v)
88    return model
89
90def randomize(p, v):
91    """
92    Randomizing parameter.
93
94    Guess the parameter type from name.
95    """
96    if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')):
97        return v
98    elif any(s in p for s in ('theta','phi','psi')):
99        # orientation in [-180,180], orientation pd in [0,45]
100        if p.endswith('_pd'):
101            return 45*np.random.rand()
102        else:
103            return 360*np.random.rand() - 180
104    elif 'sld' in p:
105        # sld in in [-0.5,10]
106        return 10.5*np.random.rand() - 0.5
107    elif p.endswith('_pd'):
108        # length pd in [0,1]
109        return np.random.rand()
110    else:
111        # values from 0 to 2*x for all other parameters
112        return 2*np.random.rand()*(v if v != 0 else 1)
113
114def randomize_model(pars, seed=None):
115    if seed is None:
116        seed = np.random.randint(1e9)
117    np.random.seed(seed)
118    # Note: the sort guarantees order of calls to random number generator
119    pars = dict((p,randomize(p,v)) for p,v in sorted(pars.items()))
120
121    return pars, seed
122
123def constrain_pars(model_definition, pars):
124    """
125    Restrict parameters to valid values.
126    """
127    name = model_definition.name
128    if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']:
129        pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius']
130    if name == 'barbell' and pars['bell_radius'] < pars['radius']:
131        pars['radius'],pars['bell_radius'] = pars['bell_radius'],pars['radius']
132
133    # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff)
134    if name == 'guinier':
135        #q_max = 0.2  # mid q maximum
136        q_max = 1.0  # high q maximum
137        rg_max = np.sqrt(90*np.log(10) + 3*np.log(pars['scale']))/q_max
138        pars['rg'] = min(pars['rg'],rg_max)
139
140def parlist(pars):
141    return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items()))
142
143def suppress_pd(pars):
144    """
145    Suppress theta_pd for now until the normalization is resolved.
146
147    May also suppress complete polydispersity of the model to test
148    models more quickly.
149    """
150    for p in pars:
151        if p.endswith("_pd"): pars[p] = 0
152
153def eval_sasview(model_definition, pars, data, Nevals=1):
154    # importing sas here so that the error message will be that sas failed to
155    # import rather than the more obscure smear_selection not imported error
156    import sas
157    from sas.models.qsmearing import smear_selection
158    model = sasview_model(model_definition, **pars)
159    smearer = smear_selection(data, model=model)
160    value = None  # silence the linter
161    toc = tic()
162    for _ in range(max(Nevals, 1)):  # make sure there is at least one eval
163        if hasattr(data, 'qx_data'):
164            q = np.sqrt(data.qx_data**2 + data.qy_data**2)
165            index = ((~data.mask) & (~np.isnan(data.data))
166                     & (q >= data.qmin) & (q <= data.qmax))
167            if smearer is not None:
168                smearer.model = model  # because smear_selection has a bug
169                smearer.accuracy = data.accuracy
170                smearer.set_index(index)
171                value = smearer.get_value()
172            else:
173                value = model.evalDistribution([data.qx_data[index], data.qy_data[index]])
174        else:
175            value = model.evalDistribution(data.x)
176            if smearer is not None:
177                value = smearer(value)
178    average_time = toc()*1000./Nevals
179    return value, average_time
180
181def eval_opencl(model_definition, pars, data, dtype='single', Nevals=1,
182                cutoff=0., fast=False):
183    try:
184        model = core.load_model(model_definition, dtype=dtype,
185                                platform="ocl", fast=fast)
186    except Exception as exc:
187        print(exc)
188        print("... trying again with single precision")
189        model = core.load_model(model_definition, dtype='single',
190                                platform="ocl", fast=fast)
191    calculator = DirectModel(data, model, cutoff=cutoff)
192    value = None  # silence the linter
193    toc = tic()
194    for _ in range(max(Nevals, 1)):  # force at least one eval
195        value = calculator(**pars)
196    average_time = toc()*1000./Nevals
197    return value, average_time
198
199
200def eval_ctypes(model_definition, pars, data, dtype='double', Nevals=1, cutoff=0.):
201    model = core.load_model(model_definition, dtype=dtype, platform="dll")
202    calculator = DirectModel(data, model, cutoff=cutoff)
203    value = None  # silence the linter
204    toc = tic()
205    for _ in range(max(Nevals, 1)):  # force at least one eval
206        value = calculator(**pars)
207    average_time = toc()*1000./Nevals
208    return value, average_time
209
210
211def make_data(qmax, is2D, Nq=128, resolution=0.0, accuracy='Low', view='log'):
212    if is2D:
213        data = empty_data2D(np.linspace(-qmax, qmax, Nq), resolution=resolution)
214        data.accuracy = accuracy
215        set_beam_stop(data, 0.004)
216        index = ~data.mask
217    else:
218        if view == 'log':
219            qmax = math.log10(qmax)
220            q = np.logspace(qmax-3, qmax, Nq)
221        else:
222            q = np.linspace(0.001*qmax, qmax, Nq)
223        data = empty_data1D(q, resolution=resolution)
224        index = slice(None, None)
225    return data, index
226
227def compare(name, pars, Ncomp, Nbase, opts, set_pars):
228    model_definition = core.load_model_definition(name)
229
230    view = ('linear' if '-linear' in opts
231            else 'log' if '-log' in opts
232            else 'q4' if '-q4' in opts
233            else 'log')
234
235    opt_values = dict(split
236                      for s in opts for split in ((s.split('='),))
237                      if len(split) == 2)
238    # Sort out data
239    qmax = (10.0 if '-exq' in opts
240            else 1.0 if '-highq' in opts
241            else 0.2 if '-midq' in opts
242            else 0.05)
243    Nq = int(opt_values.get('-Nq', '128'))
244    res = float(opt_values.get('-res', '0'))
245    accuracy = opt_values.get('-accuracy', 'Low')
246    is2D = "-2d" in opts
247    data, index = make_data(qmax, is2D, Nq, res, accuracy, view=view)
248
249
250    # modelling accuracy is determined by dtype and cutoff
251    dtype = ('longdouble' if '-quad' in opts
252             else 'double' if '-double' in opts
253             else 'half' if '-half' in opts
254             else 'single')
255    cutoff = float(opt_values.get('-cutoff','1e-5'))
256    fast = "-fast" in opts and dtype is 'single'
257
258    # randomize parameters
259    #pars.update(set_pars)  # set value before random to control range
260    if '-random' in opts or '-random' in opt_values:
261        seed = int(opt_values['-random']) if '-random' in opt_values else None
262        pars, seed = randomize_model(pars, seed=seed)
263        print("Randomize using -random=%i"%seed)
264    pars.update(set_pars)  # set value after random to control value
265    constrain_pars(model_definition, pars)
266    constrain_new_to_old(model_definition, pars)
267
268    # parameter selection
269    if '-mono' in opts:
270        suppress_pd(pars)
271    if '-pars' in opts:
272        print("pars "+str(parlist(pars)))
273
274    # Base calculation
275    if 0:
276        from sasmodels.models import sphere as target
277        base_name = target.name
278        base, base_time = eval_ctypes(target, pars, data,
279                dtype='longdouble', cutoff=0., Nevals=Ncomp)
280    elif Nbase > 0 and "-ctypes" in opts and "-sasview" in opts:
281        try:
282            base, base_time = eval_sasview(model_definition, pars, data, Ncomp)
283            base_name = "sasview"
284            #print("base/sasview", (base-pars['background'])/(comp-pars['background']))
285            print("sasview t=%.1f ms, intensity=%.0f"%(base_time, sum(base)))
286            #print("sasview",comp)
287        except ImportError:
288            traceback.print_exc()
289            Nbase = 0
290    elif Nbase > 0:
291        base, base_time = eval_opencl(model_definition, pars, data,
292                dtype=dtype, cutoff=cutoff, Nevals=Nbase, fast=fast)
293        base_name = "ocl"
294        print("opencl t=%.1f ms, intensity=%.0f"%(base_time, sum(base)))
295        #print("base " + base)
296        #print(max(base), min(base))
297
298    # Comparison calculation
299    if Ncomp > 0 and "-ctypes" in opts:
300        comp, comp_time = eval_ctypes(model_definition, pars, data,
301                dtype=dtype, cutoff=cutoff, Nevals=Ncomp)
302        comp_name = "ctypes"
303        print("ctypes t=%.1f ms, intensity=%.0f"%(comp_time, sum(comp)))
304    elif Ncomp > 0:
305        try:
306            comp, comp_time = eval_sasview(model_definition, pars, data, Ncomp)
307            comp_name = "sasview"
308            #print("base/sasview", (base-pars['background'])/(comp-pars['background']))
309            print("sasview t=%.1f ms, intensity=%.0f"%(comp_time, sum(comp)))
310            #print("sasview",comp)
311        except ImportError:
312            traceback.print_exc()
313            Ncomp = 0
314
315    # Compare, but only if computing both forms
316    if Nbase > 0 and Ncomp > 0:
317        #print("speedup %.2g"%(comp_time/base_time))
318        #print("max |base/comp|", max(abs(base/comp)), "%.15g"%max(abs(base)), "%.15g"%max(abs(comp)))
319        #comp *= max(base/comp)
320        resid = (base - comp)
321        relerr = resid/comp
322        #bad = (relerr>1e-4)
323        #print(relerr[bad],comp[bad],base[bad],data.qx_data[bad],data.qy_data[bad])
324        _print_stats("|%s-%s|"%(base_name,comp_name)+(" "*(3+len(comp_name))), resid)
325        _print_stats("|(%s-%s)/%s|"%(base_name,comp_name,comp_name), relerr)
326
327    # Plot if requested
328    if '-noplot' in opts: return
329    import matplotlib.pyplot as plt
330    if Ncomp > 0:
331        if Nbase > 0: plt.subplot(131)
332        plot_theory(data, comp, view=view, plot_data=False)
333        plt.title("%s t=%.1f ms"%(comp_name,comp_time))
334        #cbar_title = "log I"
335    if Nbase > 0:
336        if Ncomp > 0: plt.subplot(132)
337        plot_theory(data, base, view=view, plot_data=False)
338        plt.title("%s t=%.1f ms"%(base_name,base_time))
339        #cbar_title = "log I"
340    if Ncomp > 0 and Nbase > 0:
341        plt.subplot(133)
342        if '-abs' in opts:
343            err,errstr,errview = resid, "abs err", "linear"
344        else:
345            err,errstr,errview = abs(relerr), "rel err", "log"
346        #err,errstr = base/comp,"ratio"
347        plot_theory(data, None, resid=err, view=errview, plot_data=False)
348        plt.title("max %s = %.3g"%(errstr, max(abs(err))))
349        #cbar_title = errstr if errview=="linear" else "log "+errstr
350    #if is2D:
351    #    h = plt.colorbar()
352    #    h.ax.set_title(cbar_title)
353
354    if Ncomp > 0 and Nbase > 0 and '-hist' in opts:
355        plt.figure()
356        v = relerr
357        v[v==0] = 0.5*np.min(np.abs(v[v!=0]))
358        plt.hist(np.log10(np.abs(v)), normed=1, bins=50);
359        plt.xlabel('log10(err), err = | F(q) single - F(q) double| / | F(q) double |');
360        plt.ylabel('P(err)')
361        plt.title('Comparison of single and double precision models for %s'%name)
362
363    plt.show()
364
365def _print_stats(label, err):
366    sorted_err = np.sort(abs(err))
367    p50 = int((len(err)-1)*0.50)
368    p98 = int((len(err)-1)*0.98)
369    data = [
370        "max:%.3e"%sorted_err[-1],
371        "median:%.3e"%sorted_err[p50],
372        "98%%:%.3e"%sorted_err[p98],
373        "rms:%.3e"%np.sqrt(np.mean(err**2)),
374        "zero-offset:%+.3e"%np.mean(err),
375        ]
376    print(label+"  ".join(data))
377
378
379
380# ===========================================================================
381#
382USAGE="""
383usage: compare.py model [Nopencl] [Nsasview] [options...] [key=val]
384
385Compare the speed and value for a model between the SasView original and the
386OpenCL rewrite.
387
388model is the name of the model to compare (see below).
389Nopencl is the number of times to run the OpenCL model (default=5)
390Nsasview is the number of times to run the Sasview model (default=1)
391
392Options (* for default):
393
394    -plot*/-noplot plots or suppress the plot of the model
395    -half/-single*/-double/-quad/-fast sets the calculation precision
396    -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0
397    -Nq=128 sets the number of Q points in the data set
398    -1d*/-2d computes 1d or 2d data
399    -preset*/-random[=seed] preset or random parameters
400    -mono/-poly* force monodisperse/polydisperse
401    -ctypes/-sasview* selects gpu:cpu, gpu:sasview, or sasview:cpu if both
402    -cutoff=1e-5* cutoff value for including a point in polydispersity
403    -pars/-nopars* prints the parameter set or not
404    -abs/-rel* plot relative or absolute error
405    -linear/-log/-q4 intensity scaling
406    -hist/-nohist* plot histogram of relative error
407    -res=0 sets the resolution width dQ/Q if calculating with resolution
408    -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh
409
410Key=value pairs allow you to set specific values to any of the model
411parameters.
412
413Available models:
414"""
415
416
417NAME_OPTIONS = set([
418    'plot', 'noplot',
419    'half', 'single', 'double', 'quad', 'fast',
420    'lowq', 'midq', 'highq', 'exq',
421    '2d', '1d',
422    'preset', 'random',
423    'poly', 'mono',
424    'sasview', 'ctypes',
425    'nopars', 'pars',
426    'rel', 'abs',
427    'linear', 'log', 'q4',
428    'hist', 'nohist',
429    ])
430VALUE_OPTIONS = [
431    # Note: random is both a name option and a value option
432    'cutoff', 'random', 'Nq', 'res', 'accuracy',
433    ]
434
435def columnize(L, indent="", width=79):
436    column_width = max(len(w) for w in L) + 1
437    num_columns = (width - len(indent)) // column_width
438    num_rows = len(L) // num_columns
439    L = L + [""] * (num_rows*num_columns - len(L))
440    columns = [L[k*num_rows:(k+1)*num_rows] for k in range(num_columns)]
441    lines = [" ".join("%-*s"%(column_width, entry) for entry in row)
442             for row in zip(*columns)]
443    output = indent + ("\n"+indent).join(lines)
444    return output
445
446
447def get_demo_pars(model_definition):
448    info = generate.make_info(model_definition)
449    pars = dict((p[0],p[2]) for p in info['parameters'])
450    pars.update(info['demo'])
451    return pars
452
453def main():
454    opts = [arg for arg in sys.argv[1:] if arg.startswith('-')]
455    popts = [arg for arg in sys.argv[1:] if not arg.startswith('-') and '=' in arg]
456    args = [arg for arg in sys.argv[1:] if not arg.startswith('-') and '=' not in arg]
457    models = "\n    ".join("%-15s"%v for v in MODELS)
458    if len(args) == 0:
459        print(USAGE)
460        print(columnize(MODELS, indent="  "))
461        sys.exit(1)
462    if args[0] not in MODELS:
463        print("Model %r not available. Use one of:\n    %s"%(args[0],models))
464        sys.exit(1)
465    if len(args) > 3:
466        print("expected parameters: model Nopencl Nsasview")
467
468    invalid = [o[1:] for o in opts
469               if o[1:] not in NAME_OPTIONS
470                  and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)]
471    if invalid:
472        print("Invalid options: %s"%(", ".join(invalid)))
473        sys.exit(1)
474
475    # Get demo parameters from model definition, or use default parameters
476    # if model does not define demo parameters
477    name = args[0]
478    model_definition = core.load_model_definition(name)
479    pars = get_demo_pars(model_definition)
480
481    Ncomp = int(args[1]) if len(args) > 1 else 5
482    Nbase = int(args[2]) if len(args) > 2 else 1
483
484    # Fill in default polydispersity parameters
485    pds = set(p.split('_pd')[0] for p in pars if p.endswith('_pd'))
486    for p in pds:
487        if p+"_pd_nsigma" not in pars: pars[p+"_pd_nsigma"] = 3
488        if p+"_pd_type" not in pars: pars[p+"_pd_type"] = "gaussian"
489
490    # Fill in parameters given on the command line
491    set_pars = {}
492    for arg in popts:
493        k,v = arg.split('=',1)
494        if k not in pars:
495            # extract base name without distribution
496            s = set(p.split('_pd')[0] for p in pars)
497            print("%r invalid; parameters are: %s"%(k,", ".join(sorted(s))))
498            sys.exit(1)
499        set_pars[k] = float(v) if not v.endswith('type') else v
500
501    compare(name, pars, Ncomp, Nbase, opts, set_pars)
502
503if __name__ == "__main__":
504    main()
Note: See TracBrowser for help on using the repository browser.