source: sasmodels/sasmodels/compare_many.py @ 40a87fa

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 40a87fa was 40a87fa, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

lint and latex cleanup

  • Property mode set to 100755
File size: 8.9 KB
Line 
1#!/usr/bin/env python
2"""
3Program to compare results from many random parameter sets for a given model.
4
5The result is a comma separated value (CSV) table that can be redirected
6from standard output into a file and loaded into a spreadsheet.
7
8The models are compared for each parameter set and if the difference is
9greater than expected for that precision, the parameter set is labeled
10as bad and written to the output, along with the random seed used to
11generate that parameter value.  This seed can be used with :mod:`compare`
12to reload and display the details of the model.
13"""
14from __future__ import print_function
15
16import sys
17import traceback
18
19import numpy as np  # type: ignore
20
21from . import core
22from .compare import (randomize_pars, suppress_pd, make_data,
23                      make_engine, get_pars, columnize,
24                      constrain_pars, constrain_new_to_old)
25
26MODELS = core.list_models()
27
28def calc_stats(target, value, index):
29    """
30    Calculate statistics between the target value and the computed value.
31
32    *target* and *value* are the vectors being compared, with the
33    difference normalized by *target* to get relative error.  Only
34    the elements listed in *index* are used, though index may be
35    and empty slice defined by *slice(None, None)*.
36
37    Returns:
38
39        *maxrel* the maximum relative difference
40
41        *rel95* the relative difference with the 5% biggest differences ignored
42
43        *maxabs* the maximum absolute difference for the 5% biggest differences
44
45        *maxval* the maximum value for the 5% biggest differences
46    """
47    resid = abs(value-target)[index]
48    relerr = resid/target[index]
49    sorted_rel_index = np.argsort(relerr)
50    #p90 = int(len(relerr)*0.90)
51    p95 = int(len(relerr)*0.95)
52    maxrel = np.max(relerr)
53    rel95 = relerr[sorted_rel_index[p95]]
54    maxabs = np.max(resid[sorted_rel_index[p95:]])
55    maxval = np.max(value[sorted_rel_index[p95:]])
56    return maxrel, rel95, maxabs, maxval
57
58def print_column_headers(pars, parts):
59    """
60    Generate column headers for the differences and for the parameters,
61    and print them to standard output.
62    """
63    stats = list('Max rel err|95% rel err|Max abs err above 90% rel|Max value above 90% rel'.split('|'))
64    groups = ['']
65    for p in parts:
66        groups.append(p)
67        groups.extend(['']*(len(stats)-1))
68    groups.append("Parameters")
69    columns = ['Seed'] + stats*len(parts) +  list(sorted(pars.keys()))
70    print(','.join('"%s"'%c for c in groups))
71    print(','.join('"%s"'%c for c in columns))
72
73# Target 'good' value for various precision levels.
74PRECISION = {
75    'fast': 1e-3,
76    'half': 1e-3,
77    'single': 5e-5,
78    'double': 5e-14,
79    'single!': 5e-5,
80    'double!': 5e-14,
81    'quad!': 5e-18,
82    'sasview': 5e-14,
83}
84def compare_instance(name, data, index, N=1, mono=True, cutoff=1e-5,
85                     base='sasview', comp='double'):
86    r"""
87    Compare the model under different calculation engines.
88
89    *name* is the name of the model.
90
91    *data* is the data object giving $q, \Delta q$ calculation points.
92
93    *index* is the active set of points.
94
95    *N* is the number of comparisons to make.
96
97    *cutoff* is the polydispersity weight cutoff to make the calculation
98    a little bit faster.
99
100    *base* and *comp* are the names of the calculation engines to compare.
101    """
102
103    is_2d = hasattr(data, 'qx_data')
104    model_info = core.load_model_info(name)
105    pars = get_pars(model_info, use_demo=True)
106    header = ('\n"Model","%s","Count","%d","Dimension","%s"'
107              % (name, N, "2D" if is_2d else "1D"))
108    if not mono: header += ',"Cutoff",%g'%(cutoff,)
109    print(header)
110
111    if is_2d:
112        if not model_info.parameters.has_2d:
113            print(',"1-D only"')
114            return
115
116    # Some not very clean macros for evaluating the models and checking the
117    # results.  They freely use variables from the current scope, even some
118    # which have not been defined yet, complete with abuse of mutable lists
119    # to allow them to update values in the current scope since nonlocal
120    # declarations are not available in python 2.7.
121    def try_model(fn, pars):
122        """
123        Return the model evaluated at *pars*.  If there is an exception,
124        print it and return NaN of the right shape.
125        """
126        try:
127            result = fn(**pars)
128        except Exception:
129            traceback.print_exc()
130            print("when comparing %s for %d"%(name, seed))
131            if hasattr(data, 'qx_data'):
132                result = np.NaN*data.data
133            else:
134                result = np.NaN*data.x
135        return result
136    def check_model(pars):
137        """
138        Run the two calculators against *pars*, returning statistics
139        on the differences.  See :func:`calc_stats` for the list of stats.
140        """
141        base_value = try_model(calc_base, pars)
142        comp_value = try_model(calc_comp, pars)
143        stats = calc_stats(base_value, comp_value, index)
144        max_diff[0] = max(max_diff[0], stats[0])
145        good[0] = good[0] and (stats[0] < expected)
146        return list(stats)
147
148
149    try:
150        calc_base = make_engine(model_info, data, base, cutoff)
151        calc_comp = make_engine(model_info, data, comp, cutoff)
152    except Exception as exc:
153        #raise
154        print('"Error: %s"'%str(exc).replace('"', "'"))
155        print('"good","%d of %d","max diff",%g' % (0, N, np.NaN))
156        return
157    expected = max(PRECISION[base], PRECISION[comp])
158
159    num_good = 0
160    first = True
161    max_diff = [0]
162    for k in range(N):
163        print("%s %d"%(name, k), file=sys.stderr)
164        seed = np.random.randint(1e6)
165        pars_i = randomize_pars(model_info, pars, seed)
166        constrain_pars(model_info, pars_i)
167        constrain_new_to_old(model_info, pars_i)
168        if mono:
169            pars_i = suppress_pd(pars_i)
170
171        good = [True]
172        columns = check_model(pars_i)
173        columns += [v for _, v in sorted(pars_i.items())]
174        if first:
175            labels = [" vs. ".join((calc_base.engine, calc_comp.engine))]
176            print_column_headers(pars_i, labels)
177            first = False
178        if good[0]:
179            num_good += 1
180        else:
181            print(("%d,"%seed)+','.join("%s"%v for v in columns))
182    print('"good","%d of %d","max diff",%g'%(num_good, N, max_diff[0]))
183
184
185def print_usage():
186    """
187    Print the command usage string.
188    """
189    print("usage: compare_many.py MODEL COUNT (1dNQ|2dNQ) (CUTOFF|mono) (single|double|quad)")
190
191
192def print_models():
193    """
194    Print the list of available models in columns.
195    """
196    print(columnize(MODELS, indent="  "))
197
198
199def print_help():
200    """
201    Print usage string, the option description and the list of available models.
202    """
203    print_usage()
204    print("""\
205
206MODEL is the model name of the model or "all" for all the models
207in alphabetical order.
208
209COUNT is the number of randomly generated parameter sets to try. A value
210of "10000" is a reasonable check for monodisperse models, or "100" for
211polydisperse models.   For a quick check, use "100" and "5" respectively.
212
213NQ is the number of Q values to calculate.  If it starts with "1d", then
214it is a 1-dimensional problem, with log spaced Q points from 1e-3 to 1.0.
215If it starts with "2d" then it is a 2-dimensional problem, with linearly
216spaced points Q points from -1.0 to 1.0 in each dimension. The usual
217values are "1d100" for 1-D and "2d32" for 2-D.
218
219CUTOFF is the cutoff value to use for the polydisperse distribution. Weights
220below the cutoff will be ignored.  Use "mono" for monodisperse models.  The
221choice of polydisperse parameters, and the number of points in the distribution
222is set in compare.py defaults for each model.
223
224PRECISION is the floating point precision to use for comparisons.  If two
225precisions are given, then compare one to the other, ignoring sasview.
226
227Available models:
228""")
229    print_models()
230
231def main():
232    """
233    Main program.
234    """
235    if len(sys.argv) not in (6, 7):
236        print_help()
237        sys.exit(1)
238
239    model = sys.argv[1]
240    if not (model in MODELS) and (model != "all"):
241        print('Bad model %s.  Use "all" or one of:'%model)
242        print_models()
243        sys.exit(1)
244    try:
245        count = int(sys.argv[2])
246        is2D = sys.argv[3].startswith('2d')
247        assert sys.argv[3][1] == 'd'
248        Nq = int(sys.argv[3][2:])
249        mono = sys.argv[4] == 'mono'
250        cutoff = float(sys.argv[4]) if not mono else 0
251        base = sys.argv[5]
252        comp = sys.argv[6] if len(sys.argv) > 6 else "sasview"
253    except Exception:
254        traceback.print_exc()
255        print_usage()
256        sys.exit(1)
257
258    data, index = make_data({'qmax':1.0, 'is2d':is2D, 'nq':Nq, 'res':0.,
259                             'accuracy': 'Low', 'view':'log', 'zero': False})
260    model_list = [model] if model != "all" else MODELS
261    for model in model_list:
262        compare_instance(model, data, index, N=count, mono=mono,
263                         cutoff=cutoff, base=base, comp=comp)
264
265if __name__ == "__main__":
266    #from .compare import push_seed
267    #with push_seed(1): main()
268    main()
Note: See TracBrowser for help on using the repository browser.