source: sasmodels/sasmodels/compare.py @ d0fdba2

Last change on this file since d0fdba2 was d0fdba2, checked in by Paul Kienzle <pkienzle@…>, 12 months ago

document default paths in sascomp

  • 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 kerneldll
42from . import kernelcl
43from . import weights
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
47
48# pylint: disable=unused-import
49try:
50    from typing import Optional, Dict, Any, Callable, Tuple
51except ImportError:
52    pass
53else:
54    from .modelinfo import ModelInfo, Parameter, ParameterSet
55    from .data import Data
56    Calculator = Callable[[float], np.ndarray]
57# pylint: enable=unused-import
58
59USAGE = """
60usage: sascomp model [options...] [key=val]
61
62Generate and compare SAS models.  If a single model is specified it shows
63a plot of that model.  Different models can be compared, or the same model
64with different parameters.  The same model with the same parameters can
65be compared with different calculation engines to see the effects of precision
66on the resultant values.
67
68model or model1,model2 are the names of the models to compare (see below).
69
70Options (* for default):
71
72    === data generation ===
73    -data="path" uses q, dq from the data file
74    -noise=0 sets the measurement error dI/I
75    -res=0 sets the resolution width dQ/Q if calculating with resolution
76    -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0
77    -q=min:max alternative specification of qrange
78    -nq=128 sets the number of Q points in the data set
79    -1d*/-2d computes 1d or 2d data
80    -zero indicates that q=0 should be included
81
82    === model parameters ===
83    -preset*/-random[=seed] preset or random parameters
84    -sets=n generates n random datasets with the seed given by -random=seed
85    -pars/-nopars* prints the parameter set or not
86    -default/-demo* use demo vs default parameters
87    -sphere[=150] set up spherical integration over theta/phi using n points
88
89    === calculation options ===
90    -mono*/-poly force monodisperse or allow polydisperse random parameters
91    -cutoff=1e-5* cutoff value for including a point in polydispersity
92    -magnetic/-nonmagnetic* suppress magnetism
93    -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh
94    -neval=1 sets the number of evals for more accurate timing
95    -ngauss=0 overrides the number of points in the 1-D gaussian quadrature
96
97    === precision options ===
98    -engine=default uses the default calcution precision
99    -single/-double/-half/-fast sets an OpenCL calculation engine
100    -single!/-double!/-quad! sets an OpenMP calculation engine
101
102    === plotting ===
103    -plot*/-noplot plots or suppress the plot of the model
104    -linear/-log*/-q4 intensity scaling on plots
105    -hist/-nohist* plot histogram of relative error
106    -abs/-rel* plot relative or absolute error
107    -title="note" adds note to the plot title, after the model name
108    -weights shows weights plots for the polydisperse parameters
109    -profile shows the sld profile if the model has a plottable sld profile
110
111    === output options ===
112    -edit starts the parameter explorer
113    -help/-html shows the model docs instead of running the model
114
115    === environment variables ===
116    -DSAS_MODELPATH=~/.sasmodels/custom_models sets path to custom models
117    -DSAS_WEIGHTS_PATH=~/.sasview/weights sets path to custom distributions
118    -DSAS_OPENCL=vendor:device|none sets the target OpenCL 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=0 set to 1 to turn on OpenMP for the DLLs
122    -DSAS_DLL_PATH=~/.sasmodels/compiled_models sets the DLL cache
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    # In practice, the model specific random generator will override.
356    par = model_info.parameters[name]
357    if len(par.limits) > 2:  # choice list
358        return np.random.randint(len(par.limits))
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.startswith('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.startswith(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('M0:'+p.id, 0.),
554            mphi=pars.get('mphi:'+p.id, 0.),
555            mtheta=pars.get('mtheta:'+p.id, 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.startswith("M0:"):
622                pars[p] = 0
623    else:
624        any_mag = False
625        first_mag = None
626        for p in pars:
627            if p.startswith("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('!') and not kernelcl.use_opencl():
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            weights.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'  : True,
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        data = load_data(os.path.expanduser(opts['datafile']))
1292    else:
1293        # Hack around the fact that make_data doesn't take a pair of resolutions
1294        res = opts['res']
1295        opts['res'] = res[0]
1296        data0, _ = make_data(opts)
1297        if res[0] != res[1]:
1298            opts['res'] = res[1]
1299            data1, _ = make_data(opts)
1300        else:
1301            data1 = data0
1302        opts['res'] = res
1303        data = data0, data1
1304
1305    base = make_engine(model_info[0], data[0], opts['engine'][0],
1306                       opts['cutoff'][0], opts['ngauss'][0])
1307    if comparison:
1308        comp = make_engine(model_info[1], data[1], opts['engine'][1],
1309                           opts['cutoff'][1], opts['ngauss'][1])
1310    else:
1311        comp = None
1312
1313    # pylint: disable=bad-whitespace
1314    # Remember it all
1315    opts.update({
1316        'data'      : data,
1317        'name'      : names,
1318        'info'      : model_info,
1319        'engines'   : [base, comp],
1320        'values'    : values,
1321    })
1322    # pylint: enable=bad-whitespace
1323
1324    # Set the integration parameters to the half sphere
1325    if opts['sphere'] > 0:
1326        set_spherical_integration_parameters(opts, opts['sphere'])
1327
1328    return opts
1329
1330def set_spherical_integration_parameters(opts, steps):
1331    # type: (Dict[str, Any], int) -> None
1332    """
1333    Set integration parameters for spherical integration over the entire
1334    surface in theta-phi coordinates.
1335    """
1336    # Set the integration parameters to the half sphere
1337    opts['values'].extend([
1338        #'theta=90',
1339        'theta_pd=%g'%(90/np.sqrt(3)),
1340        'theta_pd_n=%d'%steps,
1341        'theta_pd_type=rectangle',
1342        #'phi=0',
1343        'phi_pd=%g'%(180/np.sqrt(3)),
1344        'phi_pd_n=%d'%(2*steps),
1345        'phi_pd_type=rectangle',
1346        #'background=0',
1347    ])
1348    if 'psi' in opts['info'][0].parameters:
1349        opts['values'].extend([
1350            #'psi=0',
1351            'psi_pd=%g'%(180/np.sqrt(3)),
1352            'psi_pd_n=%d'%(2*steps),
1353            'psi_pd_type=rectangle',
1354        ])
1355
1356def parse_pars(opts, maxdim=np.inf):
1357    # type: (Dict[str, Any], float) -> Tuple[Dict[str, float], Dict[str, float]]
1358    """
1359    Generate a parameter set.
1360
1361    The default values come from the model, or a randomized model if a seed
1362    value is given.  Next, evaluate any parameter expressions, constraining
1363    the value of the parameter within and between models.  If *maxdim* is
1364    given, limit parameters with units of Angstrom to this value.
1365
1366    Returns a pair of parameter dictionaries for base and comparison models.
1367    """
1368    model_info, model_info2 = opts['info']
1369
1370    # Get demo parameters from model definition, or use default parameters
1371    # if model does not define demo parameters
1372    pars = get_pars(model_info, opts['use_demo'])
1373    pars2 = get_pars(model_info2, opts['use_demo'])
1374    pars2.update((k, v) for k, v in pars.items() if k in pars2)
1375    # randomize parameters
1376    #pars.update(set_pars)  # set value before random to control range
1377    if opts['seed'] > -1:
1378        pars = randomize_pars(model_info, pars)
1379        limit_dimensions(model_info, pars, maxdim)
1380        if model_info != model_info2:
1381            pars2 = randomize_pars(model_info2, pars2)
1382            limit_dimensions(model_info2, pars2, maxdim)
1383            # Share values for parameters with the same name
1384            for k, v in pars.items():
1385                if k in pars2:
1386                    pars2[k] = v
1387        else:
1388            pars2 = pars.copy()
1389        constrain_pars(model_info, pars)
1390        constrain_pars(model_info2, pars2)
1391    pars = suppress_pd(pars, opts['mono'])
1392    pars2 = suppress_pd(pars2, opts['mono'])
1393    pars = suppress_magnetism(pars, not opts['magnetic'])
1394    pars2 = suppress_magnetism(pars2, not opts['magnetic'])
1395
1396    # Fill in parameters given on the command line
1397    presets = {}
1398    presets2 = {}
1399    for arg in opts['values']:
1400        k, v = arg.split('=', 1)
1401        if k not in pars and k not in pars2:
1402            # extract base name without polydispersity info
1403            s = set(p.split('_pd')[0] for p in pars)
1404            print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s))))
1405            return None
1406        v1, v2 = v.split(PAR_SPLIT, 2) if PAR_SPLIT in v else (v, v)
1407        if v1 and k in pars:
1408            presets[k] = float(v1) if isnumber(v1) else v1
1409        if v2 and k in pars2:
1410            presets2[k] = float(v2) if isnumber(v2) else v2
1411
1412    # If pd given on the command line, default pd_n to 35
1413    for k, v in list(presets.items()):
1414        if k.endswith('_pd'):
1415            presets.setdefault(k+'_n', 35.)
1416    for k, v in list(presets2.items()):
1417        if k.endswith('_pd'):
1418            presets2.setdefault(k+'_n', 35.)
1419
1420    # Evaluate preset parameter expressions
1421    # Note: need to replace ':' with '_' in parameter names and expressions
1422    # in order to support math on magnetic parameters.
1423    context = MATH.copy()
1424    context['np'] = np
1425    context.update((k.replace(':', '_'), v) for k, v in pars.items())
1426    context.update((k, v) for k, v in presets.items() if isinstance(v, float))
1427    #for k,v in sorted(context.items()): print(k, v)
1428    for k, v in presets.items():
1429        if not isinstance(v, float) and not k.endswith('_type'):
1430            presets[k] = eval(v.replace(':', '_'), context)
1431    context.update(presets)
1432    context.update((k.replace(':', '_'), v) for k, v in presets2.items() if isinstance(v, float))
1433    for k, v in presets2.items():
1434        if not isinstance(v, float) and not k.endswith('_type'):
1435            presets2[k] = eval(v.replace(':', '_'), context)
1436
1437    # update parameters with presets
1438    pars.update(presets)  # set value after random to control value
1439    pars2.update(presets2)  # set value after random to control value
1440    #import pprint; pprint.pprint(model_info)
1441
1442    # Hack to load user-defined distributions; run through all parameters
1443    # and make sure any pd_type parameter is a defined distribution.
1444    if (any(p.endswith('pd_type') and v not in weights.MODELS
1445            for p, v in pars.items()) or
1446        any(p.endswith('pd_type') and v not in weights.MODELS
1447            for p, v in pars2.items())):
1448       weights.load_weights()
1449
1450    if opts['show_pars']:
1451        if model_info.name != model_info2.name or pars != pars2:
1452            print("==== %s ====="%model_info.name)
1453            print(str(parlist(model_info, pars, opts['is2d'])))
1454            print("==== %s ====="%model_info2.name)
1455            print(str(parlist(model_info2, pars2, opts['is2d'])))
1456        else:
1457            print(str(parlist(model_info, pars, opts['is2d'])))
1458
1459    return pars, pars2
1460
1461def show_docs(opts):
1462    # type: (Dict[str, Any]) -> None
1463    """
1464    show html docs for the model
1465    """
1466    from .generate import make_html
1467    from . import rst2html
1468
1469    info = opts['info'][0]
1470    html = make_html(info)
1471    path = os.path.dirname(info.filename)
1472    url = "file://" + path.replace("\\", "/")[2:] + "/"
1473    rst2html.view_html_wxapp(html, url)
1474
1475def explore(opts):
1476    # type: (Dict[str, Any]) -> None
1477    """
1478    explore the model using the bumps gui.
1479    """
1480    import wx  # type: ignore
1481    from bumps.names import FitProblem  # type: ignore
1482    from bumps.gui.app_frame import AppFrame  # type: ignore
1483    from bumps.gui import signal
1484
1485    is_mac = "cocoa" in wx.version()
1486    # Create an app if not running embedded
1487    app = wx.App() if wx.GetApp() is None else None
1488    model = Explore(opts)
1489    problem = FitProblem(model)
1490    frame = AppFrame(parent=None, title="explore", size=(1000, 700))
1491    if not is_mac:
1492        frame.Show()
1493    frame.panel.set_model(model=problem)
1494    frame.panel.Layout()
1495    frame.panel.aui.Split(0, wx.TOP)
1496    def _reset_parameters(event):
1497        model.revert_values()
1498        signal.update_parameters(problem)
1499    frame.Bind(wx.EVT_TOOL, _reset_parameters, frame.ToolBar.GetToolByPos(1))
1500    if is_mac:
1501        frame.Show()
1502    # If running withing an app, start the main loop
1503    if app:
1504        app.MainLoop()
1505
1506class Explore(object):
1507    """
1508    Bumps wrapper for a SAS model comparison.
1509
1510    The resulting object can be used as a Bumps fit problem so that
1511    parameters can be adjusted in the GUI, with plots updated on the fly.
1512    """
1513    def __init__(self, opts):
1514        # type: (Dict[str, Any]) -> None
1515        from bumps.cli import config_matplotlib  # type: ignore
1516        from . import bumps_model
1517        config_matplotlib()
1518        self.opts = opts
1519        opts['pars'] = list(opts['pars'])
1520        p1, p2 = opts['pars']
1521        m1, m2 = opts['info']
1522        self.fix_p2 = m1 != m2 or p1 != p2
1523        model_info = m1
1524        pars, pd_types = bumps_model.create_parameters(model_info, **p1)
1525        # Initialize parameter ranges, fixing the 2D parameters for 1D data.
1526        if not opts['is2d']:
1527            for p in model_info.parameters.user_parameters({}, is2d=False):
1528                for ext in ['', '_pd', '_pd_n', '_pd_nsigma']:
1529                    k = p.name+ext
1530                    v = pars.get(k, None)
1531                    if v is not None:
1532                        v.range(*parameter_range(k, v.value))
1533        else:
1534            for k, v in pars.items():
1535                v.range(*parameter_range(k, v.value))
1536
1537        self.pars = pars
1538        self.starting_values = dict((k, v.value) for k, v in pars.items())
1539        self.pd_types = pd_types
1540        self.limits = None
1541
1542    def revert_values(self):
1543        # type: () -> None
1544        """
1545        Restore starting values of the parameters.
1546        """
1547        for k, v in self.starting_values.items():
1548            self.pars[k].value = v
1549
1550    def model_update(self):
1551        # type: () -> None
1552        """
1553        Respond to signal that model parameters have been changed.
1554        """
1555        pass
1556
1557    def numpoints(self):
1558        # type: () -> int
1559        """
1560        Return the number of points.
1561        """
1562        return len(self.pars) + 1  # so dof is 1
1563
1564    def parameters(self):
1565        # type: () -> Any   # Dict/List hierarchy of parameters
1566        """
1567        Return a dictionary of parameters.
1568        """
1569        return self.pars
1570
1571    def nllf(self):
1572        # type: () -> float
1573        """
1574        Return cost.
1575        """
1576        # pylint: disable=no-self-use
1577        return 0.  # No nllf
1578
1579    def plot(self, view='log'):
1580        # type: (str) -> None
1581        """
1582        Plot the data and residuals.
1583        """
1584        pars = dict((k, v.value) for k, v in self.pars.items())
1585        pars.update(self.pd_types)
1586        self.opts['pars'][0] = pars
1587        if not self.fix_p2:
1588            self.opts['pars'][1] = pars
1589        result = run_models(self.opts)
1590        limits = plot_models(self.opts, result, limits=self.limits)
1591        if self.limits is None:
1592            vmin, vmax = limits
1593            self.limits = vmax*1e-7, 1.3*vmax
1594            import pylab
1595            pylab.clf()
1596            plot_models(self.opts, result, limits=self.limits)
1597
1598
1599def main(*argv):
1600    # type: (*str) -> None
1601    """
1602    Main program.
1603    """
1604    opts = parse_opts(argv)
1605    if opts is not None:
1606        if opts['seed'] > -1:
1607            print("Randomize using -random=%i"%opts['seed'])
1608            np.random.seed(opts['seed'])
1609        if opts['html']:
1610            show_docs(opts)
1611        elif opts['explore']:
1612            opts['pars'] = parse_pars(opts)
1613            if opts['pars'] is None:
1614                return
1615            explore(opts)
1616        else:
1617            compare(opts)
1618
1619if __name__ == "__main__":
1620    main(*sys.argv[1:])
Note: See TracBrowser for help on using the repository browser.