{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":53666,"databundleVersionId":6589269,"sourceType":"competition"}],"dockerImageVersionId":30558,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# データセットの作成","metadata":{}},{"cell_type":"markdown","source":"## ライブラリ","metadata":{}},{"cell_type":"code","source":"! pip install japanize_matplotlib --quiet","metadata":{"execution":{"iopub.status.busy":"2023-11-27T06:36:09.518887Z","iopub.execute_input":"2023-11-27T06:36:09.519296Z","iopub.status.idle":"2023-11-27T06:36:27.853506Z","shell.execute_reply.started":"2023-11-27T06:36:09.519262Z","shell.execute_reply":"2023-11-27T06:36:27.852172Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport pyarrow as pa\nimport pyarrow.parquet as pq\nimport matplotlib.pyplot as plt\nimport japanize_matplotlib\nimport os\nimport gc\nimport time\nimport IPython.display\nimport pandas.tseries.offsets as offsets\n\nfrom scipy.stats import norm\nfrom pathlib import Path\nfrom tqdm import tqdm","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-11-27T06:36:27.855555Z","iopub.execute_input":"2023-11-27T06:36:27.855900Z","iopub.status.idle":"2023-11-27T06:36:28.875203Z","shell.execute_reply.started":"2023-11-27T06:36:27.855872Z","shell.execute_reply":"2023-11-27T06:36:28.873552Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"HOME = Path('/kaggle')\nORIGIN = HOME / 'input/child-mind-institute-detect-sleep-states'\nWORKING = HOME / 'working'\n\nEXP = ['a_b_var', 'a_a_var', 'e_b_mean', 'e_a_mean', 'pd']\nOBJ = ['state']\n\nSTEP_OF_WARMUP = 12 * 60 * 2\nSTEP_OF_RATE = 3\nSTEP_OF_STRIDE = 12\nSTEP_OF_MOVING = 12 * 60 * 6\n\nSEC_OF_DAY  = 60 * 60 * 24\nSEC_OF_HOUR = 60 * 60\n\nSTEP_OF_DAY   = 12 * 60 * 24\nSTEP_OF_HOUR  = 12 * 60\nSTEP_OF_30MIN = 12 * 30\n\nEVENT_MAP = {'onset':0, 'wakeup':1}","metadata":{"execution":{"iopub.status.busy":"2023-11-27T06:36:28.876507Z","iopub.execute_input":"2023-11-27T06:36:28.876958Z","iopub.status.idle":"2023-11-27T06:36:28.884689Z","shell.execute_reply.started":"2023-11-27T06:36:28.876927Z","shell.execute_reply":"2023-11-27T06:36:28.882920Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def graphs(df, id_map, start=0, end=0):\n    if not end:\n        df = df[df['id_map']==id_map].copy().reset_index(drop=True)\n    else:\n        df = df[df['id_map']==id_map].copy().reset_index(drop=True).iloc[start:end]\n    steps = df[df['time'].apply(lambda x: str(x)[-8:]=='00:00:00')]['step']\n    flag = len(df) != df['step'].iat[-1]\n\n    _, axes = plt.subplots(2, 1, figsize=(20,6))\n    axes[0].set_title(id_map)\n    for ax, col in zip(axes.ravel(), ['anglez', 'enmo']):\n        if 'state' in df.columns:\n            if col == 'anglez':\n                df_state = df['state'].apply(lambda x: 50 if x==1 else -50 if x==0 else 0)\n            else:\n                df_state = df['state'].apply(lambda x: 2 if x==1 else 0 if x==0 else 1)\n        X = df.index if flag else df['step']\n        ax.plot(X, df[col], linewidth=0.5, marker='', color='0.5', label=col)\n        if 'state' in df.columns:\n            ax.plot(X, df_state, linewidth=1, marker='', color='b', label='state')\n        for step in steps:\n            X = df[df['step']==step].index[0] if flag else step\n            ax.axvline(x=X, linewidth=1, color='r', alpha=0.5)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-11-27T06:36:28.888354Z","iopub.execute_input":"2023-11-27T06:36:28.891685Z","iopub.status.idle":"2023-11-27T06:36:28.906153Z","shell.execute_reply.started":"2023-11-27T06:36:28.891615Z","shell.execute_reply":"2023-11-27T06:36:28.903578Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 分析用データ成形","metadata":{}},{"cell_type":"code","source":"%%time\n\ndf = pq.read_table(ORIGIN / 'train_series.parquet').to_pandas()\ndf_test = pq.read_table(ORIGIN / 'test_series.parquet').to_pandas()\ndf_event = pd.read_csv(ORIGIN / 'train_events.csv')\ndf_sample = pd.read_csv(ORIGIN / 'sample_submission.csv', index_col=0)","metadata":{"execution":{"iopub.status.busy":"2023-11-27T06:36:28.907859Z","iopub.execute_input":"2023-11-27T06:36:28.909356Z","iopub.status.idle":"2023-11-27T06:37:33.196943Z","shell.execute_reply.started":"2023-11-27T06:36:28.909187Z","shell.execute_reply":"2023-11-27T06:37:33.195440Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_metric = [12, 36, 60, 90, 120, 150, 180, 240, 300, 360]\ndf_metric = pd.DataFrame([[y, y/12, (10-x)/10] for x,y in enumerate(df_metric)], columns=['step_under','minutes_under','metric_score'])","metadata":{"execution":{"iopub.status.busy":"2023-11-27T06:37:33.198202Z","iopub.execute_input":"2023-11-27T06:37:33.198600Z","iopub.status.idle":"2023-11-27T06:37:33.211846Z","shell.execute_reply.started":"2023-11-27T06:37:33.198569Z","shell.execute_reply":"2023-11-27T06:37:33.207538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.info()","metadata":{"execution":{"iopub.status.busy":"2023-11-27T06:37:33.213752Z","iopub.execute_input":"2023-11-27T06:37:33.214188Z","iopub.status.idle":"2023-11-27T06:37:33.248096Z","shell.execute_reply.started":"2023-11-27T06:37:33.214152Z","shell.execute_reply":"2023-11-27T06:37:33.244996Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\ndf_id = pd.DataFrame([[y, x] for x, y in enumerate(df['series_id'].unique())], columns=['series_id', 'id_map'])\ndf_id['id_map'] = df_id['id_map'].astype('int16')","metadata":{"execution":{"iopub.status.busy":"2023-11-27T06:37:33.250371Z","iopub.execute_input":"2023-11-27T06:37:33.250743Z","iopub.status.idle":"2023-11-27T06:37:41.613967Z","shell.execute_reply.started":"2023-11-27T06:37:33.250713Z","shell.execute_reply":"2023-11-27T06:37:41.612473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_id","metadata":{"execution":{"iopub.status.busy":"2023-11-27T06:37:41.616521Z","iopub.execute_input":"2023-11-27T06:37:41.617258Z","iopub.status.idle":"2023-11-27T06:37:41.642438Z","shell.execute_reply.started":"2023-11-27T06:37:41.617198Z","shell.execute_reply":"2023-11-27T06:37:41.640304Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def data_pre(df):\n    df['time'] = df['timestamp'].apply(lambda x: x[:-5]).astype('datetime64[ns]')\n    df = df.merge(df_id, on='series_id', how=\"left\").drop(['series_id', 'timestamp'], axis=1)\n    df.info()\n    return df","metadata":{"execution":{"iopub.status.busy":"2023-11-27T06:37:41.649440Z","iopub.execute_input":"2023-11-27T06:37:41.649809Z","iopub.status.idle":"2023-11-27T06:37:41.658658Z","shell.execute_reply.started":"2023-11-27T06:37:41.649769Z","shell.execute_reply":"2023-11-27T06:37:41.656452Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\ndf = data_pre(df)\nprint()","metadata":{"execution":{"iopub.status.busy":"2023-11-27T06:37:41.660182Z","iopub.execute_input":"2023-11-27T06:37:41.660507Z","iopub.status.idle":"2023-11-27T06:39:18.000460Z","shell.execute_reply.started":"2023-11-27T06:37:41.660484Z","shell.execute_reply":"2023-11-27T06:39:17.998466Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.head()","metadata":{"execution":{"iopub.status.busy":"2023-11-27T06:39:18.002459Z","iopub.execute_input":"2023-11-27T06:39:18.002830Z","iopub.status.idle":"2023-11-27T06:39:18.025919Z","shell.execute_reply.started":"2023-11-27T06:39:18.002806Z","shell.execute_reply":"2023-11-27T06:39:18.024824Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\n# 欠損値は就寝/起床をセットで扱う\nnull_list = [558, 1425, 3879, 13553, 14105]\ndf_event.loc[df_event.index.isin(null_list), 'step'] = np.nan\ndf_event.loc[df_event.index.isin(null_list), 'timestamp'] = np.nan\n\ndf_event.loc[pd.isna(df_event['timestamp']), 'step'] = -1\ndf_event.loc[pd.isna(df_event['timestamp']), 'timestamp'] = '1999-12-31T23:59:59-0400'\ndf_event['event'] = df_event['event'].replace(EVENT_MAP).astype('int8')\ndf_event['step'] = df_event['step'].astype('int32')\ndf_event = data_pre(df_event)\nprint()","metadata":{"execution":{"iopub.status.busy":"2023-11-27T06:39:18.028866Z","iopub.execute_input":"2023-11-27T06:39:18.029261Z","iopub.status.idle":"2023-11-27T06:39:18.080144Z","shell.execute_reply.started":"2023-11-27T06:39:18.029230Z","shell.execute_reply":"2023-11-27T06:39:18.079075Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\n# dfにeventカラムを追加：'onset':0, 'wakeup':1, 'other':2\ndf_tmp = df_event[['id_map', 'step', 'event']]\ndf = pd.merge(df, df_tmp, on=['id_map', 'step'], how='left').fillna(2)\ndf['event'] = df['event'].astype('int8')\ndf.info()","metadata":{"execution":{"iopub.status.busy":"2023-11-27T06:39:18.081720Z","iopub.execute_input":"2023-11-27T06:39:18.083401Z","iopub.status.idle":"2023-11-27T06:40:11.302039Z","shell.execute_reply.started":"2023-11-27T06:39:18.083362Z","shell.execute_reply":"2023-11-27T06:40:11.300534Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def moving_var(x_array, term=None, w=STEP_OF_MOVING):\n    m_init = np.mean(x_array[:w]); mean_list = [m_init] #平均の初期値\n    v_init = np.var(x_array[:w]); var_list = [v_init] #分散の初期値\n    for i in range(len(x_array)-w): # 更新式\n        mean_list.append(mean_list[-1] + (-x_array[i] +x_array[i+w])/w)\n        var_list.append(var_list[-1]+mean_list[-2]**2 +(x_array[i+w]**2 - x_array[i]**2)/w  -mean_list[-1]**2)\n    if term == 'before':\n        mean_list = mean_list + [mean_list[-1] for _ in range(w-1)]        \n        var_list = var_list + [var_list[-1] for _ in range(w-1)]\n    if term == 'after':\n        mean_list = [mean_list[0] for _ in range(w-1)] + mean_list        \n        var_list = [var_list[0] for _ in range(w-1)] + var_list\n    return mean_list, var_list","metadata":{"execution":{"iopub.status.busy":"2023-11-27T06:40:11.304004Z","iopub.execute_input":"2023-11-27T06:40:11.304414Z","iopub.status.idle":"2023-11-27T06:40:11.315071Z","shell.execute_reply.started":"2023-11-27T06:40:11.304385Z","shell.execute_reply":"2023-11-27T06:40:11.313464Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"a_b_var, a_a_var = np.array([]), np.array([])\ne_b_mean, e_a_mean = np.array([]), np.array([])\nwith tqdm(df_id['id_map']) as t:\n    for id_map in t:\n        a = df[df['id_map']==id_map]['anglez'].values\n        e = df[df['id_map']==id_map]['enmo'].values                       \n        a_b_var = np.concatenate([a_b_var, moving_var(a, 'before')[1]])\n        a_a_var = np.concatenate([a_a_var, moving_var(a, 'after')[1]])\n        e_b_mean = np.concatenate([e_b_mean, moving_var(e, 'before')[0]])\n        e_a_mean = np.concatenate([e_a_mean, moving_var(e, 'after')[0]])\ndf = df.assign(a_b_var=a_b_var)\ndf = df.assign(a_a_var=a_a_var)\ndf = df.assign(e_b_mean=e_b_mean)\ndf = df.assign(e_a_mean=e_a_mean)\ndel a_b_var, a_a_var, e_b_mean, e_a_mean, a, e; gc.collect()\nprint()\ndf.info()","metadata":{"execution":{"iopub.status.busy":"2023-11-27T06:40:11.317681Z","iopub.execute_input":"2023-11-27T06:40:11.318073Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def time_to_sec(x):\n    return x.timestamp() % SEC_OF_DAY\n\ndef time_to_step(x):\n    return (x.timestamp() % SEC_OF_DAY) / 5","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n\ndf_event['sec'] = -1\ndf_event.loc[(df_event['event']==0)&(df_event['step']!=-1), 'sec'] = df_event['time'].apply(lambda x: time_to_sec(x) + (0 if x.hour>12 else SEC_OF_DAY))\ndf_event.loc[(df_event['event']==1)&(df_event['step']!=-1), 'sec'] = df_event['time'].apply(time_to_sec)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_event.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 学習データ全体の統計情報を df_stat に集約\ndf_stat = {}\n\ndf_tmp = df_event[(df_event['step']!=-1) & df_event['event']]['sec']\ndf_stat['w_mean'] = round(df_tmp.mean() / 5)\ndf_stat['w_std']  = round(df_tmp.std() / 5)\ndf_stat['w_var']  = round(df_tmp.var() / 25)\ndf_tmp = df_event[(df_event['step']!=-1) & ~df_event['event']]['sec']\ndf_stat['o_mean'] = round(df_tmp.mean() / 5)\ndf_stat['o_std']  = round(df_tmp.std() / 5)\ndf_stat['o_var']  = round(df_tmp.var() / 25)\n\nprint('average')\nprint(f\"wakeup:  {df_stat['w_mean']//STEP_OF_HOUR}:{round((df_stat['w_mean']%STEP_OF_HOUR)/12)} (std {df_stat['w_std']//12}min)\")\nprint(f\"onset : {df_stat['o_mean']//STEP_OF_HOUR}:{round((df_stat['o_mean']%STEP_OF_HOUR)/12)} (std {df_stat['o_std']//12}min)\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"## ID別の欠損値を含む期間 length と欠損値を除いた期間 night\ndf_id = df_id.assign(length=[len(df_event[df_event['id_map']==x])/2 for x in df_event['id_map'].unique()])\ndf_id = df_id.assign(night=[len(df_event[(df_event['id_map']==x) & (~df_event['step'])])/2 for x in df_event['id_map'].unique()])\ndf_id['length'] = df_id['length'].astype('int8')\ndf_id['night'] = df_id['night'].astype('int8')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 学習データ全体の length と night の平均\ndf_stat['l_mean'] = round(df_id['length'].mean())\ndf_stat['n_mean'] = round(df_id['night'].mean())\ndf_stat","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def calc_stat(df_id, df_event=df_event, df_stat=df_stat):\n    for event in ['wakeup', 'onset']:\n        e = event[0]\n        v = df_stat[f'{e}_var'] * (df_stat['n_mean']-1)\n        dic = {'mean':f'{e}_mean', 'std':f'{e}_std'}\n        if event == 'wakeup':\n            df_tmp = df_event[(df_event['step']!=-1) & df_event['event']][['id_map','sec']]\n        elif event == 'onset':\n            df_tmp = df_event[(df_event['step']!=-1) & ~df_event['event']][['id_map','sec']]\n        n = len(df_tmp) - 1\n        df_tmp = df_tmp.groupby('id_map').agg(['mean', 'var', 'count'])['sec']\n        df_tmp['mean'] = df_tmp['mean'] / 5\n        df_tmp['var'] = df_tmp['var'] / 25\n        df_tmp['std'] = ((df_tmp['var'] * (df_tmp['count']-1) + v) / (df_tmp['count'] + df_stat['n_mean'] - 1)) ** 0.5\n        df_tmp = df_tmp.drop(['var','count'], axis=1).apply(round).rename(columns=dic)\n        df_id = df_id.merge(df_tmp, on='id_map', how='left').fillna(-1)\n        df_id.loc[df_id['night']==0, f'{e}_mean'] = df_stat[f'{e}_mean']\n        df_id.loc[df_id['night']<=1, f'{e}_std'] = df_id[f'{e}_std'].max()\n        df_id[f'{e}_mean'] = df_id[f'{e}_mean'].astype('int32')\n        df_id[f'{e}_std'] = df_id[f'{e}_std'].astype('int16')\n    return df_id","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ID別に状態変化の確率密度を求めるために必要な統計量を求める\ndf_id = calc_stat(df_id)\ndf_id.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 24時間周期で anglez と enmo が全く同じ値を取る区間を特定する。\n# 正しく記録できていないケースが多い。\n\ndef null_range(id_map, df):\n    r = []\n    if len(df) < STEP_OF_DAY:\n        return r\n    df['a_s'] = df['anglez'].shift(STEP_OF_DAY)\n    df['e_s'] = df['enmo'].shift(STEP_OF_DAY)\n    df['flag'] = (df['anglez']==df['a_s']) & (df['enmo']==df['e_s'])\n    running = df['flag'].groupby(df['flag'].diff().ne(0).cumsum()).cumsum()\n    i = len(df)-1\n    if running[i] >= STEP_OF_DAY:\n        s = running[-running[i]:].index[0]\n        e = i + 1\n        r.append(range(s,e))\n        i = s - 1\n    while i >= 0:\n        if running[i] < STEP_OF_DAY:\n            i -= 1\n        else:\n            e = i + 1\n            s = e - STEP_OF_DAY + 1\n            r.append(range(s,e))\n            i = s - 1\n    return r","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ID別の24時間周期の anglez と enmo の平均/移動分散を求めて、df_mean に集約\nwith tqdm(df_id['id_map']) as t:  \n    for id_map in t:\n        r_list = null_range(id_map, df[df['id_map']==id_map][['anglez','enmo']].reset_index(drop=True))\n        df_tmp = df[df['id_map']==id_map]\n        for r in r_list:\n            df_tmp = df_tmp[~df_tmp['step'].isin(r)]\n        df_tmp['time'] = df_tmp['time'].apply(time_to_step).astype('int16')\n        IPython.display.clear_output()\n        df_tmp = df_tmp[['time','id_map','anglez','enmo','a_b_var','a_a_var','e_b_mean','e_a_mean']].groupby('time').mean()\n        df_mean = pd.concat([df_mean, df_tmp]) if id_map else df_tmp\n    df_mean['id_map'] = df_mean['id_map'].astype('int16')\n    df_mean = df_mean.reset_index()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def pdf_handle(id_map, state='wakeup', df_id=df_id): # PDF: Probability Density Function\n    mean = df_id.loc[df_id['id_map']==id_map, f'{state[0]}_mean'].iat[0]\n    std = df_id.loc[df_id['id_map']==id_map, f'{state[0]}_std'].iat[0]\n    pd = norm.pdf(np.arange(STEP_OF_DAY * 3), mean, std)\n    pd = pd[:STEP_OF_DAY] + pd[STEP_OF_DAY:STEP_OF_DAY*2] + pd[STEP_OF_DAY*2:]\n    return pd","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# ID別の24時間周期のstepに対する各状態（wakeup, onset）の確率密度\nwith tqdm(df_id['id_map']) as t:\n    w_pd, o_pd = np.array([]), np.array([])\n    for i in t:\n        w_pd = np.concatenate([w_pd, pdf_handle(i, 'wakeup')])\n        o_pd = np.concatenate([o_pd, pdf_handle(i, 'onset')])    \n    df_mean = df_mean.assign(w_pd=w_pd)\n    df_mean = df_mean.assign(o_pd=o_pd)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"state = np.array([])\nend = df_id.index[-1]\nfor i in df_id.index:\n    l = 0\n    if df_id['o_mean'].iat[i] > STEP_OF_DAY:\n        l = df_id['o_mean'].iat[i] - STEP_OF_DAY\n        state = np.concatenate([state, np.ones(l)])\n    l = df_id['w_mean'].iat[i] - l\n    state = np.concatenate([state, np.zeros(l)])\n    if df_id['o_mean'].iat[i] < STEP_OF_DAY:\n        l = df_id['o_mean'].iat[i] - df_id['w_mean'].iat[i]\n    else:\n        l = STEP_OF_DAY - df_id['w_mean'].iat[i]\n    state = np.concatenate([state, np.ones(l)])\n    if df_id['o_mean'].iat[i] < STEP_OF_DAY:\n        l = STEP_OF_DAY - df_id['o_mean'].iat[i]\n        state = np.concatenate([state, np.zeros(l)])\ndf_mean = df_mean.assign(state=state)\ndf_mean['state'] = df_mean['state'].astype('int8')\ndf_mean['event'] = df_mean['state'] - df_mean['state'].shift(1).fillna(0)\ndf_mean['event'] = df_mean['event'].apply(lambda x: 0 if x==-1 else 1 if x==1 else 2).astype('int8')\ndf_mean.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 状態'state'カラムの追加\n# 欠損値期間は直前の起床後1時間から+24hの期間: state -1\n\ndef state_handle(id_map, length, df):\n    if len(df) == sum(df['step']==-1):\n        return np.full(length, -1)\n    \n    state, s = np.array([]), 1\n    start, count = 0, 0\n    \n    for i in range(len(df)):\n        if df['step'].iat[i] == -1:\n            if i == len(df)-1: # 最後が欠損区間の場合\n                end = min(length, start+STEP_OF_HOUR)\n                l = end - start\n                state = np.concatenate([state, np.full(l, s)])\n                l = length - end\n                state = np.concatenate([state, np.full(l, -1)])\n            else:\n                count += 0.5\n        else:\n            end = df['step'].iat[i]\n            if count: # 欠損値区間の直後\n                if start != 0:\n                    l = STEP_OF_HOUR\n                    state = np.concatenate([state, np.full(l, s)])\n                    l = STEP_OF_DAY * int(count)\n                    state = np.concatenate([state, np.full(l, -1)])\n                    l = end - len(state)\n                    state = np.concatenate([state, np.full(l, s)])\n                else: # 最初から欠損区間の場合\n                    l = max(0, end-STEP_OF_HOUR*6)\n                    state = np.concatenate([state, np.full(l, -1)])\n                    l = min(end, STEP_OF_HOUR*6)\n                    state = np.concatenate([state, np.full(l, s)])\n            else:\n                l = end - start\n                state = np.concatenate([state, np.full(l, s)])\n            s = df['event'].iat[i]\n            start = end\n            if i == len(df) - 1:\n                l = min(length-end, STEP_OF_HOUR * 6)\n                state = np.concatenate([state, np.full(l, s)])\n                if l == STEP_OF_HOUR * 6: # 末尾の余剰部分は-1\n                    l = length - end - STEP_OF_HOUR * 6\n                    state = np.concatenate([state, np.full(l, -1)])\n            count = 0\n    if len(state) != length:\n        print(f'error: {id_map}, {len(state)}, {length}')\n    return state","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"state = np.array([])\nwith tqdm(df_id['id_map']) as t:\n    for id_map in t:\n        df_tmp = df[df['id_map']==id_map][['step','event']]\n        state = np.concatenate([state, state_handle(id_map, len(df_tmp), df_event[df_event['id_map']==id_map])])\ndf = df.assign(state=state)\ndf['state'] = df['state'].astype('int8')\ndel df_tmp, state; gc.collect()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.to_parquet('custom_train_series.parquet')\ndf_event.to_csv('custom_train_events.csv')\ndf_id.to_csv('custom_train_stat_groupby_id.csv')\ndf_mean.to_csv('custom_train_series_mean.csv')\ndf_metric.to_csv('metric.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# コンペの概要","metadata":{}},{"cell_type":"markdown","source":"## 内容\n* 子供の手に付けた計測器の情報（角度と加速度）から睡眠/起床の変化点を予測","metadata":{}},{"cell_type":"markdown","source":"## 睡眠の定義  \n* 1回の睡眠時間は30分以上\n* 夜間の最も長い睡眠状態だった時間を睡眠時間とする\n* 睡眠時間として計測されるのは基本的に24時間あたり1回\n* 上記を満たせば同日で 01:00～06:00 と 19:00～23:30 が睡眠時間となるケースもある","metadata":{}},{"cell_type":"markdown","source":"## 評価指標\n* 下記 df_metric の内容にて各イベントでスコアを加算し、その平均が最終的なスコアとなる\n* submission.csvのscoreカラムの大きい順に評価されていく\n* 欠損値（正しく装着していないと思われる期間）に対してどのように予測すれば良いかが調査不足","metadata":{}},{"cell_type":"code","source":"df_metric","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 提出時のファイル名は submission.csv\ndf_sample","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# データの分析","metadata":{}},{"cell_type":"markdown","source":"## 時系列データ","metadata":{}},{"cell_type":"markdown","source":"### 時系列データの特徴 \n* データ数：1億2794万（元）\n* 被験者数：277名","metadata":{}},{"cell_type":"markdown","source":"### 時系列データのカラム情報\n* 時間　　：5秒間隔\n* ステップ：同一人物の時間の経過に対応\n* 角度　　：-90°～90°\n* 加速度　：標準化された0以上の値","metadata":{}},{"cell_type":"code","source":"df.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 変化点データ","metadata":{}},{"cell_type":"markdown","source":"### 変化点データのカラム\n* 夜の日数：同一人物にとって計測開始後何日目の夜か\n* イベント：睡眠 or 起床\n* 時間　　：イベント発生時の時間\n* ステップ：イベント発生時のステップ","metadata":{}},{"cell_type":"code","source":"df_event.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### ID単位での変化点について\n* 時系列データ自体の長さの違い、変化点データの欠損値の有無、により'night'の数はID単位で異なる","metadata":{}},{"cell_type":"code","source":"print(f\"変化点データが全て欠損値のID: {len(df_id[df_id['night']==0])}\\n\")\n_, axes = plt.subplots(1, 2, figsize=(20,6))\naxes = axes.ravel()\naxes[0].hist(df_id['length'], bins=20)\naxes[0].set_title('ID別 night の長さ（元）')\naxes[1].hist(df_id['night'], bins=20)\naxes[1].set_title('ID別 night の長さ（欠損値の削除後）')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 変化点の確率密度の例","metadata":{}},{"cell_type":"code","source":"id_map = 16\nplt.figure(figsize=(20,5))\nplt.title('24時間での各状態（wakeup, onset）の確率密度')\nplt.plot(df_mean[df_mean['id_map']==id_map].reset_index()['w_pd'], color='b')\nplt.plot(df_mean[df_mean['id_map']==id_map].reset_index()['o_pd'], color='r')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## データの加工","metadata":{}},{"cell_type":"markdown","source":"### 時系列データの不要部分の削除\n* 変化点が欠損値の期間を削除  \n* 前後の余剰部分を削除","metadata":{}},{"cell_type":"code","source":"_, axes = plt.subplots(1, 2, figsize=(20,6))\naxes = axes.ravel()\naxes[0].hist(df.groupby('id_map').max()['step'], bins=20)\naxes[0].set_title('ID別データ数（元）')\naxes[1].hist(df[df['state']==-1].groupby('id_map').count()['step'], bins=20)\naxes[1].set_title('ID別データ数（加工後）')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 時系列の元データから train/val/test の分離","metadata":{}},{"cell_type":"code","source":"def add_pd(df, id_map, df_mean=df_mean):\n    df_tmp = df_mean[df_mean['id_map']==id_map][['time', 'w_pd', 'o_pd']].rename(columns={'time':'t'})\n    df['t'] = df['time'].apply(time_to_step).astype('int64')\n    df = df.merge(df_tmp, on='t').sort_values('step').drop('t', axis=1)\n    return df\n\ndef add_warmup(df, id_map, warmup, col, df_mean=df_mean):\n    return df","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def train_test_extract(id_map, warmup, n_val=2, n_test=2, df=df, df_id=df_id):\n    df = df[(df['id_map']==id_map) & (df['state']!=-1)].copy() if 'state' in df.columns else df.copy()\n    df = add_pd(df, id_map).reset_index(drop=True)\n    \n    night   = df_id[df_id['id_map']==id_map]['night'].iat[0]\n    n_train = night - n_test - n_val\n    steps = df_event[(df_event['id_map']==id_map) & (df_event['step']!=-1)]['step']\n\n    if n_test:\n        end    = steps.iat[(n_train+n_val)*2]\n        start  = steps.iat[(n_train+n_val)*2+1]\n        s_test = (end + start) // 2\n    else:\n        s_test = df['step'].iat[-1] + 1\n        \n    if n_val:\n        end   = steps.iat[n_train*2]\n        start = steps.iat[n_train*2+1]\n        s_val = (end + start) // 2\n    else:\n        s_val = s_test\n\n    col = ['id_map', 'step', 'time', 'anglez', 'enmo', 'state']\n    df = df[col]\n\n    df_tmp = df[df['step']<s_val].reset_index(drop=True)\n    train = add_warmup(df_tmp, id_map, warmup, col)\n    df_tmp = df[(df['step']>=s_val) & (df['step']<s_test)].reset_index(drop=True)\n    val = add_warmup(df_tmp, id_map, warmup, col) if n_val else df.head(0)\n    df_tmp = df[df['step']>=s_test].reset_index(drop=True)\n    test  = add_warmup(df_tmp, id_map, warmup, col) if n_test else df.head(0)\n    \n    IPython.display.clear_output()\n    print(f'night: {night}')\n    print(f'train: {n_train}')\n    print(f'val  : {n_val}')\n    print(f'test : {n_test}')\n    \n    return train, val, test","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"id_map = 84\nwarmup = 0\ntrain, val, test = train_test_extract(id_map, warmup)\ndf_tmp = df[df['id_map']==id_map]\n\nprint('\\noriginal\\n')\ngraphs(df_tmp, id_map)\nprint('\\ncustom\\n')\ngraphs(df_tmp[df_tmp['state']!=-1], id_map)\nprint('\\ntrain\\n')\ngraphs(train, id_map)\nprint('\\nval\\n')\ngraphs(val, id_map)\nprint('\\ntest\\n')\ngraphs(test, id_map)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}