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
« 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
11import uncertainties
13from larch import (Group, Make_CallArgs, parse_group_args, isgroup)
14from larch.math import index_of, index_nearest, realimag, remove_dups
16from .xafsutils import ETOK, TINY_ENERGY, set_xafsGroup
17from .xafsft import ftwindow, xftf_fast
18from .pre_edge import find_e0, pre_edge
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
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:]))
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
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]
83 Output arrays are written to the provided group.
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)
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
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
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)
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
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)
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
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)]
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))
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
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)
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)
200 best, covar, _infodict, errmsg, ier = lsout
201 final_coefs = coefs[:]
202 final_coefs[:nspl] = best[:]
203 final_coefs[nspl:] = best[-1]
205 chisqr = ((_resid(best, *userargs))**2).sum()
206 redchi = chisqr / (2*irbkg+2*nclamp - nspl)
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
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
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
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)
237 if calc_uncertainties and covar is not None:
238 autobk_delta_chi(group, err_sigma=err_sigma)
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
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]
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])
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]
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)
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
287## version of autobk using lmfit
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:]))
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
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]
344 Output arrays are written to the provided group.
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()
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)
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
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
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
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
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
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)]
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))
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]
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))
452 initbkg, initchi = spline_eval(kraw[:iemax-iek0+1], mu[iek0:iemax+1],
453 knots, coefs, order, kout)
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))
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
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
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))
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))
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)
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)
522 cvals[i] = cval0 - cerrs[i]
523 bkg2, chi2 = spline_eval(_k, _m, knots, cvals, order, kout)
525 cvals[i] = cval0
526 jac_chi[i] = (chi1 - chi2) / (2*cerrs[i])
527 jac_bkg[i] = (bkg1 - bkg2) / (2*cerrs[i])
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]
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)
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