Repository: StingraySoftware/stingray

Issue:
Upper limits: add case with norm different from Leahy
See powerspectrum.py:339 

Discussion:
pep8speaks: Hello @matteobachetti! Thanks for opening this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:

* In the file [`stingray/powerspectrum.py`](https://github.com/StingraySoftware/stingray/blob/658917176fc7a7036fccb1002a36e2b78ef18126/stingray/powerspectrum.py):

> [Line 395:101](https://github.com/StingraySoftware/stingray/blob/658917176fc7a7036fccb1002a36e2b78ef18126/stingray/powerspectrum.py#L395): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (101 > 100 characters)





Patch:
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 18a72819..d2474331 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -337,6 +337,10 @@ class Powerspectrum(Crossspectrum):
 
         See `stingray.stats.power_upper_limit`, `stingray.stats.amplitude_upper_limit`
         for more information.
+        
+        The formula used to calculate the upper limit assumes the Leahy normalization.
+        If the periodogram is in another normalization, we will internally convert
+        it to Leahy before calculating the upper limit. 
 
         Parameters
         ----------
@@ -367,20 +371,15 @@ class Powerspectrum(Crossspectrum):
         >>> pds.nphots = 30000
         >>> pds.modulation_upper_limit(fmin=2, fmax=5, c=0.99)
         0.1016...
-        >>> pds.norm = "frac"
-        >>> pds.modulation_upper_limit(2, 5, 0.99)
-        Traceback (most recent call last):
-         ...
-        ValueError: Modulation upper limit is only available in the Leahy...
         """
 
-        # TODO: add case with norm different from Leahy
+        pds = self
         if self.norm != 'leahy':
-            raise ValueError("Modulation upper limit is only available in the Leahy normalization")
+            pds = self.to_norm("leahy")
 
-        freq = self.freq
+        freq = pds.freq
         fnyq = np.max(freq)
-        power = self.power
+        power = pds.power
         freq_mask = freq > 0
         if fmin is not None or fmax is not None:
             if fmin is not None:
@@ -393,10 +392,11 @@ class Powerspectrum(Crossspectrum):
         maximum_val = np.argmax(power)
         nyq_ratio = freq[maximum_val] / fnyq
 
-        # I multiply by M because the formulas from Vaughan+94 treat summed powers, while here
-        # we have averaged powers.
+        # I multiply by M because the formulas from Vaughan+94 treat summed
+        # powers, while here we have averaged powers.
         return amplitude_upper_limit(
-            power[maximum_val] * self.m, self.nphots, n=self.m, c=c, nyq_ratio=nyq_ratio, fft_corr=True)
+            power[maximum_val] * pds.m, pds.nphots, n=pds.m, c=c,
+            nyq_ratio=nyq_ratio, fft_corr=True)
 
     @staticmethod
     def from_time_array(times, dt, segment_size=None, gti=None, norm="frac",
@@ -467,7 +467,7 @@ class Powerspectrum(Crossspectrum):
         silent : bool, default False
             Silence the progress bars
         gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
-            Additional, optional, Good Time intervals, that get interesected with the 
+            Additional, optional, Good Time intervals, that get interesected with the
             GTIs of the input object.
         """
 
@@ -506,7 +506,7 @@ class Powerspectrum(Crossspectrum):
         silent : bool, default False
             Silence the progress bars
         gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
-            Additional, optional, Good Time intervals, that get interesected with the 
+            Additional, optional, Good Time intervals, that get interesected with the
             GTIs of the input object.
         """
 
@@ -1176,7 +1176,7 @@ def powerspectrum_from_events(events, dt, segment_size=None, norm="frac",
     silent : bool, default False
         Silence the progress bars
     gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
-        Additional, optional, Good Time intervals, that get interesected with the 
+        Additional, optional, Good Time intervals, that get interesected with the
         GTIs of the input object.
 
     Returns
@@ -1221,7 +1221,7 @@ def powerspectrum_from_lightcurve(lc, segment_size=None, norm="frac",
     silent : bool, default False
         Silence the progress bars
     gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
-        Additional, optional, Good Time intervals, that get interesected with the 
+        Additional, optional, Good Time intervals, that get interesected with the
         GTIs of the input object.
 
     Returns
@@ -1277,7 +1277,7 @@ def powerspectrum_from_lc_iterable(iter_lc, dt, segment_size=None, norm="frac",
     silent : bool, default False
         Silence the progress bars
     gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
-        Good Time intervals. 
+        Good Time intervals.
 
     Returns
     -------


Test patch:
diff --git a/stingray/tests/test_powerspectrum.py b/stingray/tests/test_powerspectrum.py
index c8c82248..2becfd42 100644
--- a/stingray/tests/test_powerspectrum.py
+++ b/stingray/tests/test_powerspectrum.py
@@ -48,6 +48,16 @@ class TestAveragedPowerspectrumEvents(object):
 
         assert np.isclose(acs_comm.power.std(), acs.power.std(), rtol=0.1)
 
+    @pytest.mark.parametrize("norm", ["frac", "leahy", "none", "abs"])
+    def test_modulation_upper_limit(self, norm):
+        val = 70
+        unnorm_val = 70 * self.leahy_pds.nphots / 2
+        pds = copy.deepcopy(self.leahy_pds)
+        pds.power[25] = val
+        pds.unnorm_power[25] = unnorm_val
+        pds_norm = pds.to_norm(norm)
+        assert np.isclose(pds_norm.modulation_upper_limit(2, 5, 0.99), 0.5412103, atol=1e-4)
+
     def test_legacy_equivalent(self):
         leahy_pds = AveragedPowerspectrum(
             self.lc, segment_size=self.segment_size, dt=self.dt, norm="leahy", silent=True, legacy=True)