source: sasmodels/sasmodels/compare.py @ 5770493

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

sascomp: plot 2d models on separate plots when multiple parameter sets (fix)

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