From cd6eb6696a2eb98708ec42448425d067269d554c Mon Sep 17 00:00:00 2001 From: Hiroki Yoneda Date: Thu, 6 Jun 2024 16:14:58 +0200 Subject: [PATCH] Major changes on DataLoader MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The following methods are added: ・event #getter ・keys_bkg_models() ・bkg_model(key) ・exposure_map #getter ・model_axes #getter ・data_axes #getter ・calc_expectation(model) ・calc_T_product(dataspace_histogram) ・calc_bkg_model_product(key, dataspace_histogram) ・calc_likelihood(expectation) --- cosipy/image_deconvolution/data_loader.py | 534 ++++++++++------------ 1 file changed, 247 insertions(+), 287 deletions(-) diff --git a/cosipy/image_deconvolution/data_loader.py b/cosipy/image_deconvolution/data_loader.py index 9e1d87f2..bd5df99f 100644 --- a/cosipy/image_deconvolution/data_loader.py +++ b/cosipy/image_deconvolution/data_loader.py @@ -1,8 +1,10 @@ -import warnings import numpy as np from tqdm.autonotebook import tqdm import astropy.units as u +import logging +logger = logging.getLogger(__name__) + from histpy import Histogram, Axes from cosipy.response import FullDetectorResponse @@ -12,40 +14,73 @@ class DataLoader(object): """ A class to manage data for image analysis, - namely event data, background model, response, coordsys conversion matrix. + namely event data, background models, response, coordsys conversion matrix. Ideally, these data should be input directly to ImageDeconvolution class, but considering their data formats are not fixed, this class is introduced. The purpose of this class is to check the consistency between input data and calculate intermediate files etc. In the future, this class may be removed or hidden in ImageDeconvolution class. """ - def __init__(self): - self.event_dense = None - self.bkg_dense = None - self.full_detector_response = None - self.coordsys_conv_matrix = None - self.mask = None + def __init__(self, name = None): + self._name = name + + # must assign them in DataLoader.load + self._event = None # histpy.Histogram (dense) + self._bkg_models = {} # a dictionary of histpy.Histogram (dense) + self._image_response = None # histpy.Histogram (dense) + self._exposure_map = None + + self._data_axes = None # histpy.Axes + self._model_axes = None # histpy.Axes + + # None if using Galactic CDS, mandotary if using local CDS + self._coordsys_conv_matrix = None + + # optional + self.is_miniDC2_format = False #should be removed in the future - self.is_miniDC2_format = False + @property + def name(self): + return self._name - self.response_on_memory = False + @property + def event(self): + return self._event - self.exposure_map = None + @property + def exposure_map(self): + return self._exposure_map + + @property + def model_axes(self): + return self._model_axes + + @property + def data_axes(self): + return self._data_axes + + def keys_bkg_models(self): + return list(self._bkg_models.keys()) + + def bkg_model(self, key): + return self._bkg_models[key] @classmethod - def load(cls, event_binned_data, bkg_binned_data, rsp, coordsys_conv_matrix, is_miniDC2_format = False): + def load(cls, name, event_binned_data, dict_bkg_binned_data, rsp, coordsys_conv_matrix = None, is_miniDC2_format = False): """ Load data Parameters ---------- + name : str + The name of data event_binned_data : :py:class:`histpy.Histogram` Event histogram - bkg_binned_data : :py:class:`histpy.Histogram` - Background model + dict_bkg_binned_data : dict + Background models as {background_model_name: :py:class:`histpy.Histogram`} rsp : :py:class:`histpy.Histogram` or :py:class:`cosipy.response.FullDetectorResponse` Response - coordsys_conv_matrix : :py:class:`cosipy.image_deconvolution.CoordsysConversionMatrix` + coordsys_conv_matrix : :py:class:`cosipy.image_deconvolution.CoordsysConversionMatrix`, default False Coordsys conversion matrix is_miniDC2_format : bool, default False Whether the file format is for mini-DC2. It will be removed in the future. @@ -53,243 +88,92 @@ def load(cls, event_binned_data, bkg_binned_data, rsp, coordsys_conv_matrix, is_ Returns ------- :py:class:`cosipy.image_deconvolution.DataLoader` - DataLoader instance containing the input data set + An instance of DataLoader containing the input data set """ - new = cls() + new = cls(name) - new.event_dense = event_binned_data.to_dense() + new._event = event_binned_data.to_dense() - new.bkg_dense = bkg_binned_data.to_dense() + new._bkg_models = dict_bkg_binned_data - new.full_detector_response = rsp + for key in new._bkg_models: + if new._bkg_models[key].is_sparse: + new._bkg_models[key] = new._bkg_models[key].to_dense() - new.coordsys_conv_matrix = coordsys_conv_matrix + new._coordsys_conv_matrix = coordsys_conv_matrix new.is_miniDC2_format = is_miniDC2_format - return new - - @classmethod - def load_from_filepath(cls, event_hdf5_filepath = None, event_yaml_filepath = None, - bkg_hdf5_filepath = None, bkg_yaml_filepath = None, - rsp_filepath = None, ccm_filepath = None, - is_miniDC2_format = False): - """ - Load data from file pathes - - Parameters - ---------- - event_hdf5_filepath : str or None, default None - File path of HDF5 file for event histogram. - event_yaml_filepath : str or None, default None - File path of yaml file to read the HDF5 file. - bkg_hdf5_filepath : str or None, default None - File path of HDF5 file for background model. - bkg_yaml_filepath : str or None, default None - File path of yaml file to read the HDF5 file. - rsp_filepath : str or None, default None - File path of the response matrix. - ccm_filepath : str or None, default None - File path of the coordsys conversion matrix. - is_miniDC2_format : bool, default False - Whether the file format is for mini-DC2. should be removed in the future. - - Returns - ------- - :py:class:`cosipy.image_deconvolution.DataLoader` - DataLoader instance containing the input data set - """ - - new = cls() - - new.set_event_from_filepath(event_hdf5_filepath, event_yaml_filepath) - - new.set_bkg_from_filepath(bkg_hdf5_filepath, bkg_yaml_filepath) + if isinstance(rsp, FullDetectorResponse): + logger.info('Loading the response matrix onto your computer memory...') + new._load_full_detector_response_on_memory(rsp, is_miniDC2_format) + logger.info('Finished') + elif isinstance(rsp, Histogram): + new._image_response = rsp - new.set_rsp_from_filepath(rsp_filepath) - - new.set_ccm_from_filepath(ccm_filepath) - - new.is_miniDC2_format = is_miniDC2_format - - return new - - def set_event_from_filepath(self, hdf5_filepath, yaml_filepath): - """ - Load event data from file pathes - - Parameters - ---------- - hdf5_filepath : str - File path of HDF5 file for event histogram. - yaml_filepath : str - File path of yaml file to read the HDF5 file. - """ - - self._event_hdf5_filepath = hdf5_filepath - self._event_yaml_filepath = yaml_filepath - - print(f'... loading event from {hdf5_filepath} and {yaml_filepath}') - - event = BinnedData(self._event_yaml_filepath) - event.load_binned_data_from_hdf5(self._event_hdf5_filepath) - - self.event_dense = event.binned_data.to_dense() - - print("... Done ...") - - def set_bkg_from_filepath(self, hdf5_filepath, yaml_filepath): - """ - Load background model from file pathes - - Parameters - ---------- - hdf5_filepath : str - File path of HDF5 file for background model. - yaml_filepath : str - File path of yaml file to read the HDF5 file. - """ - - self._bkg_hdf5_filepath = hdf5_filepath - self._bkg_yaml_filepath = yaml_filepath - - print(f'... loading background from {hdf5_filepath} and {yaml_filepath}') - - bkg = BinnedData(self._bkg_yaml_filepath) - bkg.load_binned_data_from_hdf5(self._bkg_hdf5_filepath) - - self.bkg_dense = bkg.binned_data.to_dense() - - print("... Done ...") - - def set_rsp_from_filepath(self, filepath): - """ - Load response matrix from file pathes - - Parameters - ---------- - filepath : str - File path of the response matrix. - """ - - self._rsp_filepath = filepath - - print(f'... loading full detector response from {filepath}') - - self.full_detector_response = FullDetectorResponse.open(self._rsp_filepath) - - print("... Done ...") - - def set_ccm_from_filepath(self, filepath): - """ - Load coordsys conversion matrix from file pathes - - Parameters - ---------- - filepath : str - File path of the coordsys conversion matrix. - """ - - self._ccm_filepath = filepath + # We modify the axes in event, bkg_models, response. This is only for DC2. + result = new._modify_axes() + if result == False: + logger.warning('Please rerun after checking the axes of the input histograms') + return - print(f'... loading coordsys conversion matrix from {filepath}') - - self.coordsys_conv_matrix = CoordsysConversionMatrix.open(self._ccm_filepath) - - print("... Done ...") - - def _check_file_registration(self): - """ - Check whether files are loaded. - - Returns - ------- - bool - True if all required files are loaded. - """ - - print(f"... checking the file registration ...") - - if self.event_dense and self.bkg_dense \ - and self.full_detector_response and self.coordsys_conv_matrix: - - print(f" --> pass") - return True - - return False - - def _check_axis_consistency(self): - """ - Check whether the axes of event/background/response are consistent with each other. - - Returns - ------- - bool - True if their axes are consistent. - """ + new._data_axes = new._event.axes - print(f"... checking the axis consistency ...") - - # check the axes of the event/background files - if self.coordsys_conv_matrix.binning_method == 'Time': - axis_name = ['Time', 'Em', 'Phi', 'PsiChi'] - - elif self.coordsys_conv_matrix.binning_method == 'ScAtt': - axis_name = ['ScAtt', 'Em', 'Phi', 'PsiChi'] - - for name in axis_name: - if not self.event_dense.axes[name] == self.bkg_dense.axes[name]: - print(f"Warning: the axis {name} is not consistent between the event and background!") - return False + if new._coordsys_conv_matrix is None: + axes = [new._image_response.axes['NuLambda'], new._image_response.axes['Ei']] + axes[0].label = 'lb' + # The gamma-ray direction of pre-computed response in DC2 is in the galactic coordinate, not in the local coordinate. + # Actually, it is labeled as 'NuLambda'. So I replace it with 'lb'. + new._model_axes = Axes(axes) + else: + new._model_axes = Axes([new._coordsys_conv_matrix.axes['lb'], new._image_response.axes['Ei']]) - # check the axes of the event/response files - axis_name = ['Em', 'Phi', 'PsiChi'] - - for name in axis_name: - if not self.event_dense.axes[name] == self.full_detector_response.axes[name]: - print(f"Warning: the axis {name} is not consistent between the event and response!") - return False + new._calc_exposure_map() - print(f" --> pass") - return True + return new def _modify_axes(self): """ Modify the axes of data. This method will be removed in the future. """ - warnings.warn("Note that _modify_axes() in DataLoader was implemented for a temporary use. It will be removed in the future.", FutureWarning) - warnings.warn("Make sure to perform _modify_axes() only once after the data are loaded.") + logger.warning("Note that _modify_axes() in DataLoader was implemented for a temporary use. It will be removed in the future.") - if self.coordsys_conv_matrix.binning_method == 'Time': + if self._coordsys_conv_matrix is None: + axis_name = ['Em', 'Phi', 'PsiChi'] + + elif self._coordsys_conv_matrix.binning_method == 'Time': axis_name = ['Time', 'Em', 'Phi', 'PsiChi'] - elif self.coordsys_conv_matrix.binning_method == 'ScAtt': + elif self._coordsys_conv_matrix.binning_method == 'ScAtt': axis_name = ['ScAtt', 'Em', 'Phi', 'PsiChi'] for name in axis_name: - print(f"... checking the axis {name} of the event and background files...") + logger.info(f"... checking the axis {name} of the event and background files...") - event_edges, event_unit = self.event_dense.axes[name].edges, self.event_dense.axes[name].unit - bkg_edges, bkg_unit = self.bkg_dense.axes[name].edges, self.bkg_dense.axes[name].unit + event_edges, event_unit = self._event.axes[name].edges, self._event.axes[name].unit - if np.all(event_edges == bkg_edges): - print(f" --> pass (edges)") - else: - print(f"Warning: the edges of the axis {name} are not consistent between the event and background!") - print(f" event : {event_edges}") - print(f" background : {bkg_edges}") - return False + for key in self._bkg_models: - if event_unit == bkg_unit: - print(f" --> pass (unit)") - else: - print(f"Warning: the unit of the axis {name} are not consistent between the event and background!") - print(f" event : {event_unit}") - print(f" background : {bkg_unit}") - return False + bkg_edges, bkg_unit = self._bkg_models[key].axes[name].edges, self._bkg_models[key].axes[name].unit + + if np.all(event_edges == bkg_edges): + logger.info(f" --> pass (edges)") + else: + logger.warning(f"Warning: the edges of the axis {name} are not consistent between the event and the background model {key}!") + logger.warning(f" event : {event_edges}") + logger.warning(f" background : {bkg_edges}") + return False + + if event_unit == bkg_unit: + logger.info(f" --> pass (unit)") + else: + logger.warning(f"Warning: the unit of the axis {name} are not consistent between the event and the background model {key}!") + logger.warning(f" event : {event_unit}") + logger.warning(f" background : {bkg_unit}") + return False # check the axes of the event/response files. # Note that currently (2023-08-29) no unit is stored in the binned data. So only the edges are compared. This should be modified in the future. @@ -298,81 +182,85 @@ def _modify_axes(self): for name in axis_name: - print(f"...checking the axis {name} of the event and response files...") + logger.info(f"...checking the axis {name} of the event and response files...") - event_edges, event_unit = self.event_dense.axes[name].edges, self.event_dense.axes[name].unit - response_edges, response_unit = self.full_detector_response.axes[name].edges, self.full_detector_response.axes[name].unit + event_edges, event_unit = self._event.axes[name].edges, self._event.axes[name].unit + response_edges, response_unit = self._image_response.axes[name].edges, self._image_response.axes[name].unit if type(response_edges) == u.quantity.Quantity and self.is_miniDC2_format == True: response_edges = response_edges.value if np.all(event_edges == response_edges): - print(f" --> pass (edges)") + logger.info(f" --> pass (edges)") else: - print(f"Warning: the edges of the axis {name} are not consistent between the event and background!") - print(f" event : {event_edges}") - print(f" response : {response_edges}") + logger.warning(f"Warning: the edges of the axis {name} are not consistent between the event and background!") + logger.warning(f" event : {event_edges}") + logger.warning(f" response : {response_edges}") return False - axes_cds = Axes([self.event_dense.axes[0], \ - self.full_detector_response.axes["Em"], \ - self.full_detector_response.axes["Phi"], \ - self.full_detector_response.axes["PsiChi"]]) + if self._coordsys_conv_matrix is None: + axes_cds = Axes([self._image_response.axes["Em"], \ + self._image_response.axes["Phi"], \ + self._image_response.axes["PsiChi"]]) + else: + axes_cds = Axes([self._event.axes[0], \ + self._image_response.axes["Em"], \ + self._image_response.axes["Phi"], \ + self._image_response.axes["PsiChi"]]) - self.event_dense = Histogram(axes_cds, unit = self.event_dense.unit, contents = self.event_dense.contents) + self._event = Histogram(axes_cds, unit = self._event.unit, contents = self._event.contents) - self.bkg_dense = Histogram(axes_cds, unit = self.bkg_dense.unit, contents = self.bkg_dense.contents) + for key in self._bkg_models: + bkg_model = self._bkg_models[key] + self._bkg_models[key] = Histogram(axes_cds, unit = bkg_model.unit, contents = bkg_model.contents) + del bkg_model - print(f"The axes in the event and background files are redefined. Now they are consistent with those of the response file.") + logger.info(f"The axes in the event and background files are redefined. Now they are consistent with those of the response file.") + + return True - def load_full_detector_response_on_memory(self): + def _load_full_detector_response_on_memory(self, full_detector_response, is_miniDC2_format): """ Load a response file on the computer memory. """ - axes_image_response = [self.full_detector_response.axes["NuLambda"], self.full_detector_response.axes["Ei"], - self.full_detector_response.axes["Em"], self.full_detector_response.axes["Phi"], self.full_detector_response.axes["PsiChi"]] + axes_image_response = [full_detector_response.axes["NuLambda"], full_detector_response.axes["Ei"], + full_detector_response.axes["Em"], full_detector_response.axes["Phi"], full_detector_response.axes["PsiChi"]] - self.image_response_dense = Histogram(axes_image_response, unit = self.full_detector_response.unit) + self._image_response = Histogram(axes_image_response, unit = full_detector_response.unit) - nside = self.full_detector_response.axes["NuLambda"].nside - npix = self.full_detector_response.axes["NuLambda"].npix + nside = full_detector_response.axes["NuLambda"].nside + npix = full_detector_response.axes["NuLambda"].npix - if self.is_miniDC2_format: + if is_miniDC2_format: for ipix in tqdm(range(npix)): - self.image_response_dense[ipix] = np.sum(self.full_detector_response[ipix].to_dense(), axis = (4,5)) #Ei, Em, Phi, ChiPsi + self._image_response[ipix] = np.sum(full_detector_response[ipix].to_dense(), axis = (4,5)) #Ei, Em, Phi, ChiPsi else: - contents = self.full_detector_response._file['DRM']['CONTENTS'][:] - self.image_response_dense[:] = contents * self.full_detector_response.unit + contents = full_detector_response._file['DRM']['CONTENTS'][:] + self._image_response[:] = contents * full_detector_response.unit - self.response_on_memory = True - - def calc_exposure_map(self): + def _calc_exposure_map(self): """ Calculate exposure_map, which is an intermidiate matrix used in RL algorithm. """ - print("... (DataLoader) calculating a projected image response ...") - - self.exposure_map = Histogram([ self.coordsys_conv_matrix.axes["lb"], self.full_detector_response.axes["Ei"] ], - unit = self.full_detector_response.unit * self.coordsys_conv_matrix.unit) - - if self.response_on_memory == False: - self.load_full_detector_response_on_memory() - - self.exposure_map[:] = np.tensordot( np.sum(self.coordsys_conv_matrix, axis = (0)), - np.sum(self.image_response_dense, axis = (2,3,4)), - axes = ([1], [0]) ) * self.full_detector_response.unit * self.coordsys_conv_matrix.unit - # [Time/ScAtt, lb, NuLambda] -> [lb, NuLambda] - # [NuLambda, Ei, Em, Phi, PsiChi] -> [NuLambda, Ei] - # [lb, NuLambda] x [NuLambda, Ei] -> [lb, Ei] + logger.info("Calculating an exposure map...") + + if self._coordsys_conv_matrix is None: + self._exposure_map = Histogram(self._model_axes, unit = self._image_response.unit) + self._exposure_map[:] = np.sum(self._image_response.contents, axis = (2,3,4)) + else: + self._exposure_map = Histogram(self._model_axes, unit = self._image_response.unit * self._coordsys_conv_matrix.unit) + self._exposure_map[:] = np.tensordot(np.sum(self._coordsys_conv_matrix, axis = (0)), + np.sum(self._image_response, axis = (2,3,4)), + axes = ([1], [0]) ) * self._image_response.unit * self._coordsys_conv_matrix.unit + # [Time/ScAtt, lb, NuLambda] -> [lb, NuLambda] + # [NuLambda, Ei, Em, Phi, PsiChi] -> [NuLambda, Ei] + # [lb, NuLambda] x [NuLambda, Ei] -> [lb, Ei] - if np.any(self.exposure_map.contents == 0): - print("... There are zero-exposure pixels. Preparing a mask to ignore them ...") - self.mask = Histogram(self.exposure_map.axes, \ - contents = self.exposure_map.contents > 0) + logger.info("Finished...") - def calc_expectation(self, model_map, bkg_norm = 1.0, almost_zero = 1e-12): + def calc_expectation(self, model_map, dict_bkg_norm = None, almost_zero = 1e-12): """ Calculate expected counts from a given model map. @@ -380,6 +268,8 @@ def calc_expectation(self, model_map, bkg_norm = 1.0, almost_zero = 1e-12): ---------- model_map : :py:class:`cosipy.image_deconvolution.ModelMap` Model map + dict_bkg_norm : dict, default None + background normalization for each background model, e.g, {'albedo': 1} almost_zero : float, default 1e-12 In order to avoid zero components in extended count histogram, a tiny offset is introduced. It should be small enough not to effect statistics. @@ -398,36 +288,106 @@ def calc_expectation(self, model_map, bkg_norm = 1.0, almost_zero = 1e-12): # However it is likely that it will have such an axis in the future in order to consider background variability depending on time and pointign direction etc. # Then, the implementation here will not work. Thus, keep in mind that we need to modify it once the response format is fixed. - expectation = Histogram(self.event_dense.axes) - - map_rotated = np.tensordot(self.coordsys_conv_matrix.contents, model_map.contents, axes = ([1], [0])) - # ['Time/ScAtt', 'lb', 'NuLambda'] x ['lb', 'Ei'] -> [Time/ScAtt, NuLambda, Ei] - map_rotated *= self.coordsys_conv_matrix.unit * model_map.unit - map_rotated *= model_map.axes['lb'].pixarea() - # data.coordsys_conv_matrix.contents is sparse, so the unit should be restored. - # the unit of map_rotated is 1/cm2 ( = s * 1/cm2/s/sr * sr) - - expectation[:] = np.tensordot( map_rotated, self.image_response_dense.contents, axes = ([1,2], [0,1])) - expectation += self.bkg_dense * bkg_norm + expectation = Histogram(self.data_axes) + + if self._coordsys_conv_matrix is None: + expectation[:] = np.tensordot( model_map.contents, self._image_response.contents, axes = ([0,1],[0,1])) * model_map.axes['lb'].pixarea() + # ['lb', 'Ei'] x [NuLambda(lb), Ei, Em, Phi, PsiChi] -> [Em, Phi, PsiChi] + else: + map_rotated = np.tensordot(self._coordsys_conv_matrix.contents, model_map.contents, axes = ([1], [0])) + # ['Time/ScAtt', 'lb', 'NuLambda'] x ['lb', 'Ei'] -> [Time/ScAtt, NuLambda, Ei] + map_rotated *= self._coordsys_conv_matrix.unit * model_map.unit + map_rotated *= model_map.axes['lb'].pixarea() + # data.coordsys_conv_matrix.contents is sparse, so the unit should be restored. + # the unit of map_rotated is 1/cm2 ( = s * 1/cm2/s/sr * sr) + expectation[:] = np.tensordot( map_rotated, self._image_response.contents, axes = ([1,2], [0,1])) + # [Time/ScAtt, NuLambda, Ei] x [NuLambda, Ei, Em, Phi, PsiChi] -> [Time/ScAtt, Em, Phi, PsiChi] + + if dict_bkg_norm is not None: + for key in self.keys_bkg_models(): + expectation += self.bkg_model(key) * dict_bkg_norm[key] expectation += almost_zero return expectation - def calc_T_product(self, dataspace_matrix): + def calc_T_product(self, dataspace_histogram): + """ + Calculate the product of the input histogram with the transonse matrix of the response function. + Let R_{ij}, H_{i} be the response matrix and dataspace_histogram, respectively. + Note that i is the index for the data space, and j is for the model space. + In this method, \sum_{j} H{i} R_{ij}, namely, R^{T} H is calculated. + + Parameters + ---------- + dataspace_histogram: :py:class:`histpy.Histogram` + Its axes must be the same as self.data_axes + + Returns + ------- + :py:class:`histpy.Histogram` + The product with self.model_axes + """ + # TODO: currently, dataspace_histogram is assumed to be a dense. hist_unit = self.exposure_map.unit - if dataspace_matrix.unit is not None: - hist_unit *= dataspace_matrix.unit + if dataspace_histogram.unit is not None: + hist_unit *= dataspace_histogram.unit - hist = Histogram(self.exposure_map.axes, \ - unit = hist_unit) + hist = Histogram(self.model_axes, unit = hist_unit) - _ = np.tensordot(dataspace_matrix.contents, self.image_response_dense.contents, axes = ([1,2,3], [2,3,4])) + if self._coordsys_conv_matrix is None: + hist[:] = np.tensordot(dataspace_histogram.contents, self._image_response.contents, axes = ([0,1,2], [2,3,4])) + # [Em, Phi, PsiChi] x [NuLambda (lb), Ei, Em, Phi, PsiChi] -> [NuLambda (lb), Ei] + else: + _ = np.tensordot(dataspace_histogram.contents, self._image_response.contents, axes = ([1,2,3], [2,3,4])) # [Time/ScAtt, Em, Phi, PsiChi] x [NuLambda, Ei, Em, Phi, PsiChi] -> [Time/ScAtt, NuLambda, Ei] - hist[:] = np.tensordot(self.coordsys_conv_matrix.contents, _, axes = ([0,2], [0,1])) \ - * _.unit * self.coordsys_conv_matrix.unit + hist[:] = np.tensordot(self._coordsys_conv_matrix.contents, _, axes = ([0,2], [0,1])) \ + * _.unit * self._coordsys_conv_matrix.unit # [Time/ScAtt, lb, NuLambda] x [Time/ScAtt, NuLambda, Ei] -> [lb, Ei] # note that coordsys_conv_matrix is the sparse, so the unit should be recovered. return hist + + def calc_bkg_model_product(self, key, dataspace_histogram): + """ + Calculate the product of the input histogram with the background model. + Let B_{i}, H_{i} be the background model and dataspace_histogram, respectively. + In this method, \sum_{i} B_{i} H_{i} is calculated. + + Parameters + ---------- + key: str + Background model name + dataspace_histogram: :py:class:`histpy.Histogram` + its axes must be the same as self.data_axes + + Returns + ------- + flaot + """ + # TODO: currently, dataspace_histogram is assumed to be a dense. + + if self._coordsys_conv_matrix is None: + + return np.tensordot(dataspace_histogram.contents, self.bkg_model(key).contents, axes = ([0,1,2], [0,1,2])) + + return np.tensordot(dataspace_histogram.contents, self.bkg_model(key).contents, axes = ([0,1,2,3], [0,1,2,3])) + + def calc_loglikelihood(self, expectation): + """ + Calculate log-likelihood from given expected counts or model/expectation. + + Parameters + ---------- + expectation : :py:class:`histpy.Histogram` + Expected count histogram. + + Returns + ------- + float + Log-likelood + """ + loglikelood = np.sum( self.event * np.log(expectation) ) - np.sum(expectation) + + return loglikelood