source: sasmodels/sasmodels/compare.py @ 2b2d431

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

play with resolution defined by fixed dtheta,dlambda

  • Property mode set to 100755
File size: 54.0 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], float) -> 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        # TODO: provide command line control of lambda and Delta lambda/lambda
670        #L, dLoL = 5, 0.14/np.sqrt(6)  # wavelength and 14% triangular FWHM
671        L, dLoL = 0, 0
672        data = empty_data1D(q, resolution=res, L=L, dL=L*dLoL)
673        index = slice(None, None)
674    return data, index
675
676DTYPE_MAP = {
677    'half': '16',
678    'fast': 'fast',
679    'single': '32',
680    'double': '64',
681    'quad': '128',
682    'f16': '16',
683    'f32': '32',
684    'f64': '64',
685    'float16': '16',
686    'float32': '32',
687    'float64': '64',
688    'float128': '128',
689    'longdouble': '128',
690}
691def eval_opencl(model_info, data, dtype='single', cutoff=0.):
692    # type: (ModelInfo, Data, str, float) -> Calculator
693    """
694    Return a model calculator using the OpenCL calculation engine.
695    """
696
697def eval_ctypes(model_info, data, dtype='double', cutoff=0.):
698    # type: (ModelInfo, Data, str, float) -> Calculator
699    """
700    Return a model calculator using the DLL calculation engine.
701    """
702    model = core.build_model(model_info, dtype=dtype, platform="dll")
703    calculator = DirectModel(data, model, cutoff=cutoff)
704    calculator.engine = "OMP%s"%DTYPE_MAP[str(model.dtype)]
705    return calculator
706
707def make_engine(model_info, data, dtype, cutoff, ngauss=0):
708    # type: (ModelInfo, Data, str, float) -> Calculator
709    """
710    Generate the appropriate calculation engine for the given datatype.
711
712    Datatypes with '!' appended are evaluated using external C DLLs rather
713    than OpenCL.
714    """
715    if ngauss:
716        set_integration_size(model_info, ngauss)
717
718    if dtype != "default" and not dtype.endswith('!') and not kernelcl.use_opencl():
719        raise RuntimeError("OpenCL not available " + kernelcl.OPENCL_ERROR)
720
721    model = core.build_model(model_info, dtype=dtype, platform="ocl")
722    calculator = DirectModel(data, model, cutoff=cutoff)
723    engine_type = calculator._model.__class__.__name__.replace('Model', '').upper()
724    bits = calculator._model.dtype.itemsize*8
725    precision = "fast" if getattr(calculator._model, 'fast', False) else str(bits)
726    calculator.engine = "%s[%s]" % (engine_type, precision)
727    return calculator
728
729def _show_invalid(data, theory):
730    # type: (Data, np.ma.ndarray) -> None
731    """
732    Display a list of the non-finite values in theory.
733    """
734    if not theory.mask.any():
735        return
736
737    if hasattr(data, 'x'):
738        bad = zip(data.x[theory.mask], theory[theory.mask])
739        print("   *** ", ", ".join("I(%g)=%g"%(x, y) for x, y in bad))
740
741
742def compare(opts, limits=None, maxdim=np.inf):
743    # type: (Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float]
744    """
745    Preform a comparison using options from the command line.
746
747    *limits* are the limits on the values to use, either to set the y-axis
748    for 1D or to set the colormap scale for 2D.  If None, then they are
749    inferred from the data and returned. When exploring using Bumps,
750    the limits are set when the model is initially called, and maintained
751    as the values are adjusted, making it easier to see the effects of the
752    parameters.
753
754    *maxdim* is the maximum value for any parameter with units of Angstrom.
755    """
756    for k in range(opts['sets']):
757        if k > 0:
758            # print a separate seed for each dataset for better reproducibility
759            new_seed = np.random.randint(1000000)
760            print("=== Set %d uses -random=%i ==="%(k+1, new_seed))
761            np.random.seed(new_seed)
762        opts['pars'] = parse_pars(opts, maxdim=maxdim)
763        if opts['pars'] is None:
764            return
765        result = run_models(opts, verbose=True)
766        if opts['plot']:
767            if opts['is2d'] and k > 0:
768                import matplotlib.pyplot as plt
769                plt.figure()
770            limits = plot_models(opts, result, limits=limits, setnum=k)
771        if opts['show_weights']:
772            base, _ = opts['engines']
773            base_pars, _ = opts['pars']
774            model_info = base._kernel.info
775            dim = base._kernel.dim
776            plot_weights(model_info, get_mesh(model_info, base_pars, dim=dim))
777    if opts['plot']:
778        import matplotlib.pyplot as plt
779        plt.show()
780    return limits
781
782def run_models(opts, verbose=False):
783    # type: (Dict[str, Any]) -> Dict[str, Any]
784    """
785    Process a parameter set, return calculation results and times.
786    """
787
788    base, comp = opts['engines']
789    base_n, comp_n = opts['count']
790    base_pars, comp_pars = opts['pars']
791    base_data, comp_data = opts['data']
792
793    comparison = comp is not None
794
795    base_time = comp_time = None
796    base_value = comp_value = resid = relerr = None
797
798    # Base calculation
799    try:
800        base_raw, base_time = time_calculation(base, base_pars, base_n)
801        base_value = np.ma.masked_invalid(base_raw)
802        if verbose:
803            print("%s t=%.2f ms, intensity=%.0f"
804                  % (base.engine, base_time, base_value.sum()))
805        _show_invalid(base_data, base_value)
806    except ImportError:
807        traceback.print_exc()
808
809    # Comparison calculation
810    if comparison:
811        try:
812            comp_raw, comp_time = time_calculation(comp, comp_pars, comp_n)
813            comp_value = np.ma.masked_invalid(comp_raw)
814            if verbose:
815                print("%s t=%.2f ms, intensity=%.0f"
816                      % (comp.engine, comp_time, comp_value.sum()))
817            _show_invalid(base_data, comp_value)
818        except ImportError:
819            traceback.print_exc()
820
821    # Compare, but only if computing both forms
822    if comparison:
823        resid = (base_value - comp_value)
824        relerr = resid/np.where(comp_value != 0., abs(comp_value), 1.0)
825        if verbose:
826            _print_stats("|%s-%s|"
827                         % (base.engine, comp.engine) + (" "*(3+len(comp.engine))),
828                         resid)
829            _print_stats("|(%s-%s)/%s|"
830                         % (base.engine, comp.engine, comp.engine),
831                         relerr)
832
833    return dict(base_value=base_value, comp_value=comp_value,
834                base_time=base_time, comp_time=comp_time,
835                resid=resid, relerr=relerr)
836
837
838def _print_stats(label, err):
839    # type: (str, np.ma.ndarray) -> None
840    # work with trimmed data, not the full set
841    sorted_err = np.sort(abs(err.compressed()))
842    if len(sorted_err) == 0:
843        print(label + "  no valid values")
844        return
845
846    p50 = int((len(sorted_err)-1)*0.50)
847    p98 = int((len(sorted_err)-1)*0.98)
848    data = [
849        "max:%.3e"%sorted_err[-1],
850        "median:%.3e"%sorted_err[p50],
851        "98%%:%.3e"%sorted_err[p98],
852        "rms:%.3e"%np.sqrt(np.mean(sorted_err**2)),
853        "zero-offset:%+.3e"%np.mean(sorted_err),
854        ]
855    print(label+"  "+"  ".join(data))
856
857
858def plot_models(opts, result, limits=None, setnum=0):
859    # type: (Dict[str, Any], Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float]
860    """
861    Plot the results from :func:`run_model`.
862    """
863    import matplotlib.pyplot as plt
864
865    base_value, comp_value = result['base_value'], result['comp_value']
866    base_time, comp_time = result['base_time'], result['comp_time']
867    resid, relerr = result['resid'], result['relerr']
868
869    have_base, have_comp = (base_value is not None), (comp_value is not None)
870    base, comp = opts['engines']
871    base_data, comp_data = opts['data']
872    use_data = (opts['datafile'] is not None) and (have_base ^ have_comp)
873
874    # Plot if requested
875    view = opts['view']
876    #view = 'log'
877    if limits is None:
878        vmin, vmax = np.inf, -np.inf
879        if have_base:
880            vmin = min(vmin, base_value.min())
881            vmax = max(vmax, base_value.max())
882        if have_comp:
883            vmin = min(vmin, comp_value.min())
884            vmax = max(vmax, comp_value.max())
885        limits = vmin, vmax
886
887    if have_base:
888        if have_comp:
889            plt.subplot(131)
890        plot_theory(base_data, base_value, view=view, use_data=use_data, limits=limits)
891        plt.title("%s t=%.2f ms"%(base.engine, base_time))
892        #cbar_title = "log I"
893    if have_comp:
894        if have_base:
895            plt.subplot(132)
896        if not opts['is2d'] and have_base:
897            plot_theory(comp_data, base_value, view=view, use_data=use_data, limits=limits)
898        plot_theory(comp_data, comp_value, view=view, use_data=use_data, limits=limits)
899        plt.title("%s t=%.2f ms"%(comp.engine, comp_time))
900        #cbar_title = "log I"
901    if have_base and have_comp:
902        plt.subplot(133)
903        if not opts['rel_err']:
904            err, errstr, errview = resid, "abs err", "linear"
905        else:
906            err, errstr, errview = abs(relerr), "rel err", "log"
907            if (err == 0.).all():
908                errview = 'linear'
909        if 0:  # 95% cutoff
910            sorted_err = np.sort(err.flatten())
911            cutoff = sorted_err[int(sorted_err.size*0.95)]
912            err[err > cutoff] = cutoff
913        #err,errstr = base/comp,"ratio"
914        # Note: base_data only since base and comp have same q values (though
915        # perhaps different resolution), and we are plotting the difference
916        # at each q
917        plot_theory(base_data, None, resid=err, view=errview, use_data=use_data)
918        plt.xscale('log' if view == 'log' and not opts['is2d'] else 'linear')
919        plt.legend(['P%d'%(k+1) for k in range(setnum+1)], loc='best')
920        plt.title("max %s = %.3g"%(errstr, abs(err).max()))
921        #cbar_title = errstr if errview=="linear" else "log "+errstr
922    #if is2D:
923    #    h = plt.colorbar()
924    #    h.ax.set_title(cbar_title)
925    fig = plt.gcf()
926    extra_title = ' '+opts['title'] if opts['title'] else ''
927    fig.suptitle(":".join(opts['name']) + extra_title)
928
929    if have_base and have_comp and opts['show_hist']:
930        plt.figure()
931        v = relerr
932        v[v == 0] = 0.5*np.min(np.abs(v[v != 0]))
933        plt.hist(np.log10(np.abs(v)), normed=1, bins=50)
934        plt.xlabel('log10(err), err = |(%s - %s) / %s|'
935                   % (base.engine, comp.engine, comp.engine))
936        plt.ylabel('P(err)')
937        plt.title('Distribution of relative error between calculation engines')
938
939    return limits
940
941
942# ===========================================================================
943#
944
945# Set of command line options.
946# Normal options such as -plot/-noplot are specified as 'name'.
947# For options such as -nq=500 which require a value use 'name='.
948#
949OPTIONS = [
950    # Plotting
951    'plot', 'noplot', 'weights',
952    'linear', 'log', 'q4',
953    'rel', 'abs',
954    'hist', 'nohist',
955    'title=',
956
957    # Data generation
958    'data=', 'noise=', 'res=', 'nq=', 'q=',
959    'lowq', 'midq', 'highq', 'exq', 'zero',
960    '2d', '1d',
961
962    # Parameter set
963    'preset', 'random', 'random=', 'sets=',
964    'demo', 'default',  # TODO: remove demo/default
965    'nopars', 'pars',
966    'sphere', 'sphere=', # integrate over a sphere in 2d with n points
967
968    # Calculation options
969    'poly', 'mono', 'cutoff=',
970    'magnetic', 'nonmagnetic',
971    'accuracy=', 'ngauss=',
972    'neval=',  # for timing...
973
974    # Precision options
975    'engine=',
976    'half', 'fast', 'single', 'double', 'single!', 'double!', 'quad!',
977
978    # Output options
979    'help', 'html', 'edit',
980    ]
981
982NAME_OPTIONS = (lambda: set(k for k in OPTIONS if not k.endswith('=')))()
983VALUE_OPTIONS = (lambda: [k[:-1] for k in OPTIONS if k.endswith('=')])()
984
985
986def columnize(items, indent="", width=79):
987    # type: (List[str], str, int) -> str
988    """
989    Format a list of strings into columns.
990
991    Returns a string with carriage returns ready for printing.
992    """
993    column_width = max(len(w) for w in items) + 1
994    num_columns = (width - len(indent)) // column_width
995    num_rows = len(items) // num_columns
996    items = items + [""] * (num_rows * num_columns - len(items))
997    columns = [items[k*num_rows:(k+1)*num_rows] for k in range(num_columns)]
998    lines = [" ".join("%-*s"%(column_width, entry) for entry in row)
999             for row in zip(*columns)]
1000    output = indent + ("\n"+indent).join(lines)
1001    return output
1002
1003
1004def get_pars(model_info, use_demo=False):
1005    # type: (ModelInfo, bool) -> ParameterSet
1006    """
1007    Extract demo parameters from the model definition.
1008    """
1009    # Get the default values for the parameters
1010    pars = {}
1011    for p in model_info.parameters.call_parameters:
1012        parts = [('', p.default)]
1013        if p.polydisperse:
1014            parts.append(('_pd', 0.0))
1015            parts.append(('_pd_n', 0))
1016            parts.append(('_pd_nsigma', 3.0))
1017            parts.append(('_pd_type', "gaussian"))
1018        for ext, val in parts:
1019            if p.length > 1:
1020                dict(("%s%d%s" % (p.id, k, ext), val)
1021                     for k in range(1, p.length+1))
1022            else:
1023                pars[p.id + ext] = val
1024
1025    # Plug in values given in demo
1026    if use_demo and model_info.demo:
1027        pars.update(model_info.demo)
1028    return pars
1029
1030INTEGER_RE = re.compile("^[+-]?[1-9][0-9]*$")
1031def isnumber(s):
1032    # type: (str) -> bool
1033    """Return True if string contains an int or float"""
1034    match = FLOAT_RE.match(s)
1035    isfloat = (match and not s[match.end():])
1036    return isfloat or INTEGER_RE.match(s)
1037
1038# For distinguishing pairs of models for comparison
1039# key-value pair separator =
1040# shell characters  | & ; <> $ % ' " \ # `
1041# model and parameter names _
1042# parameter expressions - + * / . ( )
1043# path characters including tilde expansion and windows drive ~ / :
1044# not sure about brackets [] {}
1045# maybe one of the following @ ? ^ ! ,
1046PAR_SPLIT = ','
1047def parse_opts(argv):
1048    # type: (List[str]) -> Dict[str, Any]
1049    """
1050    Parse command line options.
1051    """
1052    MODELS = core.list_models()
1053    flags = [arg for arg in argv
1054             if arg.startswith('-')]
1055    values = [arg for arg in argv
1056              if not arg.startswith('-') and '=' in arg]
1057    positional_args = [arg for arg in argv
1058                       if not arg.startswith('-') and '=' not in arg]
1059    models = "\n    ".join("%-15s"%v for v in MODELS)
1060    if len(positional_args) == 0:
1061        print(USAGE)
1062        print("\nAvailable models:")
1063        print(columnize(MODELS, indent="  "))
1064        return None
1065
1066    invalid = [o[1:] for o in flags
1067               if o[1:] not in NAME_OPTIONS
1068               and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)]
1069    if invalid:
1070        print("Invalid options: %s"%(", ".join(invalid)))
1071        return None
1072
1073    name = positional_args[-1]
1074
1075    # pylint: disable=bad-whitespace,C0321
1076    # Interpret the flags
1077    opts = {
1078        'plot'      : True,
1079        'view'      : 'log',
1080        'is2d'      : False,
1081        'qmin'      : None,
1082        'qmax'      : 0.05,
1083        'nq'        : 128,
1084        'res'       : '0.0',
1085        'noise'     : 0.0,
1086        'accuracy'  : 'Low',
1087        'cutoff'    : '0.0',
1088        'seed'      : -1,  # default to preset
1089        'mono'      : True,
1090        # Default to magnetic a magnetic moment is set on the command line
1091        'magnetic'  : False,
1092        'show_pars' : False,
1093        'show_hist' : False,
1094        'rel_err'   : True,
1095        'explore'   : False,
1096        'use_demo'  : True,
1097        'zero'      : False,
1098        'html'      : False,
1099        'title'     : None,
1100        'datafile'  : None,
1101        'sets'      : 0,
1102        'engine'    : 'default',
1103        'count'     : '1',
1104        'show_weights' : False,
1105        'sphere'    : 0,
1106        'ngauss'    : '0',
1107    }
1108    for arg in flags:
1109        if arg == '-noplot':    opts['plot'] = False
1110        elif arg == '-plot':    opts['plot'] = True
1111        elif arg == '-linear':  opts['view'] = 'linear'
1112        elif arg == '-log':     opts['view'] = 'log'
1113        elif arg == '-q4':      opts['view'] = 'q4'
1114        elif arg == '-1d':      opts['is2d'] = False
1115        elif arg == '-2d':      opts['is2d'] = True
1116        elif arg == '-exq':     opts['qmax'] = 10.0
1117        elif arg == '-highq':   opts['qmax'] = 1.0
1118        elif arg == '-midq':    opts['qmax'] = 0.2
1119        elif arg == '-lowq':    opts['qmax'] = 0.05
1120        elif arg == '-zero':    opts['zero'] = True
1121        elif arg.startswith('-nq='):       opts['nq'] = int(arg[4:])
1122        elif arg.startswith('-q='):
1123            opts['qmin'], opts['qmax'] = [float(v) for v in arg[3:].split(':')]
1124        elif arg.startswith('-res='):      opts['res'] = arg[5:]
1125        elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:])
1126        elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:])
1127        elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:]
1128        elif arg.startswith('-cutoff='):   opts['cutoff'] = arg[8:]
1129        elif arg.startswith('-title='):    opts['title'] = arg[7:]
1130        elif arg.startswith('-data='):     opts['datafile'] = arg[6:]
1131        elif arg.startswith('-engine='):   opts['engine'] = arg[8:]
1132        elif arg.startswith('-neval='):    opts['count'] = arg[7:]
1133        elif arg.startswith('-ngauss='):   opts['ngauss'] = arg[8:]
1134        elif arg.startswith('-random='):
1135            opts['seed'] = int(arg[8:])
1136            opts['sets'] = 0
1137        elif arg == '-random':
1138            opts['seed'] = np.random.randint(1000000)
1139            opts['sets'] = 0
1140        elif arg.startswith('-sphere'):
1141            opts['sphere'] = int(arg[8:]) if len(arg) > 7 else 150
1142            opts['is2d'] = True
1143        elif arg == '-preset':  opts['seed'] = -1
1144        elif arg == '-mono':    opts['mono'] = True
1145        elif arg == '-poly':    opts['mono'] = False
1146        elif arg == '-magnetic':       opts['magnetic'] = True
1147        elif arg == '-nonmagnetic':    opts['magnetic'] = False
1148        elif arg == '-pars':    opts['show_pars'] = True
1149        elif arg == '-nopars':  opts['show_pars'] = False
1150        elif arg == '-hist':    opts['show_hist'] = True
1151        elif arg == '-nohist':  opts['show_hist'] = False
1152        elif arg == '-rel':     opts['rel_err'] = True
1153        elif arg == '-abs':     opts['rel_err'] = False
1154        elif arg == '-half':    opts['engine'] = 'half'
1155        elif arg == '-fast':    opts['engine'] = 'fast'
1156        elif arg == '-single':  opts['engine'] = 'single'
1157        elif arg == '-double':  opts['engine'] = 'double'
1158        elif arg == '-single!': opts['engine'] = 'single!'
1159        elif arg == '-double!': opts['engine'] = 'double!'
1160        elif arg == '-quad!':   opts['engine'] = 'quad!'
1161        elif arg == '-edit':    opts['explore'] = True
1162        elif arg == '-demo':    opts['use_demo'] = True
1163        elif arg == '-default': opts['use_demo'] = False
1164        elif arg == '-weights': opts['show_weights'] = True
1165        elif arg == '-html':    opts['html'] = True
1166        elif arg == '-help':    opts['html'] = True
1167    # pylint: enable=bad-whitespace,C0321
1168
1169    # Magnetism forces 2D for now
1170    if opts['magnetic']:
1171        opts['is2d'] = True
1172
1173    # Force random if sets is used
1174    if opts['sets'] >= 1 and opts['seed'] < 0:
1175        opts['seed'] = np.random.randint(1000000)
1176    if opts['sets'] == 0:
1177        opts['sets'] = 1
1178
1179    # Create the computational engines
1180    if opts['qmin'] is None:
1181        opts['qmin'] = 0.001*opts['qmax']
1182
1183    comparison = any(PAR_SPLIT in v for v in values)
1184
1185    if PAR_SPLIT in name:
1186        names = name.split(PAR_SPLIT, 2)
1187        comparison = True
1188    else:
1189        names = [name]*2
1190    try:
1191        model_info = [core.load_model_info(k) for k in names]
1192    except ImportError as exc:
1193        print(str(exc))
1194        print("Could not find model; use one of:\n    " + models)
1195        return None
1196
1197    if PAR_SPLIT in opts['ngauss']:
1198        opts['ngauss'] = [int(k) for k in opts['ngauss'].split(PAR_SPLIT, 2)]
1199        comparison = True
1200    else:
1201        opts['ngauss'] = [int(opts['ngauss'])]*2
1202
1203    if PAR_SPLIT in opts['engine']:
1204        opts['engine'] = opts['engine'].split(PAR_SPLIT, 2)
1205        comparison = True
1206    else:
1207        opts['engine'] = [opts['engine']]*2
1208
1209    if PAR_SPLIT in opts['count']:
1210        opts['count'] = [int(k) for k in opts['count'].split(PAR_SPLIT, 2)]
1211        comparison = True
1212    else:
1213        opts['count'] = [int(opts['count'])]*2
1214
1215    if PAR_SPLIT in opts['cutoff']:
1216        opts['cutoff'] = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)]
1217        comparison = True
1218    else:
1219        opts['cutoff'] = [float(opts['cutoff'])]*2
1220
1221    if PAR_SPLIT in opts['res']:
1222        opts['res'] = [float(k) for k in opts['res'].split(PAR_SPLIT, 2)]
1223        comparison = True
1224    else:
1225        opts['res'] = [float(opts['res'])]*2
1226
1227    if opts['datafile'] is not None:
1228        data = load_data(os.path.expanduser(opts['datafile']))
1229    else:
1230        # Hack around the fact that make_data doesn't take a pair of resolutions
1231        res = opts['res']
1232        opts['res'] = res[0]
1233        data0, _ = make_data(opts)
1234        if res[0] != res[1]:
1235            opts['res'] = res[1]
1236            data1, _ = make_data(opts)
1237        else:
1238            data1 = data0
1239        opts['res'] = res
1240        data = data0, data1
1241
1242    base = make_engine(model_info[0], data[0], opts['engine'][0],
1243                       opts['cutoff'][0], opts['ngauss'][0])
1244    if comparison:
1245        comp = make_engine(model_info[1], data[1], opts['engine'][1],
1246                           opts['cutoff'][1], opts['ngauss'][1])
1247    else:
1248        comp = None
1249
1250    # pylint: disable=bad-whitespace
1251    # Remember it all
1252    opts.update({
1253        'data'      : data,
1254        'name'      : names,
1255        'info'      : model_info,
1256        'engines'   : [base, comp],
1257        'values'    : values,
1258    })
1259    # pylint: enable=bad-whitespace
1260
1261    # Set the integration parameters to the half sphere
1262    if opts['sphere'] > 0:
1263        set_spherical_integration_parameters(opts, opts['sphere'])
1264
1265    return opts
1266
1267def set_spherical_integration_parameters(opts, steps):
1268    # type: (Dict[str, Any], int) -> None
1269    """
1270    Set integration parameters for spherical integration over the entire
1271    surface in theta-phi coordinates.
1272    """
1273    # Set the integration parameters to the half sphere
1274    opts['values'].extend([
1275        #'theta=90',
1276        'theta_pd=%g'%(90/np.sqrt(3)),
1277        'theta_pd_n=%d'%steps,
1278        'theta_pd_type=rectangle',
1279        #'phi=0',
1280        'phi_pd=%g'%(180/np.sqrt(3)),
1281        'phi_pd_n=%d'%(2*steps),
1282        'phi_pd_type=rectangle',
1283        #'background=0',
1284    ])
1285    if 'psi' in opts['info'][0].parameters:
1286        opts['values'].extend([
1287            #'psi=0',
1288            'psi_pd=%g'%(180/np.sqrt(3)),
1289            'psi_pd_n=%d'%(2*steps),
1290            'psi_pd_type=rectangle',
1291        ])
1292
1293def parse_pars(opts, maxdim=np.inf):
1294    # type: (Dict[str, Any], float) -> Tuple[Dict[str, float], Dict[str, float]]
1295    """
1296    Generate a parameter set.
1297
1298    The default values come from the model, or a randomized model if a seed
1299    value is given.  Next, evaluate any parameter expressions, constraining
1300    the value of the parameter within and between models.  If *maxdim* is
1301    given, limit parameters with units of Angstrom to this value.
1302
1303    Returns a pair of parameter dictionaries for base and comparison models.
1304    """
1305    model_info, model_info2 = opts['info']
1306
1307    # Get demo parameters from model definition, or use default parameters
1308    # if model does not define demo parameters
1309    pars = get_pars(model_info, opts['use_demo'])
1310    pars2 = get_pars(model_info2, opts['use_demo'])
1311    pars2.update((k, v) for k, v in pars.items() if k in pars2)
1312    # randomize parameters
1313    #pars.update(set_pars)  # set value before random to control range
1314    if opts['seed'] > -1:
1315        pars = randomize_pars(model_info, pars)
1316        limit_dimensions(model_info, pars, maxdim)
1317        if model_info != model_info2:
1318            pars2 = randomize_pars(model_info2, pars2)
1319            limit_dimensions(model_info2, pars2, maxdim)
1320            # Share values for parameters with the same name
1321            for k, v in pars.items():
1322                if k in pars2:
1323                    pars2[k] = v
1324        else:
1325            pars2 = pars.copy()
1326        constrain_pars(model_info, pars)
1327        constrain_pars(model_info2, pars2)
1328    pars = suppress_pd(pars, opts['mono'])
1329    pars2 = suppress_pd(pars2, opts['mono'])
1330    pars = suppress_magnetism(pars, not opts['magnetic'])
1331    pars2 = suppress_magnetism(pars2, not opts['magnetic'])
1332
1333    # Fill in parameters given on the command line
1334    presets = {}
1335    presets2 = {}
1336    for arg in opts['values']:
1337        k, v = arg.split('=', 1)
1338        if k not in pars and k not in pars2:
1339            # extract base name without polydispersity info
1340            s = set(p.split('_pd')[0] for p in pars)
1341            print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s))))
1342            return None
1343        v1, v2 = v.split(PAR_SPLIT, 2) if PAR_SPLIT in v else (v, v)
1344        if v1 and k in pars:
1345            presets[k] = float(v1) if isnumber(v1) else v1
1346        if v2 and k in pars2:
1347            presets2[k] = float(v2) if isnumber(v2) else v2
1348
1349    # If pd given on the command line, default pd_n to 35
1350    for k, v in list(presets.items()):
1351        if k.endswith('_pd'):
1352            presets.setdefault(k+'_n', 35.)
1353    for k, v in list(presets2.items()):
1354        if k.endswith('_pd'):
1355            presets2.setdefault(k+'_n', 35.)
1356
1357    # Evaluate preset parameter expressions
1358    # Note: need to replace ':' with '_' in parameter names and expressions
1359    # in order to support math on magnetic parameters.
1360    context = MATH.copy()
1361    context['np'] = np
1362    context.update((k.replace(':', '_'), v) for k, v in pars.items())
1363    context.update((k, v) for k, v in presets.items() if isinstance(v, float))
1364    #for k,v in sorted(context.items()): print(k, v)
1365    for k, v in presets.items():
1366        if not isinstance(v, float) and not k.endswith('_type'):
1367            presets[k] = eval(v.replace(':', '_'), context)
1368    context.update(presets)
1369    context.update((k.replace(':', '_'), v) for k, v in presets2.items() if isinstance(v, float))
1370    for k, v in presets2.items():
1371        if not isinstance(v, float) and not k.endswith('_type'):
1372            presets2[k] = eval(v.replace(':', '_'), context)
1373
1374    # update parameters with presets
1375    pars.update(presets)  # set value after random to control value
1376    pars2.update(presets2)  # set value after random to control value
1377    #import pprint; pprint.pprint(model_info)
1378
1379    if opts['show_pars']:
1380        if model_info.name != model_info2.name or pars != pars2:
1381            print("==== %s ====="%model_info.name)
1382            print(str(parlist(model_info, pars, opts['is2d'])))
1383            print("==== %s ====="%model_info2.name)
1384            print(str(parlist(model_info2, pars2, opts['is2d'])))
1385        else:
1386            print(str(parlist(model_info, pars, opts['is2d'])))
1387
1388    return pars, pars2
1389
1390def show_docs(opts):
1391    # type: (Dict[str, Any]) -> None
1392    """
1393    show html docs for the model
1394    """
1395    from .generate import make_html
1396    from . import rst2html
1397
1398    info = opts['info'][0]
1399    html = make_html(info)
1400    path = os.path.dirname(info.filename)
1401    url = "file://" + path.replace("\\", "/")[2:] + "/"
1402    rst2html.view_html_wxapp(html, url)
1403
1404def explore(opts):
1405    # type: (Dict[str, Any]) -> None
1406    """
1407    explore the model using the bumps gui.
1408    """
1409    import wx  # type: ignore
1410    from bumps.names import FitProblem  # type: ignore
1411    from bumps.gui.app_frame import AppFrame  # type: ignore
1412    from bumps.gui import signal
1413
1414    is_mac = "cocoa" in wx.version()
1415    # Create an app if not running embedded
1416    app = wx.App() if wx.GetApp() is None else None
1417    model = Explore(opts)
1418    problem = FitProblem(model)
1419    frame = AppFrame(parent=None, title="explore", size=(1000, 700))
1420    if not is_mac:
1421        frame.Show()
1422    frame.panel.set_model(model=problem)
1423    frame.panel.Layout()
1424    frame.panel.aui.Split(0, wx.TOP)
1425    def _reset_parameters(event):
1426        model.revert_values()
1427        signal.update_parameters(problem)
1428    frame.Bind(wx.EVT_TOOL, _reset_parameters, frame.ToolBar.GetToolByPos(1))
1429    if is_mac:
1430        frame.Show()
1431    # If running withing an app, start the main loop
1432    if app:
1433        app.MainLoop()
1434
1435class Explore(object):
1436    """
1437    Bumps wrapper for a SAS model comparison.
1438
1439    The resulting object can be used as a Bumps fit problem so that
1440    parameters can be adjusted in the GUI, with plots updated on the fly.
1441    """
1442    def __init__(self, opts):
1443        # type: (Dict[str, Any]) -> None
1444        from bumps.cli import config_matplotlib  # type: ignore
1445        from . import bumps_model
1446        config_matplotlib()
1447        self.opts = opts
1448        opts['pars'] = list(opts['pars'])
1449        p1, p2 = opts['pars']
1450        m1, m2 = opts['info']
1451        self.fix_p2 = m1 != m2 or p1 != p2
1452        model_info = m1
1453        pars, pd_types = bumps_model.create_parameters(model_info, **p1)
1454        # Initialize parameter ranges, fixing the 2D parameters for 1D data.
1455        if not opts['is2d']:
1456            for p in model_info.parameters.user_parameters({}, is2d=False):
1457                for ext in ['', '_pd', '_pd_n', '_pd_nsigma']:
1458                    k = p.name+ext
1459                    v = pars.get(k, None)
1460                    if v is not None:
1461                        v.range(*parameter_range(k, v.value))
1462        else:
1463            for k, v in pars.items():
1464                v.range(*parameter_range(k, v.value))
1465
1466        self.pars = pars
1467        self.starting_values = dict((k, v.value) for k, v in pars.items())
1468        self.pd_types = pd_types
1469        self.limits = None
1470
1471    def revert_values(self):
1472        # type: () -> None
1473        """
1474        Restore starting values of the parameters.
1475        """
1476        for k, v in self.starting_values.items():
1477            self.pars[k].value = v
1478
1479    def model_update(self):
1480        # type: () -> None
1481        """
1482        Respond to signal that model parameters have been changed.
1483        """
1484        pass
1485
1486    def numpoints(self):
1487        # type: () -> int
1488        """
1489        Return the number of points.
1490        """
1491        return len(self.pars) + 1  # so dof is 1
1492
1493    def parameters(self):
1494        # type: () -> Any   # Dict/List hierarchy of parameters
1495        """
1496        Return a dictionary of parameters.
1497        """
1498        return self.pars
1499
1500    def nllf(self):
1501        # type: () -> float
1502        """
1503        Return cost.
1504        """
1505        # pylint: disable=no-self-use
1506        return 0.  # No nllf
1507
1508    def plot(self, view='log'):
1509        # type: (str) -> None
1510        """
1511        Plot the data and residuals.
1512        """
1513        pars = dict((k, v.value) for k, v in self.pars.items())
1514        pars.update(self.pd_types)
1515        self.opts['pars'][0] = pars
1516        if not self.fix_p2:
1517            self.opts['pars'][1] = pars
1518        result = run_models(self.opts)
1519        limits = plot_models(self.opts, result, limits=self.limits)
1520        if self.limits is None:
1521            vmin, vmax = limits
1522            self.limits = vmax*1e-7, 1.3*vmax
1523            import pylab
1524            pylab.clf()
1525            plot_models(self.opts, result, limits=self.limits)
1526
1527
1528def main(*argv):
1529    # type: (*str) -> None
1530    """
1531    Main program.
1532    """
1533    opts = parse_opts(argv)
1534    if opts is not None:
1535        if opts['seed'] > -1:
1536            print("Randomize using -random=%i"%opts['seed'])
1537            np.random.seed(opts['seed'])
1538        if opts['html']:
1539            show_docs(opts)
1540        elif opts['explore']:
1541            opts['pars'] = parse_pars(opts)
1542            if opts['pars'] is None:
1543                return
1544            explore(opts)
1545        else:
1546            compare(opts)
1547
1548if __name__ == "__main__":
1549    main(*sys.argv[1:])
Note: See TracBrowser for help on using the repository browser.