{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":81933,"databundleVersionId":9643020,"sourceType":"competition"},{"sourceId":7453542,"sourceType":"datasetVersion","datasetId":921302}],"dockerImageVersionId":30761,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip -q install /kaggle/input/pytorchtabnet/pytorch_tabnet-4.1.0-py3-none-any.whl","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:38:05.418991Z","iopub.execute_input":"2024-12-01T15:38:05.420180Z","iopub.status.idle":"2024-12-01T15:38:50.052228Z","shell.execute_reply.started":"2024-12-01T15:38:05.420132Z","shell.execute_reply":"2024-12-01T15:38:50.051029Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport random\nimport pickle\n\nimport warnings\nfrom functools import partial\nfrom pathlib import Path\n\nimport numpy as np\nimport torch\nfrom torch import nn\nimport torch.optim as optim\nimport torch.nn.functional as F\nimport optuna\nimport pandas as pd\nimport polars as pl\nimport polars.selectors as cs\nfrom polars.testing import assert_frame_equal\nimport lightgbm as lgb\nimport xgboost as xgb\nfrom xgboost import XGBRegressor\nfrom catboost import CatBoostRegressor, MultiTargetCustomMetric\nfrom numpy.typing import ArrayLike, NDArray\nfrom sklearn.base import BaseEstimator, RegressorMixin\nfrom sklearn.metrics import cohen_kappa_score\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.model_selection import StratifiedKFold\nfrom sklearn.impute import SimpleImputer, KNNImputer\nfrom sklearn.base import clone\nfrom scipy.optimize import minimize\n\nfrom pytorch_tabnet.callbacks import Callback\nfrom pytorch_tabnet.tab_model import TabNetRegressor\nfrom pytorch_tabnet.metrics import Metric\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:38:50.054441Z","iopub.execute_input":"2024-12-01T15:38:50.054774Z","iopub.status.idle":"2024-12-01T15:38:56.630550Z","shell.execute_reply.started":"2024-12-01T15:38:50.054744Z","shell.execute_reply":"2024-12-01T15:38:56.629449Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# CFG\nclass CFG:\n    N_FOLD = 5\n    random_seed=42\n\nIS_TEST = False\n    ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:38:56.631810Z","iopub.execute_input":"2024-12-01T15:38:56.632314Z","iopub.status.idle":"2024-12-01T15:38:56.637692Z","shell.execute_reply.started":"2024-12-01T15:38:56.632281Z","shell.execute_reply":"2024-12-01T15:38:56.636387Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# set random seed\ndef seed_everything(seed: int):\n    \n    random.seed(seed)\n    os.environ['PYTHONHASHSEED'] = str(seed)\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n    torch.cuda.manual_seed(seed)\n    torch.backends.cudnn.deterministic = True\n    torch.backends.cudnn.benchmark = True\n    \nseed_everything(seed=CFG.random_seed)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:38:56.638887Z","iopub.execute_input":"2024-12-01T15:38:56.639224Z","iopub.status.idle":"2024-12-01T15:38:56.654861Z","shell.execute_reply.started":"2024-12-01T15:38:56.639175Z","shell.execute_reply":"2024-12-01T15:38:56.653514Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Feature\nFEATURE_COLS = [\n 'Basic_Demos-Age',\n 'CGAS-CGAS_Score',\n 'Physical-BMI',\n 'Physical-Height',\n 'Physical-Weight',\n 'Physical-Waist_Circumference',\n 'Physical-Diastolic_BP',\n 'Physical-HeartRate',\n 'Physical-Systolic_BP',\n 'Fitness_Endurance-Max_Stage',\n 'Fitness_Endurance-Time_Mins',\n 'Fitness_Endurance-Time_Sec',\n 'FGC-FGC_CU',\n 'FGC-FGC_CU_Zone',\n 'FGC-FGC_GSND',\n 'FGC-FGC_GSND_Zone',\n 'FGC-FGC_GSD',\n 'FGC-FGC_GSD_Zone',\n 'FGC-FGC_PU',\n 'FGC-FGC_PU_Zone',\n 'FGC-FGC_SRL',\n 'FGC-FGC_SRL_Zone',\n 'FGC-FGC_SRR',\n 'FGC-FGC_SRR_Zone',\n 'FGC-FGC_TL',\n 'FGC-FGC_TL_Zone',\n 'BIA-BIA_Activity_Level_num',\n 'BIA-BIA_BMC',\n 'BIA-BIA_BMI',\n 'BIA-BIA_BMR',\n 'BIA-BIA_DEE',\n 'BIA-BIA_ECW',\n 'BIA-BIA_FFM',\n 'BIA-BIA_FFMI',\n 'BIA-BIA_FMI',\n 'BIA-BIA_Fat',\n 'BIA-BIA_Frame_num',\n 'BIA-BIA_ICW',\n 'BIA-BIA_LDM',\n 'BIA-BIA_LST',\n 'BIA-BIA_SMM',\n 'BIA-BIA_TBW',\n 'PAQ_A-PAQ_A_Total',\n 'PAQ_C-PAQ_C_Total',\n 'SDS-SDS_Total_Raw',\n 'SDS-SDS_Total_T',\n 'PreInt_EduHx-computerinternet_hoursday',\n]\n\nTIME_SERIES_COLS = [\n    'non_wear_percentage',\n     'unique_days',\n     'count_enmo',\n     'count_anglez',\n     'count_light',\n     'count_battery_voltage',\n     'null_count_enmo',\n     'null_count_anglez',\n     'null_count_light',\n     'null_count_battery_voltage',\n     'mean_enmo',\n     'mean_anglez',\n     'mean_light',\n     'mean_battery_voltage',\n     'std_enmo',\n     'std_anglez',\n     'std_light',\n     'std_battery_voltage',\n     'min_enmo',\n     'min_anglez',\n     'min_light',\n     'min_battery_voltage',\n     '25%_enmo',\n     '25%_anglez',\n     '25%_light',\n     '25%_battery_voltage',\n     '50%_enmo',\n     '50%_anglez',\n     '50%_light',\n     '50%_battery_voltage',\n     '75%_enmo',\n     '75%_anglez',\n     '75%_light',\n     '75%_battery_voltage',\n     'max_enmo',\n     'max_anglez',\n     'max_light',\n     'max_battery_voltage',\n     'no_motion_total_duration_median',\n     'no_motion_total_duration_max',\n     'no_motion_total_duration_std',\n     'no_motion_count_periods_median',\n     'no_motion_count_periods_max',\n     'no_motion_count_periods_std',\n     'std_across_hours_median_day',\n     'std_across_hours_max_day',\n     'std_across_hours_std_day',\n     'peak_hour_median_day',\n     'peak_hour_max_day',\n     'peak_hour_std_day',\n     'entropy_median_day',\n     'entropy_max_day',\n     'entropy_std_day',\n     'circadian_rhythm_std_across_hours_median_night',\n     'circadian_rhythm_std_across_hours_max_night',\n     'circadian_rhythm_std_across_hours_std_night',\n     'circadian_rhythm_peak_hour_median_night',\n     'circadian_rhythm_peak_hour_max_night',\n     'circadian_rhythm_peak_hour_std_night',\n     'circadian_rhythm_entropy_median_night',\n     'circadian_rhythm_entropy_max_night',\n     'circadian_rhythm_entropy_std_night',\n     'physical_activity_total_duration_median',\n     'physical_activity_total_duration_max',\n     'physical_activity_total_duration_std',\n     'physical_activity_count_periods_median',\n     'physical_activity_count_periods_max',\n     'physical_activity_count_periods_std'\n]\n\n# Target\nPCIAT_COLS = [\n    f\"PCIAT-PCIAT_{i+1:02}\" for i in range(20)\n]\nTOTAL_COLS = [\n    \"PCIAT-PCIAT_Total\",    \n]\nSII_COLS = [\n    'sii',\n]\nTARGET_COLS = (\n    PCIAT_COLS + \n    # TOTAL_COLS + \n    SII_COLS\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:38:56.659261Z","iopub.execute_input":"2024-12-01T15:38:56.659693Z","iopub.status.idle":"2024-12-01T15:38:56.669348Z","shell.execute_reply.started":"2024-12-01T15:38:56.659663Z","shell.execute_reply":"2024-12-01T15:38:56.668312Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"DATA_DIR = Path(\"/kaggle/input/child-mind-institute-problematic-internet-use\")\n\ntrain = pl.read_csv(DATA_DIR / 'train.csv')\ntest = pl.read_csv(DATA_DIR / 'test.csv')\ndata_dict = pl.read_csv(DATA_DIR / 'data_dictionary.csv')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:38:56.670759Z","iopub.execute_input":"2024-12-01T15:38:56.671169Z","iopub.status.idle":"2024-12-01T15:38:56.792918Z","shell.execute_reply.started":"2024-12-01T15:38:56.671139Z","shell.execute_reply":"2024-12-01T15:38:56.792040Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Feature Engineering","metadata":{}},{"cell_type":"markdown","source":"## Table Features\n","metadata":{}},{"cell_type":"code","source":"def feature_engineering(df: pl.DataFrame) -> pl.DataFrame:\n    df = df.with_columns(\n        (pl.col(\"Physical-BMI\") * pl.col(\"Basic_Demos-Age\")).alias(\"BMI_Age\"),\n        (pl.col(\"PreInt_EduHx-computerinternet_hoursday\") * pl.col(\"Basic_Demos-Age\")).alias(\"Internet_Hours_Age\"),\n        (pl.col(\"Physical-BMI\") * pl.col(\"PreInt_EduHx-computerinternet_hoursday\")).alias(\"BMI_Internet_Hours\"),\n        (pl.col(\"BIA-BIA_Fat\") / pl.col(\"BIA-BIA_BMI\")).alias(\"BFP_BMI\"),\n        (pl.col(\"BIA-BIA_FFMI\") / pl.col(\"BIA-BIA_Fat\")).alias(\"FFMI_BFP\"),\n        (pl.col(\"BIA-BIA_FMI\") / pl.col(\"BIA-BIA_Fat\")).alias(\"FMI_BFP\"),\n        (pl.col(\"BIA-BIA_LST\") / pl.col(\"BIA-BIA_TBW\")).alias(\"LST_TBW\"),\n        (pl.col(\"BIA-BIA_Fat\") * pl.col(\"BIA-BIA_BMR\")).alias(\"BFP_BMR\"),\n        (pl.col(\"BIA-BIA_Fat\") * pl.col(\"BIA-BIA_DEE\")).alias(\"BFP_DEE\"),\n        (pl.col(\"BIA-BIA_DEE\") / pl.col(\"Physical-Weight\")).alias(\"DEE_Weight\"),\n        (pl.col(\"BIA-BIA_SMM\") / pl.col(\"Physical-Height\")).alias(\"SMM_Height\"),\n        (pl.col(\"BIA-BIA_SMM\") / pl.col(\"BIA-BIA_FMI\")).alias(\"Muscle_to_Fat\"),\n        (pl.col(\"BIA-BIA_TBW\") / pl.col(\"Physical-Weight\")).alias(\"Hydration_Status\"),\n        (pl.col(\"BIA-BIA_ICW\") / pl.col(\"BIA-BIA_TBW\")).alias(\"ICW_TBW\"),\n    )\n    \n    return df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:38:56.794300Z","iopub.execute_input":"2024-12-01T15:38:56.794731Z","iopub.status.idle":"2024-12-01T15:38:56.804230Z","shell.execute_reply.started":"2024-12-01T15:38:56.794686Z","shell.execute_reply":"2024-12-01T15:38:56.803149Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train = feature_engineering(train)\ntest = feature_engineering(test)\n\ntrain.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:38:56.805599Z","iopub.execute_input":"2024-12-01T15:38:56.806035Z","iopub.status.idle":"2024-12-01T15:38:56.905982Z","shell.execute_reply.started":"2024-12-01T15:38:56.805986Z","shell.execute_reply":"2024-12-01T15:38:56.905080Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Time Series Features","metadata":{}},{"cell_type":"code","source":"def describe_feature(file_path, participant_id):\n    \"\"\" describeで計算できる統計量の特徴量 \"\"\"\n\n    data = pl.read_parquet(file_path)\n    non_wear_percentage = (data['non-wear_flag'].sum() / len(data)) * 100\n    worn_data = data.filter(\n        pl.col('non-wear_flag') == 0\n    )\n\n    describe_cols = [\n        \"enmo\",\n        \"anglez\",\n        \"light\",\n        \"battery_voltage\",\n    ]\n\n    described_worn_data = worn_data.select(describe_cols).describe()\n    described_worn_data_dicts = described_worn_data.to_dicts()\n\n    result_dict = {\n        \"id\": participant_id.replace('id=', ''), \n        \"non_wear_percentage\": non_wear_percentage,    \n    }\n\n    # nuniqueに関する処理\n    result_dict[\"unique_days\"] = worn_data['relative_date_PCIAT'].n_unique()\n\n    # describeに関する処理\n    for described_worn_data_dict in described_worn_data_dicts:\n        statistic = described_worn_data_dict[\"statistic\"]\n        \n        d = {f\"{statistic}_{i}\" : v for i, v in described_worn_data_dict.items() if i != \"statistic\"}\n        result_dict.update(d)       \n\n    return result_dict","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:38:56.907332Z","iopub.execute_input":"2024-12-01T15:38:56.907693Z","iopub.status.idle":"2024-12-01T15:38:56.915419Z","shell.execute_reply.started":"2024-12-01T15:38:56.907653Z","shell.execute_reply":"2024-12-01T15:38:56.914407Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def no_motion_feature(file_path, participant_id):\n    \"\"\" no_motion 特徴量 \"\"\"\n    data = pl.read_parquet(file_path)\n    worn_data = data.filter(\n        pl.col('non-wear_flag') == 0\n    )    \n    \n    # recalculate time difference between rows and measurement_after_gap flag\n    worn_data = worn_data.with_columns(\n        (pl.col(\"time_of_day\")  / 1e9 / 3600).alias(\"time_of_day_hours\"), #nanoseconds to hours,\n    ).with_columns(\n        (pl.col(\"relative_date_PCIAT\") + (pl.col('time_of_day_hours') / 24)).alias(\"day_time\")\n    )\n\n    worn_data = worn_data.with_columns(\n        (pl.col(\"day_time\").diff() * 86400).round(0).alias(\"time_diff\")\n    ).with_columns(\n        (pl.col(\"time_diff\") > 5).alias(\"measurement_after_gap\") #expected_diff\n    )\n\n    worn_data = worn_data.with_columns(\n        no_motion = pl.col('enmo') == 0\n    ).with_columns([\n        ((pl.col(\"no_motion\") != pl.col(\"no_motion\").shift(1)) |\n         pl.col(\"measurement_after_gap\")).cum_sum().alias(\"motion_group\")\n    ])\n\n    no_motion_periods = worn_data.filter(\"no_motion\").group_by(\n        \"motion_group\"\n    ).agg(\n        pl.col(\"day_time\").min().alias(\"no_motion_day_time_min\"), \n        pl.col(\"day_time\").max().alias(\"no_motion_day_time_max\")\n    )\n\n    no_motion_periods = no_motion_periods.with_columns(\n        (\n            (\n                pl.col('no_motion_day_time_max') - pl.col('no_motion_day_time_min')\n            ) * 86400\n        ).round(0).cast(int).alias(\"duration_sec\")\n    )\n\n    no_motion_periods = no_motion_periods.with_columns(\n        pl.col('no_motion_day_time_min').cast(int).alias(\"day\")\n    )\n    daily_stats = no_motion_periods.group_by(\"day\").agg(\n        pl.col(\"duration_sec\").sum().alias(\"total_duration\"), \n        pl.col(\"duration_sec\").count().alias(\"count_periods\")\n    )\n\n    features = daily_stats.select(\n        pl.col(\"total_duration\").median().alias(\"no_motion_total_duration_median\"), \n        pl.col(\"total_duration\").max().alias(\"no_motion_total_duration_max\"),\n        pl.col(\"total_duration\").std().alias(\"no_motion_total_duration_std\"),\n        pl.col(\"count_periods\").median().alias(\"no_motion_count_periods_median\"), \n        pl.col(\"count_periods\").max().alias(\"no_motion_count_periods_max\"),\n        pl.col(\"count_periods\").std().alias(\"no_motion_count_periods_std\"),\n    ).with_columns(\n        pl.lit(participant_id.replace('id=', '')).alias(\"id\")\n    )\n    \n    return features","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:38:56.917028Z","iopub.execute_input":"2024-12-01T15:38:56.917796Z","iopub.status.idle":"2024-12-01T15:38:56.932281Z","shell.execute_reply.started":"2024-12-01T15:38:56.917751Z","shell.execute_reply":"2024-12-01T15:38:56.931120Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def circadian_rhythm_feature(file_path, participant_id):\n    \"\"\" circadian_rhythm_feature \"\"\"\n\n    data = pl.read_parquet(file_path)\n    worn_data = data.filter(\n        pl.col('non-wear_flag') == 0\n    ) \n\n    day_start_hour = 8\n    day_end_hour = 21\n    # recalculate time difference between rows and measurement_after_gap flag    \n    worn_data = worn_data.with_columns(\n        (pl.col(\"time_of_day\")  / 1e9 / 3600).cast(int).alias(\"time_of_day_hours\"),\n        pl.col(\"relative_date_PCIAT\").cast(int).alias(\"relative_date_PCIAT\"),        \n    ).with_columns(\n        pl.when(\n            (pl.col(\"time_of_day_hours\") >= day_start_hour)\n            & (pl.col(\"time_of_day_hours\") < day_end_hour)\n        ).then(pl.lit(\"day\"))\n        .otherwise(pl.lit(\"night\")).alias(\"day_period\")\n    )\n\n    hourly_activity = worn_data.group_by(\n        \"relative_date_PCIAT\",\n        \"time_of_day_hours\",\n        \"day_period\",\n    ).agg(\n        pl.col(\"enmo\").mean().alias(\"enmo_mean\"), \n        pl.col(\"enmo\").max().alias(\"enmo_max\"),\n    )\n    \n    features = hourly_activity.sort(\n        \"enmo_mean\", descending=True\n    ).group_by(['relative_date_PCIAT', 'day_period']).agg(\n        pl.col(\"enmo_mean\").std().alias(\"std_across_hours\"),\n        pl.col(\"time_of_day_hours\").first().alias(\"peak_hour\"),\n        (\n            -(pl.col(\"enmo_mean\") / pl.col(\"enmo_mean\").sum() * np.log(pl.col(\"enmo_mean\") / pl.col(\"enmo_mean\").sum() + 1e-9)).sum()\n        ).alias(\"entropy\")\n    )\n\n    day_features = features.filter(pl.col(\"day_period\") == \"day\").select(\n        pl.col(\"std_across_hours\").median().alias(\"std_across_hours_median\"), \n        pl.col(\"std_across_hours\").max().alias(\"std_across_hours_max\"),\n        pl.col(\"std_across_hours\").std().alias(\"std_across_hours_std\"),\n        pl.col(\"peak_hour\").median().alias(\"peak_hour_median\"), \n        pl.col(\"peak_hour\").max().alias(\"peak_hour_max\"),\n        pl.col(\"peak_hour\").std().alias(\"peak_hour_std\"),\n        pl.col(\"entropy\").median().alias(\"entropy_median\"), \n        pl.col(\"entropy\").max().alias(\"entropy_max\"),\n        pl.col(\"entropy\").std().alias(\"entropy_std\"),        \n    ).with_columns(\n        pl.all().name.suffix(\"_day\")\n    ).select(\n        pl.selectors.ends_with(\"_day\")\n    )\n    \n    night_features = features.filter(pl.col(\"day_period\") == \"night\").select(\n        pl.col(\"std_across_hours\").median().alias(\"circadian_rhythm_std_across_hours_median\"),\n        pl.col(\"std_across_hours\").max().alias(\"circadian_rhythm_std_across_hours_max\"),\n        pl.col(\"std_across_hours\").std().alias(\"circadian_rhythm_std_across_hours_std\"),\n        pl.col(\"peak_hour\").median().alias(\"circadian_rhythm_peak_hour_median\"),\n        pl.col(\"peak_hour\").max().alias(\"circadian_rhythm_peak_hour_max\"),\n        pl.col(\"peak_hour\").std().alias(\"circadian_rhythm_peak_hour_std\"),\n        pl.col(\"entropy\").median().alias(\"circadian_rhythm_entropy_median\"),\n        pl.col(\"entropy\").max().alias(\"circadian_rhythm_entropy_max\"),\n        pl.col(\"entropy\").std().alias(\"circadian_rhythm_entropy_std\"),\n    ).with_columns(\n        pl.all().name.suffix(\"_night\")\n    ).select(\n        pl.selectors.ends_with(\"_night\")\n    )\n\n    feature_df = pl.concat(\n        [day_features, night_features], how=\"horizontal\"\n    ).with_columns(\n        pl.lit(participant_id.replace('id=', '')).alias(\"id\")\n    )\n    \n    return feature_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:38:56.933742Z","iopub.execute_input":"2024-12-01T15:38:56.934523Z","iopub.status.idle":"2024-12-01T15:38:56.951302Z","shell.execute_reply.started":"2024-12-01T15:38:56.934476Z","shell.execute_reply":"2024-12-01T15:38:56.950257Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def merge_mvpa_groups(df, allowed_gap=60, merge_gap=60):\n    # 1. 最後の is_mvpa の時刻を ffill と shift で取得\n    df = df.with_columns(\n        pl.when(pl.col(\"is_mvpa\")).then(pl.col(\"day_time\"))\n        .forward_fill()\n        .shift(1).alias(\"last_mvpa_time\")    \n    )\n\n    # 2. is_mvpa が True の場合の時間差を秒単位に変換\n    df = df.with_columns(\n        (\n            (pl.col(\"day_time\") - pl.col(\"last_mvpa_time\")) * 86400\n        ).round(0).alias(\"mvpa_time_diff\")\n    )\n\n    # 3. グループを作成するための条件式\n    df = df.with_columns(\n        (\n            (pl.col(\"is_mvpa\") != pl.col(\"is_mvpa\").shift(1)) |\n            (pl.col(\"mvpa_time_diff\") >= allowed_gap)\n        ).cum_sum().alias(\"mvpa_group\")\n    ).with_columns(\n        ((pl.col(\"mvpa_group\") != pl.col(\"mvpa_group\").shift(1)) & pl.col(\"is_mvpa\")).alias(\"is_mvpa_start\")\n    )\n\n\n    # 5. マージするグループのインクリメント条件\n    df = df.with_columns(\n        (pl.col(\"is_mvpa_start\") & (\n            (pl.col(\"mvpa_time_diff\") >= merge_gap) | pl.col(\"last_mvpa_time\").is_null()\n        )).alias(\"group_increment\")\n    )\n\n    # 6. 累積和で最終的なマージグループを決定し、is_mvpa が False の部分には NaN を設定\n    df = df.with_columns(\n        pl.col(\"group_increment\").cum_sum().alias(\"merged_group\")\n    )\n    df = df.with_columns(\n        (pl.when(pl.col(\"is_mvpa\"))\n        .then(pl.col(\"merged_group\"))\n        .otherwise(None)).alias(\"mvpa_merged_group\")\n    )\n\n    # 新しい DataFrame に列を追加して返す\n    return df\n\ndef physical_activity_feature(file_path, participant_id):\n    \"\"\" circadian_rhythm_feature \"\"\"\n\n    data = pl.read_parquet(file_path)\n    worn_data = data.filter(\n        pl.col('non-wear_flag') == 0\n    ) \n\n    # recalculate time difference between rows and measurement_after_gap flag\n    worn_data = worn_data.with_columns(\n        (pl.col(\"time_of_day\")  / 1e9 / 3600).alias(\"time_of_day_hours\"), #nanoseconds to hours,\n    ).with_columns(\n        (pl.col(\"relative_date_PCIAT\") + (pl.col('time_of_day_hours') / 24)).alias(\"day_time\")\n    )\n\n    worn_data = worn_data.with_columns(\n        (pl.col(\"day_time\").diff() * 86400).round(0).alias(\"time_diff\")\n    ).with_columns(\n        (pl.col(\"time_diff\") > 5).alias(\"measurement_after_gap\") #expected_diff\n    )\n    \n    window_size = 12 \n\n    worn_data = worn_data.with_columns(\n        pl.col(\"measurement_after_gap\").cum_sum().alias(\"motion_group\")        \n    )\n    \n    # Polarsでのグループごとのローリング平均の計算\n    worn_data = worn_data.with_columns(\n        pl.col(\"enmo\")\n        .rolling_mean(window_size)  # ローリング平均\n        .over(\"motion_group\")      # グループごとに適用\n        .alias(\"smoothed_enmo\")\n    )\n        \n    mvpa_threshold = 0.1\n    merge_gap = 60\n\n    worn_data = worn_data.with_columns(\n        (pl.col('enmo') > mvpa_threshold).alias(\"is_mvpa\")\n    )\n    worn_data = merge_mvpa_groups(worn_data)    \n\n    mvpa_periods = worn_data.filter(\n        pl.col(\"is_mvpa\")\n    ).group_by(\"mvpa_merged_group\").agg(\n        pl.col(\"day_time\").min().alias(\"day_time_min\"),\n        pl.col(\"day_time\").max().alias(\"day_time_max\"),\n    ).with_columns(\n        (\n            (pl.col(\"day_time_max\") - pl.col(\"day_time_min\")) * 86400\n        ).alias(\"duration_sec\")\n    ).filter(\n        pl.col(\"duration_sec\") >= 60\n    ).with_columns(\n        (pl.col(\"duration_sec\") / 60).alias(\"duration_min\")\n    ).sort(by='duration_sec')\n\n\n    mvpa_periods = mvpa_periods.with_columns(\n        pl.col(\"duration_min\").cast(int).alias('day')\n    )\n\n    daily_stats = mvpa_periods.group_by('day').agg(\n        pl.col(\"duration_sec\").sum().alias(\"total_duration\"),\n        pl.col(\"duration_sec\").count().alias(\"count_periods\"),\n    )\n\n    feature_df = daily_stats.select(\n        pl.col(\"total_duration\").median().alias(\"physical_activity_total_duration_median\"), \n        pl.col(\"total_duration\").max().alias(\"physical_activity_total_duration_max\"),\n        pl.col(\"total_duration\").std().alias(\"physical_activity_total_duration_std\"),\n        pl.col(\"count_periods\").median().alias(\"physical_activity_count_periods_median\"),\n        pl.col(\"count_periods\").max().alias(\"physical_activity_count_periods_max\"),\n        pl.col(\"count_periods\").std().alias(\"physical_activity_count_periods_std\"),\n\n    ).with_columns(\n        pl.lit(participant_id.replace('id=', '')).alias(\"id\")\n    )\n    \n    return feature_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:38:56.952913Z","iopub.execute_input":"2024-12-01T15:38:56.953353Z","iopub.status.idle":"2024-12-01T15:38:56.971906Z","shell.execute_reply.started":"2024-12-01T15:38:56.953309Z","shell.execute_reply":"2024-12-01T15:38:56.970711Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_dir = '/kaggle/input/child-mind-institute-problematic-internet-use/series_train.parquet'\ntest_dir = '/kaggle/input/child-mind-institute-problematic-internet-use/series_test.parquet'","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:38:56.973301Z","iopub.execute_input":"2024-12-01T15:38:56.973720Z","iopub.status.idle":"2024-12-01T15:38:56.988712Z","shell.execute_reply.started":"2024-12-01T15:38:56.973686Z","shell.execute_reply":"2024-12-01T15:38:56.987433Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# train\n\nresults = []\nfor participant_id in os.listdir(train_dir):\n    file_path = os.path.join(train_dir, participant_id, 'part-0.parquet')\n    result = describe_feature(file_path, participant_id)\n    results.append(result)\n\ndescribe_feature_df = pl.DataFrame(results)\n\nresults = []\nfor participant_id in os.listdir(train_dir):\n    file_path = os.path.join(train_dir, participant_id, 'part-0.parquet')\n    result = no_motion_feature(file_path, participant_id)\n    results.append(result)\n\nno_motion_feature_df = pl.concat(results)\n\nresults = []\nfor participant_id in os.listdir(train_dir):\n    file_path = os.path.join(train_dir, participant_id, 'part-0.parquet')\n    result = circadian_rhythm_feature(file_path, participant_id)\n    results.append(result)\n\ncircadian_rhythm_feature_df = pl.concat(results)\n\n\nresults = []\nfor participant_id in os.listdir(train_dir):\n    file_path = os.path.join(train_dir, participant_id, 'part-0.parquet')\n    result = physical_activity_feature(file_path, participant_id)\n    results.append(result)\n\nphysical_activity_feature_df = pl.concat(results)\n\ntrain_ts = train.select(\"id\").join(\n    describe_feature_df, on=\"id\", how=\"inner\"\n).join(\n    no_motion_feature_df, on=\"id\", how=\"inner\"\n).join(\n    circadian_rhythm_feature_df, on=\"id\", how=\"inner\"\n).join(\n    physical_activity_feature_df, on=\"id\", how=\"inner\"\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:38:56.993421Z","iopub.execute_input":"2024-12-01T15:38:56.994143Z","iopub.status.idle":"2024-12-01T15:42:13.029270Z","shell.execute_reply.started":"2024-12-01T15:38:56.994092Z","shell.execute_reply":"2024-12-01T15:42:13.027124Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# test\n\nresults = []\nfor participant_id in os.listdir(test_dir):\n    file_path = os.path.join(test_dir, participant_id, 'part-0.parquet')\n    result = describe_feature(file_path, participant_id)\n    results.append(result)\n\ndescribe_feature_df = pl.DataFrame(results)\n\nresults = []\nfor participant_id in os.listdir(test_dir):\n    file_path = os.path.join(test_dir, participant_id, 'part-0.parquet')\n    result = no_motion_feature(file_path, participant_id)\n    results.append(result)\n\nno_motion_feature_df = pl.concat(results)\n\nresults = []\nfor participant_id in os.listdir(test_dir):\n    file_path = os.path.join(test_dir, participant_id, 'part-0.parquet')\n    result = circadian_rhythm_feature(file_path, participant_id)\n    results.append(result)\n\ncircadian_rhythm_feature_df = pl.concat(results)\n\n\nresults = []\nfor participant_id in os.listdir(test_dir):\n    file_path = os.path.join(test_dir, participant_id, 'part-0.parquet')\n    result = physical_activity_feature(file_path, participant_id)\n    results.append(result)\n\nphysical_activity_feature_df = pl.concat(results)\n\ntest_ts = test.select(\"id\").join(\n    describe_feature_df, on=\"id\", how=\"inner\"\n).join(\n    no_motion_feature_df, on=\"id\", how=\"inner\"\n).join(\n    circadian_rhythm_feature_df, on=\"id\", how=\"inner\"\n).join(\n    physical_activity_feature_df, on=\"id\", how=\"inner\"\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:13.031788Z","iopub.execute_input":"2024-12-01T15:42:13.032210Z","iopub.status.idle":"2024-12-01T15:42:13.329507Z","shell.execute_reply.started":"2024-12-01T15:42:13.032171Z","shell.execute_reply":"2024-12-01T15:42:13.328321Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Attach Fold","metadata":{}},{"cell_type":"code","source":"# Fill unknown sii as \"-1\", \n# because null is not supported by stratifiedkfold function.\n\ntrain = train.with_columns(\n    pl.col(\"sii\").fill_null(-1).alias(\"sii\")\n)\n\n# Attach fold number to train df\n\nskf = StratifiedKFold(\n    n_splits = CFG.N_FOLD, \n    random_state = CFG.random_seed,\n    shuffle=True,\n)\n\ndfs = []\nval_indices = []\nfor fold, (_, test_index) in enumerate(skf.split(train, train[\"sii\"])):\n    _train = train[test_index].with_columns(\n        pl.lit(fold).alias(\"fold\")\n    )\n    dfs.append(_train)\n    val_indices.append(test_index)\n\ntrain = pl.concat(dfs)\ntrain","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:13.331380Z","iopub.execute_input":"2024-12-01T15:42:13.331840Z","iopub.status.idle":"2024-12-01T15:42:13.380728Z","shell.execute_reply.started":"2024-12-01T15:42:13.331790Z","shell.execute_reply":"2024-12-01T15:42:13.379485Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Join Features","metadata":{}},{"cell_type":"code","source":"train = train.join(train_ts, how=\"left\", on='id')\ntest = test.join(test_ts, how=\"left\", on='id')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:13.382298Z","iopub.execute_input":"2024-12-01T15:42:13.382630Z","iopub.status.idle":"2024-12-01T15:42:13.394524Z","shell.execute_reply.started":"2024-12-01T15:42:13.382600Z","shell.execute_reply":"2024-12-01T15:42:13.393536Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 列の共通部分を見つけて、欠損している箇所はnullで埋めるconcat\ntrain_test = pl.concat([train, test], how=\"diagonal\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:13.395943Z","iopub.execute_input":"2024-12-01T15:42:13.396379Z","iopub.status.idle":"2024-12-01T15:42:13.403852Z","shell.execute_reply.started":"2024-12-01T15:42:13.396335Z","shell.execute_reply":"2024-12-01T15:42:13.402666Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Cast string columns to categorical\ntrain_test = train_test.with_columns(cs.string().cast(pl.Categorical).fill_null(\"NAN\"))\ntrain = train_test[: train.height]\ntest = train_test[train.height :]\n\n# ignore rows with null values in TARGET_COLS\ntrain_without_null = train_test.drop_nulls(subset=TARGET_COLS)\n\n# X = train_without_null.select(FEATURE_COLS + TIME_SERIES_COLS)\n# X_test = test.select(FEATURE_COLS + TIME_SERIES_COLS)\n# y = train_without_null.select(PCIAT_COLS)\n\n# y_sii = train_without_null.get_column(\"sii\").to_numpy()  # ground truth\n\ncat_features = train_without_null.select(\n    FEATURE_COLS + TIME_SERIES_COLS\n).select(\n    cs.categorical()\n).columns","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:13.405379Z","iopub.execute_input":"2024-12-01T15:42:13.405844Z","iopub.status.idle":"2024-12-01T15:42:13.437139Z","shell.execute_reply.started":"2024-12-01T15:42:13.405796Z","shell.execute_reply":"2024-12-01T15:42:13.436116Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Model define and train\n\n","metadata":{}},{"cell_type":"markdown","source":"## Functions\n","metadata":{}},{"cell_type":"code","source":"def convert_pciat_pciat_k_to_sii(output: np.ndarray):\n    \"\"\" 20個のpciatのカラムの値からsiiのラベルを計算する \"\"\"\n    column_sum = np.sum(output[:, :20], axis=1)\n    boundaries = [30+10e-8, 50, 80]  # 境界値\n\n    sii = np.digitize(column_sum, boundaries)\n    return sii\n\ndef convert_pciat_pciat_total_to_sii(output: np.ndarray):\n    boundaries = [30+10e-8, 50, 80]  # 境界値\n    sii = np.digitize(output, boundaries)\n    return sii\n\ndef convert_value_to_sii(output: np.ndarray):\n    boundaries = [0.5, 1.5, 2.5]  # 境界値\n    sii = np.digitize(output, boundaries)\n    return sii\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:13.438621Z","iopub.execute_input":"2024-12-01T15:42:13.439157Z","iopub.status.idle":"2024-12-01T15:42:13.447144Z","shell.execute_reply.started":"2024-12-01T15:42:13.439103Z","shell.execute_reply":"2024-12-01T15:42:13.445974Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## XGBoost\n","metadata":{}},{"cell_type":"code","source":"# setting parameters\nXGB_Params = {\n    'learning_rate': 0.05,\n    'max_depth': 7,\n    'n_estimators': 1 if IS_TEST else 10000,\n    'subsample': 0.8,\n    'colsample_bytree': 0.8,\n    'reg_alpha': 1,  # Increased from 0.1\n    'reg_lambda': 5,  # Increased from 1\n    'random_state': CFG.random_seed,\n}\n\nclass XGBoostSiiPredicter:\n    def quadratic_kappa_coefficient(self, output, target):\n        # Calculate Quadratic Weighted Kappa\n        return cohen_kappa_score(output, target, weights='quadratic')\n    \n    def qwk(self, dtrain: xgb.DMatrix, predt: np.ndarray) -> tuple[str, float]:\n        output = pred_sii = np.clip(predt[:, -1], 0, 3).round().astype(int)\n        convert_pciat_pciat_k_to_sii(predt)\n        \n        target = dtrain.reshape(predt.shape)[:, -1]\n        # target = convert_pciat_pciat_k_to_sii(targets)\n        # Calculate and return the QWK score\n        return -self.quadratic_kappa_coefficient(output, target)\n\n    # def sii_rmse(self, dtrain: xgb.DMatrix, predt: np.ndarray) -> tuple[str, float]:\n    #     pred_sii = predt[:, -1]\n    #     target = dtrain.reshape(predt.shape)[:, -1]\n        \n    #     # Calculate and return the QWK score\n    #     return mean_squared_error(target, pred_sii, squared=True)\n\n    \n    def train(self, train_df: pl.DataFrame):\n\n        preds = []\n\n        for fold in range(CFG.N_FOLD):\n            print(f\"# XGBoost Fold {fold}\")\n            train = train_df.filter(\n                pl.col(\"fold\") != fold\n            )\n            val = train_df.filter(\n                pl.col(\"fold\") == fold\n            )\n            \n            X_train = train[FEATURE_COLS + TIME_SERIES_COLS].to_pandas()\n            X_val   = val[FEATURE_COLS + TIME_SERIES_COLS].to_pandas()\n            X_train[cat_features] = X_train[cat_features].astype(\"category\")\n            X_val[cat_features]   = X_val[cat_features].astype(\"category\")\n            \n            y_train = train[TARGET_COLS].to_pandas()\n            y_val   = val[TARGET_COLS].to_pandas()\n            \n            Xy = xgb.DMatrix(X_train, y_train)\n            evalXy = xgb.DMatrix(X_val, y_val)\n            \n            # train model\n            model = XGBRegressor(\n                **XGB_Params,\n                enable_categorical=True,\n                tree_method=\"hist\",\n                multi_strategy=\"multi_output_tree\",\n                objective=\"reg:squarederror\",\n                eval_metric=self.qwk,\n                early_stopping_rounds=100,\n            )\n            model.fit(\n                X=X_train, \n                y=y_train, \n                eval_set=[(X_val, y_val)],\n                verbose=10\n            )\n            \n            file_name = f\"./xgboost_fold_{fold}.json\"\n            model.save_model(file_name)\n            \n            pciat_pred = model.predict(X_val)\n            df = pl.DataFrame(pciat_pred, schema=[name+\"_xgb\" for name in TARGET_COLS]).with_columns(\n                (val[\"id\"]).alias(\"id\"),\n                (val[\"sii\"]).alias(\"sii\"),\n            )\n            preds.append(df)\n\n        y_pred = pl.concat(preds)\n\n        # print cv score\n        y_pred_sii = convert_value_to_sii(y_pred[\"sii_xgb\"].to_numpy())\n        y_sii = y_pred[\"sii\"]\n\n        qwk = cohen_kappa_score(y_sii, y_pred_sii, weights=\"quadratic\")\n        print(f\"Cross-Validated QWK Score: {qwk}\")\n\n        return y_pred\n    \n    def predict(self, test_df: pl.DataFrame):\n        _test_df = test_df[FEATURE_COLS + TIME_SERIES_COLS].to_pandas()\n        _test_df[cat_features] = _test_df[cat_features].astype(\"category\")\n        \n        X = xgb.DMatrix(_test_df)\n        preds: list[NDArray[np.int_]] = []\n\n        for fold in range(CFG.N_FOLD):\n            model = xgb.Booster()\n            model.load_model(f\"xgboost_fold_{fold}.json\")\n            \n            pred = model.predict(X)\n            preds.append(pred)\n\n        return np.mean(preds, axis=0)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:13.448521Z","iopub.execute_input":"2024-12-01T15:42:13.448868Z","iopub.status.idle":"2024-12-01T15:42:13.467166Z","shell.execute_reply.started":"2024-12-01T15:42:13.448838Z","shell.execute_reply":"2024-12-01T15:42:13.466274Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"xgb_regressor = XGBoostSiiPredicter()\nxgb_val_df = xgb_regressor.train(train_without_null)\nxgb_test_pred = xgb_regressor.predict(test)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:13.468437Z","iopub.execute_input":"2024-12-01T15:42:13.468793Z","iopub.status.idle":"2024-12-01T15:42:23.266516Z","shell.execute_reply.started":"2024-12-01T15:42:13.468760Z","shell.execute_reply":"2024-12-01T15:42:23.264993Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## CatBoost\n","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# class MultiTargetQWK(MultiTargetCustomMetric):\n#     def get_final_error(self, error, weight):\n#         return np.sum(error)  # / np.sum(weight)\n\n#     def is_max_optimal(self):\n#         # if True, the bigger the better\n#         return True\n\n#     def evaluate(self, approxes, targets, weight):\n#         # approxes: 予測値 (shape: [ターゲット数, サンプル数])\n#         # targets: 実際の値 (shape: [ターゲット数, サンプル数])\n#         # weight: サンプルごとの重み (Noneも可)\n\n#         preds = np.array(np.array([a.tolist() for a in approxes[:20]])).T\n#         _targets = np.array(np.array([a.tolist() for a in targets[:20]])).T\n\n#         pred = convert_pciat_pciat_k_to_sii(preds)\n#         target = convert_pciat_pciat_k_to_sii(_targets)\n\n#         qwk = cohen_kappa_score(pred, target, weights='quadratic')\n#         return qwk, 1\n\n#     def get_custom_metric_name(self):\n#         return \"MultiTargetQWK\"\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:23.267967Z","iopub.execute_input":"2024-12-01T15:42:23.268330Z","iopub.status.idle":"2024-12-01T15:42:23.274254Z","shell.execute_reply.started":"2024-12-01T15:42:23.268295Z","shell.execute_reply":"2024-12-01T15:42:23.273020Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class SiiQWK(MultiTargetCustomMetric):\n    def get_final_error(self, error, weight):\n        return np.sum(error)  # / np.sum(weight)\n\n    def is_max_optimal(self):\n        # if True, the bigger the better\n        return True\n\n    def evaluate(self, approxes, targets, weight):\n        # approxes: 予測値 (shape: [ターゲット数, サンプル数])\n        # targets: 実際の値 (shape: [ターゲット数, サンプル数])\n        # weight: サンプルごとの重み (Noneも可)\n        predt = np.array([a.tolist() for a in approxes[-1]])\n        pred = np.clip(predt, 0, 3).round().astype(int)\n        \n        target = np.array([a.tolist() for a in targets[-1]])\n\n        qwk = cohen_kappa_score(pred, target, weights='quadratic')\n        return qwk, 1\n\n    def get_custom_metric_name(self):\n        return \"SiiQWK\"\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:23.276114Z","iopub.execute_input":"2024-12-01T15:42:23.276542Z","iopub.status.idle":"2024-12-01T15:42:23.291990Z","shell.execute_reply.started":"2024-12-01T15:42:23.276500Z","shell.execute_reply":"2024-12-01T15:42:23.290810Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"CAT_params = dict(\n    loss_function=\"MultiRMSE\",\n    eval_metric=SiiQWK(),\n    iterations=1 if IS_TEST else 10000,\n    learning_rate=0.02,\n    depth=7,\n    early_stopping_rounds=100,\n)\n\nclass CatBoostSiiPredicter:\n    \n    def train(self, train_df: pl.DataFrame):\n\n        preds = []\n\n        for fold in range(CFG.N_FOLD):\n            print(f\"# CatBoost Fold {fold}\")\n            train = train_df.filter(\n                pl.col(\"fold\") != fold\n            )\n            val = train_df.filter(\n                pl.col(\"fold\") == fold\n            )\n            \n            X_train = train[FEATURE_COLS + TIME_SERIES_COLS]\n            X_val = val[FEATURE_COLS + TIME_SERIES_COLS]\n            y_train = train[TARGET_COLS]\n            y_val = val[TARGET_COLS]\n             \n            # train model\n            model = CatBoostRegressor(**CAT_params)\n            model.fit(\n                X_train.to_pandas(),\n                y_train.to_pandas(),\n                eval_set=(X_val.to_pandas(), y_val.to_pandas()),\n                cat_features=cat_features,\n                verbose=50,\n            )\n            \n            model.save_model(f\"catboost_fold_{fold}.json\")\n            \n            pciat_pred = model.predict(X_val.to_pandas())\n            df = pl.DataFrame(pciat_pred, schema=[name+\"_cat\" for name in TARGET_COLS]).with_columns(\n                (val[\"id\"]).alias(\"id\"),\n                (val[\"sii\"]).alias(\"sii\"),\n            )\n            preds.append(df)\n\n        y_pred = pl.concat(preds)\n\n        # print cv score\n        # print cv score\n        y_pred_sii = convert_value_to_sii(y_pred[\"sii_cat\"].to_numpy())\n        y_sii = y_pred[\"sii\"]\n\n        qwk = cohen_kappa_score(y_sii, y_pred_sii, weights=\"quadratic\")\n        \n        print(f\"[CatBoost] CV QWK Score: {qwk}\")\n\n        return y_pred\n    \n    def predict(self, test_df: pl.DataFrame):\n        X = test_df[FEATURE_COLS + TIME_SERIES_COLS].to_pandas()\n        preds: list[NDArray[np.int_]] = []\n\n        for fold in range(CFG.N_FOLD):\n            model = CatBoostRegressor(**CAT_params)\n            model.load_model(f\"catboost_fold_{fold}.json\")            \n            \n            pred = model.predict(X)\n            preds.append(pred)\n\n        return np.mean(preds, axis=0)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:23.294126Z","iopub.execute_input":"2024-12-01T15:42:23.294665Z","iopub.status.idle":"2024-12-01T15:42:23.311255Z","shell.execute_reply.started":"2024-12-01T15:42:23.294607Z","shell.execute_reply":"2024-12-01T15:42:23.310109Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"cat_regressor = CatBoostSiiPredicter()\ncat_val_df = cat_regressor.train(train_without_null)\ncat_test_pred = cat_regressor.predict(test)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:23.312775Z","iopub.execute_input":"2024-12-01T15:42:23.313252Z","iopub.status.idle":"2024-12-01T15:42:28.675261Z","shell.execute_reply.started":"2024-12-01T15:42:23.313205Z","shell.execute_reply":"2024-12-01T15:42:28.674104Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Tabnet","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_pd = train.to_pandas()\ntest_pd = test.to_pandas()\n\nimputer = KNNImputer(n_neighbors=5)\nnumeric_cols = train_pd.select_dtypes(include=['float32', 'float64', 'int32', 'int64']).columns.to_list()\n\n# targetに関するカラムはimputeから除外\nfor col in TARGET_COLS:\n    numeric_cols.remove(col)\n\n\n# 無限大を NaN に置き換え\ntrain_pd[numeric_cols] = train_pd[numeric_cols].replace([np.inf, -np.inf], np.nan)\ntest_pd[numeric_cols] = test_pd[numeric_cols].replace([np.inf, -np.inf], np.nan)\n\nimputer.fit(pd.concat([\n    train_pd[numeric_cols],\n    test_pd[numeric_cols],    \n]))\n\n# knnで欠損値を補完\nimputed_data = imputer.transform(pd.concat([\n    train_pd[numeric_cols],\n    test_pd[numeric_cols],    \n]))\ntrain_imputed_data = imputed_data[:len(train_pd), ]\ntest_imputed_data = imputed_data[len(train_pd):, ]\n\ntrain_imputed = pd.DataFrame(train_imputed_data, columns=numeric_cols)\ntest_imputed = pd.DataFrame(test_imputed_data, columns=numeric_cols)\n\n\nfor col in train_pd.columns:\n    if col not in numeric_cols:\n        train_imputed[col] = train_pd[col]\n\nfor col in test.columns:\n    if col not in numeric_cols:\n        test_imputed[col] = test_pd[col]\n\ntrain_pd = train_imputed\ntest_pd = test_imputed","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:28.676636Z","iopub.execute_input":"2024-12-01T15:42:28.676962Z","iopub.status.idle":"2024-12-01T15:42:44.471966Z","shell.execute_reply.started":"2024-12-01T15:42:28.676931Z","shell.execute_reply":"2024-12-01T15:42:44.471104Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# すべて0であることを確認\nassert train_pd[FEATURE_COLS].isnull().all().sum() == 0, \"train has null value\"\nassert test_pd[FEATURE_COLS].isnull().all().sum() == 0,  \"test has null value\"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:44.473565Z","iopub.execute_input":"2024-12-01T15:42:44.474321Z","iopub.status.idle":"2024-12-01T15:42:44.483261Z","shell.execute_reply.started":"2024-12-01T15:42:44.474275Z","shell.execute_reply":"2024-12-01T15:42:44.482117Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# https://discuss.pytorch.org/t/rmse-loss-function/16540\n# https://github.com/dreamquark-ai/tabnet/blob/develop/customizing_example.ipynb\n\ndef MCRMSELoss(y_pred, y_true):\n    \"\"\"\n    Dummy example similar to using default torch.nn.functional.cross_entropy\n    \"\"\"\n    num_scored=21\n    score = 0\n    criterion = nn.MSELoss()\n\n    for i in range(num_scored):\n        score += (\n            torch.sqrt(criterion(y_pred[:, i], y_true[:, i])) / num_scored\n        )\n\n    return score","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:44.484758Z","iopub.execute_input":"2024-12-01T15:42:44.485211Z","iopub.status.idle":"2024-12-01T15:42:44.494201Z","shell.execute_reply.started":"2024-12-01T15:42:44.485164Z","shell.execute_reply":"2024-12-01T15:42:44.493090Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# QWK metric\n# ref: https://www.kaggle.com/code/mawanda/qwk-metric-and-loss-in-pytorch#Implementation-of-metric-function\n\nclass QWK(Metric):\n    def __init__(self):\n        self._name = \"QWK\"\n        self._maximize = True\n\n    def __call__(self, y_true, y_pred):\n        # Convert predictions and true labels to numpy arrays\n        # output = y_pred.cpu().numpy()\n        # target = y_true.cpu().numpy()\n\n        output = np.clip(y_pred[:, -1], 0, 3).round().astype(int)\n        target = y_true[:, -1]\n\n        # Calculate and return the QWK score\n        return self.quadratic_kappa_coefficient(output, target)\n    \n    def quadratic_kappa_coefficient(self, output, target):\n        # Calculate Quadratic Weighted Kappa\n        return cohen_kappa_score(output, target, weights='quadratic')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:44.495518Z","iopub.execute_input":"2024-12-01T15:42:44.496087Z","iopub.status.idle":"2024-12-01T15:42:44.511782Z","shell.execute_reply.started":"2024-12-01T15:42:44.495998Z","shell.execute_reply":"2024-12-01T15:42:44.510582Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class TabNetWrapper(BaseEstimator, RegressorMixin):\n    def __init__(self, **kwargs):\n        self.model = TabNetRegressor(**kwargs)\n        self.kwargs = kwargs\n        self.imputer = SimpleImputer(strategy='median')\n        \n    def fit(self, X_train, X_valid, y_train, y_valid, best_model_path):        \n        \n        # Train TabNet model\n        history = self.model.fit(\n            X_train=X_train,\n            y_train=y_train,\n            eval_set=[(X_valid, y_valid)],\n            eval_name=['val'],\n            eval_metric=[QWK],\n            max_epochs=1 if IS_TEST else 2000,\n            patience=100,\n            batch_size=1024,\n            virtual_batch_size=128,\n            num_workers=0,\n            drop_last=False,\n            loss_fn=MCRMSELoss,\n            callbacks=[\n                TabNetPretrainedModelCheckpoint(\n                    filepath=best_model_path,\n                    monitor='valid_mse',\n                    mode='max',\n                    save_best_only=True,\n                    verbose=10\n                )\n            ]\n        )\n\n        # Load the best model\n        if os.path.exists(best_model_path):\n            self.model.load_model(best_model_path)\n        \n        self.model.save_model(best_model_path)\n        \n        return self\n    \n    def predict(self, X):\n        if hasattr(X, 'values'):\n            X = X.values\n            \n        pred = self.model.predict(X)\n\n        return pred # PCIAT-PCIAT_kとsiiの両方返す\n\n    def load_model(self, model_path):\n        self.model.load_model(model_path)\n\n    def __deepcopy__(self, memo):\n        # Add deepcopy support for scikit-learn\n        cls = self.__class__\n        result = cls.__new__(cls)\n        memo[id(self)] = result\n        for k, v in self.__dict__.items():\n            setattr(result, k, deepcopy(v, memo))\n        return result","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:44.513429Z","iopub.execute_input":"2024-12-01T15:42:44.513944Z","iopub.status.idle":"2024-12-01T15:42:44.528101Z","shell.execute_reply.started":"2024-12-01T15:42:44.513892Z","shell.execute_reply":"2024-12-01T15:42:44.526920Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# TabNet hyperparameters\nTabNet_Params = {\n    'n_d': 64,              # Width of the decision prediction layer\n    'n_a': 64,              # Width of the attention embedding for each step\n    'n_steps': 5,           # Number of steps in the architecture\n    'gamma': 1.5,           # Coefficient for feature selection regularization\n    'n_independent': 2,     # Number of independent GLU layer in each GLU block\n    'n_shared': 2,          # Number of shared GLU layer in each GLU block\n    'lambda_sparse': 1e-4,  # Sparsity regularization\n    'optimizer_fn': torch.optim.Adam,\n    'optimizer_params': dict(lr=2e-2, weight_decay=1e-5),\n    'mask_type': 'entmax',\n    'scheduler_params': dict(mode=\"max\", patience=10, min_lr=1e-5, factor=0.5),\n    'scheduler_fn': torch.optim.lr_scheduler.ReduceLROnPlateau,\n    'verbose': 10,\n    'device_name': 'cuda' if torch.cuda.is_available() else 'cpu'\n}\n\nclass TabNetPretrainedModelCheckpoint(Callback):\n    def __init__(self, filepath, monitor='val_loss', mode='min', \n                 save_best_only=True, verbose=1):\n        super().__init__()  # Initialize parent class\n        self.filepath = filepath\n        self.monitor = monitor\n        self.mode = mode\n        self.save_best_only = save_best_only\n        self.verbose = verbose\n        self.best = float('inf') if mode == 'min' else -float('inf')\n        \n    def on_train_begin(self, logs=None):\n        self.model = self.trainer  # Use trainer itself as model\n        \n    def on_epoch_end(self, epoch, logs=None):\n        logs = logs or {}\n        current = logs.get(self.monitor)\n        if current is None:\n            return\n        \n        # Check if current metric is better than best\n        if (self.mode == 'min' and current < self.best) or \\\n           (self.mode == 'max' and current > self.best):\n            if self.verbose:\n                print(f'\\nEpoch {epoch}: {self.monitor} improved from {self.best:.4f} to {current:.4f}')\n            self.best = current\n            if self.save_best_only:\n                self.model.save_model(self.filepath)  # Save the entire model","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:44.529730Z","iopub.execute_input":"2024-12-01T15:42:44.530245Z","iopub.status.idle":"2024-12-01T15:42:44.552356Z","shell.execute_reply.started":"2024-12-01T15:42:44.530193Z","shell.execute_reply":"2024-12-01T15:42:44.551286Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class Trainer:\n    def train(\n        self, \n        train_df: pd.DataFrame, \n    ):\n        preds = []\n\n        for fold in range(CFG.N_FOLD):\n            print(f\"# TabNet Fold {fold}\")\n            train = train_df[train_df[\"fold\"] != fold]\n            val   = train_df[train_df[\"fold\"] == fold]\n            \n            X_train = train[FEATURE_COLS + TIME_SERIES_COLS]\n            X_val = val[FEATURE_COLS + TIME_SERIES_COLS]\n            y_train = train[TARGET_COLS]\n            y_val = val[TARGET_COLS]\n\n            if hasattr(X_val, 'values'):\n                X_train = X_train.values\n                X_val = X_val.values\n            \n            if hasattr(y_val, 'values'):\n                y_train = y_train.values\n                y_val = y_val.values\n\n            # train model\n            model = TabNetWrapper(**TabNet_Params)\n\n            best_model_path = f\"tabnet_fold_{fold}.pt\"\n            model.fit(X_train, X_val, y_train, y_val, best_model_path)\n            \n            pciat_pred = model.predict(X_val)\n\n            df = pl.DataFrame(pciat_pred, schema=[name+\"_tabnet\" for name in TARGET_COLS]).with_columns(\n                (pl.from_pandas(val)[\"id\"]).alias(\"id\"),\n                (pl.from_pandas(val)[\"sii\"]).alias(\"sii\"),\n            )\n            preds.append(df)\n\n        y_pred = pl.concat(preds)\n    \n        # print cv score\n        y_pred_sii = convert_value_to_sii(y_pred.select(\"sii_tabnet\").to_numpy())\n        y_sii = y_pred[\"sii\"]\n\n        qwk = cohen_kappa_score(y_sii, y_pred_sii, weights=\"quadratic\")\n        print(f\"[TabNet] CV QWK Score: {qwk}\")\n\n        return y_pred\n\n    def predict(self, test: pd.DataFrame):\n        X = test[FEATURE_COLS + TIME_SERIES_COLS]\n        preds: list[NDArray[np.int_]] = []\n\n        for fold in range(CFG.N_FOLD):\n            model = TabNetWrapper(**TabNet_Params)\n            \n            best_model_path = f\"tabnet_fold_{fold}.pt\"\n            model.load_model(best_model_path + \".zip\")\n\n            pred = model.predict(X)\n            preds.append(pred)\n\n        return np.mean(preds, axis=0)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:44.553800Z","iopub.execute_input":"2024-12-01T15:42:44.554170Z","iopub.status.idle":"2024-12-01T15:42:44.572711Z","shell.execute_reply.started":"2024-12-01T15:42:44.554136Z","shell.execute_reply":"2024-12-01T15:42:44.571578Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"tabnet_regressor = Trainer()\ntabnet_val_df = tabnet_regressor.train(train_pd.dropna(subset=TARGET_COLS))\ntabnet_test_pred = tabnet_regressor.predict(test_pd)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:44.574309Z","iopub.execute_input":"2024-12-01T15:42:44.574859Z","iopub.status.idle":"2024-12-01T15:42:53.029239Z","shell.execute_reply.started":"2024-12-01T15:42:44.574802Z","shell.execute_reply":"2024-12-01T15:42:53.028297Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# y_pred_sii = convert_value_to_sii(tabnet_val_df.select(\"sii_tabnet\").to_numpy())\n# y_sii = tabnet_val_df[\"sii\"]\n\n# qwk = cohen_kappa_score(y_sii, y_pred_sii, weights=\"quadratic\")\n# print(f\"[TabNet] CV QWK Score: {qwk}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:53.030739Z","iopub.execute_input":"2024-12-01T15:42:53.031514Z","iopub.status.idle":"2024-12-01T15:42:53.036542Z","shell.execute_reply.started":"2024-12-01T15:42:53.031467Z","shell.execute_reply":"2024-12-01T15:42:53.035386Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Check predict results","metadata":{}},{"cell_type":"code","source":"xgb_test_pred[:, -1]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:53.037840Z","iopub.execute_input":"2024-12-01T15:42:53.038163Z","iopub.status.idle":"2024-12-01T15:42:53.051201Z","shell.execute_reply.started":"2024-12-01T15:42:53.038131Z","shell.execute_reply":"2024-12-01T15:42:53.049994Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"cat_test_pred[:, -1]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:53.052649Z","iopub.execute_input":"2024-12-01T15:42:53.053117Z","iopub.status.idle":"2024-12-01T15:42:53.066886Z","shell.execute_reply.started":"2024-12-01T15:42:53.053044Z","shell.execute_reply":"2024-12-01T15:42:53.065539Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"tabnet_test_pred[:, -1]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:53.068423Z","iopub.execute_input":"2024-12-01T15:42:53.068882Z","iopub.status.idle":"2024-12-01T15:42:53.085206Z","shell.execute_reply.started":"2024-12-01T15:42:53.068836Z","shell.execute_reply":"2024-12-01T15:42:53.084020Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Stacking\n","metadata":{}},{"cell_type":"markdown","source":"## XGBoost","metadata":{}},{"cell_type":"code","source":"# for train\ntrain_with_pseudo_label = train.join(\n    cat_val_df,\n    on=[\"id\", \"sii\"],\n    how=\"left\"\n).join(\n    xgb_val_df,\n    on=[\"id\", \"sii\"],\n    how=\"left\"\n).join(\n    tabnet_val_df,\n    on=[\"id\", \"sii\"],\n    how=\"left\"\n)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:53.086833Z","iopub.execute_input":"2024-12-01T15:42:53.087219Z","iopub.status.idle":"2024-12-01T15:42:53.103220Z","shell.execute_reply.started":"2024-12-01T15:42:53.087187Z","shell.execute_reply":"2024-12-01T15:42:53.102113Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# for test\ncat_test = pl.DataFrame(\n    cat_test_pred, \n    schema=[name+\"_cat\" for name in TARGET_COLS]\n).with_columns(\n    (test[\"id\"]).alias(\"id\"),\n)\nxgb_test = pl.DataFrame(\n    xgb_test_pred, \n    schema=[name+\"_xgb\" for name in TARGET_COLS]\n).with_columns(\n    (test[\"id\"]).alias(\"id\"),\n)\ntabnet_test = pl.DataFrame(\n    tabnet_test_pred, \n    schema=[name+\"_tabnet\" for name in TARGET_COLS]\n).with_columns(\n    (test[\"id\"]).alias(\"id\"),\n)\n\ntest_with_pseudo_label = test.join(\n    cat_test,\n    on=[\"id\"],\n    how=\"left\"\n).join(\n    xgb_test,\n    on=[\"id\"],\n    how=\"left\"\n).join(\n    tabnet_test,\n    on=[\"id\"],\n    how=\"left\"\n)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:53.104983Z","iopub.execute_input":"2024-12-01T15:42:53.105377Z","iopub.status.idle":"2024-12-01T15:42:53.122024Z","shell.execute_reply.started":"2024-12-01T15:42:53.105344Z","shell.execute_reply":"2024-12-01T15:42:53.120628Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# setting parameters\nXGB_Stacking_Params = {\n    'learning_rate': 0.05,\n    'max_depth': 7,\n    'n_estimators': 1 if IS_TEST else 10000,\n    'subsample': 0.8,\n    'colsample_bytree': 0.8,\n    'reg_alpha': 1,  # Increased from 0.1\n    'reg_lambda': 5,  # Increased from 1\n    'random_state': CFG.random_seed,\n}\n\nclass XGBStackingPredicter:\n    def __init__(self):\n        self.xgb_pseudo_cols = [name+\"_xgb\" for name in TARGET_COLS]\n        self.cat_pseudo_cols = [name+\"_cat\" for name in TARGET_COLS]\n        self.tabnet_pseudo_cols = [name+\"_tabnet\" for name in TARGET_COLS]\n        self.pseudo_cols = self.xgb_pseudo_cols + self.cat_pseudo_cols + self.tabnet_pseudo_cols\n\n    def quadratic_kappa_coefficient(self, output, target):\n        # Calculate Quadratic Weighted Kappa\n        return cohen_kappa_score(output, target, weights='quadratic')\n    \n    def qwk(self, dtrain: xgb.DMatrix, predt: np.ndarray) -> tuple[str, float]:\n        output = pred_sii = np.clip(predt[:, -1], 0, 3).round().astype(int)\n        convert_pciat_pciat_k_to_sii(predt)\n        \n        target = dtrain.reshape(predt.shape)[:, -1]\n        # target = convert_pciat_pciat_k_to_sii(targets)\n        # Calculate and return the QWK score\n        return -self.quadratic_kappa_coefficient(output, target)\n    \n    def train(self, train_df: pl.DataFrame):\n\n        preds = []\n\n        for fold in range(CFG.N_FOLD):\n            print(f\"## Stacking XGB fold {fold} ####\")\n\n            train = train_df.filter(\n                pl.col(\"fold\") != fold\n            )\n            val = train_df.filter(\n                pl.col(\"fold\") == fold\n            )\n            \n            X_train = train[FEATURE_COLS + TIME_SERIES_COLS + self.pseudo_cols].to_pandas()\n            X_val = val[FEATURE_COLS + TIME_SERIES_COLS + self.pseudo_cols].to_pandas()\n            X_train[cat_features] = X_train[cat_features].astype(\"category\")\n            X_val[cat_features]   = X_val[cat_features].astype(\"category\")\n\n            y_train = train[TARGET_COLS].to_pandas()\n            y_val   = val[TARGET_COLS].to_pandas()\n            \n            Xy = xgb.DMatrix(X_train, y_train)\n            evalXy = xgb.DMatrix(X_val, y_val)\n            \n            # train model\n            model = XGBRegressor(\n                **XGB_Stacking_Params,\n                enable_categorical=True,\n                tree_method=\"hist\",\n                multi_strategy=\"multi_output_tree\",\n                objective=\"reg:squarederror\",\n                eval_metric=self.qwk,\n                early_stopping_rounds=100,\n            )\n            model.fit(\n                X=X_train, \n                y=y_train, \n                eval_set=[(X_val, y_val)],\n                verbose=10\n            )\n            \n            file_name = f\"./xgboost_stacking_fold_{fold}.json\"\n            model.save_model(file_name)\n            \n            pciat_pred = model.predict(X_val)\n            df = pl.DataFrame(pciat_pred, schema=[name+\"_stacking_xgb\" for name in TARGET_COLS]).with_columns(\n                (val[\"id\"]).alias(\"id\"),\n                (val[\"sii\"]).alias(\"sii\"),\n            )\n            preds.append(df)\n\n        y_pred = pl.concat(preds)\n\n        # print cv score\n        y_pred_sii = convert_value_to_sii(y_pred[\"sii_stacking_xgb\"].to_numpy())\n        y_sii = y_pred[\"sii\"]\n\n        qwk = cohen_kappa_score(y_sii, y_pred_sii, weights=\"quadratic\")\n        print(f\"Cross-Validated QWK Score: {qwk}\")\n\n        return y_pred\n    \n    def predict(self, test_df: pl.DataFrame):\n        _test_df = test_df[FEATURE_COLS + TIME_SERIES_COLS + self.pseudo_cols].to_pandas()\n        _test_df[cat_features] = _test_df[cat_features].astype(\"category\")\n        \n        X = xgb.DMatrix(_test_df)\n        preds: list[NDArray[np.int_]] = []\n\n        for fold in range(CFG.N_FOLD):\n            model = xgb.Booster()\n            model.load_model(f\"xgboost_stacking_fold_{fold}.json\")\n            \n            pred = model.predict(X)\n            preds.append(pred)\n\n        return np.mean(preds, axis=0)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:53.129606Z","iopub.execute_input":"2024-12-01T15:42:53.129984Z","iopub.status.idle":"2024-12-01T15:42:53.149154Z","shell.execute_reply.started":"2024-12-01T15:42:53.129951Z","shell.execute_reply":"2024-12-01T15:42:53.147879Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"stacking_model = XGBStackingPredicter()\nxgb_stacking_val_df = stacking_model.train(train_with_pseudo_label.drop_nulls(subset=TARGET_COLS))\nxgb_stacking_test_pred = stacking_model.predict(test_with_pseudo_label)[:, -1]\n# xgb_stacking_test_sii = convert_value_to_sii(xgb_stacking_test_pred)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:42:53.150605Z","iopub.execute_input":"2024-12-01T15:42:53.150960Z","iopub.status.idle":"2024-12-01T15:43:09.129175Z","shell.execute_reply.started":"2024-12-01T15:42:53.150927Z","shell.execute_reply":"2024-12-01T15:43:09.128257Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Catboost","metadata":{}},{"cell_type":"code","source":"class CatBoostStackingPredicter:\n    def __init__(self):\n        self.xgb_pseudo_cols = [name+\"_xgb\" for name in TARGET_COLS]\n        self.cat_pseudo_cols = [name+\"_cat\" for name in TARGET_COLS]\n        self.tabnet_pseudo_cols = [name+\"_tabnet\" for name in TARGET_COLS]\n        self.pseudo_cols = self.xgb_pseudo_cols + self.cat_pseudo_cols + self.tabnet_pseudo_cols\n    \n    def train(self, train_df: pl.DataFrame):\n\n        preds = []\n\n        for fold in range(CFG.N_FOLD):\n            print(f\"# CatBoost Fold {fold}\")\n            train = train_df.filter(\n                pl.col(\"fold\") != fold\n            )\n            val = train_df.filter(\n                pl.col(\"fold\") == fold\n            )\n            \n            X_train = train[FEATURE_COLS + TIME_SERIES_COLS + self.pseudo_cols]\n            X_val = val[FEATURE_COLS + TIME_SERIES_COLS + self.pseudo_cols]\n            y_train = train[TARGET_COLS]\n            y_val = val[TARGET_COLS]\n             \n            # train model\n            model = CatBoostRegressor(**CAT_params)\n            model.fit(\n                X_train.to_pandas(),\n                y_train.to_pandas(),\n                eval_set=(X_val.to_pandas(), y_val.to_pandas()),\n                cat_features=cat_features,\n                verbose=50,\n            )\n            \n            model.save_model(f\"catboost_stacking_fold_{fold}.json\")\n\n            pciat_pred = model.predict(X_val.to_pandas())\n            df = pl.DataFrame(pciat_pred, schema=[name+\"_stacking_cat\" for name in TARGET_COLS]).with_columns(\n                (val[\"id\"]).alias(\"id\"),\n                (val[\"sii\"]).alias(\"sii\"),\n            )\n            preds.append(df)\n\n        y_pred = pl.concat(preds)\n        \n        # print cv score\n        y_pred_sii = convert_value_to_sii(y_pred[\"sii_stacking_cat\"].to_numpy())\n        y_sii = y_pred[\"sii\"]\n\n        qwk = cohen_kappa_score(y_sii, y_pred_sii, weights=\"quadratic\")\n\n        print(f\"[CatBoost] CV QWK Score: {qwk}\")\n\n        return y_pred\n    \n    def predict(self, test_df: pl.DataFrame):\n        X = test_df[FEATURE_COLS + TIME_SERIES_COLS + self.pseudo_cols].to_pandas()\n        preds: list[NDArray[np.int_]] = []\n\n        for fold in range(CFG.N_FOLD):\n            model = CatBoostRegressor(**CAT_params)\n            model.load_model(f\"catboost_stacking_fold_{fold}.json\")   \n\n            pred = model.predict(X)\n            preds.append(pred)\n            \n        return np.mean(preds, axis=0)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:44:48.374516Z","iopub.execute_input":"2024-12-01T15:44:48.374998Z","iopub.status.idle":"2024-12-01T15:44:48.388442Z","shell.execute_reply.started":"2024-12-01T15:44:48.374963Z","shell.execute_reply":"2024-12-01T15:44:48.387279Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"cat_stacking_model = CatBoostStackingPredicter()\ncat_stacking_val_df = cat_stacking_model.train(train_with_pseudo_label.drop_nulls(subset=TARGET_COLS))\ncat_stacking_test_pred = cat_stacking_model.predict(test_with_pseudo_label)[:, -1]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:44:48.829974Z","iopub.execute_input":"2024-12-01T15:44:48.830938Z","iopub.status.idle":"2024-12-01T15:44:52.370738Z","shell.execute_reply.started":"2024-12-01T15:44:48.830898Z","shell.execute_reply":"2024-12-01T15:44:52.369578Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Tabnet","metadata":{}},{"cell_type":"code","source":"# for train\ntrain_pd_with_pseudo_label = train_pd.merge(\n    cat_val_df.to_pandas(),\n    on=[\"id\", \"sii\"],\n    how=\"left\"\n).merge(\n    xgb_val_df.to_pandas(),\n    on=[\"id\", \"sii\"],\n    how=\"left\"\n).merge(\n    tabnet_val_df.to_pandas(),\n    on=[\"id\", \"sii\"],\n    how=\"left\"\n)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:45:08.918613Z","iopub.execute_input":"2024-12-01T15:45:08.919073Z","iopub.status.idle":"2024-12-01T15:45:08.984385Z","shell.execute_reply.started":"2024-12-01T15:45:08.919019Z","shell.execute_reply":"2024-12-01T15:45:08.983160Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# for test\n\ntest_pd_with_pseudo_label = test_pd.merge(\n    cat_test.to_pandas(),\n    on=[\"id\"],\n    how=\"left\"\n).merge(\n    xgb_test.to_pandas(),\n    on=[\"id\"],\n    how=\"left\"\n).merge(\n    tabnet_test.to_pandas(),\n    on=[\"id\"],\n    how=\"left\"\n)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:45:09.405593Z","iopub.execute_input":"2024-12-01T15:45:09.405986Z","iopub.status.idle":"2024-12-01T15:45:09.442447Z","shell.execute_reply.started":"2024-12-01T15:45:09.405953Z","shell.execute_reply":"2024-12-01T15:45:09.441324Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"class TabnetStackingTrainer:\n    def __init__(self):\n        self.xgb_pseudo_cols = [name+\"_xgb\" for name in TARGET_COLS]\n        self.cat_pseudo_cols = [name+\"_cat\" for name in TARGET_COLS]\n        self.tabnet_pseudo_cols = [name+\"_tabnet\" for name in TARGET_COLS]\n        self.pseudo_cols = self.xgb_pseudo_cols + self.cat_pseudo_cols + self.tabnet_pseudo_cols\n\n    def train(\n        self, \n        train_df: pd.DataFrame, \n    ):\n        preds = []\n\n        for fold in range(CFG.N_FOLD):\n            print(f\"# TabNet Fold {fold}\")\n            train = train_df[train_df[\"fold\"] != fold]\n            val   = train_df[train_df[\"fold\"] == fold]\n            \n            X_train = train[FEATURE_COLS + TIME_SERIES_COLS + self.pseudo_cols].astype(float)\n            X_val = val[FEATURE_COLS + TIME_SERIES_COLS + self.pseudo_cols].astype(float)\n            y_train = train[TARGET_COLS].astype(float)\n            y_val = val[TARGET_COLS].astype(float)\n\n            if hasattr(X_val, 'values'):\n                X_train = X_train.values\n                X_val = X_val.values\n            \n            if hasattr(y_val, 'values'):\n                y_train = y_train.values\n                y_val = y_val.values\n\n            # train model\n            model = TabNetWrapper(**TabNet_Params)\n\n            best_model_path = f\"tabnet_stacking_fold_{fold}.pt\"\n            model.fit(X_train, X_val, y_train, y_val, best_model_path)\n            \n            pciat_pred = model.predict(X_val)\n\n            df = pl.DataFrame(pciat_pred, schema=[name+\"_stacking_tabnet\" for name in TARGET_COLS]).with_columns(\n                pl.from_pandas(val[\"id\"]).alias(\"id\"),\n                pl.from_pandas(val[\"sii\"]).alias(\"sii\"),\n            )\n            preds.append(df)\n\n        y_pred = pl.concat(preds)\n\n        # print cv score\n        y_pred_sii = convert_value_to_sii(y_pred[\"sii_stacking_tabnet\"].to_numpy())\n        y_sii = y_pred[\"sii\"]\n        qwk = cohen_kappa_score(y_sii, y_pred_sii, weights=\"quadratic\")\n\n        print(f\"[TabNet] CV QWK Score: {qwk}\")\n\n        return y_pred\n\n        y_pred = pl.concat(preds)\n\n    def predict(self, test: pd.DataFrame):\n        X = test[FEATURE_COLS + TIME_SERIES_COLS + self.pseudo_cols].astype(float)\n        preds: list[NDArray[np.int_]] = []\n\n        for fold in range(CFG.N_FOLD):\n            model = TabNetWrapper(**TabNet_Params)\n            \n            best_model_path = f\"tabnet_stacking_fold_{fold}.pt\"\n            model.load_model(best_model_path + \".zip\")\n\n            pred = model.predict(X)\n            preds.append(pred)\n\n        return np.mean(preds, axis=0)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:47:01.053938Z","iopub.execute_input":"2024-12-01T15:47:01.054504Z","iopub.status.idle":"2024-12-01T15:47:01.072261Z","shell.execute_reply.started":"2024-12-01T15:47:01.054450Z","shell.execute_reply":"2024-12-01T15:47:01.071127Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"tabnet_stacking_model = TabnetStackingTrainer()\ntabnet_stacking_val_df = tabnet_stacking_model.train(train_pd_with_pseudo_label.dropna(subset=TARGET_COLS))\ntabnet_stacking_test_pred = tabnet_stacking_model.predict(test_pd_with_pseudo_label)[:, -1]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:47:01.272071Z","iopub.execute_input":"2024-12-01T15:47:01.272849Z","iopub.status.idle":"2024-12-01T15:47:08.618014Z","shell.execute_reply.started":"2024-12-01T15:47:01.272807Z","shell.execute_reply":"2024-12-01T15:47:08.616813Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"stacking_val_df = np.mean([xgb_stacking_val_df[\"sii_stacking_xgb\", \"sii\"].to_numpy(), cat_stacking_val_df[\"sii_stacking_cat\", \"sii\"].to_numpy(), tabnet_stacking_val_df[\"sii_stacking_tabnet\", \"sii\"].to_numpy()], axis=0)\nstacking_pred_df = np.mean([xgb_stacking_test_pred, cat_stacking_test_pred, tabnet_stacking_test_pred], axis=0)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:48:55.309221Z","iopub.execute_input":"2024-12-01T15:48:55.310281Z","iopub.status.idle":"2024-12-01T15:48:55.317691Z","shell.execute_reply.started":"2024-12-01T15:48:55.310241Z","shell.execute_reply":"2024-12-01T15:48:55.316311Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"stacking_pred_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:49:16.783347Z","iopub.execute_input":"2024-12-01T15:49:16.783764Z","iopub.status.idle":"2024-12-01T15:49:16.790836Z","shell.execute_reply.started":"2024-12-01T15:49:16.783729Z","shell.execute_reply":"2024-12-01T15:49:16.789720Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Optimizer\n","metadata":{}},{"cell_type":"code","source":"class OptimizedRounder:\n    \"\"\"\n    A class for optimizing the rounding of continuous predictions into discrete class labels using Optuna.\n    The optimization process maximizes the Quadratic Weighted Kappa score by learning thresholds that separate\n    continuous predictions into class intervals.\n\n    Args:\n        n_classes (int): The number of discrete class labels.\n        n_trials (int, optional): The number of trials for the Optuna optimization. Defaults to 100.\n\n    Attributes:\n        n_classes (int): The number of discrete class labels.\n        labels (NDArray[np.int_]): An array of class labels from 0 to `n_classes - 1`.\n        n_trials (int): The number of optimization trials.\n        metric (Callable): The Quadratic Weighted Kappa score metric used for optimization.\n        thresholds (List[float]): The optimized thresholds learned after calling `fit()`.\n\n    Methods:\n        fit(y_pred: NDArray[np.float_], y_true: NDArray[np.int_]) -> None:\n            Fits the rounding thresholds based on continuous predictions and ground truth labels.\n\n            Args:\n                y_pred (NDArray[np.float_]): Continuous predictions that need to be rounded.\n                y_true (NDArray[np.int_]): Ground truth class labels.\n\n            Returns:\n                None\n\n        predict(y_pred: NDArray[np.float_]) -> NDArray[np.int_]:\n            Predicts discrete class labels by rounding continuous predictions using the fitted thresholds.\n            `fit()` must be called before `predict()`.\n\n            Args:\n                y_pred (NDArray[np.float_]): Continuous predictions to be rounded.\n\n            Returns:\n                NDArray[np.int_]: Predicted class labels.\n\n        _normalize(y: NDArray[np.float_]) -> NDArray[np.float_]:\n            Normalizes the continuous values to the range [0, `n_classes - 1`].\n\n            Args:\n                y (NDArray[np.float_]): Continuous values to be normalized.\n\n            Returns:\n                NDArray[np.float_]: Normalized values.\n\n    References:\n        - This implementation uses Optuna for threshold optimization.\n        - Quadratic Weighted Kappa is used as the evaluation metric.\n    \"\"\"\n\n    def __init__(self, n_classes: int, n_trials: int = 100):\n        self.n_classes = n_classes\n        self.labels = np.arange(n_classes)\n        self.n_trials = n_trials\n        self.metric = partial(cohen_kappa_score, weights=\"quadratic\")\n\n    def fit(self, y_pred: NDArray[np.float_], y_true: NDArray[np.int_]) -> None:\n        y_pred = self._normalize(y_pred)\n\n        def objective(trial: optuna.Trial) -> float:\n            thresholds = []\n            for i in range(self.n_classes - 1):\n                low = max(thresholds) if i > 0 else min(self.labels)\n                high = max(self.labels)\n                th = trial.suggest_float(f\"threshold_{i}\", low, high)\n                thresholds.append(th)\n            try:\n                y_pred_rounded = np.digitize(y_pred, thresholds)\n            except ValueError:\n                return -100\n            return self.metric(y_true, y_pred_rounded)\n\n        optuna.logging.disable_default_handler()\n        study = optuna.create_study(direction=\"maximize\")\n        study.optimize(\n            objective,\n            n_trials=self.n_trials,\n        )\n        self.thresholds = [study.best_params[f\"threshold_{i}\"] for i in range(self.n_classes - 1)]\n\n    def predict(self, y_pred: NDArray[np.float_]) -> NDArray[np.int_]:\n        assert hasattr(self, \"thresholds\"), \"fit() must be called before predict()\"\n        y_pred = self._normalize(y_pred)\n        return np.digitize(y_pred, self.thresholds)\n\n    def _normalize(self, y: NDArray[np.float_]) -> NDArray[np.float_]:\n        # normalize y_pred to [0, n_classes - 1]\n        return (y - y.min()) / (y.max() - y.min()) * (self.n_classes - 1)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:55:36.663213Z","iopub.execute_input":"2024-12-01T15:55:36.664079Z","iopub.status.idle":"2024-12-01T15:55:36.679011Z","shell.execute_reply.started":"2024-12-01T15:55:36.664017Z","shell.execute_reply":"2024-12-01T15:55:36.677585Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"optimizer = OptimizedRounder(n_classes=4, n_trials=300)\ny_pred = stacking_val_df[:, 0]\ny_true = stacking_val_df[:, 1]\noptimizer.fit(y_pred, y_true)\ny_pred_rounded = optimizer.predict(y_pred)\n\n# Calculate QWK\nqwk = cohen_kappa_score(y_true, y_pred_rounded, weights=\"quadratic\")\nprint(f\"Cross-Validated QWK Score: {qwk}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:55:36.862814Z","iopub.execute_input":"2024-12-01T15:55:36.863767Z","iopub.status.idle":"2024-12-01T15:55:43.045377Z","shell.execute_reply.started":"2024-12-01T15:55:36.863714Z","shell.execute_reply":"2024-12-01T15:55:43.044156Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"test_pred_rounded = optimizer.predict(stacking_pred_df)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:55:44.305934Z","iopub.execute_input":"2024-12-01T15:55:44.306368Z","iopub.status.idle":"2024-12-01T15:55:44.312042Z","shell.execute_reply.started":"2024-12-01T15:55:44.306331Z","shell.execute_reply":"2024-12-01T15:55:44.310901Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Submission","metadata":{}},{"cell_type":"code","source":"# test_pred_sii = convert_pciat_pciat_k_to_sii(test_pred)\n\ntest.select(\"id\").with_columns(\n    pl.Series(\"sii\", pl.Series(\"sii\", test_pred_rounded)),\n).write_csv(\"submission.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:55:45.461359Z","iopub.execute_input":"2024-12-01T15:55:45.461822Z","iopub.status.idle":"2024-12-01T15:55:45.469035Z","shell.execute_reply.started":"2024-12-01T15:55:45.461785Z","shell.execute_reply":"2024-12-01T15:55:45.467712Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"test.select(\"id\").with_columns(\n    pl.Series(\"sii\", pl.Series(\"sii\", test_pred_rounded)),\n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-12-01T15:55:46.132724Z","iopub.execute_input":"2024-12-01T15:55:46.133198Z","iopub.status.idle":"2024-12-01T15:55:46.142919Z","shell.execute_reply.started":"2024-12-01T15:55:46.133161Z","shell.execute_reply":"2024-12-01T15:55:46.141562Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}