source: sasmodels/sasmodels/compare.py @ a8a1d48

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

Merge branch 'beta_approx' into ticket-1157

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