source: sasmodels/sasmodels/compare.py @ a557a99

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since a557a99 was a557a99, checked in by wojciech, 5 years ago

A few fixes to make code cleaner and add PYOPENCL_CTX control

  • Property mode set to 100755
File size: 35.2 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3"""
4Program to compare models using different compute engines.
5
6This program lets you compare results between OpenCL and DLL versions
7of the code and between precision (half, fast, single, double, quad),
8where fast precision is single precision using native functions for
9trig, etc., and may not be completely IEEE 754 compliant.  This lets
10make sure that the model calculations are stable, or if you need to
11tag the model as double precision only.
12
13Run using ./compare.sh (Linux, Mac) or compare.bat (Windows) in the
14sasmodels root to see the command line options.
15
16Note that there is no way within sasmodels to select between an
17OpenCL CPU device and a GPU device, but you can do so by setting the
18PYOPENCL_CTX environment variable ahead of time.  Start a python
19interpreter and enter::
20
21    import pyopencl as cl
22    cl.create_some_context()
23
24This will prompt you to select from the available OpenCL devices
25and tell you which string to use for the PYOPENCL_CTX variable.
26On Windows you will need to remove the quotes.
27"""
28
29from __future__ import print_function
30
31import sys
32import math
33import datetime
34import traceback
35
36import numpy as np  # type: ignore
37
38from . import core
39from . import kerneldll
40from .data import plot_theory, empty_data1D, empty_data2D
41from .direct_model import DirectModel
42from .convert import revert_name, revert_pars, constrain_new_to_old
43
44try:
45    from typing import Optional, Dict, Any, Callable, Tuple
46except:
47    pass
48else:
49    from .modelinfo import ModelInfo, Parameter, ParameterSet
50    from .data import Data
51    Calculator = Callable[[float], np.ndarray]
52
53USAGE = """
54usage: compare.py model N1 N2 [options...] [key=val]
55
56Compare the speed and value for a model between the SasView original and the
57sasmodels rewrite.
58
59model is the name of the model to compare (see below).
60N1 is the number of times to run sasmodels (default=1).
61N2 is the number times to run sasview (default=1).
62
63Options (* for default):
64
65    -plot*/-noplot plots or suppress the plot of the model
66    -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0
67    -nq=128 sets the number of Q points in the data set
68    -zero indicates that q=0 should be included
69    -1d*/-2d computes 1d or 2d data
70    -preset*/-random[=seed] preset or random parameters
71    -mono/-poly* force monodisperse/polydisperse
72    -cutoff=1e-5* cutoff value for including a point in polydispersity
73    -pars/-nopars* prints the parameter set or not
74    -abs/-rel* plot relative or absolute error
75    -linear/-log*/-q4 intensity scaling
76    -hist/-nohist* plot histogram of relative error
77    -res=0 sets the resolution width dQ/Q if calculating with resolution
78    -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh
79    -edit starts the parameter explorer
80    -default/-demo* use demo vs default parameters
81
82Any two calculation engines can be selected for comparison:
83
84    -single/-double/-half/-fast sets an OpenCL calculation engine
85    -single!/-double!/-quad! sets an OpenMP calculation engine
86    -sasview sets the sasview calculation engine
87
88The default is -single -sasview.  Note that the interpretation of quad
89precision depends on architecture, and may vary from 64-bit to 128-bit,
90with 80-bit floats being common (1e-19 precision).
91
92Key=value pairs allow you to set specific values for the model parameters.
93"""
94
95# Update docs with command line usage string.   This is separate from the usual
96# doc string so that we can display it at run time if there is an error.
97# lin
98__doc__ = (__doc__  # pylint: disable=redefined-builtin
99           + """
100Program description
101-------------------
102
103"""
104           + USAGE)
105
106kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True
107
108# CRUFT python 2.6
109if not hasattr(datetime.timedelta, 'total_seconds'):
110    def delay(dt):
111        """Return number date-time delta as number seconds"""
112        return dt.days * 86400 + dt.seconds + 1e-6 * dt.microseconds
113else:
114    def delay(dt):
115        """Return number date-time delta as number seconds"""
116        return dt.total_seconds()
117
118
119class push_seed(object):
120    """
121    Set the seed value for the random number generator.
122
123    When used in a with statement, the random number generator state is
124    restored after the with statement is complete.
125
126    :Parameters:
127
128    *seed* : int or array_like, optional
129        Seed for RandomState
130
131    :Example:
132
133    Seed can be used directly to set the seed::
134
135        >>> from numpy.random import randint
136        >>> push_seed(24)
137        <...push_seed object at...>
138        >>> print(randint(0,1000000,3))
139        [242082    899 211136]
140
141    Seed can also be used in a with statement, which sets the random
142    number generator state for the enclosed computations and restores
143    it to the previous state on completion::
144
145        >>> with push_seed(24):
146        ...    print(randint(0,1000000,3))
147        [242082    899 211136]
148
149    Using nested contexts, we can demonstrate that state is indeed
150    restored after the block completes::
151
152        >>> with push_seed(24):
153        ...    print(randint(0,1000000))
154        ...    with push_seed(24):
155        ...        print(randint(0,1000000,3))
156        ...    print(randint(0,1000000))
157        242082
158        [242082    899 211136]
159        899
160
161    The restore step is protected against exceptions in the block::
162
163        >>> with push_seed(24):
164        ...    print(randint(0,1000000))
165        ...    try:
166        ...        with push_seed(24):
167        ...            print(randint(0,1000000,3))
168        ...            raise Exception()
169        ...    except Exception:
170        ...        print("Exception raised")
171        ...    print(randint(0,1000000))
172        242082
173        [242082    899 211136]
174        Exception raised
175        899
176    """
177    def __init__(self, seed=None):
178        # type: (Optional[int]) -> None
179        self._state = np.random.get_state()
180        np.random.seed(seed)
181
182    def __enter__(self):
183        # type: () -> None
184        pass
185
186    def __exit__(self, exc_type, exc_value, traceback):
187        # type: (Any, BaseException, Any) -> None
188        # TODO: better typing for __exit__ method
189        np.random.set_state(self._state)
190
191def tic():
192    # type: () -> Callable[[], float]
193    """
194    Timer function.
195
196    Use "toc=tic()" to start the clock and "toc()" to measure
197    a time interval.
198    """
199    then = datetime.datetime.now()
200    return lambda: delay(datetime.datetime.now() - then)
201
202
203def set_beam_stop(data, radius, outer=None):
204    # type: (Data, float, float) -> None
205    """
206    Add a beam stop of the given *radius*.  If *outer*, make an annulus.
207
208    Note: this function does not require sasview
209    """
210    if hasattr(data, 'qx_data'):
211        q = np.sqrt(data.qx_data**2 + data.qy_data**2)
212        data.mask = (q < radius)
213        if outer is not None:
214            data.mask |= (q >= outer)
215    else:
216        data.mask = (data.x < radius)
217        if outer is not None:
218            data.mask |= (data.x >= outer)
219
220
221def parameter_range(p, v):
222    # type: (str, float) -> Tuple[float, float]
223    """
224    Choose a parameter range based on parameter name and initial value.
225    """
226    # process the polydispersity options
227    if p.endswith('_pd_n'):
228        return 0., 100.
229    elif p.endswith('_pd_nsigma'):
230        return 0., 5.
231    elif p.endswith('_pd_type'):
232        raise ValueError("Cannot return a range for a string value")
233    elif any(s in p for s in ('theta', 'phi', 'psi')):
234        # orientation in [-180,180], orientation pd in [0,45]
235        if p.endswith('_pd'):
236            return 0., 45.
237        else:
238            return -180., 180.
239    elif p.endswith('_pd'):
240        return 0., 1.
241    elif 'sld' in p:
242        return -0.5, 10.
243    elif p == 'background':
244        return 0., 10.
245    elif p == 'scale':
246        return 0., 1.e3
247    elif v < 0.:
248        return 2.*v, -2.*v
249    else:
250        return 0., (2.*v if v > 0. else 1.)
251
252
253def _randomize_one(model_info, p, v):
254    # type: (ModelInfo, str, float) -> float
255    # type: (ModelInfo, str, str) -> str
256    """
257    Randomize a single parameter.
258    """
259    if any(p.endswith(s) for s in ('_pd', '_pd_n', '_pd_nsigma', '_pd_type')):
260        return v
261
262    # Find the parameter definition
263    for par in model_info.parameters.call_parameters:
264        if par.name == p:
265            break
266    else:
267        raise ValueError("unknown parameter %r"%p)
268
269    if len(par.limits) > 2:  # choice list
270        return np.random.randint(len(par.limits))
271
272    limits = par.limits
273    if not np.isfinite(limits).all():
274        low, high = parameter_range(p, v)
275        limits = (max(limits[0], low), min(limits[1], high))
276
277    return np.random.uniform(*limits)
278
279
280def randomize_pars(model_info, pars, seed=None):
281    # type: (ModelInfo, ParameterSet, int) -> ParameterSet
282    """
283    Generate random values for all of the parameters.
284
285    Valid ranges for the random number generator are guessed from the name of
286    the parameter; this will not account for constraints such as cap radius
287    greater than cylinder radius in the capped_cylinder model, so
288    :func:`constrain_pars` needs to be called afterward..
289    """
290    with push_seed(seed):
291        # Note: the sort guarantees order `of calls to random number generator
292        random_pars = dict((p, _randomize_one(model_info, p, v))
293                           for p, v in sorted(pars.items()))
294    return random_pars
295
296def constrain_pars(model_info, pars):
297    # type: (ModelInfo, ParameterSet) -> None
298    """
299    Restrict parameters to valid values.
300
301    This includes model specific code for models such as capped_cylinder
302    which need to support within model constraints (cap radius more than
303    cylinder radius in this case).
304
305    Warning: this updates the *pars* dictionary in place.
306    """
307    name = model_info.id
308    # if it is a product model, then just look at the form factor since
309    # none of the structure factors need any constraints.
310    if '*' in name:
311        name = name.split('*')[0]
312
313    if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']:
314        pars['radius'], pars['cap_radius'] = pars['cap_radius'], pars['radius']
315    if name == 'barbell' and pars['bell_radius'] < pars['radius']:
316        pars['radius'], pars['bell_radius'] = pars['bell_radius'], pars['radius']
317
318    # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff)
319    if name == 'guinier':
320        #q_max = 0.2  # mid q maximum
321        q_max = 1.0  # high q maximum
322        rg_max = np.sqrt(90*np.log(10) + 3*np.log(pars['scale']))/q_max
323        pars['rg'] = min(pars['rg'], rg_max)
324
325    if name == 'rpa':
326        # Make sure phi sums to 1.0
327        if pars['case_num'] < 2:
328            pars['Phi1'] = 0.
329            pars['Phi2'] = 0.
330        elif pars['case_num'] < 5:
331            pars['Phi1'] = 0.
332        total = sum(pars['Phi'+c] for c in '1234')
333        for c in '1234':
334            pars['Phi'+c] /= total
335
336def parlist(model_info, pars, is2d):
337    # type: (ModelInfo, ParameterSet, bool) -> str
338    """
339    Format the parameter list for printing.
340    """
341    lines = []
342    parameters = model_info.parameters
343    for p in parameters.user_parameters(pars, is2d):
344        fields = dict(
345            value=pars.get(p.id, p.default),
346            pd=pars.get(p.id+"_pd", 0.),
347            n=int(pars.get(p.id+"_pd_n", 0)),
348            nsigma=pars.get(p.id+"_pd_nsgima", 3.),
349            pdtype=pars.get(p.id+"_pd_type", 'gaussian'),
350            relative_pd=p.relative_pd,
351        )
352        lines.append(_format_par(p.name, **fields))
353    return "\n".join(lines)
354
355    #return "\n".join("%s: %s"%(p, v) for p, v in sorted(pars.items()))
356
357def _format_par(name, value=0., pd=0., n=0, nsigma=3., pdtype='gaussian',
358                relative_pd=False):
359    # type: (str, float, float, int, float, str) -> str
360    line = "%s: %g"%(name, value)
361    if pd != 0.  and n != 0:
362        if relative_pd:
363            pd *= value
364        line += " +/- %g  (%d points in [-%g,%g] sigma %s)"\
365                % (pd, n, nsigma, nsigma, pdtype)
366    return line
367
368def suppress_pd(pars):
369    # type: (ParameterSet) -> ParameterSet
370    """
371    Suppress theta_pd for now until the normalization is resolved.
372
373    May also suppress complete polydispersity of the model to test
374    models more quickly.
375    """
376    pars = pars.copy()
377    for p in pars:
378        if p.endswith("_pd_n"): pars[p] = 0
379    return pars
380
381def eval_sasview(model_info, data):
382    # type: (Modelinfo, Data) -> Calculator
383    """
384    Return a model calculator using the pre-4.0 SasView models.
385    """
386    # importing sas here so that the error message will be that sas failed to
387    # import rather than the more obscure smear_selection not imported error
388    import sas
389    import sas.models
390    from sas.models.qsmearing import smear_selection
391    from sas.models.MultiplicationModel import MultiplicationModel
392
393    def get_model_class(name):
394        # type: (str) -> "sas.models.BaseComponent"
395        #print("new",sorted(_pars.items()))
396        __import__('sas.models.' + name)
397        ModelClass = getattr(getattr(sas.models, name, None), name, None)
398        if ModelClass is None:
399            raise ValueError("could not find model %r in sas.models"%name)
400        return ModelClass
401
402    # WARNING: ugly hack when handling model!
403    # Sasview models with multiplicity need to be created with the target
404    # multiplicity, so we cannot create the target model ahead of time for
405    # for multiplicity models.  Instead we store the model in a list and
406    # update the first element of that list with the new multiplicity model
407    # every time we evaluate.
408
409    # grab the sasview model, or create it if it is a product model
410    if model_info.composition:
411        composition_type, parts = model_info.composition
412        if composition_type == 'product':
413            P, S = [get_model_class(revert_name(p))() for p in parts]
414            model = [MultiplicationModel(P, S)]
415        else:
416            raise ValueError("sasview mixture models not supported by compare")
417    else:
418        old_name = revert_name(model_info)
419        if old_name is None:
420            raise ValueError("model %r does not exist in old sasview"
421                            % model_info.id)
422        ModelClass = get_model_class(old_name)
423        model = [ModelClass()]
424
425    # build a smearer with which to call the model, if necessary
426    smearer = smear_selection(data, model=model)
427    if hasattr(data, 'qx_data'):
428        q = np.sqrt(data.qx_data**2 + data.qy_data**2)
429        index = ((~data.mask) & (~np.isnan(data.data))
430                 & (q >= data.qmin) & (q <= data.qmax))
431        if smearer is not None:
432            smearer.model = model  # because smear_selection has a bug
433            smearer.accuracy = data.accuracy
434            smearer.set_index(index)
435            def _call_smearer():
436                smearer.model = model[0]
437                return smearer.get_value()
438            theory = _call_smearer
439        else:
440            theory = lambda: model[0].evalDistribution([data.qx_data[index],
441                                                        data.qy_data[index]])
442    elif smearer is not None:
443        theory = lambda: smearer(model[0].evalDistribution(data.x))
444    else:
445        theory = lambda: model[0].evalDistribution(data.x)
446
447    def calculator(**pars):
448        # type: (float, ...) -> np.ndarray
449        """
450        Sasview calculator for model.
451        """
452        oldpars = revert_pars(model_info, pars)
453        # For multiplicity models, create a model with the correct multiplicity
454        control = oldpars.pop("CONTROL", None)
455        if control is not None:
456            # sphericalSLD has one fewer multiplicity.  This update should
457            # happen in revert_pars, but it hasn't been called yet.
458            model[0] = ModelClass(control)
459        # paying for parameter conversion each time to keep life simple, if not fast
460        #print("sasview pars",oldpars)
461        for k, v in oldpars.items():
462            name_attr = k.split('.')  # polydispersity components
463            if len(name_attr) == 2:
464                model[0].dispersion[name_attr[0]][name_attr[1]] = v
465            else:
466                model[0].setParam(k, v)
467        return theory()
468
469    calculator.engine = "sasview"
470    return calculator
471
472DTYPE_MAP = {
473    'half': '16',
474    'fast': 'fast',
475    'single': '32',
476    'double': '64',
477    'quad': '128',
478    'f16': '16',
479    'f32': '32',
480    'f64': '64',
481    'longdouble': '128',
482}
483def eval_opencl(model_info, data, dtype='single', cutoff=0.):
484    # type: (ModelInfo, Data, str, float) -> Calculator
485    """
486    Return a model calculator using the OpenCL calculation engine.
487    """
488    if not core.HAVE_OPENCL:
489        raise RuntimeError("OpenCL not available")
490    model = core.build_model(model_info, dtype=dtype, platform="ocl")
491    calculator = DirectModel(data, model, cutoff=cutoff)
492    calculator.engine = "OCL%s"%DTYPE_MAP[dtype]
493    return calculator
494
495def eval_ctypes(model_info, data, dtype='double', cutoff=0.):
496    # type: (ModelInfo, Data, str, float) -> Calculator
497    """
498    Return a model calculator using the DLL calculation engine.
499    """
500    model = core.build_model(model_info, dtype=dtype, platform="dll")
501    calculator = DirectModel(data, model, cutoff=cutoff)
502    calculator.engine = "OMP%s"%DTYPE_MAP[dtype]
503    return calculator
504
505def time_calculation(calculator, pars, evals=1):
506    # type: (Calculator, ParameterSet, int) -> Tuple[np.ndarray, float]
507    """
508    Compute the average calculation time over N evaluations.
509
510    An additional call is generated without polydispersity in order to
511    initialize the calculation engine, and make the average more stable.
512    """
513    # initialize the code so time is more accurate
514    if evals > 1:
515        calculator(**suppress_pd(pars))
516    toc = tic()
517    # make sure there is at least one eval
518    value = calculator(**pars)
519    for _ in range(evals-1):
520        value = calculator(**pars)
521    average_time = toc()*1000. / evals
522    #print("I(q)",value)
523    return value, average_time
524
525def make_data(opts):
526    # type: (Dict[str, Any]) -> Tuple[Data, np.ndarray]
527    """
528    Generate an empty dataset, used with the model to set Q points
529    and resolution.
530
531    *opts* contains the options, with 'qmax', 'nq', 'res',
532    'accuracy', 'is2d' and 'view' parsed from the command line.
533    """
534    qmax, nq, res = opts['qmax'], opts['nq'], opts['res']
535    if opts['is2d']:
536        q = np.linspace(-qmax, qmax, nq)  # type: np.ndarray
537        data = empty_data2D(q, resolution=res)
538        data.accuracy = opts['accuracy']
539        set_beam_stop(data, 0.0004)
540        index = ~data.mask
541    else:
542        if opts['view'] == 'log' and not opts['zero']:
543            qmax = math.log10(qmax)
544            q = np.logspace(qmax-3, qmax, nq)
545        else:
546            q = np.linspace(0.001*qmax, qmax, nq)
547        if opts['zero']:
548            q = np.hstack((0, q))
549        data = empty_data1D(q, resolution=res)
550        index = slice(None, None)
551    return data, index
552
553def make_engine(model_info, data, dtype, cutoff):
554    # type: (ModelInfo, Data, str, float) -> Calculator
555    """
556    Generate the appropriate calculation engine for the given datatype.
557
558    Datatypes with '!' appended are evaluated using external C DLLs rather
559    than OpenCL.
560    """
561    if dtype == 'sasview':
562        return eval_sasview(model_info, data)
563    elif dtype.endswith('!'):
564        return eval_ctypes(model_info, data, dtype=dtype[:-1], cutoff=cutoff)
565    else:
566        return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff)
567
568def _show_invalid(data, theory):
569    # type: (Data, np.ma.ndarray) -> None
570    """
571    Display a list of the non-finite values in theory.
572    """
573    if not theory.mask.any():
574        return
575
576    if hasattr(data, 'x'):
577        bad = zip(data.x[theory.mask], theory[theory.mask])
578        print("   *** ", ", ".join("I(%g)=%g"%(x, y) for x, y in bad))
579
580
581def compare(opts, limits=None):
582    # type: (Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float]
583    """
584    Preform a comparison using options from the command line.
585
586    *limits* are the limits on the values to use, either to set the y-axis
587    for 1D or to set the colormap scale for 2D.  If None, then they are
588    inferred from the data and returned. When exploring using Bumps,
589    the limits are set when the model is initially called, and maintained
590    as the values are adjusted, making it easier to see the effects of the
591    parameters.
592    """
593    n_base, n_comp = opts['n1'], opts['n2']
594    pars = opts['pars']
595    data = opts['data']
596
597    # silence the linter
598    base = opts['engines'][0] if n_base else None
599    comp = opts['engines'][1] if n_comp else None
600    base_time = comp_time = None
601    base_value = comp_value = resid = relerr = None
602
603    # Base calculation
604    if n_base > 0:
605        try:
606            base_raw, base_time = time_calculation(base, pars, n_base)
607            base_value = np.ma.masked_invalid(base_raw)
608            print("%s t=%.2f ms, intensity=%.0f"
609                  % (base.engine, base_time, base_value.sum()))
610            _show_invalid(data, base_value)
611        except ImportError:
612            traceback.print_exc()
613            n_base = 0
614
615    # Comparison calculation
616    if n_comp > 0:
617        try:
618            comp_raw, comp_time = time_calculation(comp, pars, n_comp)
619            comp_value = np.ma.masked_invalid(comp_raw)
620            print("%s t=%.2f ms, intensity=%.0f"
621                  % (comp.engine, comp_time, comp_value.sum()))
622            _show_invalid(data, comp_value)
623        except ImportError:
624            traceback.print_exc()
625            n_comp = 0
626
627    # Compare, but only if computing both forms
628    if n_base > 0 and n_comp > 0:
629        resid = (base_value - comp_value)
630        relerr = resid/np.where(comp_value != 0., abs(comp_value), 1.0)
631        _print_stats("|%s-%s|"
632                     % (base.engine, comp.engine) + (" "*(3+len(comp.engine))),
633                     resid)
634        _print_stats("|(%s-%s)/%s|"
635                     % (base.engine, comp.engine, comp.engine),
636                     relerr)
637
638    # Plot if requested
639    if not opts['plot'] and not opts['explore']: return
640    view = opts['view']
641    import matplotlib.pyplot as plt
642    if limits is None:
643        vmin, vmax = np.Inf, -np.Inf
644        if n_base > 0:
645            vmin = min(vmin, base_value.min())
646            vmax = max(vmax, base_value.max())
647        if n_comp > 0:
648            vmin = min(vmin, comp_value.min())
649            vmax = max(vmax, comp_value.max())
650        limits = vmin, vmax
651
652    if n_base > 0:
653        if n_comp > 0: plt.subplot(131)
654        plot_theory(data, base_value, view=view, use_data=False, limits=limits)
655        plt.title("%s t=%.2f ms"%(base.engine, base_time))
656        #cbar_title = "log I"
657    if n_comp > 0:
658        if n_base > 0: plt.subplot(132)
659        plot_theory(data, comp_value, view=view, use_data=False, limits=limits)
660        plt.title("%s t=%.2f ms"%(comp.engine, comp_time))
661        #cbar_title = "log I"
662    if n_comp > 0 and n_base > 0:
663        plt.subplot(133)
664        if not opts['rel_err']:
665            err, errstr, errview = resid, "abs err", "linear"
666        else:
667            err, errstr, errview = abs(relerr), "rel err", "log"
668        #err,errstr = base/comp,"ratio"
669        plot_theory(data, None, resid=err, view=errview, use_data=False)
670        if view == 'linear':
671            plt.xscale('linear')
672        plt.title("max %s = %.3g"%(errstr, abs(err).max()))
673        #cbar_title = errstr if errview=="linear" else "log "+errstr
674    #if is2D:
675    #    h = plt.colorbar()
676    #    h.ax.set_title(cbar_title)
677    fig = plt.gcf()
678    fig.suptitle(opts['name'])
679
680    if n_comp > 0 and n_base > 0 and '-hist' in opts:
681        plt.figure()
682        v = relerr
683        v[v == 0] = 0.5*np.min(np.abs(v[v != 0]))
684        plt.hist(np.log10(np.abs(v)), normed=1, bins=50)
685        plt.xlabel('log10(err), err = |(%s - %s) / %s|'
686                   % (base.engine, comp.engine, comp.engine))
687        plt.ylabel('P(err)')
688        plt.title('Distribution of relative error between calculation engines')
689
690    if not opts['explore']:
691        plt.show()
692
693    return limits
694
695def _print_stats(label, err):
696    # type: (str, np.ma.ndarray) -> None
697    # work with trimmed data, not the full set
698    sorted_err = np.sort(abs(err.compressed()))
699    p50 = int((len(sorted_err)-1)*0.50)
700    p98 = int((len(sorted_err)-1)*0.98)
701    data = [
702        "max:%.3e"%sorted_err[-1],
703        "median:%.3e"%sorted_err[p50],
704        "98%%:%.3e"%sorted_err[p98],
705        "rms:%.3e"%np.sqrt(np.mean(sorted_err**2)),
706        "zero-offset:%+.3e"%np.mean(sorted_err),
707        ]
708    print(label+"  "+"  ".join(data))
709
710
711
712# ===========================================================================
713#
714NAME_OPTIONS = set([
715    'plot', 'noplot',
716    'half', 'fast', 'single', 'double',
717    'single!', 'double!', 'quad!', 'sasview',
718    'lowq', 'midq', 'highq', 'exq', 'zero',
719    '2d', '1d',
720    'preset', 'random',
721    'poly', 'mono',
722    'nopars', 'pars',
723    'rel', 'abs',
724    'linear', 'log', 'q4',
725    'hist', 'nohist',
726    'edit',
727    'demo', 'default',
728    ])
729VALUE_OPTIONS = [
730    # Note: random is both a name option and a value option
731    'cutoff', 'random', 'nq', 'res', 'accuracy',
732    ]
733
734def columnize(items, indent="", width=79):
735    # type: (List[str], str, int) -> str
736    """
737    Format a list of strings into columns.
738
739    Returns a string with carriage returns ready for printing.
740    """
741    column_width = max(len(w) for w in items) + 1
742    num_columns = (width - len(indent)) // column_width
743    num_rows = len(items) // num_columns
744    items = items + [""] * (num_rows * num_columns - len(items))
745    columns = [items[k*num_rows:(k+1)*num_rows] for k in range(num_columns)]
746    lines = [" ".join("%-*s"%(column_width, entry) for entry in row)
747             for row in zip(*columns)]
748    output = indent + ("\n"+indent).join(lines)
749    return output
750
751
752def get_pars(model_info, use_demo=False):
753    # type: (ModelInfo, bool) -> ParameterSet
754    """
755    Extract demo parameters from the model definition.
756    """
757    # Get the default values for the parameters
758    pars = {}
759    for p in model_info.parameters.call_parameters:
760        parts = [('', p.default)]
761        if p.polydisperse:
762            parts.append(('_pd', 0.0))
763            parts.append(('_pd_n', 0))
764            parts.append(('_pd_nsigma', 3.0))
765            parts.append(('_pd_type', "gaussian"))
766        for ext, val in parts:
767            if p.length > 1:
768                dict(("%s%d%s" % (p.id, k, ext), val)
769                     for k in range(1, p.length+1))
770            else:
771                pars[p.id + ext] = val
772
773    # Plug in values given in demo
774    if use_demo:
775        pars.update(model_info.demo)
776    return pars
777
778
779def parse_opts():
780    # type: () -> Dict[str, Any]
781    """
782    Parse command line options.
783    """
784    MODELS = core.list_models()
785    flags = [arg for arg in sys.argv[1:]
786             if arg.startswith('-')]
787    values = [arg for arg in sys.argv[1:]
788              if not arg.startswith('-') and '=' in arg]
789    args = [arg for arg in sys.argv[1:]
790            if not arg.startswith('-') and '=' not in arg]
791    models = "\n    ".join("%-15s"%v for v in MODELS)
792    if len(args) == 0:
793        print(USAGE)
794        print("\nAvailable models:")
795        print(columnize(MODELS, indent="  "))
796        sys.exit(1)
797    if len(args) > 3:
798        print("expected parameters: model N1 N2")
799
800    name = args[0]
801    try:
802        model_info = core.load_model_info(name)
803    except ImportError as exc:
804        print(str(exc))
805        print("Could not find model; use one of:\n    " + models)
806        sys.exit(1)
807
808    invalid = [o[1:] for o in flags
809               if o[1:] not in NAME_OPTIONS
810               and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)]
811    if invalid:
812        print("Invalid options: %s"%(", ".join(invalid)))
813        sys.exit(1)
814
815
816    # pylint: disable=bad-whitespace
817    # Interpret the flags
818    opts = {
819        'plot'      : True,
820        'view'      : 'log',
821        'is2d'      : False,
822        'qmax'      : 0.05,
823        'nq'        : 128,
824        'res'       : 0.0,
825        'accuracy'  : 'Low',
826        'cutoff'    : 0.0,
827        'seed'      : -1,  # default to preset
828        'mono'      : False,
829        'show_pars' : False,
830        'show_hist' : False,
831        'rel_err'   : True,
832        'explore'   : False,
833        'use_demo'  : True,
834        'zero'      : False,
835    }
836    engines = []
837    for arg in flags:
838        if arg == '-noplot':    opts['plot'] = False
839        elif arg == '-plot':    opts['plot'] = True
840        elif arg == '-linear':  opts['view'] = 'linear'
841        elif arg == '-log':     opts['view'] = 'log'
842        elif arg == '-q4':      opts['view'] = 'q4'
843        elif arg == '-1d':      opts['is2d'] = False
844        elif arg == '-2d':      opts['is2d'] = True
845        elif arg == '-exq':     opts['qmax'] = 10.0
846        elif arg == '-highq':   opts['qmax'] = 1.0
847        elif arg == '-midq':    opts['qmax'] = 0.2
848        elif arg == '-lowq':    opts['qmax'] = 0.05
849        elif arg == '-zero':    opts['zero'] = True
850        elif arg.startswith('-nq='):       opts['nq'] = int(arg[4:])
851        elif arg.startswith('-res='):      opts['res'] = float(arg[5:])
852        elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:]
853        elif arg.startswith('-cutoff='):   opts['cutoff'] = float(arg[8:])
854        elif arg.startswith('-random='):   opts['seed'] = int(arg[8:])
855        elif arg == '-random':  opts['seed'] = np.random.randint(1000000)
856        elif arg == '-preset':  opts['seed'] = -1
857        elif arg == '-mono':    opts['mono'] = True
858        elif arg == '-poly':    opts['mono'] = False
859        elif arg == '-pars':    opts['show_pars'] = True
860        elif arg == '-nopars':  opts['show_pars'] = False
861        elif arg == '-hist':    opts['show_hist'] = True
862        elif arg == '-nohist':  opts['show_hist'] = False
863        elif arg == '-rel':     opts['rel_err'] = True
864        elif arg == '-abs':     opts['rel_err'] = False
865        elif arg == '-half':    engines.append(arg[1:])
866        elif arg == '-fast':    engines.append(arg[1:])
867        elif arg == '-single':  engines.append(arg[1:])
868        elif arg == '-double':  engines.append(arg[1:])
869        elif arg == '-single!': engines.append(arg[1:])
870        elif arg == '-double!': engines.append(arg[1:])
871        elif arg == '-quad!':   engines.append(arg[1:])
872        elif arg == '-sasview': engines.append(arg[1:])
873        elif arg == '-edit':    opts['explore'] = True
874        elif arg == '-demo':    opts['use_demo'] = True
875        elif arg == '-default':    opts['use_demo'] = False
876    # pylint: enable=bad-whitespace
877
878    if len(engines) == 0:
879        engines.extend(['single', 'double'])
880    elif len(engines) == 1:
881        if engines[0][0] == 'double':
882            engines.append('single')
883        else:
884            engines.append('double')
885    elif len(engines) > 2:
886        del engines[2:]
887
888    n1 = int(args[1]) if len(args) > 1 else 1
889    n2 = int(args[2]) if len(args) > 2 else 1
890    use_sasview = any(engine == 'sasview' and count > 0
891                      for engine, count in zip(engines, [n1, n2]))
892
893    # Get demo parameters from model definition, or use default parameters
894    # if model does not define demo parameters
895    pars = get_pars(model_info, opts['use_demo'])
896
897
898    # Fill in parameters given on the command line
899    presets = {}
900    for arg in values:
901        k, v = arg.split('=', 1)
902        if k not in pars:
903            # extract base name without polydispersity info
904            s = set(p.split('_pd')[0] for p in pars)
905            print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s))))
906            sys.exit(1)
907        presets[k] = float(v) if not k.endswith('type') else v
908
909    # randomize parameters
910    #pars.update(set_pars)  # set value before random to control range
911    if opts['seed'] > -1:
912        pars = randomize_pars(model_info, pars, seed=opts['seed'])
913        print("Randomize using -random=%i"%opts['seed'])
914    if opts['mono']:
915        pars = suppress_pd(pars)
916    pars.update(presets)  # set value after random to control value
917    #import pprint; pprint.pprint(model_info)
918    constrain_pars(model_info, pars)
919    if use_sasview:
920        constrain_new_to_old(model_info, pars)
921    if opts['show_pars']:
922        print(str(parlist(model_info, pars, opts['is2d'])))
923
924    # Create the computational engines
925    data, _ = make_data(opts)
926    if n1:
927        base = make_engine(model_info, data, engines[0], opts['cutoff'])
928    else:
929        base = None
930    if n2:
931        comp = make_engine(model_info, data, engines[1], opts['cutoff'])
932    else:
933        comp = None
934
935    # pylint: disable=bad-whitespace
936    # Remember it all
937    opts.update({
938        'name'      : name,
939        'def'       : model_info,
940        'n1'        : n1,
941        'n2'        : n2,
942        'presets'   : presets,
943        'pars'      : pars,
944        'data'      : data,
945        'engines'   : [base, comp],
946    })
947    # pylint: enable=bad-whitespace
948
949    return opts
950
951def explore(opts):
952    # type: (Dict[str, Any]) -> None
953    """
954    Explore the model using the Bumps GUI.
955    """
956    import wx  # type: ignore
957    from bumps.names import FitProblem  # type: ignore
958    from bumps.gui.app_frame import AppFrame  # type: ignore
959
960    problem = FitProblem(Explore(opts))
961    is_mac = "cocoa" in wx.version()
962    app = wx.App()
963    frame = AppFrame(parent=None, title="explore")
964    if not is_mac: frame.Show()
965    frame.panel.set_model(model=problem)
966    frame.panel.Layout()
967    frame.panel.aui.Split(0, wx.TOP)
968    if is_mac: frame.Show()
969    app.MainLoop()
970
971class Explore(object):
972    """
973    Bumps wrapper for a SAS model comparison.
974
975    The resulting object can be used as a Bumps fit problem so that
976    parameters can be adjusted in the GUI, with plots updated on the fly.
977    """
978    def __init__(self, opts):
979        # type: (Dict[str, Any]) -> None
980        from bumps.cli import config_matplotlib  # type: ignore
981        from . import bumps_model
982        config_matplotlib()
983        self.opts = opts
984        model_info = opts['def']
985        pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars'])
986        # Initialize parameter ranges, fixing the 2D parameters for 1D data.
987        if not opts['is2d']:
988            for p in model_info.parameters.user_parameters(is2d=False):
989                for ext in ['', '_pd', '_pd_n', '_pd_nsigma']:
990                    k = p.name+ext
991                    v = pars.get(k, None)
992                    if v is not None:
993                        v.range(*parameter_range(k, v.value))
994        else:
995            for k, v in pars.items():
996                v.range(*parameter_range(k, v.value))
997
998        self.pars = pars
999        self.pd_types = pd_types
1000        self.limits = None
1001
1002    def numpoints(self):
1003        # type: () -> int
1004        """
1005        Return the number of points.
1006        """
1007        return len(self.pars) + 1  # so dof is 1
1008
1009    def parameters(self):
1010        # type: () -> Any   # Dict/List hierarchy of parameters
1011        """
1012        Return a dictionary of parameters.
1013        """
1014        return self.pars
1015
1016    def nllf(self):
1017        # type: () -> float
1018        """
1019        Return cost.
1020        """
1021        # pylint: disable=no-self-use
1022        return 0.  # No nllf
1023
1024    def plot(self, view='log'):
1025        # type: (str) -> None
1026        """
1027        Plot the data and residuals.
1028        """
1029        pars = dict((k, v.value) for k, v in self.pars.items())
1030        pars.update(self.pd_types)
1031        self.opts['pars'] = pars
1032        limits = compare(self.opts, limits=self.limits)
1033        if self.limits is None:
1034            vmin, vmax = limits
1035            self.limits = vmax*1e-7, 1.3*vmax
1036
1037
1038def main():
1039    # type: () -> None
1040    """
1041    Main program.
1042    """
1043    opts = parse_opts()
1044    if opts['explore']:
1045        explore(opts)
1046    else:
1047        compare(opts)
1048
1049if __name__ == "__main__":
1050    main()
Note: See TracBrowser for help on using the repository browser.