Skip to content

Commit

Permalink
fix(background_models.py): implementing a smooth method
Browse files Browse the repository at this point in the history 
implementing a smooth method for smooth function. This allows to override the default smoothing method. Adding Gaussian2DKernel and convolve_fft to imports

#67
  • Loading branch information
@steob92
steob92 committed Dec 17, 2024
1 parent 2d46ac4 commit f17a55b
Show file tree
Hide file tree
Showing 6 changed files with 101 additions and 63 deletions.
2 changes: 1 addition & 1 deletion Dockerfile
View file
Original file line number Diff line number Diff line change
@@ -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


Expand Down
16 changes: 8 additions & 8 deletions gammapy-tools.def
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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 .
Expand All @@ -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
Expand Down
97 changes: 60 additions & 37 deletions gammapy_tools/make_background/background_models.py
View file
Original file line number Diff line number Diff line change
@@ -1,19 +1,19 @@
from copy import deepcopy

import numpy as np
from scipy.ndimage import gaussian_filter
from multiprocess import Pool

# Astropy stuff


# 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

Expand All @@ -22,6 +22,7 @@

from typing import Optional


class BackgroundModelEstimator:
"""
Largely based off https://docs.gammapy.org/0.18.2/tutorials/background_model.html
Expand Down Expand Up @@ -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
----------
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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)
Expand All @@ -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()
Expand All @@ -390,18 +404,20 @@ 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)
Parameters
----------
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
Expand All @@ -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
44 changes: 29 additions & 15 deletions gammapy_tools/make_background/make_background.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -132,48 +132,62 @@ 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,
smooth=config["background_selection"]["smooth"],
smooth_sigma=smooth_sigma,
njobs=config["config"]["njobs"],
default_exclusion=default_exclusion,
)
)
else:
estimator = BackgroundModelEstimator(
energy,
offset,
smooth=config["background_selection"]["smooth"],
smooth_sigma=smooth_sigma,
njobs=config["config"]["njobs"],
)
)

estimator.run(observations)
# Check if a background currently exists
Expand Down Expand Up @@ -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,
)
Expand Down
1 change: 1 addition & 0 deletions gammapy_tools/make_background/prepare_data.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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?"
Expand Down
4 changes: 2 additions & 2 deletions pyproject.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down

0 comments on commit f17a55b

Please sign in to comment.