|
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 |
|
|
|
|
|
IAA-CSIC /
MUTIS
Loading history...
