tests fixed

kalman_detector/efficiency.py

@@ -134,15 +134,9 @@ def eff_plot(
     niters = 10000
     nchans = 336
     corr_len = 300
-    freqs = np.arange(0, 335) + 1104
+    freqs = np.arange(nchans) + 1104
     template = utils.simulate_gaussian_signal(nchans, corr_len, complex_process=True)
 
-    figure = plt.figure(figsize=(5, 5.5), dpi=200)
-    grid = figure.add_gridspec(left=0.05, right=0.95, bottom=0.1, top=0.95)
-    inner_grid = grid[0, 0].subgridspec(
-        nrows=2, ncols=1, hspace=0.24, height_ratios=(2, 1)
-    )
-    ax_prof = plt.subplot(inner_grid[1, 0])
-    ax_eff = plt.subplot(inner_grid[0, 0])
+    fig, (ax_eff, ax_prof) = plt.subplots(2, 1, height_ratios=(2, 1), figsize=(5, 5.5))
     results = eff_plot(template, freqs, ax_eff, ax_prof, niters=niters)
     plt.show()

kalman_detector/main.py

@@ -117,7 +117,7 @@ def get_best_significance(self, spec: np.ndarray) -> float:
         Returns
         -------
         float
-            The best kalman significance
+            The best kalman significance in logsf units.
         """
         sigs, scores = self.get_significance(spec)
         return np.min(sigs)
@@ -346,7 +346,7 @@ def secondary_spectrum_cumulative_chi2_score(
     Parameters
     ----------
     spec : numpy.ndarray
-        1D array of the observed spectrum I(f) of the candidate FRB.
+        1d array of the observed spectrum.
     spec_std : numpy.ndarray
         1D array of the standard deviation of the observed spectrum.
     mask_tol : float, optional
@@ -355,7 +355,7 @@ def secondary_spectrum_cumulative_chi2_score(
     Returns
     -------
     float
-        Maximum value of the significance. max_f0 of ln(P(sig|H1,ff0)/P(sig|H0))
+        Best significance in logsf units (max_f0 of ln(P(sig|H1,ff0)/P(sig|H0)))
 
     Notes
     -----
@@ -375,4 +375,4 @@ def secondary_spectrum_cumulative_chi2_score(
     score_arr = np.cumsum(fft_sig[1:] / (len(signal) / 2))
     dof_arr = 2 * np.arange(1, len(score_arr) + 1)
     significance_arr = stats.chi2.logsf(score_arr, dof_arr)
-    return np.max(-significance_arr)
+    return np.min(significance_arr)

kalman_detector/utils.py

@@ -133,9 +133,12 @@ def simulate_gaussian_signal(
     rng = np.random.default_rng()
     kernel = np.exp(-np.linspace(-nchans / corr_len, nchans / corr_len, nchans) ** 2)
     kernel /= np.dot(kernel, kernel) ** 0.5
-    base_process = rng.normal(0, 1 / np.sqrt(corr_len), nchans)
     if complex_process:
-        base_process += 1.0j * rng.normal(0, 1 / np.sqrt(corr_len), nchans)
+        base_process = rng.normal(0, 1 / np.sqrt(corr_len), nchans) + 1.0j * rng.normal(
+            0, 1 / np.sqrt(corr_len), nchans
+        )
+    else:
+        base_process = rng.normal(0, 1 / np.sqrt(corr_len), nchans)
     signal = np.abs(np.fft.ifft(np.fft.fft(kernel) * np.fft.fft(base_process)))
     signal -= np.mean(signal)
     return signal / np.dot(signal, signal) ** 0.5

setup.cfg

@@ -97,14 +97,19 @@ omit =
 	*__init__.py
 	*tests*
 	*docs*
+    kalman_detector/efficiency.py
 
 [coverage:report]
 show_missing = True
 ignore_errors = True
-#fail_under = 85
-exclude_lines =
-	raise AssertionError
-	raise NotImplementedError
+exclude_also =
+    def __repr__
+    raise AssertionError
+    raise NotImplementedError
+    if __name__ == .__main__.:
+    if TYPE_CHECKING:
+    class .*\bProtocol\):
+    @(abc\.)?abstractmethod
     @jit
 
 [coverage:paths]

tests/test_kalman.py

@@ -1,6 +1,5 @@
 import pytest
 import numpy as np
-from scipy import stats
 from pathlib import Path
 from numpy.polynomial import Polynomial
 
@@ -9,45 +8,35 @@
     KalmanDistribution,
     secondary_spectrum_cumulative_chi2_score,
 )
+from kalman_detector import utils
 
 
 class TestKalmanDetector(object):
-    def test_sig_ts_float(self):
-        sig_ts = 0.1
-        std_vec = np.arange(0.1, 1, 0.01)
-        kalman = KalmanDetector(std_vec, sig_ts)
-        np.testing.assert_array_equal(kalman.transit_sigmas, np.array([sig_ts]))
-
-    def test_sig_ts_list(self):
-        sig_ts = [0.1, 0.2]
-        std_vec = np.arange(0.1, 1, 0.01)
-        kalman = KalmanDetector(std_vec, sig_ts)
-        np.testing.assert_array_equal(kalman.transit_sigmas, np.array(sig_ts))
-
-    def test_sig_ts_fail(self):
-        sig_ts = {"a": 0.1, "b": 0.2}
-        std_vec = np.arange(0.1, 1, 0.01)
-        with pytest.raises(ValueError):
-            kalman = KalmanDetector(std_vec, sig_ts)
-
-    def test_sig_ts_nan_fail(self):
-        sig_ts = [0.1, np.nan]
-        std_vec = np.arange(0.1, 1, 0.01)
-        with pytest.raises(ValueError):
-            kalman = KalmanDetector(std_vec, sig_ts)
+    def test_q_par_float(self):
+        q_par = 0.1
+        std_vec = np.ones(336)
+        kalman = KalmanDetector(std_vec, q_par**2)
+        np.testing.assert_array_equal(kalman.transit_sigmas, np.array([q_par]))
+
+    def test_q_par_list(self):
+        q_par = [0.1, 0.2]
+        q_par_sq = [qq**2 for qq in q_par]
+        std_vec = np.ones(336)
+        kalman = KalmanDetector(std_vec, q_par_sq)
+        np.testing.assert_array_equal(kalman.transit_sigmas, np.array(q_par))
 
     def test_initialization_fail(self):
         std_vec = np.zeros(1024)
         with pytest.raises(ValueError):
-            kalman = KalmanDetector(std_vec)
+            KalmanDetector(std_vec)
 
     def test_prepare_fits(self):
-        sig_ts = [0.1, 0.2]
+        q_par = [0.1, 0.2]
         std_vec = np.arange(0.1, 1, 0.01)
-        kalman = KalmanDetector(std_vec, sig_ts)
+        kalman = KalmanDetector(std_vec, q_par)
         kalman.prepare_fits(ntrials=1000)
         assert isinstance(kalman.distributions[0], KalmanDistribution)
-        np.testing.assert_equal(len(kalman.distributions), len(sig_ts))
+        np.testing.assert_equal(len(kalman.distributions), len(q_par))
 
     def test_get_significance(self):
         nchans = 128
@@ -57,8 +46,9 @@ def test_get_significance(self):
         kalman = KalmanDetector(std_vec)
         kalman.prepare_fits(ntrials=1000)
         spectrum = rng.normal(target, std_vec)
-        sig_kalman = kalman.get_significance(spectrum)
-        assert sig_kalman <= 0
+        sigs, scores = kalman.get_significance(spectrum)
+        np.testing.assert_equal(len(sigs), len(scores))
+        np.testing.assert_array_less(sigs, 0.1)
 
     def test_get_significance_fail(self):
         nchans = 128
@@ -68,134 +58,80 @@ def test_get_significance_fail(self):
         kalman = KalmanDetector(std_vec)
         kalman.prepare_fits(ntrials=1000)
         with pytest.raises(ValueError):
-            sig_kalman = kalman.get_significance(spectrum)
+            kalman.get_significance(spectrum)
 
-    def test_sensitivity_up_down(self):
-        nchans = 800
-        signal = np.concatenate([np.ones(10), np.zeros(10)] * 40)
-        std_noise = np.sin(0.1 * np.arange(nchans)) + 2
-        optimal_snr = np.sum(signal**2 / std_noise**2) ** 0.5
-        snr_sum = np.sum(signal / std_noise**2) / np.sum(1 / std_noise**2) ** 0.5
+    def test_get_best_significance(self):
+        nchans = 336
+        corr_len = 300
+        ntrials = 1000
+        target_snr = 20
 
+        scores_arr = []
+        template = utils.simulate_gaussian_signal(nchans, corr_len, complex_process=True)
+        spec_mean = np.zeros_like(template)
+        spec_std = np.ones_like(template)
+        kalman = KalmanDetector(spec_std)
+        kalman.prepare_fits(ntrials=1000)
         rng = np.random.default_rng()
-        kalman = KalmanDetector(std_noise)
-        kalman.prepare_fits(ntrials=10000)
-        sig_kalman = []
-        sig_naive = []
-        for i in range(1000):
-            spectrum = signal + rng.normal(0, std_noise)
-            snr_naive = (
-                np.sum(spectrum / std_noise**2) / np.sum(1 / std_noise**2) ** 0.5
-            )
-            sig_naive.append(stats.norm.logsf(snr_naive))
-            sig_kalman.append(kalman.get_significance(spectrum))
-
-        sig_naive_mean = np.mean(sig_naive)
-        improvement = (np.array(sig_naive) + np.array(sig_kalman)) / stats.norm.logsf(
-            optimal_snr
-        )
-        np.testing.assert_almost_equal(
-            sig_naive_mean, stats.norm.logsf(snr_sum), decimal=0
-        )
-        np.testing.assert_almost_equal(np.mean(improvement), 0.6, decimal=1)
-
-    def test_sin(self):
-        nchans = 400
-        rng = np.random.default_rng()
-        sig_list = [
-            rng.random()
-            * np.sin(
-                rng.random() * 2 * np.pi + rng.random(nchans) * np.arange(nchans) / 4
-            )
-            for i in range(4)
-        ]
-        signal = np.abs((1.5 + np.sum(sig_list, axis=0)) * 0.6)
-        std_noise = np.sin(0.1 * np.arange(nchans)) + 2
-        optimal_snr = np.sum(signal**2 / std_noise**2) ** 0.5
-        snr_sum = np.sum(signal / std_noise**2) / np.sum(1 / std_noise**2) ** 0.5
-
-        kalman = KalmanDetector(std_noise)
-        kalman.prepare_fits(ntrials=10000)
-        sig_kalman = []
-        sig_naive = []
-        optimal_sig_list = []
-
-        for i in range(1000):
-            spectrum = signal + rng.normal(0, std_noise)
-            snr_naive = (
-                np.sum(spectrum / std_noise**2) / np.sum(1 / std_noise**2) ** 0.5
-            )
-            sig_naive.append(stats.norm.logsf(snr_naive))
-            optimal_sig_list.append(
-                stats.norm.logsf(
-                    np.sum(signal * spectrum / std_noise**2)
-                    / np.sum(signal**2 / std_noise**2) ** 0.5
-                )
-            )
-            sig_kalman.append(kalman.get_significance(spectrum))
-
-        sig_naive_mean = np.mean(sig_naive)
-        improvement = (np.array(sig_naive) + np.array(sig_kalman)) / stats.norm.logsf(
-            optimal_snr
-        )
-        np.testing.assert_almost_equal(
-            sig_naive_mean, stats.norm.logsf(snr_sum), decimal=0
-        )
-        np.testing.assert_almost_equal(np.mean(improvement), 0.8, decimal=1)
+        for _ in range(ntrials):
+            spec = target_snr * template + rng.normal(spec_mean, spec_std, len(template))
+            best_sig = kalman.get_best_significance(spec)
+            scores_arr.append(best_sig)
+        score_mean = utils.get_snr_from_logsf(np.mean(scores_arr))
+        np.testing.assert_allclose(score_mean, target_snr, rtol=0.2)
 
 
 class TestKalmanDistribution(object):
     def test_initialization(self):
-        std_vec = np.arange(0.1, 1, 0.01)
-        t_sig = 0.1
+        sigma_arr = np.arange(0.1, 1, 0.01)
+        sig_eta = 0.1
         mask_tol = 0.1
         ntrials = 1000
-        kdist = KalmanDistribution(std_vec, t_sig, ntrials=ntrials, mask_tol=mask_tol)
+        kdist = KalmanDistribution(sigma_arr, sig_eta, ntrials=ntrials, mask_tol=mask_tol)
         np.testing.assert_equal(kdist.mask_tol, mask_tol)
-        np.testing.assert_equal(len(kdist.mask), len(std_vec))
-        np.testing.assert_equal(kdist.t_sig, t_sig)
+        np.testing.assert_equal(len(kdist.mask), len(sigma_arr))
+        np.testing.assert_equal(kdist.sig_eta, sig_eta)
         np.testing.assert_equal(kdist.ntrials, ntrials)
 
     def test_initialization_fail(self):
-        std_vec = np.zeros(1024)
+        sigma_arr = np.zeros(1024)
         with pytest.raises(ValueError):
-            kdist = KalmanDistribution(std_vec, 0.1)
+            KalmanDistribution(sigma_arr, 0.1)
 
     def test_plot_diagnostic(self, tmpfile):
         outfile_path = Path(f"{tmpfile}.pdf")
-        std_vec = np.arange(0.1, 1, 0.01)
-        kdist = KalmanDistribution(std_vec, 0.01, ntrials=1000)
-        fig = kdist.plot_diagnostic(logy=True, outfile=outfile_path)
+        sigma_arr = np.arange(0.1, 1, 0.01)
+        kdist = KalmanDistribution(sigma_arr, 0.01, ntrials=1000)
+        kdist.plot_diagnostic(logy=True, outfile=outfile_path)
         assert outfile_path.is_file()
         outfile_path.unlink()
 
     def test_polyfit(self):
-        std_vec = np.arange(0.1, 1, 0.01)
-        kdist = KalmanDistribution(std_vec, 0.01, ntrials=1000)
+        sigma_arr = np.arange(0.1, 1, 0.01)
+        kdist = KalmanDistribution(sigma_arr, 0.01, ntrials=1000)
         assert isinstance(kdist.polyfit, Polynomial)
 
     def test_str(self):
-        std_vec = np.arange(0.1, 1, 0.01)
-        kdist = KalmanDistribution(std_vec, 0.01, ntrials=1000)
+        sigma_arr = np.arange(0.1, 1, 0.01)
+        kdist = KalmanDistribution(sigma_arr, 0.01, ntrials=1000)
         assert str(kdist).startswith("KalmanDistribution")
 
-    def test_repr(self):
-        std_vec = np.arange(0.1, 1, 0.01)
-        kdist = KalmanDistribution(std_vec, 0.01, ntrials=1000)
-        assert repr(kdist).startswith("KalmanDistribution")
-
 
 class TestSecondarySpectrym(object):
     def test_secondary_spectrum_cumulative_chi2_score(self):
-        nchans = 128
+        nchans = 336
+        corr_len = 300
         ntrials = 1000
+        target_snr = 20
+
         scores_arr = []
-        target = 5
+        template = utils.simulate_gaussian_signal(nchans, corr_len, complex_process=True)
+        spec_mean = np.zeros_like(template)
+        spec_std = np.ones_like(template)
         rng = np.random.default_rng()
         for _ in range(ntrials):
-            std_vec = rng.normal(1, 0.1, size=nchans)
-            spectrum = rng.normal(target, std_vec)
-            score = secondary_spectrum_cumulative_chi2_score(spectrum, std_vec)
+            spec = target_snr * template + rng.normal(spec_mean, spec_std, len(template))
+            score = secondary_spectrum_cumulative_chi2_score(spec, spec_std)
             scores_arr.append(score)
-        new_score_mean = np.mean(scores_arr)
-        np.testing.assert_almost_equal(new_score_mean, target, decimal=0)
+        score_mean = utils.get_snr_from_logsf(np.mean(scores_arr))
+        np.testing.assert_allclose(score_mean, target_snr, rtol=0.2)