Coverage for src / lstautorta / hdf5_data_check.py: 0%
159 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-08-10 11:56 +0000
1import datetime
2import time
3from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser
4from pathlib import Path
5from threading import Lock
6from typing import Annotated, NamedTuple
8import numpy as np
9import pandas as pd
10import seaborn as sns
11from annotated_types import Ge
12from matplotlib import pyplot as plt
13from matplotlib.backends.backend_pdf import PdfPages
14from scipy import ndimage
15from tqdm import tqdm
17from lstautorta.config.hdf5_data_check import HDFDataCheckConfiguration
18from lstautorta.utils.hdf5 import pd_read_hdf_with_retry
21class FileProvenance(NamedTuple):
22 obs_id: Annotated[int, Ge(0)]
23 line_idx: Annotated[int, Ge(0)]
24 worker_node: str
27# def parse_file_provenance(cta_data_file_path: str, process_idx_to_nodelist: Dict[int, str]) -> FileProvenance:
28# # get the string after "..._id_", and split on "_" to get only this part of the file path.
29# obs_id = int(cta_data_file_path.split("obs_id_", 1)[1].split("_", 1)[0])
30# line_idx = int(cta_data_file_path.split("line_idx_", 1)[1].split("_", 1)[0])
31# return FileProvenance(obs_id=obs_id, line_idx=line_idx, worker_node=process_idx_to_nodelist[line_idx])
34def parse_file_provenance(cta_data_file_path: Path, process_idx_to_nodelist: dict[int, str]) -> FileProvenance:
35 # dl1_v06_${obs_id}_${run_id}_${tel_id}_${processIndex}_${threadIndex}_${fileIndex}.h5
36 split_ = cta_data_file_path.stem.split("_")
37 obs_id = int(split_[3])
38 line_idx = int(split_[5])
39 return FileProvenance(obs_id=obs_id, line_idx=line_idx, worker_node=process_idx_to_nodelist[line_idx])
42def add_debuging_fields_dl1_df(
43 dl1_df: pd.DataFrame,
44 obs_id: int,
45 line_idx: int,
46 worker_node: str,
47 worker_node_categories: list[str],
48 start=None,
49 end=None,
50 step=None,
51) -> None:
52 # Some events have a trigger time "from the future" that mess with event_id_diff
53 # so set the event id diff of all of these future events and the following event as well to 0
54 # prepend 1st value to get diff = 0 for 1st elem
55 event_id_diffs = np.diff(dl1_df["event_id"], prepend=dl1_df["event_id"][0])
56 event_id_diffs[
57 ndimage.binary_dilation(dl1_df["trigger_time"] > pd.to_datetime("21000101"), np.array([0, 1, 1]))
58 ] = 0
59 dl1_df.loc[slice(start, end, step), "event_id_diff"] = event_id_diffs
60 dl1_df.loc[slice(start, end, step), "obs_id"] = obs_id
61 dl1_df.loc[slice(start, end, step), "line_idx"] = line_idx
62 dl1_df.loc[slice(start, end, step), "worker_node"] = pd.Categorical(
63 [worker_node] * dl1_df.shape[0], categories=worker_node_categories
64 )
67def timestamp_to_datetime(df: pd.DataFrame, columns: list[str]):
68 for col in columns:
69 df[col] = pd.to_datetime(df[col], unit="s")
72class HDFDataCheck:
73 def __init__(
74 self,
75 data_check_config: HDFDataCheckConfiguration,
76 process_idx_to_nodelist: dict[int, str],
77 ):
78 self.conf = data_check_config
80 self.last_plot_datetime = datetime.datetime.now()
81 self.plot_refresh_interval_timedelta = datetime.timedelta(seconds=data_check_config.plot_refresh_interval_s)
83 self.process_idx_to_nodelist = process_idx_to_nodelist
85 self.data_lock = Lock()
86 self.hdf_data = None
87 self.current_length = 0
89 def init_db(self, dtypes: pd.Series):
90 # Allocate the big df that will accumulate all the data to plot.
91 # size should be more than enough to store an entire run data
92 # we first create an empty df with the right length
93 # then we create the columns with a default value, according to the dtype
94 # (pandas automatically fills missing values with NaNs, which changes the type of the column to float
95 # if we were to try to create the df with only the dtype and an example data.)
96 self.hdf_data = pd.DataFrame(index=np.arange(self.conf.size_storage_dataframe))
97 for col in dtypes.index:
98 col_dtype = dtypes[col]
99 if isinstance(col_dtype, pd.CategoricalDtype):
100 self.hdf_data[col] = pd.Categorical(
101 [col_dtype.categories[0]] * self.hdf_data.shape[0], categories=col_dtype.categories
102 )
103 elif np.issubdtype(col_dtype, np.datetime64):
104 self.hdf_data[col] = np.datetime64(0, "ns")
105 elif np.issubdtype(col_dtype, bool):
106 self.hdf_data[col] = False
107 elif np.issubdtype(col_dtype, np.integer):
108 self.hdf_data[col] = col_dtype.type(np.iinfo(col_dtype).max)
109 elif np.issubdtype(col_dtype, np.inexact):
110 self.hdf_data[col] = col_dtype.type(np.nan)
111 else:
112 raise ValueError(
113 f"Can not set default value of column {col} with dtype {col_dtype}! Unsupported dtype!"
114 )
116 def read_hdf5_data(self, hdf5_file_path: Path):
117 file_provenance = parse_file_provenance(hdf5_file_path, self.process_idx_to_nodelist)
118 # Read the hdf5 data directly in pandas DataFrame
119 # Alternatively, we could read with pytables and then copy the columns in the big df with table.col("col_name")
120 # but moving to pandas straight away allows to perform all operations on df
121 file_df = pd_read_hdf_with_retry(
122 hdf5_file_path,
123 key=self.conf.data_path_in_hdf5,
124 mode="r",
125 nb_tries=self.conf.hdf5_open_nb_retries,
126 retry_wait_time_s=self.conf.hdf5_open_wait_time_s,
127 retry_on_os_error=True,
128 )
129 if self.conf.column_to_datetime is not None:
130 timestamp_to_datetime(file_df, self.conf.column_to_datetime)
131 add_debuging_fields_dl1_df(
132 file_df,
133 file_provenance.obs_id,
134 file_provenance.line_idx,
135 file_provenance.worker_node,
136 list(self.process_idx_to_nodelist.values()),
137 )
139 file_nb_events = file_df.shape[0]
141 with self.data_lock:
142 # If the big df is not initialized: do it now that we know the dtype
143 if self.hdf_data is None:
144 self.init_db(file_df.dtypes)
146 # copy the file df in the big df and update length. Do not use loc but iloc !!
147 self.hdf_data.iloc[self.current_length : self.current_length + file_nb_events] = file_df
149 self.current_length += file_nb_events
151 def data_check_new_file(self, hdf5_file_path: Path, png_dir_path: Path, plots_pdf_path: Path):
152 self.read_hdf5_data(hdf5_file_path)
154 if datetime.datetime.now() - self.last_plot_datetime > self.plot_refresh_interval_timedelta:
155 self.last_plot_datetime = datetime.datetime.now()
156 with self.data_lock:
157 self.plot(png_dir_path, plots_pdf_path)
159 def plot(self, png_dir_path: Path, plots_pdf_path: Path):
160 sns.set_theme()
162 png_dir_path.mkdir(exist_ok=True)
164 # legend order
165 legend_order = sorted(self.hdf_data.loc[:, "worker_node"].dtype.categories)
166 # masks
167 mask_event_quality = self.hdf_data["event_quality"] == 0
168 mask_is_good_event = self.hdf_data["is_good_event"] == 1
169 mask_future_events = self.hdf_data["trigger_time"] <= pd.to_datetime("21000101")
171 # TODO:
172 # - repair histogram plots with numpy_x fcts or numpy_y
173 # - apply masks on plots that have e31 values ...
174 # - find out why seaborn is so slow (do only histograms ?)
175 # - implement tests
177 with PdfPages(plots_pdf_path) as pdf_pages:
178 for plt_png_path, plt_conf in self.conf.plots_config.items():
179 fig, ax = plt.subplots(1, 1, figsize=(12, 8))
181 # select data based on masks and current length
182 fig_mask = np.ones((self.hdf_data.shape[0],), dtype=bool)
183 fig_mask[self.current_length :] = 0 # do not select anything over current length
184 if plt_conf.mask_event_quality:
185 fig_mask &= mask_event_quality
186 if plt_conf.mask_future_events:
187 fig_mask &= mask_future_events
188 if plt_conf.mask_is_good_event:
189 fig_mask &= mask_is_good_event
190 data = self.hdf_data.loc[fig_mask, :]
191 # Apply groupby operation if requested
192 if plt_conf.groupby is not None:
193 data = getattr(
194 data.groupby(
195 [
196 item if isinstance(item, str) else pd.Grouper(**item)
197 for item in plt_conf.groupby.groupby
198 ]
199 ),
200 plt_conf.groupby.computation[0],
201 **(plt_conf.groupby.kwargs if plt_conf.groupby.kwargs is not None else {}),
202 )(**(plt_conf.groupby.computation[1] if plt_conf.groupby.computation[1] is not None else {}))
204 # Apply function on x or y data if requested
205 for field, fct_conf in zip(["x", "y"], [plt_conf.np_function_x, plt_conf.np_function_y]):
206 if fct_conf is not None:
207 try:
208 plt_conf.kwargs[field] = getattr(np, fct_conf[0])(
209 data.loc[:, plt_conf.kwargs[field]],
210 **(fct_conf[1] if fct_conf[1] is not None else {}),
211 )
212 except AttributeError as e:
213 raise ValueError(
214 f"Could not apply function {fct_conf[0]} with kwargs {fct_conf[1]} on {plt_conf.kwargs[field]}"
215 ) from e
217 if plt_conf.kind == "plot":
218 plot_fct = sns.lineplot
219 # do not sort unless user requested it (lineplot sorted is quite confusing)
220 plt_conf.kwargs["sort"] = plt_conf.kwargs.get("sort", False)
221 plt_conf.kwargs["style"] = plt_conf.kwargs.get("style", "worker_node")
222 plt_conf.kwargs["style_order"] = plt_conf.kwargs.get("style_order", legend_order)
223 elif plt_conf.kind == "hist":
224 plot_fct = sns.histplot
225 plt_conf.kwargs["multiple"] = plt_conf.kwargs.get("multiple", "dodge")
226 # plt_conf.kwargs["align"] = plt_conf.kwargs.get("align", "center") # not supported by seaborn which already supplies the 'align' arg
227 plt_conf.kwargs["shrink"] = plt_conf.kwargs.get("shrink", 0.8)
228 else:
229 raise ValueError(f"Plot kind {plt_conf.kind} not supported")
230 # TODO: fix histo intensity
231 # fix xticks time: show second
232 # fix histogram column placement
233 plot_fct(
234 **plt_conf.kwargs,
235 data=data,
236 hue="worker_node",
237 hue_order=legend_order,
238 ax=ax,
239 )
240 if plt_conf.title is not None:
241 ax.set_title(plt_conf.title)
242 if plt_conf.xlabel is not None:
243 ax.set_xlabel(plt_conf.xlabel)
244 if plt_conf.ylabel is not None:
245 ax.set_ylabel(plt_conf.ylabel)
247 if len(ax.get_xticklabels()[0].get_text()) >= 6:
248 ax.set_xticks(ax.get_xticks(), ax.get_xticklabels(), rotation=45, ha="right")
250 fig.tight_layout()
251 fig.savefig(png_dir_path / plt_png_path, bbox_inches="tight")
252 pdf_pages.savefig(fig)
253 plt.close(fig)
255 def write_df_and_reset_data(self, data_dir: Path, df_write_path: Path):
256 self.hdf_data.loc[: self.current_length].to_hdf(df_write_path, key=self.conf.write_df_store_key, mode="w")
257 self.current_length = 0
260def main():
261 parser = ArgumentParser(
262 description="Data check plotting util of LST AUTO RTA: "
263 "run this entrypoint to reproduce offline the plots generated during the night",
264 formatter_class=ArgumentDefaultsHelpFormatter,
265 )
266 parser.add_argument(
267 "-c",
268 "--data_check_config",
269 type=str,
270 required=True,
271 dest="config",
272 help="Path to the data check configuration file",
273 )
274 parser.add_argument(
275 "-d",
276 "--data_dir",
277 type=str,
278 required=True,
279 dest="data_dir",
280 help="Path to the directory containing the data files, to run the data check on",
281 )
282 parser.add_argument(
283 "-n",
284 "--night_nodes",
285 type=str,
286 nargs="+",
287 dest="nodes",
288 help="List of worker nodes used during the night. Eg: cp12 cp15 cp42 cp58",
289 )
290 parser.add_argument(
291 "-p",
292 "--png_dir_path",
293 type=str,
294 required=True,
295 dest="png_dir_path",
296 help="Path to the directory where to write the plots as png images.",
297 )
298 parser.add_argument(
299 "-f", "--pdf_path", type=str, required=True, dest="pdf_path", help="Path to the output pdf with the plots."
300 )
301 args = parser.parse_args()
303 start_time = time.process_time()
304 with open(Path(args.config)) as config_f:
305 config = HDFDataCheckConfiguration.model_validate_json(config_f.read())
306 config_loaded_time = time.process_time()
307 print(f"Configuration loading time: {config_loaded_time - start_time:.2f}s")
309 data_checker = HDFDataCheck(config, {i: n for i, n in enumerate(args.nodes)})
310 instantiation_time = time.process_time()
311 print(f"Instantiation time: {instantiation_time - config_loaded_time:.2f}s")
313 file_list = [
314 p for p in Path(args.data_dir).glob("*") if p.is_file() and p.suffix == ".h5" and p.stem.startswith("dl1_v06")
315 ]
316 data_checker.read_hdf5_data(file_list[0])
317 init_and_load_first_file_time = time.process_time()
318 print(f"Big df init and load 1st file time: {init_and_load_first_file_time - instantiation_time:.2f}s")
320 for data_file in tqdm(file_list[1:]):
321 data_checker.read_hdf5_data(data_file)
322 load_all_files_time = time.process_time()
323 print(f"Load remaining files: {load_all_files_time - init_and_load_first_file_time:.2f}s")
325 data_checker.plot(Path(args.png_dir_path), Path(args.pdf_path))
326 plot_time = time.process_time()
327 print(f"Plot time: {plot_time - load_all_files_time:.2f}s")
330if __name__ == "__main__":
331 main()