source: sasmodels/sasmodels/compare.py @ 7609046

Last change on this file since 7609046 was 7609046, checked in by GitHub <noreply@…>, 14 months ago

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

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