|
1
|
|
|
# Licensed under a 3-clause BSD style license - see LICENSE |
|
2
|
|
|
"""Synthetic generation of light curves.""" |
|
3
|
|
|
|
|
4
|
|
|
import logging |
|
5
|
|
|
|
|
6
|
|
|
import matplotlib.pyplot as plt |
|
|
|
|
|
|
7
|
|
|
import numpy as np |
|
|
|
|
|
|
8
|
|
|
import scipy.optimize as scipy_optimize |
|
|
|
|
|
|
9
|
|
|
import scipy.special as scipy_special |
|
|
|
|
|
|
10
|
|
|
import scipy.stats as scipy_stats |
|
|
|
|
|
|
11
|
|
|
from matplotlib.offsetbox import AnchoredText |
|
|
|
|
|
|
12
|
|
|
|
|
13
|
|
|
from mutis.lib.signal import * |
|
|
|
|
|
|
14
|
|
|
|
|
15
|
|
|
__all__ = ["Signal"] |
|
16
|
|
|
|
|
17
|
|
|
log = logging.getLogger(__name__) |
|
18
|
|
|
|
|
19
|
|
|
|
|
20
|
|
|
|
|
21
|
|
|
class Signal: |
|
22
|
|
|
"""Analysis and generation of a signal. |
|
23
|
|
|
|
|
24
|
|
|
Description goes here. |
|
25
|
|
|
|
|
26
|
|
|
Parameters |
|
27
|
|
|
---------- |
|
28
|
|
|
times : :class:`numpy.ndarray` or :class:`pandas.Series` |
|
29
|
|
|
Values of the time axis. |
|
30
|
|
|
values : :class:`numpy.ndarray` or :class:`pandas.Series` |
|
31
|
|
|
Values of the signal axis. |
|
32
|
|
|
fgen : :class:`str` |
|
33
|
|
|
Method used to generate the synthetic signal. |
|
34
|
|
|
""" |
|
35
|
|
|
|
|
36
|
|
|
def __init__(self, times, values, fgen): |
|
37
|
|
|
self.times = np.array(times) |
|
38
|
|
|
self.values = np.array(values) |
|
39
|
|
|
self.fgen = fgen |
|
40
|
|
|
self.synth = None |
|
41
|
|
|
|
|
42
|
|
|
# TODO make attributes below specific of OU method / not the entire class |
|
|
|
|
|
|
43
|
|
|
self.theta = None |
|
44
|
|
|
self.mu = None |
|
|
|
|
|
|
45
|
|
|
self.sigma = None |
|
46
|
|
|
|
|
47
|
|
|
def gen_synth(self, samples): |
|
48
|
|
|
"""Description goes here.""" |
|
49
|
|
|
|
|
50
|
|
|
self.synth = np.empty((samples, self.times.size)) |
|
51
|
|
|
for n in range(samples): |
|
|
|
|
|
|
52
|
|
|
if self.fgen == "lc_gen_samp": |
|
53
|
|
|
self.synth[n] = lc_gen_samp(self.values) |
|
54
|
|
|
elif self.fgen == "lc_gen_psd_nft": |
|
55
|
|
|
self.synth[n] = lc_gen_psd_nft(self.times, self.values) |
|
56
|
|
|
elif self.fgen == "lc_gen_psd_lombscargle": |
|
57
|
|
|
self.synth[n] = lc_gen_psd_lombscargle(self.times, self.values) |
|
58
|
|
|
elif self.fgen == "lc_gen_psd_fft": |
|
59
|
|
|
self.synth[n] = lc_gen_psd_fft(self.values) |
|
60
|
|
|
elif self.fgen == "lc_gen_ou": |
|
61
|
|
|
if self.theta is None or self.mu is None or self.sigma is None: |
|
62
|
|
|
raise Exception("You need to set the parameters for the signal") |
|
63
|
|
|
self.synth[n] = lc_gen_ou(self.theta, self.mu, self.sigma, self.times) |
|
64
|
|
|
else: |
|
65
|
|
|
raise Exception(f"Unknown fgen method {self.fgen}") |
|
66
|
|
|
|
|
67
|
|
|
def OU_fit(self, bins=None, rang=None, a=1e-5, b=100): |
|
|
|
|
|
|
68
|
|
|
"""Description goes here.""" |
|
69
|
|
|
|
|
70
|
|
|
# TODO: make a generic fit method |
|
|
|
|
|
|
71
|
|
|
res = dict() |
|
72
|
|
|
|
|
73
|
|
|
# estimate sigma |
|
74
|
|
|
try: |
|
75
|
|
|
# dy = np.diff(self.values) |
|
76
|
|
|
# dt = np.diff(self.times) |
|
77
|
|
|
sigma_est = (np.nanmean(np.diff(self.values) ** 2 / self.values[:-1] ** 2 / np.diff(self.times))) ** 0.5 |
|
|
|
|
|
|
78
|
|
|
res["sigma_est"] = sigma_est |
|
79
|
|
|
except Exception as e: |
|
|
|
|
|
|
80
|
|
|
raise Exception(f"Could not estimate sigma: {e}") |
|
81
|
|
|
|
|
82
|
|
|
# plot histogram |
|
83
|
|
|
if bins is None: |
|
84
|
|
|
bins = np.int(self.values.size ** 0.5 / 1.5) # bins='auto' |
|
85
|
|
|
if rang is None: |
|
86
|
|
|
rang = (np.percentile(self.values, 0), np.percentile(self.values, 99)) |
|
87
|
|
|
p, x = np.histogram(self.values, density=True, bins=bins, range=rang) # bins='sqrt') |
|
|
|
|
|
|
88
|
|
|
x = (x + np.roll(x, -1))[:-1] / 2.0 |
|
|
|
|
|
|
89
|
|
|
|
|
90
|
|
|
plt.subplots() |
|
91
|
|
|
plt.hist(self.values, density=True, alpha=0.75, bins=bins, range=rang) |
|
92
|
|
|
plt.plot(x, p, "r-", alpha=0.5) |
|
93
|
|
|
anchored_text = AnchoredText( |
|
94
|
|
|
f"mean {np.mean(self.values):.2g} \n " |
|
95
|
|
|
"median {np.median(self.values):.2g} \n " |
|
96
|
|
|
"mode {scipy_stats.mode(self.values)[0][0]:.2g} \n " |
|
97
|
|
|
"std {np.std(self.values):.2g} \n " |
|
98
|
|
|
"var {np.var(self.values):.2g}", |
|
99
|
|
|
loc="upper right", |
|
100
|
|
|
) |
|
101
|
|
|
plt.gca().add_artist(anchored_text) |
|
102
|
|
|
|
|
103
|
|
|
# curve_fit |
|
104
|
|
|
try: |
|
105
|
|
|
popt, pcov = scipy_optimize.curve_fit(f=self.pdf, xdata=x, ydata=p) |
|
106
|
|
|
# print('curve_fit: (l, mu)') |
|
107
|
|
|
# print('popt: ') |
|
108
|
|
|
# print(popt) |
|
109
|
|
|
# print('pcov: ') |
|
110
|
|
|
# print(np.sqrt(np.diag(pcov))) |
|
111
|
|
|
x_c = np.linspace(1e-5, 1.1 * np.max(x), 1000) |
|
112
|
|
|
plt.plot(x_c, self.pdf(x_c, *popt), "k-", label="curve_fit", alpha=0.8) |
|
113
|
|
|
res["curve_fit"] = (popt, np.sqrt(np.diag(pcov))) |
|
114
|
|
|
except Exception as e: |
|
|
|
|
|
|
115
|
|
|
raise Exception(f"Some error fitting with curve_fit {e}") |
|
116
|
|
|
|
|
117
|
|
|
# TODO: place outside as a helper class |
|
|
|
|
|
|
118
|
|
|
# and understand how it is used |
|
119
|
|
|
# is it mu really needed ? |
|
120
|
|
|
# mu here is different than self.mu |
|
121
|
|
|
# fit pdf with MLE |
|
122
|
|
|
class OU(scipy_stats.rv_continuous): |
|
|
|
|
|
|
123
|
|
|
def _pdf(self, x, l, mu): |
|
|
|
|
|
|
124
|
|
|
return (l * mu) ** (1 + l) / scipy_special.gamma(1 + l) * np.exp(-l * mu / x) / x ** (l + 2) |
|
|
|
|
|
|
125
|
|
|
|
|
126
|
|
|
try: |
|
127
|
|
|
fit = OU(a=a, b=np.percentile(self.values, b)).fit(self.values, 1, 1, floc=0, fscale=1) |
|
128
|
|
|
# print('MLE fit: (l, mu)') |
|
129
|
|
|
# print(fit) |
|
130
|
|
|
x_c = np.linspace(0, 1.1 * np.max(x), 1000) |
|
131
|
|
|
plt.plot(x_c, self.pdf(x_c, fit[0], fit[1]), "k-.", label="MLE", alpha=0.8) |
|
132
|
|
|
res["MLE_fit"] = fit[:-2] |
|
133
|
|
|
except Exception as e: |
|
|
|
|
|
|
134
|
|
|
raise Exception(f"Some error fitting with MLE {e}") |
|
135
|
|
|
plt.legend(loc="lower right") |
|
136
|
|
|
plt.show() |
|
137
|
|
|
|
|
138
|
|
|
# estimate theta |
|
139
|
|
|
res["th_est1"] = fit[0] * sigma_est ** 2 / 2 |
|
140
|
|
|
res["th_est2"] = popt[0] * sigma_est ** 2 / 2 |
|
141
|
|
|
|
|
142
|
|
|
return res |
|
143
|
|
|
|
|
144
|
|
|
def OU_check_gen(self, theta, mu, sigma): |
|
|
|
|
|
|
145
|
|
|
"""Description goes here.""" |
|
146
|
|
|
|
|
147
|
|
|
# TODO: make a generic check_gen method |
|
|
|
|
|
|
148
|
|
|
|
|
149
|
|
|
t, y = self.times, self.values |
|
|
|
|
|
|
150
|
|
|
y2 = lc_gen_ou(theta, mu, sigma, self.times, scale=np.std(self.values), loc=np.mean(self.values)) |
|
|
|
|
|
|
151
|
|
|
|
|
152
|
|
|
# plot the two signals |
|
153
|
|
|
fig, ax = plt.subplots() |
|
|
|
|
|
|
154
|
|
|
ax.plot(t, y, "b-", label="orig", lw=0.5, alpha=0.8) |
|
155
|
|
|
ax2 = ax.twinx() |
|
156
|
|
|
ax2.plot(t, y2, "r-", label="gen", lw=0.5, alpha=0.8) |
|
157
|
|
|
plt.show() |
|
158
|
|
|
|
|
159
|
|
|
# plot their histogram |
|
160
|
|
|
fig, ax = plt.subplots() |
|
|
|
|
|
|
161
|
|
|
bins = "auto" # bins = np.int(y.size**0.5/1.5) # |
|
162
|
|
|
rang = (np.percentile(y, 0), np.percentile(y, 99)) |
|
163
|
|
|
ax.hist(y, density=True, color="b", alpha=0.4, bins=bins, range=rang) |
|
164
|
|
|
ax2 = ax.twinx() |
|
165
|
|
|
bins = "auto" # bins = np.int(y.size**0.5/1.5) # |
|
166
|
|
|
rang = (np.percentile(y2, 0), np.percentile(y2, 99)) |
|
167
|
|
|
ax2.hist(y2, density=True, color="r", alpha=0.4, bins=bins, range=rang) |
|
168
|
|
|
plt.show() |
|
169
|
|
|
|
|
170
|
|
|
# plot their PSD |
|
171
|
|
|
fig, ax = plt.subplots() |
|
|
|
|
|
|
172
|
|
|
ax.psd(y, color="b", lw=1, alpha=0.5) |
|
173
|
|
|
ax2 = ax.twinx() |
|
174
|
|
|
ax2.psd(y2, color="r", lw=1, alpha=0.5) |
|
175
|
|
|
plt.show() |
|
176
|
|
|
|
|
177
|
|
|
def PSD_check_gen(self, fgen=None, ax=None): |
|
|
|
|
|
|
178
|
|
|
"""Description goes here.""" |
|
179
|
|
|
|
|
180
|
|
|
# TODO: make a generic check_gen method |
|
|
|
|
|
|
181
|
|
|
|
|
182
|
|
|
if fgen is None: |
|
183
|
|
|
fgen = self.fgen |
|
184
|
|
|
if ax is None: |
|
185
|
|
|
ax = plt.gca() |
|
186
|
|
|
|
|
187
|
|
|
if fgen == "lc_gen_psd_nft": |
|
188
|
|
|
s2 = lc_gen_psd_nft(self.times, self.values) |
|
|
|
|
|
|
189
|
|
|
elif fgen == "lc_gen_psd_lombscargle": |
|
190
|
|
|
s2 = lc_gen_psd_lombscargle(self.times, self.values) |
|
|
|
|
|
|
191
|
|
|
elif fgen == "lc_gen_psd_fft": |
|
192
|
|
|
s2 = lc_gen_psd_fft(self.values) |
|
|
|
|
|
|
193
|
|
|
else: |
|
194
|
|
|
raise Exception("No valid fgen specified") |
|
195
|
|
|
|
|
196
|
|
|
ax.plot(self.times, self.values, "b-", label="orig", lw=0.5, alpha=0.8) |
|
197
|
|
|
ax.plot(self.times, s2, "r-", label="gen", lw=0.5, alpha=0.8) |
|
198
|
|
|
ax.legend() |
|
199
|
|
|
|
|
200
|
|
|
@staticmethod |
|
201
|
|
|
def pdf(xx, ll, mu): |
|
|
|
|
|
|
202
|
|
|
"""Fit pdf as a curve.""" |
|
203
|
|
|
return (ll * mu) ** (1 + ll) / scipy_special.gamma(1 + ll) * np.exp(-ll * mu / xx) / xx ** (ll + 2) |
|
|
|
|
|
|
204
|
|
|
|
IAA-CSIC /
MUTIS
Loading history...