Attachment 'characterize_variability.py'
Download 1 import pyfits as pf
2 import numpy as np
3 import scipy.optimize as opt
4 from scipy.stats.distributions import chi2
5 ###units 1000 if ksec, 1 if sec!
6 s2u = 1.
7 def heaviside(x):
8 return 0.5 * (np.sign(x) + 1)
9
10 def compute_variability_parameters(rate, rate_err, time, tstart):
11 # read the values in the binned light curve
12 # (time is measured from tstart)
13
14 BestfitPar = BestfitParameters(time, rate, rate_err)
15
16 # extract the chi2 test statistics with its associated P-value
17 # for a rate model linear in time
18 BestfitPar.FitLinear()
19 # for a rate model quadratic in time
20 BestfitPar.FitFlare()
21 # for a constant rate model with an eclipse
22 BestfitPar.FitEclipse()
23 return BestfitPar
24
25 class BestfitParameters:
26 def __init__(self, X=[], Y=[], Y_err=[]):
27 self.N_points = len(X)
28
29 self.__X = X
30 self.__Y = Y
31 self.__Y_err = Y_err
32
33 self.Linear = []
34 self.Flare = []
35 self.Eclipse = []
36
37 self.Hdu = 0
38 self.Hdr = []
39
40
41 def FitLinear(self):
42 self.Linear = compute_chi2_linear(self.__Y, self.__Y_err, self.__X)
43
44 def FitFlare(self):
45 self.Flare = compute_chi2_flare(self.__Y, self.__Y_err, \
46 self.__X)
47
48 def FitEclipse(self):
49 self.Eclipse = compute_chi2_eclipse(self.__Y, self.__Y_err, \
50 self.__X)
51
52
53 def compute_chi2_linear(rate, rate_err, time):
54 n_points = np.size(rate)
55 epsilon = 1.e-8
56 chi2_TS = -1.; chi2_prob = -1.;
57 cons = 0.; cons_err = -1.
58 lin = 0.; lin_err = -1.
59 chi2_dof = n_points - 2
60 if chi2_dof <= 0:
61 print 'linear fit notice: not enough bins, skip'
62 return [chi2_TS, chi2_prob, chi2_dof, cons, cons_err, lin, lin_err]
63 # extract the chi2 test statistics with its associated P-value
64 drate = [np.sqrt(epsilon + d ** 2.) for d in rate_err]
65 params = opt.leastsq(residual_linear, x0=[0.,0.], \
66 args=(rate, drate, time), full_output=True)
67 # get parameters and errors only if leastsq converged
68 if params[1] is not None and params[4] > 0 and params[4] < 5:
69 if params[1][0][0] >= 0. and params[1][1][1] >= 0.:
70 cons = params[0][0]
71 cons_err = np.sqrt(params[1][0][0])
72 lin = params[0][1] * s2u
73 lin_err = np.sqrt(params[1][1][1]) * s2u
74 offsets = residual_linear(params[0], rate, drate, time)
75 chi2_TS = np.inner(offsets, offsets)
76 chi2_prob = chi2.sf(chi2_TS, chi2_dof)
77 else:
78 print 'linear fit error: negative values for errors squared'
79 else:
80 print 'linear fit warning: leastsq did not converge'
81 if chi2_TS >= 0.:
82 print 'linear fit: Slope=%.3g+/-%.2g Const=%.3g+/-%.2g chi2(%d)=%.1f' % \
83 (lin, lin_err, cons, cons_err, chi2_dof, chi2_TS)
84 return [chi2_TS, chi2_prob, chi2_dof, cons, cons_err, lin, lin_err]
85
86
87
88 def compute_chi2_flare(rate, rate_err, time):
89 n_points = np.size(rate)
90 epsilon = 1.e-8
91 chi2_TS = -1.; chi2_prob = -1.;
92 norm = 0.; norm_err = -1.
93 t_burst = 0.; t_burst_err = -1.
94 t_decay = 0.; t_dec_err = -1.
95 cons = 0.; cons_err = -1.
96 chi2_dof = n_points - 4
97 if chi2_dof <= 0:
98 print 'flare fit notice: not enough bins, skip'
99 return [chi2_TS, chi2_prob, chi2_dof, \
100 cons, cons_err, norm, norm_err, t_burst, t_burst_err, t_decay, t_dec_err]
101 # extract the chi2 test statistics with its associated P-value
102 drate = [np.sqrt(epsilon + d ** 2.) for d in rate_err]
103 constr = get_flare_constraints(time)
104 # get a first guess of the global minimum by brute force
105 params = locate_flare(rate, drate, time, constr)
106 # follow up only if there is a flare
107 if params[1] is not None and params[0][0] > 0.:
108 params = opt.leastsq(residual_flare, x0=params[0], \
109 args=(rate, drate, time, constr), full_output=True)
110 # get parameters and errors only if leastsq converged
111 if params[1] is not None and params[4] > 0 and params[4] < 5:
112 if params[1][0][0] >= 0. and params[1][1][1] >= 0. \
113 and params[1][2][2] >= 0. and params[1][3][3] >= 0.:
114 norm = params[0][0]
115 norm_err = np.sqrt(params[1][0][0])
116 t_burst = params[0][1] / s2u
117 t_burst_err = np.sqrt(params[1][1][1]) / s2u
118 t_decay = params[0][2] / s2u
119 t_dec_err = np.sqrt(params[1][2][2]) / s2u
120 cons = params[0][3]
121 cons_err = np.sqrt(params[1][3][3])
122 chi2_TS = res_square_flare(params[0], rate, drate, time, constr)
123 chi2_prob = chi2.sf(chi2_TS, chi2_dof)
124 else:
125 print 'flare fit error: negative values for errors squared'
126 else:
127 print 'flare fit warning: leastsq did not converge'
128 else:
129 print 'flare fit warning: brute did not converge'
130 if chi2_TS >= 0.:
131 print 'flare fit: Norm=%.3g+/-%.2g Time_Burst=%.3g+/-%.2g Time_Decay=%.3g+/-%.2g Const=%.3g+/-%.2g chi2(%d)=%.1f' % \
132 (norm, norm_err, t_burst, t_burst_err, t_decay, t_dec_err, cons, cons_err, chi2_dof, chi2_TS)
133 return [chi2_TS, chi2_prob, chi2_dof, \
134 cons, cons_err, norm, norm_err, t_burst, t_burst_err, t_decay, t_dec_err]
135
136 def compute_chi2_eclipse(rate, rate_err, time):
137 n_points = np.size(rate)
138 epsilon = 1.e-8
139 chi2_TS = -1.; chi2_prob = -1.;
140 t_init = 0.; t_init_err = -1.
141 t_exit = 0.; t_exit_err = -1.
142 drop = 0.; drop_err = -1.
143 cons = 0.; cons_err = -1.
144 chi2_dof = n_points - 4
145 if chi2_dof <= 0:
146 print 'eclipse fit notice: not enough bins, skip'
147 return [chi2_TS, chi2_prob, chi2_dof, \
148 cons, cons_err, drop, drop_err, t_init, t_init_err, t_exit, t_exit_err]
149 # extract the chi2 test statistics with its associated P-value
150 drate = [np.sqrt(epsilon + d ** 2.) for d in rate_err]
151 constr = get_eclipse_constraints(time)
152 # get a first guess of the global minimum by brute force
153 params = locate_eclipse(rate, drate, time, constr)
154 # follow up only if there is a flare
155 if params[1] is not None and params[0][0] > 0.:
156 params = opt.leastsq(residual_eclipse, x0=params[0], \
157 args=(rate, drate, time, constr), full_output=True)
158 # get parameters and errors only if leastsq converged
159 if params[1] is not None and params[4] > 0 and params[4] < 5:
160 if params[1][0][0] >= 0. and params[1][1][1] >= 0. \
161 and params[1][2][2] >= 0. and params[1][3][3] >= 0.:
162 drop = params[0][0]
163 drop_err = np.sqrt(params[1][0][0])
164 t_init = params[0][1] / s2u
165 t_init_err = np.sqrt(params[1][1][1]) / s2u
166 t_exit = params[0][2] / s2u
167 t_exit_err = np.sqrt(params[1][2][2]) / s2u
168 cons = params[0][3]
169 cons_err = np.sqrt(params[1][3][3])
170 chi2_TS = res_square_eclipse(params[0], rate, drate, time, constr)
171 chi2_prob = chi2.sf(chi2_TS, chi2_dof)
172 else:
173 print 'eclipse fit error: negative values for errors squared'
174 else:
175 print 'eclipse fit warning: leastsq did not converge'
176 else:
177 print 'eclipse fit warning: brute did not converge'
178 if chi2_TS >= 0.:
179 print 'eclipse fit: Drop=%.3g+/-%.2g Time_Init=%.3g+/-%.2g Time_Exit=%.3g+/-%.2g Const=%.3g+/-%.2g chi2(%d)=%.1f' % \
180 (drop, drop_err, t_init, t_init_err, t_exit, t_exit_err, cons, cons_err, chi2_dof, chi2_TS)
181 return [chi2_TS, chi2_prob, chi2_dof, \
182 cons, cons_err, drop, drop_err, t_init, t_init_err, t_exit, t_exit_err]
183 def locate_eclipse(rate, drate, time, constr):
184 n_points = np.size(rate)
185 n_steps = 12
186 ti_min, te_max, dtime = constr
187 tirange = slice(ti_min, te_max-dtime, 2.*(te_max-ti_min)/n_points)
188 terange = slice(ti_min+dtime, te_max, 2.*(te_max-ti_min)/n_points)
189 r_min = np.min(rate)
190 r_max = np.max(rate)
191 r_ave = np.mean(rate)
192 crange = slice(r_ave, r_max, (r_max-r_ave)/float(n_steps))
193 drange = slice(0., r_max-r_min, (r_max-r_min)/float(n_steps))
194 return opt.brute(res_square_eclipse, \
195 ranges=(drange, tirange, terange, crange), \
196 args=(rate, drate, time, constr), full_output=True)
197
198 def get_eclipse_constraints(time):
199 dtime = (time[1] - time[0]) / 2.
200 if time[-1] - time[0] > dtime:
201 ti_min = time[0] + dtime
202 te_max = time[-1] - dtime
203 else:
204 tave = (time[-1] + time[0]) / 2.
205 ti_min = tave - dtime / 2.
206 te_max = tave + dtime / 2.
207 return ti_min, te_max, dtime
208
209 def bind_eclipse_params(params, constr):
210 ti_min, te_max, dt = constr
211 d, ti, te, c = params
212 if ti < ti_min:
213 ti = ti_min
214 elif ti > te_max - dt:
215 ti = te_max - dt
216 if d < 0.:
217 d = 0.
218 if te < ti + dt:
219 te = ti + dt
220 elif te > te_max:
221 te = te_max
222 return d, ti, te, c
223
224 def residual_eclipse(params, y, dy, x, constr):
225 d, ti, te, c = bind_eclipse_params(params, constr)
226 dt = 2. * constr[2]
227 # center and duration of the eclipse
228 tc = (ti + te) / 2.
229 dur = te - ti
230 # this way, also partially eclipsed bins are considered
231 y_mod = c - d * np.clip(((dt+dur)/2. - np.fabs(x-tc)) / dt, 0., 1.)
232
233 return (y - y_mod) / dy
234
235 def res_square_eclipse(params, y, dy, x, constr):
236 residuals = residual_eclipse(params, y, dy, x, constr)
237 return np.inner(residuals, residuals)
238
239 def locate_flare(rate, drate, time, constr):
240 n_points = np.size(rate)
241 n_steps = 12
242 tb_min, tb_max, dtime = constr
243 tbrange = slice(tb_min, tb_max, 2.*(tb_max-tb_min)/n_points)
244 tdrange = slice(tb_min, (tb_min+tb_max)/2., (tb_max-tb_min)/float(n_steps))
245 r_min = np.min(rate)
246 r_max = np.max(rate)
247 r_ave = np.mean(rate)
248 crange = slice(r_min, r_ave, (r_ave-r_min)/float(n_steps))
249 nrange = slice(0., r_max-r_min, (r_max-r_min)/float(n_steps))
250 return opt.brute(res_square_flare, \
251 ranges=(nrange, tbrange, tdrange, crange), \
252 args=(rate, drate, time, constr), full_output=True)
253
254 def get_flare_constraints(time):
255 dtime = (time[1] - time[0]) / 2.
256 tb_min = (time[0] + time[1]) / 2.
257 tb_max = (time[-3] + time[-2]) / 2.
258 return tb_min, tb_max, dtime
259
260 def bind_flare_params(params, constr):
261 tb_min, tb_max, dt = constr
262 n, tb, td, c = params
263 if tb < tb_min:
264 tb = tb_min
265 elif tb > tb_max:
266 tb = tb_max
267 if n < 0.:
268 n = 0.
269 if td < dt:
270 td = dt
271 elif td > 1.e6:
272 td = 1.e6
273 return n, tb, td, c
274
275 def residual_flare(params, y, dy, x, constr):
276 n, tb, td, c = bind_flare_params(params, constr)
277 dt = 2. * constr[2]
278 y_mod = n * np.clip((x - tb) / dt + .5, 0., 1.) * \
279 np.exp(-1. * np.clip((x - tb) / td, -2., 10.)) + c
280 return (y - y_mod) / dy
281
282 def res_square_flare(params, y, dy, x, constr):
283 residuals = residual_flare(params, y, dy, x, constr)
284 return np.inner(residuals, residuals)
285
286
287 def residual_linear(params, y, dy, x):
288 c, b = params
289 return (y - (b*x + c)) / dy
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.