source: sasmodels/sasmodels/compare.py @ fbb9397

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since fbb9397 was fbb9397, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

fix handling of value limits on plots

  • Property mode set to 100755
File size: 52.3 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 os
33import math
34import datetime
35import traceback
36import re
37
38import numpy as np  # type: ignore
39
40from . import core
41from . import kerneldll
42from . import exception
43from .data import plot_theory, empty_data1D, empty_data2D, load_data
44from .direct_model import DirectModel
45from .convert import revert_name, revert_pars, constrain_new_to_old
46from .generate import FLOAT_RE
47
48try:
49    from typing import Optional, Dict, Any, Callable, Tuple
50except Exception:
51    pass
52else:
53    from .modelinfo import ModelInfo, Parameter, ParameterSet
54    from .data import Data
55    Calculator = Callable[[float], np.ndarray]
56
57USAGE = """
58usage: sascomp model [options...] [key=val]
59
60Generate and compare SAS models.  If a single model is specified it shows
61a plot of that model.  Different models can be compared, or the same model
62with different parameters.  The same model with the same parameters can
63be compared with different calculation engines to see the effects of precision
64on the resultant values.
65
66model or model1,model2 are the names of the models to compare (see below).
67
68Options (* for default):
69
70    === data generation ===
71    -data="path" uses q, dq from the data file
72    -noise=0 sets the measurement error dI/I
73    -res=0 sets the resolution width dQ/Q if calculating with resolution
74    -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0
75    -nq=128 sets the number of Q points in the data set
76    -1d*/-2d computes 1d or 2d data
77    -zero indicates that q=0 should be included
78
79    === model parameters ===
80    -preset*/-random[=seed] preset or random parameters
81    -sets=n generates n random datasets with the seed given by -random=seed
82    -pars/-nopars* prints the parameter set or not
83    -default/-demo* use demo vs default parameters
84
85    === calculation options ===
86    -mono*/-poly force monodisperse or allow polydisperse demo parameters
87    -cutoff=1e-5* cutoff value for including a point in polydispersity
88    -magnetic/-nonmagnetic* suppress magnetism
89    -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh
90    -neval=1 sets the number of evals for more accurate timing
91
92    === precision options ===
93    -engine=default uses the default calcution precision
94    -single/-double/-half/-fast sets an OpenCL calculation engine
95    -single!/-double!/-quad! sets an OpenMP calculation engine
96    -sasview sets the sasview calculation engine
97
98    === plotting ===
99    -plot*/-noplot plots or suppress the plot of the model
100    -linear/-log*/-q4 intensity scaling on plots
101    -hist/-nohist* plot histogram of relative error
102    -abs/-rel* plot relative or absolute error
103    -title="note" adds note to the plot title, after the model name
104
105    === output options ===
106    -edit starts the parameter explorer
107    -help/-html shows the model docs instead of running the model
108
109The interpretation of quad precision depends on architecture, and may
110vary from 64-bit to 128-bit, with 80-bit floats being common (1e-19 precision).
111On unix and mac you may need single quotes around the DLL computation
112engines, such as -engine='single!,double!' since !, is treated as a history
113expansion request in the shell.
114
115Key=value pairs allow you to set specific values for the model parameters.
116Key=value1,value2 to compare different values of the same parameter. The
117value can be an expression including other parameters.
118
119Items later on the command line override those that appear earlier.
120
121Examples:
122
123    # compare single and double precision calculation for a barbell
124    sascomp barbell -engine=single,double
125
126    # generate 10 random lorentz models, with seed=27
127    sascomp lorentz -sets=10 -seed=27
128
129    # compare ellipsoid with R = R_polar = R_equatorial to sphere of radius R
130    sascomp sphere,ellipsoid radius_polar=radius radius_equatorial=radius
131
132    # model timing test requires multiple evals to perform the estimate
133    sascomp pringle -engine=single,double -timing=100,100 -noplot
134"""
135
136# Update docs with command line usage string.   This is separate from the usual
137# doc string so that we can display it at run time if there is an error.
138# lin
139__doc__ = (__doc__  # pylint: disable=redefined-builtin
140           + """
141Program description
142-------------------
143
144""" + USAGE)
145
146kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True
147
148# list of math functions for use in evaluating parameters
149MATH = dict((k,getattr(math, k)) for k in dir(math) if not k.startswith('_'))
150
151# CRUFT python 2.6
152if not hasattr(datetime.timedelta, 'total_seconds'):
153    def delay(dt):
154        """Return number date-time delta as number seconds"""
155        return dt.days * 86400 + dt.seconds + 1e-6 * dt.microseconds
156else:
157    def delay(dt):
158        """Return number date-time delta as number seconds"""
159        return dt.total_seconds()
160
161
162class push_seed(object):
163    """
164    Set the seed value for the random number generator.
165
166    When used in a with statement, the random number generator state is
167    restored after the with statement is complete.
168
169    :Parameters:
170
171    *seed* : int or array_like, optional
172        Seed for RandomState
173
174    :Example:
175
176    Seed can be used directly to set the seed::
177
178        >>> from numpy.random import randint
179        >>> push_seed(24)
180        <...push_seed object at...>
181        >>> print(randint(0,1000000,3))
182        [242082    899 211136]
183
184    Seed can also be used in a with statement, which sets the random
185    number generator state for the enclosed computations and restores
186    it to the previous state on completion::
187
188        >>> with push_seed(24):
189        ...    print(randint(0,1000000,3))
190        [242082    899 211136]
191
192    Using nested contexts, we can demonstrate that state is indeed
193    restored after the block completes::
194
195        >>> with push_seed(24):
196        ...    print(randint(0,1000000))
197        ...    with push_seed(24):
198        ...        print(randint(0,1000000,3))
199        ...    print(randint(0,1000000))
200        242082
201        [242082    899 211136]
202        899
203
204    The restore step is protected against exceptions in the block::
205
206        >>> with push_seed(24):
207        ...    print(randint(0,1000000))
208        ...    try:
209        ...        with push_seed(24):
210        ...            print(randint(0,1000000,3))
211        ...            raise Exception()
212        ...    except Exception:
213        ...        print("Exception raised")
214        ...    print(randint(0,1000000))
215        242082
216        [242082    899 211136]
217        Exception raised
218        899
219    """
220    def __init__(self, seed=None):
221        # type: (Optional[int]) -> None
222        self._state = np.random.get_state()
223        np.random.seed(seed)
224
225    def __enter__(self):
226        # type: () -> None
227        pass
228
229    def __exit__(self, exc_type, exc_value, traceback):
230        # type: (Any, BaseException, Any) -> None
231        # TODO: better typing for __exit__ method
232        np.random.set_state(self._state)
233
234def tic():
235    # type: () -> Callable[[], float]
236    """
237    Timer function.
238
239    Use "toc=tic()" to start the clock and "toc()" to measure
240    a time interval.
241    """
242    then = datetime.datetime.now()
243    return lambda: delay(datetime.datetime.now() - then)
244
245
246def set_beam_stop(data, radius, outer=None):
247    # type: (Data, float, float) -> None
248    """
249    Add a beam stop of the given *radius*.  If *outer*, make an annulus.
250
251    Note: this function does not require sasview
252    """
253    if hasattr(data, 'qx_data'):
254        q = np.sqrt(data.qx_data**2 + data.qy_data**2)
255        data.mask = (q < radius)
256        if outer is not None:
257            data.mask |= (q >= outer)
258    else:
259        data.mask = (data.x < radius)
260        if outer is not None:
261            data.mask |= (data.x >= outer)
262
263
264def parameter_range(p, v):
265    # type: (str, float) -> Tuple[float, float]
266    """
267    Choose a parameter range based on parameter name and initial value.
268    """
269    # process the polydispersity options
270    if p.endswith('_pd_n'):
271        return 0., 100.
272    elif p.endswith('_pd_nsigma'):
273        return 0., 5.
274    elif p.endswith('_pd_type'):
275        raise ValueError("Cannot return a range for a string value")
276    elif any(s in p for s in ('theta', 'phi', 'psi')):
277        # orientation in [-180,180], orientation pd in [0,45]
278        if p.endswith('_pd'):
279            return 0., 45.
280        else:
281            return -180., 180.
282    elif p.endswith('_pd'):
283        return 0., 1.
284    elif 'sld' in p:
285        return -0.5, 10.
286    elif p == 'background':
287        return 0., 10.
288    elif p == 'scale':
289        return 0., 1.e3
290    elif v < 0.:
291        return 2.*v, -2.*v
292    else:
293        return 0., (2.*v if v > 0. else 1.)
294
295
296def _randomize_one(model_info, name, value):
297    # type: (ModelInfo, str, float) -> float
298    # type: (ModelInfo, str, str) -> str
299    """
300    Randomize a single parameter.
301    """
302    # Set the amount of polydispersity/angular dispersion, but by default pd_n
303    # is zero so there is no polydispersity.  This allows us to turn on/off
304    # pd by setting pd_n, and still have randomly generated values
305    if name.endswith('_pd'):
306        par = model_info.parameters[name[:-3]]
307        if par.type == 'orientation':
308            # Let oriention variation peak around 13 degrees; 95% < 42 degrees
309            return 180*np.random.beta(2.5, 20)
310        else:
311            # Let polydispersity peak around 15%; 95% < 0.4; max=100%
312            return np.random.beta(1.5, 7)
313
314    # pd is selected globally rather than per parameter, so set to 0 for no pd
315    # In particular, when multiple pd dimensions, want to decrease the number
316    # of points per dimension for faster computation
317    if name.endswith('_pd_n'):
318        return 0
319
320    # Don't mess with distribution type for now
321    if name.endswith('_pd_type'):
322        return 'gaussian'
323
324    # type-dependent value of number of sigmas; for gaussian use 3.
325    if name.endswith('_pd_nsigma'):
326        return 3.
327
328    # background in the range [0.01, 1]
329    if name == 'background':
330        return 10**np.random.uniform(-2, 0)
331
332    # scale defaults to 0.1% to 30% volume fraction
333    if name == 'scale':
334        return 10**np.random.uniform(-3, -0.5)
335
336    # If it is a list of choices, pick one at random with equal probability
337    # In practice, the model specific random generator will override.
338    par = model_info.parameters[name]
339    if len(par.limits) > 2:  # choice list
340        return np.random.randint(len(par.limits))
341
342    # If it is a fixed range, pick from it with equal probability.
343    # For logarithmic ranges, the model will have to override.
344    if np.isfinite(par.limits).all():
345        return np.random.uniform(*par.limits)
346
347    # If the paramter is marked as an sld use the range of neutron slds
348    # TODO: ought to randomly contrast match a pair of SLDs
349    if par.type == 'sld':
350        return np.random.uniform(-0.5, 12)
351
352    # Limit magnetic SLDs to a smaller range, from zero to iron=5/A^2
353    if par.name.startswith('M0:'):
354        return np.random.uniform(0, 5)
355
356    # Guess at the random length/radius/thickness.  In practice, all models
357    # are going to set their own reasonable ranges.
358    if par.type == 'volume':
359        if ('length' in par.name or
360                'radius' in par.name or
361                'thick' in par.name):
362            return 10**np.random.uniform(2, 4)
363
364    # In the absence of any other info, select a value in [0, 2v], or
365    # [-2|v|, 2|v|] if v is negative, or [0, 1] if v is zero.  Mostly the
366    # model random parameter generators will override this default.
367    low, high = parameter_range(par.name, value)
368    limits = (max(par.limits[0], low), min(par.limits[1], high))
369    return np.random.uniform(*limits)
370
371def _random_pd(model_info, pars):
372    pd = [p for p in model_info.parameters.kernel_parameters if p.polydisperse]
373    pd_volume = []
374    pd_oriented = []
375    for p in pd:
376        if p.type == 'orientation':
377            pd_oriented.append(p.name)
378        elif p.length_control is not None:
379            n = int(pars.get(p.length_control, 1) + 0.5)
380            pd_volume.extend(p.name+str(k+1) for k in range(n))
381        elif p.length > 1:
382            pd_volume.extend(p.name+str(k+1) for k in range(p.length))
383        else:
384            pd_volume.append(p.name)
385    u = np.random.rand()
386    n = len(pd_volume)
387    if u < 0.01 or n < 1:
388        pass  # 1% chance of no polydispersity
389    elif u < 0.86 or n < 2:
390        pars[np.random.choice(pd_volume)+"_pd_n"] = 35
391    elif u < 0.99 or n < 3:
392        choices = np.random.choice(len(pd_volume), size=2)
393        pars[pd_volume[choices[0]]+"_pd_n"] = 25
394        pars[pd_volume[choices[1]]+"_pd_n"] = 10
395    else:
396        choices = np.random.choice(len(pd_volume), size=3)
397        pars[pd_volume[choices[0]]+"_pd_n"] = 25
398        pars[pd_volume[choices[1]]+"_pd_n"] = 10
399        pars[pd_volume[choices[2]]+"_pd_n"] = 5
400    if pd_oriented:
401        pars['theta_pd_n'] = 20
402        if np.random.rand() < 0.1:
403            pars['phi_pd_n'] = 5
404        if np.random.rand() < 0.1:
405            if any(p.name == 'psi' for p in model_info.parameters.kernel_parameters):
406                #print("generating psi_pd_n")
407                pars['psi_pd_n'] = 5
408
409    ## Show selected polydispersity
410    #for name, value in pars.items():
411    #    if name.endswith('_pd_n') and value > 0:
412    #        print(name, value, pars.get(name[:-5], 0), pars.get(name[:-2], 0))
413
414
415def randomize_pars(model_info, pars):
416    # type: (ModelInfo, ParameterSet) -> ParameterSet
417    """
418    Generate random values for all of the parameters.
419
420    Valid ranges for the random number generator are guessed from the name of
421    the parameter; this will not account for constraints such as cap radius
422    greater than cylinder radius in the capped_cylinder model, so
423    :func:`constrain_pars` needs to be called afterward..
424    """
425    # Note: the sort guarantees order of calls to random number generator
426    random_pars = dict((p, _randomize_one(model_info, p, v))
427                       for p, v in sorted(pars.items()))
428    if model_info.random is not None:
429        random_pars.update(model_info.random())
430    _random_pd(model_info, random_pars)
431    return random_pars
432
433
434def constrain_pars(model_info, pars):
435    # type: (ModelInfo, ParameterSet) -> None
436    """
437    Restrict parameters to valid values.
438
439    This includes model specific code for models such as capped_cylinder
440    which need to support within model constraints (cap radius more than
441    cylinder radius in this case).
442
443    Warning: this updates the *pars* dictionary in place.
444    """
445    # TODO: move the model specific code to the individual models
446    name = model_info.id
447    # if it is a product model, then just look at the form factor since
448    # none of the structure factors need any constraints.
449    if '*' in name:
450        name = name.split('*')[0]
451
452    # Suppress magnetism for python models (not yet implemented)
453    if callable(model_info.Iq):
454        pars.update(suppress_magnetism(pars))
455
456    if name == 'barbell':
457        if pars['radius_bell'] < pars['radius']:
458            pars['radius'], pars['radius_bell'] = pars['radius_bell'], pars['radius']
459
460    elif name == 'capped_cylinder':
461        if pars['radius_cap'] < pars['radius']:
462            pars['radius'], pars['radius_cap'] = pars['radius_cap'], pars['radius']
463
464    elif name == 'guinier':
465        # Limit guinier to an Rg such that Iq > 1e-30 (single precision cutoff)
466        # I(q) = A e^-(Rg^2 q^2/3) > e^-(30 ln 10)
467        # => ln A - (Rg^2 q^2/3) > -30 ln 10
468        # => Rg^2 q^2/3 < 30 ln 10 + ln A
469        # => Rg < sqrt(90 ln 10 + 3 ln A)/q
470        #q_max = 0.2  # mid q maximum
471        q_max = 1.0  # high q maximum
472        rg_max = np.sqrt(90*np.log(10) + 3*np.log(pars['scale']))/q_max
473        pars['rg'] = min(pars['rg'], rg_max)
474
475    elif name == 'pearl_necklace':
476        if pars['radius'] < pars['thick_string']:
477            pars['radius'], pars['thick_string'] = pars['thick_string'], pars['radius']
478        pass
479
480    elif name == 'rpa':
481        # Make sure phi sums to 1.0
482        if pars['case_num'] < 2:
483            pars['Phi1'] = 0.
484            pars['Phi2'] = 0.
485        elif pars['case_num'] < 5:
486            pars['Phi1'] = 0.
487        total = sum(pars['Phi'+c] for c in '1234')
488        for c in '1234':
489            pars['Phi'+c] /= total
490
491def parlist(model_info, pars, is2d):
492    # type: (ModelInfo, ParameterSet, bool) -> str
493    """
494    Format the parameter list for printing.
495    """
496    lines = []
497    parameters = model_info.parameters
498    magnetic = False
499    magnetic_pars = []
500    for p in parameters.user_parameters(pars, is2d):
501        if any(p.id.startswith(x) for x in ('M0:', 'mtheta:', 'mphi:')):
502            continue
503        if p.id.startswith('up:'):
504            magnetic_pars.append("%s=%s"%(p.id, pars.get(p.id, p.default)))
505            continue
506        fields = dict(
507            value=pars.get(p.id, p.default),
508            pd=pars.get(p.id+"_pd", 0.),
509            n=int(pars.get(p.id+"_pd_n", 0)),
510            nsigma=pars.get(p.id+"_pd_nsgima", 3.),
511            pdtype=pars.get(p.id+"_pd_type", 'gaussian'),
512            relative_pd=p.relative_pd,
513            M0=pars.get('M0:'+p.id, 0.),
514            mphi=pars.get('mphi:'+p.id, 0.),
515            mtheta=pars.get('mtheta:'+p.id, 0.),
516        )
517        lines.append(_format_par(p.name, **fields))
518        magnetic = magnetic or fields['M0'] != 0.
519    if magnetic and magnetic_pars:
520        lines.append(" ".join(magnetic_pars))
521    return "\n".join(lines)
522
523    #return "\n".join("%s: %s"%(p, v) for p, v in sorted(pars.items()))
524
525def _format_par(name, value=0., pd=0., n=0, nsigma=3., pdtype='gaussian',
526                relative_pd=False, M0=0., mphi=0., mtheta=0.):
527    # type: (str, float, float, int, float, str) -> str
528    line = "%s: %g"%(name, value)
529    if pd != 0.  and n != 0:
530        if relative_pd:
531            pd *= value
532        line += " +/- %g  (%d points in [-%g,%g] sigma %s)"\
533                % (pd, n, nsigma, nsigma, pdtype)
534    if M0 != 0.:
535        line += "  M0:%.3f  mtheta:%.1f  mphi:%.1f" % (M0, mtheta, mphi)
536    return line
537
538def suppress_pd(pars, suppress=True):
539    # type: (ParameterSet) -> ParameterSet
540    """
541    If suppress is True complete eliminate polydispersity of the model to test
542    models more quickly.  If suppress is False, make sure at least one
543    parameter is polydisperse, setting the first polydispersity parameter to
544    15% if no polydispersity is given (with no explicit demo parameters given
545    in the model, there will be no default polydispersity).
546    """
547    pars = pars.copy()
548    #print("pars=", pars)
549    if suppress:
550        for p in pars:
551            if p.endswith("_pd_n"):
552                pars[p] = 0
553    else:
554        any_pd = False
555        first_pd = None
556        for p in pars:
557            if p.endswith("_pd_n"):
558                pd = pars.get(p[:-2], 0.)
559                any_pd |= (pars[p] != 0 and pd != 0.)
560                if first_pd is None:
561                    first_pd = p
562        if not any_pd and first_pd is not None:
563            if pars[first_pd] == 0:
564                pars[first_pd] = 35
565            if first_pd[:-2] not in pars or pars[first_pd[:-2]] == 0:
566                pars[first_pd[:-2]] = 0.15
567    return pars
568
569def suppress_magnetism(pars, suppress=True):
570    # type: (ParameterSet) -> ParameterSet
571    """
572    If suppress is True complete eliminate magnetism of the model to test
573    models more quickly.  If suppress is False, make sure at least one sld
574    parameter is magnetic, setting the first parameter to have a strong
575    magnetic sld (8/A^2) at 60 degrees (with no explicit demo parameters given
576    in the model, there will be no default magnetism).
577    """
578    pars = pars.copy()
579    if suppress:
580        for p in pars:
581            if p.startswith("M0:"):
582                pars[p] = 0
583    else:
584        any_mag = False
585        first_mag = None
586        for p in pars:
587            if p.startswith("M0:"):
588                any_mag |= (pars[p] != 0)
589                if first_mag is None:
590                    first_mag = p
591        if not any_mag and first_mag is not None:
592            pars[first_mag] = 8.
593    return pars
594
595def eval_sasview(model_info, data):
596    # type: (Modelinfo, Data) -> Calculator
597    """
598    Return a model calculator using the pre-4.0 SasView models.
599    """
600    # importing sas here so that the error message will be that sas failed to
601    # import rather than the more obscure smear_selection not imported error
602    import sas
603    import sas.models
604    from sas.models.qsmearing import smear_selection
605    from sas.models.MultiplicationModel import MultiplicationModel
606    from sas.models.dispersion_models import models as dispersers
607
608    def get_model_class(name):
609        # type: (str) -> "sas.models.BaseComponent"
610        #print("new",sorted(_pars.items()))
611        __import__('sas.models.' + name)
612        ModelClass = getattr(getattr(sas.models, name, None), name, None)
613        if ModelClass is None:
614            raise ValueError("could not find model %r in sas.models"%name)
615        return ModelClass
616
617    # WARNING: ugly hack when handling model!
618    # Sasview models with multiplicity need to be created with the target
619    # multiplicity, so we cannot create the target model ahead of time for
620    # for multiplicity models.  Instead we store the model in a list and
621    # update the first element of that list with the new multiplicity model
622    # every time we evaluate.
623
624    # grab the sasview model, or create it if it is a product model
625    if model_info.composition:
626        composition_type, parts = model_info.composition
627        if composition_type == 'product':
628            P, S = [get_model_class(revert_name(p))() for p in parts]
629            model = [MultiplicationModel(P, S)]
630        else:
631            raise ValueError("sasview mixture models not supported by compare")
632    else:
633        old_name = revert_name(model_info)
634        if old_name is None:
635            raise ValueError("model %r does not exist in old sasview"
636                            % model_info.id)
637        ModelClass = get_model_class(old_name)
638        model = [ModelClass()]
639    model[0].disperser_handles = {}
640
641    # build a smearer with which to call the model, if necessary
642    smearer = smear_selection(data, model=model)
643    if hasattr(data, 'qx_data'):
644        q = np.sqrt(data.qx_data**2 + data.qy_data**2)
645        index = ((~data.mask) & (~np.isnan(data.data))
646                 & (q >= data.qmin) & (q <= data.qmax))
647        if smearer is not None:
648            smearer.model = model  # because smear_selection has a bug
649            smearer.accuracy = data.accuracy
650            smearer.set_index(index)
651            def _call_smearer():
652                smearer.model = model[0]
653                return smearer.get_value()
654            theory = _call_smearer
655        else:
656            theory = lambda: model[0].evalDistribution([data.qx_data[index],
657                                                        data.qy_data[index]])
658    elif smearer is not None:
659        theory = lambda: smearer(model[0].evalDistribution(data.x))
660    else:
661        theory = lambda: model[0].evalDistribution(data.x)
662
663    def calculator(**pars):
664        # type: (float, ...) -> np.ndarray
665        """
666        Sasview calculator for model.
667        """
668        oldpars = revert_pars(model_info, pars)
669        # For multiplicity models, create a model with the correct multiplicity
670        control = oldpars.pop("CONTROL", None)
671        if control is not None:
672            # sphericalSLD has one fewer multiplicity.  This update should
673            # happen in revert_pars, but it hasn't been called yet.
674            model[0] = ModelClass(control)
675        # paying for parameter conversion each time to keep life simple, if not fast
676        for k, v in oldpars.items():
677            if k.endswith('.type'):
678                par = k[:-5]
679                if v == 'gaussian': continue
680                cls = dispersers[v if v != 'rectangle' else 'rectangula']
681                handle = cls()
682                model[0].disperser_handles[par] = handle
683                try:
684                    model[0].set_dispersion(par, handle)
685                except Exception:
686                    exception.annotate_exception("while setting %s to %r"
687                                                 %(par, v))
688                    raise
689
690
691        #print("sasview pars",oldpars)
692        for k, v in oldpars.items():
693            name_attr = k.split('.')  # polydispersity components
694            if len(name_attr) == 2:
695                par, disp_par = name_attr
696                model[0].dispersion[par][disp_par] = v
697            else:
698                model[0].setParam(k, v)
699        return theory()
700
701    calculator.engine = "sasview"
702    return calculator
703
704DTYPE_MAP = {
705    'half': '16',
706    'fast': 'fast',
707    'single': '32',
708    'double': '64',
709    'quad': '128',
710    'f16': '16',
711    'f32': '32',
712    'f64': '64',
713    'float16': '16',
714    'float32': '32',
715    'float64': '64',
716    'float128': '128',
717    'longdouble': '128',
718}
719def eval_opencl(model_info, data, dtype='single', cutoff=0.):
720    # type: (ModelInfo, Data, str, float) -> Calculator
721    """
722    Return a model calculator using the OpenCL calculation engine.
723    """
724    if not core.HAVE_OPENCL:
725        raise RuntimeError("OpenCL not available")
726    model = core.build_model(model_info, dtype=dtype, platform="ocl")
727    calculator = DirectModel(data, model, cutoff=cutoff)
728    calculator.engine = "OCL%s"%DTYPE_MAP[str(model.dtype)]
729    return calculator
730
731def eval_ctypes(model_info, data, dtype='double', cutoff=0.):
732    # type: (ModelInfo, Data, str, float) -> Calculator
733    """
734    Return a model calculator using the DLL calculation engine.
735    """
736    model = core.build_model(model_info, dtype=dtype, platform="dll")
737    calculator = DirectModel(data, model, cutoff=cutoff)
738    calculator.engine = "OMP%s"%DTYPE_MAP[str(model.dtype)]
739    return calculator
740
741def time_calculation(calculator, pars, evals=1):
742    # type: (Calculator, ParameterSet, int) -> Tuple[np.ndarray, float]
743    """
744    Compute the average calculation time over N evaluations.
745
746    An additional call is generated without polydispersity in order to
747    initialize the calculation engine, and make the average more stable.
748    """
749    # initialize the code so time is more accurate
750    if evals > 1:
751        calculator(**suppress_pd(pars))
752    toc = tic()
753    # make sure there is at least one eval
754    value = calculator(**pars)
755    for _ in range(evals-1):
756        value = calculator(**pars)
757    average_time = toc()*1000. / evals
758    #print("I(q)",value)
759    return value, average_time
760
761def make_data(opts):
762    # type: (Dict[str, Any]) -> Tuple[Data, np.ndarray]
763    """
764    Generate an empty dataset, used with the model to set Q points
765    and resolution.
766
767    *opts* contains the options, with 'qmax', 'nq', 'res',
768    'accuracy', 'is2d' and 'view' parsed from the command line.
769    """
770    qmax, nq, res = opts['qmax'], opts['nq'], opts['res']
771    if opts['is2d']:
772        q = np.linspace(-qmax, qmax, nq)  # type: np.ndarray
773        data = empty_data2D(q, resolution=res)
774        data.accuracy = opts['accuracy']
775        set_beam_stop(data, 0.0004)
776        index = ~data.mask
777    else:
778        if opts['view'] == 'log' and not opts['zero']:
779            qmax = math.log10(qmax)
780            q = np.logspace(qmax-3, qmax, nq)
781        else:
782            q = np.linspace(0.001*qmax, qmax, nq)
783        if opts['zero']:
784            q = np.hstack((0, q))
785        data = empty_data1D(q, resolution=res)
786        index = slice(None, None)
787    return data, index
788
789def make_engine(model_info, data, dtype, cutoff):
790    # type: (ModelInfo, Data, str, float) -> Calculator
791    """
792    Generate the appropriate calculation engine for the given datatype.
793
794    Datatypes with '!' appended are evaluated using external C DLLs rather
795    than OpenCL.
796    """
797    if dtype == 'sasview':
798        return eval_sasview(model_info, data)
799    elif dtype is None or not dtype.endswith('!'):
800        return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff)
801    else:
802        return eval_ctypes(model_info, data, dtype=dtype[:-1], cutoff=cutoff)
803
804def _show_invalid(data, theory):
805    # type: (Data, np.ma.ndarray) -> None
806    """
807    Display a list of the non-finite values in theory.
808    """
809    if not theory.mask.any():
810        return
811
812    if hasattr(data, 'x'):
813        bad = zip(data.x[theory.mask], theory[theory.mask])
814        print("   *** ", ", ".join("I(%g)=%g"%(x, y) for x, y in bad))
815
816
817def compare(opts, limits=None):
818    # type: (Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float]
819    """
820    Preform a comparison using options from the command line.
821
822    *limits* are the limits on the values to use, either to set the y-axis
823    for 1D or to set the colormap scale for 2D.  If None, then they are
824    inferred from the data and returned. When exploring using Bumps,
825    the limits are set when the model is initially called, and maintained
826    as the values are adjusted, making it easier to see the effects of the
827    parameters.
828    """
829    for k in range(opts['sets']):
830        opts['pars'] = parse_pars(opts)
831        if opts['pars'] is None:
832            return
833        result = run_models(opts, verbose=True)
834        if opts['plot']:
835            limits = plot_models(opts, result, limits=limits, setnum=k)
836    if opts['plot']:
837        import matplotlib.pyplot as plt
838        plt.show()
839    return limits
840
841def run_models(opts, verbose=False):
842    # type: (Dict[str, Any]) -> Dict[str, Any]
843
844    base, comp = opts['engines']
845    base_n, comp_n = opts['count']
846    base_pars, comp_pars = opts['pars']
847    data = opts['data']
848
849    comparison = comp is not None
850
851    base_time = comp_time = None
852    base_value = comp_value = resid = relerr = None
853
854    # Base calculation
855    try:
856        base_raw, base_time = time_calculation(base, base_pars, base_n)
857        base_value = np.ma.masked_invalid(base_raw)
858        if verbose:
859            print("%s t=%.2f ms, intensity=%.0f"
860                  % (base.engine, base_time, base_value.sum()))
861        _show_invalid(data, base_value)
862    except ImportError:
863        traceback.print_exc()
864
865    # Comparison calculation
866    if comparison:
867        try:
868            comp_raw, comp_time = time_calculation(comp, comp_pars, comp_n)
869            comp_value = np.ma.masked_invalid(comp_raw)
870            if verbose:
871                print("%s t=%.2f ms, intensity=%.0f"
872                      % (comp.engine, comp_time, comp_value.sum()))
873            _show_invalid(data, comp_value)
874        except ImportError:
875            traceback.print_exc()
876
877    # Compare, but only if computing both forms
878    if comparison:
879        resid = (base_value - comp_value)
880        relerr = resid/np.where(comp_value != 0., abs(comp_value), 1.0)
881        if verbose:
882            _print_stats("|%s-%s|"
883                         % (base.engine, comp.engine) + (" "*(3+len(comp.engine))),
884                         resid)
885            _print_stats("|(%s-%s)/%s|"
886                         % (base.engine, comp.engine, comp.engine),
887                         relerr)
888
889    return dict(base_value=base_value, comp_value=comp_value,
890                base_time=base_time, comp_time=comp_time,
891                resid=resid, relerr=relerr)
892
893
894def _print_stats(label, err):
895    # type: (str, np.ma.ndarray) -> None
896    # work with trimmed data, not the full set
897    sorted_err = np.sort(abs(err.compressed()))
898    if len(sorted_err) == 0.:
899        print(label + "  no valid values")
900        return
901
902    p50 = int((len(sorted_err)-1)*0.50)
903    p98 = int((len(sorted_err)-1)*0.98)
904    data = [
905        "max:%.3e"%sorted_err[-1],
906        "median:%.3e"%sorted_err[p50],
907        "98%%:%.3e"%sorted_err[p98],
908        "rms:%.3e"%np.sqrt(np.mean(sorted_err**2)),
909        "zero-offset:%+.3e"%np.mean(sorted_err),
910        ]
911    print(label+"  "+"  ".join(data))
912
913
914def plot_models(opts, result, limits=None, setnum=0):
915    # type: (Dict[str, Any], Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float]
916    import matplotlib.pyplot as plt
917
918    base_value, comp_value = result['base_value'], result['comp_value']
919    base_time, comp_time = result['base_time'], result['comp_time']
920    resid, relerr = result['resid'], result['relerr']
921
922    have_base, have_comp = (base_value is not None), (comp_value is not None)
923    base, comp = opts['engines']
924    data = opts['data']
925    use_data = (opts['datafile'] is not None) and (have_base ^ have_comp)
926
927    # Plot if requested
928    view = opts['view']
929    if limits is None:
930        vmin, vmax = np.inf, -np.inf
931        if have_base:
932            vmin = min(vmin, base_value.min())
933            vmax = max(vmax, base_value.max())
934        if have_comp:
935            vmin = min(vmin, comp_value.min())
936            vmax = max(vmax, comp_value.max())
937        limits = vmin, vmax
938
939    if have_base:
940        if have_comp:
941            plt.subplot(131)
942        plot_theory(data, base_value, view=view, use_data=use_data, limits=limits)
943        plt.title("%s t=%.2f ms"%(base.engine, base_time))
944        #cbar_title = "log I"
945    if have_comp:
946        if have_base:
947            plt.subplot(132)
948        if not opts['is2d'] and have_base:
949            plot_theory(data, base_value, view=view, use_data=use_data, limits=limits)
950        plot_theory(data, comp_value, view=view, use_data=use_data, limits=limits)
951        plt.title("%s t=%.2f ms"%(comp.engine, comp_time))
952        #cbar_title = "log I"
953    if have_base and have_comp:
954        plt.subplot(133)
955        if not opts['rel_err']:
956            err, errstr, errview = resid, "abs err", "linear"
957        else:
958            err, errstr, errview = abs(relerr), "rel err", "log"
959        if 0:  # 95% cutoff
960            sorted = np.sort(err.flatten())
961            cutoff = sorted[int(sorted.size*0.95)]
962            err[err > cutoff] = cutoff
963        #err,errstr = base/comp,"ratio"
964        plot_theory(data, None, resid=err, view=errview, use_data=use_data)
965        if view == 'linear':
966            plt.xscale('linear')
967        plt.title("max %s = %.3g"%(errstr, abs(err).max()))
968        #cbar_title = errstr if errview=="linear" else "log "+errstr
969    #if is2D:
970    #    h = plt.colorbar()
971    #    h.ax.set_title(cbar_title)
972    fig = plt.gcf()
973    extra_title = ' '+opts['title'] if opts['title'] else ''
974    fig.suptitle(":".join(opts['name']) + extra_title)
975
976    if have_base and have_comp and opts['show_hist']:
977        plt.figure()
978        v = relerr
979        v[v == 0] = 0.5*np.min(np.abs(v[v != 0]))
980        plt.hist(np.log10(np.abs(v)), normed=1, bins=50)
981        plt.xlabel('log10(err), err = |(%s - %s) / %s|'
982                   % (base.engine, comp.engine, comp.engine))
983        plt.ylabel('P(err)')
984        plt.title('Distribution of relative error between calculation engines')
985
986    return limits
987
988
989# ===========================================================================
990#
991
992# Set of command line options.
993# Normal options such as -plot/-noplot are specified as 'name'.
994# For options such as -nq=500 which require a value use 'name='.
995#
996OPTIONS = [
997    # Plotting
998    'plot', 'noplot',
999    'linear', 'log', 'q4',
1000    'rel', 'abs',
1001    'hist', 'nohist',
1002    'title=',
1003
1004    # Data generation
1005    'data=', 'noise=', 'res=',
1006    'nq=', 'lowq', 'midq', 'highq', 'exq', 'zero',
1007    '2d', '1d',
1008
1009    # Parameter set
1010    'preset', 'random', 'random=', 'sets=',
1011    'demo', 'default',  # TODO: remove demo/default
1012    'nopars', 'pars',
1013
1014    # Calculation options
1015    'poly', 'mono', 'cutoff=',
1016    'magnetic', 'nonmagnetic',
1017    'accuracy=',
1018    'neval=',  # for timing...
1019
1020    # Precision options
1021    'engine=',
1022    'half', 'fast', 'single', 'double', 'single!', 'double!', 'quad!',
1023    'sasview',  # TODO: remove sasview 3.x support
1024
1025    # Output options
1026    'help', 'html', 'edit',
1027    ]
1028
1029NAME_OPTIONS = set(k for k in OPTIONS if not k.endswith('='))
1030VALUE_OPTIONS = [k[:-1] for k in OPTIONS if k.endswith('=')]
1031
1032
1033def columnize(items, indent="", width=79):
1034    # type: (List[str], str, int) -> str
1035    """
1036    Format a list of strings into columns.
1037
1038    Returns a string with carriage returns ready for printing.
1039    """
1040    column_width = max(len(w) for w in items) + 1
1041    num_columns = (width - len(indent)) // column_width
1042    num_rows = len(items) // num_columns
1043    items = items + [""] * (num_rows * num_columns - len(items))
1044    columns = [items[k*num_rows:(k+1)*num_rows] for k in range(num_columns)]
1045    lines = [" ".join("%-*s"%(column_width, entry) for entry in row)
1046             for row in zip(*columns)]
1047    output = indent + ("\n"+indent).join(lines)
1048    return output
1049
1050
1051def get_pars(model_info, use_demo=False):
1052    # type: (ModelInfo, bool) -> ParameterSet
1053    """
1054    Extract demo parameters from the model definition.
1055    """
1056    # Get the default values for the parameters
1057    pars = {}
1058    for p in model_info.parameters.call_parameters:
1059        parts = [('', p.default)]
1060        if p.polydisperse:
1061            parts.append(('_pd', 0.0))
1062            parts.append(('_pd_n', 0))
1063            parts.append(('_pd_nsigma', 3.0))
1064            parts.append(('_pd_type', "gaussian"))
1065        for ext, val in parts:
1066            if p.length > 1:
1067                dict(("%s%d%s" % (p.id, k, ext), val)
1068                     for k in range(1, p.length+1))
1069            else:
1070                pars[p.id + ext] = val
1071
1072    # Plug in values given in demo
1073    if use_demo and model_info.demo:
1074        pars.update(model_info.demo)
1075    return pars
1076
1077INTEGER_RE = re.compile("^[+-]?[1-9][0-9]*$")
1078def isnumber(str):
1079    match = FLOAT_RE.match(str)
1080    isfloat = (match and not str[match.end():])
1081    return isfloat or INTEGER_RE.match(str)
1082
1083# For distinguishing pairs of models for comparison
1084# key-value pair separator =
1085# shell characters  | & ; <> $ % ' " \ # `
1086# model and parameter names _
1087# parameter expressions - + * / . ( )
1088# path characters including tilde expansion and windows drive ~ / :
1089# not sure about brackets [] {}
1090# maybe one of the following @ ? ^ ! ,
1091PAR_SPLIT = ','
1092def parse_opts(argv):
1093    # type: (List[str]) -> Dict[str, Any]
1094    """
1095    Parse command line options.
1096    """
1097    MODELS = core.list_models()
1098    flags = [arg for arg in argv
1099             if arg.startswith('-')]
1100    values = [arg for arg in argv
1101              if not arg.startswith('-') and '=' in arg]
1102    positional_args = [arg for arg in argv
1103                       if not arg.startswith('-') and '=' not in arg]
1104    models = "\n    ".join("%-15s"%v for v in MODELS)
1105    if len(positional_args) == 0:
1106        print(USAGE)
1107        print("\nAvailable models:")
1108        print(columnize(MODELS, indent="  "))
1109        return None
1110
1111    invalid = [o[1:] for o in flags
1112               if o[1:] not in NAME_OPTIONS
1113               and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)]
1114    if invalid:
1115        print("Invalid options: %s"%(", ".join(invalid)))
1116        return None
1117
1118    name = positional_args[-1]
1119
1120    # pylint: disable=bad-whitespace
1121    # Interpret the flags
1122    opts = {
1123        'plot'      : True,
1124        'view'      : 'log',
1125        'is2d'      : False,
1126        'qmax'      : 0.05,
1127        'nq'        : 128,
1128        'res'       : 0.0,
1129        'noise'     : 0.0,
1130        'accuracy'  : 'Low',
1131        'cutoff'    : '0.0',
1132        'seed'      : -1,  # default to preset
1133        'mono'      : True,
1134        # Default to magnetic a magnetic moment is set on the command line
1135        'magnetic'  : False,
1136        'show_pars' : False,
1137        'show_hist' : False,
1138        'rel_err'   : True,
1139        'explore'   : False,
1140        'use_demo'  : True,
1141        'zero'      : False,
1142        'html'      : False,
1143        'title'     : None,
1144        'datafile'  : None,
1145        'sets'      : 0,
1146        'engine'    : 'default',
1147        'evals'     : '1',
1148    }
1149    for arg in flags:
1150        if arg == '-noplot':    opts['plot'] = False
1151        elif arg == '-plot':    opts['plot'] = True
1152        elif arg == '-linear':  opts['view'] = 'linear'
1153        elif arg == '-log':     opts['view'] = 'log'
1154        elif arg == '-q4':      opts['view'] = 'q4'
1155        elif arg == '-1d':      opts['is2d'] = False
1156        elif arg == '-2d':      opts['is2d'] = True
1157        elif arg == '-exq':     opts['qmax'] = 10.0
1158        elif arg == '-highq':   opts['qmax'] = 1.0
1159        elif arg == '-midq':    opts['qmax'] = 0.2
1160        elif arg == '-lowq':    opts['qmax'] = 0.05
1161        elif arg == '-zero':    opts['zero'] = True
1162        elif arg.startswith('-nq='):       opts['nq'] = int(arg[4:])
1163        elif arg.startswith('-res='):      opts['res'] = float(arg[5:])
1164        elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:])
1165        elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:])
1166        elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:]
1167        elif arg.startswith('-cutoff='):   opts['cutoff'] = arg[8:]
1168        elif arg.startswith('-random='):   opts['seed'] = int(arg[8:])
1169        elif arg.startswith('-title='):    opts['title'] = arg[7:]
1170        elif arg.startswith('-data='):     opts['datafile'] = arg[6:]
1171        elif arg.startswith('-engine='):   opts['engine'] = arg[8:]
1172        elif arg.startswith('-neval='):    opts['evals'] = arg[7:]
1173        elif arg == '-random':  opts['seed'] = np.random.randint(1000000)
1174        elif arg == '-preset':  opts['seed'] = -1
1175        elif arg == '-mono':    opts['mono'] = True
1176        elif arg == '-poly':    opts['mono'] = False
1177        elif arg == '-magnetic':       opts['magnetic'] = True
1178        elif arg == '-nonmagnetic':    opts['magnetic'] = False
1179        elif arg == '-pars':    opts['show_pars'] = True
1180        elif arg == '-nopars':  opts['show_pars'] = False
1181        elif arg == '-hist':    opts['show_hist'] = True
1182        elif arg == '-nohist':  opts['show_hist'] = False
1183        elif arg == '-rel':     opts['rel_err'] = True
1184        elif arg == '-abs':     opts['rel_err'] = False
1185        elif arg == '-half':    opts['engine'] = 'half'
1186        elif arg == '-fast':    opts['engine'] = 'fast'
1187        elif arg == '-single':  opts['engine'] = 'single'
1188        elif arg == '-double':  opts['engine'] = 'double'
1189        elif arg == '-single!': opts['engine'] = 'single!'
1190        elif arg == '-double!': opts['engine'] = 'double!'
1191        elif arg == '-quad!':   opts['engine'] = 'quad!'
1192        elif arg == '-sasview': opts['engine'] = 'sasview'
1193        elif arg == '-edit':    opts['explore'] = True
1194        elif arg == '-demo':    opts['use_demo'] = True
1195        elif arg == '-default': opts['use_demo'] = False
1196        elif arg == '-html':    opts['html'] = True
1197        elif arg == '-help':    opts['html'] = True
1198    # pylint: enable=bad-whitespace
1199
1200    # Magnetism forces 2D for now
1201    if opts['magnetic']:
1202        opts['is2d'] = True
1203
1204    # Force random if sets is used
1205    if opts['sets'] >= 1 and opts['seed'] < 0:
1206        opts['seed'] = np.random.randint(1000000)
1207    if opts['sets'] == 0:
1208        opts['sets'] = 1
1209
1210    # Create the computational engines
1211    if opts['datafile'] is not None:
1212        data = load_data(os.path.expanduser(opts['datafile']))
1213    else:
1214        data, _ = make_data(opts)
1215
1216    comparison = any(PAR_SPLIT in v for v in values)
1217    if PAR_SPLIT in name:
1218        names = name.split(PAR_SPLIT, 2)
1219        comparison = True
1220    else:
1221        names = [name]*2
1222    try:
1223        model_info = [core.load_model_info(k) for k in names]
1224    except ImportError as exc:
1225        print(str(exc))
1226        print("Could not find model; use one of:\n    " + models)
1227        return None
1228
1229    if PAR_SPLIT in opts['engine']:
1230        engine_types = opts['engine'].split(PAR_SPLIT, 2)
1231        comparison = True
1232    else:
1233        engine_types = [opts['engine']]*2
1234
1235    if PAR_SPLIT in opts['evals']:
1236        evals = [int(k) for k in opts['evals'].split(PAR_SPLIT, 2)]
1237        comparison = True
1238    else:
1239        evals = [int(opts['evals'])]*2
1240
1241    if PAR_SPLIT in opts['cutoff']:
1242        cutoff = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)]
1243        comparison = True
1244    else:
1245        cutoff = [float(opts['cutoff'])]*2
1246
1247    base = make_engine(model_info[0], data, engine_types[0], cutoff[0])
1248    if comparison:
1249        comp = make_engine(model_info[1], data, engine_types[1], cutoff[1])
1250    else:
1251        comp = None
1252
1253    # pylint: disable=bad-whitespace
1254    # Remember it all
1255    opts.update({
1256        'data'      : data,
1257        'name'      : names,
1258        'def'       : model_info,
1259        'count'     : evals,
1260        'engines'   : [base, comp],
1261        'values'    : values,
1262    })
1263    # pylint: enable=bad-whitespace
1264
1265    return opts
1266
1267def parse_pars(opts):
1268    model_info, model_info2 = opts['def']
1269
1270    # Get demo parameters from model definition, or use default parameters
1271    # if model does not define demo parameters
1272    pars = get_pars(model_info, opts['use_demo'])
1273    pars2 = get_pars(model_info2, opts['use_demo'])
1274    pars2.update((k, v) for k, v in pars.items() if k in pars2)
1275    # randomize parameters
1276    #pars.update(set_pars)  # set value before random to control range
1277    if opts['seed'] > -1:
1278        pars = randomize_pars(model_info, pars)
1279        if model_info != model_info2:
1280            pars2 = randomize_pars(model_info2, pars2)
1281            # Share values for parameters with the same name
1282            for k, v in pars.items():
1283                if k in pars2:
1284                    pars2[k] = v
1285        else:
1286            pars2 = pars.copy()
1287        constrain_pars(model_info, pars)
1288        constrain_pars(model_info2, pars2)
1289    pars = suppress_pd(pars, opts['mono'])
1290    pars2 = suppress_pd(pars2, opts['mono'])
1291    pars = suppress_magnetism(pars, not opts['magnetic'])
1292    pars2 = suppress_magnetism(pars2, not opts['magnetic'])
1293
1294    # Fill in parameters given on the command line
1295    presets = {}
1296    presets2 = {}
1297    for arg in opts['values']:
1298        k, v = arg.split('=', 1)
1299        if k not in pars and k not in pars2:
1300            # extract base name without polydispersity info
1301            s = set(p.split('_pd')[0] for p in pars)
1302            print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s))))
1303            return None
1304        v1, v2 = v.split(PAR_SPLIT, 2) if PAR_SPLIT in v else (v,v)
1305        if v1 and k in pars:
1306            presets[k] = float(v1) if isnumber(v1) else v1
1307        if v2 and k in pars2:
1308            presets2[k] = float(v2) if isnumber(v2) else v2
1309
1310    # If pd given on the command line, default pd_n to 35
1311    for k, v in list(presets.items()):
1312        if k.endswith('_pd'):
1313            presets.setdefault(k+'_n', 35.)
1314    for k, v in list(presets2.items()):
1315        if k.endswith('_pd'):
1316            presets2.setdefault(k+'_n', 35.)
1317
1318    # Evaluate preset parameter expressions
1319    context = MATH.copy()
1320    context['np'] = np
1321    context.update(pars)
1322    context.update((k, v) for k, v in presets.items() if isinstance(v, float))
1323    for k, v in presets.items():
1324        if not isinstance(v, float) and not k.endswith('_type'):
1325            presets[k] = eval(v, context)
1326    context.update(presets)
1327    context.update((k, v) for k, v in presets2.items() if isinstance(v, float))
1328    for k, v in presets2.items():
1329        if not isinstance(v, float) and not k.endswith('_type'):
1330            presets2[k] = eval(v, context)
1331
1332    # update parameters with presets
1333    pars.update(presets)  # set value after random to control value
1334    pars2.update(presets2)  # set value after random to control value
1335    #import pprint; pprint.pprint(model_info)
1336
1337    if opts['show_pars']:
1338        if model_info.name != model_info2.name or pars != pars2:
1339            print("==== %s ====="%model_info.name)
1340            print(str(parlist(model_info, pars, opts['is2d'])))
1341            print("==== %s ====="%model_info2.name)
1342            print(str(parlist(model_info2, pars2, opts['is2d'])))
1343        else:
1344            print(str(parlist(model_info, pars, opts['is2d'])))
1345
1346    return pars, pars2
1347
1348def show_docs(opts):
1349    # type: (Dict[str, Any]) -> None
1350    """
1351    show html docs for the model
1352    """
1353    import os
1354    from .generate import make_html
1355    from . import rst2html
1356
1357    info = opts['def'][0]
1358    html = make_html(info)
1359    path = os.path.dirname(info.filename)
1360    url = "file://"+path.replace("\\","/")[2:]+"/"
1361    rst2html.view_html_qtapp(html, url)
1362
1363def explore(opts):
1364    # type: (Dict[str, Any]) -> None
1365    """
1366    explore the model using the bumps gui.
1367    """
1368    import wx  # type: ignore
1369    from bumps.names import FitProblem  # type: ignore
1370    from bumps.gui.app_frame import AppFrame  # type: ignore
1371    from bumps.gui import signal
1372
1373    is_mac = "cocoa" in wx.version()
1374    # Create an app if not running embedded
1375    app = wx.App() if wx.GetApp() is None else None
1376    model = Explore(opts)
1377    problem = FitProblem(model)
1378    frame = AppFrame(parent=None, title="explore", size=(1000, 700))
1379    if not is_mac:
1380        frame.Show()
1381    frame.panel.set_model(model=problem)
1382    frame.panel.Layout()
1383    frame.panel.aui.Split(0, wx.TOP)
1384    def reset_parameters(event):
1385        model.revert_values()
1386        signal.update_parameters(problem)
1387    frame.Bind(wx.EVT_TOOL, reset_parameters, frame.ToolBar.GetToolByPos(1))
1388    if is_mac: frame.Show()
1389    # If running withing an app, start the main loop
1390    if app:
1391        app.MainLoop()
1392
1393class Explore(object):
1394    """
1395    Bumps wrapper for a SAS model comparison.
1396
1397    The resulting object can be used as a Bumps fit problem so that
1398    parameters can be adjusted in the GUI, with plots updated on the fly.
1399    """
1400    def __init__(self, opts):
1401        # type: (Dict[str, Any]) -> None
1402        from bumps.cli import config_matplotlib  # type: ignore
1403        from . import bumps_model
1404        config_matplotlib()
1405        self.opts = opts
1406        opts['pars'] = list(opts['pars'])
1407        p1, p2 = opts['pars']
1408        m1, m2 = opts['def']
1409        self.fix_p2 = m1 != m2 or p1 != p2
1410        model_info = m1
1411        pars, pd_types = bumps_model.create_parameters(model_info, **p1)
1412        # Initialize parameter ranges, fixing the 2D parameters for 1D data.
1413        if not opts['is2d']:
1414            for p in model_info.parameters.user_parameters({}, is2d=False):
1415                for ext in ['', '_pd', '_pd_n', '_pd_nsigma']:
1416                    k = p.name+ext
1417                    v = pars.get(k, None)
1418                    if v is not None:
1419                        v.range(*parameter_range(k, v.value))
1420        else:
1421            for k, v in pars.items():
1422                v.range(*parameter_range(k, v.value))
1423
1424        self.pars = pars
1425        self.starting_values = dict((k, v.value) for k, v in pars.items())
1426        self.pd_types = pd_types
1427        self.limits = None
1428
1429    def revert_values(self):
1430        for k, v in self.starting_values.items():
1431            self.pars[k].value = v
1432
1433    def model_update(self):
1434        pass
1435
1436    def numpoints(self):
1437        # type: () -> int
1438        """
1439        Return the number of points.
1440        """
1441        return len(self.pars) + 1  # so dof is 1
1442
1443    def parameters(self):
1444        # type: () -> Any   # Dict/List hierarchy of parameters
1445        """
1446        Return a dictionary of parameters.
1447        """
1448        return self.pars
1449
1450    def nllf(self):
1451        # type: () -> float
1452        """
1453        Return cost.
1454        """
1455        # pylint: disable=no-self-use
1456        return 0.  # No nllf
1457
1458    def plot(self, view='log'):
1459        # type: (str) -> None
1460        """
1461        Plot the data and residuals.
1462        """
1463        pars = dict((k, v.value) for k, v in self.pars.items())
1464        pars.update(self.pd_types)
1465        self.opts['pars'][0] = pars
1466        if not self.fix_p2:
1467            self.opts['pars'][1] = pars
1468        result = run_models(self.opts)
1469        limits = plot_models(self.opts, result, limits=self.limits)
1470        if self.limits is None:
1471            vmin, vmax = limits
1472            self.limits = vmax*1e-7, 1.3*vmax
1473            import pylab; pylab.clf()
1474            plot_models(self.opts, result, limits=self.limits)
1475
1476
1477def main(*argv):
1478    # type: (*str) -> None
1479    """
1480    Main program.
1481    """
1482    opts = parse_opts(argv)
1483    if opts is not None:
1484        if opts['seed'] > -1:
1485            print("Randomize using -random=%i"%opts['seed'])
1486            np.random.seed(opts['seed'])
1487        if opts['html']:
1488            show_docs(opts)
1489        elif opts['explore']:
1490            opts['pars'] = parse_pars(opts)
1491            if opts['pars'] is None:
1492                return
1493            explore(opts)
1494        else:
1495            compare(opts)
1496
1497if __name__ == "__main__":
1498    main(*sys.argv[1:])
Note: See TracBrowser for help on using the repository browser.