1
|
|
|
""" Flare analysis """ |
2
|
|
|
|
3
|
|
|
import logging |
4
|
|
|
|
5
|
|
|
import numpy as np |
|
|
|
|
6
|
|
|
import scipy as sp |
|
|
|
|
7
|
|
|
import pandas as pd |
|
|
|
|
8
|
|
|
|
9
|
|
|
import matplotlib as mplt |
|
|
|
|
10
|
|
|
import matplotlib.pyplot as plt |
|
|
|
|
11
|
|
|
|
12
|
|
|
import sys |
|
|
|
|
13
|
|
|
import os |
|
|
|
|
14
|
|
|
import glob |
|
|
|
|
15
|
|
|
|
16
|
|
|
import re |
|
|
|
|
17
|
|
|
from datetime import datetime |
|
|
|
|
18
|
|
|
|
19
|
|
|
from astropy.stats import bayesian_blocks |
|
|
|
|
20
|
|
|
|
21
|
|
|
import warnings |
|
|
|
|
22
|
|
|
|
23
|
|
|
import mutis.flares.flare as flare |
24
|
|
|
|
25
|
|
|
log = logging.getLogger(__name__) |
26
|
|
|
|
27
|
|
|
class BayesianBlocks: |
28
|
|
|
""" Return a Bayesian Block representation of the signal |
29
|
|
|
|
|
|
|
|
30
|
|
|
Attributes |
31
|
|
|
---------- |
32
|
|
|
edges: list |
33
|
|
|
list of edges defining the blocks. |
34
|
|
|
values: list |
35
|
|
|
list of values defining the height of the blocks. |
36
|
|
|
signal: mutis.Signal |
37
|
|
|
the mutis.Signal() instance that this BayesianBlocks() represents. |
38
|
|
|
inflare: list |
39
|
|
|
list of boolean representing whether a block is in flare or not. |
40
|
|
|
""" |
41
|
|
|
|
42
|
|
|
def __init__(self, signal, p=0.1): |
43
|
|
|
edges = bayesian_blocks(signal.times, signal.values, signal.dvalues, fitness='measures', p0=p) |
|
|
|
|
44
|
|
|
|
|
|
|
|
45
|
|
|
values = list() |
46
|
|
|
for i in range( len(edges)-1 ): |
|
|
|
|
47
|
|
|
msk = (edges[i] < signal.times) & (signal.times < edges[i+1]) |
48
|
|
|
if np.sum(msk) > 0: |
49
|
|
|
value = np.average( signal.values[msk], weights=1/signal.dvalues[msk] ) |
|
|
|
|
50
|
|
|
values.append( value ) |
|
|
|
|
51
|
|
|
else: |
52
|
|
|
values.append( 0 ) |
|
|
|
|
53
|
|
|
|
|
|
|
|
54
|
|
|
self.edges = np.array(edges) |
55
|
|
|
self.values = np.array(values) |
56
|
|
|
self.signal = signal |
57
|
|
|
|
|
|
|
|
58
|
|
|
self.inflare = None |
59
|
|
|
|
|
|
|
|
60
|
|
|
def plot(self, style='color.strict', ax=None, **kwargs): |
|
|
|
|
61
|
|
|
""" Plot the bayesian block representation |
62
|
|
|
|
|
|
|
|
63
|
|
|
E.g: |
64
|
|
|
>>> signal = mutis.Signal(data['jyear'], data['Flux'], data['Flux_err']) |
65
|
|
|
>>> BayesianBlocks(signal).plot(color='k') |
66
|
|
|
|
67
|
|
|
""" |
68
|
|
|
|
|
|
|
|
69
|
|
|
ax = plt.gca() if ax is None else ax |
70
|
|
|
|
|
|
|
|
71
|
|
|
#self.signal.plot() |
72
|
|
|
|
|
|
|
|
73
|
|
|
x = self.edges |
|
|
|
|
74
|
|
|
y = self.values[np.r_[0:len(self.values),-1]] |
|
|
|
|
75
|
|
|
|
76
|
|
|
|
|
|
|
|
77
|
|
|
if self.inflare is None or style == 'none': |
78
|
|
|
|
|
|
|
|
79
|
|
|
ax.step(x, y, 'k', where='post', **kwargs) |
80
|
|
|
|
|
|
|
|
81
|
|
|
else: |
82
|
|
|
|
|
|
|
|
83
|
|
|
if style == 'area': |
84
|
|
|
|
|
|
|
|
85
|
|
|
ax.step(x, y, 'k', where='post') |
86
|
|
|
for i in range(len(self.inflare)): |
87
|
|
|
if self.inflare[i]: |
88
|
|
|
ax.axvspan(x[i], x[i+1], facecolor='r', edgecolor=None, alpha=0.2, **kwargs) |
89
|
|
|
|
|
|
|
|
90
|
|
|
elif style == 'color.simple': |
91
|
|
|
|
|
|
|
|
92
|
|
|
for i in range(len(self.edges)-2): |
93
|
|
|
ax.step(self.edges[[i,i+1]], self.values[[i,i+1]], 'r' if self.inflare[i] else 'k', where='post', **kwargs) |
|
|
|
|
94
|
|
|
ax.step(self.edges[[-2,-1]], self.values[[-1,-1]], 'r' if self.inflare[-1] else 'k', where='post', **kwargs) |
|
|
|
|
95
|
|
|
|
|
|
|
|
96
|
|
|
elif style == 'color.loose': |
97
|
|
|
|
|
|
|
|
98
|
|
|
for i in range(len(self.edges)-2): |
99
|
|
|
ax.step(self.edges[[i,i+1]], self.values[[i,i]], 'r' if self.inflare[i] else 'k', where='post', **kwargs) |
|
|
|
|
100
|
|
|
ax.step(self.edges[[i+1,i+1]], self.values[[i,i+1]], 'r' if self.inflare[i] or self.inflare[i+1] else 'k', where='post', **kwargs) |
|
|
|
|
101
|
|
|
ax.step(self.edges[[-2,-1]], self.values[[-1,-1]], 'r' if self.inflare[-1] else 'k', where='post', **kwargs) |
|
|
|
|
102
|
|
|
|
|
|
|
|
103
|
|
|
elif style == 'color.strict': |
104
|
|
|
|
|
|
|
|
105
|
|
|
for i in range(len(self.edges)-2): |
106
|
|
|
ax.step(self.edges[[i,i+1]], self.values[[i,i]], 'r' if self.inflare[i] else 'k', where='post', **kwargs) |
|
|
|
|
107
|
|
|
ax.step(self.edges[[i+1,i+1]], self.values[[i,i+1]], 'r' if self.inflare[i] and self.inflare[i+1] else 'k', where='post', **kwargs) |
|
|
|
|
108
|
|
|
ax.step(self.edges[[-2,-1]], self.values[[-1,-1]], 'r' if self.inflare[-1] else 'k', where='post', **kwargs) |
|
|
|
|
109
|
|
|
|
|
|
|
|
110
|
|
|
else: |
111
|
|
|
raise Exception("Unknown style. Available styles are 'none', 'area', 'color.simple', 'color.strict', 'color.loose'.") |
|
|
|
|
112
|
|
|
|
113
|
|
|
|
|
|
|
|
114
|
|
|
|
|
|
|
|
115
|
|
|
def get_flares(self, thresh=1): |
116
|
|
|
""" Get a list of flares following the algorithm proposed in |
|
|
|
|
117
|
|
|
[Meyer, Scargle, Blandford (2019)] |
118
|
|
|
(https://iopscience.iop.org/article/10.3847/1538-4357/ab1651/pdf): |
119
|
|
|
|
|
|
|
|
120
|
|
|
```There is no generally accepted consensus on the best way to |
121
|
|
|
determine which data points belong to a flaring state and which |
122
|
|
|
characterize the quiescent level. Nalewajko (2013)suggested |
123
|
|
|
the simple definition that a flare is a continuous time interval |
124
|
|
|
associated with a flux peak in which the flux is larger than half |
125
|
|
|
the peak flux value. This definition is intuitive, however, and it |
126
|
|
|
is unclear how to treat overlapping flares and identify flux |
127
|
|
|
peaks in an objective way. Here we use a simple two-step |
128
|
|
|
procedure tailored to the block representation: (1)identify a |
129
|
|
|
block that is higher than both the previous and subsequent |
130
|
|
|
blocks as a peak, and (2)proceed downward from the peak in |
131
|
|
|
both directions as long as the blocks are successively lower.``` |
132
|
|
|
""" |
133
|
|
|
di_max = 5 |
134
|
|
|
|
|
|
|
|
135
|
|
|
# get list of local maxima |
136
|
|
|
# bad beheaviour at 0 and -1: |
137
|
|
|
#imaxL = sp.signal.argrelextrema(self.values, np.greater)[0] |
138
|
|
|
imaxL = list() |
|
|
|
|
139
|
|
|
for i in range(0,len(self.values)): |
|
|
|
|
140
|
|
|
if i == 0 and thresh*self.values[0] > self.values[1]: |
141
|
|
|
imaxL.append(i) |
142
|
|
|
elif i == len(self.values)-1 and self.values[i-1] < thresh*self.values[i]: |
143
|
|
|
imaxL.append(i) |
144
|
|
|
elif self.values[i-1] < thresh*self.values[i] > self.values[i+1]: |
145
|
|
|
imaxL.append(i) |
146
|
|
|
imaxL = np.array(imaxL) |
|
|
|
|
147
|
|
|
|
|
|
|
|
148
|
|
|
|
|
|
|
|
149
|
|
|
inflareL = np.full(self.values.shape, False) |
|
|
|
|
150
|
|
|
|
|
|
|
|
151
|
|
|
# get list of local maxima over a threshold (flare peaks) |
152
|
|
|
for imax in imaxL: |
153
|
|
|
inflare = True |
154
|
|
|
|
|
|
|
|
155
|
|
|
if imax == 0: |
156
|
|
|
pass |
157
|
|
|
elif not (self.values[imax-1] < thresh*self.values[imax]): |
|
|
|
|
158
|
|
|
inflare = False |
159
|
|
|
|
|
|
|
|
160
|
|
|
|
|
|
|
|
161
|
|
|
if imax == len(self.values)-1: |
162
|
|
|
pass |
163
|
|
|
elif not (thresh*self.values[imax] > self.values[imax+1]): |
|
|
|
|
164
|
|
|
inflare = False |
165
|
|
|
|
|
|
|
|
166
|
|
|
inflareL[imax] = inflare |
167
|
|
|
|
|
|
|
|
168
|
|
|
# extend the flare to adyacent blocks |
169
|
|
|
for imax in np.argwhere(inflareL): |
170
|
|
|
|
|
|
|
|
171
|
|
|
if imax == 0: |
172
|
|
|
pass |
173
|
|
|
else: |
174
|
|
|
di = 1 |
|
|
|
|
175
|
|
|
stop = False |
176
|
|
|
while not stop and di < di_max: |
177
|
|
|
if imax-di-1 < 0: |
178
|
|
|
break |
179
|
|
|
|
|
|
|
|
180
|
|
|
if self.values[imax-di-1] < thresh*self.values[imax-di]: |
181
|
|
|
inflareL[imax-di] = True |
182
|
|
|
else: |
183
|
|
|
stop = True |
184
|
|
|
|
|
|
|
|
185
|
|
|
|
|
|
|
|
186
|
|
|
di = di + 1 |
|
|
|
|
187
|
|
|
|
|
|
|
|
188
|
|
|
if imax == len(self.values)-1: |
189
|
|
|
pass |
190
|
|
|
else: |
191
|
|
|
di = 1 |
|
|
|
|
192
|
|
|
stop = False |
193
|
|
|
while not stop and di < di_max: |
194
|
|
|
if imax+di+1 > len(self.values)-1: |
195
|
|
|
break |
196
|
|
|
|
|
|
|
|
197
|
|
|
if thresh*self.values[imax+di] > self.values[imax+di+1]: |
198
|
|
|
inflareL[imax+di] = True |
199
|
|
|
else: |
200
|
|
|
stop = True |
201
|
|
|
|
|
|
|
|
202
|
|
|
|
|
|
|
|
203
|
|
|
di = di + 1 |
|
|
|
|
204
|
|
|
|
|
|
|
|
205
|
|
|
|
|
|
|
|
206
|
|
|
self.inflare = np.array(inflareL) |
207
|
|
|
|
|
|
|
|
208
|
|
|
|
209
|
|
|
def get_flare_list(self): |
210
|
|
|
""" Join all flares into a list of mutis.flares.Flare() """ |
211
|
|
|
|
212
|
|
|
groups = np.split(np.arange(len(self.inflare)), |
|
|
|
|
213
|
|
|
np.where((np.abs(np.diff(np.asarray(self.inflare, dtype=int))) == 1))[0]+1) |
|
|
|
|
214
|
|
|
|
|
|
|
|
215
|
|
|
# (groups contains the groups of indices of self.inflare/self.values |
|
|
|
|
216
|
|
|
# that correspond to flares or no flares) |
217
|
|
|
|
|
|
|
|
218
|
|
|
#print(self.inflare) |
219
|
|
|
#print([list(grp) for grp in groups]) |
220
|
|
|
|
|
|
|
|
221
|
|
|
flareL = list() |
|
|
|
|
222
|
|
|
for i, grp in enumerate(groups): |
|
|
|
|
223
|
|
|
if np.all(self.inflare[grp]): |
224
|
|
|
#print(grp) |
225
|
|
|
tstart = self.edges[grp[0]] |
226
|
|
|
tstop = self.edges[grp[-1]+1] # flare ends when no flare begins |
227
|
|
|
|
|
|
|
|
228
|
|
|
flareL.append(flare.Flare(tstart,tstop)) |
|
|
|
|
229
|
|
|
|
|
|
|
|
230
|
|
|
#print(flareL) |
231
|
|
|
|
|
|
|
|
232
|
|
|
return flareL |
|
|
|
|