1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
|
from __future__ import print_function
import sys
import math
from optparse import OptionParser
try:
import matplotlib
matplotlib.use('WXAgg')
import pylab
HASPYLAB = True
except ImportError:
HASPYLAB = False
if 'nose' in arg:
HASPYLAB = False
from lmfit import Parameters, minimize
from NISTModels import Models, ReadNistData
def ndig(a, b):
"precision for NIST values"
return round(-math.log10((abs(abs(a)-abs(b)) +1.e-15)/ abs(b)))
def Compare_NIST_Results(DataSet, myfit, params, NISTdata):
print(' ======================================')
print(' %s: ' % DataSet)
print(' | Parameter Name | Value Found | Certified Value | # Matching Digits |')
print(' |----------------+----------------+------------------+-------------------|')
val_dig_min = 200
err_dig_min = 200
for i in range(NISTdata['nparams']):
parname = 'b%i' % (i+1)
par = params[parname]
thisval = par.value
certval = NISTdata['cert_values'][i]
vdig = ndig(thisval, certval)
pname = (parname + ' value ' + ' '*14)[:14]
print(' | %s | % -.7e | % -.7e | %2i |' % (pname, thisval, certval, vdig))
val_dig_min = min(val_dig_min, vdig)
thiserr = par.stderr
certerr = NISTdata['cert_stderr'][i]
if thiserr is not None and myfit.errorbars:
edig = ndig(thiserr, certerr)
ename = (parname + ' stderr' + ' '*14)[:14]
print(' | %s | % -.7e | % -.7e | %2i |' % (ename, thiserr, certerr, edig))
err_dig_min = min(err_dig_min, edig)
print(' |----------------+----------------+------------------+-------------------|')
sumsq = NISTdata['sum_squares']
try:
chi2 = myfit.chisqr
print(' | Sum of Squares | %.7e | %.7e | %2i |' % (chi2, sumsq,
ndig(chi2, sumsq)))
except:
pass
print(' |----------------+----------------+------------------+-------------------|')
if not myfit.errorbars:
print(' | * * * * COULD NOT ESTIMATE UNCERTAINTIES * * * * |')
err_dig_min = 0
if err_dig_min < 199:
print(' Worst agreement: %i digits for value, %i digits for error ' % (val_dig_min, err_dig_min))
else:
print(' Worst agreement: %i digits' % (val_dig_min))
return val_dig_min
def NIST_Test(DataSet, method='leastsq', start='start2', plot=True):
NISTdata = ReadNistData(DataSet)
resid, npar, dimx = Models[DataSet]
y = NISTdata['y']
x = NISTdata['x']
params = Parameters()
for i in range(npar):
pname = 'b%i' % (i+1)
cval = NISTdata['cert_values'][i]
cerr = NISTdata['cert_stderr'][i]
pval1 = NISTdata[start][i]
params.add(pname, value=pval1)
myfit = minimize(resid, params, method=method, args=(x,), kws={'y':y})
digs = Compare_NIST_Results(DataSet, myfit, params, NISTdata)
if plot and HASPYLAB:
fit = -resid(params, x, )
pylab.plot(x, y, 'ro')
pylab.plot(x, fit, 'k+-')
pylab.show()
return digs > 2
modelnames = []
ms = ''
for d in sorted(Models.keys()):
ms = ms + ' %s ' % d
if len(ms) > 55:
modelnames.append(ms)
ms = ' '
modelnames.append(ms)
modelnames = '\n'.join(modelnames)
usage = """
=== Test Fit to NIST StRD Models ===
usage:
------
python fit_NIST.py [options] Model Start
where Start is one of 'start1','start2' or 'cert', for different
starting values, and Model is one of
%s
if Model = 'all', all models and starting values will be run.
options:
--------
-m name of fitting method. One of:
leastsq, nelder, powell, lbfgsb, bfgs,
tnc, cobyla, slsqp, cg, newto-cg
leastsq (Levenberg-Marquardt) is the default
""" % modelnames
############################
parser = OptionParser(usage=usage, prog="fit-NIST.py")
parser.add_option("-m", "--method", dest="method", metavar='METH',
default='leastsq', help="set method name, default = 'leastsq'")
(opts, args) = parser.parse_args()
dset = ''
start = 'start2'
if len(args) > 0:
dset = args[0]
if len(args) > 1:
start = args[1]
if dset.lower() == 'all':
tpass = 0
tfail = 0
failures = []
dsets = sorted(Models.keys())
for dset in dsets:
for start in ('start1', 'start2', 'cert'):
if NIST_Test(dset, method=opts.method, start=start, plot=False):
tpass += 1
else:
tfail += 1
failures.append(" %s (starting at '%s')" % (dset, start))
print('--------------------------------------')
print(' Fit Method: %s ' % opts.method)
print(' Final Results: %i pass, %i fail.' % (tpass, tfail))
print(' Tests Failed for:\n %s' % '\n '.join(failures))
print('--------------------------------------')
elif dset not in Models:
print(usage)
else:
NIST_Test(dset, method=opts.method, start=start, plot=True)
|
