Feature/issue 177/output sib spt windows (#225)

* Initial changes to return SIB/SPT

* More small changes

* output sib/spt

Use the nonwear_utils package to fix timing for sib/spt
Update orchestrator +tests to output data

* Move cleanup to own function in analytics

* Renaming utils module

* Update documentation

* Update wristpy_tutorial.md

This line is missing from the tutorial

* Changes for review

Remove abstract sleepdetector class as we have no plans to extend to new sleep detection algorithms in the immediate future
Add new sleep parameters data class to store internal sleep parameters to be used in post-proc

* Remove tuple returns and run function twice

* Update test_analytics.py

No need for List in 3.10+

Author: Asanto32

README.md

@@ -36,7 +36,7 @@ The main processing pipeline of the wristpy module can be described as follows:
 - ***Data imputation*** In the special case when dealing with the Actigraph `idle_sleep_mode == enabled`, the gaps in acceleration are filled in after calibration, to avoid biasing the calibration phase.
 - **Metrics Calculation**: Calculates various activity metrics on the calibrated data, namely ENMO (Euclidean norm, minus one), MAD (mean amplitude deviation) <sup>1</sup>, Actigraph activity counts<sup>2</sup>, MIMS (monitor-independent movement summary) unit <sup>3</sup>, and angle-Z (angle of acceleration relative to the *x-y* axis).
 - **Non-wear detection**: We find periods of non-wear based on the acceleration data. Specifically, the standard deviation of the acceleration values in a given time window, along each axis, is used as a threshold to decide `wear` or `not wear`. Additionally, we can use the temperature sensor, when avaia\lable, to augment the acceleration data. This is used in the CTA (combined temperature and acceleration) algorithm <sup>4</sup>, and in the `skdh` DETACH algorithm <sup>5</sup>. Furthermore, ensemble classification of non-wear periods is possible by providing a list (of any length) of non-wear algorithm options.
-- **Sleep Detection**: Using the HDCZ<sup>6</sup> and HSPT<sup>7</sup> algorithms to analyze changes in arm angle we are able to find periods of sleep. We find the sleep onset-wakeup times for all sleep windows detected. Any sleep periods that overlap with detected non-wear times are removed, and any remaining sleep periods shorter than 15 minutes (default value) are removed.
+- **Sleep Detection**: Using the HDCZ<sup>6</sup> and HSPT<sup>7</sup> algorithms to analyze changes in arm angle we are able to find periods of sleep. We find the sleep onset-wakeup times for all sleep windows detected. Any sleep periods that overlap with detected non-wear times are removed, and any remaining sleep periods shorter than 15 minutes (default value) are removed. Additionally, the SIB (sustained inactivity bouts) and the SPT (sleep period time) windows are provided as part of the output to aid in sleep metric post-processing.
 - **Physical activity levels**: Using the chosen physical activity metric (aggregated into time bins, 5 second default) we compute activity levels into the following categories: [`inactive`, `light`, `moderate`, `vigorous`]. The threshold values can be defined by the user, while the default values are chosen based on the specific activity metric and the values found in the literature  <sup>8-10</sup>.
 - **Data output**: The output results can be saved in `.csv` or `.parquet` data formats, with the run-time configuration parameters saved in a `.json` dictionary.
 
@@ -102,8 +102,11 @@ results = orchestrator.run(
 physical_activity_metric = results.physical_activity_metric
 anglez = results.anglez
 physical_activity_levels = results.physical_activity_levels
-nonwear_array = results.nonwear_epoch
-sleep_windows = results.sleep_windows_epoch
+nonwear_array = results.nonwear_status
+sleep_windows = results.sleep_status
+sib_periods = results.sib_periods
+spt_periods = results.spt_periods
+
 ```
 #### Running entire directories:
 ```Python
@@ -131,8 +134,11 @@ subject1 = results_dict['subject1']
 physical_activity_metric = subject1.physical_activity_metric
 anglez = subject1.anglez
 physical_activity_levels = subject1.physical_activity_levels
-nonwear_array = subject1.nonwear_epoch
-sleep_windows = subject1.sleep_windows_epoch
+nonwear_array = subject1.nonwear_status
+sleep_windows = subject1.sleep_status
+sib_periods = subject1.sib_periods
+spt_periods = subject1.spt_periods
+
 ```
 
 ### Using Wristpy Through Docker

docs/wristpy_tutorial.md

@@ -96,6 +96,8 @@ We can also view and process these outputs from the saved `.csv` output file:
 import polars as pl
 import matplotlib.pyplot as plt
 
+output_results = pl.read_csv('path/to/save/file_name.csv', try_parse_dates=True)
+
 activity_mapping = {
     "inactive": 0,
     "light": 1,

src/wristpy/core/orchestrator.py

@@ -13,7 +13,7 @@
     calibration,
     idle_sleep_mode_imputation,
     metrics,
-    nonwear_utils,
+    processing_utils,
 )
 
 logger = config.get_logger()
@@ -292,15 +292,6 @@ def _run_file(
     if output is not None:
         writers.OrchestratorResults.validate_output(output=output)
 
-    parameters_dictionary = {
-        "thresholds": list(thresholds),
-        "calibrator": calibrator,
-        "epoch_length": epoch_length,
-        "activity_metric": activity_metric,
-        "nonwear_algorithm": list(nonwear_algorithm),
-        "input_file": str(input),
-    }
-
     if calibrator is not None and calibrator not in ["ggir", "gradient"]:
         msg = (
             f"Invalid calibrator: {calibrator}. Choose: 'ggir', 'gradient'. "
@@ -342,17 +333,14 @@ def _run_file(
         dynamic_range=watch_data.dynamic_range,
     )
 
-    sleep_detector = analytics.GgirSleepDetection(anglez)
-    sleep_windows = sleep_detector.run_sleep_detection()
-
-    nonwear_array = nonwear_utils.get_nonwear_measurements(
+    nonwear_array = processing_utils.get_nonwear_measurements(
         calibrated_acceleration=calibrated_acceleration,
         temperature=watch_data.temperature,
         non_wear_algorithms=nonwear_algorithm,
     )
 
-    nonwear_epoch = nonwear_utils.nonwear_array_cleanup(
-        nonwear_array=nonwear_array,
+    nonwear_epoch = processing_utils.synchronize_measurements(
+        data_measurement=nonwear_array,
         reference_measurement=activity_measurement,
         epoch_length=epoch_length,
     )
@@ -361,14 +349,40 @@ def _run_file(
         activity_measurement, thresholds
     )
 
+    sleep_detector = analytics.GgirSleepDetection(anglez)
+    sleep_parameters = sleep_detector.run_sleep_detection()
     sleep_array = analytics.sleep_cleanup(
-        sleep_windows=sleep_windows, nonwear_measurement=nonwear_epoch
+        sleep_windows=sleep_parameters.sleep_windows, nonwear_measurement=nonwear_epoch
+    )
+    spt_windows = analytics.sleep_bouts_cleanup(
+        sleep_parameter=sleep_parameters.spt_windows,
+        nonwear_measurement=nonwear_epoch,
+        time_reference_measurement=activity_measurement,
+        epoch_length=epoch_length,
+    )
+    sib_periods = analytics.sleep_bouts_cleanup(
+        sleep_parameter=sleep_parameters.sib_periods,
+        nonwear_measurement=nonwear_epoch,
+        time_reference_measurement=activity_measurement,
+        epoch_length=epoch_length,
     )
+
+    parameters_dictionary = {
+        "thresholds": list(thresholds),
+        "calibrator": calibrator,
+        "epoch_length": epoch_length,
+        "activity_metric": activity_metric,
+        "nonwear_algorithm": list(nonwear_algorithm),
+        "input_file": str(input),
+    }
+
     results = writers.OrchestratorResults(
         physical_activity_metric=activity_measurement,
         anglez=anglez,
         physical_activity_levels=physical_activity_levels,
         sleep_status=sleep_array,
+        sib_periods=sib_periods,
+        spt_periods=spt_windows,
         nonwear_status=nonwear_epoch,
         processing_params=parameters_dictionary,
     )

src/wristpy/io/writers/writers.py

@@ -23,6 +23,8 @@ class OrchestratorResults(pydantic.BaseModel):
     physical_activity_levels: models.Measurement
     nonwear_status: models.Measurement
     sleep_status: models.Measurement
+    sib_periods: models.Measurement
+    spt_periods: models.Measurement
     processing_params: Optional[Dict[str, Any]] = None
 
     def save_results(self, output: pathlib.Path) -> None:

src/wristpy/processing/analytics.py

@@ -1,6 +1,5 @@
 """Calculate sleep onset and wake up times."""
 
-import abc
 import datetime
 from dataclasses import dataclass
 from typing import List, Tuple, Union
@@ -9,6 +8,7 @@
 import polars as pl
 
 from wristpy.core import computations, config, models
+from wristpy.processing import processing_utils
 
 logger = config.get_logger()
 
@@ -26,27 +26,23 @@ class SleepWindow:
     wakeup: Union[datetime.datetime, List]
 
 
-class AbstractSleepDetector(abc.ABC):
-    """Abstract class defining the interface for sleep detection algorithms."""
-
-    @abc.abstractmethod
-    def __init__(self, anglez: models.Measurement) -> None:
-        """Initialization function for the sleep detection algorithm.
-
-        Must contain the anglez data as an input.
-        """
-        pass
+@dataclass
+class SleepParameters:
+    """Dataclass to store sleep parameters used to compute sleep metrics.
 
-    @abc.abstractmethod
-    def run_sleep_detection(self) -> List[SleepWindow]:
-        """Sleep Detector must contain a run_sleep_detection function.
+    Attributes:
+        sleep_windows: a list of SleepWindow objects, each representing a sleep period
+            with an onset and wakeup time.
+        spt_windows: a Measurement object with the sleep period guider windows.
+        sib_periods: a Measurement object with the sustained inactivity bouts.
+    """
 
-        The function must return a list of SleepWindow objects.
-        """
-        pass
+    sleep_windows: List[SleepWindow]
+    spt_windows: models.Measurement
+    sib_periods: models.Measurement
 
 
-class GgirSleepDetection(AbstractSleepDetector):
+class GgirSleepDetection:
     """Sleep Detection algorithm based on the GGIR method.
 
     This class implements the GGIR method for sleep detection. The method uses the
@@ -68,7 +64,9 @@ def __init__(
         """
         self.anglez = anglez
 
-    def run_sleep_detection(self) -> List[SleepWindow]:
+    def run_sleep_detection(
+        self,
+    ) -> SleepParameters:
         """Run the GGIR sleep detection.
 
         This algorithm uses the angle-z data to first find potential sleep periods
@@ -77,8 +75,10 @@ def run_sleep_detection(self) -> List[SleepWindow]:
         the SPT windows and SIB periods.
 
         Returns:
-            A list of SleepWindow instances, each instance contains a sleep onset/wakeup
-            time pair.
+            A SleepParameters instance containing the underlying sleep parameters:
+                - sleep_windows
+                - spt_window
+                - sib_periods
         """
         logger.debug("Beginning sleep detection.")
         spt_window = self._spt_window(self.anglez)
@@ -91,7 +91,11 @@ def run_sleep_detection(self) -> List[SleepWindow]:
         logger.debug(
             "Sleep detection complete. Windows detected: %s", len(sleep_onset_wakeup)
         )
-        return sleep_onset_wakeup
+        return SleepParameters(
+            sleep_windows=sleep_onset_wakeup,
+            spt_windows=spt_window,
+            sib_periods=sib_periods,
+        )
 
     def _spt_window(
         self, anglez_data: models.Measurement, threshold: float = 0.2
@@ -263,7 +267,7 @@ def _find_periods(
 
     This is a helper function to return the periods in the format of
     List [start_of_period, end_of_period], it is used in the
-    GGIRSleepDetection class and when removing non-wear periods from sleep windows.
+    GGIRSleepDetection class.
 
     Args:
         window_measurement: the Measurement instance, intended to be
@@ -424,6 +428,43 @@ def sleep_cleanup(
     return models.Measurement(time=sleep.time, measurements=cleaned_sleep)
 
 
+def sleep_bouts_cleanup(
+    sleep_parameter: models.Measurement,
+    nonwear_measurement: models.Measurement,
+    time_reference_measurement: models.Measurement,
+    epoch_length: float,
+) -> models.Measurement:
+    """This function will synchronize and filter the SPT and SIB windows.
+
+    The time sychrnoization is based on the time_reference_measurement, while the
+    filtering is based on the nonwear_measurement.
+
+    Args:
+        sleep_parameter: The sleep parameter measurement data, which contains
+            either the SPT and SIB windows.
+        nonwear_measurement: The nonwear measurement data used for reference time
+            stamps and to remove overlaps with periods of sleep.
+        time_reference_measurement: The time reference measurement data used for
+            reference time stamps.
+        epoch_length: The epoch length in seconds, used for resampling the data.
+
+    Returns:
+        A tuple of two Measurement instances with the cleaned SPT and SIB data.
+    """
+    logger.debug("Starting the sleep bouts cleanup.")
+    sleep_parameter_sync = processing_utils.synchronize_measurements(
+        data_measurement=sleep_parameter,
+        reference_measurement=time_reference_measurement,
+        epoch_length=epoch_length,
+    )
+    sleep_parameter_sync.measurements = np.logical_and(
+        sleep_parameter_sync.measurements,
+        np.logical_not(nonwear_measurement.measurements.astype(bool)),
+    )
+
+    return sleep_parameter_sync
+
+
 def _sleep_windows_as_measurement(
     ref_measurement_time: pl.Series, sleep_windows: List[SleepWindow]
 ) -> models.Measurement:

src/wristpy/processing/nonwear_utils.py …c/wristpy/processing/processing_utils.pysrc/wristpy/processing/nonwear_utils.py renamed to src/wristpy/processing/processing_utils.py src/wristpy/processing/nonwear_utils.py renamed to src/wristpy/processing/processing_utils.py

@@ -1,4 +1,4 @@
-"""This module contains helper functions for aggregated nonwear detection outputs."""
+"""This module contains helper functions for general processing of timeseries data."""
 
 import datetime
 from typing import Literal, Sequence, Union
@@ -160,22 +160,23 @@ def get_nonwear_measurements(
     return results[0]
 
 
-def nonwear_array_cleanup(
-    nonwear_array: models.Measurement,
+def synchronize_measurements(
+    data_measurement: models.Measurement,
     reference_measurement: models.Measurement,
     epoch_length: float = 5.0,
 ) -> models.Measurement:
-    """This function is used to match the nonwear array to a reference Measurement.
+    """This function is used to match a Measurement object to a reference Measurement.
 
-    This function ensures that the nonwear array and reference Measurement times
-    are synced up. This is accomplished by first resampling the nonwear array to
+    This function ensures that a Measurement object and reference Measurement times
+    are synced up. This is accomplished by first resampling a Measurement object to
     the specified temporal resolution.
-    It also ensures that the nonwear array is a binary array, where 1 indicates
+    It also ensures that a Measurement object is a binary array, where 1 indicates
     nonwear and 0 indicates wear.
-    It then truncates the nonwear array to match the reference Measurement time points.
+    It then truncates a Measurement object to match the reference
+    Measurement time points.
 
     Args:
-        nonwear_array: The nonwear array to clean up.
+        data_measurement: The Measurement object that requires time syncing.
         reference_measurement: The reference measurement to use for resampling.
         epoch_length: The temporal resolution of the output, in seconds.
             Defaults to 5.0.
@@ -185,7 +186,7 @@ def nonwear_array_cleanup(
         returned as a boolean (True == nonwear).
     """
     time_fix_nonwear = _time_fix(
-        nonwear_array, reference_measurement.time[-1], reference_measurement.time[0]
+        data_measurement, reference_measurement.time[-1], reference_measurement.time[0]
     )
     resampled_nonwear = computations.resample(time_fix_nonwear, epoch_length)
     binary_nonwear = np.where(resampled_nonwear.measurements >= 0.5, 1, 0)

tests/unit/test_analytics.py

@@ -2,7 +2,6 @@
 
 import datetime
 import math
-from typing import List
 
 import numpy as np
 import polars as pl
@@ -160,8 +159,10 @@ def test_run_sleep_detection(sleep_detection: analytics.GgirSleepDetection) -> N
     """Test the full sleep detection process."""
     result = sleep_detection.run_sleep_detection()
 
-    assert result == []
-    assert isinstance(result, List)
+    assert result.sleep_windows == []
+    assert isinstance(result.sleep_windows, list)
+    assert isinstance(result.spt_windows, models.Measurement)
+    assert isinstance(result.sib_periods, models.Measurement)
 
 
 def test_physical_activity_thresholds() -> None:

tests/unit/test_orchestrator.py

@@ -29,6 +29,8 @@ def dummy_results() -> writers.OrchestratorResults:
         physical_activity_levels=dummy_measure,
         nonwear_status=dummy_measure,
         sleep_status=dummy_measure,
+        sib_periods=dummy_measure,
+        spt_periods=dummy_measure,
     )
 
     return dummy_results