Repository: StingraySoftware/HENDRICS

Issue:
Use new API for covariance spectrum in HENvarenergy

Discussion:

Patch:
diff --git a/hendrics/fold.py b/hendrics/fold.py
index e585006..26dfbc9 100644
--- a/hendrics/fold.py
+++ b/hendrics/fold.py
@@ -276,7 +276,6 @@ def get_TOAs_from_events(
         # BAD!BAD!BAD!
         # [[Pay attention to time reference here.
         # We are folding wrt pepoch, and calculating TOAs wrt start]]
-
         toa, toaerr = get_TOA(
             profile,
             1 / frequency_derivatives[0],
diff --git a/hendrics/fspec.py b/hendrics/fspec.py
index fb33ade..9cd183e 100644
--- a/hendrics/fspec.py
+++ b/hendrics/fspec.py
@@ -145,7 +145,7 @@ def calc_pds(
         return
 
     ftype, data = get_file_type(lcfile)
-    mjdref = data.instr
+    mjdref = data.mjdref
     instr = data.instr
 
     lc_data = _format_lc_data(data, ftype, bintime=bintime, fftlen=fftlen)
@@ -551,7 +551,6 @@ def main(args=None):
 
     args = check_negative_numbers_in_args(args)
     args = parser.parse_args(args)
-    log.info("Starting")
 
     if args.debug:
         args.loglevel = "DEBUG"
diff --git a/hendrics/io.py b/hendrics/io.py
index ad29548..280b5bb 100644
--- a/hendrics/io.py
+++ b/hendrics/io.py
@@ -326,6 +326,8 @@ def read_from_netcdf(fname):
             to_save = values[0]
         else:
             to_save = values
+        if isinstance(to_save, (str, bytes)) and to_save.startswith("__bool_"):
+            to_save = eval(to_save.replace("__bool__", ""))
         out[k] = to_save
 
     rootgrp.close()
@@ -482,7 +484,7 @@ def save_lcurve(lcurve, fname, lctype="Lightcurve"):
     elif hasattr(lcurve, "e_interval") and lcurve.e_interval is not None:
         out["e_interval"] = lcurve.e_interval
     if hasattr(lcurve, "use_pi"):
-        out["use_pi"] = int(lcurve.use_pi)
+        out["use_pi"] = lcurve.use_pi
 
     if hasattr(lcurve, "instr") and lcurve.instr is not None:
         out["instr"] = lcurve.instr.lower()
@@ -542,7 +544,7 @@ def load_lcurve(fname):
     if "e_interval" in list(data.keys()):
         lcurve.e_interval = data["e_interval"]
     if "use_pi" in list(data.keys()):
-        lcurve.use_pi = bool(data["use_pi"])
+        lcurve.use_pi = data["use_pi"]
     if "header" in list(data.keys()):
         lcurve.header = data["header"]
     if "base" in list(data.keys()):
@@ -607,42 +609,32 @@ def save_pds(cpds, fname, save_all=False):
     if not hasattr(cpds, "instr"):
         outdata["instr"] = "unknown"
 
-    if hasattr(cpds, "show_progress"):
-        outdata["show_progress"] = "T" if cpds.show_progress else "F"
-
-    if hasattr(cpds, "amplitude"):
-        outdata["amplitude"] = int(cpds.amplitude)
+    for attr in ['show_progress', 'amplitude']:
+        if hasattr(cpds, attr):
+            outdata[attr] = getattr(cpds, attr)
 
     outdir = fname.replace(HEN_FILE_EXTENSION, "")
     if save_all:
         mkdir_p(outdir)
 
-    if "lc1" in outdata:
-        if save_all and isinstance(cpds.lc1, Lightcurve):
-            save_lcurve(
-                cpds.lc1, os.path.join(outdir, "__lc1__" + HEN_FILE_EXTENSION)
-            )
-        outdata.pop("lc1")
-    if "lc2" in outdata:
-        if save_all and isinstance(cpds.lc2, Lightcurve):
-            save_lcurve(
-                cpds.lc2, os.path.join(outdir, "__lc2__" + HEN_FILE_EXTENSION)
-            )
-        outdata.pop("lc2")
-    if "pds1" in outdata:
-        if save_all:
-            save_pds(
-                cpds.pds1,
-                os.path.join(outdir, "__pds1__" + HEN_FILE_EXTENSION),
-            )
-        outdata.pop("pds1")
-    if "pds2" in outdata:
-        if save_all:
-            save_pds(
-                cpds.pds2,
-                os.path.join(outdir, "__pds2__" + HEN_FILE_EXTENSION),
-            )
-        outdata.pop("pds2")
+    for attr in ['lc1', 'lc2', 'pds1', 'pds2']:
+        if save_all and hasattr(cpds, attr):
+            value = getattr(cpds, attr)
+
+            outf = f"__{attr}__" + HEN_FILE_EXTENSION
+            if "lc" in attr and isinstance(value, Lightcurve):
+                save_lcurve(
+                    value,
+                    os.path.join(outdir, outf)
+                )
+            elif "pds" in attr and isinstance(value, Crossspectrum):
+                save_pds(
+                    value,
+                    os.path.join(outdir, outf),
+                    save_all=False
+                )
+        outdata.pop(attr, None)
+
     if "cs_all" in outdata:
         if save_all:
             for i, c in enumerate(cpds.cs_all):
@@ -694,12 +686,6 @@ def load_pds(fname, nosub=False):
     for key in data.keys():
         setattr(cpds, key, data[key])
 
-    if "amplitude" in list(data.keys()):
-        cpds.amplitude = bool(data["amplitude"])
-
-    if "show_progress" in list(data.keys()):
-        cpds.show_progress = data["show_progress"] == "T"
-
     outdir = fname.replace(HEN_FILE_EXTENSION, "")
     modelfiles = glob.glob(
         os.path.join(outdir, fname.replace(HEN_FILE_EXTENSION, "__mod*__.p"))
@@ -847,8 +833,10 @@ def _save_data_nc(struct, fname, kind="data"):
 
     for k in struct.keys():
         var = struct[k]
-
+        if isinstance(var, bool):
+            var = f"__bool__{var}"
         probe = var
+
         if isinstance(var, Iterable) and len(var) >= 1:
             probe = var[0]
 
diff --git a/hendrics/phaseogram.py b/hendrics/phaseogram.py
index ce587ae..cead714 100644
--- a/hendrics/phaseogram.py
+++ b/hendrics/phaseogram.py
@@ -647,13 +647,24 @@ class InteractivePhaseogram(BasePhaseogram):
         dfreq, dfdot, dfddot = self._read_sliders()
         freqs = [self.freq - dfreq, self.fdot - dfdot, self.fddot - dfddot]
         folding_length = np.median(np.diff(self.times))
-        template = gaussian_filter1d(
-            np.sum(np.nan_to_num(self.unnorm_phaseogr), axis=1), sigma=1
-        )
-        template = (
-            np.roll(template[: template.size // 2], -np.argmax(template))
-            / self.nt
-        )
+        nbin = self.nph
+        template_raw = np.sum(np.nan_to_num(self.unnorm_phaseogr), axis=1)
+        template = gaussian_filter1d(template_raw, sigma=1)
+
+        if self.nph < 64:
+            warnings.warn("TOA calculation is not robust if the "
+                          "number of bins is < 64. Oversampling.")
+            nbin = 64
+            phases = np.arange(self.nph * 2 + 1) / (self.nph * 2)
+            fun = interp1d(phases,
+                           np.concatenate((template, [template[0]])),
+                           kind='cubic')
+            new_phases = np.arange(0, nbin * 2) / (nbin * 2)
+            template = fun(new_phases)
+
+        template = template[nbin // 2: nbin // 2 + template.size // 2]
+
+        template = (np.roll(template, -np.argmax(template)) / self.nt)
 
         toa, toaerr = get_TOAs_from_events(
             self.ev_times,
@@ -871,17 +882,24 @@ def run_interactive_phaseogram(
     events = load_events(event_file)
     if emin is not None or emax is not None:
         events, elabel = filter_energy(events, emin, emax)
-    try:
+
+    position = name = None
+
+    if hasattr(events, 'header') and events.header is not None:
         header = Header.fromstring(events.header)
-        position = SkyCoord(
-            header["RA_OBJ"],
-            header["DEC_OBJ"],
-            unit="deg",
-            frame=header["RADECSYS"].lower(),
-        )
-        name = header["OBJECT"]
-    except (KeyError, AttributeError):
-        position = name = None
+
+        try:
+            position = SkyCoord(
+                header["RA_OBJ"],
+                header["DEC_OBJ"],
+                unit="deg",
+                frame=header["RADECSYS"].lower(),
+            )
+        except (KeyError, AttributeError):
+            position = None
+
+        if 'OBJECT' in header:
+            name = header["OBJECT"]
 
     pepoch_mjd = pepoch
     if pepoch is None:
diff --git a/hendrics/varenergy.py b/hendrics/varenergy.py
index 8fca90e..259d67b 100644
--- a/hendrics/varenergy.py
+++ b/hendrics/varenergy.py
@@ -7,7 +7,10 @@ import warnings
 from astropy import log
 from astropy.logger import AstropyUserWarning
 import numpy as np
-from stingray.varenergyspectrum import RmsEnergySpectrum
+from stingray.varenergyspectrum import (
+    RmsEnergySpectrum,
+    _decode_energy_specification,
+)
 from stingray.varenergyspectrum import LagEnergySpectrum
 
 # from stingray.covariancespectrum import AveragedCovariancespectrum
@@ -179,33 +182,26 @@ def main(args=None):
                 filelist.append(out2)
 
             if args.covariance:
-                try:
-                    from stingray.varenergyspectrum import (
-                        CovarianceEnergySpectrum,
-                    )
-                except Exception:
-                    warnings.warn(
-                        "This version of Stingray does not implement "
-                        "the correct version of Covariance Spectrum.",
-                        AstropyUserWarning,
-                    )
-                    continue
-                cov = CovarianceEnergySpectrum(
+                from stingray.covariancespectrum import (
+                    AveragedCovariancespectrum,
+                )
+
+                energies = _decode_energy_specification(energy_spec)
+                energies = list(zip(energies[:-1], energies[1:]))
+                cov = AveragedCovariancespectrum(
                     events,
-                    args.freq_interval,
-                    energy_spec,
-                    args.ref_band,
+                    dt=args.bintime,
+                    band_interest=energies,
+                    ref_band_interest=args.ref_band,
                     segment_size=args.segment_size,
-                    bin_time=args.bintime,
-                    events2=events2,
-                    use_pi=args.use_pi,
                 )
-                start_energy = np.asarray(cov.energy_intervals)[:, 0]
-                stop_energy = np.asarray(cov.energy_intervals)[:, 1]
+
+                start_energy = np.asarray(energies)[:, 0]
+                stop_energy = np.asarray(energies)[:, 1]
                 out2 = hen_root(fname) + "_cov" + ".qdp"
                 save_as_qdp(
-                    [start_energy, stop_energy, cov.spectrum],
-                    [None, None, cov.spectrum_error],
+                    [start_energy, stop_energy, cov.covar],
+                    [None, None, cov.covar_error],
                     filename=out2,
                 )
                 filelist.append(out2)
diff --git a/setup.cfg b/setup.cfg
index a1d8b18..31862bc 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -112,6 +112,25 @@ astropy_header = true
 doctest_plus = enabled
 text_file_format = rst
 addopts = --doctest-rst
+filterwarnings =
+    ignore:SIMON says. Errorbars on cross spectra are not thoroughly tested.:UserWarning
+    ignore:SIMON says. Stingray only uses poisson err_dist at the moment.:UserWarning
+    ignore:SIMON says. Looks like your lightcurve:UserWarning
+    ignore:Matplotlib is currently using agg, which is a:UserWarning
+    ignore:Using or importing the ABCs from 'collections':DeprecationWarning
+    ignore:unclosed file:ResourceWarning
+    ignore:numpy.ufunc size changed:RuntimeWarning
+    ignore:numpy.ndarray size changed:RuntimeWarning
+    ignore:invalid value encountered:RuntimeWarning
+    ignore:divide by zero encountered:RuntimeWarning
+    ignore:underflow encountered:RuntimeWarning
+    ignore:overflow encountered:RuntimeWarning
+    ignore:table path was not set via the path= argument; using default path:UserWarning
+    ignore:Beware! For cpds and derivatives, I assume that:UserWarning
+    ignore:SIMON says. Lightcurves have different statistics:UserWarning
+    ignore: SIMON says. Coherence is not ideal:UserWarning
+    ignore:SIMON says. The excess variance in the:UserWarning
+    ignore:SIMON says. Number of segments used in averaging:UserWarning
 
 [coverage:run]
 omit =


Test patch:
diff --git a/hendrics/tests/test_a_complete_run.py b/hendrics/tests/test_a_complete_run.py
index b469ae3..c09f7bf 100644
--- a/hendrics/tests/test_a_complete_run.py
+++ b/hendrics/tests/test_a_complete_run.py
@@ -49,6 +49,10 @@ log.setLevel("DEBUG")
 # log.basicConfig(filename='HEN.log', level=log.DEBUG, filemode='w')
 
 
+def find_file_pattern_in_dir(pattern, directory):
+    return glob.glob(os.path.join(directory, pattern))
+
+
 class TestFullRun(object):
     """Test how command lines work.
 
@@ -185,6 +189,29 @@ class TestFullRun(object):
         out = hen.base.hen_root(fname) + "_rms" + ".qdp"
         os.path.exists(out)
 
+    def test_save_varen_rms(self):
+        fname = self.ev_fileAcal
+        hen.varenergy.main(
+            [
+                fname,
+                "-f",
+                "0",
+                "100",
+                "--energy-values",
+                "0.3",
+                "12",
+                "5",
+                "lin",
+                "--covariance",
+                "-b",
+                "0.5",
+                "--segment-size",
+                "128",
+            ]
+        )
+        out = hen.base.hen_root(fname) + "_cov" + ".qdp"
+        os.path.exists(out)
+
     def test_save_varen_lag(self):
         fname = self.ev_fileAcal
         hen.varenergy.main(
@@ -213,8 +240,9 @@ class TestFullRun(object):
         command = ("{0} -b 100 -e {1} {2} {2} {3}").format(
             self.ev_fileA, 3, 5, 10
         )
-        with pytest.raises(ValueError) as excinfo:
-            hen.colors.main(command.split())
+        with catch_warnings():
+            with pytest.raises(ValueError) as excinfo:
+                hen.colors.main(command.split())
 
         assert "No energy information is present " in str(excinfo.value)
 
@@ -307,13 +335,17 @@ class TestFullRun(object):
             ]
         )
 
+    def test_all_files_get_read(self):
+        # Test that HENreadfile works with all file types
+        fnames = find_file_pattern_in_dir(
+            "*" + HEN_FILE_EXTENSION, self.datadir
+        )
+        io.main(fnames)
+
     @classmethod
     def teardown_class(self):
         """Test a full run of the scripts (command lines)."""
 
-        def find_file_pattern_in_dir(pattern, directory):
-            return glob.glob(os.path.join(directory, pattern))
-
         patterns = [
             "*monol_test*" + HEN_FILE_EXTENSION,
             "*lcurve*" + HEN_FILE_EXTENSION,
diff --git a/hendrics/tests/test_fspec.py b/hendrics/tests/test_fspec.py
index 6cf3f87..4c5ec5a 100644
--- a/hendrics/tests/test_fspec.py
+++ b/hendrics/tests/test_fspec.py
@@ -633,7 +633,7 @@ model = models.Const1D()
             "*.qdp",
             "*.inf",
             "*_cpds" + HEN_FILE_EXTENSION,
-            "*_ev" + HEN_FILE_EXTENSION
+            "*_ev" + HEN_FILE_EXTENSION,
         ]
 
         file_list = []
diff --git a/hendrics/tests/test_io.py b/hendrics/tests/test_io.py
index 8822130..496f1cb 100644
--- a/hendrics/tests/test_io.py
+++ b/hendrics/tests/test_io.py
@@ -1,5 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 
+from astropy.io.fits import Header
 from stingray.events import EventList
 from stingray.lightcurve import Lightcurve
 from stingray.powerspectrum import Powerspectrum, AveragedPowerspectrum
@@ -139,6 +140,7 @@ class TestIO:
         )
         events.energy = np.array([3.0, 4.0, 5.0])
         events.mission = "nustar"
+        events.header = Header().tostring()
         save_events(events, self.dum)
         events2 = load_events(self.dum)
         assert np.allclose(events.time, events2.time)
@@ -146,6 +148,7 @@ class TestIO:
         assert np.allclose(events.mjdref, events2.mjdref)
         assert np.allclose(events.gti, events2.gti)
         assert np.allclose(events.energy, events2.energy)
+        assert events.header == events2.header
         assert events2.mission == events.mission
 
     def test_load_and_save_lcurve(self):
@@ -169,12 +172,13 @@ class TestIO:
         pds.power = np.random.poisson(30, 15)
         pds.mjdref = 54385.3254923845
         pds.gti = np.longdouble([[-0.5, 3.5]])
+        pds.show_progress = True
+        pds.amplitude = False
 
         save_pds(pds, self.dum)
         pds2 = load_pds(self.dum)
-        assert np.allclose(pds.gti, pds2.gti)
-        assert np.allclose(pds.mjdref, pds2.mjdref)
-        assert pds.m == pds2.m
+        for attr in ['gti', 'mjdref', 'm', 'show_progress', 'amplitude']:
+            assert np.allclose(getattr(pds, attr), getattr(pds2, attr))
 
     def test_load_and_save_xps(self):
         lcurve1 = Lightcurve(
diff --git a/hendrics/tests/test_phaseogram.py b/hendrics/tests/test_phaseogram.py
index 3335730..8565be4 100644
--- a/hendrics/tests/test_phaseogram.py
+++ b/hendrics/tests/test_phaseogram.py
@@ -1,3 +1,4 @@
+from astropy.io.fits import Header
 from stingray.lightcurve import Lightcurve
 from stingray.events import EventList
 import numpy as np
@@ -32,7 +33,21 @@ class TestPhaseogram:
         events.mjdref = 57000.0
         cls.event_times = events.time
         cls.dum = "events" + HEN_FILE_EXTENSION
+        cls.dum_nohead = "events_nohead" + HEN_FILE_EXTENSION
+        cls.dum_info = "events_info" + HEN_FILE_EXTENSION
+        save_events(events, cls.dum_nohead)
+
+        header = Header()
+        events.header = header.tostring()
+
         save_events(events, cls.dum)
+        header = Header()
+        header["OBJECT"] = "BUBU"
+        header["RADECSYS"] = "FK5"
+        header["RA_OBJ"] = 0.4
+        header["DEC_OBJ"] = 30.4
+        events.header = header.tostring()
+        save_events(events, cls.dum_info)
 
         curdir = os.path.abspath(os.path.dirname(__file__))
         cls.datadir = os.path.join(curdir, "data")
@@ -70,8 +85,9 @@ class TestPhaseogram:
         # Defaults to 2 harmonics
         assert efperiod.N == 2
 
-    def test_phaseogram_input_periodogram(self):
-        evfile = self.dum
+    @pytest.mark.parametrize('label', ('', '_info', '_nohead'))
+    def test_phaseogram_input_periodogram(self, label):
+        evfile = getattr(self, "dum" + label)
         main_phaseogram(
             [
                 evfile,
@@ -132,7 +148,10 @@ class TestPhaseogram:
         ip = run_interactive_phaseogram(evfile, 9.9, test=True, nbin=16, nt=8)
         ip.update(1)
         ip.recalculate(1)
-        ip.toa(1)
+        with pytest.warns(UserWarning) as record:
+            ip.toa(1)
+        assert np.any(["TOA calculation is not robust" in r.message.args[0]
+                       for r in record])
         ip.reset(1)
         ip.fdot = 2
         f, fdot, fddot = ip.get_values()
@@ -148,6 +167,8 @@ class TestPhaseogram:
 
         par = hen_root(evfile) + ".par"
         with open(par, "a") as fobj:
+            print("F0  9.9", file=fobj)
+            print("PEPOCH 57000", file=fobj)
             print("BINARY BT", file=fobj)
             print("PB  1e20", file=fobj)
             print("A1  0", file=fobj)
@@ -161,7 +182,10 @@ class TestPhaseogram:
         )
         ip.update(1)
         ip.recalculate(1)
-        ip.toa(1)
+        with pytest.warns(UserWarning) as record:
+            ip.toa(1)
+        assert np.any(["TOA calculation is not robust" in r.message.args[0]
+                       for r in record])
         ip.reset(1)
         ip.fdot = 2
         f, fdot, fddot = ip.get_values()