From 07bbe247eda5162342c7f5a482feca39147b0e16 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 23 Apr 2024 15:36:34 +0200 Subject: [PATCH 01/14] Fix operation on segment when result is NaN --- stingray/base.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/stingray/base.py b/stingray/base.py index a70fed444..da51305f3 100644 --- a/stingray/base.py +++ b/stingray/base.py @@ -2646,24 +2646,27 @@ def analyze_segments(self, func, segment_size, fraction_step=1, **kwargs): stop = np.searchsorted(self.time, stop_times) results = [] + good = np.ones(start.size, dtype=bool) + for i, (st, sp, tst, tsp) in enumerate(zip(start, stop, start_times, stop_times)): if sp - st <= 1: warnings.warn( f"Segment {i} ({tst}--{tsp}) has one data point or less. Skipping it." ) - res = np.nan - else: - lc_filt = self[st:sp] - lc_filt.gti = np.asarray([[tst, tsp]]) + good[i] = False + continue + lc_filt = self[st:sp] + lc_filt.gti = np.asarray([[tst, tsp]]) - res = func(lc_filt, **kwargs) + res = func(lc_filt, **kwargs) results.append(res) results = np.array(results) if len(results.shape) == 2: results = [results[:, i] for i in range(results.shape[1])] - return start_times, stop_times, results + + return start_times[good], stop_times[good], results def analyze_by_gti(self, func, fraction_step=1, **kwargs): """Analyze the light curve with any function, on a GTI-by-GTI base. From 235273e91f050a0a8ca962b7c288316ed8ad0309 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 23 Apr 2024 15:41:20 +0200 Subject: [PATCH 02/14] Add changelog --- docs/changes/822.bugfix.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 docs/changes/822.bugfix.rst diff --git a/docs/changes/822.bugfix.rst b/docs/changes/822.bugfix.rst new file mode 100644 index 000000000..06f526c55 --- /dev/null +++ b/docs/changes/822.bugfix.rst @@ -0,0 +1 @@ +Fix case when analyze_segments has an invalid segment From d5244db279b48dc9f2be7998d26f9dd8168bca6d Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 23 Apr 2024 16:33:43 +0200 Subject: [PATCH 03/14] Fix test for new behavior --- stingray/tests/test_base.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/stingray/tests/test_base.py b/stingray/tests/test_base.py index 37128486d..b909cfa73 100644 --- a/stingray/tests/test_base.py +++ b/stingray/tests/test_base.py @@ -1500,9 +1500,7 @@ def func(x): # I do not specify the segment_size, which means results will be calculated per-GTI with pytest.warns(UserWarning, match="has one data point or less."): - _, _, results = ts.analyze_segments(func, segment_size=None) - # the first GTI contains only one bin, the result will be invalid - assert np.isnan(results[0]) + ts.analyze_segments(func, segment_size=None) def test_analyze_segments_by_gti(self): ts = StingrayTimeseries(time=np.arange(11), dt=1, gti=[[-0.5, 5.5], [6.5, 10.5]]) From f66afd708863b6d4067621a2ce2c7468c75cf991 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Wed, 24 Apr 2024 08:50:22 +0200 Subject: [PATCH 04/14] Avoid divisionbyzero error --- stingray/events.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/stingray/events.py b/stingray/events.py index ff98ba0eb..21dead519 100644 --- a/stingray/events.py +++ b/stingray/events.py @@ -797,6 +797,8 @@ def color(ev): en1_ct = np.count_nonzero(mask1) en2_ct = np.count_nonzero(mask2) + if en1_ct == 0 or en2_ct == 0: + return np.nan, np.nan color = en2_ct / en1_ct color_err = color * (np.sqrt(en1_ct) / en1_ct + np.sqrt(en2_ct) / en2_ct) return color, color_err From 7187ad461a9ea2bfc4f6474f6c2c93f78d3e8e7d Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Wed, 24 Apr 2024 11:38:25 +0200 Subject: [PATCH 05/14] Warn when bad data happen, and test for it --- docs/notebooks | 2 +- stingray/events.py | 1 + stingray/tests/test_events.py | 6 ++++++ 3 files changed, 8 insertions(+), 1 deletion(-) diff --git a/docs/notebooks b/docs/notebooks index 0e34107ac..2db8e50be 160000 --- a/docs/notebooks +++ b/docs/notebooks @@ -1 +1 @@ -Subproject commit 0e34107ac872db8ed9b4acfcc75b97edf07f3ccb +Subproject commit 2db8e50be46e3d5cddebb2c15647c45d63de4a92 diff --git a/stingray/events.py b/stingray/events.py index 21dead519..5ad0b5ed4 100644 --- a/stingray/events.py +++ b/stingray/events.py @@ -798,6 +798,7 @@ def color(ev): en2_ct = np.count_nonzero(mask2) if en1_ct == 0 or en2_ct == 0: + warnings.warn("No counts in one of the energy ranges. Returning NaN") return np.nan, np.nan color = en2_ct / en1_ct color_err = color * (np.sqrt(en1_ct) / en1_ct + np.sqrt(en2_ct) / en2_ct) diff --git a/stingray/tests/test_events.py b/stingray/tests/test_events.py index 3f909edb6..3344367d9 100644 --- a/stingray/tests/test_events.py +++ b/stingray/tests/test_events.py @@ -644,6 +644,12 @@ def test_colors(self): assert np.allclose(start, np.arange(10) * 10000) assert np.allclose(stop, np.arange(1, 11) * 10000) + def test_colors_missing_energies(self): + events = copy.deepcopy(self.events) + events.filter_energy_range([0, 3], inplace=True) + with pytest.warns(UserWarning, match="No counts in one of the energy ranges"): + events.get_color_evolution([[0, 3], [4, 6]], 10000) + def test_colors_no_segment(self): start, stop, colors, color_errs = self.events.get_color_evolution([[0, 3], [4, 6]]) # 50000 / 50000 = 1 From 3781a18a4c8108bf60e51871e43dbd20af1dc018 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Fri, 26 Apr 2024 13:33:08 +0200 Subject: [PATCH 06/14] Fix analyze_segments for function with multiple return values --- stingray/base.py | 32 +++++++++++++++++++++----------- 1 file changed, 21 insertions(+), 11 deletions(-) diff --git a/stingray/base.py b/stingray/base.py index da51305f3..e901150ed 100644 --- a/stingray/base.py +++ b/stingray/base.py @@ -2606,8 +2606,10 @@ def analyze_segments(self, func, segment_size, fraction_step=1, **kwargs): Lower time boundaries of all time segments. stop_times : array upper time boundaries of all segments. - result : array of N elements - The result of ``func`` for each segment of the light curve + result : list of N elements + The result of ``func`` for each segment of the light curve. If the function + returns multiple outputs, they are returned as a list of arrays. + If a given interval has not enough data for a calculation, ``None`` is returned. Examples -------- @@ -2646,27 +2648,35 @@ def analyze_segments(self, func, segment_size, fraction_step=1, **kwargs): stop = np.searchsorted(self.time, stop_times) results = [] - good = np.ones(start.size, dtype=bool) + n_outs = 1 for i, (st, sp, tst, tsp) in enumerate(zip(start, stop, start_times, stop_times)): if sp - st <= 1: warnings.warn( f"Segment {i} ({tst}--{tsp}) has one data point or less. Skipping it." ) - good[i] = False + results.append(None) continue lc_filt = self[st:sp] lc_filt.gti = np.asarray([[tst, tsp]]) res = func(lc_filt, **kwargs) results.append(res) - - results = np.array(results) - - if len(results.shape) == 2: - results = [results[:, i] for i in range(results.shape[1])] - - return start_times[good], stop_times[good], results + if isinstance(res, Iterable): + n_outs = len(res) + + # If the function returns multiple outputs, we need to separate them + + if n_outs > 1: + outs = [[] for _ in range(n_outs)] + for res in results: + for i in range(n_outs): + if res is None: + outs[i] = None + else: + outs[i].append(res[i]) + results = outs + return start_times, stop_times, results def analyze_by_gti(self, func, fraction_step=1, **kwargs): """Analyze the light curve with any function, on a GTI-by-GTI base. From 6a7ed2efc11b30a47402087362259cb4db99d83e Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Fri, 26 Apr 2024 13:43:48 +0200 Subject: [PATCH 07/14] Try to make the output a numpy array --- stingray/base.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/stingray/base.py b/stingray/base.py index e901150ed..f5740bc53 100644 --- a/stingray/base.py +++ b/stingray/base.py @@ -2676,6 +2676,13 @@ def analyze_segments(self, func, segment_size, fraction_step=1, **kwargs): else: outs[i].append(res[i]) results = outs + + # Try to transform into a (possibly multi-dimensional) numpy array + try: + results = np.array(results) + except ValueError: # pragma: no cover + pass + return start_times, stop_times, results def analyze_by_gti(self, func, fraction_step=1, **kwargs): From d80b70b1f9b109012829e520ef67ca61763bb59c Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Fri, 26 Apr 2024 13:47:17 +0200 Subject: [PATCH 08/14] Verify kind of iterable --- stingray/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stingray/base.py b/stingray/base.py index f5740bc53..34384f73b 100644 --- a/stingray/base.py +++ b/stingray/base.py @@ -2662,7 +2662,7 @@ def analyze_segments(self, func, segment_size, fraction_step=1, **kwargs): res = func(lc_filt, **kwargs) results.append(res) - if isinstance(res, Iterable): + if isinstance(res, Iterable) and not isinstance(res, str): n_outs = len(res) # If the function returns multiple outputs, we need to separate them From d8a3e6f91c3c2f55f50c64f53bc235d0b55eaec7 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Fri, 26 Apr 2024 13:57:41 +0200 Subject: [PATCH 09/14] Test properly and be more robust in case of missing values --- stingray/base.py | 2 +- stingray/tests/test_base.py | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/stingray/base.py b/stingray/base.py index 34384f73b..2d71f1b91 100644 --- a/stingray/base.py +++ b/stingray/base.py @@ -2672,7 +2672,7 @@ def analyze_segments(self, func, segment_size, fraction_step=1, **kwargs): for res in results: for i in range(n_outs): if res is None: - outs[i] = None + outs[i].append(None) else: outs[i].append(res[i]) results = outs diff --git a/stingray/tests/test_base.py b/stingray/tests/test_base.py index b909cfa73..89ab1552a 100644 --- a/stingray/tests/test_base.py +++ b/stingray/tests/test_base.py @@ -1492,11 +1492,14 @@ def test_estimate_segment_size_lower_dt(self): assert ts.estimate_segment_size(100, min_samples=40) == 8.0 - def test_analyze_segments_bad_intv(self): + @pytest.mark.parametrize("n_outs", [0, 1, 2, 3]) + def test_analyze_segments_bad_intv(self, n_outs): ts = StingrayTimeseries(time=np.arange(10), dt=1, gti=[[-0.5, 0.5], [1.5, 10.5]]) def func(x): - return np.size(x.time) + if n_outs == 0: + return np.size(x.time) + return [np.size(x.time) for _ in range(n_outs)] # I do not specify the segment_size, which means results will be calculated per-GTI with pytest.warns(UserWarning, match="has one data point or less."): From 54d26bc22f052d1e7d253a5bd994f26f1157f5cf Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Fri, 26 Apr 2024 14:06:21 +0200 Subject: [PATCH 10/14] Solve old numpy warning --- setup.cfg | 1 + 1 file changed, 1 insertion(+) diff --git a/setup.cfg b/setup.cfg index 82852c01f..9a480a5e4 100644 --- a/setup.cfg +++ b/setup.cfg @@ -116,6 +116,7 @@ filterwarnings = ignore:.*datetime.datetime.utcfromtimestamp.*:DeprecationWarning ignore:.*__array_wrap__ must accept context and return_scalar arguments:DeprecationWarning ignore:.*Pyarrow: + ignore:.*Creating an ndarray from ragged nested sequences: ;addopts = --disable-warnings From 8658a14f940ef5e64debdb90dafef7285b9c12db Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 7 May 2024 09:50:42 +0200 Subject: [PATCH 11/14] Fixing warning --- stingray/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stingray/base.py b/stingray/base.py index 2d71f1b91..1c8c1583e 100644 --- a/stingray/base.py +++ b/stingray/base.py @@ -2653,7 +2653,7 @@ def analyze_segments(self, func, segment_size, fraction_step=1, **kwargs): for i, (st, sp, tst, tsp) in enumerate(zip(start, stop, start_times, stop_times)): if sp - st <= 1: warnings.warn( - f"Segment {i} ({tst}--{tsp}) has one data point or less. Skipping it." + f"Segment {i} ({tst}--{tsp}) has one data point or less. " ) results.append(None) continue From f5c18b3fbfcb64e123a402653c13b1b0823e908b Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 7 May 2024 10:06:33 +0200 Subject: [PATCH 12/14] Skip invalid intervals --- stingray/base.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/stingray/base.py b/stingray/base.py index 1c8c1583e..4c50d0dcf 100644 --- a/stingray/base.py +++ b/stingray/base.py @@ -2653,9 +2653,8 @@ def analyze_segments(self, func, segment_size, fraction_step=1, **kwargs): for i, (st, sp, tst, tsp) in enumerate(zip(start, stop, start_times, stop_times)): if sp - st <= 1: warnings.warn( - f"Segment {i} ({tst}--{tsp}) has one data point or less. " + f"Segment {i} ({tst}--{tsp}) has one data point or less. Skipping it " ) - results.append(None) continue lc_filt = self[st:sp] lc_filt.gti = np.asarray([[tst, tsp]]) @@ -2671,10 +2670,7 @@ def analyze_segments(self, func, segment_size, fraction_step=1, **kwargs): outs = [[] for _ in range(n_outs)] for res in results: for i in range(n_outs): - if res is None: - outs[i].append(None) - else: - outs[i].append(res[i]) + outs[i].append(res[i]) results = outs # Try to transform into a (possibly multi-dimensional) numpy array From 4330c2f3c46d2e879d8211261bd1f934f0b16ddd Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 7 May 2024 10:07:47 +0200 Subject: [PATCH 13/14] Add skip information [docs only] --- stingray/base.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/stingray/base.py b/stingray/base.py index 4c50d0dcf..fff9b135f 100644 --- a/stingray/base.py +++ b/stingray/base.py @@ -2579,6 +2579,8 @@ def estimate_segment_size(self, min_counts=None, min_samples=None, even_sampling def analyze_segments(self, func, segment_size, fraction_step=1, **kwargs): """Analyze segments of the light curve with any function. + Intervals with less than two data points are skipped. + Parameters ---------- func : function From ac86944c645ca485e159f13c744ea72fdeeb0c9f Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 7 May 2024 10:42:17 +0200 Subject: [PATCH 14/14] Fix docstring --- stingray/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stingray/base.py b/stingray/base.py index fff9b135f..3e499e5b7 100644 --- a/stingray/base.py +++ b/stingray/base.py @@ -2579,7 +2579,7 @@ def estimate_segment_size(self, min_counts=None, min_samples=None, even_sampling def analyze_segments(self, func, segment_size, fraction_step=1, **kwargs): """Analyze segments of the light curve with any function. - Intervals with less than two data points are skipped. + Intervals with less than one data point are skipped. Parameters ----------