Coverage for larch/xafs/autobk.py: 49%

324 statements  

« prev     ^ index     » next       coverage.py v7.6.0, created at 2024-10-16 21:04 +0000

1#!/usr/bin/env python 

2import sys 

3import time 

4import numpy as np 

5from scipy.interpolate import splrep, splev, UnivariateSpline 

6from scipy.stats import t 

7from scipy.special import erf 

8from scipy.optimize import leastsq 

9from lmfit import Parameter, Parameters, minimize, fit_report 

10 

11import uncertainties 

12 

13from larch import (Group, Make_CallArgs, parse_group_args, isgroup) 

14from larch.math import index_of, index_nearest, realimag, remove_dups 

15 

16from .xafsutils import ETOK, TINY_ENERGY, set_xafsGroup 

17from .xafsft import ftwindow, xftf_fast 

18from .pre_edge import find_e0, pre_edge 

19 

20FMT_COEF = 'coef_%2.2i' 

21NFEV = 0 

22def spline_eval(kraw, mu, knots, coefs, order, kout): 

23 """eval bkg(kraw) and chi(k) for knots, coefs, order""" 

24 bkg = splev(kraw, [knots, coefs, order]) 

25 chi = UnivariateSpline(kraw, (mu-bkg), s=0)(kout) 

26 return bkg, chi 

27 

28def _resid(vcoefs, ncoef, kraw, mu, chi_std, knots, order, kout, 

29 ftwin, nfft, irbkg, nclamp, clamp_lo, clamp_hi): 

30 global NFEV 

31 NFEV += 1 

32 nspl = len(vcoefs) 

33 coefs = np.ones(ncoef)*vcoefs[-1] 

34 coefs[:nspl] = vcoefs 

35 bkg, chi = spline_eval(kraw, mu, knots, coefs, order, kout) 

36 if chi_std is not None: 

37 chi = chi - chi_std 

38 out = realimag(xftf_fast(chi*ftwin, nfft=nfft)[:irbkg]) 

39 if nclamp == 0: 

40 return out 

41 scale = 1.0 + 100*(out*out).mean() 

42 return np.concatenate((out, 

43 abs(clamp_lo)*scale*chi[:nclamp], 

44 abs(clamp_hi)*scale*chi[-nclamp:])) 

45 

46 

47@Make_CallArgs(["energy" ,"mu"]) 

48def autobk(energy, mu=None, group=None, rbkg=1, nknots=None, e0=None, ek0=None, 

49 edge_step=None, kmin=0, kmax=None, kweight=1, dk=0.1, 

50 win='hanning', k_std=None, chi_std=None, nfft=2048, kstep=0.05, 

51 pre_edge_kws=None, nclamp=3, clamp_lo=0, clamp_hi=1, 

52 calc_uncertainties=False, err_sigma=1, _larch=None, **kws): 

53 """Use Autobk algorithm to remove XAFS background 

54 

55 Parameters: 

56 ----------- 

57 energy: 1-d array of x-ray energies, in eV, or group 

58 mu: 1-d array of mu(E) 

59 group: output group (and input group for e0 and edge_step). 

60 rbkg: distance (in Ang) for chi(R) above 

61 which the signal is ignored. Default = 1. 

62 e0: edge energy, in eV. (deprecated: use ek0) 

63 ek0: edge energy, in eV. If None, it will be determined. 

64 edge_step: edge step. If None, it will be determined. 

65 pre_edge_kws: keyword arguments to pass to pre_edge() 

66 nknots: number of knots in spline. If None, it will be determined. 

67 kmin: minimum k value [0] 

68 kmax: maximum k value [full data range]. 

69 kweight: k weight for FFT. [1] 

70 dk: FFT window window parameter. [0.1] 

71 win: FFT window function name. ['hanning'] 

72 nfft: array size to use for FFT [2048] 

73 kstep: k step size to use for FFT [0.05] 

74 k_std: optional k array for standard chi(k). 

75 chi_std: optional chi array for standard chi(k). 

76 nclamp: number of energy end-points for clamp [3] 

77 clamp_lo: weight of low-energy clamp [0] 

78 clamp_hi: weight of high-energy clamp [1] 

79 calc_uncertaintites: Flag to calculate uncertainties in 

80 mu_0(E) and chi(k) [True] 

81 err_sigma: sigma level for uncertainties in mu_0(E) and chi(k) [1] 

82 

83 Output arrays are written to the provided group. 

84 

85 Follows the 'First Argument Group' convention. 

86 """ 

