{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Exploring and Visualizing CMI Sleep Dataset \nIn this notebook I perform an exploratory analaysis and visualization of the Child Mind Institute sleep monitoring dataset, as part of the [Detect Sleep States](https://www.kaggle.com/competitions/child-mind-institute-detect-sleep-states/overview) competition.","metadata":{}},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport polars as pl\nimport datetime \nfrom tqdm import tqdm\n\nimport plotly.express as px\nfrom plotly.subplots import make_subplots\nimport plotly.graph_objects as go","metadata":{"execution":{"iopub.status.busy":"2023-10-02T17:46:38.059307Z","iopub.execute_input":"2023-10-02T17:46:38.059668Z","iopub.status.idle":"2023-10-02T17:46:39.435728Z","shell.execute_reply.started":"2023-10-02T17:46:38.059639Z","shell.execute_reply":"2023-10-02T17:46:39.434609Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Importing data\nSince the time series dataset is very large, we will use [`polars`](https://www.pola.rs) in this notebook as opposed to the more conventional `pandas` library. We begin by importing the datasets and extracting the time data from the 'timestamp' field. Wherever possible we will cast our columns into efficient datatypes to save on memory.","metadata":{}},{"cell_type":"code","source":"# Importing data \ndt_transforms = [\n    pl.col('timestamp').str.to_datetime(), \n    (pl.col('timestamp').str.to_datetime().dt.year()-2000).cast(pl.UInt8).alias('year'), \n    pl.col('timestamp').str.to_datetime().dt.month().cast(pl.UInt8).alias('month'),\n    pl.col('timestamp').str.to_datetime().dt.day().cast(pl.UInt8).alias('day'), \n    pl.col('timestamp').str.to_datetime().dt.hour().cast(pl.UInt8).alias('hour')\n]\n\ndata_transforms = [\n    pl.col('anglez').cast(pl.Int8), # Casting anglez to 8 bit integer\n    (pl.col('enmo')*1000).cast(pl.UInt16), # Convert enmo to 16 bit uint\n]\n\ntrain_series = pl.scan_parquet('/kaggle/input/child-mind-institute-detect-sleep-states/train_series.parquet').with_columns(\n    dt_transforms + data_transforms\n    )\n\ntrain_events = pl.read_csv('/kaggle/input/child-mind-institute-detect-sleep-states/train_events.csv').with_columns(\n    dt_transforms\n    )\n\ntest_series = pl.scan_parquet('/kaggle/input/child-mind-institute-detect-sleep-states/test_series.parquet').with_columns(\n    dt_transforms + data_transforms\n    )\n\n# Getting series ids as a list for convenience\nseries_ids = train_events['series_id'].unique(maintain_order=True).to_list()","metadata":{"execution":{"iopub.status.busy":"2023-10-02T17:46:39.533205Z","iopub.execute_input":"2023-10-02T17:46:39.534041Z","iopub.status.idle":"2023-10-02T17:46:39.838616Z","shell.execute_reply.started":"2023-10-02T17:46:39.534005Z","shell.execute_reply":"2023-10-02T17:46:39.837387Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Visualizing time series\nLet's create a function to visualize given series in a way which displays our two main datapoints, `anglez` and `enmo`, along with the onset and wakeup events which we would like to predict - or equivalently, combining the two, the sleeping periods vs. the wakeful periods of the time series. Our visualization should also be able to handle missing data by marking days with missing values.\n\nTo begin, let's find out whether the number of onset events always equals the number of wakeup events. In principle these numbers should always be the same.","metadata":{}},{"cell_type":"code","source":"onset_counts = train_events.filter(pl.col('event')=='onset').group_by('series_id').count().sort('series_id')['count']\nwakeup_counts = train_events.filter(pl.col('event')=='wakeup').group_by('series_id').count().sort('series_id')['count']\n\ncounts = pl.DataFrame({'series_id':sorted(series_ids), 'onset_counts':onset_counts, 'wakeup_counts':wakeup_counts})\ncount_mismatches = counts.filter(counts['onset_counts'] != counts['wakeup_counts'])\n\nprint(f'Series where the number of onsets is not the same as the number of wakeups: {count_mismatches.shape[0]}')","metadata":{"execution":{"iopub.status.busy":"2023-10-02T17:46:40.280952Z","iopub.execute_input":"2023-10-02T17:46:40.281515Z","iopub.status.idle":"2023-10-02T17:46:40.319697Z","shell.execute_reply.started":"2023-10-02T17:46:40.281483Z","shell.execute_reply":"2023-10-02T17:46:40.318965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So there are four series where there is a mismatch. Let's remove these from our training set.","metadata":{}},{"cell_type":"code","source":"train_series = train_series.filter(~pl.col('series_id').is_in(count_mismatches['series_id']))\ntrain_events = train_events.filter(~pl.col('series_id').is_in(count_mismatches['series_id']))\n\n# Updating series_ids list\nseries_ids = train_events.drop_nulls()['series_id'].unique(maintain_order=True).to_list()","metadata":{"execution":{"iopub.status.busy":"2023-10-02T17:46:41.068598Z","iopub.execute_input":"2023-10-02T17:46:41.069489Z","iopub.status.idle":"2023-10-02T17:46:41.088305Z","shell.execute_reply.started":"2023-10-02T17:46:41.069444Z","shell.execute_reply":"2023-10-02T17:46:41.086963Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We also note the presence of missing labels in some of the series. That is, for some series there are days in which no onsets or wakeups are recorded, even though the time series data indicate periods where the participant likely was sleeping. Let's see how frequently this appears in our data. ","metadata":{}},{"cell_type":"code","source":"missing_data = train_events.group_by('series_id', maintain_order=True).agg(pl.col('step').map_elements(lambda x: x.is_null().any()).alias('MissingValues'))\nprint(f\"Series with missing data: {missing_data['MissingValues'].sum()}/{train_events['series_id'].n_unique()}\")","metadata":{"execution":{"iopub.status.busy":"2023-10-02T17:46:41.905298Z","iopub.execute_input":"2023-10-02T17:46:41.905940Z","iopub.status.idle":"2023-10-02T17:46:41.945947Z","shell.execute_reply.started":"2023-10-02T17:46:41.905895Z","shell.execute_reply":"2023-10-02T17:46:41.944804Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Unfortunately a large majority of series contain days with missing data. We will take this into consideration when training our model.","metadata":{}},{"cell_type":"markdown","source":"Now that we can be certain that the onsets and wakeup numbers match, we can define a plotting function that will display our two primary datapoints - anglez and enmo - on a time series for a given series_id, with the sleep stages highlighted. For completeness, we will also flag any days where there are missing or null values. ","metadata":{}},{"cell_type":"code","source":"def getdata(i, every=1) : \n    '''\n    Get data for a given sample series. \n    '''\n    \n    idx = series_ids[i] if type(i) == int else i\n    series = train_series.filter((pl.col('series_id') == idx) & (pl.col('step') % every == 0)).collect()\n    events = train_events.filter(pl.col('series_id') == idx)\n    \n    return idx, series, events\n\ndef plotdata(idx, series, events) :\n    '''\n    Plots a time series of anglez and enmo for the i-th series, using a \n    given frequency of samples, with wakeup and onset events marked. \n    every=1 uses every sample at 5 minute intervals, and every=12 gives \n    minute-to-minute time series. \n    '''\n    \n    fig = make_subplots(specs=[[{'secondary_y': True}]])\n    \n    # Plotting time series data\n    fig.add_trace(go.Scatter(x=series['timestamp'], y=series['anglez'], mode='lines', name='anglez'), secondary_y=False)\n    fig.add_trace(go.Scatter(x=series['timestamp'], y=series['enmo'], mode='lines', name='enmo'), secondary_y=True)\n    \n    # Plotting sleeping periods\n    onsets = events.filter((pl.col('event') == 'onset') & (pl.col('timestamp') != None))['timestamp'].to_list()\n    wakeups = events.filter((pl.col('event') == 'wakeup') & (pl.col('timestamp') != None))['timestamp'].to_list()\n    for onset, wakeup in zip(onsets, wakeups) :\n        fig.add_vrect(\n            x0=onset, x1=wakeup,\n            fillcolor='gray', opacity=0.2,\n            layer='below', line_width=0\n        )\n        \n    # Add dummy trace for the legend\n    fig.add_trace(go.Scatter(\n        x=[None],\n        y=[None],\n        mode='lines',\n        line=dict(color='gray', width=0),\n        fill='tozeroy',\n        fillcolor='gray',\n        opacity=0.2,\n        name=\"Sleeping periods\"\n    ))\n        \n    # Plotting days with missing data\n    start_date = events['timestamp'].min().replace(hour=0, minute=0, second=0, microsecond=0)\n    null_nights = events.filter(pl.col('step').is_null())['night'].unique().to_list()\n    for n in null_nights : \n        fig.add_vrect(\n            x0=start_date + datetime.timedelta(days=n-1), x1=start_date + datetime.timedelta(days=n),\n            fillcolor='violet', opacity=0.2,\n            layer='below', line_width=0\n        )\n    \n    # Add dummy trace for the legend\n    fig.add_trace(go.Scatter(\n        x=[None],\n        y=[None],\n        mode='lines',\n        line=dict(color='violet', width=0),\n        fill='tozeroy',\n        fillcolor='violet',\n        opacity=0.2,\n        name=\"Missing data\"\n    ))\n        \n    fig.update_layout(\n        title_text=f\"Accelerometer Data<br><sub>ID: {idx}</sub>\", \n        #template = 'ggplot2',\n        xaxis=dict(showgrid=False), yaxis=dict(showgrid=False), yaxis2=dict(showgrid=False)\n    )\n    \n    fig.show()\n    \n#data = getdata(10)\n\n#plotdata(*data)","metadata":{"execution":{"iopub.status.busy":"2023-10-02T17:46:43.093996Z","iopub.execute_input":"2023-10-02T17:46:43.094739Z","iopub.status.idle":"2023-10-02T17:46:43.108091Z","shell.execute_reply.started":"2023-10-02T17:46:43.094684Z","shell.execute_reply":"2023-10-02T17:46:43.107335Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"After plotting a few series samples, we can see some very clear patterns that tend to distinguish sleeping periods from wakeful periods. \n- During wakeful periods `anglez` varies constantly, rarely ever staying at the same value for consecutive measurements. During sleeping periods there tend to be periods of low to zero variability in `anglez` where the person will change positions and then remain there for many consecutive measurements. \n- During wakeful periods `enmo` also has much higher variability and reaches much greater values, while during sleeping periods `enmo` values stay very close to 0.\n- Days with missing data tend to have much smaller variations in `anglez`, as well as `enmo`. \n\n","metadata":{}},{"cell_type":"markdown","source":"## Feature Engineering","metadata":{}},{"cell_type":"markdown","source":"Let's start by looking at the series without any null values to look for features correlated with sleep/wakefullness, without being impacted by missing data.","metadata":{}},{"cell_type":"code","source":"i = 0\n\nfull_data_ids = missing_data.filter(pl.col('MissingValues')==False)['series_id'].to_list()\nidx, series, events = getdata(full_data_ids[i])\nplotdata(idx, series, events)\n\nonsets = events.filter((pl.col('event') == 'onset') & (pl.col('timestamp') != None))['timestamp'].to_list()\nwakeups = events.filter((pl.col('event') == 'wakeup') & (pl.col('timestamp') != None))['timestamp'].to_list()","metadata":{"execution":{"iopub.status.busy":"2023-10-02T17:46:44.873434Z","iopub.execute_input":"2023-10-02T17:46:44.874213Z","iopub.status.idle":"2023-10-02T17:46:58.641723Z","shell.execute_reply.started":"2023-10-02T17:46:44.874170Z","shell.execute_reply":"2023-10-02T17:46:58.640762Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's put together some features to see if we can find any that strongly correlate with sleeping phases. ","metadata":{}},{"cell_type":"code","source":"features = [sum([(onset <= pl.col('timestamp')) & (pl.col('timestamp') <= wakeup) for onset, wakeup in zip(onsets, wakeups)]).alias('asleep')]\nfeature_cols = ['asleep']\n\nrolling_durs = [1, 30, 60*3] # Rolling aggregate durations, in minutes\n\nfor var in ['anglez', 'enmo'] :\n    features += [pl.col(var).alias(var), \n                 pl.col(var).diff().fill_null(0).abs().alias(f'{var}_1v')\n                ]\n    feature_cols += [var, f'{var}_1v']\n    \n    # Rolling aggregates\n    for dur in rolling_durs :\n        features += [\n            pl.col(var).rolling_mean(dur * 12, center=True, min_periods=1).alias(f'{var}_{dur}m_mean'), # Rolling mean\n            pl.col(var).rolling_max(dur * 12, center=True, min_periods=1).alias(f'{var}_{dur}m_max'), # Rolling max\n            \n            # Rolling total variations\n            pl.col(var).diff().fill_null(0).abs().rolling_mean(dur * 12, center=True, min_periods=1).alias(f'{var}_1v_{dur}m_mean'),\n            pl.col(var).diff().fill_null(0).abs().rolling_max(dur * 12, center=True, min_periods=1).alias(f'{var}_1v_{dur}m_max')\n        ]\n\n        feature_cols += [\n            f'{var}_{dur}m_mean', f'{var}_{dur}m_max', f'{var}_1v_{dur}m_mean', f'{var}_1v_{dur}m_max'\n        ]\n        \n# Lets also include enmo x anglez features\nfeatures += [\n    (pl.col('anglez') * pl.col('enmo')).alias('anglez.enmo'), \n    (pl.col('anglez') * pl.col('enmo')).alias('anglez.enmo').diff().fill_null(0).abs().alias('anglez.enmo_1v')\n]\n\nfeature_cols += [\n    'anglez.enmo', 'anglez.enmo_1v'\n]\n\nfor dur in rolling_durs :\n    features += [\n        (pl.col('anglez') * pl.col('enmo')).rolling_mean(dur * 12, center=True, min_periods=1).alias(f'anglez.enmo_{dur}m_mean'), # Rolling mean\n        (pl.col('anglez') * pl.col('enmo')).rolling_max(dur * 12, center=True, min_periods=1).alias(f'anglez.enmo_{dur}m_max'), # Rolling max\n\n        # Rolling total variations\n        (pl.col('anglez') * pl.col('enmo')).diff().fill_null(0).abs().rolling_mean(dur * 12, center=True, min_periods=1).alias(f'anglez.enmo_1v_{dur}m_mean'),\n        (pl.col('anglez') * pl.col('enmo')).diff().fill_null(0).abs().rolling_max(dur * 12, center=True, min_periods=1).alias(f'anglez.enmo_1v_{dur}m_max')\n    ]\n\n    feature_cols += [\n        f'anglez.enmo_{dur}m_mean', f'anglez.enmo_{dur}m_max', f'anglez.enmo_1v_{dur}m_mean', f'anglez.enmo_1v_{dur}m_max'\n    ]\n\ncorr_df = series.with_columns(features)[feature_cols].corr()","metadata":{"execution":{"iopub.status.busy":"2023-10-02T18:12:20.821279Z","iopub.execute_input":"2023-10-02T18:12:20.821633Z","iopub.status.idle":"2023-10-02T18:12:21.350561Z","shell.execute_reply.started":"2023-10-02T18:12:20.821605Z","shell.execute_reply":"2023-10-02T18:12:21.349363Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = series.with_columns(features)","metadata":{"execution":{"iopub.status.busy":"2023-10-02T18:24:41.168889Z","iopub.execute_input":"2023-10-02T18:24:41.169758Z","iopub.status.idle":"2023-10-02T18:24:41.390340Z","shell.execute_reply.started":"2023-10-02T18:24:41.169715Z","shell.execute_reply":"2023-10-02T18:24:41.389470Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cols.sort()\ncols","metadata":{"execution":{"iopub.status.busy":"2023-10-02T18:40:03.783852Z","iopub.execute_input":"2023-10-02T18:40:03.784396Z","iopub.status.idle":"2023-10-02T18:40:03.791067Z","shell.execute_reply.started":"2023-10-02T18:40:03.784353Z","shell.execute_reply":"2023-10-02T18:40:03.789940Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(r'Z-Angle \\times ENMO')","metadata":{"execution":{"iopub.status.busy":"2023-10-02T18:43:53.653343Z","iopub.execute_input":"2023-10-02T18:43:53.653683Z","iopub.status.idle":"2023-10-02T18:43:53.659692Z","shell.execute_reply.started":"2023-10-02T18:43:53.653656Z","shell.execute_reply":"2023-10-02T18:43:53.658482Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from plotly.subplots import make_subplots\n\nN = (len(feature_cols)-1)//3\n\n#fig = make_subplots(cols=1, rows=3)\nvariables = ['Z-Angle', 'ENMO', 'Z-Angle × ENMO']\nfor c in range(3) :\n    fig = go.Figure()\n    cols = feature_cols[:1] + feature_cols[c*N + 1:(c+1)*N + 1]\n    fig.add_trace(\n        go.Heatmap(\n            x = cols,\n            y = cols,\n            z = df[cols].corr(),\n            text=df[cols].corr(),\n            texttemplate='%{text:.2f}',\n            zmin=-1,zmax=1\n        )\n    )\n    fig.update_layout(title_text=f\"{variables[c]} correlation map<br><sub>ID: {idx}<sub>\")\n    fig.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-02T18:47:03.694233Z","iopub.execute_input":"2023-10-02T18:47:03.694618Z","iopub.status.idle":"2023-10-02T18:47:04.258261Z","shell.execute_reply.started":"2023-10-02T18:47:03.694584Z","shell.execute_reply":"2023-10-02T18:47:04.257125Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"features = [sum([(onset <= pl.col('timestamp')) & (pl.col('timestamp') <= wakeup) for onset, wakeup in zip(onsets, wakeups)]).alias('asleep')]\nfeature_cols = ['asleep']\n\nrolling_durs = [1, 30, 60*3] # Rolling aggregate durations, in minutes\n\nfor var in ['anglez', 'enmo'] :\n    features += [pl.col(var).alias(var), pl.col(var).diff().abs().alias(f'{var}_1v')]\n    feature_cols += [var, f'{var}_1v']\n    \n    # Rolling aggregates\n    for dur in rolling_durs :\n        features += [\n            pl.col(var).rolling_mean(dur * 12, center=True, min_periods=1).alias(f'{var}_{dur}m_mean'), # Rolling mean\n            pl.col(var).rolling_max(dur * 12, center=True, min_periods=1).alias(f'{var}_{dur}m_max'), # Rolling max\n            \n            # Rolling total variations\n            pl.col(var).diff().abs().rolling_mean(dur * 12, center=True, min_periods=1).alias(f'{var}_1v_{dur}m_mean'),\n            pl.col(var).diff().abs().rolling_max(dur * 12, center=True, min_periods=1).alias(f'{var}_1v_{dur}m_max')\n        ]\n\n        feature_cols += [\n            f'{var}_{dur}m_mean', f'{var}_{dur}m_max', f'{var}_1v_{dur}m_mean', f'{var}_1v_{dur}m_max'\n        ]\n        \n# Lets also include enmo x anglez features\nfeatures += [\n    (pl.col('anglez') * pl.col('enmo')).alias('anglez.enmo'), \n    (pl.col('anglez') * pl.col('enmo')).alias('anglez.enmo').diff().abs().alias('anglez.enmo_1v')\n]\n\nfeature_cols += [\n    'anglez.enmo', 'anglez.enmo_1v'\n]\n\nfor dur in rolling_durs :\n    features += [\n        (pl.col('anglez') * pl.col('enmo')).rolling_mean(dur * 12, center=True, min_periods=1).alias(f'anglez.enmo_{dur}m_mean'), # Rolling mean\n        (pl.col('anglez') * pl.col('enmo')).rolling_max(dur * 12, center=True, min_periods=1).alias(f'anglez.enmo_{dur}m_max'), # Rolling max\n\n        # Rolling total variations\n        (pl.col('anglez') * pl.col('enmo')).diff().abs().rolling_mean(dur * 12, center=True, min_periods=1).alias(f'anglez.enmo_1v_{dur}m_mean'),\n        (pl.col('anglez') * pl.col('enmo')).diff().abs().rolling_max(dur * 12, center=True, min_periods=1).alias(f'anglez.enmo_1v_{dur}m_max')\n    ]\n\n    feature_cols += [\n        f'anglez.enmo_{dur}m_mean', f'anglez.enmo_{dur}m_max', f'anglez.enmo_1v_{dur}m_mean', f'anglez.enmo_1v_{dur}m_max'\n    ]\n\ncorr_df = series.with_columns(features)\n\n\nfig = go.Figure()\nfig.add_trace(\n    go.Heatmap(\n        x = cols,\n        y = cols,\n        z = corr_df[cols].drop_nulls().corr(),\n        text=corr_df[cols].drop_nulls().corr(),\n        texttemplate='%{text:.2f}'\n    )\n)\nfig.update_layout(title_text=f\"Feature correlation map<br><sub>ID: {idx}<sub>\")\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-02T17:46:58.643034Z","iopub.execute_input":"2023-10-02T17:46:58.643680Z","iopub.status.idle":"2023-10-02T17:46:59.097407Z","shell.execute_reply.started":"2023-10-02T17:46:58.643651Z","shell.execute_reply":"2023-10-02T17:46:59.096526Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Some observations: \n- The longer 30 minute rolling windows outperformed both the minute rolling window and the individual datapoints by a wide margin.\n- The feature correlated most closely with sleeping is the 30 minute rolling average of the first order variation of anglez, with a 75% negative correlation.\n- Both the enmo 30 minute rolling average and the enmo 30 minute rolling first order variation offered strong negative correlations. \n- The rolling 30 minute multiplicative combitionation of anglez and enmo offered a strong positive correlation as well. ","metadata":{}},{"cell_type":"markdown","source":"Let's see how these four correlations hold up across the other sample series with full data.","metadata":{}},{"cell_type":"code","source":"cols = ['asleep', 'anglez_1v_r30', 'enmo_1v_r30', 'enmo_r30', 'anglezxenmo_r30']\nsleep_corr = pl.DataFrame()\n\nfor idx in tqdm(full_data_ids) : \n    idx, series, events = getdata(idx)\n\n    onsets = events.filter((pl.col('event') == 'onset') & (pl.col('timestamp') != None))['timestamp'].to_list()\n    wakeups = events.filter((pl.col('event') == 'wakeup') & (pl.col('timestamp') != None))['timestamp'].to_list()\n    \n    corr_df = series.with_columns(\n        sum([(onset <= pl.col('timestamp')) & (pl.col('timestamp') <= wakeup) for onset, wakeup in zip(onsets, wakeups)]).alias('asleep'), \n        pl.col('enmo').rolling_mean(12*30, center=True).alias('enmo_r30'), # Rolling mean of anglez for 30 minutes\n        (pl.col('enmo') * pl.col('anglez')).rolling_mean(12*30, center=True).alias('anglezxenmo_r30'), # Rolling mean of anglez times enmo for 30 minutes\n        (pl.col('anglez').shift(-1) - pl.col('anglez')).abs().rolling_mean(12*30, center=True).alias('anglez_1v_r30'), # 1st order variation of anglez for 30 min\n        (pl.col('enmo').shift(-1) - pl.col('enmo')).abs().rolling_mean(12*30, center=True).alias('enmo_1v_r30') # First order varation of enmo for 30 minutes\n    )\n    \n    sleep_corr = sleep_corr.vstack(corr_df[cols].drop_nulls().corr()[0].drop('asleep').with_columns(pl.lit(idx).alias('ID')))","metadata":{"execution":{"iopub.status.busy":"2023-09-13T23:46:49.089175Z","iopub.execute_input":"2023-09-13T23:46:49.089561Z","iopub.status.idle":"2023-09-13T23:53:27.511579Z","shell.execute_reply.started":"2023-09-13T23:46:49.089533Z","shell.execute_reply":"2023-09-13T23:53:27.510414Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_cols = [col for col in cols if col != 'asleep']\n\nfig = go.Figure()\nfor col in plot_cols : \n    fig.add_trace(go.Violin(y=sleep_corr[col],\n                            name=col,\n                            box_visible=True,\n                            meanline_visible=True))\n    \nfig.update_layout(title_text='Feature correlation with sleep across all full-data time series')\n    \nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-09-13T23:59:47.421967Z","iopub.execute_input":"2023-09-13T23:59:47.422393Z","iopub.status.idle":"2023-09-13T23:59:47.441259Z","shell.execute_reply.started":"2023-09-13T23:59:47.422351Z","shell.execute_reply":"2023-09-13T23:59:47.440077Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Our rolling average anglez $\\times$ enmo feature's correlation with sleep does not appear to be robust across the other time series, so we will discard it and move forward with the remaining three features, whose correlations do appear to be robust. ","metadata":{}},{"cell_type":"markdown","source":"### Exploiting periodicity","metadata":{}},{"cell_type":"code","source":"px.histogram(train_events, x='hour', color='event', opacity=0.5, barmode='overlay')","metadata":{"execution":{"iopub.status.busy":"2023-09-14T01:30:39.484454Z","iopub.execute_input":"2023-09-14T01:30:39.485832Z","iopub.status.idle":"2023-09-14T01:30:39.576781Z","shell.execute_reply.started":"2023-09-14T01:30:39.485789Z","shell.execute_reply":"2023-09-14T01:30:39.575885Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training Random Forest Model","metadata":{}},{"cell_type":"code","source":"import torch\nfrom torch.utils.data import Dataset, DataLoader, TensorDataset\nimport torch.nn as nn\nimport torch.nn.functional as F\n\nfrom sklearn.preprocessing import StandardScaler","metadata":{"execution":{"iopub.status.busy":"2023-09-16T19:53:27.711893Z","iopub.execute_input":"2023-09-16T19:53:27.712350Z","iopub.status.idle":"2023-09-16T19:53:31.522405Z","shell.execute_reply.started":"2023-09-16T19:53:27.712319Z","shell.execute_reply":"2023-09-16T19:53:31.521130Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"i = 10\nidx = series_ids[i]\nchunk_df = train_series.filter(pl.col('series_id') == idx).collect().sort('timestamp')","metadata":{"execution":{"iopub.status.busy":"2023-09-14T02:29:22.932244Z","iopub.execute_input":"2023-09-14T02:29:22.933109Z","iopub.status.idle":"2023-09-14T02:29:34.176646Z","shell.execute_reply.started":"2023-09-14T02:29:22.933066Z","shell.execute_reply":"2023-09-14T02:29:34.175407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"i = 10\nidx = series_ids[i]\n\nchunk_df = train_series.filter(pl.col('series_id') == idx).collect().sort('timestamp').with_columns( \n    pl.col('enmo').rolling_mean(12*30, center=True).alias('enmo_r30'), # Rolling mean of anglez for 30 minutes\n    pl.col('enmo').diff().abs().rolling_mean(12*30, center=True).alias('enmo_1v_r30'), # First order varation of enmo for 30 minutes\n    pl.col('anglez').diff().abs().rolling_mean(12*30, center=True).alias('anglez_1v_r30') # 1st order variation of anglez for 30 min\n)","metadata":{"execution":{"iopub.status.busy":"2023-09-14T02:37:53.951947Z","iopub.execute_input":"2023-09-14T02:37:53.952387Z","iopub.status.idle":"2023-09-14T02:38:05.015684Z","shell.execute_reply.started":"2023-09-14T02:37:53.952356Z","shell.execute_reply":"2023-09-14T02:38:05.014589Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"chunk_df ","metadata":{"execution":{"iopub.status.busy":"2023-09-14T02:38:07.619704Z","iopub.execute_input":"2023-09-14T02:38:07.620177Z","iopub.status.idle":"2023-09-14T02:38:07.645199Z","shell.execute_reply.started":"2023-09-14T02:38:07.620141Z","shell.execute_reply":"2023-09-14T02:38:07.643677Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class SleepDataset(Dataset) :\n    def __init__(self, series_data, events_data, chunk_size):\n        self.series_data = series_data\n        self.events_data = events_data\n        self.chunk_size = chunk_size\n        self.series_ids = events_data.select('series_id').unique(maintain_order=True).to_list()\n        self.scaler = StandardScaler()\n\n    def __len__(self):\n        return len(self.series_ids)\n\n    def __getitem__(self, i):\n        idx = self.series_ids[i]\n        \n        chunk_df = self.series_data.filter(pl.col('series_id') == idx).collect().sort('timestamp').with_columns( \n            pl.col('enmo').rolling_mean(12*30, center=True).alias('enmo_r30'), # Rolling mean of anglez for 30 minutes\n            (pl.col('anglez').shift(-1) - pl.col('anglez')).abs().rolling_mean(12*30, center=True).alias('anglez_1v_r30'), # 1st order variation of anglez for 30 min\n            (pl.col('enmo').shift(-1) - pl.col('enmo')).abs().rolling_mean(12*30, center=True).alias('enmo_1v_r30') # First order varation of enmo for 30 minutes\n        )\n        \n        series_tensor = torch.tensor(chunk_df)\n        \n        # Determine the range for the current chunk\n        start_idx = idx * self.chunk_size\n        end_idx = min(start_idx + self.chunk_size, self.total_rows)\n        \n        # Fetch the chunk from the Polars DataFrame\n        chunk_df = self.pl_lazy_df.filter(pl.col(\"index_column_name\").between(start_idx, end_idx)).collect()\n\n        # Convert the chunk to appropriate tensors (Modify this based on your data and model)\n        data_tensor = torch.tensor(chunk_df[\"data_column_name\"].to_numpy())\n        label_tensor = torch.tensor(chunk_df[\"label_column_name\"].to_numpy())\n        \n        return data_tensor, label_tensor","metadata":{},"execution_count":null,"outputs":[]}]}