Skip to content

GitLab

Commit 002f9c29 authored by JOSSOUD Olivier's avatar JOSSOUD Olivier
Browse files

Transition. Gaussian fit. Works but not very effective.

parent 7f5b5592
Hide whitespace changes
Inline Side-by-side
Showing with 137 additions and 49 deletions
cfatools/tsanalyser/transition.py
View file @ 002f9c29
......@@ -154,62 +154,150 @@ def get_transition_duration(df: pd.DataFrame,
return duration_sec
def __detect_transition__(df: pd.DataFrame, col_value_name: str = "value") -> pd.DataFrame:
#
# /!\ DOES NOT WORK!
#
df["diff"] = df[col_value_name].diff()
def get_transition_duration_gauss_sigma(df: pd.DataFrame,
first_datetime: datetime.datetime,
last_datetime: datetime.datetime,
col_value_name: str = "value",
debug_plot: bool = False) -> float:
# Restrict the data to the transition area
df = df[(df.index >= first_datetime) & (df.index <= last_datetime)]
df = df.copy()
sp = csaps.UnivariateCubicSmoothingSpline(df.index, df[col_value_name], smooth=1e-29)
df["smooth"] = sp(df.index)
first_date = datetime.datetime(2019, 11, 4, 9, 26, 47, tzinfo=datetime.timezone.utc)
last_date = datetime.datetime(2019, 11, 4, 9, 36, 4, tzinfo=datetime.timezone.utc)
df = df[(df.index >= first_date) & (df.index < last_date)]
timestamp = [(dt.to_pydatetime() - datetime.datetime(1970, 1, 1, tzinfo=datetime.timezone.utc)).total_seconds() for dt in df.index]
model = piecewise(timestamp, df[col_value_name], min_stop_frac=0.03)
plt.plot(df.index, df[col_value_name], '-')
for segment in model.segments:
plt.axvline(x=datetime.datetime.fromtimestamp(segment.start_t, tz=datetime.timezone.utc), color="r")
plt.axvline(x=datetime.datetime.fromtimestamp(segment.end_t, tz=datetime.timezone.utc), color="b")
df["deriv"] = df[col_value_name].diff() / df.index.to_series(keep_tz=True).diff().dt.total_seconds()
df["deriv_smooth"] = df["smooth"].diff() / df.index.to_series(keep_tz=True).diff().dt.total_seconds()
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
def gauss(x, H, A, x0, sigma):
return H + A * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))
def gauss_fit(x, y):
mean = sum(x * y) / sum(y)
sigma = np.sqrt(sum(y * (x - mean) ** 2) / sum(y))
popt, pcov = curve_fit(gauss, x, y, p0=[min(y), max(y), mean, sigma])
return popt
# generate simulated data
xdata = df.index.view(int) / 10 ** 9
# ydata_perfect = gauss(xdata, 20, 5, 6, 1)
ydata = df["deriv_smooth"].tolist()
xdata = xdata[1:]
ydata = ydata[1:]
# ydata = [-y for y in ydata]
H, A, x0, sigma = gauss_fit(xdata, ydata)
FWHM = 2.35482 * sigma
print('The offset of the gaussian baseline is', H)
print('The center of the gaussian fit is', x0)
print('The sigma of the gaussian fit is', sigma)
print('The maximum intensity of the gaussian fit is', H + A)
print('The Amplitude of the gaussian fit is', A)
print('The FWHM of the gaussian fit is', FWHM)
plt.plot(xdata, ydata, 'ko', label='data')
# plt.plot(xdata, ydata_perfect, '-k', label='data (without_noise)')
plt.plot(xdata, gauss(xdata, *gauss_fit(xdata, ydata)), '--r', label='fit')
plt.legend()
plt.title('Gaussian fit, $f(x) = A e^{(-(x-x_0)^2/(2sigma^2))}$, sigma=' + "{:.2f}".format(sigma) + "s")
plt.xlabel('Time (s)')
plt.ylabel('Delta_18_16 deriv')
plt.show()
base_level = df[df["datetime"] >= start_datetime].iloc[0]["rol10_mean"]
current_plateau_end = df[(df["datetime"] >= start_datetime) & (abs(df["diff"]) > threshold)].iloc[0]["datetime"]
next_plateau_start = df[(df["datetime"] >= current_plateau_end) & (abs(df["diff"]) < threshold)].iloc[0]["datetime"]
next_plateau_end = df[(df["datetime"] >= next_plateau_start) & (abs(df["diff"]) > threshold)].iloc[0]["datetime"]
subset_df = df.loc[(df["datetime"] >= start_datetime) & (df["datetime"] <= next_plateau_end)]
np.random.seed(1234)
x = np.linspace(-5., 5., 25)
y = np.exp(-(x / 2.5) ** 2) + (np.random.rand(25) - 0.2) * 0.3
sp = csaps.UnivariateCubicSmoothingSpline(x, y, smooth=0.85)
if debug_plot:
plt.plot(df.index, df[col_value_name], 'g-')
plt.plot(df.index, df["smooth"], 'r-')
# plt.axvline(x=transition_time, color="b")
# plt.title('Transition time =' + transition_time.strftime("%H:%M:%S"))
plt.show()
xs = np.linspace(x[0], x[-1], 150)
ys = sp(xs)
plt.plot(x, y, 'o', xs, ys, '-')
for segment in model.segments:
plt.axvline(x=segment.start_t)
plt.axvline(x=segment.end_t)
plt.show()
# def __detect_transition__(df: pd.DataFrame, col_value_name: str = "value") -> pd.DataFrame:
# #
# # /!\ DOES NOT WORK!
# #
# df["diff"] = df[col_value_name].diff()
# sp = csaps.UnivariateCubicSmoothingSpline(df.index, df[col_value_name], smooth=1e-29)
# df["smooth"] = sp(df.index)
#
# first_date = datetime.datetime(2019, 11, 4, 9, 26, 47, tzinfo=datetime.timezone.utc)
# last_date = datetime.datetime(2019, 11, 4, 9, 36, 4, tzinfo=datetime.timezone.utc)
# df = df[(df.index >= first_date) & (df.index < last_date)]
#
# timestamp = [(dt.to_pydatetime() - datetime.datetime(1970, 1, 1, tzinfo=datetime.timezone.utc)).total_seconds() for dt in df.index]
# model = piecewise(timestamp, df[col_value_name], min_stop_frac=0.03)
#
# plt.plot(df.index, df[col_value_name], '-')
# for segment in model.segments:
# plt.axvline(x=datetime.datetime.fromtimestamp(segment.start_t, tz=datetime.timezone.utc), color="r")
# plt.axvline(x=datetime.datetime.fromtimestamp(segment.end_t, tz=datetime.timezone.utc), color="b")
# plt.show()
#
# base_level = df[df["datetime"] >= start_datetime].iloc[0]["rol10_mean"]
# current_plateau_end = df[(df["datetime"] >= start_datetime) & (abs(df["diff"]) > threshold)].iloc[0]["datetime"]
# next_plateau_start = df[(df["datetime"] >= current_plateau_end) & (abs(df["diff"]) < threshold)].iloc[0]["datetime"]
# next_plateau_end = df[(df["datetime"] >= next_plateau_start) & (abs(df["diff"]) > threshold)].iloc[0]["datetime"]
# subset_df = df.loc[(df["datetime"] >= start_datetime) & (df["datetime"] <= next_plateau_end)]
#
# np.random.seed(1234)
#
# x = np.linspace(-5., 5., 25)
# y = np.exp(-(x / 2.5) ** 2) + (np.random.rand(25) - 0.2) * 0.3
#
# sp = csaps.UnivariateCubicSmoothingSpline(x, y, smooth=0.85)
#
# xs = np.linspace(x[0], x[-1], 150)
# ys = sp(xs)
#
# plt.plot(x, y, 'o', xs, ys, '-')
# for segment in model.segments:
# plt.axvline(x=segment.start_t)
# plt.axvline(x=segment.end_t)
# plt.show()
if __name__ == '__main__':
log = LogReader("/home/dryas/cfa")
df = log.get_timeseries("20191104_picarro_delay", "PICARRO_periodic", "Delta_18_16")
get_transition_time(df,
datetime.datetime(2019, 11, 4, 9, 26, 47, tzinfo=datetime.timezone.utc),
datetime.datetime(2019, 11, 4, 9, 36, 24, tzinfo=datetime.timezone.utc),
debug_plot=True)
get_transition_duration(df,
datetime.datetime(2019, 11, 4, 9, 26, 47, tzinfo=datetime.timezone.utc),
datetime.datetime(2019, 11, 4, 9, 30, 17, tzinfo=datetime.timezone.utc),
datetime.datetime(2019, 11, 4, 9, 33, 0, tzinfo=datetime.timezone.utc),
datetime.datetime(2019, 11, 4, 9, 36, 24, tzinfo=datetime.timezone.utc),
debug_plot=True)
detect_plateau(df)
detect_transition(df)
\ No newline at end of file
from cfatools.logreader.dataset import DatasetReader
log = DatasetReader("/homel/ojossoud/_temp/data_cfa", "20210910_calib_isotopes")
df = log.get_timeseries("PICARRO_periodic", "Delta_18_16")
# get_transition_duration_gauss_sigma(df,
# datetime.datetime(2021, 9, 10, 6, 49, 33, tzinfo=datetime.timezone.utc),
# datetime.datetime(2021, 9, 10, 6, 57, 14, tzinfo=datetime.timezone.utc),
# debug_plot=True
# )
# get_transition_duration_gauss_sigma(df,
# datetime.datetime(2021, 9, 10, 7, 3, 0, tzinfo=datetime.timezone.utc),
# datetime.datetime(2021, 9, 10, 7, 15, 41, tzinfo=datetime.timezone.utc),
# debug_plot=True
# )
get_transition_duration_gauss_sigma(df,
datetime.datetime(2021, 9, 10, 7, 24, 3, tzinfo=datetime.timezone.utc),
datetime.datetime(2021, 9, 10, 7, 26, 43, tzinfo=datetime.timezone.utc),
debug_plot=True
)
get_transition_duration_gauss_sigma(df,
datetime.datetime(2021, 9, 10, 7, 22, 33, tzinfo=datetime.timezone.utc),
datetime.datetime(2021, 9, 10, 7, 28, 13, tzinfo=datetime.timezone.utc),
debug_plot=True
)
# get_transition_time(df,
# datetime.datetime(2019, 11, 4, 9, 26, 47, tzinfo=datetime.timezone.utc),
# datetime.datetime(2019, 11, 4, 9, 36, 24, tzinfo=datetime.timezone.utc),
# debug_plot=True)
# get_transition_duration(df,
# datetime.datetime(2019, 11, 4, 9, 26, 47, tzinfo=datetime.timezone.utc),
# datetime.datetime(2019, 11, 4, 9, 30, 17, tzinfo=datetime.timezone.utc),
# datetime.datetime(2019, 11, 4, 9, 33, 0, tzinfo=datetime.timezone.utc),
# datetime.datetime(2019, 11, 4, 9, 36, 24, tzinfo=datetime.timezone.utc),
# debug_plot=True)
# detect_plateau(df)
# detect_transition(df)
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment