Merge branch 'ctapipe_0.17' of github.com:cta-observatory/cta-lstchai…

…n into ctapipe_0.17_new_write_pedestal_ff_events_to_h5

lstchain/high_level/interpolate.py

@@ -13,13 +13,13 @@
     create_psf_table_hdu,
 )
 from pyirf.interpolation import (
+    GridDataInterpolator,
     interpolate_effective_area_per_energy_and_fov,
     interpolate_energy_dispersion,
     interpolate_psf_table,
-    # interpolate_rad_max,
+    interpolate_rad_max,
 )
 from scipy.spatial import Delaunay, distance, QhullError
-from scipy.interpolate import griddata
 
 log = logging.getLogger(__name__)
 
@@ -318,67 +318,38 @@ def interpolate_gh_cuts(
     method="linear",
 ):
     """
-    Interpolates a grid of GH CUTS tables to a target-point.
+    Interpolates a grid of GH_CUTS tables to a target-point. Same as pyirf's
+    interpolate_rad_max function.
+
     Wrapper around scipy.interpolate.griddata [1].
 
     Parameters
     ----------
     gh_cuts: numpy.ndarray, shape=(N, M, ...)
-        Gammaness-cuts for all combinations of grid-points, energy and
-        fov_offset.
-        Shape (N:n_grid_points, M:n_energy_bins, n_fov_offset_bins)
+        Gammaness-cuts for all combinations of grid-points, like energy.
+        Shape (N:n_grid_points, M:n_energy_bins)
     grid_points: numpy.ndarray, shape=(N, O)
         Array of the N O-dimensional morphing parameter values corresponding
         to the N input templates.
     target_point: numpy.ndarray, shape=(O)
         Value for which the interpolation is performed (target point)
-    method: 'linear’, ‘nearest’, ‘cubic’
+    method: 'linear', 'nearest', 'cubic'
         Interpolation method for scipy.interpolate.griddata [1].
         Defaults to 'linear'.
 
     Returns
     -------
     gh_cuts_interp: numpy.ndarray, shape=(1, M, ...)
-        Gammaness-cuts for the target grid-point,
-        shape (1, M:n_energy_bins, n_fov_offset_bins)
-    """
-
-    return griddata(grid_points, gh_cuts, target_point, method=method)
-
-
-def interpolate_rad_max(
-    rad_max,
-    grid_points,
-    target_point,
-    method="linear",
-):
-    """
-    Interpolates a grid of RAD_MAX tables to a target-point.
-    Wrapper around scipy.interpolate.griddata [1].
+        Gammaness-cuts for the target grid-point, shape (1, M:n_energy_bins)
 
-    Parameters
+    References
     ----------
-    rad_max: numpy.ndarray, shape=(N, M, ...)
-        Theta-cuts for all combinations of grid-points, energy and
-        fov_offset.
-        Shape (N:n_grid_points, M:n_energy_bins, n_fov_offset_bins)
-    grid_points: numpy.ndarray, shape=(N, O)
-        Array of the N O-dimensional morphing parameter values corresponding
-        to the N input templates.
-    target_point: numpy.ndarray, shape=(O)
-        Value for which the interpolation is performed (target point)
-    method: 'linear’, ‘nearest’, ‘cubic’
-        Interpolation method for scipy.interpolate.griddata [1].
-        Defaults to 'linear'.
-
-    Returns
-    -------
-    rad_max_interp: numpy.ndarray, shape=(1, M, ...)
-        Theta-cuts for the target grid-point,
-        shape (1, M:n_energy_bins, n_fov_offset_bins)
+    .. [1] https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html
     """
+    interp = GridDataInterpolator(grid_points=grid_points, params=gh_cuts)
+    gh_cuts_interp = interp(target_point, method=method)
 
-    return griddata(grid_points, rad_max, target_point, method=method)
+    return gh_cuts_interp
 
 
 def interpolate_irf(irfs, data_pars, interp_method="linear"):