87 msg = sys.stdout.write 

88 if _larch is not None: 

89 msg = _larch.writer.write 

90 if 'kw' in kws: 

91 kweight = kws.pop('kw') 

92 if len(kws) > 0: 

93 msg('Unrecognized arguments for autobk():\n') 

94 msg(' %s\n' % (', '.join(kws.keys()))) 

95 return 

96 energy, mu, group = parse_group_args(energy, members=('energy', 'mu'), 

97 defaults=(mu,), group=group, 

98 fcn_name='autobk') 

99 if len(energy.shape) > 1: 

100 energy = energy.squeeze() 

101 if len(mu.shape) > 1: 

102 mu = mu.squeeze() 

103 energy = remove_dups(energy, tiny=TINY_ENERGY) 

104 # if e0 or edge_step are not specified, get them, either from the 

105 # passed-in group or from running pre_edge() 

106 group = set_xafsGroup(group, _larch=_larch) 

107 

108 if edge_step is None and isgroup(group, 'edge_step'): 

109 edge_step = group.edge_step 

110 if e0 is not None and ek0 is None: # command-line e0 still valid 

111 ek0 = e0 

112 if ek0 is None and isgroup(group, 'ek0'): 

113 ek0 = group.ek0 

114 if ek0 is None and isgroup(group, 'e0'): 

115 ek0 = group.e0 

116 

117 if ek0 is not None and (ek0 < energy.min() or ek0 > energy.max()): 

118 ek0 = None 

119 if ek0 is None or edge_step is None: 

120 # need to run pre_edge: 

121 pre_kws = dict(nnorm=None, nvict=0, pre1=None, 

122 pre2=None, norm1=None, norm2=None) 

123 if pre_edge_kws is not None: 

124 pre_kws.update(pre_edge_kws) 

125 pre_edge(energy, mu, group=group, _larch=_larch, **pre_kws) 

126 if ek0 is None: 

127 ek0 = group.e0 

128 if edge_step is None: 

129 edge_step = group.edge_step 

130 if ek0 is None or edge_step is None: 

131 msg('autobk() could not determine ek0 or edge_step!: trying running pre_edge first\n') 

132 return 

133 

134 # get array indices for rkbg and ek0: irbkg, iek0 

135 iek0 = index_of(energy, ek0) 

136 rgrid = np.pi/(kstep*nfft) 

137 rbkg = max(rbkg, 2*rgrid) 

138 

139 # save ungridded k (kraw) and grided k (kout) 

140 # and ftwin (*k-weighting) for FT in residual 

141 enpe = energy[iek0:] - ek0 

142 kraw = np.sign(enpe)*np.sqrt(ETOK*abs(enpe)) 

143 if kmax is None: 

144 kmax = max(kraw) 

145 else: 

146 kmax = max(0, min(max(kraw), kmax)) 

147 kout = kstep * np.arange(int(1.01+kmax/kstep), dtype='float64') 

148 iemax = min(len(energy), 2+index_of(energy, ek0+kmax*kmax/ETOK)) - 1 

149 

150 # interpolate provided chi(k) onto the kout grid 

151 if chi_std is not None and k_std is not None: 

152 chi_std = np.interp(kout, k_std, chi_std) 

153 # pre-load FT window 

154 ftwin = kout**kweight * ftwindow(kout, xmin=kmin, xmax=kmax, 

155 window=win, dx=dk, dx2=dk) 

156 # calc k-value and initial guess for y-values of spline params 

157 nspl = 1 + int(2*rbkg*(kmax-kmin)/np.pi) 

158 irbkg = int(1 + (nspl-1)*np.pi/(2*rgrid*(kmax-kmin))) 

