diff --git a/cosipy/ts_map/fast_norm_fit.py b/cosipy/ts_map/fast_norm_fit.py index d4df39c0d..ead16df18 100644 --- a/cosipy/ts_map/fast_norm_fit.py +++ b/cosipy/ts_map/fast_norm_fit.py @@ -16,8 +16,8 @@ class FastNormFit: TS(N) = 2 \\sum_i \\left( \\frac{\\log P(d_i; b_i+N e_i)}{\\log P(d_i; b_i)} \\right) - where :math:`P(d; \lambda)` is the Poisson probability of - obtaining :math:`d` count where :math:`\lambda` is expected on + where :math:`P(d; \\lambda)` is the Poisson probability of + obtaining :math:`d` count where :math:`\\lambda` is expected on average; :math:`b` is the estimated number of background counts; :math:`N` is the normalization; and :math:`e` is the expected excess -i.e. signal- per normalization unit -i.e. the number of @@ -37,7 +37,7 @@ class FastNormFit: .. note:: Because of the Poisson probability, :math:`TS(N)` is only - well-defined for :math:`N \geq 1`. By default, + well-defined for :math:`N \\geq 1`. By default, underfluctuations are set to :math:`TS(N=0) = 0`. For cases when there is benefit in letting the normalization float to negative values, you can use `allow_negative`, but in that diff --git a/cosipy/ts_map/moc_ts_fit.py b/cosipy/ts_map/moc_ts_fit.py index 322f698c9..e3afbb509 100644 --- a/cosipy/ts_map/moc_ts_fit.py +++ b/cosipy/ts_map/moc_ts_fit.py @@ -169,7 +169,7 @@ def fit(self, max_nside, spectrum, energy_channel = None, lowest nside used in map strategy : MOCTSMap.Strategy subclass, optional strategy to use in selecting pixels to refine. If None, - default to TopKStrategy with k=8 + default to PaddingStrategy(ContainmentStrategy(0.999)). Returns ------- @@ -197,7 +197,7 @@ def refine(pix): return res if strategy is None: - self.strategy = self.TopKStrategy(k=8) + self.strategy = self.PaddingStrategy(self.ContainmentStrategy(0.999)) else: self.strategy = strategy diff --git a/tests/ts_map/test_fast_ts_map.py b/tests/ts_map/test_fast_ts_map.py index 8db588c3c..c4b41dca7 100644 --- a/tests/ts_map/test_fast_ts_map.py +++ b/tests/ts_map/test_fast_ts_map.py @@ -136,6 +136,48 @@ def test_moc_ts_fit(): ts_results = ts.fit(max_nside = 2, energy_channel = [2,3], spectrum = spectrum, cpu_cores = 1) + ts_values, pixels = ts_results + assert all(pixels == [ + 16, 20, 24, 28, 32, + 36, 40, 44, 48, 52, + 56, 60, 17, 21, 25, + 29, 33, 37, 41, 45, + 49, 53, 57, 61, 18, + 22, 26, 30, 34, 38, + 42, 46, 50, 54, 58, + 62, 19, 23, 27, 31, + 35, 39, 43, 47, 51, + 55, 59, 63 + ]) + + assert np.allclose(ts_values, [ + 40.31750179, 39.40582836, 37.39229509, 39.78630473, 40.39347596, + 40.10805455, 38.79974817, 38.86985166, 40.14551663, 39.92706709, + 39.19653532, 40.07420192, 40.07720833, 38.76345776, 37.46243839, + 40.19657905, 40.41047825, 39.9701239, 37.23577505, 39.28060583, + 40.2803205, 39.83138091, 39.25707762, 39.90376762, 40.20425492, + 39.93311807, 37.35958207, 39.07431591, 40.3217545, 40.1935334, + 38.98305396, 37.24284051, 40.41652632, 39.61022258, 39.22641583, + 40.09740223, 39.81314166, 39.49561083, 37.43747447, 39.65452285, + 39.940159, 39.61014066, 37.24223294, 37.35636565, 40.65108076, + 39.92533792, 37.24385822, 40.2989865 + ]) + + ts.plot_ts(*ts_results, + skycoord = SkyCoord(l=0, b=0, unit=u.deg, frame="galactic")) + + ts.plot_ts(*ts_results, containment = 0.9, save_plot = True, + save_dir = "", save_name = "ts_map.png") + + assert Path("ts_map.png").exists() + + os.remove("ts_map.png") + + # test top-k strategy + ts_results = ts.fit(max_nside = 2, energy_channel = [2,3], + spectrum = spectrum, cpu_cores = 1, + strategy=MOCTSMap.TopKStrategy(k=8)) + ts_values, pixels = ts_results assert all(pixels == [ 6, 10, 11, 14, 16, @@ -159,22 +201,15 @@ def test_moc_ts_fit(): 40.2989865 ]) - ts.plot_ts(*ts_results, - skycoord = SkyCoord(l=0, b=0, unit=u.deg, frame="galactic")) - - ts.plot_ts(*ts_results, containment = 0.9, save_plot = True, - save_dir = "", save_name = "ts_map.png") - - assert Path("ts_map.png").exists() - - os.remove("ts_map.png") - # test containment strategy ts_results = ts.fit(max_nside = 2, energy_channel = [2,3], spectrum = spectrum, strategy=MOCTSMap.ContainmentStrategy(0.9)) ts_values, pixels = ts_results + print(pixels) + print(ts_values) + assert all(pixels == [ 16, 20, 24, 28, 32, 36, 40, 44, 48, 52,