source: sasmodels/sasmodels/compare.py @ 82c299f

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

rpa model
  • Property mode set to 100755
File size: 22.9 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 parameter_range(p, v):
69    """
70    Choose a parameter range based on parameter name and initial value.
71    """
72    if p.endswith('_pd_n'):
73        return [0, 100]
74    elif p.endswith('_pd_nsigma'):
75        return [0, 5]
76    elif p.endswith('_pd_type'):
77        return v
78    elif any(s in p for s in ('theta','phi','psi')):
79        # orientation in [-180,180], orientation pd in [0,45]
80        if p.endswith('_pd'):
81            return [0,45]
82        else:
83            return [-180, 180]
84    elif 'sld' in p:
85        return [-0.5, 10]
86    elif p.endswith('_pd'):
87        return [0, 1]
88    elif p in ['background', 'scale']:
89        return [0, 1e3]
90    else:
91        return [0, (2*v if v>0 else 1)]
92
93def _randomize_one(p, v):
94    """
95    Randomizing parameter.
96    """
97    if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')):
98        return v
99    else:
100        return np.random.uniform(*parameter_range(p, v))
101
102def randomize_pars(pars, seed=None):
103    np.random.seed(seed)
104    # Note: the sort guarantees order `of calls to random number generator
105    pars = dict((p,_randomize_one(p,v))
106                for p,v in sorted(pars.items()))
107    return pars
108
109def constrain_pars(model_definition, pars):
110    """
111    Restrict parameters to valid values.
112    """
113    name = model_definition.name
114    if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']:
115        pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius']
116    if name == 'barbell' and pars['bell_radius'] < pars['radius']:
117        pars['radius'],pars['bell_radius'] = pars['bell_radius'],pars['radius']
118
119    # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff)
120    if name == 'guinier':
121        #q_max = 0.2  # mid q maximum
122        q_max = 1.0  # high q maximum
123        rg_max = np.sqrt(90*np.log(10) + 3*np.log(pars['scale']))/q_max
124        pars['rg'] = min(pars['rg'],rg_max)
125
126    if name == 'rpa':
127        # Make sure phi sums to 1.0
128        if pars['case_num'] < 2:
129            pars['Phia'] = 0.
130            pars['Phib'] = 0.
131        elif pars['case_num'] < 5:
132            pars['Phia'] = 0.
133        total = sum(pars['Phi'+c] for c in 'abcd')
134        for c in 'abcd':
135            pars['Phi'+c] /= total
136
137def parlist(pars):
138    return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items()))
139
140def suppress_pd(pars):
141    """
142    Suppress theta_pd for now until the normalization is resolved.
143
144    May also suppress complete polydispersity of the model to test
145    models more quickly.
146    """
147    pars = pars.copy()
148    for p in pars:
149        if p.endswith("_pd_n"): pars[p] = 0
150    return pars
151
152def eval_sasview(model_definition, data):
153    # importing sas here so that the error message will be that sas failed to
154    # import rather than the more obscure smear_selection not imported error
155    import sas
156    from sas.models.qsmearing import smear_selection
157
158    # convert model parameters from sasmodel form to sasview form
159    #print("old",sorted(pars.items()))
160    modelname, pars = revert_model(model_definition, {})
161    #print("new",sorted(pars.items()))
162    sas = __import__('sas.models.'+modelname)
163    ModelClass = getattr(getattr(sas.models,modelname,None),modelname,None)
164    if ModelClass is None:
165        raise ValueError("could not find model %r in sas.models"%modelname)
166    model = ModelClass()
167    smearer = smear_selection(data, model=model)
168
169    if hasattr(data, 'qx_data'):
170        q = np.sqrt(data.qx_data**2 + data.qy_data**2)
171        index = ((~data.mask) & (~np.isnan(data.data))
172                 & (q >= data.qmin) & (q <= data.qmax))
173        if smearer is not None:
174            smearer.model = model  # because smear_selection has a bug
175            smearer.accuracy = data.accuracy
176            smearer.set_index(index)
177            theory = lambda: smearer.get_value()
178        else:
179            theory = lambda: model.evalDistribution([data.qx_data[index], data.qy_data[index]])
180    elif smearer is not None:
181        theory = lambda: smearer(model.evalDistribution(data.x))
182    else:
183        theory = lambda: model.evalDistribution(data.x)
184
185    def calculator(**pars):
186        # paying for parameter conversion each time to keep life simple, if not fast
187        _, pars = revert_model(model_definition, pars)
188        for k,v in pars.items():
189            parts = k.split('.')  # polydispersity components
190            if len(parts) == 2:
191                model.dispersion[parts[0]][parts[1]] = v
192            else:
193                model.setParam(k, v)
194        return theory()
195
196    calculator.engine = "sasview"
197    return calculator
198
199DTYPE_MAP = {
200    'half': '16',
201    'fast': 'fast',
202    'single': '32',
203    'double': '64',
204    'quad': '128',
205    'f16': '16',
206    'f32': '32',
207    'f64': '64',
208    'longdouble': '128',
209}
210def eval_opencl(model_definition, data, dtype='single', cutoff=0.):
211    try:
212        model = core.load_model(model_definition, dtype=dtype, platform="ocl")
213    except Exception as exc:
214        print(exc)
215        print("... trying again with single precision")
216        dtype = 'single'
217        model = core.load_model(model_definition, dtype=dtype, platform="ocl")
218    calculator = DirectModel(data, model, cutoff=cutoff)
219    calculator.engine = "OCL%s"%DTYPE_MAP[dtype]
220    return calculator
221
222def eval_ctypes(model_definition, data, dtype='double', cutoff=0.):
223    if dtype=='quad':
224        dtype = 'longdouble'
225    model = core.load_model(model_definition, dtype=dtype, platform="dll")
226    calculator = DirectModel(data, model, cutoff=cutoff)
227    calculator.engine = "OMP%s"%DTYPE_MAP[dtype]
228    return calculator
229
230def time_calculation(calculator, pars, Nevals=1):
231    # initialize the code so time is more accurate
232    value = calculator(**suppress_pd(pars))
233    toc = tic()
234    for _ in range(max(Nevals, 1)):  # make sure there is at least one eval
235        value = calculator(**pars)
236    average_time = toc()*1000./Nevals
237    return value, average_time
238
239def make_data(opts):
240    qmax, nq, res = opts['qmax'], opts['nq'], opts['res']
241    if opts['is2d']:
242        data = empty_data2D(np.linspace(-qmax, qmax, nq), resolution=res)
243        data.accuracy = opts['accuracy']
244        set_beam_stop(data, 0.004)
245        index = ~data.mask
246    else:
247        if opts['view'] == 'log':
248            qmax = math.log10(qmax)
249            q = np.logspace(qmax-3, qmax, nq)
250        else:
251            q = np.linspace(0.001*qmax, qmax, nq)
252        data = empty_data1D(q, resolution=res)
253        index = slice(None, None)
254    return data, index
255
256def make_engine(model_definition, data, dtype, cutoff):
257    if dtype == 'sasview':
258        return eval_sasview(model_definition, data)
259    elif dtype.endswith('!'):
260        return eval_ctypes(model_definition, data, dtype=dtype[:-1],
261                           cutoff=cutoff)
262    else:
263        return eval_opencl(model_definition, data, dtype=dtype,
264                           cutoff=cutoff)
265
266def compare(opts, limits=None):
267    Nbase, Ncomp = opts['N1'], opts['N2']
268    pars = opts['pars']
269    data = opts['data']
270
271    # Base calculation
272    if Nbase > 0:
273        base = opts['engines'][0]
274        try:
275            base_value, base_time = time_calculation(base, pars, Nbase)
276            print("%s t=%.1f ms, intensity=%.0f"%(base.engine, base_time, sum(base_value)))
277        except ImportError:
278            traceback.print_exc()
279            Nbase = 0
280
281    # Comparison calculation
282    if Ncomp > 0:
283        comp = opts['engines'][1]
284        try:
285            comp_value, comp_time = time_calculation(comp, pars, Ncomp)
286            print("%s t=%.1f ms, intensity=%.0f"%(comp.engine, comp_time, sum(comp_value)))
287        except ImportError:
288            traceback.print_exc()
289            Ncomp = 0
290
291    # Compare, but only if computing both forms
292    if Nbase > 0 and Ncomp > 0:
293        #print("speedup %.2g"%(comp_time/base_time))
294        #print("max |base/comp|", max(abs(base_value/comp_value)), "%.15g"%max(abs(base_value)), "%.15g"%max(abs(comp_value)))
295        #comp *= max(base_value/comp_value)
296        resid = (base_value - comp_value)
297        relerr = resid/comp_value
298        _print_stats("|%s - %s|"%(base.engine,comp.engine)+(" "*(3+len(comp.engine))), resid)
299        _print_stats("|(%s - %s) / %s|"%(base.engine,comp.engine,comp.engine), relerr)
300
301    # Plot if requested
302    if not opts['plot'] and not opts['explore']: return
303    view = opts['view']
304    import matplotlib.pyplot as plt
305    if limits is None:
306        vmin, vmax = np.Inf, -np.Inf
307        if Nbase > 0:
308            vmin = min(vmin, min(base_value))
309            vmax = max(vmax, max(base_value))
310        if Ncomp > 0:
311            vmin = min(vmin, min(comp_value))
312            vmax = max(vmax, max(comp_value))
313        limits = vmin, vmax
314
315    if Nbase > 0:
316        if Ncomp > 0: plt.subplot(131)
317        plot_theory(data, base_value, view=view, plot_data=False, limits=limits)
318        plt.title("%s t=%.1f ms"%(base.engine, base_time))
319        #cbar_title = "log I"
320    if Ncomp > 0:
321        if Nbase > 0: plt.subplot(132)
322        plot_theory(data, comp_value, view=view, plot_data=False, limits=limits)
323        plt.title("%s t=%.1f ms"%(comp.engine,comp_time))
324        #cbar_title = "log I"
325    if Ncomp > 0 and Nbase > 0:
326        plt.subplot(133)
327        if '-abs' in opts:
328            err,errstr,errview = resid, "abs err", "linear"
329        else:
330            err,errstr,errview = abs(relerr), "rel err", "log"
331        #err,errstr = base/comp,"ratio"
332        plot_theory(data, None, resid=err, view=errview, plot_data=False)
333        plt.title("max %s = %.3g"%(errstr, max(abs(err))))
334        #cbar_title = errstr if errview=="linear" else "log "+errstr
335    #if is2D:
336    #    h = plt.colorbar()
337    #    h.ax.set_title(cbar_title)
338
339    if Ncomp > 0 and Nbase > 0 and '-hist' in opts:
340        plt.figure()
341        v = relerr
342        v[v==0] = 0.5*np.min(np.abs(v[v!=0]))
343        plt.hist(np.log10(np.abs(v)), normed=1, bins=50);
344        plt.xlabel('log10(err), err = |(%s - %s) / %s|'%(base.engine, comp.engine, comp.engine));
345        plt.ylabel('P(err)')
346        plt.title('Distribution of relative error between calculation engines')
347
348    if not opts['explore']:
349        plt.show()
350
351    return limits
352
353def _print_stats(label, err):
354    sorted_err = np.sort(abs(err))
355    p50 = int((len(err)-1)*0.50)
356    p98 = int((len(err)-1)*0.98)
357    data = [
358        "max:%.3e"%sorted_err[-1],
359        "median:%.3e"%sorted_err[p50],
360        "98%%:%.3e"%sorted_err[p98],
361        "rms:%.3e"%np.sqrt(np.mean(err**2)),
362        "zero-offset:%+.3e"%np.mean(err),
363        ]
364    print(label+"  ".join(data))
365
366
367
368# ===========================================================================
369#
370USAGE="""
371usage: compare.py model N1 N2 [options...] [key=val]
372
373Compare the speed and value for a model between the SasView original and the
374sasmodels rewrite.
375
376model is the name of the model to compare (see below).
377N1 is the number of times to run sasmodels (default=1).
378N2 is the number times to run sasview (default=1).
379
380Options (* for default):
381
382    -plot*/-noplot plots or suppress the plot of the model
383    -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0
384    -nq=128 sets the number of Q points in the data set
385    -1d*/-2d computes 1d or 2d data
386    -preset*/-random[=seed] preset or random parameters
387    -mono/-poly* force monodisperse/polydisperse
388    -cutoff=1e-5* cutoff value for including a point in polydispersity
389    -pars/-nopars* prints the parameter set or not
390    -abs/-rel* plot relative or absolute error
391    -linear/-log*/-q4 intensity scaling
392    -hist/-nohist* plot histogram of relative error
393    -res=0 sets the resolution width dQ/Q if calculating with resolution
394    -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh
395    -edit starts the parameter explorer
396
397Any two calculation engines can be selected for comparison:
398
399    -single/-double/-half/-fast sets an OpenCL calculation engine
400    -single!/-double!/-quad! sets an OpenMP calculation engine
401    -sasview sets the sasview calculation engine
402
403The default is -single -sasview.  Note that the interpretation of quad
404precision depends on architecture, and may vary from 64-bit to 128-bit,
405with 80-bit floats being common (1e-19 precision).
406
407Key=value pairs allow you to set specific values for the model parameters.
408
409Available models:
410"""
411
412
413NAME_OPTIONS = set([
414    'plot', 'noplot',
415    'half', 'fast', 'single', 'double',
416    'single!', 'double!', 'quad!', 'sasview',
417    'lowq', 'midq', 'highq', 'exq',
418    '2d', '1d',
419    'preset', 'random',
420    'poly', 'mono',
421    'nopars', 'pars',
422    'rel', 'abs',
423    'linear', 'log', 'q4',
424    'hist', 'nohist',
425    'edit',
426    ])
427VALUE_OPTIONS = [
428    # Note: random is both a name option and a value option
429    'cutoff', 'random', 'nq', 'res', 'accuracy',
430    ]
431
432def columnize(L, indent="", width=79):
433    column_width = max(len(w) for w in L) + 1
434    num_columns = (width - len(indent)) // column_width
435    num_rows = len(L) // num_columns
436    L = L + [""] * (num_rows*num_columns - len(L))
437    columns = [L[k*num_rows:(k+1)*num_rows] for k in range(num_columns)]
438    lines = [" ".join("%-*s"%(column_width, entry) for entry in row)
439             for row in zip(*columns)]
440    output = indent + ("\n"+indent).join(lines)
441    return output
442
443
444def get_demo_pars(model_definition):
445    info = generate.make_info(model_definition)
446    # Get the default values for the parameters
447    pars = dict((p[0],p[2]) for p in info['parameters'])
448
449    # Fill in default values for the polydispersity parameters
450    for p in info['parameters']:
451        if p[4] in ('volume', 'orientation'):
452            pars[p[0]+'_pd'] = 0.0
453            pars[p[0]+'_pd_n'] = 0
454            pars[p[0]+'_pd_nsigma'] = 3.0
455            pars[p[0]+'_pd_type'] = "gaussian"
456
457    # Plug in values given in demo
458    pars.update(info['demo'])
459    return pars
460
461def parse_opts():
462    flags = [arg for arg in sys.argv[1:] if arg.startswith('-')]
463    values = [arg for arg in sys.argv[1:] if not arg.startswith('-') and '=' in arg]
464    args = [arg for arg in sys.argv[1:] if not arg.startswith('-') and '=' not in arg]
465    models = "\n    ".join("%-15s"%v for v in MODELS)
466    if len(args) == 0:
467        print(USAGE)
468        print(columnize(MODELS, indent="  "))
469        sys.exit(1)
470    if args[0] not in MODELS:
471        print("Model %r not available. Use one of:\n    %s"%(args[0],models))
472        sys.exit(1)
473    if len(args) > 3:
474        print("expected parameters: model Nopencl Nsasview")
475
476    invalid = [o[1:] for o in flags
477               if o[1:] not in NAME_OPTIONS
478                  and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)]
479    if invalid:
480        print("Invalid options: %s"%(", ".join(invalid)))
481        sys.exit(1)
482
483
484    # Interpret the flags
485    opts = {
486        'plot'      : True,
487        'view'      : 'log',
488        'is2d'      : False,
489        'qmax'      : 0.05,
490        'nq'        : 128,
491        'res'       : 0.0,
492        'accuracy'  : 'Low',
493        'cutoff'    : 1e-5,
494        'seed'      : -1,  # default to preset
495        'mono'      : False,
496        'show_pars' : False,
497        'show_hist' : False,
498        'rel_err'   : True,
499        'explore'   : False,
500    }
501    engines = []
502    for arg in flags:
503        if arg == '-noplot':    opts['plot'] = False
504        elif arg == '-plot':    opts['plot'] = True
505        elif arg == '-linear':  opts['view'] = 'linear'
506        elif arg == '-log':     opts['view'] = 'log'
507        elif arg == '-q4':      opts['view'] = 'q4'
508        elif arg == '-1d':      opts['is2d'] = False
509        elif arg == '-2d':      opts['is2d'] = True
510        elif arg == '-exq':     opts['qmax'] = 10.0
511        elif arg == '-highq':   opts['qmax'] = 1.0
512        elif arg == '-midq':    opts['qmax'] = 0.2
513        elif arg == '-loq':     opts['qmax'] = 0.05
514        elif arg.startswith('-nq='):       opts['nq'] = int(arg[4:])
515        elif arg.startswith('-res='):      opts['res'] = float(arg[5:])
516        elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:]
517        elif arg.startswith('-cutoff='):   opts['cutoff'] = float(arg[8:])
518        elif arg.startswith('-random='):   opts['seed'] = int(arg[8:])
519        elif arg == '-random':  opts['seed'] = np.random.randint(1e6)
520        elif arg == '-preset':  opts['seed'] = -1
521        elif arg == '-mono':    opts['mono'] = True
522        elif arg == '-poly':    opts['mono'] = False
523        elif arg == '-pars':    opts['show_pars'] = True
524        elif arg == '-nopars':  opts['show_pars'] = False
525        elif arg == '-hist':    opts['show_hist'] = True
526        elif arg == '-nohist':  opts['show_hist'] = False
527        elif arg == '-rel':     opts['rel_err'] = True
528        elif arg == '-abs':     opts['rel_err'] = False
529        elif arg == '-half':    engines.append(arg[1:])
530        elif arg == '-fast':    engines.append(arg[1:])
531        elif arg == '-single':  engines.append(arg[1:])
532        elif arg == '-double':  engines.append(arg[1:])
533        elif arg == '-single!': engines.append(arg[1:])
534        elif arg == '-double!': engines.append(arg[1:])
535        elif arg == '-quad!':   engines.append(arg[1:])
536        elif arg == '-sasview': engines.append(arg[1:])
537        elif arg == '-edit':    opts['explore'] = True
538
539    if len(engines) == 0:
540        engines.extend(['single','sasview'])
541    elif len(engines) == 1:
542        if engines[0][0] != 'sasview':
543            engines.append('sasview')
544        else:
545            engines.append('single')
546    elif len(engines) > 2:
547        del engines[2:]
548
549    name = args[0]
550    model_definition = core.load_model_definition(name)
551
552    N1 = int(args[1]) if len(args) > 1 else 1
553    N2 = int(args[2]) if len(args) > 2 else 1
554
555    # Get demo parameters from model definition, or use default parameters
556    # if model does not define demo parameters
557    pars = get_demo_pars(model_definition)
558
559    # Fill in parameters given on the command line
560    presets = {}
561    for arg in values:
562        k,v = arg.split('=',1)
563        if k not in pars:
564            # extract base name without polydispersity info
565            s = set(p.split('_pd')[0] for p in pars)
566            print("%r invalid; parameters are: %s"%(k,", ".join(sorted(s))))
567            sys.exit(1)
568        presets[k] = float(v) if not k.endswith('type') else v
569
570    # randomize parameters
571    #pars.update(set_pars)  # set value before random to control range
572    if opts['seed'] > -1:
573        pars = randomize_pars(pars, seed=opts['seed'])
574        print("Randomize using -random=%i"%opts['seed'])
575    if opts['mono']:
576        pars = suppress_pd(pars)
577    pars.update(presets)  # set value after random to control value
578    constrain_pars(model_definition, pars)
579    constrain_new_to_old(model_definition, pars)
580    if opts['show_pars']:
581        print("pars " + str(parlist(pars)))
582
583    # Create the computational engines
584    data, _index = make_data(opts)
585    if N1:
586        base = make_engine(model_definition, data, engines[0], opts['cutoff'])
587    else:
588        base = None
589    if N2:
590        comp = make_engine(model_definition, data, engines[1], opts['cutoff'])
591    else:
592        comp = None
593
594    # Remember it all
595    opts.update({
596        'name'      : name,
597        'def'       : model_definition,
598        'N1'        : N1,
599        'N2'        : N2,
600        'presets'   : presets,
601        'pars'      : pars,
602        'data'      : data,
603        'engines'   : [base, comp],
604    })
605
606    return opts
607
608def main():
609    opts = parse_opts()
610    if opts['explore']:
611        explore(opts)
612    else:
613        compare(opts)
614
615def explore(opts):
616    import wx
617    from bumps.names import FitProblem
618    from bumps.gui.app_frame import AppFrame
619
620    problem = FitProblem(Explore(opts))
621    isMac = "cocoa" in wx.version()
622    app = wx.App()
623    frame = AppFrame(parent=None, title="explore")
624    if not isMac: frame.Show()
625    frame.panel.set_model(model=problem)
626    frame.panel.Layout()
627    frame.panel.aui.Split(0, wx.TOP)
628    if isMac: frame.Show()
629    app.MainLoop()
630
631class Explore(object):
632    """
633    Return a bumps wrapper for a SAS model comparison.
634    """
635    def __init__(self, opts):
636        from bumps.cli import config_matplotlib
637        import bumps_model
638        config_matplotlib()
639        self.opts = opts
640        info = generate.make_info(opts['def'])
641        pars, pd_types = bumps_model.create_parameters(info, **opts['pars'])
642        if not opts['is2d']:
643            active = [base + ext
644                      for base in info['partype']['pd-1d']
645                      for ext in ['','_pd','_pd_n','_pd_nsigma']]
646            active.extend(info['partype']['fixed-1d'])
647            for k in active:
648                v = pars[k]
649                v.range(*parameter_range(k, v.value))
650        else:
651            for k, v in pars.items():
652                v.range(*parameter_range(k, v.value))
653
654        self.pars = pars
655        self.pd_types = pd_types
656        self.limits = None
657
658    def numpoints(self):
659        """
660        Return the number of points
661        """
662        return len(self.pars) + 1  # so dof is 1
663
664    def parameters(self):
665        """
666        Return a dictionary of parameters
667        """
668        return self.pars
669
670    def nllf(self):
671        return 0.  # No nllf
672
673    def plot(self, view='log'):
674        """
675        Plot the data and residuals.
676        """
677        pars = dict((k, v.value) for k,v in self.pars.items())
678        pars.update(self.pd_types)
679        self.opts['pars'] = pars
680        limits = compare(self.opts, limits=self.limits)
681        if self.limits is None:
682            vmin, vmax = limits
683            vmax = 1.3*vmax
684            vmin = vmax*1e-7
685            self.limits = vmin, vmax
686
687
688if __name__ == "__main__":
689    main()
Note: See TracBrowser for help on using the repository browser.