159 if nknots is not None: 

160 nspl = nknots 

161 nspl = max(5, min(128, nspl)) 

162 spl_y, spl_k = np.ones(nspl), np.zeros(nspl) 

163 

164 for i in range(nspl): 

165 q = kmin + i*(kmax-kmin)/(nspl - 1) 

166 ik = index_nearest(kraw, q) 

167 i1 = min(len(kraw)-1, ik + 5) 

168 i2 = max(0, ik - 5) 

169 spl_k[i] = kraw[ik] 

170 spl_y[i] = (2*mu[ik+iek0] + mu[i1+iek0] + mu[i2+iek0] ) / 4.0 

171 

172 order = 3 

173 qmin, qmax = spl_k[0], spl_k[nspl-1] 

174 knots = [spl_k[0] - 1.e-4*(order-i) for i in range(order)] 

175 

176 for i in range(order, nspl): 

177 knots.append((i-order)*(qmax - qmin)/(nspl-order+1)) 

178 qlast = knots[-1] 

179 for i in range(order+1): 

180 knots.append(qlast + 1.e-4*(i+1)) 

181 

182 # coefs = [mu[index_nearest(energy, ek0 + q**2/ETOK)] for q in knots] 

183 knots, coefs, order = splrep(spl_k, spl_y, k=order) 

184 coefs[nspl:] = coefs[nspl-1] 

185 ncoefs = len(coefs) 

186 kraw_ = kraw[:iemax-iek0+1] 

187 mu_ = mu[iek0:iemax+1] 

188 initbkg, initchi = spline_eval(kraw_, mu_, knots, coefs, order, kout) 

189 global NFEV 

190 NFEV = 0 

191 

192 vcoefs = 1.0*coefs[:nspl] 

193 userargs = (len(coefs), kraw_, mu_, chi_std, knots, order, kout, 

194 ftwin, nfft, irbkg, nclamp, clamp_lo, clamp_hi) 

195 

196 lsout = leastsq(_resid, vcoefs, userargs, maxfev=2000*(ncoefs+1), 

197 gtol=0.0, ftol=1.e-6, xtol=1.e-6, epsfcn=1.e-6, 

198 full_output=1, col_deriv=0, factor=100, diag=None) 

199 

200 best, covar, _infodict, errmsg, ier = lsout 

201 final_coefs = coefs[:] 

202 final_coefs[:nspl] = best[:] 

203 final_coefs[nspl:] = best[-1] 

204 

205 chisqr = ((_resid(best, *userargs))**2).sum() 

206 redchi = chisqr / (2*irbkg+2*nclamp - nspl) 

207 

208 coefs_std = np.array([np.sqrt(redchi*covar[i, i]) for i in range(nspl)]) 

209 bkg, chi = spline_eval(kraw[:iemax-iek0+1], mu[iek0:iemax+1], 

210 knots, final_coefs, order, kout) 

211 obkg = mu[:]*1.0 

212 obkg[iek0:iek0+len(bkg)] = bkg 

213 

214 # outputs to group 

215 group = set_xafsGroup(group, _larch=_larch) 

216 group.bkg = obkg 

217 group.chie = (mu-obkg)/edge_step 

218 group.k = kout 

219 group.chi = chi/edge_step 

220 group.ek0 = ek0 

221 group.rbkg = rbkg 

222 

223 knots_y = np.array([coefs[i] for i in range(nspl)]) 

224 init_bkg = mu[:]*1.0 

225 init_bkg[iek0:iek0+len(bkg)] = initbkg 

226 # now fill in 'autobk_details' group 

227 

228 group.autobk_details = Group(kmin=kmin, kmax=kmax, irbkg=irbkg, 

229 nknots=len(spl_k), knots=knots, order=order, 

230 init_knots_y=spl_y, nspl=nspl, 

231 init_chi=initchi/edge_step, coefs=final_coefs, 

232 coefs_std=coefs_std, iek0=iek0, iemax=iemax, 

233 ek0=ek0, covar=covar, chisqr=chisqr, 

234 redchi=redchi, init_bkg=init_bkg, 

235 knots_y=knots_y, kraw=kraw, mu=mu) 

