{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"dfd65b3d","cell_type":"markdown","source":"# LANL Earthquake Prediction: Competition Submission Notebook\n\nThis notebook is optimized for generating a clean, competition-ready baseline submission.\n\nDesign decisions:\n\n- Use **full training segments** for stronger baseline quality.\n- Use memory-efficient chunked processing for `train.csv`.\n- Use rich engineered features with practical runtime.\n- Use **time-aware, leakage-conscious** validation.\n- Write final file to `/kaggle/working/submission.csv`.\n","metadata":{}},{"id":"7f8efae5","cell_type":"code","source":"import os\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\n\nfrom scipy.stats import skew, kurtosis, entropy\nfrom scipy.signal import find_peaks\n\nfrom sklearn.metrics import mean_absolute_error\n\nfrom lightgbm import LGBMRegressor, early_stopping, log_evaluation\n\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom joblib import Parallel, delayed\n\nsns.set_theme(style='whitegrid')\n\nSEED = 42\nnp.random.seed(SEED)\n\nSEGMENT_SIZE = 150_000\n\nKAGGLE_TRAIN = Path('/kaggle/input/competitions/LANL-Earthquake-Prediction/train.csv')\nKAGGLE_TEST_DIR = Path('/kaggle/input/competitions/LANL-Earthquake-Prediction/test')\nKAGGLE_SAMPLE = Path('/kaggle/input/competitions/LANL-Earthquake-Prediction/sample_submission.csv')\nKAGGLE_WORKING = Path('/kaggle/working')\n\nif KAGGLE_TRAIN.exists() and KAGGLE_TEST_DIR.exists() and KAGGLE_SAMPLE.exists():\n    TRAIN_PATH = KAGGLE_TRAIN\n    TEST_DIR = KAGGLE_TEST_DIR\n    SAMPLE_SUB_PATH = KAGGLE_SAMPLE\n    WORKING_DIR = KAGGLE_WORKING\n    RUN_MODE = 'kaggle'\nelse:\n    local_data_dir = Path(os.environ.get('LANL_LOCAL_DATA_DIR', ''))\n    assert local_data_dir.exists(), (\n        'Kaggle paths are unavailable. Set LANL_LOCAL_DATA_DIR to a local copy of competition files.'\n    )\n    TRAIN_PATH = local_data_dir / 'train.csv'\n    TEST_DIR = local_data_dir / 'test'\n    SAMPLE_SUB_PATH = local_data_dir / 'sample_submission.csv'\n    WORKING_DIR = Path(os.environ.get('LANL_LOCAL_WORKING_DIR', './_local_working')).resolve()\n    RUN_MODE = 'local-mirror'\n\nfor required_path in [TRAIN_PATH, TEST_DIR, SAMPLE_SUB_PATH]:\n    assert required_path.exists(), f'Missing required competition path: {required_path}'\nWORKING_DIR.mkdir(parents=True, exist_ok=True)\n\nprint(f'Run mode: {RUN_MODE}')\nprint(f'Train path: {TRAIN_PATH}')\nprint(f'Test dir:   {TEST_DIR}')\nprint(f'Output dir: {WORKING_DIR}')","metadata":{"execution":{"iopub.status.busy":"2026-06-12T14:46:31.549900Z","iopub.execute_input":"2026-06-12T14:46:31.550290Z","iopub.status.idle":"2026-06-12T14:46:37.434131Z","shell.execute_reply.started":"2026-06-12T14:46:31.550260Z","shell.execute_reply":"2026-06-12T14:46:37.433188Z"},"trusted":true},"outputs":[],"execution_count":null},{"id":"689af86a","cell_type":"markdown","source":"## 1. Feature Engineering Toolkit\n\nImplemented feature groups:\n\n- Statistical features\n- Trend features\n- Signal/FFT features\n- Rolling-window features\n- Percentile and distribution features\n- Change-rate and peak features\n","metadata":{}},{"id":"2d5c0147","cell_type":"code","source":"def linear_trend(values: np.ndarray) -> float:\n    idx = np.arange(values.size, dtype=np.float64)\n    return float(np.polyfit(idx, values, 1)[0])\n\n\ndef rolling_tail(values: np.ndarray, window: int) -> tuple[float, float, float]:\n    if values.size <= window:\n        return float(np.mean(values)), float(np.std(values)), float(np.quantile(values, 0.5))\n    s = pd.Series(values)\n    r = s.rolling(window=window)\n    return float(r.mean().iloc[-1]), float(r.std().iloc[-1]), float(r.quantile(0.5).iloc[-1])\n\n\ndef compute_features(x: np.ndarray) -> dict[str, float]:\n    x = x.astype(np.float64)\n    abs_x = np.abs(x)\n    diff_x = np.diff(x)\n\n    q01, q05, q10, q25, q50, q75, q90, q95, q99 = np.quantile(\n        x, [0.01, 0.05, 0.10, 0.25, 0.50, 0.75, 0.90, 0.95, 0.99]\n    )\n\n    roll_1k_mean, roll_1k_std, roll_1k_q50 = rolling_tail(x, 1000)\n    roll_10k_mean, roll_10k_std, roll_10k_q50 = rolling_tail(x, 10_000)\n\n    fft_vals = np.fft.rfft(x)\n    power = np.abs(fft_vals) ** 2\n    freqs = np.fft.rfftfreq(x.size, d=1.0)\n    dominant_freq = float(freqs[int(np.argmax(power[1:]) + 1)]) if power.size > 1 else 0.0\n\n    bands = np.array_split(power, 8)\n    band_energy = {f'fft_band_energy_{i}': float(np.mean(b)) for i, b in enumerate(bands)}\n\n    hist = np.histogram(x, bins=64)[0].astype(np.float64)\n    prob = (hist + 1e-12) / np.sum(hist + 1e-12)\n\n    pos_peaks, _ = find_peaks(x, distance=80)\n    neg_peaks, _ = find_peaks(-x, distance=80)\n\n    feats = {\n        # Statistical features\n        'mean': float(np.mean(x)),\n        'std': float(np.std(x)),\n        'min': float(np.min(x)),\n        'max': float(np.max(x)),\n        'median': float(np.median(x)),\n        'q01': float(q01),\n        'q05': float(q05),\n        'q10': float(q10),\n        'q25': float(q25),\n        'q50': float(q50),\n        'q75': float(q75),\n        'q90': float(q90),\n        'q95': float(q95),\n        'q99': float(q99),\n        'abs_mean': float(np.mean(abs_x)),\n        'abs_max': float(np.max(abs_x)),\n\n        # Trend features\n        'trend': linear_trend(x),\n        'abs_trend': linear_trend(abs_x),\n\n        # Signal and spectral features\n        'signal_power': float(np.mean(x ** 2)),\n        'spectral_energy': float(np.sum(power) / power.size),\n        'spectral_power_std': float(np.std(power)),\n        'dominant_freq': dominant_freq,\n\n        # Rolling features\n        'roll_1k_mean_tail': roll_1k_mean,\n        'roll_1k_std_tail': roll_1k_std,\n        'roll_1k_q50_tail': roll_1k_q50,\n        'roll_10k_mean_tail': roll_10k_mean,\n        'roll_10k_std_tail': roll_10k_std,\n        'roll_10k_q50_tail': roll_10k_q50,\n\n        # Additional distribution and change-rate features\n        'mean_abs_change': float(np.mean(np.abs(diff_x))),\n        'std_change': float(np.std(diff_x)),\n        'nonzero_change_ratio': float(np.mean(diff_x != 0)),\n        'signal_entropy': float(entropy(prob)),\n        'skewness': float(skew(x)),\n        'kurtosis': float(kurtosis(x)),\n        'num_pos_peaks': float(len(pos_peaks)),\n        'num_neg_peaks': float(len(neg_peaks)),\n        'peak_diff': float(len(pos_peaks) - len(neg_peaks)),\n    }\n\n    feats.update(band_energy)\n    return feats","metadata":{"execution":{"iopub.status.busy":"2026-06-12T14:46:37.435735Z","iopub.execute_input":"2026-06-12T14:46:37.436004Z","iopub.status.idle":"2026-06-12T14:46:37.455896Z","shell.execute_reply.started":"2026-06-12T14:46:37.435978Z","shell.execute_reply":"2026-06-12T14:46:37.454877Z"},"trusted":true},"outputs":[],"execution_count":null},{"id":"5402a29a","cell_type":"markdown","source":"## 2. Efficient Train Feature Pipeline\n\nWe read `train.csv` in fixed-size chunks (`150,000`) so each chunk maps directly to one training segment.\n\nWe also infer event IDs from `time_to_failure` resets for leakage-aware validation.\n","metadata":{}},{"id":"c67651da","cell_type":"code","source":"def build_train_features(train_path: Path, segment_size: int) -> pd.DataFrame:\n    rows = []\n    prev_ttf = None\n    event_id = 0\n\n    reader = pd.read_csv(\n        train_path,\n        dtype={'acoustic_data': np.int16, 'time_to_failure': np.float32},\n        chunksize=segment_size,\n    )\n\n    for segment_idx, chunk in enumerate(reader):\n        if len(chunk) != segment_size:\n            continue\n\n        x = chunk['acoustic_data'].to_numpy(dtype=np.float64)\n        ttf = chunk['time_to_failure'].to_numpy(dtype=np.float64)\n\n        boundary_reset = int(prev_ttf is not None and ttf[0] > prev_ttf + 1e-8)\n        intra_reset = int(np.any(np.diff(ttf) > 1e-8))\n        event_id += boundary_reset + intra_reset\n        prev_ttf = float(ttf[-1])\n\n        feats = compute_features(x)\n        feats.update({\n            'segment_idx': segment_idx,\n            'event_id': event_id,\n            'target': float(ttf[-1]),\n            'boundary_reset': boundary_reset,\n            'intra_reset': intra_reset,\n        })\n        rows.append(feats)\n\n    return pd.DataFrame(rows)\n\n\ntrain_features = build_train_features(TRAIN_PATH, SEGMENT_SIZE)\nprint('Train feature shape:', train_features.shape)\n\ntrain_feature_path = WORKING_DIR / 'train_features.parquet'\ntrain_features.to_parquet(train_feature_path, index=False)\nprint('Saved:', train_feature_path)","metadata":{"execution":{"iopub.status.busy":"2026-06-12T14:46:37.457314Z","iopub.execute_input":"2026-06-12T14:46:37.457688Z"},"trusted":true},"outputs":[],"execution_count":null},{"id":"b134fd4b","cell_type":"markdown","source":"## 3. Time-Aware Validation\n\nRandom split is avoided.\n\nValidation logic:\n\n- Keep segments ordered by time.\n- Create contiguous validation blocks from later event IDs.\n- Purge one event block before validation to reduce leakage bleed.\n","metadata":{}},{"id":"ecf1ae21","cell_type":"code","source":"def make_purged_group_folds(df: pd.DataFrame, n_splits: int = 4, purge_groups: int = 1):\n    groups = df['event_id'].to_numpy()\n    unique_groups = pd.unique(groups)\n    blocks = [b for b in np.array_split(unique_groups, n_splits) if len(b) > 0]\n\n    folds = []\n    for fold_id in range(1, len(blocks)):\n        valid_groups = blocks[fold_id]\n        valid_start = int(np.where(unique_groups == valid_groups[0])[0][0])\n        train_end = max(0, valid_start - purge_groups)\n        train_groups = unique_groups[:train_end]\n        if len(train_groups) == 0:\n            continue\n\n        train_idx = np.where(np.isin(groups, train_groups))[0]\n        valid_idx = np.where(np.isin(groups, valid_groups))[0]\n        if len(train_idx) == 0 or len(valid_idx) == 0:\n            continue\n\n        folds.append((fold_id, train_idx, valid_idx))\n\n    if not folds:\n        raise RuntimeError('Unable to create valid purged folds.')\n    return folds\n\n\nfeature_cols = [\n    c for c in train_features.columns\n    if c not in {'target', 'segment_idx', 'event_id', 'boundary_reset', 'intra_reset'}\n]\n\nfolds = make_purged_group_folds(train_features, n_splits=5, purge_groups=1)\nprint('Number of folds:', len(folds))\nfor fold_id, tr_idx, va_idx in folds:\n    print(f'Fold {fold_id}: train={len(tr_idx)}, valid={len(va_idx)}')","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"772bddf7","cell_type":"markdown","source":"## 4. LightGBM Baseline + Alternative Spot Check\n\nPrimary model: **LightGBM Regressor** (fast and strong on tabular engineered features).\n\nAlternative spot-checks:\n\n- XGBoost (if available)\n- CatBoost (if available)\n\nThese alternatives are evaluated on the first fold only to keep runtime practical.\n","metadata":{}},{"id":"94057743","cell_type":"code","source":"cv_records = []\nfold_models = []\n\nfor fold_id, tr_idx, va_idx in folds:\n    X_tr = train_features.iloc[tr_idx][feature_cols]\n    y_tr = train_features.iloc[tr_idx]['target']\n    X_va = train_features.iloc[va_idx][feature_cols]\n    y_va = train_features.iloc[va_idx]['target']\n\n    model = LGBMRegressor(\n        objective='mae',\n        n_estimators=7000,\n        learning_rate=0.02,\n        num_leaves=96,\n        subsample=0.85,\n        colsample_bytree=0.8,\n        random_state=SEED,\n        verbosity=-1,\n    )\n\n    model.fit(\n        X_tr,\n        y_tr,\n        eval_set=[(X_va, y_va)],\n        eval_metric='l1',\n        callbacks=[early_stopping(stopping_rounds=250, verbose=False), log_evaluation(period=0)],\n    )\n\n    pred = model.predict(X_va)\n    mae = mean_absolute_error(y_va, pred)\n\n    cv_records.append({'model': 'lightgbm', 'fold': fold_id, 'mae': mae, 'n_train': len(tr_idx), 'n_valid': len(va_idx)})\n    fold_models.append(model)\n\ncv_df = pd.DataFrame(cv_records)\ncv_df","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"5e94b33e","cell_type":"code","source":"# Optional quick alternatives on first fold only\nalt_results = []\nfold_id, tr_idx, va_idx = folds[0]\nX_tr = train_features.iloc[tr_idx][feature_cols]\ny_tr = train_features.iloc[tr_idx]['target']\nX_va = train_features.iloc[va_idx][feature_cols]\ny_va = train_features.iloc[va_idx]['target']\n\n# XGBoost\ntry:\n    from xgboost import XGBRegressor\n    xgb = XGBRegressor(\n        objective='reg:absoluteerror',\n        n_estimators=1500,\n        learning_rate=0.03,\n        max_depth=8,\n        subsample=0.8,\n        colsample_bytree=0.8,\n        random_state=SEED,\n        tree_method='hist',\n    )\n    xgb.fit(X_tr, y_tr)\n    alt_results.append({'model': 'xgboost_fold1', 'mae': mean_absolute_error(y_va, xgb.predict(X_va))})\nexcept Exception as exc:\n    alt_results.append({'model': 'xgboost_fold1', 'mae': np.nan, 'note': f'skipped: {exc}'})\n\n# CatBoost\ntry:\n    from catboost import CatBoostRegressor\n    cat = CatBoostRegressor(\n        loss_function='MAE',\n        eval_metric='MAE',\n        iterations=1500,\n        learning_rate=0.03,\n        depth=8,\n        random_seed=SEED,\n        verbose=False,\n    )\n    cat.fit(X_tr, y_tr)\n    alt_results.append({'model': 'catboost_fold1', 'mae': mean_absolute_error(y_va, cat.predict(X_va))})\nexcept Exception as exc:\n    alt_results.append({'model': 'catboost_fold1', 'mae': np.nan, 'note': f'skipped: {exc}'})\n\npd.DataFrame(alt_results)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"1a5b8647","cell_type":"code","source":"print('LightGBM CV MAE mean:', round(cv_df['mae'].mean(), 6))\nprint('LightGBM CV MAE std :', round(cv_df['mae'].std(), 6))","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"70e327bd","cell_type":"markdown","source":"## 5. Feature Importance\n\nUsing the final fold model as a practical proxy for feature importance ranking.\n","metadata":{}},{"id":"27f72387","cell_type":"code","source":"importance_df = pd.DataFrame({\n    'feature': feature_cols,\n    'importance': fold_models[-1].feature_importances_,\n}).sort_values('importance', ascending=False)\n\nplt.figure(figsize=(10, 8))\nsns.barplot(data=importance_df.head(25), y='feature', x='importance', orient='h')\nplt.title('Top-25 Feature Importances (LightGBM)')\nplt.tight_layout()\nplt.show()\n\nimportance_df.head(10)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"c9ddee6f","cell_type":"markdown","source":"## 6. Test Inference Pipeline\n\nWe process all test segments with the same feature extractor.\n\nTo improve runtime, feature extraction is parallelized with joblib.\n","metadata":{}},{"id":"48ca5014","cell_type":"code","source":"def process_test_segment(file_path: Path) -> dict[str, float]:\n    seg = pd.read_csv(file_path, dtype={'acoustic_data': np.int16})\n    x = seg['acoustic_data'].to_numpy(dtype=np.float64)\n    feats = compute_features(x)\n    feats['seg_id'] = file_path.stem\n    return feats\n\n\ntest_files = sorted(TEST_DIR.glob('*.csv'))\nprint('Test files:', len(test_files))\n\ntest_rows = Parallel(n_jobs=min(8, os.cpu_count() or 2), backend='loky')(\n    delayed(process_test_segment)(fp) for fp in test_files\n)\n\ntest_features = pd.DataFrame(test_rows).sort_values('seg_id').reset_index(drop=True)\n\nX_test = test_features[feature_cols]\n\ntest_feature_path = WORKING_DIR / 'test_features.parquet'\ntest_features.to_parquet(test_feature_path, index=False)\nprint('Saved:', test_feature_path)\nprint('Test feature shape:', test_features.shape)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"5649de2b","cell_type":"markdown","source":"## 7. Final Training and Submission Generation\n\nWe train one final LightGBM model on all training segments and produce Kaggle submission format.\n","metadata":{}},{"id":"15ffa88c","cell_type":"code","source":"final_model = LGBMRegressor(\n    objective='mae',\n    n_estimators=3000,\n    learning_rate=0.02,\n    num_leaves=96,\n    subsample=0.9,\n    colsample_bytree=0.8,\n    random_state=SEED,\n    verbosity=-1,\n)\n\nfinal_model.fit(train_features[feature_cols], train_features['target'])\nfinal_pred = final_model.predict(X_test)\n\nsubmission = pd.DataFrame({\n    'seg_id': test_features['seg_id'],\n    'time_to_failure': np.clip(final_pred, a_min=0.0, a_max=None),\n})\n\nsubmission_path = WORKING_DIR / 'submission.csv'\nsubmission.to_csv(submission_path, index=False)\n\nprint('Submission saved to:', submission_path)\nsubmission.head()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"id":"82dce092","cell_type":"code","source":"assert (WORKING_DIR / 'submission.csv').exists(), 'submission.csv was not created.'\nsample = pd.read_csv(SAMPLE_SUB_PATH)\nsub = pd.read_csv(WORKING_DIR / 'submission.csv')\n\nassert list(sub.columns) == ['seg_id', 'time_to_failure'], 'Submission columns mismatch.'\nassert len(sub) == len(sample), 'Submission row count mismatch.'\nprint('Submission format validation passed.')","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}