{"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":"code","source":"import itertools\nfrom pathlib import Path\nfrom sklearn.linear_model import LinearRegression\nfrom sklearn.metrics import confusion_matrix, f1_score\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport pyarrow as pa\nimport pyarrow.parquet as pq\nimport pyarrow.dataset as pads\n\nimport time\nfrom tqdm import tqdm\nfrom tqdm.contrib.concurrent import process_map\nimport matplotlib.pyplot as plt ","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-09-25T09:57:22.101564Z","iopub.execute_input":"2023-09-25T09:57:22.101945Z","iopub.status.idle":"2023-09-25T09:57:24.227468Z","shell.execute_reply.started":"2023-09-25T09:57:22.101884Z","shell.execute_reply":"2023-09-25T09:57:24.226232Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install ruptures --no-index --find-links=/kaggle/input/ruptures-offline","metadata":{"execution":{"iopub.status.busy":"2023-09-25T09:57:24.229179Z","iopub.execute_input":"2023-09-25T09:57:24.229566Z","iopub.status.idle":"2023-09-25T09:57:35.246726Z","shell.execute_reply.started":"2023-09-25T09:57:24.229543Z","shell.execute_reply":"2023-09-25T09:57:35.245303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"WINDOW_SIZE = 361","metadata":{"execution":{"iopub.status.busy":"2023-09-25T09:57:35.248366Z","iopub.execute_input":"2023-09-25T09:57:35.248750Z","iopub.status.idle":"2023-09-25T09:57:35.254459Z","shell.execute_reply.started":"2023-09-25T09:57:35.248719Z","shell.execute_reply":"2023-09-25T09:57:35.253349Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"INPUT_DIR = Path('/kaggle/input/child-mind-institute-detect-sleep-states')\nWORKING_DIR = Path('/kaggle/working/cmi-detect-sleep-states/')\n\nTEST_SERIES = INPUT_DIR / 'test_series.parquet'\n\nTEST_DATASET = WORKING_DIR / 'test'","metadata":{"execution":{"iopub.status.busy":"2023-09-25T09:57:35.258170Z","iopub.execute_input":"2023-09-25T09:57:35.258545Z","iopub.status.idle":"2023-09-25T09:57:35.270577Z","shell.execute_reply.started":"2023-09-25T09:57:35.258515Z","shell.execute_reply":"2023-09-25T09:57:35.269786Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Setup the parquet dataset partitioned by series id\npartitioning = pads.partitioning(schema=pa.schema([pa.field('series_id', pa.string())]))\nprint('Reading dataset...')\ntry:\n    pq_ds = pq.ParquetDataset(TEST_DATASET, use_legacy_dataset=False, partitioning=partitioning)\nexcept FileNotFoundError:\n    # Create it if it doesn't exist\n    print('Dataset not found, Creating dataset grouped by series_id...')\n    pq_table = pq.read_table(TEST_SERIES)\n    pq.write_to_dataset(pq_table, root_path=TEST_DATASET, partitioning=partitioning)\n    \n    # Read it again\n    pq_ds = pq.ParquetDataset(TEST_DATASET, use_legacy_dataset=False, partitioning=partitioning)\n\n# Each series can be accessed by iterating over the fragments\ndef get_series_id_from_fragment(f):\n    pkeys = pads._get_partition_keys(f.partition_expression)\n    return pkeys.get('series_id', None)","metadata":{"execution":{"iopub.status.busy":"2023-09-25T09:57:35.271426Z","iopub.execute_input":"2023-09-25T09:57:35.271714Z","iopub.status.idle":"2023-09-25T09:57:35.352404Z","shell.execute_reply.started":"2023-09-25T09:57:35.271689Z","shell.execute_reply":"2023-09-25T09:57:35.351471Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def compute_delta(df):\n    # Extract the series - enmo\n    enmo = df['enmo'].values\n    \n    # Pad input\n    enmo = np.pad(enmo, 1, mode='edge')\n    \n    # Compute central difference\n    kernel = np.array([-0.5, 0, 0.5], dtype=np.float32)\n    enmo_diff = np.convolve(enmo, kernel, 'valid')\n    \n    return enmo_diff","metadata":{"execution":{"iopub.status.busy":"2023-09-25T09:57:35.353221Z","iopub.execute_input":"2023-09-25T09:57:35.353467Z","iopub.status.idle":"2023-09-25T09:57:35.358792Z","shell.execute_reply.started":"2023-09-25T09:57:35.353445Z","shell.execute_reply":"2023-09-25T09:57:35.357996Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"min_valid_sleep_minutes = 2\nwindow = int(min_valid_sleep_minutes * 60/5)","metadata":{"execution":{"iopub.status.busy":"2023-09-25T09:57:35.360116Z","iopub.execute_input":"2023-09-25T09:57:35.360410Z","iopub.status.idle":"2023-09-25T09:57:35.374767Z","shell.execute_reply.started":"2023-09-25T09:57:35.360386Z","shell.execute_reply":"2023-09-25T09:57:35.373426Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import random\nimport ruptures as rpt \nimport gc","metadata":{"execution":{"iopub.status.busy":"2023-09-25T09:57:35.376125Z","iopub.execute_input":"2023-09-25T09:57:35.376553Z","iopub.status.idle":"2023-09-25T09:57:35.468197Z","shell.execute_reply.started":"2023-09-25T09:57:35.376528Z","shell.execute_reply":"2023-09-25T09:57:35.466982Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"algo = rpt.KernelCPD(kernel=\"rbf\", min_size=120, jump=2)\npenalty_value = 25 # beta","metadata":{"execution":{"iopub.status.busy":"2023-09-25T09:57:35.469684Z","iopub.execute_input":"2023-09-25T09:57:35.470007Z","iopub.status.idle":"2023-09-25T09:57:35.475323Z","shell.execute_reply.started":"2023-09-25T09:57:35.469981Z","shell.execute_reply":"2023-09-25T09:57:35.474249Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def compute_stats_in_window(column, window, ):\n    rw = column.rolling(window, center=True, step=1)\n    return (rw.mean(), rw.median(), rw.std())\n\ndef generate_enmo_df(frag, series_id):\n    sub_df = frag.to_table().to_pandas()\n    sub_df[\"enmo_delta\"] = compute_delta(sub_df)\n    sub_df[\"enmo_delta_abs\"] = sub_df[\"enmo_delta\"].apply(np.abs) \n    sub_df[\"enmo_delta_abs_scaled\"] = sub_df[\"enmo_delta_abs\"] * 1e6# Scale for clarity\n    sub_df['step'] = sub_df['step'].astype(int)\n    sub_df = sub_df.sort_values(\"step\").copy().reset_index(drop=True)\n\n    sub_df['anglez'] = sub_df['anglez'] * (np.pi / 180.0)\n    sub_df['anglez_mean'], sub_df['anglez_median'], sub_df['anglez_std'] = compute_stats_in_window(sub_df['anglez'], WINDOW_SIZE)\n    \n    sub_df[\"series_id\"] = series_id\n    return sub_df\n\n\ndef generate_lr_coeffs(sub_df):\n    chunk_size = int((sub_df[\"enmo_delta_abs_scaled\"].values.shape[0]+window-1)//window)\n    chunk_size = min(chunk_size, sub_df.shape[0])\n    arrs = np.array_split(sub_df[\"enmo_delta_abs_scaled\"].values, chunk_size)\n    coeffs = []\n    lns = []\n    for y in arrs:\n        X = np.arange(y.shape[0]).reshape(-1, 1)\n        reg = LinearRegression(fit_intercept=False).fit(X, y)\n        coeffs.append(reg.coef_[0])\n        lns.append(y.shape[0])\n\n    gc.collect()\n    #     print(\"LR Coeff- Piecemeal determined\")\n    # Post process: Map splits back to steps.\n    step = 0\n    steps = []\n    for x in lns:\n        steps.append(step)\n        step += x\n\n    # Post process: LR df\n    lrdf = pd.DataFrame()\n    lrdf[\"step\"] = steps\n    lrdf[\"coeff\"] = coeffs\n\n    return lrdf","metadata":{"execution":{"iopub.status.busy":"2023-09-25T09:57:35.478172Z","iopub.execute_input":"2023-09-25T09:57:35.478476Z","iopub.status.idle":"2023-09-25T09:57:35.492895Z","shell.execute_reply.started":"2023-09-25T09:57:35.478448Z","shell.execute_reply":"2023-09-25T09:57:35.491877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def refine_windows(df, result):\n    series_id = df['series_id'].iloc[0]\n    \n    N = df.shape[0]\n\n    # find if a windows is valid\n    windows = list(itertools.pairwise([0,] + result))\n    \n    # Classify the valid_windows as night / day based on length or enmo\n    rows = []\n    for s, e in windows:\n        sub_df = df[df['step'].between(s, e)]\n        if sub_df['anglez'].std() < 0.35:\n            continue\n        \n        window_size = e - s + 1\n        window_enmo_median = sub_df['enmo_delta_abs'].median()\n        \n        print(f'[{s}, {e}] - Size: {window_size}, median: {window_enmo_median}:.3f')\n        \n        if window_enmo_median < 1e-3:\n            rows.append((series_id, s, 'onset'))\n            rows.append((series_id, e, 'wakeup'))\n\n    preds = pd.DataFrame(rows, columns=['series_id', 'step', 'event'])\n    preds.reset_index(drop=False).rename(columns={'index': 'row_id'})\n    \n    return preds","metadata":{"execution":{"iopub.status.busy":"2023-09-25T09:57:35.494147Z","iopub.execute_input":"2023-09-25T09:57:35.494458Z","iopub.status.idle":"2023-09-25T09:57:35.509589Z","shell.execute_reply.started":"2023-09-25T09:57:35.494433Z","shell.execute_reply":"2023-09-25T09:57:35.508418Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def generate_full_breakpoints(frag):\n    series_id = get_series_id_from_fragment(frag)\n\n    # 1. Prep time series\n    sub_df = generate_enmo_df(frag, series_id)\n    #     print(\"GT Prepped\")\n\n    # 3. Window-wise Linear Reg without Bias. --> 30 min\n    lrdf = generate_lr_coeffs(sub_df)\n    lrdf = lrdf.reset_index(drop=True)\n\n    # 4. Predict breakpoints\n    signal = lrdf[\"coeff\"].values\n    try:\n        result = algo.fit_predict(\n            signal,\n            pen=penalty_value\n        )\n    except Exception:\n        return pd.DataFrame([[series_id, 0, \"onset\", 0.0], [series_id, 1, \"wakeup\", 0.0]], columns=[\"series_id\", \"step\", \"event\", \"score\"])\n    # 5. Map back to enmo dataset\n    lrdf[\"predicted_change_pt\"] = 0\n    lrdf.loc[result[:-1], \"predicted_change_pt\"] = 1\n    lrdf[\"step_adjusted\"] = lrdf[\"step\"] #+180 #TODO:\n    result_raw = lrdf.loc[lrdf[\"predicted_change_pt\"]==1, \"step_adjusted\"].values\n    \n    test = refine_windows(sub_df, result_raw)\n    test[\"score\"] = 1.0\n    \n    return test","metadata":{"execution":{"iopub.status.busy":"2023-09-25T09:57:35.510701Z","iopub.execute_input":"2023-09-25T09:57:35.511091Z","iopub.status.idle":"2023-09-25T09:57:35.528480Z","shell.execute_reply.started":"2023-09-25T09:57:35.511065Z","shell.execute_reply":"2023-09-25T09:57:35.527493Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission = []\nfor frag in pq_ds.fragments:\n    test = generate_full_breakpoints(frag)\n    submission.append(test)\n\nsubmission = pd.concat(submission, ignore_index=True).reset_index(drop=False).rename(columns={\"index\": \"row_id\"})\nsubmission.to_csv(\"submission.csv\", index=False)","metadata":{"execution":{"iopub.status.busy":"2023-09-25T09:57:35.529625Z","iopub.execute_input":"2023-09-25T09:57:35.529925Z","iopub.status.idle":"2023-09-25T09:57:35.859309Z","shell.execute_reply.started":"2023-09-25T09:57:35.529877Z","shell.execute_reply":"2023-09-25T09:57:35.858223Z"},"trusted":true},"execution_count":null,"outputs":[]}]}