source: sasmodels/sasmodels/compare.py @ e2da671

Last change on this file since e2da671 was e2da671, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

Merge branch 'beta_approx' into ticket-608-user-defined-weights

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