{"cells":[{"cell_type":"code","execution_count":null,"id":"6ecd7090","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","execution":{"iopub.execute_input":"2025-10-20T02:50:15.104822Z","iopub.status.busy":"2025-10-20T02:50:15.104368Z","iopub.status.idle":"2025-10-20T02:50:26.506274Z","shell.execute_reply":"2025-10-20T02:50:26.505186Z"},"papermill":{"duration":11.415258,"end_time":"2025-10-20T02:50:26.508022","exception":false,"start_time":"2025-10-20T02:50:15.092764","status":"completed"},"tags":[]},"outputs":[],"source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport seaborn as sns   \nfrom scipy.signal import butter, filtfilt\nfrom scipy.stats import skew, kurtosis\n# import seglearn as sglearn        # For windowing and sequence modeling\nimport tsfresh     \nimport os\nfrom sklearn.preprocessing import StandardScaler\n\n\n \n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\n# import os\n# for dirname, _, filenames in os.walk('/kaggle/input'):\n#     for filename in filenames:\n#         print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session\nimport polars as pl\nimport dask.dataframe as dd\nfrom pathlib import Path"},{"cell_type":"markdown","id":"36356536","metadata":{"papermill":{"duration":0.006313,"end_time":"2025-10-20T02:50:26.521309","exception":false,"start_time":"2025-10-20T02:50:26.514996","status":"completed"},"tags":[]},"source":"# Data Exploration"},{"cell_type":"code","execution_count":null,"id":"92b6eb89","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:26.535773Z","iopub.status.busy":"2025-10-20T02:50:26.535273Z","iopub.status.idle":"2025-10-20T02:50:26.539756Z","shell.execute_reply":"2025-10-20T02:50:26.53901Z"},"papermill":{"duration":0.013279,"end_time":"2025-10-20T02:50:26.541034","exception":false,"start_time":"2025-10-20T02:50:26.527755","status":"completed"},"tags":[]},"outputs":[],"source":"# File paths for three training datasets\ndefog = Path('/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction/train/defog')\nnotype = Path('/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction/train/notype')\ntdcsfog = Path('/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction/train/tdcsfog')"},{"cell_type":"code","execution_count":null,"id":"6146f126","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:26.55664Z","iopub.status.busy":"2025-10-20T02:50:26.556321Z","iopub.status.idle":"2025-10-20T02:50:38.872267Z","shell.execute_reply":"2025-10-20T02:50:38.870852Z"},"papermill":{"duration":12.325156,"end_time":"2025-10-20T02:50:38.873698","exception":false,"start_time":"2025-10-20T02:50:26.548542","status":"completed"},"tags":[]},"outputs":[],"source":"defog_files = [f for f in os.listdir(defog) if f.endswith('.csv')]\n\n# List to store individual DataFrames\ndefog_list = []\n\nfor path in defog.glob(\"*.csv\"):\n    patient_id = path.stem  # removes .csv\n\n    df = pl.read_csv(path)\n    df = df.with_columns([\n        pl.lit(patient_id).alias(\"patient_id\")\n    ])\n    \n    defog_list.append(df)\n\ndefog_df = pl.concat(defog_list)\n# for f in defog_files:\n#     file_path = os.path.join(defog, f)\n#     df = pl.read_csv(file_path)\n#     df = df.with_columns([\n#         pl.lit(f).alias('file')  # Add filename as identifier\n#     ])\n#     defog_list.append(df)\n\n# # Concatenate into one large DataFrame\n# defog_df = pl.concat(defog_list)"},{"cell_type":"code","execution_count":null,"id":"1f96dab2","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:38.892919Z","iopub.status.busy":"2025-10-20T02:50:38.89144Z","iopub.status.idle":"2025-10-20T02:50:38.915454Z","shell.execute_reply":"2025-10-20T02:50:38.914421Z"},"papermill":{"duration":0.034594,"end_time":"2025-10-20T02:50:38.917561","exception":false,"start_time":"2025-10-20T02:50:38.882967","status":"completed"},"tags":[]},"outputs":[],"source":"defog_df.head()"},{"cell_type":"code","execution_count":null,"id":"96d0d4a5","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:38.933264Z","iopub.status.busy":"2025-10-20T02:50:38.932855Z","iopub.status.idle":"2025-10-20T02:50:52.33757Z","shell.execute_reply":"2025-10-20T02:50:52.336481Z"},"papermill":{"duration":13.414076,"end_time":"2025-10-20T02:50:52.339062","exception":false,"start_time":"2025-10-20T02:50:38.924986","status":"completed"},"tags":[]},"outputs":[],"source":"tdcsfog_files = [f for f in os.listdir(tdcsfog) if f.endswith('.csv')]\n\n# List to store individual DataFrames\ntdcsfog_list = []\n\nfor path in tdcsfog.glob(\"*.csv\"):\n    patient_id = path.stem  # removes .csv\n\n    df = pl.read_csv(path)\n    df = df.with_columns([\n        pl.lit(patient_id).alias(\"patient_id\")\n    ])\n    \n    tdcsfog_list.append(df)\n\ntdcsfog_df = pl.concat(tdcsfog_list)"},{"cell_type":"code","execution_count":null,"id":"160c7e7d","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:52.354607Z","iopub.status.busy":"2025-10-20T02:50:52.354225Z","iopub.status.idle":"2025-10-20T02:50:52.360426Z","shell.execute_reply":"2025-10-20T02:50:52.359643Z"},"papermill":{"duration":0.015544,"end_time":"2025-10-20T02:50:52.361961","exception":false,"start_time":"2025-10-20T02:50:52.346417","status":"completed"},"tags":[]},"outputs":[],"source":"tdcsfog_df.head()"},{"cell_type":"code","execution_count":null,"id":"96df1bd1","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:52.377858Z","iopub.status.busy":"2025-10-20T02:50:52.377578Z","iopub.status.idle":"2025-10-20T02:50:54.262408Z","shell.execute_reply":"2025-10-20T02:50:54.260977Z"},"papermill":{"duration":1.89533,"end_time":"2025-10-20T02:50:54.264458","exception":false,"start_time":"2025-10-20T02:50:52.369128","status":"completed"},"tags":[]},"outputs":[],"source":"print(defog_df.head())\n# print(defog_df.info())\nprint(defog_df.describe())\nprint(defog_df.shape)     # (rows, columns)\nprint(defog_df.columns)   # list of column names\nprint(defog_df.dtypes)    # list of column types"},{"cell_type":"code","execution_count":null,"id":"fa9114c9","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:54.280597Z","iopub.status.busy":"2025-10-20T02:50:54.279711Z","iopub.status.idle":"2025-10-20T02:50:55.261928Z","shell.execute_reply":"2025-10-20T02:50:55.260675Z"},"papermill":{"duration":0.992019,"end_time":"2025-10-20T02:50:55.263817","exception":false,"start_time":"2025-10-20T02:50:54.271798","status":"completed"},"tags":[]},"outputs":[],"source":"print(tdcsfog_df.head())\n# print(tdcsfog_df.info())\nprint(tdcsfog_df.shape)     # (rows, columns)\nprint(tdcsfog_df.columns)   # list of column names\nprint(tdcsfog_df.dtypes) \nprint(tdcsfog_df.describe())"},{"cell_type":"code","execution_count":null,"id":"bb0e7f1f","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:55.28073Z","iopub.status.busy":"2025-10-20T02:50:55.279751Z","iopub.status.idle":"2025-10-20T02:50:55.330132Z","shell.execute_reply":"2025-10-20T02:50:55.329164Z"},"papermill":{"duration":0.060136,"end_time":"2025-10-20T02:50:55.331587","exception":false,"start_time":"2025-10-20T02:50:55.271451","status":"completed"},"tags":[]},"outputs":[],"source":"events_df = pd.read_csv('/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction/events.csv')\nprint(events_df.head())\nprint(events_df.shape)   \nprint(events_df.columns)   \nprint(events_df.dtypes) \nprint(events_df.describe())"},{"cell_type":"code","execution_count":null,"id":"68805d6d","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:55.347476Z","iopub.status.busy":"2025-10-20T02:50:55.347076Z","iopub.status.idle":"2025-10-20T02:50:55.511647Z","shell.execute_reply":"2025-10-20T02:50:55.51058Z"},"papermill":{"duration":0.174297,"end_time":"2025-10-20T02:50:55.513386","exception":false,"start_time":"2025-10-20T02:50:55.339089","status":"completed"},"tags":[]},"outputs":[],"source":"unique_defog_patients = defog_df[\"patient_id\"].unique()\n\nprint(unique_defog_patients)"},{"cell_type":"code","execution_count":null,"id":"fe56bde2","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:55.529274Z","iopub.status.busy":"2025-10-20T02:50:55.529013Z","iopub.status.idle":"2025-10-20T02:50:56.340649Z","shell.execute_reply":"2025-10-20T02:50:56.339714Z"},"papermill":{"duration":0.821078,"end_time":"2025-10-20T02:50:56.342072","exception":false,"start_time":"2025-10-20T02:50:55.520994","status":"completed"},"tags":[]},"outputs":[],"source":"# 1. Filter your Polars DF for a single patient and convert to pandas\ndf = defog_df.filter(pl.col(\"patient_id\") == 'be9d33541d').to_pandas()\n\n# 2. Plot\nplt.figure(figsize=(15, 6))\n\n# Plot acceleration\nplt.plot(df['Time'], df['AccV'], label='AccV', alpha=0.7)\nplt.plot(df['Time'], df['AccML'], label='AccML', alpha=0.7)\nplt.plot(df['Time'], df['AccAP'], label='AccAP', alpha=0.7)\n\n# 3. Plot events\nplt.plot(df['Time'], df['StartHesitation'], label='StartHesitation', alpha=0.7)\nplt.plot(df['Time'], df['Turn'], label='Turn', alpha=0.7)\nplt.plot(df['Time'], df['Walking'], label='Walking', alpha=0.7)\n\n\n# 4. Final touches\nplt.xlabel(\"Time\")\nplt.ylabel(\"Acceleration (g)\")\nplt.title(f\"Patient: {patient_id} - Acceleration + FOG Events\")\nplt.legend(loc=\"upper right\")\nplt.grid(True)\nplt.tight_layout()\nplt.show()"},{"cell_type":"markdown","id":"70b5b6f1","metadata":{"papermill":{"duration":0.009042,"end_time":"2025-10-20T02:50:56.360729","exception":false,"start_time":"2025-10-20T02:50:56.351687","status":"completed"},"tags":[]},"source":"# Data Cleaning"},{"cell_type":"code","execution_count":null,"id":"87c66541","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:56.381329Z","iopub.status.busy":"2025-10-20T02:50:56.380423Z","iopub.status.idle":"2025-10-20T02:50:56.385903Z","shell.execute_reply":"2025-10-20T02:50:56.384904Z"},"papermill":{"duration":0.017519,"end_time":"2025-10-20T02:50:56.3876","exception":false,"start_time":"2025-10-20T02:50:56.370081","status":"completed"},"tags":[]},"outputs":[],"source":"# Data types of features \nprint(f'DEFOG DATA TYPES:\\n{defog_df.dtypes}\\n')\nprint(f'TDCSFOG DATA TYPES:\\n{tdcsfog_df.dtypes}\\n')"},{"cell_type":"code","execution_count":null,"id":"91c1f2a8","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:56.407527Z","iopub.status.busy":"2025-10-20T02:50:56.407237Z","iopub.status.idle":"2025-10-20T02:50:56.412732Z","shell.execute_reply":"2025-10-20T02:50:56.411566Z"},"papermill":{"duration":0.016992,"end_time":"2025-10-20T02:50:56.414126","exception":false,"start_time":"2025-10-20T02:50:56.397134","status":"completed"},"tags":[]},"outputs":[],"source":"print(tdcsfog_df.null_count())"},{"cell_type":"code","execution_count":null,"id":"4bd7dd3c","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:56.434794Z","iopub.status.busy":"2025-10-20T02:50:56.434544Z","iopub.status.idle":"2025-10-20T02:50:56.841512Z","shell.execute_reply":"2025-10-20T02:50:56.840243Z"},"papermill":{"duration":0.418667,"end_time":"2025-10-20T02:50:56.842952","exception":false,"start_time":"2025-10-20T02:50:56.424285","status":"completed"},"tags":[]},"outputs":[],"source":"# Convert accerlations in defog to m/s^2\nG_CONVERSION = 9.80665\ndefog_df[[\"AccV\", \"AccML\", \"AccAP\"]] *= G_CONVERSION\nprint(defog_df)"},{"cell_type":"code","execution_count":null,"id":"93d3925d","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:56.864747Z","iopub.status.busy":"2025-10-20T02:50:56.864205Z","iopub.status.idle":"2025-10-20T02:50:57.244794Z","shell.execute_reply":"2025-10-20T02:50:57.243897Z"},"papermill":{"duration":0.393553,"end_time":"2025-10-20T02:50:57.246277","exception":false,"start_time":"2025-10-20T02:50:56.852724","status":"completed"},"tags":[]},"outputs":[],"source":"# Convert the Valid and Task Columns to Integer Columns\ndef convert_valid_and_t(df):\n    df = df.with_columns(\n        pl.col(\"Valid\").cast(pl.Int8).alias(\"Valid\")\n    )\n    \n    df = df.with_columns(\n        pl.col(\"Task\").cast(pl.Int8).alias(\"Task\")\n    )\n    return df\ndefog_df = convert_valid_and_t(defog_df)\n# tdcsfog_df = convert_valid_and_t(tdcsfog_df)\n\n\nprint(defog_df)"},{"cell_type":"code","execution_count":null,"id":"ba85eb8e","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:57.26725Z","iopub.status.busy":"2025-10-20T02:50:57.266864Z","iopub.status.idle":"2025-10-20T02:50:58.558578Z","shell.execute_reply":"2025-10-20T02:50:58.557587Z"},"papermill":{"duration":1.303802,"end_time":"2025-10-20T02:50:58.56016","exception":false,"start_time":"2025-10-20T02:50:57.256358","status":"completed"},"tags":[]},"outputs":[],"source":"# Create a new column that contains the acceleration magnitude\ndef acc_magnitude(df):\n    df = df.with_columns(\n        (\n            (pl.col(\"AccV\") ** 2 + pl.col(\"AccML\") ** 2 + pl.col(\"AccAP\") ** 2).sqrt()\n        ).alias(\"Acc_MAGNITUDE\")\n    )\n\n    return df\n\ntdcsfog_df = acc_magnitude(tdcsfog_df)\ndefog_df = acc_magnitude(defog_df)\ndefog_df"},{"cell_type":"code","execution_count":null,"id":"78ba2ab1","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:50:58.639133Z","iopub.status.busy":"2025-10-20T02:50:58.638833Z","iopub.status.idle":"2025-10-20T02:51:04.043533Z","shell.execute_reply":"2025-10-20T02:51:04.042682Z"},"papermill":{"duration":5.474575,"end_time":"2025-10-20T02:51:04.04541","exception":false,"start_time":"2025-10-20T02:50:58.570835","status":"completed"},"tags":[]},"outputs":[],"source":"# Standardize acceleration per patient for each training dataframe\ndef standardize_acc_by_patient(df: pl.DataFrame):\n    acc_cols = ['AccV', 'AccML', 'AccAP']\n    for col in acc_cols:\n        df = df.with_columns(\n            (\n                (pl.col(col) - pl.col(col).mean().over(\"patient_id\")) /\n                pl.col(col).std().over(\"patient_id\")\n            ).alias(col)  # overwrite original column\n        )\n    return df\n\ntdcsfog_df = standardize_acc_by_patient(tdcsfog_df)\ndefog_df = standardize_acc_by_patient(defog_df)\ndefog_df"},{"cell_type":"code","execution_count":null,"id":"798ec8c7","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:51:04.070257Z","iopub.status.busy":"2025-10-20T02:51:04.069899Z","iopub.status.idle":"2025-10-20T02:51:04.079677Z","shell.execute_reply":"2025-10-20T02:51:04.078801Z"},"papermill":{"duration":0.024552,"end_time":"2025-10-20T02:51:04.081575","exception":false,"start_time":"2025-10-20T02:51:04.057023","status":"completed"},"tags":[]},"outputs":[],"source":"# Band-pass Filter \ndef infer_fs(time_seconds: np.ndarray) -> float:\n    dt = np.diff(np.asarray(time_seconds, dtype=float))\n    dt = dt[np.isfinite(dt) & (dt > 0)]\n    if dt.size == 0:\n        raise ValueError(\"Cannot infer sampling frequency from Time column.\")\n    return 1.0 / np.median(dt)\n\ndef design_bandpass(low_hz: float, high_hz: float, fs: float, order: int = 4):\n    nyq = fs / 2.0\n    low = max(1e-6, low_hz / nyq)\n    high = min(0.999999, high_hz / nyq)\n    if not (0 < low < high < 1):\n        raise ValueError(f\"Invalid band for fs={fs:.3f}Hz: low={low_hz}Hz, high={high_hz}Hz\")\n    b, a = butter(order, [low, high], btype=\"band\")\n    return b, a\n\ndef bandpass_series(y: pd.Series, b, a) -> np.ndarray:\n    sig = pd.to_numeric(y, errors=\"coerce\").interpolate(limit_direction=\"both\").to_numpy(float)\n    return filtfilt(b, a, sig, method=\"pad\")\n\ndef bandpass_dataframe(df: pd.DataFrame, cols=('AccV','AccML','AccAP'),\n                       low_hz=0.1, high_hz=30.0, order=4) -> pd.DataFrame:\n    out = df.copy()\n    # Only keep columns that exist\n    cols = tuple([c for c in cols if c in out.columns])\n    if len(cols) == 0:\n        return out\n\n    fs = infer_fs(out['Time'].to_numpy())\n    b, a = design_bandpass(low_hz, high_hz, fs, order)\n    for col in cols:\n        out[f\"{col}_bp\"] = bandpass_series(out[col], b, a)\n    return out\n\n"},{"cell_type":"code","execution_count":null,"id":"51241fc8","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:51:04.105432Z","iopub.status.busy":"2025-10-20T02:51:04.10513Z","iopub.status.idle":"2025-10-20T02:54:37.619172Z","shell.execute_reply":"2025-10-20T02:54:37.618223Z"},"papermill":{"duration":213.538407,"end_time":"2025-10-20T02:54:37.630332","exception":false,"start_time":"2025-10-20T02:51:04.091925","status":"completed"},"tags":[]},"outputs":[],"source":"# Apply Band-pass to all patients \ndef add_bandpass_to_all_patients(pl_df: pl.DataFrame,\n                                 cols=('AccV','AccML','AccAP'),\n                                 low_hz=0.1, high_hz=30.0, order=4) -> pl.DataFrame:\n    if \"patient_id\" not in pl_df.columns:\n        raise ValueError(\"Expected a 'patient_id' column.\")\n\n    out_chunks = []\n    # Unique patient list\n    patient_ids = pl_df.select(\"patient_id\").unique().to_series().to_list()\n\n    for pid in patient_ids:\n        g = pl_df.filter(pl.col(\"patient_id\") == pid).to_pandas()\n        # Skip tiny or malformed groups\n        if \"Time\" not in g.columns or len(g) < 5:\n            out_chunks.append(pl.from_pandas(g))  # just pass through\n            continue\n\n        try:\n            g_bp = bandpass_dataframe(g, cols=cols, low_hz=low_hz, high_hz=high_hz, order=order)\n        except Exception as e:\n            print(f\"[WARN] Skipping bandpass for patient {pid}: {e}\")\n            g_bp = g  # pass through raw if something fails\n\n        out_chunks.append(pl.from_pandas(g_bp))\n\n    return pl.concat(out_chunks, how=\"vertical_relaxed\")\n\ndefog_df_bp   = add_bandpass_to_all_patients(defog_df,   cols=('AccV','AccML','AccAP'),\n                                             low_hz=0.1, high_hz=30.0, order=4)\ntdcsfog_df_bp = add_bandpass_to_all_patients(tdcsfog_df, cols=('AccV','AccML','AccAP'),\n                                             low_hz=0.1, high_hz=30.0, order=4)\n\nprint(\"DEFOG with band-pass columns:\", [c for c in defog_df_bp.columns if c.endswith(\"_bp\")][:6], \"...\")\nprint(\"TDCSFOG with band-pass columns:\", [c for c in tdcsfog_df_bp.columns if c.endswith(\"_bp\")][:6], \"...\")"},{"cell_type":"markdown","id":"ebe15917","metadata":{"papermill":{"duration":0.010129,"end_time":"2025-10-20T02:54:37.6506","exception":false,"start_time":"2025-10-20T02:54:37.640471","status":"completed"},"tags":[]},"source":"# Plot Magnitude "},{"cell_type":"code","execution_count":null,"id":"082a56e8","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:37.673471Z","iopub.status.busy":"2025-10-20T02:54:37.673133Z","iopub.status.idle":"2025-10-20T02:54:40.328239Z","shell.execute_reply":"2025-10-20T02:54:40.327308Z"},"papermill":{"duration":2.668145,"end_time":"2025-10-20T02:54:40.329612","exception":false,"start_time":"2025-10-20T02:54:37.661467","status":"completed"},"tags":[]},"outputs":[],"source":"def add_magnitude_cols(pl_df: pl.DataFrame) -> pl.DataFrame: \n    out = pl_df.with_columns(\n        ((pl.col(\"AccV\")**2 + pl.col(\"AccML\")**2 + pl.col(\"AccAP\")**2).sqrt()).alias(\"AccMag\")\n    )\n    bp_cols = {\"AccV_bp\", \"AccML_bp\", \"AccAP_bp\"}\n    if bp_cols.issubset(set(out.columns)):\n        out = out.with_columns(\n            ((pl.col(\"AccV_bp\")**2 + pl.col(\"AccML_bp\")**2 + pl.col(\"AccAP_bp\")**2).sqrt()).alias(\"AccMag_bp\")\n        )\n    return out \n    \ndef plot_patient_mag(pl_df: pl.DataFrame, patient_id: str,\n                     time_col: str = \"Time\",\n                     show_events: bool = True,\n                     title_suffix: str = \"\"):\n    dfp = pl_df.filter(pl.col(\"patient_id\") == patient_id).to_pandas()\n    if not {\"AccMag\", time_col}.issubset(dfp.columns):\n        raise ValueError(\"Missing magnitude or time column\")\n    plt.figure(figsize=(16,6))\n    plt.plot(dfp[time_col], dfp[\"AccMag\"], label=\"|a| (raw)\", alpha=0.7)\n    if \"AccMag_bp\" in dfp.columns:\n        plt.plot(dfp[time_col], dfp[\"AccMag_bp\"], label=\"|a| (0.1–30 Hz)\", linewidth=1.6)\n    if show_events:\n        for ev in [\"StartHesitation\", \"Turn\", \"Walking\"]:\n            if ev in dfp.columns:\n                plt.plot(dfp[time_col], dfp[ev], label=ev, alpha=0.5)\n    plt.xlabel(\"Time (s)\")\n    plt.ylabel(\"Acceleration magnitude\")\n    plt.title(f\"Patient {patient_id} – Acc Magnitude {title_suffix}\")\n    plt.legend(ncol=3)\n    plt.grid(True)\n    plt.tight_layout()\n    plt.show()\n\n# Add magnitude columns to dataframes\ndefog_df_bp = add_magnitude_cols(defog_df_bp)\ntdcsfog_df_bp = add_magnitude_cols(tdcsfog_df_bp)\ndefog_df = add_magnitude_cols(defog_df)\ntdcsfog_df = add_magnitude_cols(tdcsfog_df)\n\n# Plot magnitude for one patient\nplot_patient_mag(defog_df_bp, patient_id=\"4c3aa8ea6e\", title_suffix=\"(raw vs band-pass)\")\n"},{"cell_type":"code","execution_count":null,"id":"7956a916","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:40.355756Z","iopub.status.busy":"2025-10-20T02:54:40.35544Z","iopub.status.idle":"2025-10-20T02:54:40.557397Z","shell.execute_reply":"2025-10-20T02:54:40.55633Z"},"papermill":{"duration":0.217194,"end_time":"2025-10-20T02:54:40.559224","exception":false,"start_time":"2025-10-20T02:54:40.34203","status":"completed"},"tags":[]},"outputs":[],"source":"# Create a new column that contains Time as seconds\ndef time_to_seconds(df, hertz):\n    df = df.with_columns(\n        (\n            (pl.col(\"Time\") / hertz)\n        ).alias(\"Time (seconds)\")\n    )\n\n    return df\n\ntdcsfog_df = time_to_seconds(tdcsfog_df, 128)\ndefog_df = time_to_seconds(tdcsfog_df, 100)\ndefog_df"},{"cell_type":"code","execution_count":null,"id":"1ee86132","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:40.584847Z","iopub.status.busy":"2025-10-20T02:54:40.584342Z","iopub.status.idle":"2025-10-20T02:54:40.731446Z","shell.execute_reply":"2025-10-20T02:54:40.730575Z"},"papermill":{"duration":0.161858,"end_time":"2025-10-20T02:54:40.733118","exception":false,"start_time":"2025-10-20T02:54:40.57126","status":"completed"},"tags":[]},"outputs":[],"source":"# Check for outliers from acceleration\ndef detect_outliers(df: pl.DataFrame):\n    acc_cols = ['AccV', 'AccML', 'AccAP']\n    for col in acc_cols: \n        z_col = col\n        outlier_df = df.filter(pl.col(z_col).abs() > 3.0)\n    return outlier_df\nprint(detect_outliers(defog_df))\nprint(detect_outliers(tdcsfog_df))"},{"cell_type":"markdown","id":"bbbd2ee4","metadata":{"papermill":{"duration":0.012372,"end_time":"2025-10-20T02:54:40.75728","exception":false,"start_time":"2025-10-20T02:54:40.744908","status":"completed"},"tags":[]},"source":"## Visualize Acceleration  Signals During FoG Events"},{"cell_type":"code","execution_count":null,"id":"09af5eca","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:40.783017Z","iopub.status.busy":"2025-10-20T02:54:40.78269Z","iopub.status.idle":"2025-10-20T02:54:40.940145Z","shell.execute_reply":"2025-10-20T02:54:40.938752Z"},"papermill":{"duration":0.172712,"end_time":"2025-10-20T02:54:40.942656","exception":false,"start_time":"2025-10-20T02:54:40.769944","status":"completed"},"tags":[]},"outputs":[],"source":"# Get unique patient IDs with a StartHesitation, Turn, and Walking event\n# Take a subset of 3 patients for each event\nStartHesPatients = (\n    defog_df.filter(pl.col(\"StartHesitation\") == 1)\n            .select(\"patient_id\")\n            .unique()\n            .to_series()[:3]  # take first 3\n)\nprint(f\"Patients with Start Hesitation: {StartHesPatients.to_list()}\")\n\nTurnPatients = (\n    defog_df.filter(pl.col(\"Turn\") == 1)\n            .select(\"patient_id\")\n            .unique()\n            .to_series()[:3]\n)\nprint(f\"Patients with Turn: {TurnPatients.to_list()}\")\n\nWalkingPatients = (\n    defog_df.filter(pl.col(\"Walking\") == 1)\n            .select(\"patient_id\")\n            .unique()\n            .to_series()[:3]\n)\nprint(f\"Patients with Walking: {WalkingPatients.to_list()}\")"},{"cell_type":"code","execution_count":null,"id":"10e3a488","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:40.970371Z","iopub.status.busy":"2025-10-20T02:54:40.970055Z","iopub.status.idle":"2025-10-20T02:54:41.046664Z","shell.execute_reply":"2025-10-20T02:54:41.045435Z"},"papermill":{"duration":0.093121,"end_time":"2025-10-20T02:54:41.048683","exception":false,"start_time":"2025-10-20T02:54:40.955562","status":"completed"},"tags":[]},"outputs":[],"source":"# Get unique patient IDs with a StartHesitation, Turn, and Walking event \n# (including band-pass)\nif {\"StartHesitation\",\"Turn\",\"Walking\"}.issubset(set(defog_df_bp.columns)):\n    StartHesPatients = (\n        defog_df_bp.filter(pl.col(\"StartHesitation\") == 1)\n                   .select(\"patient_id\").unique().to_series()[:3]\n    )\n    TurnPatients = (\n        defog_df_bp.filter(pl.col(\"Turn\") == 1)\n                   .select(\"patient_id\").unique().to_series()[:3]\n    )\n    WalkingPatients = (\n        defog_df_bp.filter(pl.col(\"Walking\") == 1)\n                   .select(\"patient_id\").unique().to_series()[:3]\n    )\n    print(f\"Patients with Start Hesitation: {StartHesPatients.to_list()}\")\n    print(f\"Patients with Turn: {TurnPatients.to_list()}\")\n    print(f\"Patients with Walking: {WalkingPatients.to_list()}\")"},{"cell_type":"code","execution_count":null,"id":"4676d34b","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:41.07781Z","iopub.status.busy":"2025-10-20T02:54:41.077481Z","iopub.status.idle":"2025-10-20T02:54:42.176441Z","shell.execute_reply":"2025-10-20T02:54:42.175448Z"},"papermill":{"duration":1.117996,"end_time":"2025-10-20T02:54:42.179172","exception":false,"start_time":"2025-10-20T02:54:41.061176","status":"completed"},"tags":[]},"outputs":[],"source":"# Start Hestitation\n# 1. Filter your Polars DF for a single patient and convert to pandas\ndf = defog_df.filter(pl.col(\"patient_id\") == '81262644e7').to_pandas()\n\n# 2. Plot\nplt.figure(figsize=(15, 6))\n\n# Plot acceleration\nplt.plot(df['Time'], df['AccV'], label='AccV', alpha=0.7)\nplt.plot(df['Time'], df['AccML'], label='AccML', alpha=0.7)\nplt.plot(df['Time'], df['AccAP'], label='AccAP', alpha=0.7)\n\n# 3. Plot events\nplt.plot(df['Time'], df['StartHesitation'], label='StartHesitation', alpha=0.7)\n\n\n# 4. Final touches\nplt.xlabel(\"Time\")\nplt.ylabel(\"Acceleration (g)\")\nplt.title(f\"Patient: {patient_id} - Acceleration + FOG Events\")\nplt.legend(loc=\"upper right\")\nplt.grid(True)\nplt.tight_layout()\nplt.show()\n\n\n# Start Hestitation\n# 1. Filter your Polars DF for a single patient and convert to pandas\ndf = defog_df.filter(pl.col(\"patient_id\") == '3ba3590a08').to_pandas()\n\n# 2. Plot\nplt.figure(figsize=(15, 6))\n\n# Plot acceleration\nplt.plot(df['Time'], df['AccV'], label='AccV', alpha=0.7)\nplt.plot(df['Time'], df['AccML'], label='AccML', alpha=0.7)\nplt.plot(df['Time'], df['AccAP'], label='AccAP', alpha=0.7)\n\n# 3. Plot events\nplt.plot(df['Time'], df['StartHesitation'], label='StartHesitation', alpha=0.7)\n\n\n# 4. Final touches\nplt.xlabel(\"Time\")\nplt.ylabel(\"Acceleration (g)\")\nplt.title(f\"Patient: {patient_id} - Acceleration + FOG Events\")\nplt.legend(loc=\"upper right\")\nplt.grid(True)\nplt.tight_layout()\nplt.show()\n\n\n\n# Start Hestitation\n# 1. Filter your Polars DF for a single patient and convert to pandas\ndf = defog_df.filter(pl.col(\"patient_id\") == 'd98358a75f').to_pandas()\n\n# 2. Plot\nplt.figure(figsize=(15, 6))\n\n# Plot acceleration\nplt.plot(df['Time'], df['AccV'], label='AccV', alpha=0.7)\nplt.plot(df['Time'], df['AccML'], label='AccML', alpha=0.7)\nplt.plot(df['Time'], df['AccAP'], label='AccAP', alpha=0.7)\n\n# 3. Plot events\nplt.plot(df['Time'], df['StartHesitation'], label='StartHesitation', alpha=0.7)\n\n\n# 4. Final touches\nplt.xlabel(\"Time\")\nplt.ylabel(\"Acceleration (g)\")\nplt.title(f\"Patient: {patient_id} - Acceleration + FOG Events\")\nplt.legend(loc=\"upper right\")\nplt.grid(True)\nplt.tight_layout()\nplt.show()"},{"cell_type":"markdown","id":"a78410d7","metadata":{"papermill":{"duration":0.019595,"end_time":"2025-10-20T02:54:42.218902","exception":false,"start_time":"2025-10-20T02:54:42.199307","status":"completed"},"tags":[]},"source":"# Extract Time Domain Features "},{"cell_type":"code","execution_count":null,"id":"5c8efcb7","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:42.259065Z","iopub.status.busy":"2025-10-20T02:54:42.258722Z","iopub.status.idle":"2025-10-20T02:54:42.277678Z","shell.execute_reply":"2025-10-20T02:54:42.276703Z"},"papermill":{"duration":0.04084,"end_time":"2025-10-20T02:54:42.279144","exception":false,"start_time":"2025-10-20T02:54:42.238304","status":"completed"},"tags":[]},"outputs":[],"source":"def extract_time_features_pd(\n    df: pd.DataFrame,\n    fs: float,\n    win_s: float,\n    hop_s: float,\n    signal_cols=('AccX','AccY','AccZ','GyroX','GyroY','GyroZ'),\n    label_cols=None,                 # e.g. ['StartHesitation','Turn','Walking','Event']\n    id_cols=None,                    # e.g. ['subject_id','series_id','id'] to carry through\n):\n    df = df.copy()\n    n = len(df)\n    win = int(round(fs*win_s))\n    hop = int(round(fs*hop_s))\n    if win <= 0 or hop <= 0:\n        raise ValueError(\"win_s and hop_s must be > 0\")\n\n    # Pre-pull arrays for speed\n    X = df.loc[:, signal_cols].to_numpy(dtype=float)\n\n    # Helper feature fns (safe on NaNs/empties)\n    def feats_one(w):\n        f = {}\n        # basic stats\n        f.update({f'{c}_mean': np.nanmean(w[:,i]) for i,c in enumerate(signal_cols)})\n        f.update({f'{c}_var' : np.nanvar (w[:,i]) for i,c in enumerate(signal_cols)})\n        f.update({f'{c}_std' : np.nanstd (w[:,i]) for i,c in enumerate(signal_cols)})\n        f.update({f'{c}_min' : np.nanmin (w[:,i]) for i,c in enumerate(signal_cols)})\n        f.update({f'{c}_max' : np.nanmax (w[:,i]) for i,c in enumerate(signal_cols)})\n        f.update({f'{c}_median': np.nanmedian(w[:,i]) for i,c in enumerate(signal_cols)})\n        f.update({f'{c}_iqr': np.nanpercentile(w[:,i],75)-np.nanpercentile(w[:,i],25) for i,c in enumerate(signal_cols)})\n\n        # energy & rms\n        f.update({f'{c}_energy': np.nansum(np.square(w[:,i]))/len(w) for i,c in enumerate(signal_cols)})\n        f.update({f'{c}_rms'   : np.sqrt(np.nanmean(np.square(w[:,i]))) for i,c in enumerate(signal_cols)})\n\n        # skew & kurt\n        for i,c in enumerate(signal_cols):\n            col = w[:,i]\n            f[f'{c}_skew'] = skew(col, nan_policy='omit', bias=False)\n            f[f'{c}_kurt'] = kurtosis(col, nan_policy='omit', fisher=True, bias=False)\n\n        # vector features for tri-axial groups\n        if set(['AccX','AccY','AccZ']).issubset(signal_cols):\n            ax = [signal_cols.index('AccX'), signal_cols.index('AccY'), signal_cols.index('AccZ')]\n            acc = w[:,ax]\n            mag = np.sqrt(np.sum(acc**2, axis=1))\n            f['Acc_mag_mean'] = np.nanmean(mag)\n            f['Acc_mag_std']  = np.nanstd(mag)\n            # Signal Magnitude Area (SMA)\n            f['Acc_sma'] = (np.nansum(np.abs(acc), axis=0).sum()) / len(mag)\n\n            # correlations\n            for (a,b) in [('AccX','AccY'),('AccX','AccZ'),('AccY','AccZ')]:\n                i1, i2 = signal_cols.index(a), signal_cols.index(b)\n                col1, col2 = w[:,i1], w[:,i2]\n                if np.all(np.isfinite(col1)) and np.all(np.isfinite(col2)) and len(col1) > 1:\n                    f[f'corr_{a}_{b}'] = np.corrcoef(col1, col2)[0,1]\n                else:\n                    f[f'corr_{a}_{b}'] = np.nan\n\n        if set(['GyroX','GyroY','GyroZ']).issubset(signal_cols):\n            gx = [signal_cols.index('GyroX'), signal_cols.index('GyroY'), signal_cols.index('GyroZ')]\n            gyro = w[:,gx]\n            mag = np.sqrt(np.sum(gyro**2, axis=1))\n            f['Gyro_mag_mean'] = np.nanmean(mag)\n            f['Gyro_mag_std']  = np.nanstd(mag)\n            f['Gyro_sma'] = (np.nansum(np.abs(gyro), axis=0).sum()) / len(mag)\n            for (a,b) in [('GyroX','GyroY'),('GyroX','GyroZ'),('GyroY','GyroZ')]:\n                i1, i2 = signal_cols.index(a), signal_cols.index(b)\n                col1, col2 = w[:,i1], w[:,i2]\n                if np.all(np.isfinite(col1)) and np.all(np.isfinite(col2)) and len(col1) > 1:\n                    f[f'corr_{a}_{b}'] = np.corrcoef(col1, col2)[0,1]\n                else:\n                    f[f'corr_{a}_{b}'] = np.nan\n\n        return f\n\n    rows = []\n    start_idx = 0\n    win_id = 0\n    while start_idx + win <= n:\n        end_idx = start_idx + win\n        w = X[start_idx:end_idx, :]\n        feat = feats_one(w)\n        # add time/window metadata\n        feat['win_id'] = win_id\n        feat['t_start_s'] = start_idx / fs\n        feat['t_end_s']   = (end_idx-1) / fs\n\n        # bring-through IDs from the *end* row of the window (common in DEFoG baselines)\n        if id_cols:\n            for c in id_cols:\n                feat[c] = df.iloc[end_idx-1][c]\n\n        # aggregate labels if provided (ANY>0 inside window)\n        if label_cols:\n            sub = df.iloc[start_idx:end_idx]\n            feat['label_any'] = bool((sub[label_cols].fillna(0).to_numpy() > 0).any())\n            for c in label_cols:\n                feat[f'label_{c}'] = bool((sub[c].fillna(0).to_numpy() > 0).any())\n\n        rows.append(feat)\n        win_id += 1\n        start_idx += hop\n\n    return pd.DataFrame(rows)"},{"cell_type":"code","execution_count":null,"id":"f1762c5d","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:42.318952Z","iopub.status.busy":"2025-10-20T02:54:42.318681Z","iopub.status.idle":"2025-10-20T02:54:42.328438Z","shell.execute_reply":"2025-10-20T02:54:42.327435Z"},"papermill":{"duration":0.031302,"end_time":"2025-10-20T02:54:42.33023","exception":false,"start_time":"2025-10-20T02:54:42.298928","status":"completed"},"tags":[]},"outputs":[],"source":"DATA_DIR = Path(\"/kaggle/input/tlvmc-parkinsons-freezing-gait-prediction/train\")\n\nprint(\"Available training folders:\", [p.name for p in DATA_DIR.iterdir()])"},{"cell_type":"code","execution_count":null,"id":"ca67a082","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:42.371095Z","iopub.status.busy":"2025-10-20T02:54:42.370778Z","iopub.status.idle":"2025-10-20T02:54:42.388215Z","shell.execute_reply":"2025-10-20T02:54:42.387185Z"},"papermill":{"duration":0.03992,"end_time":"2025-10-20T02:54:42.38979","exception":false,"start_time":"2025-10-20T02:54:42.34987","status":"completed"},"tags":[]},"outputs":[],"source":"tdcsfog_file = list((DATA_DIR / \"tdcsfog\").glob(\"*.csv\"))[0]\ndefog_file   = list((DATA_DIR / \"defog\").glob(\"*.csv\"))[0]\nnotype_file  = list((DATA_DIR / \"notype\").glob(\"*.csv\"))[0]\n\nprint(\"Sample files chosen:\")\nprint(\"tdcsfog:\", tdcsfog_file.name)\nprint(\"defog:\", defog_file.name)\nprint(\"notype:\", notype_file.name)"},{"cell_type":"code","execution_count":null,"id":"16a4c05b","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:42.430448Z","iopub.status.busy":"2025-10-20T02:54:42.430155Z","iopub.status.idle":"2025-10-20T02:54:42.839368Z","shell.execute_reply":"2025-10-20T02:54:42.838637Z"},"papermill":{"duration":0.431683,"end_time":"2025-10-20T02:54:42.841107","exception":false,"start_time":"2025-10-20T02:54:42.409424","status":"completed"},"tags":[]},"outputs":[],"source":"\ntdcsfog_df = pd.read_csv(tdcsfog_file)\ndefog_df   = pd.read_csv(defog_file)\nnotype_df  = pd.read_csv(notype_file)"},{"cell_type":"code","execution_count":null,"id":"f8ae1d6f","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:42.884132Z","iopub.status.busy":"2025-10-20T02:54:42.883841Z","iopub.status.idle":"2025-10-20T02:54:42.899592Z","shell.execute_reply":"2025-10-20T02:54:42.8986Z"},"papermill":{"duration":0.039881,"end_time":"2025-10-20T02:54:42.901694","exception":false,"start_time":"2025-10-20T02:54:42.861813","status":"completed"},"tags":[]},"outputs":[],"source":"import re\n\ndef tdcsfog_time_features(df, fs=128.0, win_s=2.0, hop_s=0.5):\n    \"\"\"\n    Extract time-domain features from accelerometer data in tdcsfog_df.\n    Works purely in pandas.\n    \"\"\"\n\n    # --- auto-detect accelerometer columns ---\n    cols = list(df.columns)\n    lower = {c.lower(): c for c in cols}\n    candidates = [\n        ['AccX','AccY','AccZ'],\n        ['AccelX','AccelY','AccelZ'],\n        ['acc_x','acc_y','acc_z'],\n        ['AccV','AccML','AccAP'],\n        ['accv','accml','accap'],\n    ]\n    acc_cols = None\n    for trio in candidates:\n        found = [lower.get(c.lower()) for c in trio]\n        if all(found):\n            acc_cols = found\n            break\n    if acc_cols is None:\n        acc_cols = [c for c in cols if re.search('acc', c, re.I)][:3]\n    if len(acc_cols) < 3:\n        raise KeyError(f\"Could not find 3 accelerometer columns. Found: {acc_cols}\")\n\n    # --- label columns ---\n    label_cols = [c for c in ['StartHesitation','Turn','Walking'] if c in df.columns]\n\n    # --- window setup ---\n    win = int(round(fs * win_s))\n    hop = int(round(fs * hop_s))\n    n = len(df)\n\n    X = df[acc_cols].to_numpy(dtype=float)\n    rows = []\n    start = 0\n    win_id = 0\n\n    while start + win <= n:\n        end = start + win\n        W = X[start:end, :]\n        f = {}\n\n        # per-axis features\n        for i, c in enumerate(acc_cols):\n            w = W[:, i]\n            f[f'{c}_mean']   = np.nanmean(w)\n            f[f'{c}_std']    = np.nanstd(w)\n            f[f'{c}_var']    = np.nanvar(w)\n            f[f'{c}_min']    = np.nanmin(w)\n            f[f'{c}_max']    = np.nanmax(w)\n            f[f'{c}_median'] = np.nanmedian(w)\n            q75, q25 = np.nanpercentile(w, [75, 25])\n            f[f'{c}_iqr']    = q75 - q25\n            f[f'{c}_energy'] = np.nansum(w**2) / len(w)\n            f[f'{c}_rms']    = np.sqrt(np.nanmean(w**2))\n            f[f'{c}_skew']   = skew(w, nan_policy='omit', bias=False)\n            f[f'{c}_kurt']   = kurtosis(w, nan_policy='omit', fisher=True, bias=False)\n\n        # vector magnitude features\n        mag = np.sqrt(np.sum(W**2, axis=1))\n        f['Acc_mag_mean'] = np.nanmean(mag)\n        f['Acc_mag_std']  = np.nanstd(mag)\n        f['Acc_sma']      = np.nansum(np.abs(W)) / len(W)\n\n        # correlations\n        def corr_safe(a,b):\n            if len(a) > 1 and np.isfinite(a).all() and np.isfinite(b).all():\n                return np.corrcoef(a,b)[0,1]\n            return np.nan\n        f[f'corr_{acc_cols[0]}_{acc_cols[1]}'] = corr_safe(W[:,0], W[:,1])\n        f[f'corr_{acc_cols[0]}_{acc_cols[2]}'] = corr_safe(W[:,0], W[:,2])\n        f[f'corr_{acc_cols[1]}_{acc_cols[2]}'] = corr_safe(W[:,1], W[:,2])\n\n        # labels (if exist)\n        if label_cols:\n            sub = df.iloc[start:end][label_cols].fillna(0).to_numpy()\n            f['label_any'] = bool((sub > 0).any())\n            for j, c in enumerate(label_cols):\n                f[f'label_{c}'] = bool((sub[:, j] > 0).any())\n\n        # window metadata\n        f['win_id'] = win_id\n        f['t_start_s'] = start / fs\n        f['t_end_s']   = (end - 1) / fs\n\n        rows.append(f)\n        win_id += 1\n        start += hop\n\n    return pd.DataFrame(rows)"},{"cell_type":"code","execution_count":null,"id":"19b2f73e","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:42.944202Z","iopub.status.busy":"2025-10-20T02:54:42.94323Z","iopub.status.idle":"2025-10-20T02:54:43.612622Z","shell.execute_reply":"2025-10-20T02:54:43.611813Z"},"papermill":{"duration":0.692153,"end_time":"2025-10-20T02:54:43.614382","exception":false,"start_time":"2025-10-20T02:54:42.922229","status":"completed"},"tags":[]},"outputs":[],"source":"# parameters\nFS = 128.0\nWIN = 2.0   # 2 seconds per window\nHOP = 0.5   # 0.5-second step\n\ntdcsfog_feats = tdcsfog_time_features(tdcsfog_df, fs=FS, win_s=WIN, hop_s=HOP)\n\nprint(\"Detected accelerometer columns:\", [c for c in tdcsfog_feats.columns if 'mean' in c][:3])\nprint(\"Shape:\", tdcsfog_feats.shape)\ntdcsfog_feats.head()"},{"cell_type":"markdown","id":"6b77b4cd","metadata":{"papermill":{"duration":0.019569,"end_time":"2025-10-20T02:54:43.654934","exception":false,"start_time":"2025-10-20T02:54:43.635365","status":"completed"},"tags":[]},"source":"# # Extract Frequency Domain Features \n*  The Fourier Transform is a mathematical tool that takes a complex signal and breaks it down into its individual frequency components. It tells us what frequencies make up the signal and how strong each one is.\n*  Time Domain features (mean, std, min, max, variance, median) shoes how a signal changes over time while a frequency domain feature (domiant frequency, spectral enegy, PSD) shows what frequencies make up the signal and how fast or periodic the motion is.\n*  This helps us analyze how the signal's energy is distrubuted across different frequencies, revealing rhytmic motion patterns like walking, turning, or freezing of gait.\n*  For PSD, we utilized Welch's method which smooths the FFT to give a more stable estimate of singal power over frequency. It does a great job detaching dominant movement frequencies (like steps per second)\n\nManual Frequency Features: \n* Directly linked to Parkinson's biomarkers like tremor, gait frequency, and freezing index\n* Ability to tweak freqency bands, window sizes, or filtering parameters\n  Quick to compute for large datasets "},{"cell_type":"code","execution_count":null,"id":"8e9d87c4","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:43.695813Z","iopub.status.busy":"2025-10-20T02:54:43.695467Z","iopub.status.idle":"2025-10-20T02:54:43.707751Z","shell.execute_reply":"2025-10-20T02:54:43.706863Z"},"papermill":{"duration":0.034348,"end_time":"2025-10-20T02:54:43.709173","exception":false,"start_time":"2025-10-20T02:54:43.674825","status":"completed"},"tags":[]},"outputs":[],"source":"from scipy.fft import rfft, rfftfreq      # Fast Fourier Transform functions + corresponding frequency bins for real-valued signals \nfrom scipy.signal import welch            # Power Spectral Density function (shows how signal's power varies across frequencies)\n\ndef extract_frequency_features(df, fs=128.0, win_s=2.0, hop_s=0.5, signal_cols=(('AccX','AccY','AccZ'))):\n    \"\"\"\n    Extracts frequency domain features from accelerometer data\n    Inputs: \n    - df: pandas DataFrame that has time-series sensor data \n    - fs: sampling frequency in Hz (how many readings per second)\n    - win_s: window size in seconds \n    - hop_s: step size in seconds (how far we slide each window)\n    - signal_cols: columns to use (like accelerometer axes)\n    \"\"\"\n\n    df = df.copy()                 # make a copy so we don't modify the original datafram e\n    n = len(df)                    # number of total samples (rows)\n    win = int(round(fs* win_s))    # window size in samples \n    hop = int(round(fs * hop_s))   # hop size in samples\n    rows = []                      # to store feature dictionaries for each window \n    start = 0                      # start index for the first window \n    win_id = 0                     # window counter \n\n    # Loop through the data using sliding windows \n    while start + win <= n:\n        end = start + win\n        # get window segment (subset of rows)\n        segment = df.iloc[start:end][list(signal_cols)].to_numpy(dtype=float)\n\n        # prepare frquency bins (x-axis for FFT)\n        freqs = rfftfreq(win, d=1/fs)      # rfft only returns positive frequencies \n\n        fdict = {}   # dictionary to store the feature for this window \n\n        # Loop through each sensor column (AccX, AccY, AccZ)\n        for i, col in enumerate(signal_cols): \n            signal = segment [:, i]               # extract one axis \n            signal = signal - np.nanmean(signal)  # center around 0 \n            signal = np.nan_to_num(signal)        # replace NaNs with 0s (avoid FFT errors)\n\n            # FFT: convert from time to frequency \n            fft_vals = np.abs(rfft(signal))           # absolute value of FFT (magnitude)\n            fft_power = fft_vals ** 2                 # power = magnitude squared \n\n            # PSD - Power Spectral Density \n            f_psd, psd = welch(signal, fs=fs, nperseg=min(win, len(signal)))\n\n            # Basic FFT & PSD feature summaries to help distinguish steady vs irregular movements \n            fdict[f'{col}_fft_mean'] = np.mean(fft_power)   # average FFT power\n            fdict[f'{col}_fft_max']  = np.max(fft_power)    # peak FFT power\n            fdict[f'{col}_fft_std']  = np.std(fft_power)    # variation in FFT power\n\n            fdict[f'{col}_psd_mean'] = np.mean(psd)         # average power across frequencies\n            fdict[f'{col}_psd_max']  = np.max(psd)          # maximum power (dominant peak)\n            fdict[f'{col}_psd_std']  = np.std(psd)          # variability in PSD\n\n            # Dominant Frequency \n            # frequency (in Hz) where PSD is largest\n            dom_freq = f_psd[np.argmax(psd)]\n            fdict[f'{col}_dominant_freq'] = dom_freq\n\n            # Band Power (energy in specific frequency ranges) \n            # Integrate (area under curve) power in low, mid, and high frequency ranges\n            fdict[f'{col}_bandpower_low']  = np.trapz(psd[(f_psd>=0)  & (f_psd<3)],  f_psd[(f_psd>=0)  & (f_psd<3)])   # 0–3 Hz (large body movements)\n            fdict[f'{col}_bandpower_mid']  = np.trapz(psd[(f_psd>=3)  & (f_psd<10)], f_psd[(f_psd>=3)  & (f_psd<10)])  # 3–10 Hz (normal walking frequency)\n            fdict[f'{col}_bandpower_high'] = np.trapz(psd[(f_psd>=10) & (f_psd<30)], f_psd[(f_psd>=10) & (f_psd<30)])  # 10–30 Hz (fine tremors/sudden changes)\n            \n        fdict['win_id'] = win_id                 # window number\n        fdict['t_start_s'] = start / fs          # start time (seconds)\n        fdict['t_end_s']   = (end - 1) / fs      # end time (seconds)\n\n        rows.append(fdict)    # store results\n        win_id += 1           # move to next window\n        start += hop          # slide window by hop length\n\n    # return a new DataFrame with all extracted frequency features\n    return pd.DataFrame(rows)"},{"cell_type":"code","execution_count":null,"id":"130252d8","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:43.749599Z","iopub.status.busy":"2025-10-20T02:54:43.749237Z","iopub.status.idle":"2025-10-20T02:54:43.755161Z","shell.execute_reply":"2025-10-20T02:54:43.754273Z"},"papermill":{"duration":0.027616,"end_time":"2025-10-20T02:54:43.756396","exception":false,"start_time":"2025-10-20T02:54:43.72878","status":"completed"},"tags":[]},"outputs":[],"source":"tdcsfog_df['patient_id'] = 'patient_1' \nprint(tdcsfog_df.columns)"},{"cell_type":"markdown","id":"528f1b6a","metadata":{"papermill":{"duration":0.01998,"end_time":"2025-10-20T02:54:43.796042","exception":false,"start_time":"2025-10-20T02:54:43.776062","status":"completed"},"tags":[]},"source":"# Using tsfresh "},{"cell_type":"markdown","id":"636eaba7","metadata":{"papermill":{"duration":0.019619,"end_time":"2025-10-20T02:54:43.835763","exception":false,"start_time":"2025-10-20T02:54:43.816144","status":"completed"},"tags":[]},"source":"**Frequency Domain: FFT coefficients, spectral entropy, energy ratios, etc.**\n* Time domain: mean, variance, autocorrelation, absolute energy, quantiles, etc.\n* Complex statistics: number of peaks, linear trend slopes, precentage of reoccuring datapoints, and more.\n  \n**that allows our model to**\n* Capture hidden micro-patterns in gait that's not visibilly obvious\n* Use automated feature selection later (tsfresh.select_features())\n* Complement your domain-driven features with new, data-driven ones - something revealing new insights "},{"cell_type":"code","execution_count":null,"id":"c27881de","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:43.875754Z","iopub.status.busy":"2025-10-20T02:54:43.875426Z","iopub.status.idle":"2025-10-20T02:54:43.884245Z","shell.execute_reply":"2025-10-20T02:54:43.883159Z"},"papermill":{"duration":0.0311,"end_time":"2025-10-20T02:54:43.885954","exception":false,"start_time":"2025-10-20T02:54:43.854854","status":"completed"},"tags":[]},"outputs":[],"source":"from tsfresh import extract_features\nfrom tsfresh.feature_extraction import EfficientFCParameters\n\n\ndef extract_tsfresh_features_windowed(\n    df, id_col='patient_id', time_col='Time',\n    signal_cols=('AccV','AccML','AccAP'),\n    fs=128.0, win_s=2.0\n):\n    \"\"\"\n    Extract tsfresh features per short window (fast and memory-safe).\n    Each patient signal is divided into small windows before feature extraction.\n    \"\"\"\n    win = int(fs * win_s)\n    df = df.copy()\n    all_feats = []\n\n    for pid in df[id_col].unique():\n        sub = df[df[id_col] == pid].copy().reset_index(drop=True)\n        n = len(sub)\n        print(f\"Extracting tsfresh features for {pid} ({n} samples)\")\n\n        # Create window IDs (one per window)\n        sub['window_id'] = (sub.index // win).astype(int)\n\n        # Melt into long format\n        df_long = sub.melt(\n            id_vars=[id_col, 'window_id', time_col],\n            value_vars=list(signal_cols),\n            var_name='sensor_axis',\n            value_name='value'\n        ).rename(columns={id_col: 'id', time_col: 'time'})\n\n        # Combine patient_id and window_id into one unique series id\n        df_long['id'] = df_long['id'].astype(str) + \"_\" + df_long['window_id'].astype(str)\n\n        # Clean numeric values\n        df_long['time']  = pd.to_numeric(df_long['time'], errors='coerce')\n        df_long['value'] = pd.to_numeric(df_long['value'], errors='coerce')\n        df_long = df_long.dropna(subset=['time', 'value']).sort_values(['id', 'time'])\n\n        # Extract light feature set\n        feats = extract_features(\n            df_long,\n            column_id='id',\n            column_sort='time',\n            column_value='value',\n            default_fc_parameters=EfficientFCParameters(),\n            n_jobs=0,\n            disable_progressbar=False\n        )\n\n        feats.reset_index(inplace=True)\n        feats[id_col] = pid\n        all_feats.append(feats)\n        print(f\"✓ Done {pid}\")\n\n    return pd.concat(all_feats, ignore_index=True) if all_feats else pd.DataFrame()"},{"cell_type":"code","execution_count":null,"id":"c30996f3","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:43.930924Z","iopub.status.busy":"2025-10-20T02:54:43.930421Z","iopub.status.idle":"2025-10-20T02:54:43.937469Z","shell.execute_reply":"2025-10-20T02:54:43.936342Z"},"papermill":{"duration":0.03008,"end_time":"2025-10-20T02:54:43.938886","exception":false,"start_time":"2025-10-20T02:54:43.908806","status":"completed"},"tags":[]},"outputs":[],"source":"def combine_features(manual_df, tsfresh_df):\n    \"\"\"\n    Merge manual and tsfresh features automatically, even if key columns differ.\n    \"\"\"\n    # Auto-detect join column\n    common_keys = set(manual_df.columns) & set(tsfresh_df.columns)\n    possible_keys = {'id', 'patient_id', 'window_id'}\n    join_key = list(common_keys & possible_keys)\n    \n    if not join_key:\n        print(\"[WARN] No shared key found. Using 'id' by default.\")\n        manual_df['id'] = manual_df.get('id', 'unknown')\n        tsfresh_df['id'] = tsfresh_df.get('id', 'unknown')\n        join_key = ['id']\n    else:\n        join_key = [join_key[0]]  # pick the first match\n    \n    print(f\"Merging on key: {join_key[0]}\")\n    merged = pd.merge(tsfresh_df, manual_df, how='inner', on=join_key)\n    return merged"},{"cell_type":"code","execution_count":null,"id":"285b4cf1","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:43.981898Z","iopub.status.busy":"2025-10-20T02:54:43.981624Z","iopub.status.idle":"2025-10-20T02:54:43.985339Z","shell.execute_reply":"2025-10-20T02:54:43.984572Z"},"papermill":{"duration":0.026425,"end_time":"2025-10-20T02:54:43.986671","exception":false,"start_time":"2025-10-20T02:54:43.960246","status":"completed"},"tags":[]},"outputs":[],"source":"# Sample testing\n\n# sample_df = tdcsfog_df.head(500).copy()\n# sample_df['patient_id'] = 'p_test'\n\n# tsfresh_feats = extract_tsfresh_features_windowed(\n    # sample_df,\n    # id_col='patient_id',\n    # time_col='Time',\n    # signal_cols=('AccV','AccML','AccAP'),\n    # fs=128.0,\n    # win_s=2.0\n# )"},{"cell_type":"code","execution_count":null,"id":"1fc33f09","metadata":{"execution":{"iopub.execute_input":"2025-10-20T02:54:44.028379Z","iopub.status.busy":"2025-10-20T02:54:44.027357Z","iopub.status.idle":"2025-10-20T02:54:48.757768Z","shell.execute_reply":"2025-10-20T02:54:48.756302Z"},"papermill":{"duration":4.75242,"end_time":"2025-10-20T02:54:48.759292","exception":false,"start_time":"2025-10-20T02:54:44.006872","status":"completed"},"tags":[]},"outputs":[],"source":"\n\nFS = 128.0\nWIN = 2.0\nHOP = 0.5\n\n# Manual frequency-domain features\nfreq_feats = extract_frequency_features(\n    tdcsfog_df,\n    fs=FS,\n    win_s=WIN,\n    hop_s=HOP,\n    signal_cols=('AccV','AccML','AccAP')\n)\n\n# Add patient_id to match tsfresh later\nfreq_feats['id'] = tdcsfog_df['patient_id'].iloc[0] if 'patient_id' in tdcsfog_df.columns else 'unknown'\n\nprint(\"Manual frequency features shape:\", freq_feats.shape)\n\n# tsfresh automatic features\ntsfresh_feats = extract_tsfresh_features_windowed(\n    tdcsfog_df,\n    id_col='patient_id',\n    time_col='Time',\n    signal_cols=('AccV','AccML','AccAP')\n)\nprint(\"TSFresh features shape:\", tsfresh_feats.shape)\n\n# Ensure both have lowercase 'id' column for merging\ntsfresh_feats = tsfresh_feats.rename(columns={'patient_id': 'id'})  # keep one id column\nfreq_feats = freq_feats.rename(columns={'patient_id': 'id'})        # unify naming just in case\n\n# Drop duplicates if needed\ntsfresh_feats = tsfresh_feats.loc[:, ~tsfresh_feats.columns.duplicated()]\n\n# Confirm alignment before merge\nprint(\"freq_feats ids:\", freq_feats['id'].head().tolist())\nprint(\"tsfresh_feats ids:\", tsfresh_feats['id'].head().tolist())\n\n# Merge again\nfeatures_df = combine_features(freq_feats, tsfresh_feats)\nprint(\"✅ Combined feature set shape:\", features_df.shape)"},{"cell_type":"markdown","id":"062c61d7","metadata":{},"source":"# NEW: PCA & LDA — Detailed, English-Only Addendum\n\nThis addendum extends the original notebook with **carefully annotated** sections on **Principal Component Analysis (PCA)** and **Linear Discriminant Analysis (LDA)**.  \nIt is designed to be ready-to-submit: conceptual notes, math, code, diagnostics, and clear plots (matplotlib only).\n"},{"cell_type":"markdown","id":"5fc525be","metadata":{},"source":"## 0. Environment\n\nWe import only what we need. Plots use **matplotlib** only (no seaborn). Randomness is fixed for reproducibility.\n"},{"cell_type":"code","execution_count":null,"id":"b0b8beb1","metadata":{},"outputs":[],"source":"import numpy as np\nimport matplotlib.pyplot as plt\n\nfrom sklearn import datasets\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.decomposition import PCA\nfrom sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.model_selection import StratifiedKFold, cross_val_score, train_test_split\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.metrics import accuracy_score, confusion_matrix, classification_report\n\nplt.rcParams[\"figure.figsize\"] = (6, 4)\nRANDOM_STATE = 42\nnp.random.seed(RANDOM_STATE)\n\nprint(\"PCA/LDA environment ready.\")"},{"cell_type":"markdown","id":"64e44111","metadata":{},"source":"## 1. PCA — Concept, Objective, Geometry\n\n**Goal.** Find orthogonal directions (principal components) capturing **maximum variance** of centered data \\(X \\in \\mathbb{R}^{n \\times d}\\).  \nSample covariance: \\( \\Sigma = \\tfrac{1}{n} X^\\top X \\). Solve\n\\[\n\\max_{W \\in \\mathbb{R}^{d \\times k}} \\mathrm{Tr}(W^\\top \\Sigma W) \\quad \\text{s.t. } W^\\top W = I_k,\n\\]\nyielding top-\\(k\\) eigenvectors of \\( \\Sigma \\). Equivalently, use SVD of \\(X\\) for numerical stability.\n\n**Key points**\n- **Unsupervised** (labels not used).\n- **Scale-sensitive** → standardize (zero-mean, unit-variance).\n- **Interpretability** via **loadings** (component weights).\n- High variance ≠ class separability. PCA is not discriminative by design.\n"},{"cell_type":"markdown","id":"d6f8de8c","metadata":{},"source":"### 1.1 PCA on Iris — Standardize → Project → Diagnose\n\nWe (1) standardize, (2) project to 2D with PCA, (3) visualize, (4) examine explained variance, (5) inspect loadings.\n"},{"cell_type":"code","execution_count":null,"id":"801b322d","metadata":{},"outputs":[],"source":"# Iris dataset\niris = datasets.load_iris()\nX_iris, y_iris = iris.data, iris.target\nclass_names_iris = iris.target_names\n\n# 1) Standardize\nscaler = StandardScaler()\nX_iris_std = scaler.fit_transform(X_iris)\n\n# 2) PCA to 2D\npca2 = PCA(n_components=2, random_state=RANDOM_STATE)\nX_iris_pca2 = pca2.fit_transform(X_iris_std)\n\nprint(\"Explained variance ratio (2 comps):\", np.round(pca2.explained_variance_ratio_, 4))\nprint(\"Cumulative (2 comps):\", np.round(np.sum(pca2.explained_variance_ratio_), 4))\n\n# 3) Scatter (distinct markers; no explicit colors)\nmarkers = [\"o\", \"s\", \"^\"]\nplt.figure()\nfor i, m in enumerate(markers):\n    plt.scatter(X_iris_pca2[y_iris == i, 0], X_iris_pca2[y_iris == i, 1], marker=m, label=class_names_iris[i])\nplt.xlabel(\"PC1\")\nplt.ylabel(\"PC2\")\nplt.title(\"Iris — PCA 2D Projection\")\nplt.legend()\nplt.show()\n\n# 4) Explained variance per PC (k=2)\nplt.figure()\nplt.bar([1, 2], pca2.explained_variance_ratio_)\nplt.xticks([1, 2])\nplt.ylabel(\"Explained Variance Ratio\")\nplt.title(\"Explained Variance Ratio — Iris PCA (k=2)\")\nplt.show()\n\n# 5) Loadings (components x features)\nprint(\"PCA component loadings:\\n\", np.round(pca2.components_, 4))"},{"cell_type":"markdown","id":"8fbf1e7d","metadata":{},"source":"### 1.2 Reconstruction Error vs. k & Scree/Cumulative Curves\n\n**Reconstruction error** measures information loss after projection and inverse projection.  \n**Scree plot** and **cumulative variance** guide the choice of \\(k\\).\n"},{"cell_type":"code","execution_count":null,"id":"5028ca94","metadata":{},"outputs":[],"source":"# Reconstruction MSE vs. k (Iris)\nmax_k = X_iris.shape[1]\nks = list(range(1, max_k + 1))\nmse_by_k = []\n\nfor k in ks:\n    p = PCA(n_components=k, random_state=RANDOM_STATE).fit(X_iris_std)\n    X_proj = p.transform(X_iris_std)\n    X_rec = p.inverse_transform(X_proj)\n    mse = np.mean((X_iris_std - X_rec) ** 2)\n    mse_by_k.append(mse)\n\nplt.figure()\nplt.plot(ks, mse_by_k, marker=\"o\")\nplt.xlabel(\"Number of Components (k)\")\nplt.ylabel(\"Mean Squared Reconstruction Error\")\nplt.title(\"PCA Reconstruction Error vs. k — Iris\")\nplt.show()\n\n# Scree & cumulative\npca_full = PCA().fit(X_iris_std)\nevar = pca_full.explained_variance_ratio_\ncum_evar = np.cumsum(evar)\n\nplt.figure()\nplt.plot(range(1, len(evar)+1), evar, marker=\"o\")\nplt.xlabel(\"Component index\")\nplt.ylabel(\"Explained Variance Ratio\")\nplt.title(\"Scree Plot — Iris\")\nplt.show()\n\nplt.figure()\nplt.plot(range(1, len(cum_evar)+1), cum_evar, marker=\"o\")\nplt.xlabel(\"Number of Components\")\nplt.ylabel(\"Cumulative Explained Variance\")\nplt.title(\"Cumulative Explained Variance — Iris\")\nplt.show()\n\nprint(\"Explained variance ratio:\", np.round(evar, 4))\nprint(\"Cumulative variance:\", np.round(cum_evar, 4))"},{"cell_type":"markdown","id":"d98f5d2f","metadata":{},"source":"## 2. LDA — Concept, Objective, Assumptions\n\n**Goal.** With labels, find projections that **maximize between-class scatter** while **minimizing within-class scatter**.  \nLet global mean \\( \\mu \\), class means \\( \\mu_c \\), sizes \\( n_c \\):\n\\[\nS_W = \\sum_c \\sum_{i \\in c} (x_i - \\mu_c)(x_i - \\mu_c)^\\top,\\qquad\nS_B = \\sum_c n_c (\\mu_c - \\mu)(\\mu_c - \\mu)^\\top.\n\\]\nSolve \\( S_B w = \\lambda S_W w \\). Max components \\(= \\#\\text{classes} - 1\\).\n\n**Assumptions & caveats**\n- Class-conditional normality; **shared covariance**.\n- Standardize features; watch for outliers/collinearity.\n- In high-d low-n, \\(S_W\\) can be singular → regularize or reduce \\(d\\).\n"},{"cell_type":"markdown","id":"b2d1cf81","metadata":{},"source":"### 2.1 LDA on Iris — Supervised 2D Projection\n\nWe expect tighter class separation than PCA because labels guide the projection.\n"},{"cell_type":"code","execution_count":null,"id":"88ebad93","metadata":{},"outputs":[],"source":"lda2 = LDA(n_components=2)\nX_iris_lda2 = lda2.fit_transform(X_iris_std, y_iris)\n\nplt.figure()\nfor i, m in enumerate([\"o\", \"s\", \"^\"]):\n    plt.scatter(X_iris_lda2[y_iris == i, 0], X_iris_lda2[y_iris == i, 1], marker=m, label=class_names_iris[i])\nplt.xlabel(\"LD1\")\nplt.ylabel(\"LD2\")\nplt.title(\"Iris — LDA 2D Projection\")\nplt.legend()\nplt.show()\n\nprint(\"LDA explained_variance_ratio_:\", np.round(lda2.explained_variance_ratio_, 4))"},{"cell_type":"markdown","id":"1332c5d2","metadata":{},"source":"## 3. PCA vs. LDA for Classification — Wine & Digits\n\nPipelines compared:\n1) **LR only**: Standardize → Logistic Regression  \n2) **PCA+LR**: Standardize → PCA (retain 95% variance) → LR  \n3) **LDA+LR**: Standardize → LDA (≤ classes − 1) → LR  \n\nMetric: accuracy via **Stratified 5-fold CV**.\n"},{"cell_type":"code","execution_count":null,"id":"a3ee0af7","metadata":{},"outputs":[],"source":"# Wine comparison\nwine = datasets.load_wine()\nX_wine, y_wine = wine.data, wine.target\n\ncv = StratifiedKFold(n_splits=5, shuffle=True, random_state=RANDOM_STATE)\n\npipe_lr = Pipeline([(\"scaler\", StandardScaler()),\n                    (\"clf\", LogisticRegression(max_iter=1000, random_state=RANDOM_STATE))])\n\npipe_pca_lr = Pipeline([(\"scaler\", StandardScaler()),\n                        (\"pca\", PCA(n_components=0.95, random_state=RANDOM_STATE)),\n                        (\"clf\", LogisticRegression(max_iter=1000, random_state=RANDOM_STATE))])\n\npipe_lda_lr = Pipeline([(\"scaler\", StandardScaler()),\n                        (\"lda\", LDA(n_components=min(len(np.unique(y_wine)) - 1, X_wine.shape[1]))),\n                        (\"clf\", LogisticRegression(max_iter=1000, random_state=RANDOM_STATE))])\n\nscores_lr  = cross_val_score(pipe_lr,     X_wine, y_wine, cv=cv, scoring=\"accuracy\")\nscores_pca = cross_val_score(pipe_pca_lr, X_wine, y_wine, cv=cv, scoring=\"accuracy\")\nscores_lda = cross_val_score(pipe_lda_lr, X_wine, y_wine, cv=cv, scoring=\"accuracy\")\n\nprint(\"Wine accuracy (mean ± std)\")\nprint(\"  LR only :\", np.mean(scores_lr),  \"±\", np.std(scores_lr))\nprint(\"  PCA+LR  :\", np.mean(scores_pca), \"±\", np.std(scores_pca))\nprint(\"  LDA+LR  :\", np.mean(scores_lda), \"±\", np.std(scores_lda))"},{"cell_type":"code","execution_count":null,"id":"e6a759f2","metadata":{},"outputs":[],"source":"# Digits comparison (64 features)\ndigits = datasets.load_digits()\nX_digits, y_digits = digits.data, digits.target\n\ncv = StratifiedKFold(n_splits=5, shuffle=True, random_state=RANDOM_STATE)\n\npipe_lr_d = Pipeline([(\"scaler\", StandardScaler()),\n                      (\"clf\", LogisticRegression(max_iter=2000, random_state=RANDOM_STATE))])\n\npipe_pca_lr_d = Pipeline([(\"scaler\", StandardScaler()),\n                          (\"pca\", PCA(n_components=0.95, random_state=RANDOM_STATE)),\n                          (\"clf\", LogisticRegression(max_iter=2000, random_state=RANDOM_STATE))])\n\npipe_lda_lr_d = Pipeline([(\"scaler\", StandardScaler()),\n                          (\"lda\", LDA(n_components=min(len(np.unique(y_digits)) - 1, X_digits.shape[1]))),\n                          (\"clf\", LogisticRegression(max_iter=2000, random_state=RANDOM_STATE))])\n\nscores_lr_d  = cross_val_score(pipe_lr_d,     X_digits, y_digits, cv=cv, scoring=\"accuracy\")\nscores_pca_d = cross_val_score(pipe_pca_lr_d, X_digits, y_digits, cv=cv, scoring=\"accuracy\")\nscores_lda_d = cross_val_score(pipe_lda_lr_d, X_digits, y_digits, cv=cv, scoring=\"accuracy\")\n\nprint(\"Digits accuracy (mean ± std)\")\nprint(\"  LR only :\", np.mean(scores_lr_d),  \"±\", np.std(scores_lr_d))\nprint(\"  PCA+LR  :\", np.mean(scores_pca_d), \"±\", np.std(scores_pca_d))\nprint(\"  LDA+LR  :\", np.mean(scores_lda_d), \"±\", np.std(scores_lda_d))"},{"cell_type":"markdown","id":"e1c56ceb","metadata":{},"source":"### 3.1 Fixed Train/Test Snapshot with Confusion Matrices (Wine)\n\nA single split provides concrete confusion patterns and per-class metrics.\n"},{"cell_type":"code","execution_count":null,"id":"b68f3b90","metadata":{},"outputs":[],"source":"X_tr, X_te, y_tr, y_te = train_test_split(\n    X_wine, y_wine, test_size=0.3, random_state=RANDOM_STATE, stratify=y_wine\n)\n\ndef plot_cm(cm, title):\n    plt.figure()\n    plt.imshow(cm, aspect=\"auto\")\n    plt.title(title)\n    plt.xlabel(\"Predicted\")\n    plt.ylabel(\"True\")\n    plt.colorbar()\n    plt.show()\n\n# Baseline\npipe_lr.fit(X_tr, y_tr)\ny_pred_base = pipe_lr.predict(X_te)\ncm_base = confusion_matrix(y_te, y_pred_base)\nprint(\"Accuracy (LR only):\", accuracy_score(y_te, y_pred_base))\nplot_cm(cm_base, \"Confusion Matrix — LR only\")\nprint(classification_report(y_te, y_pred_base))\n\n# PCA pipeline\npipe_pca_lr.fit(X_tr, y_tr)\ny_pred_pca = pipe_pca_lr.predict(X_te)\ncm_pca = confusion_matrix(y_te, y_pred_pca)\nprint(\"Accuracy (PCA + LR):\", accuracy_score(y_te, y_pred_pca))\nplot_cm(cm_pca, \"Confusion Matrix — PCA + LR\")\nprint(classification_report(y_te, y_pred_pca))\n\n# LDA pipeline\npipe_lda_lr.fit(X_tr, y_tr)\ny_pred_lda = pipe_lda_lr.predict(X_te)\ncm_lda = confusion_matrix(y_te, y_pred_lda)\nprint(\"Accuracy (LDA + LR):\", accuracy_score(y_te, y_pred_lda))\nplot_cm(cm_lda, \"Confusion Matrix — LDA + LR\")\nprint(classification_report(y_te, y_pred_lda))"},{"cell_type":"markdown","id":"4c519cf9","metadata":{},"source":"## 4. Practical Guidance — PCA vs. LDA\n\n**Use PCA when** you need unsupervised compression/visualization/denoising or to mitigate multicollinearity before a classifier.  \n**Use LDA when** labels are available and the primary goal is to **maximize class separability** with a linear projection.\n\n**Pitfalls**\n- PCA ignores labels; variance capture ≠ discriminative power.\n- LDA assumes shared covariance; violations hurt.\n- Always standardize; check outliers and imbalance.\n- In \\(d \\gg n\\), consider PCA→LDA or regularization.\n"}],"metadata":{"kaggle":{"accelerator":"none","dataSources":[],"dockerImageVersionId":31089,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.11.13"},"papermill":{"default_parameters":{},"duration":280.84057,"end_time":"2025-10-20T02:54:50.307475","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2025-10-20T02:50:09.466905","version":"2.6.0"}},"nbformat":4,"nbformat_minor":4}