Merge pull request #158 from StingraySoftware/use_new_stingray

Use newest stingray versions

.github/workflows/ci_test.yml

@@ -103,6 +103,10 @@ jobs:
             python: '3.10'
             tox_env: 'py310-test-alldeps-cov'
             use_remote_data: true
+          - os: macos-14
+            python: '3.10'
+            tox_env: 'py310-test-alldeps-cov'
+            use_remote_data: true
           - os: ubuntu-latest
             python: '3.10'
             tox_env: 'py310-test-devdeps'

docs/index.rst

@@ -43,14 +43,19 @@ satellite (for sure XMM-Newton and RXTE).
 What's new
 ----------
 
-HENDRICS 8.0 (under development)
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+HENDRICS 8.0.0rc1
+~~~~~~~~~~~~~~~~~
+
++ Compatible with `Stingray 2.0.0rc1 <https://github.com/StingraySoftware/stingray/releases/tag/v2.0.0rc1>`__, which introduced:
+
+    + Lomb-Scargle periodograms and cross spectra are now available
+    + Power colors
+    + Easy filling of small gaps in light curves with random data
+    + Generic timeseries (complex data, multi-dimensional data)
 
 + ``HENaccelsearch`` now has additional options for detrending, denoising and deorbiting
 + An improved Maximum likelihood algorithm as FFTFIT substitute for TOA calculation
 + NASA's IXPE added to supported missions
-+ Compatible with `Stingray 1.1 <https://github.com/StingraySoftware/stingray/releases/tag/v1.1>`__ or newer
-+ Add option to obtain weighted light curves
 + Better support of Stingray's native file formats
 
 

hendrics/binary.py

@@ -132,7 +132,8 @@ def save_events_to_binary(
             )
 
         good = (events.energy >= emin) & (events.energy < emax)
-        events.time = events.time[good]
+        events = events.apply_mask(good)
+        # events.time = events.time[good]
 
     tstop = events.gti[-1, 1]
     nbin = (tstop - tstart) / bin_time

hendrics/colors.py

@@ -5,9 +5,8 @@
 from astropy import log
 from stingray.lightcurve import Lightcurve
 import numpy as np
-from .io import HEN_FILE_EXTENSION, load_lcurve, save_lcurve
+from .io import HEN_FILE_EXTENSION, load_events, save_lcurve
 from .base import hen_root
-from .lcurve import main as henlcurve
 
 
 def colors():
@@ -27,8 +26,7 @@ def main(args=None):
         "--energies",
         nargs=4,
         required=True,
-        type=str,
-        default=None,
+        type=float,
         help="The energy boundaries in keV used to calculate "
         "the color. E.g. -e 2 3 4 6 means that the "
         "color will be calculated as 4.-6./2.-3. keV. "
@@ -46,65 +44,44 @@ def main(args=None):
     log.setLevel(args.loglevel)
 
     with log.log_to_file("HENcolors.log"):
-        option = "--energy-interval"
-        if args.use_pi:
-            option = "--pi-interval"
-
+        energies = [
+            [args.energies[0], args.energies[1]],
+            [args.energies[2], args.energies[3]],
+        ]
         if args.outfile is not None and len(files) > 1:
             raise ValueError("Specify --output only when processing " "a single file")
         for f in files:
-            henlcurve(
-                [f]
-                + [option]
-                + args.energies[:2]
-                + [
-                    "-b",
-                    str(args.bintime),
-                    "-d",
-                    ".",
-                    "-o",
-                    "lc0" + HEN_FILE_EXTENSION,
-                ]
-            )
-            lc0 = load_lcurve("lc0" + HEN_FILE_EXTENSION)
-            henlcurve(
-                [f]
-                + [option]
-                + args.energies[2:]
-                + [
-                    "-b",
-                    str(args.bintime),
-                    "-d",
-                    ".",
-                    "-o",
-                    "lc1" + HEN_FILE_EXTENSION,
-                ]
+            events = load_events(f)
+            if not args.use_pi and events.energy is None:
+                raise ValueError(
+                    "Energy information not found in file {0}. "
+                    "Use --use-pi if you want to use PI channels "
+                    "instead.".format(f)
+                )
+            h_starts, h_stops, colors, color_errs = events.get_color_evolution(
+                energy_ranges=energies, segment_size=args.bintime, use_pi=args.use_pi
             )
-            lc1 = load_lcurve("lc1" + HEN_FILE_EXTENSION)
 
-            time = lc0.time
-            counts = lc1.countrate / lc0.countrate
-            counts_err = np.sqrt(lc1.countrate_err**2 + lc0.countrate_err**2)
+            time = (h_starts + h_stops) / 2
+
             scolor = Lightcurve(
                 time=time,
-                counts=counts,
-                err=counts_err,
+                counts=colors,
+                err=color_errs,
                 input_counts=False,
                 err_dist="gauss",
-                gti=lc0.gti,
+                gti=events.gti,
                 dt=args.bintime,
+                skip_checks=True,
             )
-            del lc0
-            del lc1
-            os.unlink("lc0" + HEN_FILE_EXTENSION)
-            os.unlink("lc1" + HEN_FILE_EXTENSION)
 
             if args.outfile is None:
                 label = "_E_"
                 if args.use_pi:
                     label = "_PI_"
-                label += "{3}-{2}_over_{1}-{0}".format(*args.energies)
+                label += "{3:g}-{2:g}_over_{1:g}-{0:g}".format(*args.energies)
                 args.outfile = hen_root(f) + label + HEN_FILE_EXTENSION
             scolor.e_intervals = np.asarray([float(k) for k in args.energies])
             scolor.use_pi = args.use_pi
             save_lcurve(scolor, args.outfile, lctype="Color")
+            print(args.outfile)

hendrics/create_gti.py

@@ -58,16 +58,24 @@ def create_gti(
     # Necessary as nc variables are sometimes defined as array
     from numpy import array  # NOQA
 
-    ftype, data = get_file_type(fname, raw_data=True)
+    ftype, data = get_file_type(fname, raw_data=False)
+
+    instr = data.instr
 
-    instr = data["instr"]
     if ftype == "lc" and instr.lower() == "pca":
         warnings.warn("RXTE/PCA data; normalizing lc per no. PCUs", AstropyUserWarning)
         # If RXTE, plot per PCU count rate
-        data["counts"] /= data["nPCUs"]
-    mjdref = data["mjdref"]
+        data.counts /= data.nPCUs
+    mjdref = data.mjdref
     # Map all entries of data to local variables
-    locals().update(data)
+    array_attrs = data.array_attrs() + ["time"]
+    locals().update(zip(array_attrs, [getattr(data, attr) for attr in array_attrs]))
+    if hasattr(data, "internal_array_attrs"):
+        array_attrs = data.internal_array_attrs()
+        mod_array_attrs = [attr.replace("_", "") for attr in array_attrs]
+        locals().update(
+            zip(mod_array_attrs, [getattr(data, attr) for attr in array_attrs])
+        )
 
     good = eval(filter_expr)
 
@@ -88,31 +96,20 @@ def apply_gti(fname, gti, outname=None, minimum_length=0):
 
     File MUST have a GTI extension already, and an extension called `time`.
     """
-    ftype, data = get_file_type(fname, raw_data=True)
+    ftype, data = get_file_type(fname, raw_data=False)
 
     try:
-        datagti = data["gti"]
+        datagti = data.gti
         newgti = cross_gtis([gti, datagti])
     except Exception:  # pragma: no cover
         warnings.warn("Data have no GTI extension", AstropyUserWarning)
         newgti = gti
 
     newgti = filter_gti_by_length(newgti, minimum_length)
 
-    data["gti"] = newgti
-    good = create_gti_mask(data["time"], newgti)
-
-    data["time"] = data["time"][good]
-    if ftype == "lc":
-        data["counts"] = data["counts"][good]
-        data["counts_err"] = data["counts_err"][good]
-    elif ftype == "events":
-        for key in ["pi", "energy"]:
-            if key in data:
-                data[key] = data[key][good]
-
-        if data["instr"] == "PCA":  # pragma: no cover
-            data["PCU"] = data["PCU"][good]
+    data.gti = newgti
+    data = data.apply_gtis()
+    data._mask = None
 
     newext = "_gtifilt" + HEN_FILE_EXTENSION
     outname = _assign_value_if_none(

hendrics/exposure.py

@@ -373,11 +373,12 @@ def main(args=None):
 
         outfile = correct_lightcurve(lc_file, uf_file, outname)
 
-        outdata = load_data(outfile)
-        time = outdata["time"]
-        lc = outdata["counts"]
-        expo = outdata["expo"]
-        gti = outdata["gti"]
+        # outdata = load_data(outfile)
+        _, outdata = get_file_type(outfile, raw_data=False)
+        time = outdata.time
+        lc = outdata.counts
+        expo = outdata.expo
+        gti = outdata.gti
 
         try:
             _plot_corrected_light_curve(time, lc * expo, expo, gti, outroot)

hendrics/exvar.py

@@ -72,16 +72,16 @@ def main(args=None):
         for fname in args.files:
             lcurve = load_lcurve(fname)
             if args.norm == "fvar":
-                start, stop, res = lcurve.analyze_lc_chunks(
-                    args.chunk_length, fvar, args.fraction_step
+                start, stop, res = lcurve.analyze_segments(
+                    fvar, args.chunk_length, args.fraction_step
                 )
             elif args.norm == "excvar":
-                start, stop, res = lcurve.analyze_lc_chunks(
-                    args.chunk_length, excvar_none, args.fraction_step
+                start, stop, res = lcurve.analyze_segments(
+                    excvar_none, args.chunk_length, args.fraction_step
                 )
             elif args.norm == "norm_excvar":
-                start, stop, res = lcurve.analyze_lc_chunks(
-                    args.chunk_length, excvar_norm, args.fraction_step
+                start, stop, res = lcurve.analyze_segments(
+                    excvar_norm, args.chunk_length, args.fraction_step
                 )
             else:
                 raise ValueError("Normalization must be fvar, norm_excvar " "or excvar")

hendrics/fold.py

@@ -433,7 +433,7 @@ def get_TOAs_from_events(events, folding_length, *frequency_derivatives, **kwarg
 
         if np.any(np.isnan(pars)) or pars[0] == 0.0 or np.any(np.isnan(errs)):
             warnings.warn(
-                f"Invalid TOA in interval {start}-{stop} (idxs {startidx:stopidx}): {pars}, {errs}"
+                f"Invalid TOA in interval {start}-{stop} (idxs {startidx}:{stopidx}): {pars}, {errs}"
             )
             continue