236 

237 if calc_uncertainties and covar is not None: 

238 autobk_delta_chi(group, err_sigma=err_sigma) 

239 

240 

241def autobk_delta_chi(group, err_sigma=1): 

242 """calculate uncertainties in chi(k) and bkg(E) 

243 after running autobk 

244 """ 

245 d = getattr(group, 'autobk_details', None) 

246 if d is None or getattr(d, 'covar', None) is None: 

247 return 

248 

249 nchi = len(group.chi) 

250 nmue = d.iemax-d.iek0 + 1 

251 nspl = d.nspl 

252 jac_chi = np.zeros(nchi*nspl).reshape((nspl, nchi)) 

253 jac_bkg = np.zeros(nmue*nspl).reshape((nspl, nmue)) 

254 tcoefs = np.ones(len(d.coefs)) * d.coefs[-1] 

255 

256 step = 0.5 

257 # find derivatives by hand 

258 for i in range(nspl): 

259 b = [0, 0] 

260 c = [0, 0] 

261 for k in (0, 1): 

262 tcoefs = [1.0*d.coefs[j] for j in range(nspl)] 

263 tcoefs[i] = d.coefs[i] + (2*k-1)*step*d.coefs_std[i] 

264 b[k], c[k] = spline_eval(d.kraw[:d.iemax-d.iek0+1], 

265 d.mu[d.iek0:d.iemax+1], 

266 d.knots, tcoefs, d.order, group.k) 

267 jac_chi[i] = (c[1]- c[0])/(2*step*d.coefs_std[i]) 

268 jac_bkg[i] = (b[1]- b[0])/(2*step*d.coefs_std[i]) 

269 

270 dfchi = np.zeros(nchi) 

271 dfbkg = np.zeros(nmue) 

272 for i in range(nspl): 

273 for j in range(nspl): 

274 dfchi += jac_chi[i]*jac_chi[j]*d.covar[i,j] 

275 dfbkg += jac_bkg[i]*jac_bkg[j]*d.covar[i,j] 

276 

277 prob = 0.5*(1.0 + erf(err_sigma/np.sqrt(2.0))) 

278 dchi = t.ppf(prob, nchi-nspl) * np.sqrt(dfchi*d.redchi) 

279 dbkg = t.ppf(prob, nmue-nspl) * np.sqrt(dfbkg*d.redchi) 

280 

281 if not any(np.isnan(dchi)): 

282 group.delta_chi = dchi 

283 group.delta_bkg = 0.0*d.mu 

284 group.delta_bkg[d.iek0:d.iek0+len(dbkg)] = dbkg 

285 

286 

287## version of autobk using lmfit 

288 

289def _lmfit_resid(pars, ncoefs=1, knots=None, order=3, irbkg=1, nfft=2048, 

290 kraw=None, mu=None, kout=None, ftwin=1, kweight=1, chi_std=None, 

291 nclamp=0, clamp_lo=1, clamp_hi=1, **kws): 

292 # coefs = [getattr(pars, FMT_COEF % i) for i in range(ncoefs)] 

293 coefs = [pars[FMT_COEF % i].value for i in range(ncoefs)] 

294 bkg, chi = spline_eval(kraw, mu, knots, coefs, order, kout) 

295 if chi_std is not None: 

296 chi = chi - chi_std 

297 out = realimag(xftf_fast(chi*ftwin, nfft=nfft)[:irbkg]) 

298 if nclamp == 0: 

299 return out 

300 # spline clamps: 

301 scale = 1.0 + 100*(out*out).mean() 

302 return np.concatenate((out, 

303 abs(clamp_lo)*scale*chi[:nclamp], 

304 abs(clamp_hi)*scale*chi[-nclamp:])) 

305 

306 

307 

308@Make_CallArgs(["energy" ,"mu"]) 

