Source code for pyActigraphy.sleep.scoring_base

import numpy as np
import pandas as pd
import re
from .scoring import csm, roenneberg, sleep_midpoint, sri
from .scoring.utils import rescore
from scipy.ndimage import binary_closing, binary_opening
from ..filters import _create_inactivity_mask
from ..utils.utils import _average_daily_activity
from ..utils.utils import _activity_onset_time
from ..utils.utils import _activity_offset_time


def _td_format(td):
    return '{:02}:{:02}:{:02}'.format(
        td.components.hours,
        td.components.minutes,
        td.components.seconds
    )


def _actiware_automatic_threshold(data, scale_factor=0.88888):
    r'''Automatic Wake Threshold Value calculation


    1. Sum the activity counts for all epochs of the data set.
    2. Count the number of epochs scored as MOBILE for the data set
    (the definition of MOBILE follows).
    3. Compute the MOBILE TIME (number of epochs scored as MOBILE from step 2
    multiplied by the Epoch Length) in minutes.
    4. Compute the Auto Wake Threshold = ((sum of activity counts from step 1)
    divided by (MOBILE TIME from step 3) ) multiplied by 0.88888.

    Definition of Mobile:
    An epoch is scored as MOBILE if the number of activity counts recorded in
    that epoch is greater than or equal to the epoch length in 15-second
    intervals. For example,there are four 15-second intervals for a 1-minute
    epoch length; hence, the activity value in an epoch must be greater than,
    or equal to, four, to be scored as MOBILE.
    '''

    # Sum of activity counts
    counts_sum = data.sum()

    # Definition of the "mobile" threshold
    mobile_thr = int(data.index.freq/pd.Timedelta('15sec'))

    # Counts the number of epochs scored as "mobile"
    counts_mobile = (data.values >= mobile_thr).astype(int).sum()

    # Mobile time
    mobile_time = counts_mobile*(data.index.freq/pd.Timedelta('1min'))

    # Automatic wake threshold
    automatic_thr = (counts_sum/mobile_time)*scale_factor

    return automatic_thr


def _padded_data(data, value, periods, frequency):

    date_offset = pd.DateOffset(seconds=frequency.total_seconds())
    pad_beginning = pd.Series(
        data=value,
        index=pd.date_range(
            end=data.index[0],
            periods=periods,
            freq=date_offset,
            closed='left'
        ),
        dtype=data.dtype
    )
    pad_end = pd.Series(
        data=value,
        index=pd.date_range(
            start=data.index[-1],
            periods=periods,
            freq=date_offset,
            closed='right'
        ),
        dtype=data.dtype
    )
    return pd.concat([pad_beginning, data, pad_end])


def _ratio_sequences_of_zeroes(
    data, seq_length, n_boostrap, seed=0, with_replacement=True
):
    # Calculate a moving sum with a window of size 'seq_length'
    rolling_sum = data.rolling(seq_length).sum()
    # Set seed for reproducibility
    np.random.seed(seed)
    random_sample = np.random.choice(
        rolling_sum.values,
        size=n_boostrap*len(rolling_sum),
        replace=with_replacement
    )
    # Calculate the ratio of zero elements
    ratio = 1 - np.count_nonzero(random_sample)/len(random_sample)
    return ratio


def _estimate_zeta(data, seq_length_max, n_boostrap=100, level=0.05):
    ratios = np.fromiter((
        _ratio_sequences_of_zeroes(data, n, n_boostrap) for n in np.arange(
            1, seq_length_max+1
        )),
        np.float,
        seq_length_max
    )
    zeta_est = np.argmax(ratios < level)
    return zeta_est


def _window_convolution(x, scale, window, offset=0.0):

    return scale * np.dot(x, window) + offset


def _cole_kripke(data, scale, window, threshold):
    """Automatic scoring methods from Cole and Kripke"""

    ck = data.rolling(
        window.size, center=True
    ).apply(_window_convolution, args=(scale, window), raw=True)

    return (ck < threshold).astype(int)


def _sadeh(data, offset, weights, threshold):
    """Activity-Based Sleep-Wake Identification"""

    r = data.rolling(11, center=True)

    mean_W5 = r.mean()

    NAT = r.apply(lambda x: np.size(np.where((x > 50) & (x < 100))), raw=True)

    sd_Last6 = data.rolling(6).std()

    logAct = data.shift(-1).apply(lambda x: np.log(1+x))

    sadeh = pd.concat(
        [mean_W5, NAT, sd_Last6, logAct],
        axis=1,
        keys=['mean_W5', 'NAT', 'sd_Last6', 'logAct']
    )

    sadeh['PS'] = sadeh.apply(
        _window_convolution, axis=1, args=(1.0, weights, offset), raw=True
    )

    return (sadeh['PS'] > threshold).astype(int)


def _scripps(data, scale, window, threshold):
    """Scripps Clinic algorithm for sleep-wake scoring."""

    scripps = data.rolling(
        window.size, center=True
    ).apply(_window_convolution, args=(scale, window), raw=True)

    return (scripps < threshold).astype(int)


def _oakley(data, window, threshold):
    """Oakley's algorithm for sleep-wake scoring."""
    if threshold == 'automatic':
        threshold = _actiware_automatic_threshold(data)
    elif not np.isscalar(threshold):
        msg = "`threshold` should be a scalar or 'automatic'."
        raise ValueError(msg)

    scale = 1.
    oakley = data.rolling(
        window.size, center=True
    ).apply(_window_convolution, args=(scale, window), raw=True)

    return (oakley < threshold).astype(int)


