Repository: PlasmaPy/PlasmaPy

Issue:
Excess statistics for time series
### Feature description

For a given time series and given threshold values, find the excess statistics of the time series (total time above threshold, number of upwards crossings, average time above threshold, etc.).

### Motivation

This functionality has been developed and applied for time series analysis of intermittent fluctuations in the SOL region of fusion experiments. Some examples are listed below:
https://doi.org/10.1016/j.nme.2016.11.008
https://doi.org/10.1063/1.4947235
https://doi.org/10.1103/PhysRevE.97.012110

### Implementation strategy

Implementing this feature should be straight forward since the method is already developed and available in: https://github.com/uit-cosmo/fpp-analysis-tools/blob/main/fppanalysis/excess_statistics.py

The code still needs refactoring and appropriate tests and documentation needs to be added. I am planning to take a look at this in the upcoming days (could I be assigned to this issue?).

### Additional context

_No response_

Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1984?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) Report
Base: **98.32**% // Head: **98.33**% // Increases project coverage by **`+0.00%`** :tada:
> Coverage data is based on head [(`f34736c`)](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1984?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) compared to base [(`fc12bbd`)](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/fc12bbd286ae9db82ab313d8c380808b814a462c?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy).
> Patch coverage: 100.00% of modified lines in pull request are covered.

> :exclamation: Current head f34736c differs from pull request most recent head ad051b3. Consider uploading reports for the commit ad051b3 to get more accurate results

<details><summary>Additional details and impacted files</summary>