309def autobk_lmfit(energy, mu=None, group=None, rbkg=1, nknots=None, e0=None, ek0=None, 

310 edge_step=None, kmin=0, kmax=None, kweight=1, dk=0.1, 

311 win='hanning', k_std=None, chi_std=None, nfft=2048, kstep=0.05, 

312 pre_edge_kws=None, nclamp=3, clamp_lo=0, clamp_hi=1, 

313 calc_uncertainties=True, err_sigma=1, _larch=None, **kws): 

314 """Use Autobk algorithm to remove XAFS background 

315 

316 Parameters: 

317 ----------- 

318 energy: 1-d array of x-ray energies, in eV, or group 

319 mu: 1-d array of mu(E) 

320 group: output group (and input group for e0 and edge_step). 

321 rbkg: distance (in Ang) for chi(R) above 

322 which the signal is ignored. Default = 1. 

323 e0: edge energy, in eV. (deprecated: use ek0) 

324 ek0: edge energy, in eV. If None, it will be determined. 

325 edge_step: edge step. If None, it will be determined. 

326 pre_edge_kws: keyword arguments to pass to pre_edge() 

327 nknots: number of knots in spline. If None, it will be determined. 

328 kmin: minimum k value [0] 

329 kmax: maximum k value [full data range]. 

330 kweight: k weight for FFT. [1] 

331 dk: FFT window window parameter. [0.1] 

332 win: FFT window function name. ['hanning'] 

333 nfft: array size to use for FFT [2048] 

334 kstep: k step size to use for FFT [0.05] 

335 k_std: optional k array for standard chi(k). 

336 chi_std: optional chi array for standard chi(k). 

337 nclamp: number of energy end-points for clamp [3] 

338 clamp_lo: weight of low-energy clamp [0] 

339 clamp_hi: weight of high-energy clamp [1] 

340 calc_uncertaintites: Flag to calculate uncertainties in 

341 mu_0(E) and chi(k) [True] 

342 err_sigma: sigma level for uncertainties in mu_0(E) and chi(k) [1] 

343 

344 Output arrays are written to the provided group. 

345 

346 Follows the 'First Argument Group' convention. 

347 """ 

348 msg = sys.stdout.write 

349 if _larch is not None: 

350 msg = _larch.writer.write 

351 if 'kw' in kws: 

352 kweight = kws.pop('kw') 

353 if len(kws) > 0: 

354 msg('Unrecognized arguments for autobk_lmfit():\n') 

355 msg(' %s\n' % (', '.join(kws.keys()))) 

356 return 

357 energy, mu, group = parse_group_args(energy, members=('energy', 'mu'), 

358 defaults=(mu,), group=group, 

359 fcn_name='autobk') 

360 if len(energy.shape) > 1: 

361 energy = energy.squeeze() 

362 if len(mu.shape) > 1: 

363 mu = mu.squeeze() 

364 

365 energy = remove_dups(energy, tiny=TINY_ENERGY) 

366 # if e0 or edge_step are not specified, get them, either from the 

367 # passed-in group or from running pre_edge() 

368 group = set_xafsGroup(group, _larch=_larch) 

369 

370 if edge_step is None and isgroup(group, 'edge_step'): 

371 edge_step = group.edge_step 

372 if e0 is not None and ek0 is None: # command-line e0 still valid 

373 ek0 = e0 

374 if ek0 is None and isgroup(group, 'ek0'): 

375 ek0 = group.ek0 

376 if ek0 is None and isgroup(group, 'e0'): 

377 ek0 = group.e0 

378 

379 if ek0 is not None and (ek0 < energy.min() or ek0 > energy.max()): 

380 ek0 = None 

381 if ek0 is None or edge_step is None: 

382 # need to run pre_edge: 

383 pre_kws = dict(nnorm=None, nvict=0, pre1=None, 

384 pre2=None, norm1=None, norm2=None) 

385 if pre_edge_kws is not None: 

386 pre_kws.update(pre_edge_kws) 

387 pre_edge(energy, mu, group=group, _larch=_larch, **pre_kws) 

388 if ek0 is None: 

389 ek0 = group.e0 

390 if edge_step is None: 

391 edge_step = group.edge_step 

