Commit 76f72e00 authored by JOSSOUD Olivier's avatar JOSSOUD Olivier
Browse files

Transition. ERF fitting.

parent d8814a92
Pipeline #152097 passed with stages
in 1 minute and 38 seconds
import pandas as pd
import csaps
import numpy as np
from scipy.special import erf
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import datetime
......@@ -218,6 +220,51 @@ def get_transition_duration_gauss_sigma(df: pd.DataFrame,
plt.show()
def get_transition_duration_erf(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()
def erfunc(x, C1, C2, z0, sigma):
return (C1/2) * (1 + erf((x - z0) / (sigma * np.sqrt(2)))) + C2
xdata = df.index.view(int) / 10 ** 9
xdata = xdata - xdata[0]
xdata = xdata[1:]
ydata = df[col_value_name].tolist()
ydata = ydata[1:]
# Initial guess
C1_0 = ydata[-1] - ydata[0] # Step's amplitude
C2_0 = ydata[0] # Base line level
z0_0 = xdata[-1] - xdata[0] # Shift instant, in seconds from first xdata value
sigma_0 = 10
params, extras = curve_fit(erfunc, xdata, ydata, p0=[C1_0, C2_0, z0_0, sigma_0])
print(params)
if debug_plot:
plt.title('Erf fit, sigma=' + "{:.2f}".format(params[3]) + "s")
plt.xlabel('Time (s)')
plt.ylabel(col_value_name)
plt.plot(xdata, ydata, 'k', label="Raw")
plt.plot(xdata, erfunc(xdata, *params), label="Erf fit")
plt.axhline(y=params[1], color="yellow", linestyle="--", label="Step's plateau")
plt.axhline(y=params[1] + params[0], color="yellow", linestyle="--")
plt.axvline(x=params[2], color="red", linestyle="--", label="Step's 'center'")
plt.axvline(x=params[2]-params[3], color="green", linestyle="--", label="Step's width")
plt.axvline(x=params[2]+params[3], color="green", linestyle="--")
plt.legend()
plt.show()
return params
# def __detect_transition__(df: pd.DataFrame, col_value_name: str = "value") -> pd.DataFrame:
# #
# # /!\ DOES NOT WORK!
......@@ -266,6 +313,13 @@ if __name__ == '__main__':
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_erf(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, 6, 49, 33, tzinfo=datetime.timezone.utc),
# datetime.datetime(2021, 9, 10, 6, 57, 14, tzinfo=datetime.timezone.utc),
......@@ -276,16 +330,16 @@ if __name__ == '__main__':
# 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_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
# )
......
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