Source code for exponential fit of Green function#
\(\phantom{\xi}\)
Listing 2 Content of file
fit_green_tau.py containing the source code used to fit the time-displaced Green function with an exponential decay for determining the single particle gap. See also Section 2.7.2.1 for a derivation of the exponential decay and Appendix A.2.2.10 demonstrates how to apply the fitting function.# 1"""Fit time-displaced Green function exponentially to determine
2single particle gap."""
3# pylint: disable=invalid-name
4
5import numpy as np
6from scipy.optimize import curve_fit
7
8
9def fit_green_tau(taus, G, dG, plot=False):
10 """Fit time-displaced Green function exponentially to determine
11 single particle gap.
12
13 Parameters
14 ----------
15 taus : array-like object
16 Tau values.
17 G : array-like object
18 Green function values
19 dG : array-like object
20 Standard errors of Green functions
21 plot : bool, default=False
22 Plot G(tau).
23
24 Returns
25 -------
26 gap, dgap : floats
27 Single particle gap and standard error.
28 """
29 N_tau = len(G)
30
31 def func1(x, a, b):
32 return a*np.exp(-b*x)
33
34 def func2(x, a, b):
35 return a*np.exp(b*(x-taus[-1]))
36
37 if plot:
38 import matplotlib.pyplot as plt
39 fig, (ax, ax2) = plt.subplots(1, 2, figsize=(14, 5))
40 ax.errorbar(x=taus, y=G, yerr=dG)
41 ax.set_yscale('log')
42
43 W = N_tau//2
44 for i in range(N_tau//2):
45 if G[i]-dG[i] < 1E-7:
46 W = i
47 # print(i)
48 break
49 i1 = W//2 - 5
50 i2 = i1+10
51 popt, pcov, i1_out, i2_out, chi_squared, reduced_chi_squared, accepted = \
52 try_fit(func1, taus, G, dG, i1, i2)
53 gap1 = popt[1]
54 dgap1 = np.sqrt(pcov[1,1])
55 if plot:
56 # print('left')
57 ax.errorbar(
58 x=taus[i1_out:i2_out],
59 y=G[i1_out:i2_out],
60 yerr=dG[i1_out:i2_out]
61 )
62 title = f'{accepted} {i1_out} {i2_out} {i2_out-i1_out} {reduced_chi_squared}'
63 ax2.errorbar(x=[1], y=[gap1], yerr=[dgap1])
64
65 W = N_tau//2
66 for i in range(N_tau//2):
67 if G[-i]-dG[-i] < 1E-7:
68 W = i
69 # print(i)
70 break
71 i2 = N_tau-W//2+5
72 i1 = i2-10
73 out = try_fit(func2, taus, G, dG, i1, i2)
74 # print(out)
75 if out is not None:
76 popt, pcov, i1_out, i2_out, chi_squared, reduced_chi_squared, accepted = out
77 gap2 = popt[1]
78 dgap2 = np.sqrt(pcov[1,1])
79 gap = np.mean([gap1, gap2])
80 dgap = (np.max([gap1+dgap1, gap2+dgap2]) - np.min([gap1-dgap1, gap2-dgap2]))/2
81 if plot:
82 # print('right')
83 ax.errorbar(x=taus[i1_out:i2_out],
84 y=G[i1_out:i2_out],
85 yerr=dG[i1_out:i2_out])
86 ax.set_title(f'{title}\n{accepted} {i1_out} {i2_out} ' +
87 f'{i2_out-i1_out} {reduced_chi_squared}')
88 ax2.errorbar(x=[2], y=[gap2], yerr=[dgap2])
89 ax2.errorbar(x=[3], y=[gap], yerr=[dgap])
90 plt.show()
91 plt.close()
92
93 return gap, dgap
94
95
96def try_fit(func, x, y, dy, i1, i2, popt=None):
97 """Fit function func to data x[i1:i2], y[i1:i2], dy[i1:i2] and alternately
98 decrease i1 and increase i2 until some standards for the quality of fit are
99 no more fulfilled."""
100 this_side = 1
101 side1_fin = False
102 side2_fin = False
103 accepted = False
104
105 while True:
106 if i1 < 0 or i2 >= len(y):
107 return popt, pcov, i1, i2, chi_squared, reduced_chi_squared, accepted
108 try:
109 popt, pcov = curve_fit(func, x[i1:i2], y[i1:i2], sigma=dy[i1:i2],
110 absolute_sigma=True, p0=popt)
111 except RuntimeError as curve_fit_failed:
112 try:
113 return popt, pcov, i1, i2, chi_squared, reduced_chi_squared, accepted
114 except NameError:
115 raise Exception("curve_fit did not converge") from curve_fit_failed
116
117 # ================== Estimate quality of fit ==================
118 chi_squared = np.sum(((func(x[i1:i2], *popt)-y[i1:i2])/dy[i1:i2])**2)
119 reduced_chi_squared = (chi_squared)/(len(x[i1:i2])-len(popt))
120
121 n_ges = 0
122 n_in_h = 0
123 n_in = 0
124 n_in_4 = 0
125 n_end1 = 0
126 n_end2 = 0
127
128 for i in range(i1, i2):
129 n_ges = n_ges + 1
130 if abs( func(x[i], *popt) - y[i] ) < 0.7 * dy[i]:
131 n_in_h = n_in_h + 1
132 if abs( func(x[i], *popt) - y[i] ) < 4 * dy[i]:
133 n_in_4 = n_in_4 + 1
134 if abs( func(x[i], *popt) - y[i] ) < dy[i]:
135 n_in = n_in + 1
136
137 for i in range(i1,i1+3):
138 if abs( func(x[i], *popt) - y[i] ) < 2*dy[i]:
139 n_end1 = n_end1 + 1
140
141 for i in range(i2-3,i2):
142 if abs( func(x[i], *popt) - y[i] ) < 2*dy[i]:
143 n_end2 = n_end2 + 1
144
145 ratio_h = n_in_h / n_ges
146 ratio = n_in / n_ges
147 ratio_4 = n_in_4 / n_ges
148
149 # line = "Fitting from {} to {}: {:5.2f} {:5.2f} "
150 # "{} {:5.2f} {:5.2f} {:5.2f} {} {}"
151 # line = line.format(i1, i2, chi_squared, reduced_chi_squared,
152 # n_ges, ratio_h, ratio, ratio_4, n_end1, n_end2)
153 # print(line)
154
155 # ====== Set minimal standard for quality of fit ======
156 success = ((ratio > 0.68) and (ratio_4 > 0.99)
157 and (n_end1 > 1) and (n_end2 > 1))
158
159 if not success:
160 # print('rejected!')
161 if this_side == 1:
162 side1_fin = True
163 i1 = i1 + 1
164 elif this_side == 2:
165 side2_fin = True
166 i2 = i2 - 1
167 else:
168 # print('accepted')
169 accepted = True
170
171 if i1 == 0:
172 side1_fin = True
173 if i2 == len(x)-1:
174 side2_fin = True
175
176 if this_side == 1 and side2_fin is False:
177 i2 = i2+1
178 this_side = 2
179 elif this_side == 2 and side1_fin is False:
180 i1 = i1-1
181 this_side = 1
182 elif side2_fin is False:
183 i2 = i2+2
184 this_side = 2
185 elif side1_fin is False:
186 i1 = i1-1
187 this_side = 1
188 else:
189 return popt, pcov, i1, i2, chi_squared, reduced_chi_squared, accepted