Merge pull request #168 from StingraySoftware/safe_intv_in_read_events

Introduce safe intervals in read_events

Author: matteobachetti

.github/workflows/ci_test.yml

@@ -165,3 +165,5 @@ jobs:
       uses: codecov/codecov-action@v3
       with:
         file: ./coverage.xml
+      env:
+        CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}

docs/conf.py

@@ -28,6 +28,8 @@
 import os
 import sys
 from importlib import import_module
+import warnings
+import subprocess as sp
 
 try:
     from sphinx_astropy.conf.v1 import *  # noqa
@@ -41,13 +43,18 @@
     import tomli as tomllib
 from pathlib import Path
 
+warnings.filterwarnings("ignore", message=".*recommended numba.*")
+warnings.filterwarnings("ignore", message=".*Using pickle.*")
+warnings.filterwarnings("ignore", message=".*PINT is not installed.*")
+
 ON_RTD = os.environ.get("READTHEDOCS") == "True"
 ON_TRAVIS = os.environ.get("TRAVIS") == "true"
 
 # Grab minversion from pyproject.toml
 with (Path(__file__).parents[1] / "pyproject.toml").open("rb") as f:
     pyproject = tomllib.load(f)
 
+
 # -- General configuration ----------------------------------------------------
 
 # By default, highlight as Python 3.
@@ -203,8 +210,11 @@
 if not ON_RTD and not ON_TRAVIS:
     # scripts = dict(conf.items("options.entry_points"))["console_scripts"]
     scripts = pyproject["project"]["scripts"]
+    pwarn = os.getenv("PYTHONWARNINGS")
 
-    import subprocess as sp
+    os.environ["PYTHONWARNINGS"] = "ignore"
+    cols = os.getenv("COLUMNS")
+    os.environ["COLUMNS"] = "80"
 
     cli_file = os.path.join(os.getcwd(), "scripts", "cli.rst")
     if os.path.exists(cli_file):
@@ -213,6 +223,7 @@
             print("""======================\n""", file=fobj)
 
             for cl in sorted(scripts.keys()):
+                print(f"Writing help for {cl}")
                 if cl.startswith("MP"):
                     continue
                 print(cl, file=fobj)
@@ -227,3 +238,13 @@
                     else:
                         print("    " + l, file=fobj)
                 print(file=fobj)
+
+    if cols is not None:
+        os.environ["COLUMNS"] = cols
+    else:
+        del os.environ["COLUMNS"]
+
+    if pwarn is not None:
+        os.environ["PYTHONWARNINGS"] = pwarn
+    else:
+        del os.environ["PYTHONWARNINGS"]

docs/index.rst

@@ -43,6 +43,16 @@ satellite (for sure XMM-Newton and RXTE).
 What's new
 ----------
 
+HENDRICS 8.2 (under development)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
++ Brings all bugfixes coming with `Stingray 2.2.4 <https://github.com/StingraySoftware/stingray/releases/tag/v2.2.4>`__
++ New option for ``HENreadevents``: ``--safe-interval``, allowing to decrease GTI length by fixed amounts at the start and the end (useful, e.g., when observations show artifacts entering and exiting from occultation or SAA)
+
+Internal changes:
+
++ Improvements with docs creation, in particular to the creation of the cli.rst file.
+
 HENDRICS 8.1
 ~~~~~~~~~~~~
 

docs/scripts/cli.rst

@@ -77,27 +77,27 @@
 
 
 __all__ = [
+    "calculate_shifts",
     "check_phase_error_after_casting_to_double",
     "decide_binary_parameters",
-    "folding_orbital_search",
+    "dyn_folding_search",
     "fit",
-    "calculate_shifts",
-    "mod",
-    "shift_and_sum",
-    "z_n_fast",
-    "transient_search",
-    "plot_transient_search",
-    "search_with_qffa_step",
-    "search_with_qffa",
-    "search_with_ffa",
+    "folding_orbital_search",
     "folding_search",
-    "dyn_folding_search",
+    "h_test",
+    "main_accelsearch",
     "main_efsearch",
+    "main_z2vspf",
     "main_zsearch",
+    "mod",
+    "plot_transient_search",
+    "search_with_ffa",
+    "search_with_qffa",
+    "search_with_qffa_step",
+    "shift_and_sum",
+    "transient_search",
     "z2_vs_pf",
-    "main_z2vspf",
-    "main_accelsearch",
-    "h_test",
+    "z_n_fast",
 ]
 
 

hendrics/efsearch.py

@@ -3,6 +3,7 @@
 
 import os
 import warnings
+from collections.abc import Iterable
 
 import numpy as np
 from stingray.events import EventList