392 if ek0 is None or edge_step is None: 

393 msg('autobk() could not determine ek0 or edge_step!: trying running pre_edge first\n') 

394 return 

395 

396 # get array indices for rkbg and ek0: irbkg, iek0 

397 iek0 = index_of(energy, ek0) 

398 rgrid = np.pi/(kstep*nfft) 

399 if rbkg < 2*rgrid: rbkg = 2*rgrid 

400 

401 # save ungridded k (kraw) and grided k (kout) 

402 # and ftwin (*k-weighting) for FT in residual 

403 enpe = energy[iek0:] - ek0 

404 kraw = np.sign(enpe)*np.sqrt(ETOK*abs(enpe)) 

405 if kmax is None: 

406 kmax = max(kraw) 

407 else: 

408 kmax = max(0, min(max(kraw), kmax)) 

409 kout = kstep * np.arange(int(1.01+kmax/kstep), dtype='float64') 

410 iemax = min(len(energy), 2+index_of(energy, ek0+kmax*kmax/ETOK)) - 1 

411 

412 # interpolate provided chi(k) onto the kout grid 

413 if chi_std is not None and k_std is not None: 

414 chi_std = np.interp(kout, k_std, chi_std) 

415 # pre-load FT window 

416 ftwin = kout**kweight * ftwindow(kout, xmin=kmin, xmax=kmax, 

417 window=win, dx=dk, dx2=dk) 

418 # calc k-value and initial guess for y-values of spline params 

419 nspl = 1 + int(2*rbkg*(kmax-kmin)/np.pi) 

420 irbkg = int(1 + (nspl-1)*np.pi/(2*rgrid*(kmax-kmin))) 

421 if nknots is not None: 

422 nspl = nknots 

423 nspl = max(5, min(128, nspl)) 

424 spl_y, spl_k = np.ones(nspl), np.zeros(nspl) 

425 for i in range(nspl): 

426 q = kmin + i*(kmax-kmin)/(nspl - 1) 

427 ik = index_nearest(kraw, q) 

428 i1 = min(len(kraw)-1, ik + 5) 

429 i2 = max(0, ik - 5) 

430 spl_k[i] = kraw[ik] 

431 spl_y[i] = (2*mu[ik+iek0] + mu[i1+iek0] + mu[i2+iek0] ) / 4.0 

432 

433 order = 3 

434 qmin, qmax = spl_k[0], spl_k[nspl-1] 

435 knots = [spl_k[0] - 1.e-4*(order-i) for i in range(order)] 

436 

437 for i in range(order, nspl): 

438 knots.append((i-order)*(qmax - qmin)/(nspl-order+1)) 

439 qlast = knots[-1] 

440 for i in range(order+1): 

441 knots.append(qlast + 1.e-4*(i+1)) 

442 

443 # coefs = [mu[index_nearest(energy, ek0 + q**2/ETOK)] for q in knots] 

444 knots, coefs, order = splrep(spl_k, spl_y, k=order) 

445 coefs[nspl:] = coefs[nspl-1] 

446 

447 # set fit parameters from initial coefficients 

448 params = Parameters() 

449 for i in range(len(coefs)): 

450 params.add(name = FMT_COEF % i, value=coefs[i], vary=i<len(spl_y)) 

451 

452 initbkg, initchi = spline_eval(kraw[:iemax-iek0+1], mu[iek0:iemax+1], 

453 knots, coefs, order, kout) 

454 

455 # do fit 

456 result = minimize(_lmfit_resid, params, method='leastsq', 

457 gtol=1.e-6, ftol=1.e-6, xtol=1.e-6, epsfcn=1.e-6, 

458 kws = dict(ncoefs=len(coefs), chi_std=chi_std, 

459 knots=knots, order=order, 

460 kraw=kraw[:iemax-iek0+1], 

461 mu=mu[iek0:iemax+1], irbkg=irbkg, kout=kout, 

462 ftwin=ftwin, kweight=kweight, 

463 nfft=nfft, nclamp=nclamp, 

464 clamp_lo=clamp_lo, clamp_hi=clamp_hi)) 