```diff
@@           Coverage Diff           @@
##             main    #1984   +/-   ##
=======================================
  Coverage   98.32%   98.33%           
=======================================
  Files          96       97    +1     
  Lines        8500     8551   +51     
=======================================
+ Hits         8358     8409   +51     
  Misses        142      142           
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1984?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) | Coverage Δ | |
|---|---|---|
| [plasmapy/analysis/time\_series/\_\_init\_\_.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1984?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy#diff-cGxhc21hcHkvYW5hbHlzaXMvdGltZV9zZXJpZXMvX19pbml0X18ucHk=) | `100.00% <100.00%> (ø)` | |
| [plasmapy/analysis/time\_series/excess\_statistics.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1984?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy#diff-cGxhc21hcHkvYW5hbHlzaXMvdGltZV9zZXJpZXMvZXhjZXNzX3N0YXRpc3RpY3MucHk=) | `100.00% <100.00%> (ø)` | |

Help us with your feedback. Take ten seconds to tell us [how you rate us](https://about.codecov.io/nps?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy). Have a feature suggestion? [Share it here.](https://app.codecov.io/gh/feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy)

</details>

[:umbrella: View full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1984?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy).   
:loudspeaker: Do you have feedback about the report comment? [Let us know in this issue](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy).

gregordecristoforo: Thanks for your feedback @namurphy !

> What would you think of having something like an `ExcessStatistics` class, with methods for each of the quantities being returned? 

That sounds like a great idea! I reformatted the code accordingly. Currently all statistics are calculate when an `ExcessStatistics` object gets instantiated sine the calculations of the statistical properties are quite entangled. What do you think about the current form?

> What would you think of building on the functionality of [astropy.timeseries](https://docs.astropy.org/en/stable/timeseries/index.html)?

Do you suggest to rewrite the code so that one can instantiate a `ExcessStatistics` object with a `TimeSeries` object as an argument? Should `ExcessStatistics` be able to handle both `astropy.timeseries` and a simple sequence (as it does now) or would you suggest writing a separate function/class for it?
namurphy: Thank you for following up on this!  I think the linter issue suddenly came up because of a new release of one of our linter packages, but we just addressed that in #2046.  I think you'll need to do a `git pull` from the `main` branch.
gregordecristoforo: Thanks! That resolved the issue.
gregordecristoforo: Thanks @namurphy for your suggestions! Expanding this class to include Astropy Quantity objects sounds like an interesting idea for a future PR.

May I ask for your help one last time? Do you know the reason for the `build_docs_pins ERROR`?
namurphy: Yes, I am all too familiar with this error.  This is one example where Sphinx gives a not-very-helpful error message from a problem with our documentation build.

We just created a new [documentation troubleshooting guide](https://docs.plasmapy.org/en/latest/contributing/doc_guide.html). This paragraph is the most relevant.  I'll submit a pull request doing what's in the last sentence, which I thought I had done before.

> This warning sometimes occurs in the type specification of a [parameter](https://docs.python.org/3/glossary.html#term-parameter) in a docstring. Sphinx attempts to link words in type specifications to code objects. Type lines are intended to provide concise information about allowed types, sizes, shapes, physical types, and default values of a parameter. To resolve this warning, first move information about the meaning of a parameter from the type specification into the parameter description that begins on the following line. To expand the list of allowed words or patterns in type specifications, add a regular expression to `nitpick_ignore_regex` in [`docs/conf.py`](https://github.com/PlasmaPy/PlasmaPy/blob/main/docs/conf.py).

Patch:
diff --git a/changelog/1984.feature.rst b/changelog/1984.feature.rst
new file mode 100644
index 00000000..ce4c244c
--- /dev/null
+++ b/changelog/1984.feature.rst
@@ -0,0 +1,3 @@
+Added the `plasmapy.analysis.time_series.excess_statistics` module
+including the `~plasmapy.analysis.time_series.excess_statistics.ExcessStatistics`
+class for calculating excess statistics of time series.
diff --git a/docs/api_static/plasmapy.analysis.time_series.excess_statistics.rst b/docs/api_static/plasmapy.analysis.time_series.excess_statistics.rst
new file mode 100644
index 00000000..13a141ae
--- /dev/null
+++ b/docs/api_static/plasmapy.analysis.time_series.excess_statistics.rst
@@ -0,0 +1,8 @@
+:orphan:
+
+`plasmapy.analysis.time_series.excess_statistics`
+=================================================
+
+.. currentmodule:: plasmapy.analysis.time_series.excess_statistics
+
+.. automodapi:: plasmapy.analysis.time_series.excess_statistics
diff --git a/plasmapy/analysis/time_series/__init__.py b/plasmapy/analysis/time_series/__init__.py
index c3cb7893..696d9fa4 100644
--- a/plasmapy/analysis/time_series/__init__.py
+++ b/plasmapy/analysis/time_series/__init__.py
@@ -1,3 +1,4 @@
 """Subpackage containing routines for analyzing time series."""
 
+from plasmapy.analysis.time_series.excess_statistics import ExcessStatistics
 from plasmapy.analysis.time_series.running_moments import running_mean, running_moment
diff --git a/plasmapy/analysis/time_series/excess_statistics.py b/plasmapy/analysis/time_series/excess_statistics.py
new file mode 100644
index 00000000..5e86cbff
--- /dev/null
+++ b/plasmapy/analysis/time_series/excess_statistics.py
@@ -0,0 +1,209 @@
+"""Functionality to calculate excess statistics of time series."""
+
+__all__ = ["ExcessStatistics"]
+
+
+import numbers
+import numpy as np
+
+from collections.abc import Iterable
+
+
+class ExcessStatistics:
+    """
+    Calculate total time, number of upwards crossings, average time and
+    root-mean-square time above given thresholds of a sequence.
+
+    Parameters
+    ----------
+    signal : 1D |array_like|
+        Signal to be analyzed.
+
+    thresholds : 1D |array_like|
+        Threshold values.
+
+    time_step : int
+        Time step of ``signal``.
+
+    Raises
+    ------
+    `ValueError`
+        If ``time_step`` ≤ 0.
+
+    Example
+    -------
+    >>> from plasmapy.analysis.time_series.excess_statistics import ExcessStatistics
+    >>> signal = [0, 0, 2, 2, 0, 4]
+    >>> thresholds = [1, 3, 5]
+    >>> time_step = 1
+    >>> excess_statistics = ExcessStatistics(signal, thresholds, time_step)
+    >>> excess_statistics.total_time_above_threshold
+    [3, 1, 0]
+    >>> excess_statistics.number_of_crossings
+    [2, 1, 0]
+    >>> excess_statistics.average_times
+    [1.5, 1.0, 0]
+    >>> excess_statistics.rms_times
+    [0.5, 0.0, 0]
+    """
+
+    def __init__(self, signal, thresholds, time_step):
+        if time_step <= 0:
+            raise ValueError("time_step must be positive")
+
+        # make sure thresholds is an iterable
+        if not isinstance(thresholds, Iterable):
+            thresholds = [thresholds]
+
+        self._total_time_above_threshold = []
+        self._number_of_crossings = []
+        self._average_times = []
+        self._rms_times = []
+        self.events_per_threshold = {}
+
+        self._calculate_excess_statistics(signal, thresholds, time_step)
+
+    def _calculate_excess_statistics(self, signal, thresholds, time_step):
+        for threshold in thresholds:
+            indices_above_threshold = np.where(np.array(signal) > threshold)[0]
+
+            if len(indices_above_threshold) == 0:
+                self._times_above_threshold = []
+                self._total_time_above_threshold.append(0)
+                self._number_of_crossings.append(0)
+                self._average_times.append(0)
+                self._rms_times.append(0)
+
+            else:
+                self._total_time_above_threshold.append(
+                    time_step * len(indices_above_threshold)
+                )
+
+                distances_to_next_index = (
+                    indices_above_threshold[1:] - indices_above_threshold[:-1]
+                )
+                split_indices = np.where(distances_to_next_index != 1)[0]
+                event_lengths = np.split(distances_to_next_index, split_indices)
+
+                # set correct length for first event
+                event_lengths[0] = np.append(event_lengths[0], 1)
+
+                self._times_above_threshold = [
+                    time_step * len(event_lengths[i])
+                    for i in range(0, len(event_lengths))
+                ]
+
+                self._number_of_crossings.append(len(event_lengths))
+                if indices_above_threshold[0] == 0:
+                    # Don't count the first event if there is no crossing.
+                    self._number_of_crossings[-1] -= 1
+
+                self._average_times.append(np.mean(self._times_above_threshold))
+                self._rms_times.append(np.std(self._times_above_threshold))
+
+            self.events_per_threshold.update({threshold: self._times_above_threshold})
+
+    def hist(self, bins=32):
+        """
+        Computes the probability density function of the time above each value
+        in ``thresholds``.
+
+        Parameters
+        ----------
+        bins : int, default: 32
+            The number of bins in the estimation of the PDF above ``thresholds``.
+
+        Returns
+        -------
+        hist: 2D `~numpy.ndarray`, shape (``thresholds.size``, ``bins`` )
+            For each value in ``thresholds``, returns the estimated PDF of time
+            above threshold.
+
+        bin_centers: 2D `~numpy.ndarray`, shape (``thresholds.size``, ``bins`` )
+            Bin centers for ``hist``.
+
+        Raises
+        ------
+        `TypeError`
+            If ``bins`` is not a positive integer.
+
+        Examples
+        --------
+        >>> from plasmapy.analysis.time_series.excess_statistics import ExcessStatistics
+        >>> signal = [0, 0, 2, 0, 4]
+        >>> thresholds = [1, 3, 5]
+        >>> time_step = 1
+        >>> excess_statistics = ExcessStatistics(signal, thresholds, time_step)
+        >>> excess_statistics.hist(2)
+        (array([[0., 2.],
+           [0., 2.],
+           [0., 0.]]), array([[0.75, 1.25],
+           [0.75, 1.25],
+           [0.  , 0.  ]]))
+        """
+
+        if not isinstance(bins, numbers.Integral):
+            raise TypeError("bins must be an integer")
+
+        hist = np.zeros((len(self.events_per_threshold), bins))
+        bin_centers = np.zeros((len(self.events_per_threshold), bins))
+
+        for i, threshold in enumerate(self.events_per_threshold.keys()):
+            if len(self.events_per_threshold[threshold]) >= 1:
+                hist[i, :], bin_edges = np.histogram(
+                    self.events_per_threshold[threshold], bins=bins, density=True
+                )
+                bin_centers[i, :] = (bin_edges[1:] + bin_edges[:-1]) / 2
+        return hist, bin_centers
+
+    @property
+    def total_time_above_threshold(self):
+        """
+        Total time above threshold(s).
+
+        Returns
+        -------
+        total_time_above_threshold: 1D |array_like|
+            Total time above threshold for each value in ``thresholds``.
+        """
+
+        return self._total_time_above_threshold
+
+    @property
+    def number_of_crossings(self):
+        """
+        Total number of upwards crossings for threshold(s).
+
+        Returns
+        -------
+        number_of_crossings: 1D |array_like|
+            Total number of upwards crossings for each value in ``thresholds``.
+        """
+
+        return self._number_of_crossings
+
+    @property
+    def average_times(self):
+        """
+        Average time above threshold(s).
+
+        Returns
+        -------
+        average_times: 1D |array_like|
+            Average time above each value in ``thresholds``.
+        """
+
+        return self._average_times
+
+    @property
+    def rms_times(self):
+        """
+        Root-mean-square values of time above threshold(s).
+
+        Returns
+        -------
+        rms_times: 1D |array_like|
+            Root-mean-square values of time above each value in ``thresholds``.
+        """
+
+        return self._rms_times


Test patch:
diff --git a/plasmapy/analysis/time_series/tests/test_excess_statistics.py b/plasmapy/analysis/time_series/tests/test_excess_statistics.py
new file mode 100644
index 00000000..4a2753a4
--- /dev/null
+++ b/plasmapy/analysis/time_series/tests/test_excess_statistics.py
@@ -0,0 +1,99 @@
+"""Tests for excess_statistics.py"""
+
+import astropy.units as u
+import numpy as np
+import pytest
+
+from plasmapy.analysis.time_series.excess_statistics import ExcessStatistics
+
+
+@pytest.mark.parametrize(
+    "signal, thresholds, time_step, pdf, bins, expected",
+    [
+        (
+            [0, 2, 0],
+            1,
+            1,
+            False,
+            32,
+            ([1], [1], [1], [0]),
+        ),
+        (
+            [0, 2, 0],
+            1,
+            0.5,
+            False,
+            32,
+            ([0.5], [1], [0.5], [0]),
+        ),
+        (
+            [0, 2, 0],
+            1.5,
+            1,
+            False,
+            32,
+            ([1], [1], [1], [0]),
+        ),
+        (
+            [0, 2, 2, 0, 4, 0],
+            [1, 3, 5],
+            1,
+            False,
+            32,
+            ([3, 1, 0], [2, 1, 0], [1.5, 1, 0], [0.5, 0, 0]),
+        ),
+        (
+            np.array([2, 2, 0, 4, 0]),
+            [1, 3, 5],
+            1,
+            False,
+            32,
+            ([3, 1, 0], [1, 1, 0], [1.5, 1, 0], [0.5, 0, 0]),
+        ),
+        (
+            [0, 2, 0, 4, 0] * u.eV,
+            [1, 3, 5],
+            1,
+            False,
+            32,
+            ([2, 1, 0], [2, 1, 0], [1, 1, 0], [0, 0, 0]),
+        ),
+        (
+            [0, 2, 0, 4, 0],
+            [1, 3, 5],
+            1,
+            True,
+            2,
+            (
+                [2, 1, 0],
+                [2, 1, 0],
+                [1, 1, 0],
+                [0, 0, 0],
+                [[0, 2], [0, 2], [0, 0]],
+                [[0.75, 1.25], [0.75, 1.25], [0, 0]],
+            ),
+        ),
+    ],
+)
+def test_ExcessStatistics(signal, thresholds, time_step, pdf, bins, expected):
+    """test excess_stat function"""
+    excess_stats = ExcessStatistics(signal, thresholds, time_step)
+    assert excess_stats.total_time_above_threshold == expected[0]
+    assert excess_stats.number_of_crossings == expected[1]
+    assert excess_stats.average_times == expected[2]
+    assert excess_stats.rms_times == expected[3]
+    if pdf:
+        hist, bin_centers = excess_stats.hist(bins)
+        assert np.allclose(hist, expected[4])
+        assert np.allclose(bin_centers, expected[5])
+
+
+@pytest.mark.parametrize(
+    "signal, thresholds, time_step, bins, exception",
+    [([1, 2], 1, -1, 32, ValueError), ([1, 2], 1, 1, 1.5, TypeError)],
+)
+def test_ExcessStatistics_exception(signal, thresholds, time_step, bins, exception):
+    """test whether exception is risen"""
+    with pytest.raises(exception):
+        tmp = ExcessStatistics(signal, thresholds, time_step)
+        tmp.hist(bins)