From 9e5a232dcf62fe5b6f19b7de4fadf1f999350d32 Mon Sep 17 00:00:00 2001 From: guillaumegrolleron Date: Tue, 13 Feb 2024 10:15:29 +0100 Subject: [PATCH] Fixing multiprocessing issue and minor bugfixes (#108) * first commit to change fork to spawn multiprocessing * implementation of SPE fit with multiprocessing "spawn" * change default args * bugfix for photostat method : -SPE_result read is a generator -typo -useless line to change events_per_slice to none but it's read only * update notebooks + new one about photostat * bugfix : CamerayDisplay is expecting a sorted image with pixels id * update user scripts * remove useless line --------- Co-authored-by: guillaume.grolleron --- notebooks/tool_implementation/tuto_SPE.py | 15 +- .../tool_implementation/tuto_photostat.py | 99 +++++ notebooks/tool_implementation/tuto_tools.py | 105 ++---- .../tuto_waveforms_charges_extraction.py | 88 ++++- src/nectarchain/display/display.py | 27 +- .../calibration/gain/photostat_makers.py | 41 +- .../component/photostatistic_component.py | 2 +- .../makers/component/spe/spe_algorithm.py | 354 ++++++++++-------- src/nectarchain/makers/core.py | 10 + .../ggrolleron/gain_PhotoStat_computation.py | 22 +- .../gain_SPEfit_combined_computation.py | 21 +- .../ggrolleron/gain_SPEfit_computation.py | 28 +- .../ggrolleron/load_wfs_compute_charge.py | 23 +- 13 files changed, 509 insertions(+), 326 deletions(-) create mode 100644 notebooks/tool_implementation/tuto_photostat.py diff --git a/notebooks/tool_implementation/tuto_SPE.py b/notebooks/tool_implementation/tuto_SPE.py index f35c9c1f..19d020c5 100644 --- a/notebooks/tool_implementation/tuto_SPE.py +++ b/notebooks/tool_implementation/tuto_SPE.py @@ -29,25 +29,28 @@ ) # %% -run_number = 3936 +run_number = 3942 # %% +os.environ["NECTARCAMDATA"] +# %% # !ls -lh $NECTARCAMDATA/runs/* # %% -tool = FlatFieldSPENominalStdNectarCAMCalibrationTool( +tool = FlatFieldSPEHHVStdNectarCAMCalibrationTool( progress_bar=True, method="LocalPeakWindowSum", extractor_kwargs={"window_width": 12, "window_shift": 4}, multiproc=True, - nproc=10, + nproc=2, run_number=run_number, - max_events=10000, + max_events=1000, log_level=20, - reload_events=True, + reload_events=False, # events_per_slice = 200, - asked_pixels_id=[52, 48], + overwrite=True, + asked_pixels_id=[52, 48, 78, 94], output_path=pathlib.Path(os.environ.get("NECTARCAMDATA", "/tmp")) / "tutorials/" / f"SPEfit_{run_number}.h5", diff --git a/notebooks/tool_implementation/tuto_photostat.py b/notebooks/tool_implementation/tuto_photostat.py new file mode 100644 index 00000000..719cad81 --- /dev/null +++ b/notebooks/tool_implementation/tuto_photostat.py @@ -0,0 +1,99 @@ +# --- +# jupyter: +# jupytext: +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.14.6 +# kernelspec: +# display_name: nectarchain +# language: python +# name: python3 +# --- + +# %% [markdown] +# # Tutorial for gain computation with the Photo-statitic method + +# %% +import logging +import os +import pathlib + +logging.basicConfig( + format="%(asctime)s %(name)s %(levelname)s %(message)s", level=logging.INFO +) +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers + +import matplotlib.pyplot as plt + +from nectarchain.data.management import DataManagement +from nectarchain.makers.calibration import PhotoStatisticNectarCAMCalibrationTool +from nectarchain.makers.extractor.utils import CtapipeExtractor + +# %% +extractor_kwargs = {"window_width": 12, "window_shift": 4} + +method = "LocalPeakWindowSum" +HHV_run_number = 3942 + +Ped_run_number = 3938 +FF_run_number = 3937 + +# %% +str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str(extractor_kwargs) +path = DataManagement.find_SPE_HHV( + run_number=HHV_run_number, + method=method, + str_extractor_kwargs=str_extractor_kwargs, +) +if len(path) == 1: + log.info( + f"{path[0]} found associated to HHV run {HHV_run_number}, method {method} and extractor kwargs {str_extractor_kwargs}" + ) +else: + _text = f"no file found in $NECTARCAM_DATA/../SPEfit associated to HHV run {HHV_run_number}, method {method} and extractor kwargs {str_extractor_kwargs}" + log.error(_text) + raise FileNotFoundError(_text) + +# %% [markdown] +# WARNING : for now you can't split the event loop in slice for the Photo-statistic method, however in case of the charges havn't been computed on disk, the loop over events will only store the charge, therefore memory errors should happen rarely + +# %% +tool = PhotoStatisticNectarCAMCalibrationTool( + progress_bar=True, + run_number=FF_run_number, + Ped_run_number=Ped_run_number, + SPE_result=path[0], + method="LocalPeakWindowSum", + extractor_kwargs={"window_width": 12, "window_shift": 4}, + max_events=10000, + log_level=20, + reload_events=False, + overwrite=True, + output_path=pathlib.Path(os.environ.get("NECTARCAMDATA", "/tmp")) + / "tutorials/" + / f"Photostat_FF{FF_run_number}_Ped{Ped_run_number}.h5", +) +tool + +# %% +tool.initialize() + +# %% +tool.setup() + +# %% +tool.start() + +# %% +output = tool.finish(return_output_component=True) +output + +# %% +plt.plot(output[0].pixels_id, output[0].high_gain.T[0]) +plt.xlabel("pixels_id") +plt.ylabel("high gain") + +# %% diff --git a/notebooks/tool_implementation/tuto_tools.py b/notebooks/tool_implementation/tuto_tools.py index 595e4a11..dfe51f1b 100644 --- a/notebooks/tool_implementation/tuto_tools.py +++ b/notebooks/tool_implementation/tuto_tools.py @@ -13,17 +13,17 @@ # --- # %% [markdown] -# # Tutorial to use ctapipe tools and component +# # Tutorial to use ctapipe tools and component # %% [markdown] # ## Context # %% [markdown] -# Tool and Component are 2 modules of ctapipe, Tool is a high level module to analyse raw data (fits.fz files). This module use Component to perform computation on the raw data. Basically, we can create a class (MyTool) which inherits of Tool, where we can define 2 Component (Comp_A and Comp_B). Thus, with an instance of MyTool, we can loop over event within raw data, and for each event apply sucessively Comp_A, Comp_B. +# Tool and Component are 2 modules of ctapipe, Tool is a high level module to analyse raw data (fits.fz files). This module use Component to perform computation on the raw data. Basically, we can create a class (MyTool) which inherits of Tool, where we can define 2 Component (Comp_A and Comp_B). Thus, with an instance of MyTool, we can loop over event within raw data, and for each event apply sucessively Comp_A, Comp_B. # # A ctapipe tutorial is accessible here : https://ctapipe.readthedocs.io/en/stable/auto_examples/core/command_line_tools.html#sphx-glr-auto-examples-core-command-line-tools-py # -# You can find documentation of ctapipe Tool and Component : +# You can find documentation of ctapipe Tool and Component : # # https://ctapipe.readthedocs.io/en/stable/api-reference/tools/index.html # @@ -32,35 +32,32 @@ # # Within nectarchain, we implemented within the nectarchain.makers module both a top level Tool and Component from which all the nectarchain Component and Tool should inherit. # -# In this tutorial, we explain quickly how we can use Tool and Component to develop the nectarchain software, as an example there is the implementation of a PedestalTool which extract pedestal +# In this tutorial, we explain quickly how we can use Tool and Component to develop the nectarchain software, as an example there is the implementation of a PedestalTool which extract pedestal # %% [markdown] # ### Imports -# %% -import numpy as np -import pathlib import os +import pathlib import matplotlib.pyplot as plt -from ctapipe_io_nectarcam.containers import NectarCAMDataContainer -from ctapipe_io_nectarcam import constants - -from ctapipe.core.traits import ComponentNameList +# %% +import numpy as np from ctapipe.containers import Field -from ctapipe.core.traits import Integer from ctapipe.core import Component +from ctapipe.core.traits import ComponentNameList, Integer from ctapipe.io import HDF5TableReader +from ctapipe_io_nectarcam import constants +from ctapipe_io_nectarcam.containers import NectarCAMDataContainer - -from nectarchain.makers import EventsLoopNectarCAMCalibrationTool -from nectarchain.makers.component import NectarCAMComponent, ArrayDataComponent from nectarchain.data.container import ( - NectarCAMContainer, ArrayDataContainer, + NectarCAMContainer, TriggerMapContainer, ) +from nectarchain.makers import EventsLoopNectarCAMCalibrationTool +from nectarchain.makers.component import ArrayDataComponent, NectarCAMComponent from nectarchain.utils import ComponentUtils # %% @@ -72,13 +69,14 @@ # %% [markdown] -# The only thing to add to to fill the componentList field, which contains the names of the component to be apply on events. +# The only thing to add to to fill the componentList field, which contains the names of the component to be apply on events. # -# Then we will define a very simple component to compute the pedestal of each events. +# Then we will define a very simple component to compute the pedestal of each events. # %% [markdown] # ### Definition of container to store extracted data on disk + # %% class MyContainer(NectarCAMContainer): run_number = Field( @@ -109,6 +107,7 @@ class MyContainer(NectarCAMContainer): # %% [markdown] # ### Definition of our Component + # %% class MyComp(NectarCAMComponent): window_shift = Integer( @@ -147,7 +146,7 @@ def __call__(self, event: NectarCAMDataContainer, *args, **kwargs): #####THE JOB IS HERE###### for i, pedestal in enumerate([self.__pedestal_hg, self.__pedestal_lg]): index_peak = np.argmax(wfs[i]) - signal_start = index_peak - self.window_shiftclasses + signal_start = index_peak - self.window_shift signal_stop = index_peak + self.window_width - self.window_shift if signal_start < 0: signal_stop = self.window_width @@ -188,7 +187,8 @@ def finish(self): # ### Definition of our Tool # %% [markdown] -# Now we can define out Tool, we have just to add our component "MyComp" in the ComponentList : +# Now we can define out Tool, we have just to add our component "MyComp" in the ComponentList : + # %% def get_valid_component(): @@ -204,39 +204,6 @@ class MyTool(EventsLoopNectarCAMCalibrationTool): help="List of Component names to be apply, the order will be respected", ).tag(config=True) - ####### THIS PART IS NEEDED NOW ##### - - # This part uis needed because the component defined outside of the nectarchain - # module is not accessible from the module (with eval(componentName)). - # This issus should be fixed soon. - def __new__(cls, *args, **kwargs): - _cls = super(EventsLoopNectarCAMCalibrationTool, cls).__new__( - cls, *args, **kwargs - ) - - for componentName in _cls.componentsList: - configurable_traits = ComponentUtils.get_configurable_traits( - eval(componentName) - ) - _cls.add_traits(**configurable_traits) - _cls.aliases.update( - {key: f"{componentName}.{key}" for key in configurable_traits.keys()} - ) - return _cls - - def _get_provided_component_kwargs(self, componentName: str): - component_kwargs = ComponentUtils.get_configurable_traits(eval(componentName)) - output_component_kwargs = {} - for key in component_kwargs.keys(): - if hasattr(self, key): - output_component_kwargs[key] = getattr(self, key) - return output_component_kwargs - - ##################################################################### - - # def __init__(self,*args,**kwargs) : - # super().__init__(*args,**kwargs) - def _init_output_path(self): if self.max_events is None: filename = f"{self.name}_run{self.run_number}.h5" @@ -249,7 +216,12 @@ def _init_output_path(self): # %% tool = MyTool( - progress_bar=True, run_number=4943, max_events=500, log_level=20, window_width=14 + progress_bar=True, + run_number=4943, + max_events=500, + log_level=20, + window_width=14, + overwrite=True, ) # %% @@ -265,19 +237,19 @@ def _init_output_path(self): tool.initialize() # %% [markdown] -# Then to setup, it will in particular setup the Components : +# Then to setup, it will in particular setup the Components : # %% tool.setup() # %% [markdown] -# The following command will just start the tool and apply components looping over events +# The following command will just start the tool and apply components looping over events # %% tool.start() # %% [markdown] -# Then, we finish the tool, behind thius command the component will be finilized and will create an output container whiich will be written on disk and can be returned +# Then, we finish the tool, behind thius command the component will be finilized and will create an output container whiich will be written on disk and can be returned # %% output = tool.finish(return_output_component=True)[0] @@ -286,19 +258,19 @@ def _init_output_path(self): output # %% [markdown] -# The following file has been written : +# The following file has been written : # %% # !ls -lh $NECTARCAMDATA/tutorials # %% [markdown] -# The shape of pedestal is (n_events,n_pixels) +# The shape of pedestal is (n_events,n_pixels) # %% output.pedestal_hg.shape # %% [markdown] -# To have a look to a random pixel pedestal evolution : +# To have a look to a random pixel pedestal evolution : # %% fix, ax = plt.subplots(1, 1) @@ -313,12 +285,14 @@ def _init_output_path(self): ax.set_xlabel("time [ms]") # %% [markdown] -# If you want to load container thereafter : +# If you want to load container thereafter : # %% -container_loaded = MyContainer._container_from_hdf5( - f"{os.environ.get('NECTARCAMDATA','/tmp')}/tutorials/PedestalTutoNectarCAM_run4943_maxevents500.h5", - MyContainer, +container_loaded = next( + MyContainer._container_from_hdf5( + f"{os.environ.get('NECTARCAMDATA','/tmp')}/tutorials/PedestalTutoNectarCAM_run4943_maxevents500.h5", + MyContainer, + ) ) container_loaded.validate() container_loaded @@ -328,7 +302,7 @@ def _init_output_path(self): # %% [markdown] # An argument that are implemented in EventsLoopNectarCAMCalibrationTool is the 'event_per_slice', this argument allows to split all the events within the raw data fits.fz file in slices. It allows to, for each slice, loop over events and write container on disk. This mechanism allows to save RAM. -# The resulting hdf5 file that is written on disk , can be easily loaded thereafter. There is only one hdf5 file for the whole run, which is a mapping between slices and containers filled by computed quantity from components. +# The resulting hdf5 file that is written on disk , can be easily loaded thereafter. There is only one hdf5 file for the whole run, which is a mapping between slices and containers filled by computed quantity from components. # %% tool = MyTool( @@ -360,6 +334,7 @@ def _init_output_path(self): # container_loaded = ArrayDataContainer._container_from_hdf5(f"{os.environ.get('NECTARCAMDATA','/tmp')}/tutorials/PedestalTutoNectarCAM_run4943_maxevents2000.h5",MyContainer) # container_loaded + # %% def read_hdf5_sliced(path): container = MyContainer() diff --git a/notebooks/tool_implementation/tuto_waveforms_charges_extraction.py b/notebooks/tool_implementation/tuto_waveforms_charges_extraction.py index b77b00ea..7b4e7829 100644 --- a/notebooks/tool_implementation/tuto_waveforms_charges_extraction.py +++ b/notebooks/tool_implementation/tuto_waveforms_charges_extraction.py @@ -13,19 +13,25 @@ # --- # %% +import logging +import os + +logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") +log = logging.getLogger(__name__) +log.handlers = logging.getLogger("__main__").handlers +import pathlib + +import numpy as np +from ctapipe.containers import EventType +from ctapipe.coordinates import EngineeringCameraFrame +from ctapipe.instrument import CameraGeometry + +from nectarchain.display import ContainerDisplay from nectarchain.makers import ( ChargesNectarCAMCalibrationTool, WaveformsNectarCAMCalibrationTool, ) from nectarchain.makers.component import get_valid_component -from nectarchain.data.container import ( - ChargesContainers, - ChargesContainer, - WaveformsContainer, - WaveformsContainers, -) -from ctapipe.io import HDF5TableReader -from ctapipe.containers import EventType # %% get_valid_component() @@ -35,10 +41,38 @@ # %% tool = WaveformsNectarCAMCalibrationTool( - progress_bar=True, run_number=run_number, max_events=500, log_level=20 + progress_bar=True, + run_number=run_number, + max_events=500, + log_level=20, + overwrite=True, + output_path=pathlib.Path(os.environ.get("NECTARCAMDATA", "/tmp")) + / "tutorials/" + / f"Waveforms_{run_number}.h5", ) tool +# %% +tool.initialize() + +# %% +tool.setup() + +# %% +tool.start() + +# %% +w_output = tool.finish(return_output_component=True)[0] +w_output + +# %% +w_output.containers + +# %% + +# %% [markdown] +# ### Now for the charge extraction + # %% tool = ChargesNectarCAMCalibrationTool( progress_bar=True, @@ -47,23 +81,49 @@ run_number=run_number, max_events=500, log_level=20, + from_computed_waveforms=False, + overwrite=True, + output_path=pathlib.Path(os.environ.get("NECTARCAMDATA", "/tmp")) + / "tutorials/" + / f"Charges_{run_number}_LocalPeakWindowSum_12-4.h5", ) tool # %% tool.initialize() - -# %% tool.setup() # %% tool.start() # %% -output = tool.finish(return_output_component=True)[0] -output +c_output = tool.finish(return_output_component=True)[0] +c_output # %% -output.containers +c_output.containers + +# %% [markdown] +# ### Display + +# %% +geom = CameraGeometry.from_name("NectarCam-003").transform_to(EngineeringCameraFrame()) + +# %% +image = w_output.containers[EventType.FLATFIELD].wfs_hg.sum(axis=2)[23] +max_id = w_output.containers[EventType.FLATFIELD].pixels_id[np.argmax(image)] +max_id + +# %% +disp = ContainerDisplay.display( + w_output.containers[EventType.FLATFIELD], evt=23, geometry=geom +) +# disp.highlight_pixels(max_id,color = 'r',linewidth = 3) #to check the validity of geometry +disp.show() + +# %% +ContainerDisplay.display( + c_output.containers[EventType.FLATFIELD], evt=23, geometry=geom +) # %% diff --git a/src/nectarchain/display/display.py b/src/nectarchain/display/display.py index f9dc2f2d..9df8a759 100644 --- a/src/nectarchain/display/display.py +++ b/src/nectarchain/display/display.py @@ -10,6 +10,8 @@ log = logging.getLogger(__name__) log.handlers = logging.getLogger("__main__").handlers +import numpy as np + class ContainerDisplay(ABC): @staticmethod @@ -25,11 +27,34 @@ def display(container: ArrayDataContainer, evt, geometry, cmap="gnuplot2"): """ if isinstance(container, ChargesContainer): image = container.charges_hg + pixels_id = container.pixels_id elif isinstance(container, WaveformsContainer): image = container.wfs_hg.sum(axis=2) + pixels_id = container.pixels_id else: - log.warning("container can't be displayed") + log.error( + "container can't be displayed, must be a ChargesContainer or a WaveformsContainer" + ) + raise Exception( + "container can't be displayed, must be a ChargesContainer or a WaveformsContainer" + ) + + highlighten_pixels = np.array([], dtype=int) + if geometry.pix_id.value.shape[0] != image.shape[1]: + mask = np.array([_id in pixels_id for _id in geometry.pix_id.value]) + missing_pixels = np.array(geometry.pix_id.value[~mask], dtype=int) + + missing_values = np.empty((image.shape[0], missing_pixels.shape[0])) + missing_values.fill(np.nan) + highlighten_pixels = np.concatenate((highlighten_pixels, missing_pixels)) + + image = np.concatenate((missing_values, image), axis=1) + pixels_id = np.concatenate((missing_pixels, pixels_id)) + sort_index = [np.where(pixels_id == pix)[0][0] for pix in geometry.pix_id.value] + image = image.T[sort_index].T + disp = CameraDisplay(geometry=geometry, image=image[evt], cmap=cmap) + disp.highlight_pixels(highlighten_pixels, color="r", linewidth=2) disp.add_colorbar() return disp diff --git a/src/nectarchain/makers/calibration/gain/photostat_makers.py b/src/nectarchain/makers/calibration/gain/photostat_makers.py index 3249170c..302fac51 100644 --- a/src/nectarchain/makers/calibration/gain/photostat_makers.py +++ b/src/nectarchain/makers/calibration/gain/photostat_makers.py @@ -11,7 +11,7 @@ from ctapipe.containers import Container, EventType from ctapipe.core.traits import ComponentNameList, Integer, Path -from ....data.container import ChargesContainer +from ....data.container import ChargesContainer, ChargesContainers from ....data.container.core import NectarCAMContainer, merge_map_ArrayDataContainer from ....data.management import DataManagement from ...component import ArrayDataComponent, NectarCAMComponent @@ -30,11 +30,11 @@ class PhotoStatisticNectarCAMCalibrationTool(GainNectarCAMCalibrationTool): default_value=["PhotoStatisticNectarCAMComponent"], help="List of Component names to be apply, the order will be respected", ).tag(config=True) - run_number = Integer(help="the FF run number to be treated", default_value=-1).tag( - config=True - ) + run_number = Integer( + help="the flat-field run number to be treated", default_value=-1 + ).tag(config=True) Ped_run_number = Integer( - help="the FF run number to be treated", default_value=-1 + help="the pedestal run number to be treated", default_value=-1 ).tag(config=True) run_file = Path( @@ -54,27 +54,6 @@ class PhotoStatisticNectarCAMCalibrationTool(GainNectarCAMCalibrationTool): def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) - str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str( - self.extractor_kwargs - ) - if not (self.reload_events): - FF_files = DataManagement.find_charges( - run_number=self.run_number, - method=self.method, - str_extractor_kwargs=str_extractor_kwargs, - max_events=self.max_events, - ) - Ped_files = DataManagement.find_charges( - run_number=self.run_number, - max_events=self.max_events, - ) - - if len(FF_files) == 1 and len(Ped_files) == 1: - log.warning( - "You asked events_per_slice but you don't want to reload events and a charges file is on disk for FF and Pedestals, then events_per_slice is set to None" - ) - self.events_per_slice = None - def _init_output_path(self): str_extractor_kwargs = CtapipeExtractor.get_extractor_kwargs_str( self.extractor_kwargs @@ -115,7 +94,7 @@ def start( max_events=self.max_events, ) Ped_files = DataManagement.find_charges( - run_number=self.run_number, + run_number=self.Ped_run_number, max_events=self.max_events, ) if self.reload_events or len(FF_files) != 1 or len(Ped_files) != 1: @@ -137,9 +116,9 @@ def start( else: self.log.info(f"reading computed charge from FF file {FF_files[0]}") chargesContainers = ChargesContainer.from_hdf5(FF_files[0]) - if isinstance(chargesContainers, NectarCAMContainer): + if isinstance(chargesContainers, ChargesContainer): self.components[0]._FF_chargesContainers = chargesContainers - elif isinstance(list(chargesContainers.containers.keys())[0], EventType): + elif isinstance(chargesContainers, ChargesContainers): self.log.debug("merging along TriggerType") self.components[0]._FF_chargesContainers = merge_map_ArrayDataContainer( chargesContainers @@ -156,9 +135,9 @@ def start( self.log.info(f"reading computed charge from Ped file {Ped_files[0]}") chargesContainers = ChargesContainer.from_hdf5(Ped_files[0]) - if isinstance(chargesContainers, NectarCAMContainer): + if isinstance(chargesContainers, ChargesContainer): self.components[0]._Ped_chargesContainers = chargesContainers - elif isinstance(list(chargesContainers.containers.keys())[0], EventType): + elif isinstance(chargesContainers, ChargesContainers): self.log.debug("merging along TriggerType") self.components[ 0 diff --git a/src/nectarchain/makers/component/photostatistic_component.py b/src/nectarchain/makers/component/photostatistic_component.py index 2a1d6e64..0b73e38d 100644 --- a/src/nectarchain/makers/component/photostatistic_component.py +++ b/src/nectarchain/makers/component/photostatistic_component.py @@ -121,7 +121,7 @@ def finish(self, *args, **kwargs): self._Ped_chargesContainers = merge_map_ArrayDataContainer( self._Ped_chargesContainers ) - nectarGainSPEresult = SPEfitContainer.from_hdf5(self.SPE_result) + nectarGainSPEresult = next(SPEfitContainer.from_hdf5(self.SPE_result)) photo_stat = eval(self.PhotoStatAlgorithm).create_from_chargesContainer( FFcharge=self._FF_chargesContainers, Pedcharge=self._Ped_chargesContainers, diff --git a/src/nectarchain/makers/component/spe/spe_algorithm.py b/src/nectarchain/makers/component/spe/spe_algorithm.py index 6c656a96..c1c6607e 100644 --- a/src/nectarchain/makers/component/spe/spe_algorithm.py +++ b/src/nectarchain/makers/component/spe/spe_algorithm.py @@ -1,10 +1,12 @@ import logging +import sys logging.basicConfig(format="%(asctime)s %(name)s %(levelname)s %(message)s") log = logging.getLogger(__name__) log.handlers = logging.getLogger("__main__").handlers import copy +import multiprocessing as mp import os import time from inspect import signature @@ -44,6 +46,61 @@ ] +class ContextFit: + def __init__( + self, + _class, + minuitParameters_array: np.ndarray, + charge: np.ndarray, + counts: np.ndarray, + ): + self._class = _class + self.minuitParameters_array = minuitParameters_array + self.charge = charge + self.counts = counts + + def __enter__(self): + init_processes( + self._class, self.minuitParameters_array, self.charge, self.counts + ) + + def __exit__(self, type, value, traceback): + del globals()["_minuitParameters_array"] + del globals()["_charge"] + del globals()["_counts"] + del globals()["_chi2"] + + +def init_processes( + _class, minuitParameters_array: np.ndarray, charge: np.ndarray, counts: np.ndarray +): + """ + Initialize each process in the process pool with global variable fit_array. + """ + global _minuitParameters_array + global _charge + global _counts + global _chi2 + _minuitParameters_array = minuitParameters_array + _charge = charge + _counts = counts + + def chi2(pedestal, pp, luminosity, resolution, mean, n, pedestalWidth, index): + return _class._NG_Likelihood_Chi2( + pp, + resolution, + mean, + n, + pedestal, + pedestalWidth, + luminosity, + _charge[int(index)].data[~_charge[int(index)].mask], + _counts[int(index)].data[~_charge[int(index)].mask], + ) + + _chi2 = chi2 + + class SPEalgorithm(Component): Windows_lenght = Integer( 40, @@ -271,11 +328,14 @@ def _get_mean_gaussian_fit( ax.set_ylabel("Events", size=15) ax.legend(fontsize=7) os.makedirs( - f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/", + f"{os.environ.get('NECTARCHAIN_LOG', '/tmp')}/{os.getpid()}/figures/", exist_ok=True, ) + log.info( + f'figures of initialization of parameters will be accessible at {os.environ.get("NECTARCHAIN_LOG","/tmp")}/{os.getpid()}' + ) fig.savefig( - f"{os.environ.get('NECTARCHAIN_LOG')}/{os.getpid()}/figures/initialization_pedestal_pixel{pixel_id}_{os.getpid()}.pdf" + f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/figures/initialization_pedestal_pixel{pixel_id}_{os.getpid()}.pdf" ) fig.clf() plt.close(fig) @@ -344,32 +404,6 @@ def _get_mean_gaussian_fit( return coeff, coeff_mean -''' - def _update_table_from_parameters(self) -> None: - """ - Update the result table based on the parameters of the FlatFieldSPEMaker class. - This method adds columns to the table for each parameter and its corresponding error. - """ - - for param in self._parameters.parameters: - if not (param.name in self._results.keys()): - self._results.add_column( - Column( - data=np.empty((self.npixels), dtype=np.float64), - name=param.name, - unit=param.unit, - ) - ) - self._results.add_column( - Column( - data=np.empty((self.npixels), dtype=np.float64), - name=f"{param.name}_error", - unit=param.unit, - ) - ) -''' - - class SPEnominalalgorithm(SPEalgorithm): parameters_file = Unicode( "parameters_SPEnominal.yaml", @@ -377,8 +411,6 @@ class SPEnominalalgorithm(SPEalgorithm): help="The name of the SPE fit parameters file", ).tag(config=True) - __fit_array = None - tol = Float( 1e-1, help="The tolerance used for minuit", @@ -481,7 +513,9 @@ def _counts(self): return self.__counts # methods - def _fill_results_table_from_dict(self, dico: dict, pixels_id: np.ndarray) -> None: + def _fill_results_table_from_dict( + self, dico: dict, pixels_id: np.ndarray, return_fit_array: bool = True + ) -> None: """ Populates the results table with fit values and errors for each pixel based on the dictionary provided as input. @@ -493,37 +527,51 @@ def _fill_results_table_from_dict(self, dico: dict, pixels_id: np.ndarray) -> No None """ ########NEED TO BE OPTIMIZED!!!########### - chi2_sig = signature(__class__.cost(self._charge, self._counts)) + chi2_sig = signature(_chi2) + if return_fit_array: + fit_array = np.empty(len(pixels_id), dtype=np.object_) for i in range(len(pixels_id)): values = dico[i].get(f"values_{i}", None) errors = dico[i].get(f"errors_{i}", None) - if not ((values is None) or (errors is None)): + fit_status = dico[i].get(f"fit_status_{i}", None) + if not ((values is None) or (errors is None) or (fit_status is None)): + if return_fit_array: + fit_array[i] = fit_status index = np.argmax(self._results.pixels_id == pixels_id[i]) if len(values) != len(chi2_sig.parameters): e = Exception( "the size out the minuit output parameters values array does not fit the signature of the minimized cost function" ) - self.log.error(e, exc_info=True) + log.error(e, exc_info=True) raise e for j, key in enumerate(chi2_sig.parameters): - self._results[key][index][0] = values[j] - self._results[key][index][1] = errors[j] - self._results[key][index][2] = errors[j] - if key == "mean": - self._results.high_gain[index][0] = values[j] - self._results.high_gain[index][1] = errors[j] - self._results.high_gain[index][2] = errors[j] - self._results.is_valid[index] = True - self._results.likelihood[index] = __class__.__fit_array[i].fcn( - __class__.__fit_array[i].values + if key != "index": + self._results[key][index][0] = values[j] + self._results[key][index][1] = errors[j] + self._results[key][index][2] = errors[j] + if key == "mean": + self._results.high_gain[index][0] = values[j] + self._results.high_gain[index][1] = errors[j] + self._results.high_gain[index][2] = errors[j] + self._results.is_valid[index] = ( + fit_status["is_valid"] + and not (fit_status["has_parameters_at_limit"]) + and fit_status["has_valid_parameters"] ) + if fit_status["has_reached_call_limit"]: + self.log.warning( + f"The minuit fit for pixel {pixels_id[i]} reached the call limit" + ) + self._results.likelihood[index] = fit_status["values"] ndof = ( self._counts.data[index][~self._counts.mask[index]].shape[0] - - __class__.__fit_array[i].nfit + - fit_status["nfit"] ) self._results.p_value[index] = Statistics.chi2_pvalue( - ndof, __class__.__fit_array[i].fcn(__class__.__fit_array[i].values) + ndof, fit_status["values"] ) + if return_fit_array: + return fit_array @staticmethod def _NG_Likelihood_Chi2( @@ -554,6 +602,12 @@ def _NG_Likelihood_Chi2( Returns: float: The chi-square value. """ + if not (kwargs.get("ntotalPE", False)): + for i in range(1000): + if gammainc(i + 1, lum) < 1e-5: + ntotalPE = i + break + kwargs.update({"ntotalPE": ntotalPE}) pdf = MPE2(charge, pp, res, mu2, n, muped, sigped, lum, **kwargs) # log.debug(f"pdf : {np.sum(pdf)}") Ntot = np.sum(counts) @@ -564,49 +618,8 @@ def _NG_Likelihood_Chi2( ) # 2 times faster return Lik - @staticmethod - def cost(charge: np.ndarray, counts: np.ndarray): - """ - Defines a function called Chi2 that calculates the chi-square value using the _NG_Likelihood_Chi2 method. - Parameters: - charge (np.ndarray): An array of charge values. - counts (np.ndarray): An array of count values. - Returns: - function: The Chi2 function. - """ - - def Chi2( - pedestal: float, - pp: float, - luminosity: float, - resolution: float, - mean: float, - n: float, - pedestalWidth: float, - ): - # assert not(np.isnan(pp) or np.isnan(resolution) or np.isnan(mean) or np.isnan(n) or np.isnan(pedestal) or np.isnan(pedestalWidth) or np.isnan(luminosity)) - for i in range(1000): - if gammainc(i + 1, luminosity) < 1e-5: - ntotalPE = i - break - kwargs = {"ntotalPE": ntotalPE} - return __class__._NG_Likelihood_Chi2( - pp, - resolution, - mean, - n, - pedestal, - pedestalWidth, - luminosity, - charge, - counts, - **kwargs, - ) - - return Chi2 - # @njit(parallel=True,nopython = True) - def _make_fit_array_from_parameters( + def _make_minuitParameters_array_from_parameters( self, pixels_id: np.ndarray = None, **kwargs ) -> np.ndarray: """ @@ -624,7 +637,7 @@ def _make_fit_array_from_parameters( else: npix = len(pixels_id) - fit_array = np.empty((npix), dtype=np.object_) + minuitParameters_array = np.empty((npix), dtype=np.object_) for i, _id in enumerate(pixels_id): index = np.where(self.pixels_id == _id)[0][0] @@ -635,36 +648,20 @@ def _make_fit_array_from_parameters( pixel_id=_id, **kwargs, ) + index_parameter = Parameter(name="index", value=index, frozen=True) + parameters.append(index_parameter) + minuitParameters = UtilsMinuit.make_minuit_par_kwargs(parameters) - minuit_kwargs = { - parname: minuitParameters["values"][parname] - for parname in minuitParameters["values"] - } - self.log.info(f"creation of fit instance for pixel: {_id}") - fit_array[i] = Minuit( - __class__.cost( - self._charge[index].data[~self._charge[index].mask], - self._counts[index].data[~self._charge[index].mask], - ), - **minuit_kwargs, - ) - self.log.debug("fit created") - fit_array[i].errordef = Minuit.LIKELIHOOD - fit_array[i].strategy = 0 - fit_array[i].tol = self.tol - fit_array[i].print_level = 1 - fit_array[i].throw_nan = True - UtilsMinuit.set_minuit_parameters_limits_and_errors( - fit_array[i], minuitParameters - ) + minuitParameters_array[i] = minuitParameters + # self.log.debug(fit_array[i].values) # self.log.debug(fit_array[i].limits) # self.log.debug(fit_array[i].fixed) - return fit_array + return minuitParameters_array @staticmethod - def run_fit(i: int) -> dict: + def run_fit(i: int, tol: float) -> dict: """ Perform a fit on a specific pixel using the Minuit package. @@ -675,13 +672,43 @@ def run_fit(i: int) -> dict: dict: A dictionary containing the fit values and errors for the specified pixel. The keys are "values_i" and "errors_i", where "i" is the index of the pixel. """ + log = logging.getLogger(__name__) + log.setLevel(logging.INFO) + log.info("Starting") - __class__.__fit_array[i].migrad() - __class__.__fit_array[i].hesse() - _values = np.array([params.value for params in __class__.__fit_array[i].params]) - _errors = np.array([params.error for params in __class__.__fit_array[i].params]) + minuit_kwargs = { + parname: _minuitParameters_array[i]["values"][parname] + for parname in _minuitParameters_array[i]["values"] + } + log.info(f"creation of fit instance for pixel: {minuit_kwargs['index']}") + fit = Minuit(_chi2, **minuit_kwargs) + fit.errordef = Minuit.LIKELIHOOD + fit.strategy = 0 + fit.tol = tol + fit.print_level = 1 + fit.throw_nan = True + UtilsMinuit.set_minuit_parameters_limits_and_errors( + fit, _minuitParameters_array[i] + ) + log.info("fit created") + fit.migrad() + fit.hesse() + _values = np.array([params.value for params in fit.params]) + _errors = np.array([params.error for params in fit.params]) + _fit_status = { + "is_valid": fit.fmin.is_valid, + "has_valid_parameters": fit.fmin.has_valid_parameters, + "has_parameters_at_limit": fit.fmin.has_parameters_at_limit, + "has_reached_call_limit": fit.fmin.has_reached_call_limit, + "nfit": fit.nfit, + "values": fit.fcn(fit.values), + } log.info("Finished") - return {f"values_{i}": _values, f"errors_{i}": _errors} + return { + f"values_{i}": _values, + f"errors_{i}": _errors, + f"fit_status_{i}": _fit_status, + } def run( self, @@ -709,52 +736,69 @@ def run( self.log.warning("The asked pixels id are all out of the data") return None else: - self.log.info("creation of the fits instance array") - __class__.__fit_array = self._make_fit_array_from_parameters( + self.log.info("creation of the minuit parameters array") + minuitParameters_array = self._make_minuitParameters_array_from_parameters( pixels_id=pixels_id, display=display, **kwargs ) self.log.info("running fits") - if self.multiproc: - nproc = kwargs.get("nproc", self.nproc) - chunksize = kwargs.get( - "chunksize", - max(self.chunksize, npix // (nproc * 10)), - ) - self.log.info(f"pooling with nproc {nproc}, chunksize {chunksize}") + with ContextFit( + __class__, minuitParameters_array, self._charge, self._counts + ): + if self.multiproc: + try: + mp.set_start_method("spawn", force=True) + self.log.info("spawned multiproc") + except RuntimeError: + pass + nproc = kwargs.get("nproc", self.nproc) + chunksize = kwargs.get( + "chunksize", + max(self.chunksize, npix // (nproc * 10)), + ) + self.log.info(f"pooling with nproc {nproc}, chunksize {chunksize}") + + t = time.time() + with Pool( + nproc, + initializer=init_processes, + initargs=( + __class__, + minuitParameters_array, + self._charge, + self._counts, + ), + ) as pool: + self.log.info( + f"time to create the Pool is {time.time() - t:.2e} sec" + ) + result = pool.starmap_async( + __class__.run_fit, + [(i, self.tol) for i in range(npix)], + chunksize=chunksize, + ) + result.wait() + try: + res = result.get() + except Exception as e: + log.error(e, exc_info=True) + raise e + self.log.info( + f"total time for multiproc with starmap_async execution is {time.time() - t:.2e} sec" + ) - t = time.time() - with Pool(nproc) as pool: - result = pool.starmap_async( - __class__.run_fit, - [(i,) for i in range(npix)], - chunksize=chunksize, + else: + self.log.info("running in mono-cpu") + t = time.time() + res = [__class__.run_fit(i, self.tol) for i in range(npix)] + self.log.info( + f"time for singleproc execution is {time.time() - t:.2e} sec" ) - result.wait() - try: - res = result.get() - except Exception as e: - self.log.error(e, exc_info=True) - raise e - self.log.debug(str(res)) - self.log.info( - f"time for multiproc with starmap_async execution is {time.time() - t:.2e} sec" - ) - else: - self.log.info("running in mono-cpu") - t = time.time() - res = [__class__.run_fit(i) for i in range(npix)] - self.log.debug(res) - self.log.info( - f"time for singleproc execution is {time.time() - t:.2e} sec" + self.log.info("filling result table from fits results") + output = self._fill_results_table_from_dict( + res, pixels_id, return_fit_array=True ) - self.log.info("filling result table from fits results") - self._fill_results_table_from_dict(res, pixels_id) - - output = copy.copy(__class__.__fit_array) - __class__.__fit_array = None - if display: self.log.info("plotting") t = time.time() diff --git a/src/nectarchain/makers/core.py b/src/nectarchain/makers/core.py index dec9a9f3..8759034d 100644 --- a/src/nectarchain/makers/core.py +++ b/src/nectarchain/makers/core.py @@ -7,6 +7,7 @@ import copy import os import pathlib +from datetime import datetime import numpy as np from ctapipe.containers import Container @@ -24,6 +25,7 @@ from ctapipe.io.datawriter import DATA_MODEL_VERSION from tables.exceptions import HDF5ExtError from tqdm.auto import tqdm +from traitlets import default from ..data import DataManagement from ..data.container import LightNectarCAMEventSource @@ -48,6 +50,14 @@ class BaseNectarCAMCalibrationTool(Tool): help="show progress bar during processing", default_value=False ).tag(config=True) + @default("provenance_log") + def _default_provenance_log(self): + return f"{os.environ.get('NECTARCHAIN_LOG', '/tmp')}/{self.name}_{os.getpid()}_{datetime.now()}.provenance.log" + + @default("log_file") + def _default_log_file(self): + return f"{os.environ.get('NECTARCHAIN_LOG', '/tmp')}/{self.name}_{os.getpid()}_{datetime.now()}.log" + @staticmethod def load_run( run_number: int, max_events: int = None, run_file: str = None diff --git a/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.py b/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.py index 582fb06f..f3991be0 100644 --- a/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.py +++ b/src/nectarchain/user_scripts/ggrolleron/gain_PhotoStat_computation.py @@ -3,17 +3,6 @@ import os import sys from pathlib import Path - -# to quiet numba -os.makedirs(os.environ.get("NECTARCHAIN_LOG"), exist_ok=True) -logging.getLogger("numba").setLevel(logging.WARNING) -logging.basicConfig( - format="%(asctime)s %(name)s %(levelname)s %(message)s", - level=logging.DEBUG, - filename=f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{Path(__file__).stem}_{os.getpid()}.log", -) -log = logging.getLogger(__name__) - import argparse import copy import glob @@ -172,12 +161,19 @@ def main( if __name__ == "__main__": + import logging + + # to quiet numba + logging.getLogger("numba").setLevel(logging.WARNING) + args = parser.parse_args() kwargs = copy.deepcopy(vars(args)) kwargs["log_level"] = args.verbosity - - os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/figures") + os.makedirs( + f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/figures", + exist_ok=True, + ) logging.basicConfig( format="%(asctime)s %(name)s %(levelname)s %(message)s", force=True, diff --git a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_combined_computation.py b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_combined_computation.py index a43706bd..6cf83469 100644 --- a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_combined_computation.py +++ b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_combined_computation.py @@ -4,17 +4,6 @@ import sys import time from pathlib import Path - -# to quiet numba -os.makedirs(os.environ.get("NECTARCHAIN_LOG"), exist_ok=True) -logging.getLogger("numba").setLevel(logging.WARNING) -logging.basicConfig( - format="%(asctime)s %(name)s %(levelname)s %(message)s", - level=logging.DEBUG, - filename=f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{Path(__file__).stem}_{os.getpid()}.log", -) -log = logging.getLogger(__name__) - import argparse import copy import glob @@ -193,13 +182,19 @@ def main( if __name__ == "__main__": + import logging + + # to quiet numba + logging.getLogger("numba").setLevel(logging.WARNING) t = time.time() args = parser.parse_args() kwargs = copy.deepcopy(vars(args)) kwargs["log_level"] = args.verbosity - - os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/figures") + os.makedirs( + f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/figures", + exist_ok=True, + ) logging.basicConfig( format="%(asctime)s %(name)s %(levelname)s %(message)s", force=True, diff --git a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.py b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.py index 0a1cc23f..a24a5ce3 100644 --- a/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.py +++ b/src/nectarchain/user_scripts/ggrolleron/gain_SPEfit_computation.py @@ -1,23 +1,12 @@ +import argparse +import copy import json -import logging import os import sys import time +from multiprocessing import freeze_support from pathlib import Path -# to quiet numba -os.makedirs(os.environ.get("NECTARCHAIN_LOG"), exist_ok=True) -logging.getLogger("numba").setLevel(logging.WARNING) -logging.basicConfig( - format="%(asctime)s %(name)s %(levelname)s %(message)s", - level=logging.DEBUG, - filename=f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{Path(__file__).stem}_{os.getpid()}.log", -) -log = logging.getLogger(__name__) - -import argparse -import copy - from nectarchain.makers.calibration import ( FlatFieldSPEHHVNectarCAMCalibrationTool, FlatFieldSPEHHVStdNectarCAMCalibrationTool, @@ -187,6 +176,11 @@ def main( if __name__ == "__main__": + import logging + + # to quiet numba + logging.getLogger("numba").setLevel(logging.WARNING) + t = time.time() args = parser.parse_args() @@ -194,7 +188,10 @@ def main( kwargs["log_level"] = args.verbosity - os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/figures") + os.makedirs( + f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/figures", + exist_ok=True, + ) logging.basicConfig( format="%(asctime)s %(name)s %(levelname)s %(message)s", force=True, @@ -217,6 +214,7 @@ def main( kwargs.pop("figpath") kwargs.pop("display") kwargs.pop("HHV") + kwargs.pop("free_pp_n") log.info(f"arguments passed to main are : {kwargs}") main(log=log, **kwargs) diff --git a/src/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.py b/src/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.py index d6e9939d..7da5dd92 100644 --- a/src/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.py +++ b/src/nectarchain/user_scripts/ggrolleron/load_wfs_compute_charge.py @@ -1,18 +1,11 @@ import argparse import json -import logging import os import sys from pathlib import Path os.makedirs(os.environ.get("NECTARCHAIN_LOG"), exist_ok=True) -logging.getLogger("numba").setLevel(logging.WARNING) -logging.basicConfig( - format="%(asctime)s %(name)s %(levelname)s %(message)s", - level=logging.DEBUG, - filename=f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{Path(__file__).stem}_{os.getpid()}.log", -) -log = logging.getLogger(__name__) + import copy from nectarchain.data.container import ( @@ -78,13 +71,13 @@ "SlidingWindowMaxSum", "TwoPassWindowSum", ], - default="LocalPeakWindowSum", + default="FullWaveformSum", help="charge extractor method", type=str, ) parser.add_argument( "--extractor_kwargs", - default={"window_width": 16, "window_shift": 4}, + default={}, help="charge extractor kwargs", type=json.loads, ) @@ -162,7 +155,6 @@ def main( from_computed_waveforms=True, **charges_kwargs, ) - tool.setup() tool.start() tool.finish() @@ -171,6 +163,10 @@ def main( if __name__ == "__main__": + import logging + + # to quiet numba + logging.getLogger("numba").setLevel(logging.WARNING) # run of interest # spe_run_number = [2633,2634,3784] # ff_run_number = [2608] @@ -182,7 +178,10 @@ def main( kwargs["log_level"] = args.verbosity - os.makedirs(f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/figures") + os.makedirs( + f"{os.environ.get('NECTARCHAIN_LOG','/tmp')}/{os.getpid()}/figures", + exist_ok=True, + ) logging.basicConfig( format="%(asctime)s %(name)s %(levelname)s %(message)s", force=True,