{"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"},{"sourceId":6890980,"sourceType":"datasetVersion","datasetId":3958630}],"dockerImageVersionId":30558,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"EVALUATE = True\nDEBUG = False\nUSE_DEV_TEST = True\nNUM_SAMPLES = 10","metadata":{"execution":{"iopub.status.busy":"2023-12-06T01:49:08.899368Z","iopub.execute_input":"2023-12-06T01:49:08.899944Z","iopub.status.idle":"2023-12-06T01:49:08.917201Z","shell.execute_reply.started":"2023-12-06T01:49:08.899876Z","shell.execute_reply":"2023-12-06T01:49:08.915784Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- v1: original version - AP(local)=0.306\n- v2: bias corrected - AP(local)=0.314","metadata":{}},{"cell_type":"code","source":"from datetime import datetime as dt\nfrom datetime import timedelta\nfrom pathlib import Path\n\nimport numpy as np\nimport polars as pl\nfrom matplotlib.dates import date2num\nfrom matplotlib import pyplot as plt\nfrom tqdm.auto import tqdm\n\nplt.style.use(\"ggplot\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-12-06T01:49:10.160210Z","iopub.execute_input":"2023-12-06T01:49:10.161779Z","iopub.status.idle":"2023-12-06T01:49:10.382003Z","shell.execute_reply.started":"2023-12-06T01:49:10.161721Z","shell.execute_reply":"2023-12-06T01:49:10.380676Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATA_PATH = Path(\"/kaggle/input/child-mind-institute-detect-sleep-states\")\nTRAIN_EVENTS_PATH = DATA_PATH / \"train_events.csv\"\nTRAIN_SERIES_PATH = DATA_PATH / \"train_series.parquet\"\nTEST_SERIES_PATH = DATA_PATH / \"test_series.parquet\"\nSAMPLE_SUBMISSION_PATH = DATA_PATH / \"sample_submission.csv\"\nSLEEP_DATA_PATH = Path(\"/kaggle/input/sleep-data/sleep_time_probability.parquet\")\n\nfor path in [\n    DATA_PATH,\n    TRAIN_EVENTS_PATH,\n    TRAIN_SERIES_PATH,\n    TEST_SERIES_PATH,\n    SAMPLE_SUBMISSION_PATH,\n    SLEEP_DATA_PATH,\n]:\n    assert path.exists()","metadata":{"execution":{"iopub.status.busy":"2023-12-06T01:49:10.855063Z","iopub.execute_input":"2023-12-06T01:49:10.855667Z","iopub.status.idle":"2023-12-06T01:49:10.869104Z","shell.execute_reply.started":"2023-12-06T01:49:10.855612Z","shell.execute_reply":"2023-12-06T01:49:10.867642Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import polars as pl\n\n\ndef transform(df, night_offset=20):\n    return (\n        df.with_columns(\n            [\n                pl.col(\"timestamp\").str.slice(-5, 3).cast(pl.Int8).alias(\"tz_offset\"),\n            ]\n        )\n        .with_columns(\n            [\n                (pl.col(\"tz_offset\") == -4).alias(\"is_dst\"),\n            ]\n        )\n        .with_columns(\n            [\n                pl.col(\"timestamp\")\n                .str.strptime(pl.Datetime, format=\"%Y-%m-%dT%H:%M:%S%z\")\n                .alias(\"timestamp\"),\n            ]\n        )\n        .with_columns(\n            [\n                (pl.col(\"timestamp\").dt.year() - 2000).cast(pl.Int8).alias(\"year\"),\n                pl.col(\"timestamp\").dt.month().cast(pl.Int8).alias(\"month\"),\n                pl.col(\"timestamp\").dt.day().cast(pl.Int8).alias(\"day\"),\n                pl.col(\"timestamp\").dt.hour().cast(pl.Int8).alias(\"hour\"),\n                pl.col(\"timestamp\").dt.minute().cast(pl.Int8).alias(\"minute\"),\n                pl.col(\"timestamp\").dt.second().cast(pl.Int8).alias(\"second\"),\n                pl.col(\"timestamp\").dt.weekday().cast(pl.Int8).alias(\"weekday\"),\n            ]\n        )\n        .with_columns(  # 正午をまたいで日付を調整\n            pl.when(pl.col(\"hour\") < night_offset)\n            .then(pl.col(\"timestamp\"))\n            .otherwise(pl.col(\"timestamp\") + pl.duration(days=1))\n            .dt.date()\n            .alias(\"night_group\"),\n        )\n        .with_columns(\n            [\n                (\n                    pl.col(\"series_id\")\n                    + pl.lit(\"_\")\n                    + pl.col(\"night_group\").cast(pl.Datetime).dt.strftime(\"%Y%m%d\")\n                ).alias(\"group_id\"),\n            ]\n        )\n        .with_columns(\n            [\n                pl.col(\"timestamp\").cumcount().over(\"group_id\").alias(\"norm_step\"),\n            ]\n        )\n        .drop([\"night_group\"])\n    )\n\n\ndef transform_series(df):\n    return transform(df).with_columns(\n        [\n            (pl.col(\"enmo\") == 0).alias(\"is_enmo_clipped\"),\n        ]\n    )\n\n\ndef transform_events(df):\n    return (\n        transform(df)\n        .with_columns(\n            [\n                pl.col(\"night\").cast(pl.UInt32).alias(\"night\"),\n            ]\n        )\n        .pivot([\"step\", \"timestamp\", \"tz_offset\"], [\"series_id\", \"group_id\", \"night\"], \"event\")\n    )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-06T01:49:11.187875Z","iopub.execute_input":"2023-12-06T01:49:11.188756Z","iopub.status.idle":"2023-12-06T01:49:11.215040Z","shell.execute_reply.started":"2023-12-06T01:49:11.188706Z","shell.execute_reply":"2023-12-06T01:49:11.213614Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def add_feature(\n    df,\n    group_col=\"series_id\",\n    day_group_col=\"group_id\",\n    term1=(5 * 60) // 5,\n    term2=(30 * 60) // 5,\n    term3=(60 * 60) // 5,\n    min_threshold=0.005,\n    max_threshold=0.04,\n    center=True,\n):\n    return (\n        df.with_columns(\n            [\n                pl.col(\"anglez\").diff(1).abs().over(group_col).alias(\"anglez_diff\"),\n                pl.col(\"enmo\").diff(1).abs().over(group_col).alias(\"enmo_diff\"),\n            ]\n        )\n        .with_columns(\n            [\n                pl.col(\"anglez_diff\")\n                .rolling_median(term1, center=center)  # 5 min window\n                .over(group_col)\n                .alias(\"anglez_diff_median_5min\"),\n                pl.col(\"enmo_diff\")\n                .rolling_median(term1, center=center)  # 5 min window\n                .over(group_col)\n                .alias(\"enmo_diff_median_5min\"),\n            ]\n        )\n        .with_columns(\n            [\n                pl.col(\"anglez_diff_median_5min\")\n                .quantile(0.1)\n                .clip(min_threshold, max_threshold)\n                .over(day_group_col)\n                .alias(\"critical_threshold\")\n            ]\n        )\n        .with_columns(\n            [\n                (pl.col(\"anglez_diff_median_5min\") < pl.col(\"critical_threshold\") * 15)\n                .over(group_col)\n                .alias(\"is_static\")\n            ]\n        )\n        .with_columns(\n            [\n                pl.col(\"is_static\")\n                .cast(pl.Int32)\n                .rolling_sum(term2, center=center)\n                .over(group_col)\n                .alias(\"is_static_sum_30min\"),\n            ]\n        )\n        .with_columns(\n            [(pl.col(\"is_static_sum_30min\") == ((30 * 60) // 5)).over(group_col).alias(\"tmp\")]\n        )\n        .with_columns(\n            [\n                pl.col(\"tmp\").shift(term2 // 2).over(group_col).alias(\"tmp_left\"),\n                pl.col(\"tmp\").shift(-(term2 // 2)).over(group_col).alias(\"tmp_right\"),\n            ]\n        )\n        .with_columns(\n            [\n                (pl.col(\"tmp_left\") | pl.col(\"tmp_right\")).alias(\"is_sleep_block\"),\n            ]\n        )\n        .drop([\"tmp\", \"tmp_left\", \"tmp_right\"])\n        .with_columns([pl.col(\"is_sleep_block\").not_().alias(\"is_gap\")])\n        .with_columns(\n            [\n                pl.col(\"is_gap\")\n                .cast(pl.Int32)\n                .rolling_sum(term3, center=center)\n                .over(group_col)\n                .alias(\"gap_length\")\n            ]\n        )\n        .with_columns([(pl.col(\"gap_length\") == term3).over(group_col).alias(\"tmp\")])\n        .with_columns(\n            [\n                pl.col(\"tmp\").shift(term3 // 2).over(group_col).alias(\"tmp_left\"),\n                pl.col(\"tmp\").shift(-(term3 // 2)).over(group_col).alias(\"tmp_right\"),\n            ]\n        )\n        .with_columns(\n            [\n                (pl.col(\"tmp_left\") | pl.col(\"tmp_right\")).alias(\"is_large_gap\"),\n            ]\n        )\n        .drop([\"tmp\", \"tmp_left\", \"tmp_right\"])\n        .with_columns([pl.col(\"is_large_gap\").not_().alias(\"is_sleep_episode\")])\n        #\n        # extract longest sleep episode\n        #\n        .with_columns(\n            [\n                # extract false->true transition\n                (\n                    (\n                        pl.col(\"is_sleep_episode\")\n                        & pl.col(\"is_sleep_episode\")\n                        .shift_and_fill(pl.lit(False), periods=1)\n                        .not_()\n                    )\n                    .cumsum()\n                    .over(\"group_id\")\n                ).alias(\"sleep_episode_id\")\n            ]\n        )\n        .with_columns(\n            [\n                pl.col(\"is_sleep_episode\")\n                .sum()\n                .over([\"group_id\", \"sleep_episode_id\"])\n                .alias(\"sleep_episode_length\")\n            ]\n        )\n        .with_columns(\n            [\n                pl.col(\"sleep_episode_length\")\n                .max()\n                .over([\"group_id\"])\n                .alias(\"max_sleep_episode_length\")\n            ]\n        )\n        .with_columns(\n            [\n                (\n                    pl.col(\"is_sleep_episode\")\n                    & (pl.col(\"sleep_episode_length\") == pl.col(\"max_sleep_episode_length\"))\n                ).alias(\"is_longest_sleep_episode\")\n            ]\n        )\n    )","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-06T01:49:11.941676Z","iopub.execute_input":"2023-12-06T01:49:11.942099Z","iopub.status.idle":"2023-12-06T01:49:11.965539Z","shell.execute_reply.started":"2023-12-06T01:49:11.942067Z","shell.execute_reply":"2023-12-06T01:49:11.964088Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"Event Detection Average Precision\n\nAn average precision metric for event detection in time series and\nvideo.\n\nAttrubution: tubo213\nhttps://github.com/tubo213/kaggle-child-mind-institute-detect-sleep-states/tree/main\n\"\"\"\n\nfrom bisect import bisect_left\nfrom typing import Dict, List, Tuple\n\nimport numpy as np\nimport pandas as pd\nfrom tqdm.auto import tqdm\n\n\nclass ParticipantVisibleError(Exception):\n    pass\n\n\n# Set some placeholders for global parameters\nseries_id_column_name = \"series_id\"\ntime_column_name = \"step\"\nevent_column_name = \"event\"\nscore_column_name = \"score\"\nuse_scoring_intervals = False\ntolerances = {\n    \"onset\": [12, 36, 60, 90, 120, 150, 180, 240, 300, 360],\n    \"wakeup\": [12, 36, 60, 90, 120, 150, 180, 240, 300, 360],\n}\n\n\ndef score(\n    solution: pd.DataFrame,\n    submission: pd.DataFrame,\n    tolerances: Dict[str, List[float]],\n    series_id_column_name: str,\n    time_column_name: str,\n    event_column_name: str,\n    score_column_name: str,\n    use_scoring_intervals: bool = False,\n    verbose: bool = True,\n) -> Tuple[float, pd.DataFrame, pd.DataFrame]:\n    # Validate metric parameters\n    assert len(tolerances) > 0, \"Events must have defined tolerances.\"\n    assert set(tolerances.keys()) == set(solution[event_column_name]).difference(\n        {\"start\", \"end\"}\n    ), (\n        f\"Solution column {event_column_name} must contain the same events \"\n        \"as defined in tolerances.\"\n    )\n    assert pd.api.types.is_numeric_dtype(\n        solution[time_column_name]\n    ), f\"Solution column {time_column_name} must be of numeric type.\"\n\n    # Validate submission format\n    for column_name in [\n        series_id_column_name,\n        time_column_name,\n        event_column_name,\n        score_column_name,\n    ]:\n        if column_name not in submission.columns:\n            raise ParticipantVisibleError(\n                f\"Submission must have column '{column_name}'.\"\n            )\n\n    if not pd.api.types.is_numeric_dtype(submission[time_column_name]):\n        raise ParticipantVisibleError(\n            f\"Submission column '{time_column_name}' must be of numeric type.\"\n        )\n    if not pd.api.types.is_numeric_dtype(submission[score_column_name]):\n        raise ParticipantVisibleError(\n            f\"Submission column '{score_column_name}' must be of numeric type.\"\n        )\n\n    # Set these globally to avoid passing around a bunch of arguments\n    globals()[\"series_id_column_name\"] = series_id_column_name\n    globals()[\"time_column_name\"] = time_column_name\n    globals()[\"event_column_name\"] = event_column_name\n    globals()[\"score_column_name\"] = score_column_name\n    globals()[\"use_scoring_intervals\"] = use_scoring_intervals\n\n    return event_detection_ap(solution, submission, tolerances, verbose=verbose)\n\n\ndef find_nearest(xs: np.ndarray, value):\n    \"\"\"\n    Find the index of the closest value to x in the array xs.\n    \"\"\"\n    idx = np.searchsorted(xs, value, side=\"left\")\n    best_idx = None\n    best_error = float(\"inf\")\n    best_diff = float(\"inf\")\n\n    range_min = max(0, idx - 1)\n    range_max = min(len(xs), idx + 2)\n    for check_idx in range(\n        range_min, range_max\n    ):  # Check the exact, one before, and one after\n        error = abs(xs[check_idx] - value)\n        if error < best_error:\n            best_error = error\n            best_idx = check_idx\n            best_diff = xs[check_idx] - value\n\n    return best_idx, best_error, best_diff\n\n\ndef find_nearest_time_idx(sorted_gt_times, det_time, excluded_indices: set):\n    \"\"\"\n    search index of gt_times closest to det_time.\n\n    assumes gt_times is sorted in ascending order.\n    \"\"\"\n    # e.g. if gt_times = [0, 1, 2, 3, 4, 5] and det_time = 2.5, then idx = 3\n    sorted_gt_times = np.asarray(sorted_gt_times)\n    available_indices = np.asarray(\n        sorted(set(range(len(sorted_gt_times))) - excluded_indices), dtype=int\n    )\n    sorted_gt_times = sorted_gt_times[available_indices]\n    idx, error, diff = find_nearest(sorted_gt_times, det_time)\n    best_idx = available_indices[idx] if idx is not None else None\n\n    return best_idx, error, diff\n\n\ndef match_detections(\n    tolerance: float, ground_truths: pd.DataFrame, detections: pd.DataFrame\n) -> pd.DataFrame:\n    detections_sorted = detections.sort_values(\n        score_column_name, ascending=False\n    ).dropna()\n    is_matched = np.full_like(detections_sorted[event_column_name], False, dtype=bool)\n    diffs = np.full_like(\n        detections_sorted[event_column_name], float(\"inf\"), dtype=float\n    )\n    ground_truths_times = ground_truths.sort_values(time_column_name)[\n        time_column_name\n    ].to_list()\n    matched_gt_indices: set[int] = set()\n\n    for i, det in enumerate(detections_sorted.itertuples(index=False)):\n        det_time = getattr(det, time_column_name)\n\n        best_idx, best_error, best_diff = find_nearest_time_idx(\n            ground_truths_times, det_time, matched_gt_indices\n        )\n\n        if (best_idx is not None) and (best_error < tolerance):\n            is_matched[i] = True\n            diffs[i] = best_diff\n            matched_gt_indices.add(best_idx)\n\n    detections_sorted[\"matched\"] = is_matched\n    detections_sorted[\"diff\"] = diffs\n    return detections_sorted\n\n\ndef precision_recall_curve(\n    matches: np.ndarray, scores: np.ndarray, p: int\n) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:\n    if len(matches) == 0:\n        return [1], [0], []\n\n    # Sort matches by decreasing confidence\n    idxs = np.argsort(scores, kind=\"stable\")[::-1]\n    scores = scores[idxs]\n    matches = matches[idxs]\n\n    distinct_value_indices = np.where(np.diff(scores))[0]\n    threshold_idxs = np.r_[distinct_value_indices, matches.size - 1]\n    thresholds = scores[threshold_idxs]\n\n    # Matches become TPs and non-matches FPs as confidence threshold decreases\n    tps = np.cumsum(matches)[threshold_idxs]\n    fps = np.cumsum(~matches)[threshold_idxs]\n\n    precision = tps / (tps + fps)\n    precision[np.isnan(precision)] = 0\n    recall = (\n        tps / p\n    )  # total number of ground truths might be different than total number of matches\n\n    # Stop when full recall attained and reverse the outputs so recall is non-increasing.\n    last_ind = tps.searchsorted(tps[-1])\n    sl = slice(last_ind, None, -1)\n\n    # Final precision is 1 and final recall is 0\n    return np.r_[precision[sl], 1], np.r_[recall[sl], 0], thresholds[sl]\n\n\ndef average_precision_score(matches: np.ndarray, scores: np.ndarray, p: int) -> float:\n    precision, recall, _ = precision_recall_curve(matches, scores, p)\n    # Compute step integral\n    return -np.sum(np.diff(recall) * np.array(precision)[:-1])\n\n\ndef event_detection_ap(\n    solution: pd.DataFrame,\n    submission: pd.DataFrame,\n    tolerances: Dict[str, List[float]] = tolerances,\n    progress_bar: bool = True,\n    verbose: bool = True,\n) -> Tuple[float, pd.DataFrame, pd.DataFrame]:\n    # Ensure solution and submission are sorted properly\n    solution = solution.sort_values([series_id_column_name, time_column_name])\n    submission = submission.sort_values([series_id_column_name, time_column_name])\n\n    # Extract scoring intervals.\n    if use_scoring_intervals:\n        raise NotImplementedError(\"Scoring intervals not implemented.\")\n\n    # Extract ground-truth events.\n    ground_truths = solution.query(\"event not in ['start', 'end']\").reset_index(\n        drop=True\n    )\n\n    # Map each event class to its prevalence (needed for recall calculation)\n    class_counts = ground_truths.value_counts(event_column_name).to_dict()\n\n    # Create table for detections with a column indicating a match to a ground-truth event\n    detections = submission.assign(matched=False)\n\n    # Remove detections outside of scoring intervals\n    if use_scoring_intervals:\n        raise NotImplementedError(\"Scoring intervals not implemented.\")\n    else:\n        detections_filtered = detections\n\n    # Create table of event-class x tolerance x series_id values\n    aggregation_keys = pd.DataFrame(\n        [\n            (ev, tol, vid)\n            for ev in tolerances.keys()\n            for tol in tolerances[ev]\n            for vid in ground_truths[series_id_column_name].unique()\n        ],\n        columns=[event_column_name, \"tolerance\", series_id_column_name],\n    )\n\n    # Create match evaluation groups: event-class x tolerance x series_id\n    detections_grouped = aggregation_keys.merge(\n        detections_filtered, on=[event_column_name, series_id_column_name], how=\"left\"\n    ).groupby([event_column_name, \"tolerance\", series_id_column_name])\n    ground_truths_grouped = aggregation_keys.merge(\n        ground_truths, on=[event_column_name, series_id_column_name], how=\"left\"\n    ).groupby([event_column_name, \"tolerance\", series_id_column_name])\n\n    # Match detections to ground truth events by evaluation group\n    pbars = aggregation_keys.itertuples(index=False)\n    if progress_bar:\n        pbars = tqdm(pbars, total=len(aggregation_keys), desc=\"Matching detections\")\n    detections_matched = []\n    for key in pbars:\n        dets = detections_grouped.get_group(key)\n        gts = ground_truths_grouped.get_group(key)\n        detections_matched.append(\n            match_detections(dets[\"tolerance\"].iloc[0], gts, dets)\n        )\n    detections_matched = pd.concat(detections_matched)\n\n    # Compute AP per event x tolerance group\n    event_classes = ground_truths[event_column_name].unique()\n    ap_table = (\n        detections_matched.query(\"event in @event_classes\")\n        .groupby([event_column_name, \"tolerance\"])\n        .apply(\n            lambda group: average_precision_score(\n                group[\"matched\"].to_numpy(),\n                group[score_column_name].to_numpy(),\n                class_counts[group[event_column_name].iat[0]],\n            )\n        )\n        .reset_index()\n        .pivot(index=\"tolerance\", columns=\"event\", values=0)\n    )\n    if verbose:\n        display(ap_table)\n    # Average over tolerances, then over event classes\n    mean_ap = ap_table.mean().mean()\n\n    return mean_ap, ap_table, detections_matched","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-06T01:49:12.417902Z","iopub.execute_input":"2023-12-06T01:49:12.418299Z","iopub.status.idle":"2023-12-06T01:49:12.952197Z","shell.execute_reply.started":"2023-12-06T01:49:12.418268Z","shell.execute_reply":"2023-12-06T01:49:12.950785Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nif EVALUATE:\n    tr_events_raw = pl.read_csv(TRAIN_EVENTS_PATH)\n    tr_series = pl.read_parquet(TRAIN_SERIES_PATH)\n    sample_submission = pl.read_csv(SAMPLE_SUBMISSION_PATH)\n\n    series_ids = tr_events_raw[\"series_id\"].unique()\n    print(f\"#Events: {len(tr_events_raw)}\")\n    print(f\"#Series: {len(series_ids)}\")\n\n    if DEBUG:\n        sample_ids = series_ids.sample(NUM_SAMPLES, seed=42).to_list()\n        tr_series = tr_series.filter(pl.col(\"series_id\").is_in(sample_ids))\n        tr_events_raw = tr_events_raw.filter(pl.col(\"series_id\").is_in(sample_ids))\n\n    tr_series = transform_series(tr_series)\n    tr_series = add_feature(tr_series)\n    tr_events = transform_events(tr_events_raw)\nelse:\n    tr_series = pl.read_parquet(TEST_SERIES_PATH)\n    sample_submission = pl.read_csv(SAMPLE_SUBMISSION_PATH)\n    tr_series = add_feature(transform_series(tr_series))","metadata":{"execution":{"iopub.status.busy":"2023-12-06T01:49:12.954142Z","iopub.execute_input":"2023-12-06T01:49:12.954651Z","iopub.status.idle":"2023-12-06T01:49:54.036772Z","shell.execute_reply.started":"2023-12-06T01:49:12.954617Z","shell.execute_reply":"2023-12-06T01:49:54.035522Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_groups(df_events, df_series, group_ids, eps=1e-5):\n    target_cols = [\"anglez\", \"enmo\"]\n\n    for i, group_id in enumerate(group_ids):\n        print(f\"{i+1}/{len(group_ids)}\")\n        se = df_series.filter(pl.col(\"group_id\") == group_id)\n        ev = df_events.filter(pl.col(\"group_id\") == group_id)\n        _, axes = plt.subplots(2, 1, figsize=(24, 6))\n        display(ev)\n        print(f\"group_id={group_id}\")\n        print(se[\"critical_threshold\"][0])\n\n        for ax, target_col in zip(axes, target_cols):\n            y_min = se[target_col].min()\n            y_max = se[target_col].max()\n\n            ax.set(\n                xlabel=\"timestamp\",\n                ylabel=target_col,\n            )\n            if target_col == \"enmo\":\n                ax.set(yscale=\"log\")\n                y_min, y_max = 0.8e-5, 1.2e0\n                ax.set(ylim=(y_min, y_max))\n                offset = eps\n                ax.scatter(\n                    se[\"timestamp\"],\n                    se[\"is_static\"] * 1.0,\n                    color=\"green\",\n                    s=1,\n                    label=\"is_static\",\n                )\n                ax.scatter(\n                    se[\"timestamp\"],\n                    se[\"is_sleep_block\"] * 0.5,\n                    color=\"blue\",\n                    s=1,\n                    label=\"is_sleep_block\",\n                )\n                ax.scatter(\n                    se[\"timestamp\"],\n                    se[\"is_sleep_episode\"] * 0.2,\n                    color=\"purple\",\n                    s=1,\n                    label=\"sleep episode\",\n                )\n                ax.scatter(\n                    se[\"timestamp\"],\n                    se[\"is_longest_sleep_episode\"] * 0.1,\n                    color=\"orange\",\n                    s=1,\n                    label=\"longest sleep episode\",\n                )\n                ax.scatter(\n                    se[\"timestamp\"],\n                    se[\"enmo_diff_median_5min\"] * 5 + offset,\n                    s=1,\n                    color=\"red\",\n                    label=\"enmo_diff_median_5min (x5)\",\n                )\n            else:\n                y_min, y_max = -92, 92\n                ax.set(ylim=(y_min, y_max))\n                offset = 0\n                ax.scatter(\n                    se[\"timestamp\"],\n                    se[\"anglez_diff_median_5min\"] * 5,\n                    s=1,\n                    color=\"red\",\n                    label=\"anglez_diff_median_5min (x5)\",\n                )\n            ax.scatter(\n                se[\"timestamp\"], se[target_col] + offset, s=1, color=\"black\", label=target_col\n            )\n            ax.legend(loc=\"lower right\")\n            plt.suptitle(f\"group_id={group_id}\")\n\n            for item in ev.iter_rows():\n                data = dict(zip(ev.columns, item))\n                x_min, x_max = (\n                    data[\"timestamp_event_onset\"],\n                    data[\"timestamp_event_wakeup\"],\n                )\n                ax.fill_between(\n                    [date2num(x_min), date2num(x_max)],\n                    y1=y_min,\n                    y2=y_max,\n                    alpha=0.2,\n                    color=\"blue\",\n                    label=None,\n                )\n        plt.tight_layout()\n        plt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-12-06T01:49:54.038637Z","iopub.execute_input":"2023-12-06T01:49:54.039068Z","iopub.status.idle":"2023-12-06T01:49:54.060141Z","shell.execute_reply.started":"2023-12-06T01:49:54.039033Z","shell.execute_reply":"2023-12-06T01:49:54.058632Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if EVALUATE:\n    sample_gids = tr_series[\"group_id\"].unique().sample(NUM_SAMPLES, seed=123).to_list()\n    plot_groups(tr_events, tr_series, group_ids=sample_gids)","metadata":{"execution":{"iopub.status.busy":"2023-12-06T01:49:54.061659Z","iopub.execute_input":"2023-12-06T01:49:54.062072Z","iopub.status.idle":"2023-12-06T01:50:18.257077Z","shell.execute_reply.started":"2023-12-06T01:49:54.062037Z","shell.execute_reply":"2023-12-06T01:50:18.256055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sleep_data = pl.read_parquet(SLEEP_DATA_PATH)\npred_df = (\n    tr_series.filter(pl.col(\"is_longest_sleep_episode\"))\n    .group_by(\"group_id\")\n    .agg(\n        [\n            pl.col(\"series_id\").first().alias(\"series_id\"),\n            pl.col(\"step\").first().cast(pl.Int64).alias(\"step_event_onset\"),\n            pl.col(\"step\").last().cast(pl.Int64).alias(\"step_event_wakeup\"),\n            pl.col(\"timestamp\").first().alias(\"timestamp_event_onset\"),\n            pl.col(\"timestamp\").last().alias(\"timestamp_event_wakeup\"),\n        ]\n    )\n    .with_columns(\n        [(pl.col(\"timestamp_event_wakeup\") - pl.col(\"timestamp_event_onset\")).alias(\"duration\")]\n    )\n    .with_columns([(pl.col(\"duration\").dt.seconds() / 3600).alias(\"duration_hour\")])\n    .with_columns(\n        [\n            pl.col(\"duration_hour\").floor().alias(\"duration_start\"),\n            pl.col(\"duration_hour\").ceil().alias(\"duration_end\"),\n        ]\n    )\n    .join(sleep_data, on=[\"duration_start\", \"duration_end\"], how=\"left\")\n    .with_columns(\n        [\n            pl.col(\"probability\").fill_null(0.0),\n        ]\n    )\n)\ndisplay(pred_df.head())\nif EVALUATE:\n    (pred_df[\"duration\"].cast(pl.UInt64) / 3600 / 1e6).to_pandas().hist(bins=30)","metadata":{"execution":{"iopub.status.busy":"2023-12-06T01:50:18.259572Z","iopub.execute_input":"2023-12-06T01:50:18.260199Z","iopub.status.idle":"2023-12-06T01:50:18.920429Z","shell.execute_reply.started":"2023-12-06T01:50:18.260161Z","shell.execute_reply":"2023-12-06T01:50:18.919204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pred_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-12-06T01:50:18.922039Z","iopub.execute_input":"2023-12-06T01:50:18.922510Z","iopub.status.idle":"2023-12-06T01:50:18.933188Z","shell.execute_reply.started":"2023-12-06T01:50:18.922477Z","shell.execute_reply":"2023-12-06T01:50:18.931893Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission_raw = (\n    pl.concat(\n        [\n            pred_df.select([\"series_id\", \"step_event_onset\", \"probability\"])\n            .rename(dict(step_event_onset=\"step\"))\n            .with_columns([pl.lit(\"onset\").alias(\"event\")]),\n            pred_df.select([\"series_id\", \"step_event_wakeup\", \"probability\"])\n            .rename(dict(step_event_wakeup=\"step\"))\n            .with_columns([pl.lit(\"wakeup\").alias(\"event\")]),\n        ]\n    )\n    .rename(dict(probability=\"score\"))\n    .with_columns(pl.col(\"series_id\").cumcount().cast(pl.Int64).alias(\"row_id\"))\n    .sort([\"series_id\", \"step\", \"event\"])\n)\nsubmission_raw.head()","metadata":{"execution":{"iopub.status.busy":"2023-12-06T01:51:18.979311Z","iopub.execute_input":"2023-12-06T01:51:18.979826Z","iopub.status.idle":"2023-12-06T01:51:18.996111Z","shell.execute_reply.started":"2023-12-06T01:51:18.979786Z","shell.execute_reply":"2023-12-06T01:51:18.995015Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Bias Correction","metadata":{}},{"cell_type":"code","source":"if EVALUATE:\n    gt_events = tr_events_raw.drop_nulls()\n    submission_corrected = submission_raw.clone()\n    max_itr = 2\n    bias = dict(onset=0.0, wakeup=0.0)\n    best_iter = 0\n    best_ap = 0\n\n    tolerance_threshold = [12, 36, 60, 90, 120, 150, 180, 240, 300, 360]\n    tolerances = dict(\n        onset=tolerance_threshold,\n        wakeup=tolerance_threshold,\n    )\n    column_names = dict(\n        series_id_column_name=\"series_id\",\n        time_column_name=\"step\",\n        event_column_name=\"event\",\n        score_column_name=\"score\",\n    )\n\n    for itr in range(1, max_itr + 1):\n        print(f\"Iteration: {itr}/{max_itr}\")\n        submission_corrected = submission_corrected.with_columns(\n            [\n                pl.when(pl.col(\"event\") == \"onset\")\n                .then(pl.col(\"step\") + round(bias[\"onset\"]))\n                .otherwise(pl.col(\"step\") + round(bias[\"wakeup\"]))\n                .alias(\"step\"),\n            ]\n        )\n        mean_ap, _, detections_matched = score(\n            gt_events.drop_nulls().to_pandas(),\n            submission_corrected.to_pandas(),\n            tolerances,\n            **column_names,\n            verbose=True,\n        )\n        print(f\"AP score: {mean_ap:.3f}\")\n\n        diffs_onset = detections_matched.query(\"event == 'onset'\")[\"diff\"]\n        diffs_wakeup = detections_matched.query(\"event =='wakeup'\")[\"diff\"]\n\n        diffs_onset = diffs_onset[~np.isinf(diffs_onset)]\n        diffs_wakeup = diffs_wakeup[~np.isinf(diffs_wakeup)]\n        bias = dict(onset=diffs_onset.mean(), wakeup=diffs_wakeup.mean())\n        print(bias)\n\n        if mean_ap > best_ap:\n            best_ap = mean_ap\n            best_iter = itr\n\n        # plot bias\n        _, ax = plt.subplots()\n        diffs_onset.hist(ax=ax, bins=50, label=\"onset\", alpha=0.5)\n        diffs_wakeup.hist(ax=ax, bins=50, label=\"wakeup\", alpha=0.5)\n\n        ax.legend()\n        ax.set(xlabel=\"diff\", ylabel=\"count\")\n        plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-12-06T01:52:37.959646Z","iopub.execute_input":"2023-12-06T01:52:37.960145Z","iopub.status.idle":"2023-12-06T01:52:41.545703Z","shell.execute_reply.started":"2023-12-06T01:52:37.960111Z","shell.execute_reply":"2023-12-06T01:52:41.544213Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bias = dict(onset= 8.74, wakeup=-22.79)\nsubmission = submission_raw.with_columns(\n    [\n        pl.when(pl.col(\"event\") == \"onset\")\n        .then(pl.col(\"step\") + round(bias[\"onset\"]))\n        .otherwise(pl.col(\"step\") + round(bias[\"wakeup\"]))\n        .alias(\"step\"),\n    ]\n)","metadata":{"execution":{"iopub.status.busy":"2023-12-06T01:52:47.986299Z","iopub.execute_input":"2023-12-06T01:52:47.986801Z","iopub.status.idle":"2023-12-06T01:52:47.994664Z","shell.execute_reply.started":"2023-12-06T01:52:47.986761Z","shell.execute_reply":"2023-12-06T01:52:47.993682Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission.write_csv(\"submission.csv\")","metadata":{"execution":{"iopub.status.busy":"2023-12-06T01:52:48.401232Z","iopub.execute_input":"2023-12-06T01:52:48.402142Z","iopub.status.idle":"2023-12-06T01:52:48.412579Z","shell.execute_reply.started":"2023-12-06T01:52:48.402099Z","shell.execute_reply":"2023-12-06T01:52:48.411591Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pl.read_csv(\"submission.csv\")","metadata":{"execution":{"iopub.status.busy":"2023-12-06T01:52:49.119717Z","iopub.execute_input":"2023-12-06T01:52:49.120168Z","iopub.status.idle":"2023-12-06T01:52:49.252395Z","shell.execute_reply.started":"2023-12-06T01:52:49.120133Z","shell.execute_reply":"2023-12-06T01:52:49.251020Z"},"trusted":true},"execution_count":null,"outputs":[]}]}