source: sasmodels/sasmodels/compare.py @ aa25fc7

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

load user-defined weight functions from ~/.sasview/weights/*.py

  • Property mode set to 100755
File size: 53.4 KB
RevLine 
[8a20be5]1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
[caeb06d]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
[9cfcac8]11tag the model as double precision only.
[caeb06d]12
[9cfcac8]13Run using ./compare.sh (Linux, Mac) or compare.bat (Windows) in the
[caeb06d]14sasmodels root to see the command line options.
15
[9cfcac8]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
[caeb06d]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
[190fc2b]31import sys
[a769b54]32import os
[190fc2b]33import math
34import datetime
35import traceback
[ff1fff5]36import re
[190fc2b]37
[7ae2b7f]38import numpy as np  # type: ignore
[190fc2b]39
40from . import core
41from . import kerneldll
[3221de0]42from . import kernelcl
[aa25fc7]43from . import weights
[a769b54]44from .data import plot_theory, empty_data1D, empty_data2D, load_data
[3c24ccd]45from .direct_model import DirectModel, get_mesh
[ff31782]46from .generate import FLOAT_RE, set_integration_size
[190fc2b]47
[32398dc]48# pylint: disable=unused-import
[dd7fc12]49try:
50    from typing import Optional, Dict, Any, Callable, Tuple
[32398dc]51except ImportError:
[dd7fc12]52    pass
53else:
54    from .modelinfo import ModelInfo, Parameter, ParameterSet
55    from .data import Data
[8d62008]56    Calculator = Callable[[float], np.ndarray]
[32398dc]57# pylint: enable=unused-import
[dd7fc12]58
[caeb06d]59USAGE = """
[bb39b4a]60usage: sascomp model [options...] [key=val]
[caeb06d]61
[bb39b4a]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.
[caeb06d]67
[8c65a33]68model or model1,model2 are the names of the models to compare (see below).
[caeb06d]69
70Options (* for default):
71
[bb39b4a]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
[caeb06d]76    -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0
[ced5bd2]77    -q=min:max alternative specification of qrange
[caeb06d]78    -nq=128 sets the number of Q points in the data set
79    -1d*/-2d computes 1d or 2d data
[bb39b4a]80    -zero indicates that q=0 should be included
81
82    === model parameters ===
[caeb06d]83    -preset*/-random[=seed] preset or random parameters
[d9ec8f9]84    -sets=n generates n random datasets with the seed given by -random=seed
[caeb06d]85    -pars/-nopars* prints the parameter set or not
[98d6cfc]86    -default/-demo* use demo vs default parameters
[e3571cb]87    -sphere[=150] set up spherical integration over theta/phi using n points
[caeb06d]88
[bb39b4a]89    === calculation options ===
[e3571cb]90    -mono*/-poly force monodisperse or allow polydisperse random parameters
[bb39b4a]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
[765eb0e]94    -neval=1 sets the number of evals for more accurate timing
[2a7e20e]95    -ngauss=0 overrides the number of points in the 1-D gaussian quadrature
[caeb06d]96
[bb39b4a]97    === precision options ===
[8698a0d]98    -engine=default uses the default calcution precision
[caeb06d]99    -single/-double/-half/-fast sets an OpenCL calculation engine
100    -single!/-double!/-quad! sets an OpenMP calculation engine
101
[bb39b4a]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
[3c24ccd]108    -weights shows weights plots for the polydisperse parameters
[bb39b4a]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
[8698a0d]117engines, such as -engine='single!,double!' since !, is treated as a history
[bb39b4a]118expansion request in the shell.
[caeb06d]119
120Key=value pairs allow you to set specific values for the model parameters.
[bb39b4a]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
[8698a0d]129    sascomp barbell -engine=single,double
[bb39b4a]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
[8698a0d]138    sascomp pringle -engine=single,double -timing=100,100 -noplot
[caeb06d]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
[d15a908]144__doc__ = (__doc__  # pylint: disable=redefined-builtin
145           + """
[caeb06d]146Program description
147-------------------
148
[bb39b4a]149""" + USAGE)
[caeb06d]150
[750ffa5]151kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True
[87985ca]152
[110f69c]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()
[248561a]161
[7cf2cfd]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
[4f2478e]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()
[dd7fc12]223        ...    except Exception:
[4f2478e]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):
[dd7fc12]232        # type: (Optional[int]) -> None
[4f2478e]233        self._state = np.random.get_state()
234        np.random.seed(seed)
235
236    def __enter__(self):
[dd7fc12]237        # type: () -> None
238        pass
[4f2478e]239
[32398dc]240    def __exit__(self, exc_type, exc_value, trace):
[dd7fc12]241        # type: (Any, BaseException, Any) -> None
[4f2478e]242        np.random.set_state(self._state)
243
[7cf2cfd]244def tic():
[dd7fc12]245    # type: () -> Callable[[], float]
[7cf2cfd]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):
[dd7fc12]257    # type: (Data, float, float) -> None
[7cf2cfd]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
[8a20be5]271
[ec7e360]272def parameter_range(p, v):
[dd7fc12]273    # type: (str, float) -> Tuple[float, float]
[87985ca]274    """
[ec7e360]275    Choose a parameter range based on parameter name and initial value.
[87985ca]276    """
[8bd7b77]277    # process the polydispersity options
[ec7e360]278    if p.endswith('_pd_n'):
[dd7fc12]279        return 0., 100.
[ec7e360]280    elif p.endswith('_pd_nsigma'):
[dd7fc12]281        return 0., 5.
[ec7e360]282    elif p.endswith('_pd_type'):
[dd7fc12]283        raise ValueError("Cannot return a range for a string value")
[caeb06d]284    elif any(s in p for s in ('theta', 'phi', 'psi')):
[87985ca]285        # orientation in [-180,180], orientation pd in [0,45]
286        if p.endswith('_pd'):
[e3571cb]287            return 0., 180.
[87985ca]288        else:
[dd7fc12]289            return -180., 180.
[87985ca]290    elif p.endswith('_pd'):
[dd7fc12]291        return 0., 1.
[8bd7b77]292    elif 'sld' in p:
[dd7fc12]293        return -0.5, 10.
[eb46451]294    elif p == 'background':
[dd7fc12]295        return 0., 10.
[eb46451]296    elif p == 'scale':
[dd7fc12]297        return 0., 1.e3
298    elif v < 0.:
299        return 2.*v, -2.*v
[87985ca]300    else:
[dd7fc12]301        return 0., (2.*v if v > 0. else 1.)
[87985ca]302
[4f2478e]303
[0bdddc2]304def _randomize_one(model_info, name, value):
[dd7fc12]305    # type: (ModelInfo, str, float) -> float
306    # type: (ModelInfo, str, str) -> str
[ec7e360]307    """
[caeb06d]308    Randomize a single parameter.
[ec7e360]309    """
[31df0c9]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
[0bdddc2]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)
[8bd7b77]321
[31df0c9]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
[0bdddc2]325    if name.endswith('_pd_n'):
326        return 0
327
[31df0c9]328    # Don't mess with distribution type for now
[0bdddc2]329    if name.endswith('_pd_type'):
330        return 'gaussian'
331
[31df0c9]332    # type-dependent value of number of sigmas; for gaussian use 3.
[0bdddc2]333    if name.endswith('_pd_nsigma'):
334        return 3.
[8bd7b77]335
[31df0c9]336    # background in the range [0.01, 1]
[0bdddc2]337    if name == 'background':
[31df0c9]338        return 10**np.random.uniform(-2, 0)
[0bdddc2]339
[31df0c9]340    # scale defaults to 0.1% to 30% volume fraction
[0bdddc2]341    if name == 'scale':
[31df0c9]342        return 10**np.random.uniform(-3, -0.5)
[0bdddc2]343
[31df0c9]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.
[0bdddc2]346    par = model_info.parameters[name]
[8bd7b77]347    if len(par.limits) > 2:  # choice list
348        return np.random.randint(len(par.limits))
349
[31df0c9]350    # If it is a fixed range, pick from it with equal probability.
351    # For logarithmic ranges, the model will have to override.
[0bdddc2]352    if np.isfinite(par.limits).all():
353        return np.random.uniform(*par.limits)
354
[31df0c9]355    # If the paramter is marked as an sld use the range of neutron slds
[0f6c41c]356    # TODO: ought to randomly contrast match a pair of SLDs
[0bdddc2]357    if par.type == 'sld':
358        return np.random.uniform(-0.5, 12)
[8bd7b77]359
[0f6c41c]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
[31df0c9]364    # Guess at the random length/radius/thickness.  In practice, all models
365    # are going to set their own reasonable ranges.
[0bdddc2]366    if par.type == 'volume':
367        if ('length' in par.name or
368                'radius' in par.name or
369                'thick' in par.name):
[31df0c9]370            return 10**np.random.uniform(2, 4)
[0bdddc2]371
[31df0c9]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.
[0bdddc2]375    low, high = parameter_range(par.name, value)
376    limits = (max(par.limits[0], low), min(par.limits[1], high))
[8bd7b77]377    return np.random.uniform(*limits)
[cd3dba0]378
[109d963]379def _random_pd(model_info, pars):
[110f69c]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    """
[109d963]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:
[232bb12]399            n = int(pars.get(p.length_control, 1) + 0.5)
[109d963]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:
[4553dae]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
[109d963]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
[caeb06d]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    """
[0bdddc2]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())
[109d963]450    _random_pd(model_info, random_pars)
[dd7fc12]451    return random_pars
[cd3dba0]452
[109d963]453
[e3571cb]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:
[110f69c]462            pars[p.name] = maxdim*10**np.random.uniform(-3, 0)
[e3571cb]463
[17bbadd]464def constrain_pars(model_info, pars):
[dd7fc12]465    # type: (ModelInfo, ParameterSet) -> None
[9a66e65]466    """
467    Restrict parameters to valid values.
[caeb06d]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).
[dd7fc12]472
473    Warning: this updates the *pars* dictionary in place.
[9a66e65]474    """
[109d963]475    # TODO: move the model specific code to the individual models
[6d6508e]476    name = model_info.id
[17bbadd]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
[f72d70a]482    # Suppress magnetism for python models (not yet implemented)
483    if callable(model_info.Iq):
484        pars.update(suppress_magnetism(pars))
485
[158cee4]486    if name == 'barbell':
487        if pars['radius_bell'] < pars['radius']:
488            pars['radius'], pars['radius_bell'] = pars['radius_bell'], pars['radius']
[b514adf]489
[158cee4]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)
[48462b0]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
[b514adf]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
[caeb06d]503        pars['rg'] = min(pars['rg'], rg_max)
[cd3dba0]504
[3e8ea5d]505    elif name == 'pearl_necklace':
506        if pars['radius'] < pars['thick_string']:
507            pars['radius'], pars['thick_string'] = pars['thick_string'], pars['radius']
508
[158cee4]509    elif name == 'rpa':
[82c299f]510        # Make sure phi sums to 1.0
511        if pars['case_num'] < 2:
[8bd7b77]512            pars['Phi1'] = 0.
513            pars['Phi2'] = 0.
[82c299f]514        elif pars['case_num'] < 5:
[8bd7b77]515            pars['Phi1'] = 0.
516        total = sum(pars['Phi'+c] for c in '1234')
517        for c in '1234':
[82c299f]518            pars['Phi'+c] /= total
519
[d6850fa]520def parlist(model_info, pars, is2d):
[dd7fc12]521    # type: (ModelInfo, ParameterSet, bool) -> str
[caeb06d]522    """
523    Format the parameter list for printing.
524    """
[e3571cb]525    is2d = True
[a4a7308]526    lines = []
[6d6508e]527    parameters = model_info.parameters
[0b040de]528    magnetic = False
[97d89af]529    magnetic_pars = []
[d19962c]530    for p in parameters.user_parameters(pars, is2d):
[0b040de]531        if any(p.id.startswith(x) for x in ('M0:', 'mtheta:', 'mphi:')):
532            continue
[97d89af]533        if p.id.startswith('up:'):
534            magnetic_pars.append("%s=%s"%(p.id, pars.get(p.id, p.default)))
[0b040de]535            continue
[d19962c]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.),
[dd7fc12]541            pdtype=pars.get(p.id+"_pd_type", 'gaussian'),
[bd49c79]542            relative_pd=p.relative_pd,
[0b040de]543            M0=pars.get('M0:'+p.id, 0.),
544            mphi=pars.get('mphi:'+p.id, 0.),
545            mtheta=pars.get('mtheta:'+p.id, 0.),
[dd7fc12]546        )
[d19962c]547        lines.append(_format_par(p.name, **fields))
[0b040de]548        magnetic = magnetic or fields['M0'] != 0.
[97d89af]549    if magnetic and magnetic_pars:
550        lines.append(" ".join(magnetic_pars))
[a4a7308]551    return "\n".join(lines)
552
553    #return "\n".join("%s: %s"%(p, v) for p, v in sorted(pars.items()))
554
[bd49c79]555def _format_par(name, value=0., pd=0., n=0, nsigma=3., pdtype='gaussian',
[0b040de]556                relative_pd=False, M0=0., mphi=0., mtheta=0.):
[dd7fc12]557    # type: (str, float, float, int, float, str) -> str
[a4a7308]558    line = "%s: %g"%(name, value)
559    if pd != 0.  and n != 0:
[bd49c79]560        if relative_pd:
561            pd *= value
[a4a7308]562        line += " +/- %g  (%d points in [-%g,%g] sigma %s)"\
[dd7fc12]563                % (pd, n, nsigma, nsigma, pdtype)
[0b040de]564    if M0 != 0.:
[b76191e]565        line += "  M0:%.3f  mtheta:%.1f  mphi:%.1f" % (M0, mtheta, mphi)
[a4a7308]566    return line
[87985ca]567
[97d89af]568def suppress_pd(pars, suppress=True):
[dd7fc12]569    # type: (ParameterSet) -> ParameterSet
[87985ca]570    """
[97d89af]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).
[87985ca]576    """
[f4f3919]577    pars = pars.copy()
[4553dae]578    #print("pars=", pars)
[97d89af]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"):
[4553dae]588                pd = pars.get(p[:-2], 0.)
589                any_pd |= (pars[p] != 0 and pd != 0.)
[97d89af]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
[4553dae]595            if first_pd[:-2] not in pars or pars[first_pd[:-2]] == 0:
[97d89af]596                pars[first_pd[:-2]] = 0.15
[f4f3919]597    return pars
[87985ca]598
[97d89af]599def suppress_magnetism(pars, suppress=True):
[0b040de]600    # type: (ParameterSet) -> ParameterSet
601    """
[97d89af]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).
[0b040de]607    """
608    pars = pars.copy()
[97d89af]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.
[0b040de]623    return pars
624
[ec7e360]625
[b32dafd]626def time_calculation(calculator, pars, evals=1):
[dd7fc12]627    # type: (Calculator, ParameterSet, int) -> Tuple[np.ndarray, float]
[caeb06d]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    """
[ec7e360]634    # initialize the code so time is more accurate
[b32dafd]635    if evals > 1:
[dd7fc12]636        calculator(**suppress_pd(pars))
[216a9e1]637    toc = tic()
[dd7fc12]638    # make sure there is at least one eval
639    value = calculator(**pars)
[b32dafd]640    for _ in range(evals-1):
[7cf2cfd]641        value = calculator(**pars)
[b32dafd]642    average_time = toc()*1000. / evals
[f2f67a6]643    #print("I(q)",value)
[216a9e1]644    return value, average_time
645
[ec7e360]646def make_data(opts):
[dd7fc12]647    # type: (Dict[str, Any]) -> Tuple[Data, np.ndarray]
[caeb06d]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    """
[ced5bd2]655    qmin, qmax, nq, res = opts['qmin'], opts['qmax'], opts['nq'], opts['res']
[ec7e360]656    if opts['is2d']:
[dd7fc12]657        q = np.linspace(-qmax, qmax, nq)  # type: np.ndarray
658        data = empty_data2D(q, resolution=res)
[ec7e360]659        data.accuracy = opts['accuracy']
[376b0ee]660        set_beam_stop(data, qmin)
[87985ca]661        index = ~data.mask
[216a9e1]662    else:
[e78edc4]663        if opts['view'] == 'log' and not opts['zero']:
[ced5bd2]664            q = np.logspace(math.log10(qmin), math.log10(qmax), nq)
[b89f519]665        else:
[ced5bd2]666            q = np.linspace(qmin, qmax, nq)
[e78edc4]667        if opts['zero']:
668            q = np.hstack((0, q))
[ec7e360]669        data = empty_data1D(q, resolution=res)
[216a9e1]670        index = slice(None, None)
671    return data, index
672
[3221de0]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
[ff31782]704def make_engine(model_info, data, dtype, cutoff, ngauss=0):
[dd7fc12]705    # type: (ModelInfo, Data, str, float) -> Calculator
[caeb06d]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    """
[ff31782]712    if ngauss:
713        set_integration_size(model_info, ngauss)
714
[3221de0]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)
[d86f0fc]720    engine_type = calculator._model.__class__.__name__.replace('Model', '').upper()
[3221de0]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
[87985ca]725
[e78edc4]726def _show_invalid(data, theory):
[dd7fc12]727    # type: (Data, np.ma.ndarray) -> None
728    """
729    Display a list of the non-finite values in theory.
730    """
[e78edc4]731    if not theory.mask.any():
732        return
733
734    if hasattr(data, 'x'):
735        bad = zip(data.x[theory.mask], theory[theory.mask])
[dd7fc12]736        print("   *** ", ", ".join("I(%g)=%g"%(x, y) for x, y in bad))
[e78edc4]737
738
[e3571cb]739def compare(opts, limits=None, maxdim=np.inf):
[dd7fc12]740    # type: (Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float]
[caeb06d]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.
[e3571cb]750
751    *maxdim* is the maximum value for any parameter with units of Angstrom.
[caeb06d]752    """
[0bdddc2]753    for k in range(opts['sets']):
[5770493]754        if k > 0:
[e3571cb]755            # print a separate seed for each dataset for better reproducibility
756            new_seed = np.random.randint(1000000)
[5770493]757            print("=== Set %d uses -random=%i ==="%(k+1, new_seed))
[e3571cb]758            np.random.seed(new_seed)
759        opts['pars'] = parse_pars(opts, maxdim=maxdim)
[8f04da4]760        if opts['pars'] is None:
761            return
[0bdddc2]762        result = run_models(opts, verbose=True)
763        if opts['plot']:
[5770493]764            if opts['is2d'] and k > 0:
765                import matplotlib.pyplot as plt
766                plt.figure()
[0bdddc2]767            limits = plot_models(opts, result, limits=limits, setnum=k)
[3c24ccd]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
[aa25fc7]773            weights.plot_weights(model_info, get_mesh(model_info, base_pars, dim=dim))
[0bdddc2]774    if opts['plot']:
775        import matplotlib.pyplot as plt
776        plt.show()
[fbb9397]777    return limits
[ca9e54e]778
779def run_models(opts, verbose=False):
780    # type: (Dict[str, Any]) -> Dict[str, Any]
[110f69c]781    """
782    Process a parameter set, return calculation results and times.
783    """
[ca9e54e]784
[bb39b4a]785    base, comp = opts['engines']
786    base_n, comp_n = opts['count']
787    base_pars, comp_pars = opts['pars']
[ec7e360]788    data = opts['data']
[87985ca]789
[bb39b4a]790    comparison = comp is not None
[ca9e54e]791
[dd7fc12]792    base_time = comp_time = None
793    base_value = comp_value = resid = relerr = None
794
[4b41184]795    # Base calculation
[bb39b4a]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()
[4b41184]805
806    # Comparison calculation
[bb39b4a]807    if comparison:
[7cf2cfd]808        try:
[bb39b4a]809            comp_raw, comp_time = time_calculation(comp, comp_pars, comp_n)
[dd7fc12]810            comp_value = np.ma.masked_invalid(comp_raw)
[ca9e54e]811            if verbose:
812                print("%s t=%.2f ms, intensity=%.0f"
813                      % (comp.engine, comp_time, comp_value.sum()))
[e78edc4]814            _show_invalid(data, comp_value)
[7cf2cfd]815        except ImportError:
[5753e4e]816            traceback.print_exc()
[87985ca]817
818    # Compare, but only if computing both forms
[bb39b4a]819    if comparison:
[ec7e360]820        resid = (base_value - comp_value)
[b32dafd]821        relerr = resid/np.where(comp_value != 0., abs(comp_value), 1.0)
[ca9e54e]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()))
[e65c3ba]839    if len(sorted_err) == 0:
[ca9e54e]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
[fbb9397]855def plot_models(opts, result, limits=None, setnum=0):
[ca9e54e]856    # type: (Dict[str, Any], Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float]
[110f69c]857    """
858    Plot the results from :func:`run_model`.
859    """
[fbb9397]860    import matplotlib.pyplot as plt
861
[97d89af]862    base_value, comp_value = result['base_value'], result['comp_value']
[ca9e54e]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)
[bb39b4a]867    base, comp = opts['engines']
[ca9e54e]868    data = opts['data']
[630156b]869    use_data = (opts['datafile'] is not None) and (have_base ^ have_comp)
[87985ca]870
871    # Plot if requested
[ec7e360]872    view = opts['view']
[fbb9397]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
[013adb7]882
[ca9e54e]883    if have_base:
[bb39b4a]884        if have_comp:
885            plt.subplot(131)
[a769b54]886        plot_theory(data, base_value, view=view, use_data=use_data, limits=limits)
[af92b73]887        plt.title("%s t=%.2f ms"%(base.engine, base_time))
[ec7e360]888        #cbar_title = "log I"
[ca9e54e]889    if have_comp:
[bb39b4a]890        if have_base:
891            plt.subplot(132)
[ca9e54e]892        if not opts['is2d'] and have_base:
[a769b54]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)
[af92b73]895        plt.title("%s t=%.2f ms"%(comp.engine, comp_time))
[7cf2cfd]896        #cbar_title = "log I"
[ca9e54e]897    if have_base and have_comp:
[87985ca]898        plt.subplot(133)
[d5e650d]899        if not opts['rel_err']:
[caeb06d]900            err, errstr, errview = resid, "abs err", "linear"
[29f5536]901        else:
[caeb06d]902            err, errstr, errview = abs(relerr), "rel err", "log"
[ced5bd2]903            if (err == 0.).all():
904                errview = 'linear'
[158cee4]905        if 0:  # 95% cutoff
[110f69c]906            sorted_err = np.sort(err.flatten())
907            cutoff = sorted_err[int(sorted_err.size*0.95)]
[bb39b4a]908            err[err > cutoff] = cutoff
[4b41184]909        #err,errstr = base/comp,"ratio"
[a769b54]910        plot_theory(data, None, resid=err, view=errview, use_data=use_data)
[3bfd924]911        plt.xscale('log' if view == 'log' and not opts['is2d'] else 'linear')
[e3571cb]912        plt.legend(['P%d'%(k+1) for k in range(setnum+1)], loc='best')
[e78edc4]913        plt.title("max %s = %.3g"%(errstr, abs(err).max()))
[7cf2cfd]914        #cbar_title = errstr if errview=="linear" else "log "+errstr
915    #if is2D:
916    #    h = plt.colorbar()
917    #    h.ax.set_title(cbar_title)
[0c24a82]918    fig = plt.gcf()
[a0d75ce]919    extra_title = ' '+opts['title'] if opts['title'] else ''
[ff1fff5]920    fig.suptitle(":".join(opts['name']) + extra_title)
[ba69383]921
[ca9e54e]922    if have_base and have_comp and opts['show_hist']:
[ba69383]923        plt.figure()
[346bc88]924        v = relerr
[caeb06d]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))
[ba69383]929        plt.ylabel('P(err)')
[ec7e360]930        plt.title('Distribution of relative error between calculation engines')
[ba69383]931
[013adb7]932    return limits
933
[0763009]934
[87985ca]935# ===========================================================================
936#
[bb39b4a]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
[3c24ccd]944    'plot', 'noplot', 'weights',
[b89f519]945    'linear', 'log', 'q4',
[bb39b4a]946    'rel', 'abs',
[5d316e9]947    'hist', 'nohist',
[bb39b4a]948    'title=',
949
950    # Data generation
[ced5bd2]951    'data=', 'noise=', 'res=', 'nq=', 'q=',
952    'lowq', 'midq', 'highq', 'exq', 'zero',
[bb39b4a]953    '2d', '1d',
954
955    # Parameter set
956    'preset', 'random', 'random=', 'sets=',
957    'demo', 'default',  # TODO: remove demo/default
958    'nopars', 'pars',
[e3571cb]959    'sphere', 'sphere=', # integrate over a sphere in 2d with n points
[bb39b4a]960
961    # Calculation options
962    'poly', 'mono', 'cutoff=',
963    'magnetic', 'nonmagnetic',
[ff31782]964    'accuracy=', 'ngauss=',
[765eb0e]965    'neval=',  # for timing...
[bb39b4a]966
967    # Precision options
[8698a0d]968    'engine=',
[bb39b4a]969    'half', 'fast', 'single', 'double', 'single!', 'double!', 'quad!',
970
971    # Output options
972    'help', 'html', 'edit',
[87985ca]973    ]
974
[e65c3ba]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('=')])()
[bb39b4a]977
978
[b32dafd]979def columnize(items, indent="", width=79):
[dd7fc12]980    # type: (List[str], str, int) -> str
[caeb06d]981    """
[1d4017a]982    Format a list of strings into columns.
983
984    Returns a string with carriage returns ready for printing.
[caeb06d]985    """
[b32dafd]986    column_width = max(len(w) for w in items) + 1
[7cf2cfd]987    num_columns = (width - len(indent)) // column_width
[b32dafd]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)]
[7cf2cfd]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
[98d6cfc]997def get_pars(model_info, use_demo=False):
[dd7fc12]998    # type: (ModelInfo, bool) -> ParameterSet
[caeb06d]999    """
1000    Extract demo parameters from the model definition.
1001    """
[ec7e360]1002    # Get the default values for the parameters
[c499331]1003    pars = {}
[6d6508e]1004    for p in model_info.parameters.call_parameters:
[c499331]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:
[b32dafd]1013                dict(("%s%d%s" % (p.id, k, ext), val)
1014                     for k in range(1, p.length+1))
[c499331]1015            else:
[b32dafd]1016                pars[p.id + ext] = val
[ec7e360]1017
1018    # Plug in values given in demo
[765eb0e]1019    if use_demo and model_info.demo:
[6d6508e]1020        pars.update(model_info.demo)
[373d1b6]1021    return pars
1022
[ff1fff5]1023INTEGER_RE = re.compile("^[+-]?[1-9][0-9]*$")
[110f69c]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)
[17bbadd]1030
[8c65a33]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 @ ? ^ ! ,
[bb39b4a]1039PAR_SPLIT = ','
[424fe00]1040def parse_opts(argv):
1041    # type: (List[str]) -> Dict[str, Any]
[caeb06d]1042    """
1043    Parse command line options.
1044    """
[fc0fcd0]1045    MODELS = core.list_models()
[424fe00]1046    flags = [arg for arg in argv
[caeb06d]1047             if arg.startswith('-')]
[424fe00]1048    values = [arg for arg in argv
[caeb06d]1049              if not arg.startswith('-') and '=' in arg]
[424fe00]1050    positional_args = [arg for arg in argv
[0bdddc2]1051                       if not arg.startswith('-') and '=' not in arg]
[d547f16]1052    models = "\n    ".join("%-15s"%v for v in MODELS)
[424fe00]1053    if len(positional_args) == 0:
[7cf2cfd]1054        print(USAGE)
[caeb06d]1055        print("\nAvailable models:")
[7cf2cfd]1056        print(columnize(MODELS, indent="  "))
[424fe00]1057        return None
[87985ca]1058
[ec7e360]1059    invalid = [o[1:] for o in flags
[216a9e1]1060               if o[1:] not in NAME_OPTIONS
[d15a908]1061               and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)]
[87985ca]1062    if invalid:
[9404dd3]1063        print("Invalid options: %s"%(", ".join(invalid)))
[424fe00]1064        return None
[87985ca]1065
[bb39b4a]1066    name = positional_args[-1]
[ec7e360]1067
[e65c3ba]1068    # pylint: disable=bad-whitespace,C0321
[ec7e360]1069    # Interpret the flags
1070    opts = {
1071        'plot'      : True,
1072        'view'      : 'log',
1073        'is2d'      : False,
[ced5bd2]1074        'qmin'      : None,
[ec7e360]1075        'qmax'      : 0.05,
1076        'nq'        : 128,
1077        'res'       : 0.0,
[bb39b4a]1078        'noise'     : 0.0,
[ec7e360]1079        'accuracy'  : 'Low',
[bb39b4a]1080        'cutoff'    : '0.0',
[ec7e360]1081        'seed'      : -1,  # default to preset
[630156b]1082        'mono'      : True,
[0b040de]1083        # Default to magnetic a magnetic moment is set on the command line
[b6f10d8]1084        'magnetic'  : False,
[ec7e360]1085        'show_pars' : False,
1086        'show_hist' : False,
1087        'rel_err'   : True,
1088        'explore'   : False,
[98d6cfc]1089        'use_demo'  : True,
[dd7fc12]1090        'zero'      : False,
[234c532]1091        'html'      : False,
[a0d75ce]1092        'title'     : None,
[630156b]1093        'datafile'  : None,
[d9ec8f9]1094        'sets'      : 0,
[bb39b4a]1095        'engine'    : 'default',
[e3571cb]1096        'count'     : '1',
[3c24ccd]1097        'show_weights' : False,
[e3571cb]1098        'sphere'    : 0,
[ff31782]1099        'ngauss'    : '0',
[ec7e360]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
[ce0b154]1112        elif arg == '-lowq':    opts['qmax'] = 0.05
[e78edc4]1113        elif arg == '-zero':    opts['zero'] = True
[ec7e360]1114        elif arg.startswith('-nq='):       opts['nq'] = int(arg[4:])
[ced5bd2]1115        elif arg.startswith('-q='):
1116            opts['qmin'], opts['qmax'] = [float(v) for v in arg[3:].split(':')]
[ec7e360]1117        elif arg.startswith('-res='):      opts['res'] = float(arg[5:])
[bb39b4a]1118        elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:])
[0bdddc2]1119        elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:])
[ec7e360]1120        elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:]
[bb39b4a]1121        elif arg.startswith('-cutoff='):   opts['cutoff'] = arg[8:]
[a769b54]1122        elif arg.startswith('-title='):    opts['title'] = arg[7:]
[630156b]1123        elif arg.startswith('-data='):     opts['datafile'] = arg[6:]
[8698a0d]1124        elif arg.startswith('-engine='):   opts['engine'] = arg[8:]
[e3571cb]1125        elif arg.startswith('-neval='):    opts['count'] = arg[7:]
[ff31782]1126        elif arg.startswith('-ngauss='):   opts['ngauss'] = arg[8:]
[31eea1f]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
[e3571cb]1133        elif arg.startswith('-sphere'):
1134            opts['sphere'] = int(arg[8:]) if len(arg) > 7 else 150
1135            opts['is2d'] = True
[ec7e360]1136        elif arg == '-preset':  opts['seed'] = -1
1137        elif arg == '-mono':    opts['mono'] = True
1138        elif arg == '-poly':    opts['mono'] = False
[0b040de]1139        elif arg == '-magnetic':       opts['magnetic'] = True
1140        elif arg == '-nonmagnetic':    opts['magnetic'] = False
[ec7e360]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
[bb39b4a]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!'
[ec7e360]1154        elif arg == '-edit':    opts['explore'] = True
[98d6cfc]1155        elif arg == '-demo':    opts['use_demo'] = True
[97d89af]1156        elif arg == '-default': opts['use_demo'] = False
[3c24ccd]1157        elif arg == '-weights': opts['show_weights'] = True
[234c532]1158        elif arg == '-html':    opts['html'] = True
[630156b]1159        elif arg == '-help':    opts['html'] = True
[e65c3ba]1160    # pylint: enable=bad-whitespace,C0321
[ec7e360]1161
[97d89af]1162    # Magnetism forces 2D for now
1163    if opts['magnetic']:
1164        opts['is2d'] = True
1165
[d9ec8f9]1166    # Force random if sets is used
1167    if opts['sets'] >= 1 and opts['seed'] < 0:
[0bdddc2]1168        opts['seed'] = np.random.randint(1000000)
[d9ec8f9]1169    if opts['sets'] == 0:
1170        opts['sets'] = 1
[0bdddc2]1171
[bb39b4a]1172    # Create the computational engines
[ced5bd2]1173    if opts['qmin'] is None:
1174        opts['qmin'] = 0.001*opts['qmax']
[bb39b4a]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)
[ff31782]1181
[bb39b4a]1182    if PAR_SPLIT in name:
1183        names = name.split(PAR_SPLIT, 2)
1184        comparison = True
[ff1fff5]1185    else:
[bb39b4a]1186        names = [name]*2
[ff1fff5]1187    try:
[bb39b4a]1188        model_info = [core.load_model_info(k) for k in names]
[ff1fff5]1189    except ImportError as exc:
1190        print(str(exc))
1191        print("Could not find model; use one of:\n    " + models)
1192        return None
[87985ca]1193
[ff31782]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
[bb39b4a]1200    if PAR_SPLIT in opts['engine']:
[e3571cb]1201        opts['engine'] = opts['engine'].split(PAR_SPLIT, 2)
[bb39b4a]1202        comparison = True
1203    else:
[e3571cb]1204        opts['engine'] = [opts['engine']]*2
[0bdddc2]1205
[e3571cb]1206    if PAR_SPLIT in opts['count']:
1207        opts['count'] = [int(k) for k in opts['count'].split(PAR_SPLIT, 2)]
[bb39b4a]1208        comparison = True
[0bdddc2]1209    else:
[e3571cb]1210        opts['count'] = [int(opts['count'])]*2
[bb39b4a]1211
1212    if PAR_SPLIT in opts['cutoff']:
[e3571cb]1213        opts['cutoff'] = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)]
[bb39b4a]1214        comparison = True
[0bdddc2]1215    else:
[e3571cb]1216        opts['cutoff'] = [float(opts['cutoff'])]*2
[bb39b4a]1217
[ff31782]1218    base = make_engine(model_info[0], data, opts['engine'][0],
1219                       opts['cutoff'][0], opts['ngauss'][0])
[bb39b4a]1220    if comparison:
[ff31782]1221        comp = make_engine(model_info[1], data, opts['engine'][1],
1222                           opts['cutoff'][1], opts['ngauss'][1])
[0bdddc2]1223    else:
1224        comp = None
1225
1226    # pylint: disable=bad-whitespace
1227    # Remember it all
1228    opts.update({
1229        'data'      : data,
[bb39b4a]1230        'name'      : names,
[e3571cb]1231        'info'      : model_info,
[0bdddc2]1232        'engines'   : [base, comp],
1233        'values'    : values,
1234    })
1235    # pylint: enable=bad-whitespace
1236
[e3571cb]1237    # Set the integration parameters to the half sphere
1238    if opts['sphere'] > 0:
1239        set_spherical_integration_parameters(opts, opts['sphere'])
1240
[0bdddc2]1241    return opts
1242
[e3571cb]1243def set_spherical_integration_parameters(opts, steps):
[110f69c]1244    # type: (Dict[str, Any], int) -> None
[e3571cb]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([
[31eea1f]1251        #'theta=90',
[e3571cb]1252        'theta_pd=%g'%(90/np.sqrt(3)),
1253        'theta_pd_n=%d'%steps,
1254        'theta_pd_type=rectangle',
[31eea1f]1255        #'phi=0',
[e3571cb]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:
[a5f91a7]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        ])
[e3571cb]1268
1269def parse_pars(opts, maxdim=np.inf):
[110f69c]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    """
[e3571cb]1281    model_info, model_info2 = opts['info']
[0bdddc2]1282
[ec7e360]1283    # Get demo parameters from model definition, or use default parameters
1284    # if model does not define demo parameters
[98d6cfc]1285    pars = get_pars(model_info, opts['use_demo'])
[ff1fff5]1286    pars2 = get_pars(model_info2, opts['use_demo'])
[248561a]1287    pars2.update((k, v) for k, v in pars.items() if k in pars2)
[ff1fff5]1288    # randomize parameters
1289    #pars.update(set_pars)  # set value before random to control range
1290    if opts['seed'] > -1:
[0bdddc2]1291        pars = randomize_pars(model_info, pars)
[e3571cb]1292        limit_dimensions(model_info, pars, maxdim)
[ff1fff5]1293        if model_info != model_info2:
[0bdddc2]1294            pars2 = randomize_pars(model_info2, pars2)
[376b0ee]1295            limit_dimensions(model_info2, pars2, maxdim)
[158cee4]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
[ff1fff5]1300        else:
1301            pars2 = pars.copy()
[158cee4]1302        constrain_pars(model_info, pars)
1303        constrain_pars(model_info2, pars2)
[97d89af]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'])
[87985ca]1308
1309    # Fill in parameters given on the command line
[ec7e360]1310    presets = {}
[ff1fff5]1311    presets2 = {}
[0bdddc2]1312    for arg in opts['values']:
[d15a908]1313        k, v = arg.split('=', 1)
[ff1fff5]1314        if k not in pars and k not in pars2:
[ec7e360]1315            # extract base name without polydispersity info
[87985ca]1316            s = set(p.split('_pd')[0] for p in pars)
[d15a908]1317            print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s))))
[424fe00]1318            return None
[110f69c]1319        v1, v2 = v.split(PAR_SPLIT, 2) if PAR_SPLIT in v else (v, v)
[ff1fff5]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
[b6f10d8]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
[ff1fff5]1333    # Evaluate preset parameter expressions
[a21d889]1334    # Note: need to replace ':' with '_' in parameter names and expressions
1335    # in order to support math on magnetic parameters.
[248561a]1336    context = MATH.copy()
[fe25eda]1337    context['np'] = np
[a21d889]1338    context.update((k.replace(':', '_'), v) for k, v in pars.items())
[0bdddc2]1339    context.update((k, v) for k, v in presets.items() if isinstance(v, float))
[a21d889]1340    #for k,v in sorted(context.items()): print(k, v)
[ff1fff5]1341    for k, v in presets.items():
1342        if not isinstance(v, float) and not k.endswith('_type'):
[a21d889]1343            presets[k] = eval(v.replace(':', '_'), context)
[ff1fff5]1344    context.update(presets)
[a21d889]1345    context.update((k.replace(':', '_'), v) for k, v in presets2.items() if isinstance(v, float))
[ff1fff5]1346    for k, v in presets2.items():
1347        if not isinstance(v, float) and not k.endswith('_type'):
[a21d889]1348            presets2[k] = eval(v.replace(':', '_'), context)
[ff1fff5]1349
1350    # update parameters with presets
[ec7e360]1351    pars.update(presets)  # set value after random to control value
[ff1fff5]1352    pars2.update(presets2)  # set value after random to control value
[fcd7bbd]1353    #import pprint; pprint.pprint(model_info)
[ff1fff5]1354
[aa25fc7]1355    # Hack to load user-defined distributions; run through all parameters
1356    # and make sure any pd_type parameter is a defined distribution.
1357    if (any(p.endswith('pd_type') and v not in weights.MODELS
1358            for p, v in pars.items()) or
1359        any(p.endswith('pd_type') and v not in weights.MODELS
1360            for p, v in pars2.items())):
1361       weights.load_weights()
1362
[ec7e360]1363    if opts['show_pars']:
[0bdddc2]1364        if model_info.name != model_info2.name or pars != pars2:
[248561a]1365            print("==== %s ====="%model_info.name)
1366            print(str(parlist(model_info, pars, opts['is2d'])))
1367            print("==== %s ====="%model_info2.name)
1368            print(str(parlist(model_info2, pars2, opts['is2d'])))
1369        else:
1370            print(str(parlist(model_info, pars, opts['is2d'])))
[ec7e360]1371
[0bdddc2]1372    return pars, pars2
[ec7e360]1373
[234c532]1374def show_docs(opts):
1375    # type: (Dict[str, Any]) -> None
1376    """
1377    show html docs for the model
1378    """
[c4e3215]1379    from .generate import make_html
1380    from . import rst2html
1381
[e3571cb]1382    info = opts['info'][0]
[c4e3215]1383    html = make_html(info)
1384    path = os.path.dirname(info.filename)
[110f69c]1385    url = "file://" + path.replace("\\", "/")[2:] + "/"
[1fbadb2]1386    rst2html.view_html_wxapp(html, url)
[234c532]1387
[ec7e360]1388def explore(opts):
[dd7fc12]1389    # type: (Dict[str, Any]) -> None
[d15a908]1390    """
[234c532]1391    explore the model using the bumps gui.
[d15a908]1392    """
[7ae2b7f]1393    import wx  # type: ignore
1394    from bumps.names import FitProblem  # type: ignore
1395    from bumps.gui.app_frame import AppFrame  # type: ignore
[ca9e54e]1396    from bumps.gui import signal
[ec7e360]1397
[d15a908]1398    is_mac = "cocoa" in wx.version()
[80013a6]1399    # Create an app if not running embedded
1400    app = wx.App() if wx.GetApp() is None else None
[ca9e54e]1401    model = Explore(opts)
1402    problem = FitProblem(model)
[0bdddc2]1403    frame = AppFrame(parent=None, title="explore", size=(1000, 700))
1404    if not is_mac:
1405        frame.Show()
[ec7e360]1406    frame.panel.set_model(model=problem)
1407    frame.panel.Layout()
1408    frame.panel.aui.Split(0, wx.TOP)
[110f69c]1409    def _reset_parameters(event):
[ca9e54e]1410        model.revert_values()
1411        signal.update_parameters(problem)
[110f69c]1412    frame.Bind(wx.EVT_TOOL, _reset_parameters, frame.ToolBar.GetToolByPos(1))
[e65c3ba]1413    if is_mac:
1414        frame.Show()
[80013a6]1415    # If running withing an app, start the main loop
[0bdddc2]1416    if app:
1417        app.MainLoop()
[ec7e360]1418
1419class Explore(object):
1420    """
[d15a908]1421    Bumps wrapper for a SAS model comparison.
1422
1423    The resulting object can be used as a Bumps fit problem so that
1424    parameters can be adjusted in the GUI, with plots updated on the fly.
[ec7e360]1425    """
1426    def __init__(self, opts):
[dd7fc12]1427        # type: (Dict[str, Any]) -> None
[7ae2b7f]1428        from bumps.cli import config_matplotlib  # type: ignore
[608e31e]1429        from . import bumps_model
[ec7e360]1430        config_matplotlib()
1431        self.opts = opts
[0bdddc2]1432        opts['pars'] = list(opts['pars'])
[ca9e54e]1433        p1, p2 = opts['pars']
[e3571cb]1434        m1, m2 = opts['info']
[ca9e54e]1435        self.fix_p2 = m1 != m2 or p1 != p2
1436        model_info = m1
1437        pars, pd_types = bumps_model.create_parameters(model_info, **p1)
[21b116f]1438        # Initialize parameter ranges, fixing the 2D parameters for 1D data.
[ec7e360]1439        if not opts['is2d']:
[85fe7f8]1440            for p in model_info.parameters.user_parameters({}, is2d=False):
[303d8d6]1441                for ext in ['', '_pd', '_pd_n', '_pd_nsigma']:
[69aa451]1442                    k = p.name+ext
[303d8d6]1443                    v = pars.get(k, None)
1444                    if v is not None:
1445                        v.range(*parameter_range(k, v.value))
[ec7e360]1446        else:
[013adb7]1447            for k, v in pars.items():
[ec7e360]1448                v.range(*parameter_range(k, v.value))
1449
1450        self.pars = pars
[ca9e54e]1451        self.starting_values = dict((k, v.value) for k, v in pars.items())
[ec7e360]1452        self.pd_types = pd_types
[fbb9397]1453        self.limits = None
[ec7e360]1454
[ca9e54e]1455    def revert_values(self):
[110f69c]1456        # type: () -> None
1457        """
1458        Restore starting values of the parameters.
1459        """
[ca9e54e]1460        for k, v in self.starting_values.items():
1461            self.pars[k].value = v
1462
1463    def model_update(self):
[110f69c]1464        # type: () -> None
1465        """
1466        Respond to signal that model parameters have been changed.
1467        """
[ca9e54e]1468        pass
1469
[ec7e360]1470    def numpoints(self):
[dd7fc12]1471        # type: () -> int
[ec7e360]1472        """
[608e31e]1473        Return the number of points.
[ec7e360]1474        """
1475        return len(self.pars) + 1  # so dof is 1
1476
1477    def parameters(self):
[dd7fc12]1478        # type: () -> Any   # Dict/List hierarchy of parameters
[ec7e360]1479        """
[608e31e]1480        Return a dictionary of parameters.
[ec7e360]1481        """
1482        return self.pars
1483
1484    def nllf(self):
[dd7fc12]1485        # type: () -> float
[608e31e]1486        """
1487        Return cost.
1488        """
[d15a908]1489        # pylint: disable=no-self-use
[ec7e360]1490        return 0.  # No nllf
1491
1492    def plot(self, view='log'):
[dd7fc12]1493        # type: (str) -> None
[ec7e360]1494        """
1495        Plot the data and residuals.
1496        """
[608e31e]1497        pars = dict((k, v.value) for k, v in self.pars.items())
[ec7e360]1498        pars.update(self.pd_types)
[ff1fff5]1499        self.opts['pars'][0] = pars
[ca9e54e]1500        if not self.fix_p2:
1501            self.opts['pars'][1] = pars
1502        result = run_models(self.opts)
1503        limits = plot_models(self.opts, result, limits=self.limits)
[013adb7]1504        if self.limits is None:
1505            vmin, vmax = limits
[dd7fc12]1506            self.limits = vmax*1e-7, 1.3*vmax
[d86f0fc]1507            import pylab
1508            pylab.clf()
[ca9e54e]1509            plot_models(self.opts, result, limits=self.limits)
[87985ca]1510
1511
[424fe00]1512def main(*argv):
1513    # type: (*str) -> None
[d15a908]1514    """
1515    Main program.
1516    """
[424fe00]1517    opts = parse_opts(argv)
1518    if opts is not None:
[48462b0]1519        if opts['seed'] > -1:
1520            print("Randomize using -random=%i"%opts['seed'])
1521            np.random.seed(opts['seed'])
[234c532]1522        if opts['html']:
1523            show_docs(opts)
1524        elif opts['explore']:
[0bdddc2]1525            opts['pars'] = parse_pars(opts)
[8f04da4]1526            if opts['pars'] is None:
1527                return
[424fe00]1528            explore(opts)
1529        else:
1530            compare(opts)
[d15a908]1531
[8a20be5]1532if __name__ == "__main__":
[424fe00]1533    main(*sys.argv[1:])
Note: See TracBrowser for help on using the repository browser.