@@ -480,11 +451,14 @@ def interpolate_irf(irfs, data_pars, interp_method="linear"):
         fov_off = np.append(temp_irf["THETA_LO"][0], temp_irf["THETA_HI"][0][-1])
 
         aeff_interp = interpolate_effective_area_per_energy_and_fov(
-            effarea_list, irf_pars_sel, interp_pars_sel, method=interp_method
+            effective_area=effarea_list,
+            grid_points=irf_pars_sel,
+            target_point=interp_pars_sel,
+            method=interp_method
         )
 
         aeff_hdu_interp = create_aeff2d_hdu(
-            aeff_interp.T,
+            effective_area=aeff_interp.T[0],
             true_energy_bins=e_true,
             fov_offset_bins=fov_off,
             point_like=point_like,
@@ -510,11 +484,15 @@ def interpolate_irf(irfs, data_pars, interp_method="linear"):
         fov_off = np.append(temp_irf["THETA_LO"][0], temp_irf["THETA_HI"][0][-1])
 
         edisp_interp = interpolate_energy_dispersion(
-            e_migra, edisp_list, irf_pars_sel, interp_pars_sel, quantile_resolution=1e-3
+            migra_bins=e_migra,
+            edisps=edisp_list,
+            grid_points=irf_pars_sel,
+            target_point=interp_pars_sel,
+            quantile_resolution=1e-3
         )
 
         edisp_hdu_interp = create_energy_dispersion_hdu(
-            edisp_interp,
+            energy_dispersion=edisp_interp[0],
             true_energy_bins=e_true,
             migration_bins=e_migra,
             fov_offset_bins=fov_off,
@@ -537,7 +515,10 @@ def interpolate_irf(irfs, data_pars, interp_method="linear"):
         temp_irf = QTable.read(irfs[0], hdu="GH_CUTS")
 
         gh_cut_interp = interpolate_gh_cuts(
-            gh_cuts_list, irf_pars_sel, interp_pars_sel, method=interp_method
+            gh_cuts=gh_cuts_list,
+            grid_points=irf_pars_sel,
+            target_point=interp_pars_sel,
+            method=interp_method
         )
 
         gh_header = fits.Header()
@@ -562,10 +543,13 @@ def interpolate_irf(irfs, data_pars, interp_method="linear"):
         temp_irf = QTable.read(irfs[0], hdu="RAD_MAX")
 
         rad_max_interp = interpolate_rad_max(
-            radmax_list, irf_pars_sel, interp_pars_sel, method=interp_method
+            rad_max=radmax_list,
+            grid_points=irf_pars_sel,
+            target_point=interp_pars_sel,
+            method=interp_method
         )
 
-        temp_irf["RAD_MAX"] = rad_max_interp.T[np.newaxis, ...] * u.deg
+        temp_irf["RAD_MAX"] = rad_max_interp[0].T[np.newaxis, ...] * u.deg
 
         radmax_header = fits.Header()
         radmax_header["CREATOR"] = f"lstchain v{__version__}"
@@ -593,14 +577,14 @@ def interpolate_irf(irfs, data_pars, interp_method="linear"):
             fov_off = np.append(temp_irf["THETA_LO"][0], temp_irf["THETA_HI"][0][-1])
 
             psf_interp = interpolate_psf_table(
-                src_bins,
-                psf_list,
-                irf_pars_sel,
-                interp_pars_sel,
+                source_offset_bins=src_bins,
+                psfs=psf_list,
+                grid_points=irf_pars_sel,
+                target_point=interp_pars_sel,
                 quantile_resolution=1e-3,
             )
             psf_hdu_interp = create_psf_table_hdu(
-                psf_interp,
+                psf=psf_interp[0],
                 true_energy=e_true,
                 source_offset_bins=src_bins,
                 fov_offset_bins=fov_off,

lstchain/high_level/tests/test_interpolate.py

@@ -23,12 +23,8 @@ def test_interp_irf(simulated_irf_file, simulated_dl2_file):
     irf_file_g_3 = simulated_irf_file.parent / "irf_interp_1.fits.gz"
     irf_file_g_final = simulated_irf_file.parent / "irf_interp_final.fits.gz"
 
-    hdus_2_g = [
-        fits.PrimaryHDU(),
-    ]
-    hdus_3_g = [
-        fits.PrimaryHDU(),
-    ]
+    hdus_2_g = [fits.PrimaryHDU(), ]
+    hdus_3_g = [fits.PrimaryHDU(), ]
 
     # En-dep, point-like cuts IRF
     irf_file_en_1 = simulated_dl2_file.parent / "en_dep_cut_irf_1.fits.gz"
@@ -49,12 +45,8 @@ def test_interp_irf(simulated_irf_file, simulated_dl2_file):
         "--theta-containment=0.8",
     )
 
-    hdus_2_en = [
-        fits.PrimaryHDU(),
-    ]
-    hdus_3_en = [
-        fits.PrimaryHDU(),
-    ]
+    hdus_2_en = [fits.PrimaryHDU(), ]
+    hdus_3_en = [fits.PrimaryHDU(), ]
 
     for irf in [simulated_irf_file, irf_file_en_1]:
         # Change the effective area for different angular pointings
@@ -67,15 +59,15 @@ def test_interp_irf(simulated_irf_file, simulated_dl2_file):
         az_1 = u.Quantity(aeff_1_meta["AZ_PNT"], "deg").to_value(u.rad)
         del_1 = u.Quantity(aeff_1_meta["B_DELTA"], "deg").to_value(u.rad)
 
-        zen_2 = 2 * zen_1
-        az_2 = 2 * az_1
-        del_2 = 1.2 * del_1
+        factor_czd = 0.7
+        factor_sd = 1.2
 
-        factor_zd = np.cos(zen_2) / np.cos(zen_1)
-        factor_del = np.sin(del_2) / np.sin(del_1)
+        zen_2 = np.arccos(np.cos(zen_1) * factor_czd)
+        az_2 = az_1 * factor_czd
+        del_2 = np.arcsin(np.sin(del_1) * factor_sd)
 
-        aeff_1["EFFAREA"][0] *= factor_zd
-        aeff_2["EFFAREA"][0] *= factor_zd * factor_del
+        aeff_1["EFFAREA"][0] *= factor_czd
+        aeff_2["EFFAREA"][0] *= factor_czd * factor_sd
 
         aeff_1_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
         aeff_1_meta["AZ_PNT"] = (az_1 * 180 / np.pi, "deg")
@@ -93,8 +85,8 @@ def test_interp_irf(simulated_irf_file, simulated_dl2_file):
         edisp_2 = edisp_1.copy()
         edisp_2_meta = edisp_1_meta.copy()
 
-        edisp_1["MATRIX"][0] *= factor_zd
-        edisp_2["MATRIX"][0] *= factor_zd * factor_del
+        edisp_1["MATRIX"][0] *= factor_czd
+        edisp_2["MATRIX"][0] *= factor_czd * factor_sd
 
         edisp_1_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
         edisp_1_meta["AZ_PNT"] = (az_1 * 180 / np.pi, "deg")
@@ -119,8 +111,8 @@ def test_interp_irf(simulated_irf_file, simulated_dl2_file):
 
             mask = gh_1["cut"] < gh_1["cut"].max()
 
-            gh_1["cut"][mask] *= factor_zd
-            gh_2["cut"][mask] *= factor_zd * factor_del
+            gh_1["cut"][mask] *= factor_czd
+            gh_2["cut"][mask] *= factor_czd * factor_sd
 
             gh_1_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
             gh_1_meta["AZ_PNT"] = (az_1 * 180 / np.pi, "deg")
@@ -140,8 +132,8 @@ def test_interp_irf(simulated_irf_file, simulated_dl2_file):
 
             mask = th_1["RAD_MAX"] < th_1["RAD_MAX"].max()
 
-            th_1["RAD_MAX"][mask] *= factor_zd
-            th_2["RAD_MAX"][mask] *= factor_zd * factor_del
+            th_1["RAD_MAX"][mask] *= factor_czd
+            th_2["RAD_MAX"][mask] *= factor_czd * factor_sd
 
             th_1_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
             th_1_meta["AZ_PNT"] = (az_1 * 180 / np.pi, "deg")
@@ -176,10 +168,17 @@ def test_interp_irf(simulated_irf_file, simulated_dl2_file):
 
     irfs_g = [simulated_irf_file, irf_file_g_2, irf_file_g_3]
     irfs_en = [irf_file_en_1, irf_file_en_2, irf_file_en_3]
+
+    factor_czd_t = 1 - (1 - factor_czd) / 2  # as factor_czd < 1
+    factor_sd_t = 1 + (factor_sd - 1) / 2
+
+    zen_t = np.arccos(np.cos(zen_1) * factor_czd_t) * 180 / np.pi
+    del_t = np.arcsin(np.sin(del_1) * factor_sd_t) * 180 / np.pi
+    az_t = az_1 * factor_czd_t * 180 / np.pi
     data_pars = {
-        "ZEN_PNT": 30 * u.deg,
-        "B_DELTA": (del_1 * 0.8 * u.rad).to(u.deg),
-        "AZ_PNT": 120 * u.deg,
+        "ZEN_PNT": zen_t * u.deg,
+        "B_DELTA": del_t * u.deg,
+        "AZ_PNT": az_t * u.deg,
     }
 
     hdu_g = interpolate_irf(irfs_g, data_pars)
@@ -188,12 +187,12 @@ def test_interp_irf(simulated_irf_file, simulated_dl2_file):
     hdu_en = interpolate_irf(irfs_en, data_pars)
     hdu_en.writeto(irf_file_en_final, overwrite=True)
 
-    assert hdu_g[1].header["ZEN_PNT"] == 30
+    assert hdu_g[1].header["ZEN_PNT"] == zen_t
     assert irf_file_g_2.exists()
     assert irf_file_g_3.exists()
     assert irf_file_g_final.exists()
 
-    assert hdu_en[1].header["ZEN_PNT"] == 30
+    assert hdu_en[1].header["ZEN_PNT"] == zen_t
     assert irf_file_en_2.exists()
     assert irf_file_en_3.exists()
     assert irf_file_en_final.exists()
@@ -280,21 +279,50 @@ def test_check_delaunay_triangles(simulated_irf_file):
 
     irfs = [simulated_irf_file, irf_file_3, irf_file_4, irf_file_5]
 
+    # Fetch the grid parameters to get target points inside and outside grids.
+    czd_list = []
+    sd_list = []
+    az_list = []
+
+    for i in irfs:
+        meta_ = fits.open(i)["EFFECTIVE AREA"].header
+        czd_list.append(np.cos(meta_["ZEN_PNT"] * np.pi/180))
+        sd_list.append(np.sin(meta_["B_DELTA"] * np.pi/180))
+        az_list.append(meta_["AZ_PNT"])
+    czd_list = np.array(czd_list)
+    sd_list = np.array(sd_list)
+    az_list = np.array(az_list)
+
     # Check on target being inside or outside Delaunay simplex
-    data_pars = {"ZEN_PNT": 25 * u.deg, "B_DELTA": 45 * u.deg, "AZ_PNT": 100 * u.deg}
-    data_pars2 = {"ZEN_PNT": 58 * u.deg, "B_DELTA": 70 * u.deg, "AZ_PNT": 200 * u.deg}
+    zen_t1 = np.arccos(np.mean(czd_list)) * 180 / np.pi
+    zen_t2 = np.arccos(np.min(czd_list) - 0.1) * 180 / np.pi
+    del_t = np.arcsin(np.mean(sd_list)) * 180 / np.pi
+    az_close_t = az_list[np.where(czd_list == np.min(czd_list))[0][0]]
 
-    new_irfs = check_in_delaunay_triangle(irfs, data_pars)
-    new_irfs2 = check_in_delaunay_triangle(irfs, data_pars2)
-    new_irfs3 = check_in_delaunay_triangle(irfs, data_pars, use_nearest_irf_node=True)
-    new_irfs4 = check_in_delaunay_triangle([irfs[0]], data_pars)
+    data_pars = {
+        "ZEN_PNT": zen_t1 * u.deg,
+        "B_DELTA": del_t * u.deg,
+        "AZ_PNT": 100 * u.deg
+    }
+    data_pars2 = {
+        "ZEN_PNT": zen_t2 * u.deg,
+        "B_DELTA": del_t * u.deg,
+        "AZ_PNT": az_close_t * u.deg
+    }
+
+    new_irfs_in_grid = check_in_delaunay_triangle(irfs, data_pars)
+    new_irfs_out_grid = check_in_delaunay_triangle(irfs, data_pars2)
+    new_irfs_in_grid_near_az = check_in_delaunay_triangle(
+        irfs, data_pars2, use_nearest_irf_node=True
+    )
+    new_irfs_no_grid = check_in_delaunay_triangle([irfs[0]], data_pars)
 
-    t3 = Table.read(new_irfs3[0], hdu=1).meta
+    az_check = Table.read(new_irfs_in_grid_near_az[0], hdu=1).meta
 
-    assert len(new_irfs) == 3
-    assert len(new_irfs2) == 1
-    assert t3["ZEN_PNT"] == 20
-    assert len(new_irfs4) == 1
+    assert len(new_irfs_in_grid) == 3
+    assert len(new_irfs_out_grid) == 1
+    assert az_check["AZ_PNT"] == az_close_t
+    assert len(new_irfs_no_grid) == 1
 
 
 def test_get_nearest_az_node():
@@ -358,7 +386,6 @@ def test_get_nearest_az_node():
 def test_interpolate_gh_cuts():
     from lstchain.high_level.interpolate import interpolate_gh_cuts
 
-    # Similar function as interpolate_rad_max, hence no need for extra test
     # linear test case
     gh_cuts_1 = np.array([[0, 0], [0.1, 0], [0.2, 0.1], [0.3, 0.2]])
     gh_cuts_2 = 2 * gh_cuts_1