{"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":"# Plot events as well as time series on spectrograms\nThis notebook is based on https://www.kaggle.com/code/gvyshnya/plotting-big-data-sleeping-series/notebook","metadata":{}},{"cell_type":"code","source":"from IPython.display import HTML, Markdown\nimport time\n\nhandle = display(HTML(\"\"\"\"\"\"), display_id='html_marquee1')\ntime.sleep(2)\nhandle = display(HTML(\"\"\"\"\"\"), display_id='html_marquee1', update=True)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-10-08T01:04:51.558253Z","iopub.execute_input":"2023-10-08T01:04:51.558573Z","iopub.status.idle":"2023-10-08T01:04:53.604702Z","shell.execute_reply.started":"2023-10-08T01:04:51.558551Z","shell.execute_reply":"2023-10-08T01:04:53.603461Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install duckdb --quiet","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-10-08T01:04:53.606545Z","iopub.execute_input":"2023-10-08T01:04:53.606907Z","iopub.status.idle":"2023-10-08T01:05:03.488324Z","shell.execute_reply.started":"2023-10-08T01:04:53.606877Z","shell.execute_reply":"2023-10-08T01:05:03.486695Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Importing gc module\nimport gc\n\nimport pandas as pd\nimport numpy as np\nimport datetime as dt\n\n\nfrom scipy.signal import spectrogram\n\nimport matplotlib.pyplot as plt\nimport matplotlib.cm as cm\nimport matplotlib.dates as mdates\nimport seaborn as sns\n\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\nimport plotly.express as px\nimport plotly.offline","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-10-08T01:05:03.489748Z","iopub.execute_input":"2023-10-08T01:05:03.490057Z","iopub.status.idle":"2023-10-08T01:05:05.116979Z","shell.execute_reply.started":"2023-10-08T01:05:03.490031Z","shell.execute_reply":"2023-10-08T01:05:05.116039Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import random","metadata":{"execution":{"iopub.status.busy":"2023-10-08T01:18:19.815125Z","iopub.execute_input":"2023-10-08T01:18:19.815475Z","iopub.status.idle":"2023-10-08T01:18:19.821128Z","shell.execute_reply.started":"2023-10-08T01:18:19.815453Z","shell.execute_reply":"2023-10-08T01:18:19.819689Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import duckdb\n\nparquet_table = '/kaggle/input/child-mind-institute-detect-sleep-states/train_series.parquet'\ntesting_parquet_table = '/kaggle/input/child-mind-institute-detect-sleep-states/test_series.parquet'\ntraining_events_table = '/kaggle/input/child-mind-institute-detect-sleep-states/train_events.csv'\n\ncon = duckdb.connect()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-10-08T01:05:05.119214Z","iopub.execute_input":"2023-10-08T01:05:05.119922Z","iopub.status.idle":"2023-10-08T01:05:05.154877Z","shell.execute_reply.started":"2023-10-08T01:05:05.119891Z","shell.execute_reply":"2023-10-08T01:05:05.153473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"con.execute(f\"\"\"\n CREATE TABLE train_series AS\n SELECT * \n FROM read_parquet('{parquet_table}');\n\"\"\")\n\ncon.execute(f\"\"\"\n CREATE TABLE test_series AS\n SELECT * \n FROM read_parquet('{testing_parquet_table}');\n\"\"\")\n\n# we intentionally load every col as VARCHAR; we need to load timestamp as a varchar, to prevent the UTC time casting at load\ncon.execute(f\"\"\"\n CREATE TABLE train_events AS\n SELECT * \n FROM read_csv('{training_events_table}', ALL_VARCHAR=1, AUTO_DETECT=TRUE);\n\"\"\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-10-08T01:05:05.156234Z","iopub.execute_input":"2023-10-08T01:05:05.156526Z","iopub.status.idle":"2023-10-08T01:05:33.201656Z","shell.execute_reply.started":"2023-10-08T01:05:05.156501Z","shell.execute_reply":"2023-10-08T01:05:33.200399Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# inspired by \n# - https://www.kaggle.com/code/mitsuyasuhoshino/cmi-detectsleepstates-eda-visualization\n# - https://www.kaggle.com/code/carlmcbrideellis/zzzs-eda-utc-or-not-utc\ndef get_series_data(series_id):\n \n train_series_sql = f\"\"\"\n SELECT\n series_id,\n CAST (step AS INTEGER) AS step,\n timestamp,\n CAST (anglez AS FLOAT) AS anglez,\n CAST (enmo AS FLOAT) AS enmo\n FROM train_series WHERE series_id = '{series_id}'\n \"\"\"\n df = con.execute(train_series_sql).fetchdf()\n \n #df['date_time'] = pd.to_datetime(df['timestamp'], utc=True)\n df['date_time'] = pd.to_datetime(df['timestamp']).apply(lambda t: t.tz_localize(None))\n df = df.drop('timestamp', axis=1)\n \n train_events_sql = f\"\"\"\n SELECT \n series_id,\n CAST (night as INTEGER) as night,\n event,\n CAST (step as INTEGER) as step,\n timestamp\n FROM train_events \n WHERE series_id = '{series_id}'\n \"\"\"\n \n ev = con.execute(train_events_sql).fetchdf()\n #ev['date_time'] = pd.to_datetime(ev['timestamp'], utc=True)\n ev['date_time'] = pd.to_datetime(ev['timestamp']).apply(lambda t: t.tz_localize(None))\n ev = ev.drop('timestamp', axis=1)\n \n df2 = pd.merge(df, ev, on='date_time', how='outer')\n df2 = df2.drop(['series_id_x','series_id_y','step_x','step_y'], axis=1)\n df2['event'] = df2['event'].fillna(0)\n \n # fill values for night feature\n df2['night'][0] = 1\n df2['night'] = df2['night'].interpolate()\n df2['night'] = round(df2['night'])\n # map numeric values to event feature\n df2.loc[df2['event']=='onset','event'] = -1\n df2.loc[df2['event']=='wakeup','event'] = 1\n \n return df2\n\ndef get_test_series_data(series_id):\n test_series_sql = f\"\"\"\n SELECT\n series_id,\n CAST (step AS INTEGER) AS step,\n timestamp,\n CAST (anglez AS FLOAT) AS anglez,\n CAST (enmo AS FLOAT) AS enmo\n FROM test_series WHERE series_id = '{series_id}'\n \"\"\"\n df = con.execute(test_series_sql).fetchdf()\n \n #df['date_time'] = pd.to_datetime(df['timestamp'], utc=True)\n df['date_time'] = pd.to_datetime(df['timestamp']).apply(lambda t: t.tz_localize(None))\n df = df.drop('timestamp', axis=1)\n \n return df\n\n# inspired by https://www.kaggle.com/code/mitsuyasuhoshino/cmi-detectsleepstates-eda-visualization\ndef visualize_nights_ts(series_id):\n series_data_df = get_series_data(series_id)\n for i in set(series_data_df['night']):\n df22 = series_data_df[series_data_df['night']==i]\n df22['anglez'] = df22['anglez'] / max(abs(max(df22['anglez'])),abs(min(df22['anglez'])))\n df22['enmo'] = df22['enmo'] / max(abs(max(df22['enmo'])),abs(min(df22['enmo'])))\n df22['event'] = df22['event'] / max(abs(max(df22['event'])),abs(min(df22['event'])),1e-10)\n if df22['event'].max() == 1:\n display(Markdown('### Control chart: ' + series_id + ', night ' + str(int(i)) ))\n fig = plt.figure(figsize=(10,5))\n ax = fig.add_subplot()\n ax.plot(df22['date_time'], df22['anglez'], color='k', linewidth=0.5, label=\"anglez\")\n ax.plot(df22['date_time'], df22['enmo'], color='b', linewidth=0.5, label=\"enmo\")\n ax.plot(df22['date_time'], df22['event'], color='r', label=\"event\")\n ax.legend()\n ax.set_title(\"\".join(['Series Id : ', series_id, ' [Night# ', str(int(i)), ']']))\n ax.set_xlabel('DateTime')\n ax.set_ylabel('Normalized anglez, enmo, event [-]')\n ax.xaxis.set_major_locator(mdates.HourLocator(byhour=range(0, 24, 3)))\n ax.xaxis.set_minor_locator(mdates.HourLocator(byhour=range(0, 24, 1)))\n ax.xaxis.set_major_formatter(mdates.DateFormatter(\"%Y-%m-%d\\n%H:%M:%S\"))\n labels = ax.get_xticklabels()\n plt.setp(labels, rotation=40, fontsize=8)\n plt.show()\n \n \ndef plot_dict_nights_ts(series_id, series_data_df, night_days):\n df22 = series_data_df[series_data_df['night']==night_days]\n df22['anglez'] = df22['anglez'] / max(abs(max(df22['anglez'])),abs(min(df22['anglez'])))\n df22['enmo'] = df22['enmo'] / max(abs(max(df22['enmo'])),abs(min(df22['enmo'])))\n df22['event'] = df22['event'] / max(abs(max(df22['event'])),abs(min(df22['event'])),1e-10)\n \n plot_data_anglez = []\n plot_data_enmo = []\n plot_data_event = []\n plot_dict_anglez = dict()\n plot_dict_enmo = dict()\n plot_dict_event = dict()\n \n if df22['event'].max() == 1:\n plot_data_anglez = [df22['date_time'], df22['anglez']]\n plot_data_enmo = [df22['date_time'], df22['enmo']]\n plot_data_event = [df22['date_time'], df22['event']]\n plot_dict_anglez = dict(color='k', linewidth=0.5, label=\"anglez\")\n plot_dict_enmo = dict(color='b', linewidth=0.5, label=\"enmo\")\n plot_dict_event = dict(color='r', label=\"event\")\n \n return plot_dict_anglez, plot_dict_enmo, plot_dict_event\n ","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-10-08T01:05:33.202833Z","iopub.execute_input":"2023-10-08T01:05:33.203101Z","iopub.status.idle":"2023-10-08T01:05:33.219699Z","shell.execute_reply.started":"2023-10-08T01:05:33.203081Z","shell.execute_reply":"2023-10-08T01:05:33.218441Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# series_id = '0a96f4993bd7'\n# series_id = '08db4255286f'\n\n# \"We see the following 37 series have no NaN values whatsoever:\"\nsids = [\n '08db4255286f',\n '0a96f4993bd7',\n '0cfc06c129cc',\n '1087d7b0ff2e',\n '10f8bc1f7b07',\n '18b61dd5aae8',\n '29c75c018220',\n '31011ade7c0a',\n '3452b878e596',\n '349c5562ee2c',\n '3664fe9233f9',\n '483d6545417f',\n '55a47ff9dc8a',\n '5acc9d63b5fd',\n '5f94bb3e1bed',\n '655f19eabf1e',\n '67f5fc60e494',\n '72bbd1ac3edf',\n '76237b9406d5',\n '7822ee8fe3ec',\n '89bd631d1769',\n '8e32047cbc1f',\n '939932f1822d',\n '9ee455e4770d',\n 'a596ad0b82aa',\n 'a9a2f7fac455',\n 'a9e5f5314bcb',\n 'af91d9a50547',\n 'b364205aba43',\n 'c535634d7dcd',\n 'c6788e579967',\n 'c68260cc9e8f',\n 'ca730dbf521d',\n 'd150801f3145',\n 'd25e479ecbb7',\n 'd515236bdeec',\n 'd5e47b94477e'\n]\n# Borrowed from https://www.kaggle.com/code/carlmcbrideellis/zzzs-make-small-starter-datasets-target?scriptVersionId=143680507&cellId=13\n\n\nseries_id = random.choice(sids)\nprint(\"series_id: \", series_id)\nvisualize_nights_ts(series_id)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-10-08T01:39:10.214008Z","iopub.execute_input":"2023-10-08T01:39:10.214384Z","iopub.status.idle":"2023-10-08T01:39:25.281723Z","shell.execute_reply.started":"2023-10-08T01:39:10.214356Z","shell.execute_reply":"2023-10-08T01:39:25.280247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# inspired by https://www.kaggle.com/code/andradaolteanu/zzz-in-depth-explore-with-spectrograms-lines#3.-Spectrogram-analysis\ndef create_spectrogram(value):\n \n # defines the width of each chunk in terms of samples\n # large window size -> better frequency resolution -> poor time localization\n window_size = 256 \n # ensures that each chunk has a certain number of samples \n # that are overlapping for each chunk\n # default is 0.5*window\n overlap = 128\n # how many FFT points are desired to be computed per chunk\n # how fine-grained the frequency resolution will be\n nfft = 256 \n # sampling frequency of your signal\n # default is 1\n fs = 1.0 \n\n frequencies, times, Sxx = spectrogram(value, fs=fs, nperseg=window_size, \n noverlap=overlap, nfft=nfft)\n return frequencies, times, Sxx\n\n\n# def build_spectrogram(series_id):\n \n# series_data_df = get_series_data(series_id)\n# for i in set(series_data_df['night']):\n# spec_data = series_data_df[series_data_df['night']==i]\n# display(Markdown('### Spectrogram: ' + series_id + ', night ' + str(int(i)) ))\n \n# fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(24, 10))\n# axs = [ax1, ax2]\n# fig.suptitle(f'A night for id: {series_id}')\n \n# k=0\n\n# feat = 'anglez'\n# f, t, s = create_spectrogram(spec_data[feat].values)\n# axs[k].pcolormesh(t, f, 10 * np.log10(s), cmap=\"RdYlBu\")\n# axs[k].set_xlabel('Time (s)')\n# axs[k].set_ylabel('Frequency (Hz)')\n# axs[k].set_title(f'Spectrogram of {feat}')\n# k = k+1\n \n# feat = 'enmo'\n# f, t, s = create_spectrogram(spec_data[feat].values)\n# axs[k].pcolormesh(t, f, 10 * np.log10(s), cmap=\"RdYlBu\")\n# axs[k].set_xlabel('Time (s)')\n# axs[k].set_ylabel('Frequency (Hz)')\n# axs[k].set_title(f'Spectrogram of {feat}')\n \n# plt.show()\n \n \nclass BuildSpectrogramWithTs:\n def __init__(self):\n self.plot_data_anglez = []\n self.plot_data_enmo = []\n self.plot_data_event = []\n self.plot_dict_anglez = dict()\n self.plot_dict_enmo = dict()\n self.plot_dict_event = dict()\n \n \n def __call__(self, series_id):\n series_data_df = get_series_data(series_id)\n for i in set(series_data_df['night']):\n spec_data = series_data_df[series_data_df['night']==i]\n display(Markdown('### Spectrogram: ' + series_id + ', night ' + str(int(i)) ))\n\n # ts\n self._ts_night(spec_data, i)\n\n fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(24, 10))\n axs = [ax1, ax2]\n fig.suptitle(f'A night for id: {series_id}')\n\n k=0\n\n feat = 'anglez'\n f, t, s = create_spectrogram(spec_data[feat].values)\n axs[k].pcolormesh(t, f, 10 * np.log10(s), cmap=\"RdYlBu\")\n axs[k].set_xlabel('Time (s)')\n axs[k].set_ylabel('Frequency (Hz)')\n axs[k].set_title(f'Spectrogram of {feat}')\n # ts\n axs[k].plot(*self.plot_data_anglez, **self.plot_dict_anglez)\n axs[k].plot(*self.plot_data_event, **self.plot_dict_event)\n axs[k].legend()\n k = k+1\n\n feat = 'enmo'\n f, t, s = create_spectrogram(spec_data[feat].values)\n axs[k].pcolormesh(t, f, 10 * np.log10(s), cmap=\"RdYlBu\")\n axs[k].set_xlabel('Time (s)')\n axs[k].set_ylabel('Frequency (Hz)')\n axs[k].set_title(f'Spectrogram of {feat}')\n # ts\n axs[k].plot(*self.plot_data_enmo, **self.plot_dict_enmo)\n axs[k].plot(*self.plot_data_event, **self.plot_dict_event)\n axs[k].legend()\n\n plt.show()\n \n \n def _ts_night(self, df22, night_days):\n df22 = df22.reset_index(drop=True)\n df22['anglez'] = df22['anglez'] / max(abs(max(df22['anglez'])),abs(min(df22['anglez'])))\n df22['enmo'] = df22['enmo'] / max(abs(max(df22['enmo'])),abs(min(df22['enmo'])))\n df22['event'] = df22['event'] / max(abs(max(df22['event'])),abs(min(df22['event'])),1e-10)\n \n if df22['event'].max() == 1:\n \n # Overlap = 128\n df22 = df22[['anglez', 'enmo', 'event']]\n df22 = df22.rolling(window=128, step=128).sum().dropna()[:-1]\n \n self.plot_data_anglez = [df22.index, self._normalize(df22['anglez'])]\n self.plot_data_enmo = [df22.index, self._normalize(df22['enmo'])]\n self.plot_data_event = [df22.index, self._normalize(df22['event'])]\n self.plot_dict_anglez = dict(color='k', linewidth=0.5, label=\"anglez\")\n self.plot_dict_enmo = dict(color='b', linewidth=0.5, label=\"enmo\")\n self.plot_dict_event = dict(color='r', label=\"event\")\n \n def _normalize(self, lst):\n min_old = min(lst)\n max_old = max(lst)\n min_new = 0\n max_new = 0.5\n normalized_lst = [(min_new + ((value - min_old) * (max_new - min_new) / (max_old - min_old))) for value in lst]\n return normalized_lst","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-10-08T01:22:58.450956Z","iopub.execute_input":"2023-10-08T01:22:58.451320Z","iopub.status.idle":"2023-10-08T01:22:58.467817Z","shell.execute_reply.started":"2023-10-08T01:22:58.451294Z","shell.execute_reply":"2023-10-08T01:22:58.466105Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# build_spectrogram(series_id)\n\nbuild_spectrogram = BuildSpectrogramWithTs()\nbuild_spectrogram(series_id)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-10-08T01:39:36.440678Z","iopub.execute_input":"2023-10-08T01:39:36.441038Z","iopub.status.idle":"2023-10-08T01:39:56.006928Z","shell.execute_reply.started":"2023-10-08T01:39:36.441011Z","shell.execute_reply":"2023-10-08T01:39:56.005491Z"},"trusted":true},"execution_count":null,"outputs":[]}]}