From f17a55bdc77ee42caf2f8e06524fa56c8f3331e6 Mon Sep 17 00:00:00 2001 From: Ste O'Brien Date: Tue, 17 Dec 2024 12:23:35 -0500 Subject: [PATCH] fix(background_models.py): implementing a smooth method implementing a smooth method for smooth function. This allows to override the default smoothing method. Adding Gaussian2DKernel and convolve_fft to imports #67 --- Dockerfile | 2 +- gammapy-tools.def | 16 +-- .../make_background/background_models.py | 97 ++++++++++++------- .../make_background/make_background.py | 44 ++++++--- gammapy_tools/make_background/prepare_data.py | 1 + pyproject.toml | 4 +- 6 files changed, 101 insertions(+), 63 deletions(-) diff --git a/Dockerfile b/Dockerfile index 4a209b7..4c791fb 100644 --- a/Dockerfile +++ b/Dockerfile @@ -1,7 +1,7 @@ FROM jupyter/minimal-notebook AS base # Install gammapy -RUN mamba install gcc jupyterlab "gammapy==1.2" ipykernel --yes +RUN mamba install gcc jupyterlab "gammapy==1.3" ipykernel --yes WORKDIR /gammapy-tools diff --git a/gammapy-tools.def b/gammapy-tools.def index eddc558..e68bacb 100644 --- a/gammapy-tools.def +++ b/gammapy-tools.def @@ -10,13 +10,13 @@ Stage: build # Temporay build files ../gammapy-tools /gamma-tools/tmp_build/gammapy-tools - ./gammapy_tools/Hipparcos_MAG8_1997.dat /gamma-tools/gammapy-datasets/1.1/catalogs/Hipparcos_MAG8_1997.dat + ./gammapy_tools/Hipparcos_MAG8_1997.dat /gamma-tools/gammapy-datasets/1.3/catalogs/Hipparcos_MAG8_1997.dat %environment # Port for jupyter lab export JPORT=8000 - export GAMMAPY_DATA=/gamma-tools/gammapy-datasets/1.2/ - export GAMMAPY_DATA=/gamma-tools/gammapy-datasets/dev + export GAMMAPY_DATA=/gamma-tools/gammapy-datasets/1.3/ + #export GAMMAPY_DATA=/gamma-tools/gammapy-datasets/dev #. "/opt/conda/etc/profile.d/conda.sh" #. /opt/conda/bin/activate #export PATH=/opt/conda/bin/:$PATH @@ -36,18 +36,18 @@ Stage: build python3 -m venv /opt/gammapy_tools . /opt/gammapy_tools/bin/activate - pip install iminuit cmasher pip papermill matplotlib uproot "jupyterlab==4.0.12" notebook ipykernel ipython ipywidgets + pip install iminuit cmasher pip papermill matplotlib uproot "jupyterlab==4.0.12" notebook ipykernel ipython ipywidgets astropy #mamba install -c conda-forge iminuit cmasher pip papermill matplotlib pip "jupyterlab==4.0.12" notebook ipykernel ipython ipywidgets --yes # Install v2dl3 cd /gamma-tools - git clone https://github.com/VERITAS-Observatory/V2DL3.git + git clone --depth 1 --branch v0.6.0 https://github.com/VERITAS-Observatory/V2DL3.git cd /gamma-tools/V2DL3 # Because its an env install of requirements # Note the python version in the latest test throws issues with pytest - grep -A 100 "dependencies:" environment-eventdisplay.yml | grep "-" | grep -v "python" | awk '{print $2}' | xargs pip install + grep -A 100 "dependencies:" environment-eventdisplay.yml | grep "-" | grep -v "python" | grep -v "astropy-base" | awk '{print $2}' | xargs pip install # root_numpy throws issues too. Only VEGAS uses it mv setup.py _setup.py && grep -v "root_numpy" _setup.py > setup.py && pip install . @@ -69,8 +69,8 @@ Stage: build python -m pip cache purge - # Add the updated version of gammapy/datasets/map.py - wget https://raw.githubusercontent.com/gammapy/gammapy/main/gammapy/datasets/map.py -O /opt/gammapy_tools/lib/python3.12/site-packages/gammapy/datasets/map.py + # # Add the updated version of gammapy/datasets/map.py + # wget https://raw.githubusercontent.com/gammapy/gammapy/main/gammapy/datasets/map.py -O /opt/gammapy_tools/lib/python3.12/site-packages/gammapy/datasets/map.py %runscript cd /local_data ; jupyter-lab --port=$JPORT diff --git a/gammapy_tools/make_background/background_models.py b/gammapy_tools/make_background/background_models.py index 095a9c0..643313c 100644 --- a/gammapy_tools/make_background/background_models.py +++ b/gammapy_tools/make_background/background_models.py @@ -1,7 +1,6 @@ from copy import deepcopy import numpy as np -from scipy.ndimage import gaussian_filter from multiprocess import Pool # Astropy stuff @@ -9,11 +8,12 @@ # Gammapy stuff import gammapy -from gammapy.irf import Background2D,Background3D +from gammapy.irf import Background2D, Background3D from gammapy.maps import MapAxis from gammapy.data import Observation -from gammapy.datasets import MapDatasetEventSampler, Datasets, MapDatasetOnOff +from gammapy.datasets import MapDatasetEventSampler from scipy.interpolate import splrep, splev +from astropy.convolution import Gaussian2DKernel, convolve_fft # from gammapy.catalog import SourceCatalog3HWC, SourceCatalogGammaCat @@ -22,6 +22,7 @@ from typing import Optional + class BackgroundModelEstimator: """ Largely based off https://docs.gammapy.org/0.18.2/tutorials/background_model.html @@ -50,8 +51,10 @@ def __init__( excluded_sources (list) - list of sources to be excluded (astropy.coordinates.Skycoords) njobs (int) - Number of jobs to run - default_exclusion (float) - Default exclusion radius for sources. Defaults to None - which will use the exclusion finder to determine the exclusion radius + default_exclusion (float) - Default exclusion radius for sources. + Defaults to None + which will use the exclusion finder to determine + the exclusion radius Returns ---------- @@ -92,7 +95,7 @@ def __init__( # ) # self.hawc = SourceCatalog3HWC("$GAMMAPY_DATA/catalogs/3HWC.ecsv") - self.exclusion = ExclusionFinder(default_exclusion = default_exclusion) + self.exclusion = ExclusionFinder(default_exclusion=default_exclusion) @staticmethod def _make_bkg2d(energy: MapAxis, offset: MapAxis, unit: str) -> Background2D: @@ -343,22 +346,28 @@ def background_rate( return rate + class Background3DModelEstimator: - - def __init__(self, - energy, - fov_lon, + + def __init__( + self, + energy, + fov_lon, fov_lat, smooth: bool = False, excluded_sources: list = [], - smooth_sigma: float = 1,): + smooth_sigma: float = 1, + smooth_method: str = "gauss", + ): self.excluded_sources = excluded_sources self.smooth = smooth self.smooth_sigma = smooth_sigma + self.smooth_method = smooth_method self.counts = self._make_bkg3d(energy, fov_lon, fov_lat, unit="") self.exposure = self._make_bkg3d(energy, fov_lon, fov_lat, unit="s TeV sr") self.exclusion = ExclusionFinder() + @staticmethod def _make_bkg3d(energy, fov_lon, fov_lat, unit): return Background3D(axes=[energy, fov_lon, fov_lat], unit=unit) @@ -376,10 +385,15 @@ def fill_counts(self, obs): events = MapDatasetEventSampler.event_det_coords(obs, events) run_mask = self.exclusion.exclude_events(events.table, run_ra, run_dec) counts = np.histogramdd( - (events.energy.to('TeV')[run_mask], events.table['DETX'][run_mask], events.table['DETY'][run_mask]), - (energy.edges, fov_lon.edges, fov_lat.edges))[0] + ( + events.energy.to("TeV")[run_mask], + events.table["DETX"][run_mask], + events.table["DETY"][run_mask], + ), + (energy.edges, fov_lon.edges, fov_lat.edges), + )[0] self.counts.data += counts - + def fill_exposure(self, obs): time = obs.observation_time_duration bin_volume = self.exposure.axes.bin_volume() @@ -390,10 +404,11 @@ def background_rate(self): rate = deepcopy(self.counts) rate.quantity /= self.exposure.quantity if self.smooth: - rate = smooth(rate, sigma=self.smooth_sigma) + rate = smooth(rate, sigma=self.smooth_sigma, method=self.smooth_method) return rate -def smooth(bkg, method='poly',sigma=1): + +def smooth(bkg, method="poly", sigma=1): """Smooths background rates from BackgroundModelEstimator.background_rate (bkg input) @@ -401,7 +416,8 @@ def smooth(bkg, method='poly',sigma=1): ---------- bkg (Background2D) - 2D background to be smoothed method (str) - Method for smoothing: - "poly": 5th deg polynomial fit (same as ED) [default] + "poly": 5th deg polynomial fit + (same as ED) [default] "spline": 10 d.o.f. cubic spline fit "gauss": gaussian kernel smoothing @@ -410,37 +426,44 @@ def smooth(bkg, method='poly',sigma=1): bkg (Background2D) - Smoothed 2D background """ - - if method == 'poly': - offset = bkg.axes['offset'].center.value + + if method == "poly": + offset = bkg.axes["offset"].center.value for i in range(len(bkg.data)): - poly = np.polyfit(offset,bkg.data[i,:],deg=5) + poly = np.polyfit(offset, bkg.data[i, :], deg=5) p = np.poly1d(poly) - bkg.data[i,:] = p(offset) - elif method == 'spline': - offset = bkg.axes['offset'].center.value - smooth_offset = np.linspace(offset[0],offset[-1],10) + bkg.data[i, :] = p(offset) + elif method == "spline": + offset = bkg.axes["offset"].center.value + smooth_offset = np.linspace(offset[0], offset[-1], 10) for i in range(len(bkg.data)): - spline = splrep(offset,bkg.data[i,:],deg=5) - fit = splev(smooth_offset,spline) - bkg.data[i,:] = fit(offset) - elif method == 'gauss': - if type(bkg) == gammapy.irf.background.Background2D: + spline = splrep(offset, bkg.data[i, :], deg=5) + fit = splev(smooth_offset, spline) + bkg.data[i, :] = fit(offset) + elif method == "gauss": + if isinstance(bkg, gammapy.irf.background.Background2D): bkg = bkg.to_3d() - else: bkg_3d = bkg + # SOB... Unused? + # else: + # bkg_3d = bkg for i in range(len(bkg.data)): data = bkg.data[i, :, :] - padded_data = np.pad(data,10) + padded_data = np.pad(data, 10) window1d = np.abs(np.blackman(len(padded_data))) - window2d = np.sqrt(np.outer(window1d,window1d)) + window2d = np.sqrt(np.outer(window1d, window1d)) windowed = padded_data * window2d - gauss = Gaussian2DKernel(sigma,sigma,x_size=np.shape(windowed)[0],y_size=np.shape(windowed)[1]) + gauss = Gaussian2DKernel( + sigma, sigma, x_size=np.shape(windowed)[0], y_size=np.shape(windowed)[1] + ) - conv = convolve_fft(data,gauss,normalize_kernel=False,fft_pad=True,preserve_nan=True) + conv = convolve_fft( + data, gauss, normalize_kernel=False, fft_pad=True, preserve_nan=True + ) bkg.data[i, :, :] = conv - if type(bkg) == gammapy.irf.background.Background2D: + if isinstance(bkg, gammapy.irf.background.Background2D): return bkg.to_2d() - else: return bkg + else: + return bkg return bkg diff --git a/gammapy_tools/make_background/make_background.py b/gammapy_tools/make_background/make_background.py index 976db3e..0edc78b 100644 --- a/gammapy_tools/make_background/make_background.py +++ b/gammapy_tools/make_background/make_background.py @@ -18,7 +18,7 @@ from pyV2DL3.generateObsHduIndex import create_obs_hdu_index_file # From here -from .background_models import BackgroundModelEstimator,Background3DModelEstimator +from .background_models import BackgroundModelEstimator, Background3DModelEstimator from .background_tools import process_run, get_requested_exposure from ..utils.run_details import find_data_mimic @@ -132,32 +132,46 @@ def get_background_for_run(parms: tuple[float, dict]) -> tuple[str, list]: if "smooth_sigma" in config["background_selection"]: smooth_sigma = config["background_selection"]["smooth_sigma"] - default_exclusion = None + default_exclusion = None if "sky_map" in config: if "on_exclusion_region" in config["sky_map"]: default_exclusion = config["sky_map"]["on_exclusion_region"] - - if make_3d: fov_lon = MapAxis.from_bounds( - -config["binning"]["off_max"], config["binning"]["off_max"], nbin=config["binning"]["off_bins"], interp="lin", unit="deg",name="fov_lon" + -config["binning"]["off_max"], + config["binning"]["off_max"], + nbin=config["binning"]["off_bins"], + interp="lin", + unit="deg", + name="fov_lon", ) - fov_lat =MapAxis.from_bounds( - -config["binning"]["off_max"], config["binning"]["off_max"], nbin=config["binning"]["off_bins"], interp="lin", unit="deg",name="fov_lat" + fov_lat = MapAxis.from_bounds( + -config["binning"]["off_max"], + config["binning"]["off_max"], + nbin=config["binning"]["off_bins"], + interp="lin", + unit="deg", + name="fov_lat", ) - if default_exclusion != None: + if default_exclusion is not None: + smooth_method = ( + "gauss" + if "smooth_method" not in config["background_selection"] + else config["background_selection"]["smooth_method"] + ) estimator = Background3DModelEstimator( energy, fov_lon, fov_lat, smooth=config["background_selection"]["smooth"], smooth_sigma=smooth_sigma, + smooth_method=smooth_method, ) else: - raise Exception("3D Exclusions not implemented") + raise Exception("3D Exclusions not implemented") else: - if default_exclusion != None: + if default_exclusion is not None: estimator = BackgroundModelEstimator( energy, offset, @@ -165,7 +179,7 @@ def get_background_for_run(parms: tuple[float, dict]) -> tuple[str, list]: smooth_sigma=smooth_sigma, njobs=config["config"]["njobs"], default_exclusion=default_exclusion, - ) + ) else: estimator = BackgroundModelEstimator( energy, @@ -173,7 +187,7 @@ def get_background_for_run(parms: tuple[float, dict]) -> tuple[str, list]: smooth=config["background_selection"]["smooth"], smooth_sigma=smooth_sigma, njobs=config["config"]["njobs"], - ) + ) estimator.run(observations) # Check if a background currently exists @@ -249,14 +263,14 @@ def generate_background_from_run(parms: tuple[int, dict]) -> str: name="offset", ) - default_exclusion = None + default_exclusion = None if "sky_map" in config: if "on_exclusion_region" in config["sky_map"]: default_exclusion = config["sky_map"]["on_exclusion_region"] - estimator = BackgroundModelEstimator( - energy, offset, + energy, + offset, smooth=config["background_selection"]["smooth"], default_exclusion=default_exclusion, ) diff --git a/gammapy_tools/make_background/prepare_data.py b/gammapy_tools/make_background/prepare_data.py index 4e62f0c..25045ca 100644 --- a/gammapy_tools/make_background/prepare_data.py +++ b/gammapy_tools/make_background/prepare_data.py @@ -54,6 +54,7 @@ def prepare_dataset(config: dict, overwrite: bool = False) -> dict: try: data_store.copy_obs(obs_in_db, in_dir, overwrite=overwrite) except Exception as e: + print(e) if len(obs_in_db) == 0: raise RuntimeError( f"Observations cannot be found in {db_dir}, do these files exist?" diff --git a/pyproject.toml b/pyproject.toml index 6f40159..71166ec 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,8 +26,8 @@ requires-python = ">=3.9" dependencies = [ - # "gammapy>=1.2", - "gammapy@git+https://github.com/gammapy/gammapy", + "gammapy>=1.3", + # "gammapy@git+https://github.com/gammapy/gammapy", "pydantic", "scipy==1.11.4", "astropy",