[docs]class ScoringMixin(object): """ Mixin Class for scoring: - sleep/wake periods - ... """
[docs] def AonT(self, freq='5min', whs=12, binarize=True, threshold=4): r"""Activity onset time. Activity onset time derived from the daily activity profile. Parameters ---------- freq: str, optional Data resampling `frequency string <https://pandas.pydata.org/pandas-docs/stable/timeseries.html>`_. Default is '5min'. whs: int, optional Window half size. Default is 12. binarize: bool, optional If set to True, the data are binarized. Default is True. threshold: int, optional If binarize is set to True, data above this threshold are set to 1 and to 0 otherwise. Default is 4. Returns ------- aot: Timedelta Activity onset time. Notes ----- This automatic detection of the activity onset time is based on the daily activity profile. It returns the time point where difference between the mean activity over :math:`whs` epochs before and after this time point is maximum: .. math:: AonT = \max_{t}( \frac{\sum_{i=-whs}^{0} x_{t+i}}{\sum_{i=1}^{whs} x_{t+i}} - 1 ) with: * :math:`x_{i}` is the activity count at time :math:`i`. Examples -------- >>> import pyActigraphy >>> rawAWD = pyActigraphy.io.read_raw_awd(fpath + 'SUBJECT_01.AWD') >>> raw.AonT() Timedelta('0 days 07:15:00') >>> raw.AonT(binarize=False) Timedelta('0 days 07:05:00') """ data = self.resampled_data(freq, binarize, threshold) dailyprof = _average_daily_activity(data, cyclic=False) return _activity_onset_time(dailyprof, whs=whs)
[docs] def AoffT(self, freq='5min', whs=12, binarize=True, threshold=4): r"""Activity offset time. Activity offset time derived from the daily activity profile. Parameters ---------- freq: str, optional Data resampling `frequency string <https://pandas.pydata.org/pandas-docs/stable/timeseries.html>`_. Default is '5min'. whs: int, optional Window half size. Default is 12. binarize: bool, optional If set to True, the data are binarized. Default is True. threshold: int, optional If binarize is set to True, data above this threshold are set to 1 and to 0 otherwise. Default is 4. Returns ------- aot: Timedelta Activity offset time. Notes ----- This automatic detection of the activity offset time is based on the daily activity profile. It returns the time point where relative difference between the mean activity over :math:`whs` epochs after and before this time point is maximum: .. math:: AoffT = \max_{t}( \frac{\sum_{i=1}^{whs} x_{t+i}}{\sum_{i=-whs}^{0} x_{t+i}} - 1 ) with: * :math:`x_{i}` is the activity count at time :math:`i`. Examples -------- >>> import pyActigraphy >>> rawAWD = pyActigraphy.io.read_raw_awd(fpath + 'SUBJECT_01.AWD') >>> raw.AoffT() Timedelta('0 days 23:20:00') >>> raw.AoffT(binarize=False) Timedelta('0 days 23:05:00') """ data = self.resampled_data(freq, binarize, threshold) dailyprof = _average_daily_activity(data, cyclic=False) return _activity_offset_time(dailyprof, whs=whs)
[docs] def CK( self, settings='30sec_max_non_overlap', threshold=1.0, rescoring=True ): r"""Cole&Kripke algorithm for sleep-wake identification. Algorithm for automatic sleep scoring based on wrist activity, developped by Cole, Kripke et al [1]_. Parameters ---------- settings: str, optional Data reduction settings for which the optimal parameters have been derived. Available settings are: * "mean": mean activity per minute * "10sec_max_overlap": maximum 10-second overlapping epoch per minute * "10sec_max_non_overlap": maximum 10-second nonoverlapping epoch per minute * "30sec_max_non_overlap": maximum 30-second nonoverlapping epoch per minute Default is "30sec_max_non_overlap". threshold: float, optional Threshold value for scoring sleep/wake. Default is 1.0. rescoring: bool, optional If set to True, Webster's rescoring rules are applied [2]_. Default is True. Returns ------- ck: pandas.core.Series Time series containing the `D` scores (0: sleep, 1: wake) for each epoch. Notes ----- The output variable D of the CK algorithm is defined in [1]_ as: .. math:: D = P*( [W_{-4},\dots,W_{0},\dots,W_{+2}] \cdot [A_{-4},\dots,A_{0},\dots,A_{+2}]) with: * D < 1 == sleep, D >= 1 == wake; * P, scale factor; * :math:`W_{0},W_{-1},W_{+1},\dots`, weighting factors for the present minute, the previous minute, the following minute, etc.; * :math:`A_{0},A_{-1},A_{+1},\dots`, activity scores for the present minute, the previous minute, the following minute, etc. .. warning:: This algorithm yields scores with a period of 1 minute. Therefore, the time series returned by the CK function has a 1-min period. In the original paper, this algorithm is validated for devices configured in "zero-crossing mode". Webster's rescoring rules [2]_: 1. After at least 4 minutes scored as wake, the next 1 minute scored as sleep is rescored as wake 2. After at least 10 minutes scored as wake, the next 3 minutes scored as sleep are rescored as wake; 3. After at least 15 minutes scored as wake, the next 4 minutes scored as sleep are rescored as wake; 4. If a period of 6 minutes or less that is scored as sleep is surrounded by at least 15 minutes scored as wake, then rescore to wake; 5. If a period of 10 minutes or less that is scored as sleep is surrounded by at least 20 minutes scored as wake, then rescore to wake. References ---------- .. [1] Cole, R. J., Kripke, D. F., Gruen, W., Mullaney, D. J., & Gillin, J. C. (1992). Automatic Sleep/Wake Identification From Wrist Activity. Sleep, 15(5), 461–469. http://doi.org/10.1093/sleep/15.5.461 .. [2] Webster, J. B., Kripke, D. F., Messin, S., Mullaney, D. J., & Wyborney, G. (1982). An Activity-Based Sleep Monitor System for Ambulatory Use. Sleep, 5(4), 389–399. https://doi.org/10.1093/sleep/5.4.389 Examples -------- """ available_settings = [ "mean", "10sec_max_overlap", "10sec_max_non_overlap", "30sec_max_non_overlap" ] if settings not in available_settings: raise ValueError("CK sleep/wake identification:\n" + "Required settings: {}\n".format(settings) + "Available settings are:\n" + "\n".join(available_settings)) ck = None if settings == "mean": if pd.Timedelta(self.data.index.freq) > pd.Timedelta('60s'): raise ValueError(( "The sampling frequency of the input data does not allow" + " the use of the requested settings ({}).\n".format( settings) + "For such settings, the sampling frequency should less" + " than 60 seconds." )) else: # Compute the resampling factor as the original weights are # calculated for an average of 2-sec epochs over 60 seconds: rs_f = 30 # 60sec/2sec # Resample to 60 sec and take the sum # WARNING: valid if the input data is a sum of activity counts data_mean = self.data.resample('60s').sum()/rs_f # Define the scale and weights for this settings scale = 0.001 window = np.array( [106, 54, 58, 76, 230, 74, 67, 0, 0], np.int32 ) ck = _cole_kripke(data_mean, scale, window, threshold) elif settings == "10sec_max_overlap": if pd.Timedelta(self.data.index.freq) > pd.Timedelta('5s'): raise ValueError(( "The sampling frequency of the input data does not allow" + " the use of the requested settings ({}).\n".format( settings) + "For such settings, the sampling frequency should less" + " than 5 seconds." )) else: # Resample to 10/2 sec and then sum over a sliding window data = self.data.resample('5s').sum().rolling('10s').sum() # Resample to 60 sec and take the max of the 10 sec periods data_max = data.resample('60s').max() # Define the scale and weights for this settings scale = 0.00001 window = np.array( [50, 30, 300, 400, 1400, 500, 350, 0, 0], np.int32 ) ck = _cole_kripke(data_max, scale, window, threshold) elif settings == "10sec_max_non_overlap": if pd.Timedelta(self.data.index.freq) > pd.Timedelta('10s'): raise ValueError(( "The sampling frequency of the input data does not allow" + " the use of the requested settings ({}).\n".format( settings) + "For such settings, the sampling frequency should less" + " or equal to 10 seconds." )) else: # Resample to 10 sec and then sum data = self.data.resample('10s').sum() # Resample to 60 sec and take the max of the 10 sec periods data_max = data.resample('60s').max() # Define the scale and weights for this settings scale = 0.00001 window = np.array( [550, 378, 413, 699, 1736, 287, 300, 0, 0], np.int32 ) ck = _cole_kripke(data_max, scale, window, threshold) elif settings == "30sec_max_non_overlap": if pd.Timedelta(self.data.index.freq) > pd.Timedelta('30s'): raise ValueError(( "The sampling frequency of the input data does not allow" + " the use of the requested settings ({}).\n".format( settings) + "For such settings, the sampling frequency should less" + " or equal to 30 seconds." )) else: # Resample to 30 sec and then sum data = self.data.resample('30s').sum() # Resample to 60 sec and take the max of the 10 sec periods data_max = data.resample('60s').max() # Define the scale and weights for this settings scale = 0.0001 window = np.array( [50, 30, 14, 28, 12, 8, 50, 0, 0], np.int32 ) ck = _cole_kripke(data_max, scale, window, threshold) if rescoring: # Rescoring mask = rescore(ck.values, sleep_score=1) ck = ck*mask return ck
[docs] def Sadeh( self, offset=7.601, weights=np.array([-0.065, -1.08, -0.056, -0.703], float), threshold=0.0 ): r"""Sadeh algorithm for sleep identification Algorithm for automatic sleep scoring based on wrist activity, developped by Sadeh et al [1]_. Parameters ---------- offset: float, optional Offset parameter. Default is 7.601. weights: np.array Array of weighting factors for mean_W5, NAT, sd_Last6 and logAct. Default is [-0.065, -1.08, -0.056, -0.703]. threshold: float, optional Threshold value for scoring sleep/wake. Default is 0.0. Returns ------- Notes ----- The output variable PS of the Sadeh algorithm is defined in [2]_ as: .. math:: PS = 7.601-0.065·mean_W5-1.08·NAT-0.056·sd_Last6-0.703·logAct with: * PS >= 0 == sleep, PS < 0 == wake; * mean_W5, the average number of activity counts during the scored epoch and the window of five epochs preceding and following it; * sd_Last6, the standard deviation of the activity counts during the scored epoch and the five epochs preceding it; * NAT, the number of epochs with activity level equal to or higher than 50 but lower than 100 activity counts in a window of 11 minutes that includes the scored epoch and the five epochs preceding and following it; * logAct, the natural logarithm of the number of activity counts during the scored epoch + 1. References ---------- .. [1] Sadeh, A., Alster, J., Urbach, D., & Lavie, P. (1989). Actigraphically based automatic bedtime sleep-wake scoring: validity and clinical applications. Journal of ambulatory monitoring, 2(3), 209-216. .. [2] Sadeh, A., Sharkey, M., & Carskadon, M. A. (1994). Activity-Based Sleep-Wake Identification: An Empirical Test of Methodological Issues. Sleep, 17(3), 201–207. http://doi.org/10.1093/sleep/17.3.201 Examples -------- """ return _sadeh(self.data, offset, weights, threshold)
[docs] def Scripps( self, scale=0.204, window=np.array([ 0.0064, # b_{-10} 0.0074, # b_{-9} 0.0112, # b_{-8} 0.0112, # b_{-7} 0.0118, # b_{-6} 0.0118, # b_{-5} 0.0128, # b_{-4} 0.0188, # b_{-3} 0.0280, # b_{-2} 0.0664, # b_{-1} 0.0300, # b_{+0} 0.0112, # b_{+1} 0.0100, # b_{+2} 0.0000, # b_{+3} 0.0000, # b_{+4} 0.0000, # b_{+5} 0.0000, # b_{+6} 0.0000, # b_{+7} 0.0000, # b_{+8} 0.0000, # b_{+9} 0.0000 # b_{+10} ], float), threshold=1.0 ): r"""Scripps Clinic algorithm for sleep-wake identification. Algorithm for automatic sleep scoring based on wrist activity, developed by Kripke et al [1]_. Parameters ---------- scale: float, optional Scale parameter P Default is 0.204. window: np.array, optional Array of weighting factors :math:`W_{i}` Default values are identical to those found in the original publication [1]_. threshold: float, optional Threshold value for scoring sleep/wake. Default is 1.0. Returns ------- scripps: pandas.core.Series Time series containing the `D` scores (0: sleep, 1: wake) for each epoch. Notes ----- The output variable D of the Scripps algorithm is defined in [1]_ as: .. math:: D = P*( [W_{-10},\dots,W_{0},\dots,W_{+10}] \cdot [A_{-10},\dots,A_{0},\dots,A_{+10}]) with: * D < 1 == sleep, D >= 1 == wake; * P, scale factor; * :math:`W_{0},W_{-1},W_{+1},\dots`, weighting factors for the present epoch, the previous epoch, the following epoch, etc.; * :math:`A_{0},A_{-1},A_{+1},\dots`, activity scores for the present epoch, the previous epoch, the following epoch, etc. References ---------- .. [1] Kripke, D. F., Hahn, E. K., Grizas, A. P., Wadiak, K. H., Loving, R. T., Poceta, J. S., … Kline, L. E. (2010). Wrist actigraphic scoring for sleep laboratory patients: algorithm development. Journal of Sleep Research, 19(4), 612–619. http://doi.org/10.1111/j.1365-2869.2010.00835.x """ return _scripps(self.data, scale, window, threshold)
[docs] def Oakley( self, threshold=40 ): r"""Oakley's algorithm for sleep/wake scoring. Algorithm for automatic sleep/wake scoring based on wrist activity, developed by Oakley [1]_. Parameters ---------- threshold: float or str, optional Threshold value for scoring sleep/wake. Can be set to "automatic" (cf. Notes). Default is 40. Returns ------- oakley: pandas.core.Series Time series containing scores (1: sleep, 0: wake) for each epoch. Notes ----- The output variable O of Oakley's algorithm is defined as: .. math:: O = ( [W_{-2},W_{-1},W_{0},W_{+1},W_{+2}] \cdot [A_{-2},A_{-1},A_{0},A_{+1},A_{+2}]) with: * O <= threshold == sleep, 0 > threshold == wake; * :math:`W_{0},W_{-1},W_{+1},\dots`, weighting factors for the present epoch, the previous epoch, the following epoch, etc.; * :math:`A_{0},A_{-1},A_{+1},\dots`, activity scores for the present epoch, the previous epoch, the following epoch, etc. The current implementation of this algorithm follows the description provided in the instruction manual of the Actiwatch Communication and Sleep Analysis Software (Respironics, Inc.) [2]_ : * 15-second sampling frequency: .. math:: W_{-8} &= \ldots = W_{-5} = W_{+5} = \ldots = W_{+8} = 1/25 \\ W_{-4} &= \ldots = W_{-1} = W_{+1} = \ldots = W_{+4} = 1/5 \\ W_{0} &= 4 * 30-second sampling frequency: .. math:: W_{-4} &= W_{-3} = W_{+3} = W_{+4} = 1/25 \\ W_{-2} &= W_{-1} = W_{+1} = W_{+2} = 1/5 \\ W_{0} &= 2 * 60-second sampling frequency: .. math:: W_{-2} &= W_{+2} = 1/25 \\ W_{-1} &= W_{+1} = 1/5 \\ W_{0} &= 1 * 120-second sampling frequency: .. math:: W_{-1} &= W_{+1} = 1/8 \\ W_{0} &= 1/2 The *Automatic Wake Threshold Value* calculation is this [2]_: 1. Sum the activity counts for all epochs of the data set. 2. Count the number of epochs scored as MOBILE for the data set (the definition of MOBILE follows). 3. Compute the MOBILE TIME (number of epochs scored as MOBILE from step 2 multiplied by the Epoch Length) in minutes. 4. Compute the Auto Wake Threshold = ((sum of activity counts from step 1) divided by (MOBILE TIME from step 3)) multiplied by 0.88888. *Definition of Mobile* [2]_: An epoch is scored as MOBILE if the number of activity counts recorded in that epoch is greater than or equal to the epoch length in 15-second intervals. For example,there are four 15-second intervals for a 1-minute epoch length; hence, the activity value in an epoch must be greater than, or equal to, four, to be scored as MOBILE. References ---------- .. [1] Oakley, N.R. Validation with Polysomnography of the Sleepwatch Sleep/Wake Scoring Algorithm Used by the Actiwatch Activity Monitoring System; Technical Report; Mini-Mitter: Bend, OR, USA, 1997 .. [2] Instruction manual, Actiwatch Communication and Sleep Analysis Software (https://fccid.io/JIAAWR1/Users-Manual/USERS-MANUAL-1-920937) """ # Sampling frequency freq = self.data.index.freq.delta if freq == pd.Timedelta('15S'): window = np.array([ 0.04, # W_{-8} 0.04, # W_{-7} 0.04, # W_{-6} 0.04, # W_{-5} 0.20, # W_{-4} 0.20, # W_{-3} 0.20, # W_{-2} 0.20, # W_{-1} 4.00, # W_{+0} 0.20, # W_{+1} 0.20, # W_{+2} 0.20, # W_{+3} 0.20, # W_{+4} 0.04, # W_{+5} 0.04, # W_{+6} 0.04, # W_{+7} 0.04 # W_{+8} ], np.float) elif freq == pd.Timedelta('30S'): window = np.array([ 0.04, # W_{-4} 0.04, # W_{-3} 0.20, # W_{-2} 0.20, # W_{-1} 2.00, # W_{+0} 0.20, # W_{+1} 0.20, # W_{+2} 0.04, # W_{+3} 0.04 # W_{+4} ], np.float) elif freq == pd.Timedelta('60S'): window = np.array([ 0.04, # W_{-2} 0.20, # W_{-1} 1.00, # W_{+0} 0.20, # W_{+1} 0.04 # W_{+2} ], np.float) elif freq == pd.Timedelta('120S'): window = np.array([ 0.12, # W_{-1} 0.50, # W_{+0} 0.12, # W_{+1} ], np.float) else: raise ValueError( 'Oakley\'s algorithm is not defined for data ' + 'acquired with a sampling frequency of {}. '.format(freq) + 'Accepted frequencies are: {}'.format( ', '.join(['15sec', '30sec', '60sec', '120sec']) ) ) return _oakley(self.data, window, threshold)
def CSM( self, settings="auto", score_rest=2, score_sleep=1, binarize=False ): """Condor Sleep Model Sleep-wake scoring algorithm developed by Condor Instrument for their ActTrust devices. This algorithm works in a two-step fashion. First, it classifies all epochs as wake or rest, as function of each epoch's score. Second, using a more stringent scoring threshold, "rest" epoch are re-classified as "sleep". A relabelling mechanism using the labels of the surrounding epochs is also applied to consolidate periods of epochs labelled as rest or sleep. Parameters ---------- settings: str, optional Parameter settings for the CSM algorithm. Refers to the data acquisition frequency. Available values are: * "auto": use input data frequency. * "30s": set parameters to optimal values obtained for a 30s data acquisition frequency. * "60s": set parameters to optimal values obtained for a 60s data acquisition frequency. Default is 'auto'. score_rest: int, optional State index for epochs labelled as "rest". Default is 2. score_sleep: int, optional State index for epochs labelled as "sleep". Default is 1. binarize: bool, optional. If set to True, the state index is set to 1 if the epoch is labelled as sleep and 0 otherwise. Defautl is False. Returns ------- csm : pandas.Series Series of state indices. """ # This algorithm has been developed for ActTrust devices from # Condor Instrument. Verify if the reader has the appropriate type: if self.format != 'ATR': raise ValueError( "The CSM has been developed for ActTrust devices.\n" "It has not been validated for other devices." ) # The CSM uses the ZCMn as input data = self.ZCMn if settings == "auto": freq = data.index.freq.delta else: freq = pd.Timedelta(settings) if freq == pd.Timedelta('60s'): wa = np.array( [34.5, 133, 529, 375, 408, 400.5, 1074, 2048.5, 2424.5] ) wp = np.array( [1920, 149.5, 257.5, 125, 111.5, 120, 69, 40.5] ) elif freq == pd.Timedelta('30s'): wa = np.array( [69, 197, 730, 328, 269, 481, 528, 288, 304, 497, 1105, 1043, 1378, 2719, 2852, 1997] ) wp = np.array( [2972, 868, 269, 30, 495, 20, 39, 211, 91, 132, 203, 37, 67, 71, 81] ) else: raise NotImplementedError( "The settings for this acquistion frequency ({}) ".format( freq ) + "have not been implemented yet." ) # The overall scaling and rescoring rules are identical for the # different settings. p_rest = 0.00005 p_sleep = 0.000464 pr_rest = 0 nr_rest = 0 pr_sleep = 1 nr_sleep = 0 states = csm( data, wa=wa, wp=wp, p_rest=p_rest, p_sleep=p_sleep, pr_rest=pr_rest, nr_rest=nr_rest, pr_sleep=pr_sleep, nr_sleep=nr_sleep, score_rest=score_rest, score_sleep=score_sleep ) return (states == score_sleep).astype(int) if binarize else states
[docs] def SoD( self, freq='5min', binarize=True, bin_threshold=4, whs=4, start='12:00:00', period='5h', algo='Roenneberg', *args, **kwargs ): r"""Sleep over Daytime Quantify the volume of epochs identified as sleep over daytime (SoD), using sleep-wake scoring algorithms. Parameters ---------- freq: str, optional Resampling frequency. Default is '5min' binarize: bool, optional If set to True, the data are binarized when determining the activity onset and offset times. Only valid if start='AonT' or 'AoffT'. Default is True. bin_threshold: int, optional If binarize is set to True, data above this threshold are set to 1 and to 0 otherwise. Default is 4. whs: int, optional Window half size. Only valid if start='AonT' or 'AoffT'. Default is 4 start: str, optional Start time of the period of interest. Default: '12:00:00' Supported times: 'AonT', 'AoffT', any 'HH:MM:SS' period: str, optional Period length. Default is '5h' algo: str, optional Sleep scoring algorithm to use. Default is 'Roenneberg'. *args Variable length argument list passed to the scoring algorithm. **kwargs Arbitrary keyword arguements passed to the scoring algorithm. Returns ------- sod: pandas.core.Series Time series containing the epochs of rest (1) and activity (0) over the specified period. Examples -------- >>> import pyActigraphy >>> rawAWD = pyActigraphy.io.read_raw_awd(fpath + 'SUBJECT_01.AWD') >>> SoD = rawAWD.SoD() >>> SoD 2018-03-26 04:16:00 1 2018-03-26 04:17:00 1 2018-03-26 04:18:00 1 (...) 2018-04-05 16:59:00 0 2018-04-05 16:59:00 0 2018-04-05 17:00:00 0 Length: 3175, dtype: int64 """ # Retrieve sleep scoring function dynamically by name sleep_algo = getattr(self, algo) # Score activity sleep_scoring = sleep_algo(*args, **kwargs) # Regex pattern for HH:MM:SS time string pattern = re.compile(r"^([0-1]\d|2[0-3])(?::([0-5]\d))(?::([0-5]\d))$") if start == 'AonT': td = self.AonT( freq=freq, whs=whs, binarize=binarize, threshold=bin_threshold ) elif start == 'AoffT': td = self.AoffT( freq=freq, whs=whs, binarize=binarize, threshold=bin_threshold ) elif pattern.match(start): td = pd.Timedelta(start) else: print('Input string for start ({}) not supported.'.format(start)) return None start_time = _td_format(td) end_time = _td_format(td+pd.Timedelta(period)) sod = sleep_scoring.between_time(start_time, end_time) return sod
[docs] def fSoD( self, freq='5min', binarize=True, bin_threshold=4, whs=12, start='12:00:00', period='5h', algo='Roenneberg', *args, **kwargs ): r"""Fraction of Sleep over Daytime Fractional volume of epochs identified as sleep over daytime (SoD), using sleep-wake scoring algorithms. Parameters ---------- freq: str, optional Resampling frequency. Default is '5min' binarize: bool, optional If set to True, the data are binarized when determining the activity onset and offset times. Only valid if start='AonT' or 'AoffT'. Default is True. bin_threshold: int, optional If binarize is set to True, data above this threshold are set to 1 and to 0 otherwise. Default is 4. whs: int, optional Window half size. Default is 4 start: str, optional Start time of the period of interest. Supported times: 'AonT', 'AoffT', any 'HH:MM:SS'. Default: '12:00:00'. period: str, optional Period length. Default is '10h'. algo: str, optional Sleep scoring algorithm to use. Default is 'Roenneberg'. args Variable length argument list passed to the scoring algorithm. kwargs Arbitrary keyword arguements passed to the scoring algorithm. Returns ------- fsod: float Fraction of epochs scored as sleep, relatively to the length of the specified period. .. warning:: The value of this variable depends on the convention used by the underlying sleep scoring algorithm. The expected convention is the following: * epochs scored as 1 refer to inactivity/sleep Otherwise, this variable will actually return the fraction of epochs scored as activity. The fraction of sleep can simply be recovered by calculating (1-fSOD). Examples -------- >>> import pyActigraphy >>> rawAWD = pyActigraphy.io.read_raw_awd(fpath + 'SUBJECT_01.AWD') >>> raw.fSoD() 0.17763779527559054 >>> raw.fSoD(algo='ck') 0.23811023622047245 """ SoD = self.SoD( freq=freq, binarize=binarize, bin_threshold=bin_threshold, whs=whs, start=start, period=period, algo=algo, *args, **kwargs ) return SoD.sum()/len(SoD)
def SleepFragmentation(self): """Sleep Fragmentation is an index of restlessness during the sleep period expressed as a percentage. The higher the index, the more sleep is disrupted. ActiLife calculates three values for sleep fragmentation: movement Index, a fragmentation index, and total sleep fragmentation Index. - The Movement Index (MI) is the percentage of epochs with y-axis counts greater than zero in the sleep period. - The Fragmentation Index (FI) is the percentage of one minute periods of sleep vs. all periods of sleep during the sleep period. - The Total Sleep Fragmentation Index (SFI) is the sum of the MI and the FI""" pass
[docs] def Crespo( self, zeta=15, zeta_r=30, zeta_a=2, t=.33, alpha='8h', beta='1h', estimate_zeta=False, seq_length_max=100, verbose=False ): r"""Crespo algorithm for activity/rest identification Algorithm for automatic identification of activity-rest periods based on actigraphy, developped by Crespo et al. [1]_. Parameters ---------- zeta: int, optional Maximum number of consecutive zeroes considered valid. Default is 15. zeta_r: int, optional Maximum number of consecutive zeroes considered valid (rest). Default is 30. zeta_a: int, optional Maximum number of consecutive zeroes considered valid (active). Default is 2. t: float, optional Percentile for invalid zeroes. Default is 0.33. alpha: str, optional Average hours of sleep per night. Default is '8h'. beta: str, optional Length of the padding sequence used during the processing. Default is '1h'. estimate_zeta: bool, optional If set to True, zeta values are estimated from the distribution of ratios of the number of series of consecutive zeroes to the number of series randomly chosen from the actigraphy data. Default is False. seq_length_max: int, optional Maximal length of the aforementioned random series. Default is 100. verbose: bool, optional If set to True, print the estimated values of zeta. Default is False. Returns ------- crespo : pandas.core.Series Time series containing the estimated periods of rest (0) and activity (1). References ---------- .. [1] Crespo, C., Aboy, M., Fernández, J. R., & Mojón, A. (2012). Automatic identification of activity–rest periods based on actigraphy. Medical & Biological Engineering & Computing, 50(4), 329–340. http://doi.org/10.1007/s11517-012-0875-y Examples -------- >>> import pyActigraphy >>> rawAWD = pyActigraphy.io.read_raw_awd(fpath + 'SUBJECT_01.AWD') >>> crespo = rawAWD.Crespo() >>> crespo 2018-03-26 14:16:00 1 2018-03-26 14:17:00 0 2018-03-26 14:18:00 0 (...) 2018-04-06 08:22:00 0 2018-04-06 08:23:00 0 2018-04-06 08:24:00 1 Length: 15489, dtype: int64 """ # 1. Pre-processing # This stage produces an initial estimate of the rest-activity periods # 1.1. Signal conditioning based on empirical probability model # This step replaces sequences of more than $\zeta$ "zeroes" # with the t-percentile value of the actigraphy data # zeta = 15 if estimate_zeta: zeta = _estimate_zeta(self.raw_data, seq_length_max) if verbose: print("CRESPO: estimated zeta = {}".format(zeta)) # Determine the sequences of consecutive zeroes mask_zeta = _create_inactivity_mask(self.raw_data, zeta, 1) # Determine the user-specified t-percentile value s_t = self.raw_data.quantile(t) # Replace zeroes with the t-percentile value x = self.raw_data.copy() x[mask_zeta > 0] = s_t # Median filter window length L_w L_w = int(pd.Timedelta(alpha)/self.frequency)+1 L_w_over_2 = int((L_w-1)/2) # Pad the signal at the beginning and at the end with a sequence of # $\alpha/2$ h of elements of value $m = max(s(t))$. # # alpha_epochs = int(pd.Timedelta(alpha)/self.frequency) # alpha_epochs_half = int(alpha_epochs/2) # beta_epochs = int(pd.Timedelta(beta)/self.frequency) s_t_max = self.raw_data.max() x_p = _padded_data( self.raw_data, s_t_max, L_w_over_2, self.frequency ) # 1.2 Rank-order processing and decision logic # Apply a median filter to the $x_p$ series x_f = x_p.rolling(L_w, center=True, min_periods=L_w_over_2).median() # Rank-order thresholding # Create a series $y_1(n)$ where $y_1(n) = 1$ for $x_f(n)>p$, $0$ otw. # The threshold $p$ is the percentile of $x_f(n)$ corresponding to # $(h_s/24)\times 100\%$ p_threshold = x_f.quantile((pd.Timedelta(alpha)/pd.Timedelta('24h'))) y_1 = pd.Series(np.where(x_f > p_threshold, 1, 0), index=x_f.index) # 1.3 Morphological filtering # Morph. filter window length, L_p L_p = int(pd.Timedelta(beta)/self.frequency)+1 # Morph. filter, M_f M_f = np.ones(L_p) # Apply a morphological closing operation y_1_close = binary_closing(y_1, M_f).astype(int) # Apply a morphological opening operation y_1_close_and_open = binary_opening(y_1_close, M_f).astype(int) y_e = pd.Series(y_1_close_and_open, index=y_1.index) # 2. Processing and decision logic # This stage uses the estimates of the rest-activity periods # from the previous stage. # 2.1 Model-based data validation # Create a mask for sequences of more than $\zeta_{rest}$ zeros # during the rest periods # zeta_r = 30 # zeta_a = 2 if estimate_zeta: zeta_r = _estimate_zeta(self.data[y_e < 1], seq_length_max) zeta_a = _estimate_zeta(self.data[y_e > 0], seq_length_max) if verbose: print("CRESPO: estimated zeta@rest= {}".format(zeta_r)) print("CRESPO: estimated zeta@actv= {}".format(zeta_a)) # Find sequences of zeroes during the rest and the active periods # and mark as invalid sequences of more $\zeta_x$ zeroes. # Use y_e series as a filter for the rest periods mask_rest = _create_inactivity_mask( self.raw_data[y_e < 1], zeta_r, 1 ) # Use y_e series as a filter for the active periods mask_actv = _create_inactivity_mask( self.raw_data[y_e > 0], zeta_a, 1 ) mask = pd.concat([mask_actv, mask_rest], verify_integrity=True) # 2.2 Adaptative rank-order processing and decision logic # Replace masked values by NaN so that they are not taken into account # by the median filter. # Done before padding to avoid unaligned time series. x_nan = self.raw_data.copy() x_nan[mask < 1] = np.NaN # Pad the signal at the beginning and at the end with a sequence of 1h # of elements of value m = max(s(t)). x_sp = _padded_data( x_nan, s_t_max, L_p-1, self.frequency ) # Apply an adaptative median filter to the $x_{sp}$ series # no need to use a time-aware window as there is no time gap # in this time series by definition. x_fa = x_sp.rolling(L_w, center=True, min_periods=L_p-1).median() # The 'alpha' hour window is biased at the edges as it is not # symmetrical anymore. In the regions (start, start+alpha/2, # the median needs to be calculate by hand. # The range is start, start+alpha as the window is centered. median_start = x_sp.iloc[0:L_w].expanding( center=True ).median() median_end = x_sp.iloc[-L_w-1:-1].sort_index( ascending=False ).expanding(center=True).median()[::-1] # replace values in the original x_fa series with the new values # within the range (start, start+alpha/2) only. x_fa.iloc[0:L_w_over_2] = median_start.iloc[0:L_w_over_2] x_fa.iloc[-L_w_over_2-1:-1] = median_end.iloc[0:L_w_over_2] # restore original time range x_fa = x_fa[self.data.index[0]:self.data.index[-1]] p_threshold = x_fa.quantile((pd.Timedelta(alpha)/pd.Timedelta('24h'))) y_2 = pd.Series(np.where(x_fa > p_threshold, 1, 0), index=x_fa.index) # ### 2.3 Morphological filtering y_2_close = binary_closing( y_2, structure=np.ones(2*(L_p-1)+1) ).astype(int) y_2_open = binary_opening( y_2_close, structure=np.ones(2*(L_p-1)+1) ).astype(int) crespo = pd.Series( y_2_open, index=y_2.index) # Manual post-processing crespo.iloc[0] = 1 crespo.iloc[-1] = 1 return crespo
[docs] def Crespo_AoT( self, zeta=15, zeta_r=30, zeta_a=2, t=.33, alpha='8h', beta='1h', estimate_zeta=False, seq_length_max=100, verbose=False ): """Automatic identification of activity onset/offset times, based on the Crespo algorithm. Identification of the activity onset and offset times using the algorithm for automatic identification of activity-rest periods based on actigraphy, developped by Crespo et al. [1]_. Parameters ---------- zeta: int Maximum number of consecutive zeroes considered valid. Default is 15. zeta_r: int Maximum number of consecutive zeroes considered valid (rest). Default is 30. zeta_a: int Maximum number of consecutive zeroes considered valid (active). Default is 2. t: float Percentile for invalid zeroes. Default is 0.33. alpha: offset Average hours of sleep per night. Default is '8h'. beta: offset Length of the padding sequence used during the processing. Default is '1h'. estimate_zeta: Boolean If set to True, zeta values are estimated from the distribution of ratios of the number of series of consecutive zeroes to the number of series randomly chosen from the actigraphy data. Default is False. seq_length_max: int Maximal length of the aforementioned random series. Default is 100. verbose: If set to True, print the estimated values of zeta. Default is False. Returns ------- aot : (ndarray, ndarray) Arrays containing the estimated activity onset and offset times, respectively. References ---------- .. [1] Crespo, C., Aboy, M., Fernández, J. R., & Mojón, A. (2012). Automatic identification of activity–rest periods based on actigraphy. Medical & Biological Engineering & Computing, 50(4), 329–340. http://doi.org/10.1007/s11517-012-0875-y Examples -------- """ crespo = self.Crespo( zeta=zeta, zeta_r=zeta_r, zeta_a=zeta_a, t=t, alpha=alpha, beta=beta, estimate_zeta=estimate_zeta, seq_length_max=seq_length_max, verbose=verbose ) diff = crespo.diff(1) AonT = crespo[diff == 1].index AoffT = crespo[diff == -1].index return (AonT, AoffT)
[docs] def Roenneberg( self, trend_period='24h', min_trend_period='12h', threshold=0.15, min_seed_period='30Min', max_test_period='12h', r_consec_below='30Min', rsfreq=None ): """Automatic sleep detection. Identification of consolidated sleep episodes using the algorithm developped by Roenneberg et al. [1]_. Parameters ---------- trend_period: str, optional Time period of the rolling window used to extract the data trend. Default is '24h'. min_trend_period: str, optional Minimum time period required for the rolling window to produce a value. Values default to NaN otherwise. Default is '12h'. threshold: float, optional Fraction of the trend to use as a threshold for sleep/wake categorization. Default is '0.15' min_seed_period: str, optional Minimum time period required to identify a potential sleep onset. Default is '30Min'. max_test_period : str, optional Maximal period of the test series. Default is '12h' r_consec_below : str, optional Time range to consider, past the potential correlation peak when searching for the maximum correlation peak. Default is '30Min'. rsfreq: str, optional Resampling frequency used to evaluate the sleep periods. The final time series with rest/activity scores is returned with a frequency equal to one of the input data. If set to None, no resampling is performed. Default is None. Returns ------- rbg : pandas.core.Series Time series containing the estimated periods of rest (1) and activity (0). Notes ----- .. warning:: The performance of this algorithm has been evaluated for actigraphy data aggregated in 10-min bins [2]_. References ---------- .. [1] Roenneberg, T., Keller, L. K., Fischer, D., Matera, J. L., Vetter, C., & Winnebeck, E. C. (2015). Human Activity and Rest In Situ. In Methods in Enzymology (Vol. 552, pp. 257-283). http://doi.org/10.1016/bs.mie.2014.11.028 .. [2] Loock, A., Khan Sullivan, A., Reis, C., Paiva, T., Ghotbi, N., Pilz, L. K., Biller, A. M., Molenda, C., Vuori‐Brodowski, M. T., Roenneberg, T., & Winnebeck, E. C. (2021). Validation of the Munich Actimetry Sleep Detection Algorithm for estimating sleep–wake patterns from activity recordings. Journal of Sleep Research, April, 1–12. https://doi.org/10.1111/jsr.13371 Examples -------- """ if rsfreq is not None: rsdata = self.resampled_data(freq=rsfreq) else: rsdata = self.data rbg = roenneberg( rsdata, trend_period=trend_period, min_trend_period=min_trend_period, threshold=threshold, min_seed_period=min_seed_period, max_test_period=max_test_period, r_consec_below=r_consec_below ) if rsfreq is not None: rbg = rbg.asfreq(self.data.index.freq, method='pad') return rbg
[docs] def Roenneberg_AoT( self, trend_period='24h', min_trend_period='12h', threshold=0.15, min_seed_period='30Min', max_test_period='12h', r_consec_below='30Min', rsfreq=None ): """Automatic identification of activity onset/offset times, based on Roenneberg's algorithm. Identification of the activity onset and offset times using the algorithm for automatic identification of consolidated sleep episodes developped by Roenneberg et al. [1]_. Parameters ---------- trend_period: str, optional Time period of the rolling window used to extract the data trend. Default is '24h'. min_trend_period: str, optional Minimum time period required for the rolling window to produce a value. Values default to NaN otherwise. Default is '12h'. threshold: float, optional Fraction of the trend to use as a threshold for sleep/wake categorization. Default is '0.15' min_seed_period: str, optional Minimum time period required to identify a potential sleep onset. Default is '30Min'. max_test_period : str, optional Maximal period of the test series. Default is '12h' r_consec_below : str, optional Time range to consider, past the potential correlation peak when searching for the maximum correlation peak. Default is '30Min'. rsfreq: str, optional Resampling frequency used to evaluate the sleep periods. The final time series with rest/activity scores is returned with a frequency equal to one of the input data. If set to None, no resampling is performed. Default is None. Returns ------- aot : (ndarray, ndarray) Arrays containing the estimated activity onset and offset times, respectively. Notes ----- .. warning:: The performance of this algorithm has been evaluated for actigraphy data aggregated in 10-min bins [2]_. References ---------- .. [1] Roenneberg, T., Keller, L. K., Fischer, D., Matera, J. L., Vetter, C., & Winnebeck, E. C. (2015). Human Activity and Rest In Situ. In Methods in Enzymology (Vol. 552, pp. 257-283). http://doi.org/10.1016/bs.mie.2014.11.028 .. [2] Loock, A., Khan Sullivan, A., Reis, C., Paiva, T., Ghotbi, N., Pilz, L. K., Biller, A. M., Molenda, C., Vuori‐Brodowski, M. T., Roenneberg, T., & Winnebeck, E. C. (2021). Validation of the Munich Actimetry Sleep Detection Algorithm for estimating sleep–wake patterns from activity recordings. Journal of Sleep Research, April, 1–12. https://doi.org/10.1111/jsr.13371 Examples -------- """ rbg = self.Roenneberg( trend_period=trend_period, min_trend_period=min_trend_period, threshold=threshold, min_seed_period=min_seed_period, max_test_period=max_test_period, r_consec_below=r_consec_below, rsfreq=rsfreq ) diff = rbg.diff(1) AonT = rbg[diff == -1].index AoffT = rbg[diff == 1].index return (AonT, AoffT)
[docs] def SleepProfile( self, freq='15min', algo='Roenneberg', *args, **kwargs ): r"""Normalized sleep daily profile XXXXX Parameters ---------- freq: str, optional Resampling frequency. Default is '15min' algo: str, optional Sleep scoring algorithm to use. Default is 'Roenneberg'. *args Variable length argument list passed to the scoring algorithm. **kwargs Arbitrary keyword arguements passed to the scoring algorithm. Returns ------- sleep_prof: YYY Examples -------- >>> import pyActigraphy >>> rawAWD = pyActigraphy.io.read_raw_awd(fpath + 'SUBJECT_01.AWD') >>> raw.SleepProfile() ZZZZZZZZZZZZZZZZzZZZZZ """ # Retrieve sleep scoring function dynamically by name sleep_algo = getattr(self, algo) # Score activity sleep_scoring = sleep_algo(*args, **kwargs) # Sleep profile over 24h sleep_prof = _average_daily_activity(data=sleep_scoring, cyclic=False) # Resampled sleep profile rs_sleep_profile = sleep_prof.resample(freq).mean() return rs_sleep_profile
[docs] def SleepRegularityIndex( self, freq='15min', bin_threshold=None, algo='Roenneberg', *args, **kwargs ): r""" Sleep regularity index Likelihood that any two time-points (epoch-by-epoch) 24 hours apart are in the same sleep/wake state, across all days. This index is originally defined in [1]_ and validated in older adults in [2]_. Parameters ---------- freq: str, optional Resampling frequency. Default is '15min' bin_threshold: bool, optional If bin_threshold is not set to None, scoring data above this threshold are set to 1 and to 0 otherwise. Default is None. algo: str, optional Sleep scoring algorithm to use. Default is 'Roenneberg'. *args Variable length argument list passed to the scoring algorithm. **kwargs Arbitrary keyword arguements passed to the scoring algorithm. Returns ------- sri: float Notes ----- The sleep regularity index (SRI) is defined as: .. math:: SRI = -100 + \frac{200}{M(N-1)} \sum_{j=1}^M\sum_{i=1}^N \delta(s_{i,j}, s_{i+1,j}) with: :math:`\delta(s_{i,j}, s_{i+1,j}) = 1` if :math:`s_{i,j} = s_{i+1,j}` and 0 otherwise. References ---------- .. [1] Phillips, A. J. K., Clerx, W. M., O’Brien, C. S., Sano, A., Barger, L. K., Picard, R. W., … Czeisler, C. A. (2017). Irregular sleep/wake patterns are associated with poorer academic performance and delayed circadian and sleep/wake timing. Scientific Reports, 7(1), 1–13. https://doi.org/10.1038/s41598-017-03171-4 .. [2] Lunsford-Avery, J. R., Engelhard, M. M., Navar, A. M., & Kollins, S. H. (2018). Validation of the Sleep Regularity Index in Older Adults and Associations with Cardiometabolic Risk. Scientific Reports, 8(1), 14158. https://doi.org/10.1038/s41598-018-32402-5 Examples -------- """ # Retrieve sleep scoring function dynamically by name sleep_algo = getattr(self, algo) # Score activity sleep_scoring = sleep_algo(*args, **kwargs) # SRI sleep_regularity_index = sri(sleep_scoring, bin_threshold) return sleep_regularity_index
[docs] def SleepMidPoint( self, freq='15min', bin_threshold=None, to_td=True, algo='Roenneberg', *args, **kwargs ): r""" Sleep midpoint Center of the mean sleep periods Parameters ---------- freq: str, optional Resampling frequency. Default is '15min' bin_threshold: bool, optional If bin_threshold is not set to None, scoring data above this threshold are set to 1 and to 0 otherwise. Default is None. to_td: bool, optional If set to true, the sleep midpoint is returned as a Timedelta. Otherwise, it represents the number of minutes since midnight. algo: str, optional Sleep scoring algorithm to use. Default is 'Roenneberg'. \*args Variable length argument list passed to the scoring algorithm. \*\*kwargs Arbitrary keyword arguements passed to the scoring algorithm. Returns ------- smp: float or Timedelta Notes ----- Sleep midpoint (SMP) is an index of sleep timing and is calculated as the following [1]_: .. math:: SMP = \frac{1440}{2\pi} arctan2\left( \sum_{j=1}^M\sum_{i=1}^N s_{i,j} \times sin\left(\frac{2\pi t_i}{1440}\right), \sum_{j=1}^M\sum_{i=1}^N s_{i,j} \times cos\left(\frac{2\pi t_i}{1440}\right) \right) with: * :math:`t_j`, time of day in minutes at epoch j, * :math:`\delta(s_{i,j}, s_{i+1,j}) = 1` if :math:`s_{i,j} = s_{i+1,j}` and 0 otherwise. References ---------- .. [1] Lunsford-Avery, J. R., Engelhard, M. M., Navar, A. M., & Kollins, S. H. (2018). Validation of the Sleep Regularity Index in Older Adults and Associations with Cardiometabolic Risk. Scientific Reports, 8(1), 14158. https://doi.org/10.1038/s41598-018-32402-5 Examples -------- """ # Retrieve sleep scoring function dynamically by name sleep_algo = getattr(self, algo) # Score activity sleep_scoring = sleep_algo(*args, **kwargs) # Sleep midpoint smp = sleep_midpoint(sleep_scoring, bin_threshold) return pd.Timedelta(smp, unit='min') if to_td is True else smp