{"cells":[{"cell_type":"markdown","metadata":{},"source":"# CMI-BFRB IMU推論ノートブック v10\n\nIMU専用モデルによる推論（周波数特徴量強化版）。\n\n**Version**: v10.0.0 \n\n**モデル**: IMU Only Model with Frequency Features\n\n**改良点**:\n- 周波数領域特徴量（FFT、パワースペクトラム）\n- スペクトログラム特徴量\n- 変化点検出\n- より洗練された後処理\n\n**目標**: LB 0.70+"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"import os\nimport sys\nimport pandas as pd\nimport numpy as np\nimport polars as pl\nfrom pathlib import Path\nimport pickle\nimport warnings\nfrom scipy import signal\nfrom scipy.stats import kurtosis, skew\nwarnings.filterwarnings('ignore')\n\n# Kaggle環境でのパス設定\nif os.path.exists('/kaggle/input'):\n    DATA_DIR = Path('/kaggle/input/cmi-detect-behavior-with-sensor-data')\n    UTILS_DIR = Path('/kaggle/input/cmi-bfrb-utils-v131')\n    MODEL_DIR = Path('/kaggle/input/cmi-bfrb-imu-models-v10')\n    OUTPUT_DIR = Path('/kaggle/working')\nelse:\n    # ローカル環境（テスト用）\n    DATA_DIR = Path('../../experiments/data')\n    UTILS_DIR = Path('../../baseline/src')\n    MODEL_DIR = Path('../../models/imu_only')\n    OUTPUT_DIR = Path('../../submissions')\n\n# utilsディレクトリをパスに追加\nsys.path.insert(0, str(UTILS_DIR))\n\nprint(\"環境設定完了\")\nprint(f\"MODEL_DIR: {MODEL_DIR}\")"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# utilsのインポート\nfrom utils import (\n    load_data,\n    create_features,\n    prepare_submission,\n    ALL_GESTURES\n)\nfrom imu_features import create_imu_features\n\nprint(\"モジュール読み込み完了\")"},{"cell_type":"markdown","metadata":{},"source":"## 周波数特徴量関数の定義"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"def extract_frequency_features(data, sampling_rate=20):\n    \"\"\"\n    周波数領域の特徴量を抽出\n    \n    Args:\n        data: 時系列データ (1次元配列)\n        sampling_rate: サンプリングレート (Hz)\n    \n    Returns:\n        特徴量の辞書\n    \"\"\"\n    features = {}\n    \n    # NaN処理\n    data = np.nan_to_num(data, nan=0.0)\n    \n    if len(data) < 10:\n        # データが短すぎる場合はダミー値を返す\n        return {\n            'fft_max_freq': 0, 'fft_max_power': 0,\n            'fft_mean_power': 0, 'fft_std_power': 0,\n            'spectral_centroid': 0, 'spectral_rolloff': 0,\n            'spectral_entropy': 0\n        }\n    \n    # FFT計算\n    fft_vals = np.fft.rfft(data)\n    fft_power = np.abs(fft_vals) ** 2\n    freqs = np.fft.rfftfreq(len(data), 1/sampling_rate)\n    \n    # 最大パワーの周波数\n    max_idx = np.argmax(fft_power)\n    features['fft_max_freq'] = freqs[max_idx]\n    features['fft_max_power'] = fft_power[max_idx]\n    \n    # パワースペクトラムの統計量\n    features['fft_mean_power'] = np.mean(fft_power)\n    features['fft_std_power'] = np.std(fft_power)\n    \n    # スペクトラル重心\n    features['spectral_centroid'] = np.sum(freqs * fft_power) / np.sum(fft_power)\n    \n    # スペクトラルロールオフ（85%のエネルギーが含まれる周波数）\n    cumsum_power = np.cumsum(fft_power)\n    rolloff_idx = np.searchsorted(cumsum_power, 0.85 * cumsum_power[-1])\n    features['spectral_rolloff'] = freqs[min(rolloff_idx, len(freqs)-1)]\n    \n    # スペクトラルエントロピー\n    normalized_power = fft_power / np.sum(fft_power)\n    spectral_entropy = -np.sum(normalized_power * np.log2(normalized_power + 1e-10))\n    features['spectral_entropy'] = spectral_entropy\n    \n    return features\n\n\ndef extract_change_point_features(data, penalty=1):\n    \"\"\"\n    変化点検出による特徴量抽出\n    \n    Args:\n        data: 時系列データ\n        penalty: 変化点検出のペナルティパラメータ\n    \n    Returns:\n        特徴量の辞書\n    \"\"\"\n    features = {}\n    \n    # NaN処理\n    data = np.nan_to_num(data, nan=0.0)\n    \n    if len(data) < 10:\n        return {'n_change_points': 0, 'mean_segment_length': len(data)}\n    \n    # 簡易的な変化点検出（差分の閾値ベース）\n    diff = np.abs(np.diff(data))\n    threshold = np.mean(diff) + 2 * np.std(diff)\n    change_points = np.where(diff > threshold)[0]\n    \n    features['n_change_points'] = len(change_points)\n    \n    # セグメント長の統計\n    if len(change_points) > 0:\n        segments = np.diff(np.concatenate([[0], change_points, [len(data)]]))\n        features['mean_segment_length'] = np.mean(segments)\n    else:\n        features['mean_segment_length'] = len(data)\n    \n    return features\n\n\ndef create_enhanced_imu_features(df, is_train=True):\n    \"\"\"\n    強化されたIMU特徴量の作成（周波数特徴量含む）\n    \"\"\"\n    # 基本的なIMU特徴量\n    features = create_imu_features(df, is_train)\n    \n    # IMUカラム\n    imu_cols = ['acc_x', 'acc_y', 'acc_z', 'rot_w', 'rot_x', 'rot_y', 'rot_z']\n    \n    # 各シーケンスごとに周波数特徴量を追加\n    freq_features = []\n    change_features = []\n    \n    for seq_id in df['sequence_id'].unique():\n        seq_data = df[df['sequence_id'] == seq_id]\n        seq_freq_feat = {}\n        seq_change_feat = {}\n        \n        for col in imu_cols:\n            if col in seq_data.columns:\n                # 周波数特徴量\n                freq_feat = extract_frequency_features(seq_data[col].values)\n                for feat_name, feat_val in freq_feat.items():\n                    seq_freq_feat[f'{col}_{feat_name}'] = feat_val\n                \n                # 変化点特徴量\n                change_feat = extract_change_point_features(seq_data[col].values)\n                for feat_name, feat_val in change_feat.items():\n                    seq_change_feat[f'{col}_{feat_name}'] = feat_val\n        \n        # 加速度の大きさに対する特徴量\n        acc_magnitude = np.sqrt(\n            seq_data['acc_x'].values**2 + \n            seq_data['acc_y'].values**2 + \n            seq_data['acc_z'].values**2\n        )\n        \n        mag_freq_feat = extract_frequency_features(acc_magnitude)\n        for feat_name, feat_val in mag_freq_feat.items():\n            seq_freq_feat[f'acc_magnitude_{feat_name}'] = feat_val\n        \n        freq_features.append(seq_freq_feat)\n        change_features.append(seq_change_feat)\n    \n    # DataFrameに変換して結合\n    freq_df = pd.DataFrame(freq_features)\n    change_df = pd.DataFrame(change_features)\n    \n    # 元の特徴量と結合\n    features = pd.concat([features, freq_df, change_df], axis=1)\n    \n    # NaN処理\n    features = features.fillna(0)\n    \n    return features"},{"cell_type":"markdown","metadata":{},"source":"## モデルの読み込み"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"import lightgbm as lgb\nfrom sklearn.preprocessing import LabelEncoder\n\n# baselineモジュールのエラーを回避\nimport sys\nclass DummyModule:\n    pass\n\nsys.modules['baseline'] = DummyModule()\nsys.modules['baseline.src'] = DummyModule()\nsys.modules['baseline.src.utils'] = DummyModule()\n\n# グローバル変数として保持\nmodels = []\ntrained_feature_cols = []\nlabel_encoder = None\ntraining_results = {}\n\n# モデルファイルの読み込み\ntry:\n    model_path = MODEL_DIR / 'imu_only_model.pkl'\n    results_path = MODEL_DIR / 'imu_training_results.pkl'\n    \n    if model_path.exists() and results_path.exists():\n        with open(model_path, 'rb') as f:\n            model_data = pickle.load(f)\n        \n        with open(results_path, 'rb') as f:\n            training_results = pickle.load(f)\n        \n        # モデルとその他のデータを抽出\n        if isinstance(model_data, dict):\n            if 'models' in model_data:\n                models = model_data['models']\n                trained_feature_cols = model_data.get('feature_columns', [])\n                label_encoder = model_data.get('label_encoder', None)\n            else:\n                models = [model_data]\n        elif isinstance(model_data, list):\n            models = model_data\n        else:\n            models = [model_data]\n            \n        print(f\"モデル読み込み成功\")\n        print(f\"モデル数: {len(models)}\")\n        print(f\"Mean CV Score: {training_results.get('mean_cv_score', 'N/A')}\")\n        \nexcept Exception as e:\n    print(f\"モデル読み込みエラー: {e}\")\n    \n# ラベルエンコーダーが無い場合は作成\nif label_encoder is None:\n    print(\"\\nラベルエンコーダーを作成中...\")\n    label_encoder = LabelEncoder()\n    label_encoder.fit(ALL_GESTURES)\n    print(f\"クラス数: {len(label_encoder.classes_)}\")"},{"cell_type":"markdown","metadata":{},"source":"## 後処理関数の定義（改良版）"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# 後処理用の定数\nCONFIDENCE_THRESHOLD = 0.40  # v10では少し下げる\nBFRB_BOOST_FACTOR = 1.20    # より強いブースト\nSMOOTH_WINDOW = 5           # より大きな窓サイズ\n\n# BFRB行動の定義（正確なラベル名）\nBFRB_BEHAVIORS = [\n    'Above ear - pull hair',\n    'Cheek - pinch skin', \n    'Eyebrow - pull hair',\n    'Eyelash - pull hair',\n    'Forehead - pull hairline',\n    'Forehead - scratch',\n    'Neck - pinch skin',\n    'Neck - scratch'\n]\n\n# 混同しやすいジェスチャーのマッピング（改良版）\nGESTURE_CONFUSION_MAP = {\n    'Neck - scratch': ['Forehead - scratch', 'Scratch knee/leg skin'],\n    'Above ear - pull hair': ['Eyebrow - pull hair', 'Eyelash - pull hair'],\n    'Cheek - pinch skin': ['Neck - pinch skin', 'Pinch knee/leg skin'],\n    'Text on phone': ['Write name on leg'],  # 似た手の動き\n}\n\n\ndef apply_enhanced_post_processing(predictions, confidence_scores, frequency_features=None):\n    \"\"\"\n    強化された後処理（周波数特徴量を考慮）\n    \"\"\"\n    predictions_adjusted = predictions.copy()\n    \n    # 1. BFRB行動の確率をブースト\n    for i, gesture in enumerate(label_encoder.classes_):\n        if gesture in BFRB_BEHAVIORS:\n            boost = BFRB_BOOST_FACTOR\n            \n            # 周波数特徴量がある場合、追加の調整\n            if frequency_features is not None:\n                # 高周波成分が多い場合、細かい動きを示すので追加ブースト\n                if 'spectral_rolloff' in frequency_features and frequency_features['spectral_rolloff'] > 5:\n                    boost *= 1.05\n            \n            predictions_adjusted[:, i] *= boost\n    \n    # 2. 低信頼度予測の調整\n    for idx in range(len(predictions_adjusted)):\n        max_prob = np.max(predictions_adjusted[idx])\n        \n        if max_prob < CONFIDENCE_THRESHOLD:\n            # 上位5つの予測を考慮\n            top_5_indices = np.argsort(predictions_adjusted[idx])[-5:]\n            top_5_gestures = [label_encoder.classes_[i] for i in top_5_indices]\n            \n            # BFRB行動が上位に含まれる場合、その確率を強化\n            bfrb_count = sum(1 for g in top_5_gestures if g in BFRB_BEHAVIORS)\n            if bfrb_count >= 2:  # 複数のBFRB行動が候補にある\n                for i, gesture in enumerate(top_5_gestures):\n                    if gesture in BFRB_BEHAVIORS:\n                        predictions_adjusted[idx, top_5_indices[i]] *= 1.15\n    \n    # 3. 混同しやすいジェスチャーの調整\n    for idx in range(len(predictions_adjusted)):\n        pred_idx = np.argmax(predictions_adjusted[idx])\n        pred_gesture = label_encoder.classes_[pred_idx]\n        \n        if pred_gesture in GESTURE_CONFUSION_MAP:\n            confused_gestures = GESTURE_CONFUSION_MAP[pred_gesture]\n            for conf_gesture in confused_gestures:\n                if conf_gesture in label_encoder.classes_:\n                    conf_idx = list(label_encoder.classes_).index(conf_gesture)\n                    prob_diff = predictions_adjusted[idx, pred_idx] - predictions_adjusted[idx, conf_idx]\n                    \n                    # 確率が近い場合の調整\n                    if prob_diff < 0.15:\n                        # BFRB行動を優先\n                        if conf_gesture in BFRB_BEHAVIORS and pred_gesture not in BFRB_BEHAVIORS:\n                            predictions_adjusted[idx, conf_idx] *= 1.1\n                        elif pred_gesture in BFRB_BEHAVIORS and conf_gesture not in BFRB_BEHAVIORS:\n                            predictions_adjusted[idx, pred_idx] *= 1.1\n    \n    # 正規化\n    predictions_adjusted = predictions_adjusted / predictions_adjusted.sum(axis=1, keepdims=True)\n    \n    return predictions_adjusted"},{"cell_type":"markdown","metadata":{},"source":"## 推論関数の定義（v10）"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"def predict(sequence: pl.DataFrame, demographics: pl.DataFrame) -> str:\n    \"\"\"\n    CMI評価API用の推論関数（v10 - 周波数特徴量版）\n    \"\"\"\n    try:\n        # Polars DataFrameをPandasに変換\n        test = sequence.to_pandas() if isinstance(sequence, pl.DataFrame) else sequence\n        \n        # IMU列の定義\n        imu_cols = ['acc_x', 'acc_y', 'acc_z', 'rot_w', 'rot_x', 'rot_y', 'rot_z']\n        \n        # IMUセンサーのチェック\n        has_imu = all(col in test.columns for col in imu_cols)\n        \n        if not has_imu:\n            # IMUセンサーがない場合はデフォルト予測\n            return 'Wave hello'\n        \n        # subjectを追加\n        if demographics is not None and not demographics.is_empty():\n            demo_pd = demographics.to_pandas() if isinstance(demographics, pl.DataFrame) else demographics\n            test['subject'] = demo_pd['subject'].iloc[0] if 'subject' in demo_pd.columns else 'test_subject'\n        else:\n            test['subject'] = 'test_subject'\n        \n        # sequence_idを追加（必要に応じて）\n        if 'sequence_id' not in test.columns:\n            test['sequence_id'] = 'test_seq'\n        \n        # IMUデータの抽出\n        meta_cols = ['sequence_id', 'subject']\n        test_imu = test[meta_cols + imu_cols].copy()\n        \n        # 強化されたIMU特徴量の作成（周波数特徴量含む）\n        test_features = create_enhanced_imu_features(test_imu, is_train=False)\n        \n        # 追加の特徴量計算\n        # 1. 加速度の大きさ\n        acc_magnitude = np.sqrt(\n            test['acc_x']**2 + test['acc_y']**2 + test['acc_z']**2\n        )\n        test_features['acc_magnitude_mean'] = acc_magnitude.mean()\n        test_features['acc_magnitude_std'] = acc_magnitude.std()\n        test_features['acc_magnitude_max'] = acc_magnitude.max()\n        \n        # 2. ジャーク（加速度の変化率）\n        for axis in ['x', 'y', 'z']:\n            acc_col = f'acc_{axis}'\n            jerk = test[acc_col].diff().abs()\n            test_features[f'jerk_{axis}_mean'] = jerk.mean()\n            test_features[f'jerk_{axis}_max'] = jerk.max()\n        \n        # 3. 回転エネルギー\n        rot_energy = test['rot_w']**2 + test['rot_x']**2 + test['rot_y']**2 + test['rot_z']**2\n        test_features['rot_energy_mean'] = rot_energy.mean()\n        test_features['rot_energy_std'] = rot_energy.std()\n        \n        # 学習時の特徴量に合わせる\n        # v10では新しい特徴量も含めるため、存在しない特徴量は0で埋める\n        all_feature_cols = list(test_features.columns)\n        for col in all_feature_cols:\n            if col not in trained_feature_cols and col not in ['sequence_id', 'subject']:\n                # 新しい特徴量は使用可能\n                pass\n        \n        # 必要な特徴量だけを選択（学習時の特徴量 + 新しい特徴量）\n        available_cols = [col for col in test_features.columns if col in trained_feature_cols or col.endswith(('_fft_max_freq', '_spectral_centroid', '_n_change_points'))]\n        \n        # 学習時の特徴量で存在しないものは0で埋める\n        for col in trained_feature_cols:\n            if col not in test_features.columns:\n                test_features[col] = 0\n        \n        X_test = test_features[trained_feature_cols]\n        \n        # 各モデルで予測\n        predictions = np.zeros((1, 18))\n        \n        # モデルごとの重み（v10用に調整）\n        model_weights = [1.0, 1.15, 0.95, 1.05, 1.1]\n        \n        for i, model in enumerate(models):\n            pred = model.predict(X_test, num_iteration=model.best_iteration)\n            weight = model_weights[i] if i < len(model_weights) else 1.0\n            predictions += pred * weight\n        \n        # 重み付き平均\n        total_weight = sum(model_weights[:len(models)])\n        predictions /= total_weight\n        \n        # 信頼度の計算\n        confidence_scores = np.array([np.max(predictions[0])])\n        \n        # 周波数特徴量の情報を集約\n        freq_info = {\n            'spectral_rolloff': test_features.get('acc_magnitude_spectral_rolloff', 0).iloc[0] if 'acc_magnitude_spectral_rolloff' in test_features else 0\n        }\n        \n        # 強化された後処理の適用\n        predictions_adjusted = apply_enhanced_post_processing(\n            predictions, confidence_scores, freq_info\n        )\n        \n        # 最も確率の高いクラスを選択\n        pred_label = np.argmax(predictions_adjusted[0])\n        \n        # 非常に低い信頼度の場合の特別処理\n        if confidence_scores[0] < 0.25:\n            # 上位3つを確認\n            top_3_indices = np.argsort(predictions_adjusted[0])[-3:]\n            top_3_gestures = [label_encoder.classes_[i] for i in top_3_indices]\n            \n            # BFRB行動を優先\n            for gesture in reversed(top_3_gestures):  # 高い確率から確認\n                if gesture in BFRB_BEHAVIORS:\n                    pred_label = list(label_encoder.classes_).index(gesture)\n                    break\n        \n        # ラベル名に変換\n        pred_gesture = label_encoder.inverse_transform([pred_label])[0]\n        \n        return pred_gesture\n        \n    except Exception as e:\n        print(f\"エラー発生: {str(e)}\")\n        import traceback\n        traceback.print_exc()\n        # エラー時はデフォルト予測\n        return 'Wave hello'"},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":"# CMI評価APIを使用して推論を実行\nimport os\nfrom kaggle_evaluation.cmi_inference_server import CMIInferenceServer\n\n# 推論サーバーの作成\nprint(\"推論サーバーを作成中...\")\ninference_server = CMIInferenceServer(predict)\n\n# スコア情報\nif training_results:\n    best_cv_score = training_results.get('mean_cv_score', 'N/A')\n    bfrb_accuracy = training_results.get('bfrb_accuracy', 'N/A')\nelse:\n    best_cv_score = 'N/A'\n    bfrb_accuracy = 'N/A'\n\nprint(\"=== IMU推論準備完了 (v10.0.0) ===\")\nprint(f\"使用モデル: IMU Only Model with Frequency Features\")\nprint(f\"ベースCV: {best_cv_score}\")\nprint(f\"BFRB検出精度: {bfrb_accuracy}\")\nprint(\"\\n主要な改良点:\")\nprint(\"- 周波数領域特徴量（FFT、スペクトログラム）\")\nprint(\"- 変化点検出による行動パターン認識\")\nprint(\"- 強化された後処理（周波数情報考慮）\")\nprint(\"- より洗練されたBFRB行動の優先順位付け\")\nprint(\"\\n目標: LB 0.70+\")\n\n# 環境に応じて実行\nif os.getenv('KAGGLE_IS_COMPETITION_RERUN'):\n    print(\"\\nKaggle競技環境を検出。推論サーバーを開始します。\")\n    inference_server.serve()\nelse:\n    print(\"\\n推論サーバーを実行します...\")\n    inference_server.run_local_gateway(\n        data_paths=(\n            '/kaggle/input/cmi-detect-behavior-with-sensor-data/test.csv',\n            '/kaggle/input/cmi-detect-behavior-with-sensor-data/test_demographics.csv',\n        )\n    )"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.10.12"}},"nbformat":4,"nbformat_minor":4}