465 

466 # write final results 

467 coefs = [result.params[FMT_COEF % i].value for i in range(len(coefs))] 

468 bkg, chi = spline_eval(kraw[:iemax-iek0+1], mu[iek0:iemax+1], 

469 knots, coefs, order, kout) 

470 obkg = mu[:]*1.0 

471 obkg[iek0:iek0+len(bkg)] = bkg 

472 

473 # outputs to group 

474 group = set_xafsGroup(group, _larch=_larch) 

475 group.bkg = obkg 

476 group.chie = (mu-obkg)/edge_step 

477 group.k = kout 

478 group.chi = chi/edge_step 

479 group.ek0 = ek0 

480 group.rbkg = rbkg 

481 

482 # now fill in 'autobk_details' group 

483 details = Group(kmin=kmin, kmax=kmax, irbkg=irbkg, nknots=len(spl_k), 

484 knots_k=knots, init_knots_y=spl_y, nspl=nspl, 

485 init_chi=initchi/edge_step, report=fit_report(result)) 

486 details.init_bkg = mu[:]*1.0 

487 details.init_bkg[iek0:iek0+len(bkg)] = initbkg 

488 details.knots_y = np.array([coefs[i] for i in range(nspl)]) 

489 group.autobk_details = details 

490 for attr in ('nfev', 'redchi', 'chisqr', 'aic', 'bic', 'params'): 

491 setattr(details, attr, getattr(result, attr, None)) 

492 

493 # uncertainties in mu0 and chi: can be fairly slow. 

494 covar = getattr(result, 'covar', None) 

495 if calc_uncertainties and covar is not None: 

496 nchi = len(chi) 

497 nmue = iemax-iek0 + 1 

498 redchi = result.redchi 

499 covar = result.covar / redchi 

500 jac_chi = np.zeros(nchi*nspl).reshape((nspl, nchi)) 

501 jac_bkg = np.zeros(nmue*nspl).reshape((nspl, nmue)) 

502 

503 cvals, cerrs = [], [] 

504 for i in range(len(coefs)): 

505 par = result.params[FMT_COEF % i] 

506 cvals.append(getattr(par, 'value', 0.0)) 

507 cdel = getattr(par, 'stderr', 0.0) 

508 if cdel is None: 

509 cdel = 0.0 

510 cerrs.append(cdel/2.0) 

511 cvals = np.array(cvals) 

512 cerrs = np.array(cerrs) 

513 

514 # find derivatives by hand! 

515 _k = kraw[:nmue] 

516 _m = mu[iek0:iemax+1] 

517 for i in range(nspl): 

518 cval0 = cvals[i] 

519 cvals[i] = cval0 + cerrs[i] 

520 bkg1, chi1 = spline_eval(_k, _m, knots, cvals, order, kout) 

521 

522 cvals[i] = cval0 - cerrs[i] 

523 bkg2, chi2 = spline_eval(_k, _m, knots, cvals, order, kout) 

524 

525 cvals[i] = cval0 

526 jac_chi[i] = (chi1 - chi2) / (2*cerrs[i]) 

527 jac_bkg[i] = (bkg1 - bkg2) / (2*cerrs[i]) 

528 

529 dfchi = np.zeros(nchi) 

530 dfbkg = np.zeros(nmue) 

531 for i in range(nspl): 

532 for j in range(nspl): 

533 dfchi += jac_chi[i]*jac_chi[j]*covar[i,j] 

534 dfbkg += jac_bkg[i]*jac_bkg[j]*covar[i,j] 

535 

536 prob = 0.5*(1.0 + erf(err_sigma/np.sqrt(2.0))) 

537 dchi = t.ppf(prob, nchi-nspl) * np.sqrt(dfchi*redchi) 

538 dbkg = t.ppf(prob, nmue-nspl) * np.sqrt(dfbkg*redchi) 

539 

540 if not any(np.isnan(dchi)): 

541 group.delta_chi = dchi 

542 group.delta_bkg = 0.0*mu 

543 group.delta_bkg[iek0:iek0+len(dbkg)] = dbkg