@@ -33,6 +34,7 @@ def treat_event_file(
     additional_columns=None,
     fill_small_gaps=None,
     bin_time_for_occultations=None,
+    safe_interval=None,
 ):
     """Read data from an event file, with no external GTI information.
 
@@ -63,6 +65,9 @@ def treat_event_file(
         Create a light curve with this bin time and infer occultations not
         recorded in GTIs. The rule is that the flux drops to zero and the average
         counts per bin are significantly above 25 ct/s.
+    safe_interval : float or [float, float]
+        A safe interval to exclude at both ends (if single float) or the start
+        and the end (if pair of values) of GTIs.
     """
     # gtistring = assign_value_if_none(gtistring, "GTI,GTI0,STDGTI")
     log.info(f"Opening {filename}")
@@ -85,9 +90,9 @@ def treat_event_file(
     if bin_time_for_occultations is not None and bin_time_for_occultations > 0:
         lc = events.to_lc(bin_time_for_occultations)
         meanrate = np.median(lc.counts)
-
         if meanrate > 25:
-            good_gti = create_gti_mask(lc.time, lc.gti)
+            log.info("Filtering out occultations")
+            good_gti = create_gti_mask(lc.time, lc.gti, safe_interval=safe_interval)
             good = lc.counts > 0
             new_bad = (~good) & good_gti
             if np.any(new_bad):
@@ -96,8 +101,13 @@ def treat_event_file(
                     lc.time, (good_gti & good), safe_interval=bin_time_for_occultations
                 )
                 events.gti = gti
+    elif safe_interval is not None:
+        if not isinstance(safe_interval, Iterable):
+            safe_interval = [safe_interval, safe_interval]
+        gti[:, 0] += safe_interval[0]
+        gti[:, 1] -= safe_interval[1]
 
-    lengths = np.array([g1 - g0 for (g0, g1) in gti])
+    lengths = gti[:, 1] - gti[:, 0]
     gti = gti[lengths >= min_length]
     events.gti = gti
     detector_id = events.detector_id
@@ -150,9 +160,7 @@ def treat_event_file(
 
                 good_gti = cross_two_gtis([g], gti)
                 if noclobber and os.path.exists(outfile_local):
-                    warnings.warn(
-                        f"{outfile_local} exists, " + "and noclobber option used. Skipping"
-                    )
+                    warnings.warn(f"{outfile_local} exists, and noclobber option used. Skipping")
                     return
                 good = np.logical_and(events.time >= g[0], events.time < g[1])
                 all_good = good_det & good
@@ -453,7 +461,7 @@ def main_join(args=None):
     import argparse
 
     description = (
-        "Read a cleaned event files and saves the relevant " "information in a standard format"
+        "Read a cleaned event files and saves the relevant information in a standard format"
     )
     parser = argparse.ArgumentParser(description=description)
     parser.add_argument("files", help="Files to join", type=str, nargs="+")
@@ -526,7 +534,7 @@ def main(args=None):
     from .base import _add_default_args, check_negative_numbers_in_args
 
     description = (
-        "Read a cleaned event files and saves the relevant " "information in a standard format"
+        "Read a cleaned event files and saves the relevant information in a standard format"
     )
     parser = argparse.ArgumentParser(description=description)
     parser.add_argument("files", help="List of files", nargs="+")
@@ -586,6 +594,13 @@ def main(args=None):
         "XMM)",
         default=None,
     )
+    parser.add_argument(
+        "--safe-interval",
+        nargs=2,
+        type=float,
+        default=0,
+        help="Interval at start and stop of GTIs used" + " for filtering",
+    )
     parser.add_argument(
         "--fill-small-gaps",
         type=float,
@@ -623,6 +638,7 @@ def main(args=None):
             "fill_small_gaps": args.fill_small_gaps,
             "split_by_detector": not args.ignore_detectors,
             "bin_time_for_occultations": args.bin_time_for_occultations,
+            "safe_interval": args.safe_interval,
         }
 
         arglist = [[f, argdict] for f in files]

hendrics/read_events.py

@@ -371,6 +371,9 @@ def test_load_events_split(self):
             "monol_testB_nustar_fpmb_gti000_ev" + HEN_FILE_EXTENSION,
         )
         assert os.path.exists(new_filename)
+        with pytest.warns(UserWarning, match=".*exists, and noclobber option"):
+            read_events.main(command.split() + ["--noclobber"])
+
         command = f"{new_filename}"
         lcurve.main(command.split())
         new_filename = os.path.join(

hendrics/tests/test_read_events.py

@@ -481,3 +481,7 @@ ignore-words-list = """
     wrap-descriptions = 75
     black = true
     syntax = "numpy"
+
+
+[tool.black]
+line-length = 100