{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.11","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":101849,"databundleVersionId":12846694,"sourceType":"competition"},{"sourceId":247163208,"sourceType":"kernelVersion"}],"dockerImageVersionId":31040,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Ariel '25 Baseline: FGS1->Features->LGBM\n\n* Wow - 263GB of data!\n* Just a single transit has 135,000 32x32 FGS1 frames!\n* Hard to figure out where to start....\n\nThe trick to machine learning seems to be making a new problem look like one you already know. So - let's try some feature engineering to **turn this all into tabular data and feed it into LGBM.**\n\n*Friendly Reminder:* If re-using large parts of this work in a public notebook - **please credit where you found the code**\n\n## Approach\n\nWe focus on **just the FGS1 data** - we do ADC correction but **ignore other calibration (dead-pixels, etc.)** for now.\n\nSince we have 283 spectral targets, we train *283 separate LGBM models per CV fold.* That's **1415 LGBM models total!**\n\nMultiple signals create multiple rows of tabular data for both train and test.  Planets with multiple signals result in multiple predictions (which are averaged together).\n\n# Engineered Features\n\nI spent some quality time with Claude to come up with a pipeline extracts **175 engineered features** from the FGS1 data:\n\n| Category | Count | Purpose |\n|----------|-------|---------|\n| **Global Flux Statistics** | 19 | Overall brightness characteristics & fundamental transit depth measurement |\n| **Rolling Statistics** | 54 | Multi-timescale transit detection across 6 window sizes (5 sec to 8.3 min) |\n| **Transit Detection** | 20 | Specialized transit identification, timing, duration & shape analysis |\n| **Frequency Domain** | 16 | Signal/noise separation & periodic pattern detection via FFT & autocorrelation |\n| **Spatial Features** | 30 | Star position/shape changes during transit across 5 key timepoints |\n| **Gradient & Change** | 36 | Rate of brightness change & 4-segment temporal analysis |\n\nThese features are combined with the star info features **(Rs, Ms, Ts, Mp, e, P, sma, i)** for training the LGBM models.\n\n## LB Score Estimate + Sigma Optimization\n\nFirst, we use **residual errors** from our train process to try predict those pesky **uncertainties (sigma) columns**.\n\nThen we use the competition metric to **estimate our LB Score**.\n\nFinally, we run a **per-wavelength optimization loop to tune sigma values** - and then generate an updated LB estimate.\n* Note: Sigma optimization hasn't produced the gains it has for other notebooks.  Unclear if this is a difference in implementation - or the residual error approach is already optimal.\n\n## Timing\n\n- **Feature extraction**: ~2 hours for training data\n- **Model training**: Few minutes for 1415 models  \n- **Scoring**: ~4 hours total (includes feature extraction on train+test data)","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport lightgbm as lgb\nfrom sklearn.model_selection import KFold\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import mean_squared_error, r2_score\nfrom tqdm import tqdm\nimport time\nimport os\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport glob\n\nfrom scipy import signal, stats\nfrom scipy.stats import skew, kurtosis\nfrom scipy.optimize import minimize_scalar\nfrom scipy.stats import norm\n\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimport sys","metadata":{"trusted":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2025-07-07T02:17:21.136089Z","iopub.execute_input":"2025-07-07T02:17:21.136881Z","iopub.status.idle":"2025-07-07T02:17:27.890396Z","shell.execute_reply.started":"2025-07-07T02:17:21.136844Z","shell.execute_reply":"2025-07-07T02:17:27.889515Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Load Data","metadata":{}},{"cell_type":"code","source":"base_path = '/kaggle/input/ariel-data-challenge-2025/'\n\ntrain_df = pd.read_csv(f\"{base_path}/train.csv\")\nprint(\"Train\")\ndisplay(train_df.head(2))\n\ntrain_star_info_df = pd.read_csv(f\"{base_path}/train_star_info.csv\")\nprint(\"Train Star Info\")\ndisplay(train_star_info_df.head(2))\n\ntest_star_info_df = pd.read_csv(f\"{base_path}/test_star_info.csv\")\nprint(\"Test\")\ndisplay(test_star_info_df)","metadata":{"trusted":true,"_kg_hide-input":false,"execution":{"iopub.status.busy":"2025-07-07T02:17:27.892199Z","iopub.execute_input":"2025-07-07T02:17:27.892821Z","iopub.status.idle":"2025-07-07T02:17:28.113376Z","shell.execute_reply.started":"2025-07-07T02:17:27.892794Z","shell.execute_reply":"2025-07-07T02:17:28.112757Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# for doing quick runs (not loading all data)\nquick_run = False\nquick_run_size = 200\n\n# disable quick_run if submitting!\nif len(test_star_info_df) > 1:\n    quick_run = False\n\nif quick_run:\n    train_star_info_df=train_star_info_df.head(quick_run_size)\n    print(f\"Only loading {quick_run_size} records - this run is not a full training!\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T02:56:11.86473Z","iopub.execute_input":"2025-07-07T02:56:11.865015Z","iopub.status.idle":"2025-07-07T02:56:11.870928Z","shell.execute_reply.started":"2025-07-07T02:56:11.864997Z","shell.execute_reply":"2025-07-07T02:56:11.869937Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Load Transit Data for Demo\n\nLoads a single planet's transit data for feature engineering demonstration.","metadata":{}},{"cell_type":"code","source":"%%time\n# Load sample data for demonstrations\nplanet_id = int(train_star_info_df.iloc[0]['planet_id'])\nfile_path = f\"/kaggle/input/ariel-data-challenge-2025/train/{planet_id}/FGS1_signal_0.parquet\"\ndf = pd.read_parquet(file_path)\nsignal_data = df.values.reshape(135000, 32, 32)\n\nprint(f\"Demo data loaded for Planet ID: {planet_id}\")\nprint(f\"Signal data shape: {signal_data.shape}\")\n\nprint(\"=\"*60)","metadata":{"trusted":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2025-07-07T02:17:28.120166Z","iopub.execute_input":"2025-07-07T02:17:28.120944Z","iopub.status.idle":"2025-07-07T02:17:29.384145Z","shell.execute_reply.started":"2025-07-07T02:17:28.120916Z","shell.execute_reply":"2025-07-07T02:17:29.383345Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Apply ADC Correction\n* Converts our ADC values into physical units\n* Also - calculate total flux (need ADC correction first)","metadata":{}},{"cell_type":"code","source":"def apply_adc_correction(signal_array, instrument, adc_info_path=\"/kaggle/input/ariel-data-challenge-2025/adc_info.csv\"):\n    \"\"\"\n    Applies ADC correction to raw signal data using known gain and offset.\n\n    Args:\n        signal_array (np.ndarray): Raw uint16 signal data.\n        instrument (str): 'FGS1' or 'AIRS-CH0'.\n        adc_info_path (str): Path to adc_info.csv file.\n\n    Returns:\n        np.ndarray: Calibrated float signal array.\n    \"\"\"\n    # Load ADC values\n    adc_df = pd.read_csv(adc_info_path)\n    \n    gain_col = f\"{instrument}_adc_gain\"\n    offset_col = f\"{instrument}_adc_offset\"\n\n    if gain_col not in adc_df.columns or offset_col not in adc_df.columns:\n        raise ValueError(f\"Missing columns: {gain_col} or {offset_col} in adc_info.csv\")\n    \n    gain = adc_df.at[0, gain_col]\n    offset = adc_df.at[0, offset_col]\n\n    # Apply gain and offset\n    calibrated_signal = signal_array.astype(np.float32) * gain + offset\n    return calibrated_signal\n    \nsignal_data = apply_adc_correction(signal_data, instrument='FGS1')\ntotal_flux = np.sum(signal_data, axis=(1, 2))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T02:17:29.386378Z","iopub.execute_input":"2025-07-07T02:17:29.386629Z","iopub.status.idle":"2025-07-07T02:17:29.754865Z","shell.execute_reply.started":"2025-07-07T02:17:29.386609Z","shell.execute_reply":"2025-07-07T02:17:29.753982Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Visualize a Transit\n\n* The graph at the bottom reflects the brightness level changes across **135,000 FGS1 images**\n* Our feature engineering will try to capture details of this curve in tabular data\n* Note the \"hot\" pixels in the brightest image (a future version of this notebook will use the calibration data to fix this)","metadata":{}},{"cell_type":"code","source":"%%time\ndef visualize_transit_data(signal_data, total_flux, planet_id):    \n    # Key frame indices (quarter-points + first/last)\n    base_frames = [0, len(signal_data)//4, len(signal_data)//2, 3*len(signal_data)//4, -1]\n\n    # Find brightest and darkest frame indices\n    brightest_idx = np.argmax(total_flux)\n    darkest_idx = np.argmin(total_flux)\n\n    # Add to the frame list, avoiding duplicates\n    extra_frames = []\n    for idx in [brightest_idx, darkest_idx]:\n        if idx not in base_frames and (idx != -1 and idx != len(signal_data)-1):\n            extra_frames.append(idx)\n    frames = base_frames + extra_frames\n\n    # Set up clean styling\n    plt.style.use('default')\n    fig = plt.figure(figsize=(18, 10))\n    \n    # Top row: telescope frames\n    for i, idx in enumerate(frames):\n        ax = plt.subplot(2, len(frames), i + 1)\n        frame = signal_data[idx]\n        \n        # Contrast scaling\n        vmin, vmax = np.percentile(frame, [2, 98])\n        im = ax.imshow(frame, cmap='hot', aspect='equal', vmin=vmin, vmax=vmax)\n        \n        # Annotate\n        time_min = idx * 0.1 / 60\n        if idx == brightest_idx:\n            title = f'Brightest\\nT={time_min:.1f} min'\n        elif idx == darkest_idx:\n            title = f'Darkest\\nT={time_min:.1f} min'\n        else:\n            title = f'T={time_min:.1f} min'\n        ax.set_title(title, fontsize=11)\n        ax.set_xticks([])\n        ax.set_yticks([])\n    \n    # Bottom: Light curve\n    ax_light = plt.subplot(2, 1, 2)\n    time_hours = np.arange(len(total_flux)) * 0.1 / 3600\n    sample = slice(None, None, max(1, len(total_flux)//2000))\n\n    ax_light.plot(time_hours[sample], total_flux[sample], \n                 color='lightsteelblue', alpha=0.4, linewidth=0.5, label='Raw flux')\n\n    window = 500\n    moving_avg = pd.Series(total_flux).rolling(window, center=True).mean()\n    ax_light.plot(time_hours[sample], moving_avg.iloc[sample], \n                 color='darkblue', linewidth=3, label=f'{window}-frame average')\n\n    colors = ['red', 'orange', 'green', 'purple', 'brown', 'lime', 'black']\n    for i, idx in enumerate(frames):\n        time_point = idx * 0.1 / 3600\n        ax_light.axvline(time_point, color=colors[i % len(colors)], alpha=0.8, linewidth=1.5, linestyle='--')\n\n    ax_light.set_xlabel('Time (hours)', fontsize=12)\n    ax_light.set_ylabel('Total Flux (counts)', fontsize=12)\n    ax_light.set_title(f'Transit Light Curve - Planet {planet_id}', fontsize=14, pad=15)\n    ax_light.grid(True, alpha=0.3, linestyle='-', linewidth=0.5)\n    ax_light.legend(loc='upper right', framealpha=0.9)\n\n    smooth_flux = moving_avg.dropna()\n    flux_min, flux_max = smooth_flux.min(), smooth_flux.max()\n    flux_range = flux_max - flux_min\n    margin = flux_range * 0.1\n    ax_light.set_ylim(flux_min - margin, flux_max + margin)\n\n    plt.tight_layout(pad=2.0)\n    plt.show()\n\n    # Enhanced summary\n    duration = time_hours[-1]\n    transit_depth = np.mean(total_flux) - np.min(total_flux)\n\n    print(f\"🌟 Planet {planet_id} Transit Observation\")\n    print(f\"   Duration: {duration:.2f} hours ({len(signal_data):,} frames)\")\n    print(f\"   Brightness: {np.min(total_flux):,.0f} → {np.max(total_flux):,.0f} counts\")\n    print(f\"   Brightest Frame: {brightest_idx} | Darkest Frame: {darkest_idx}\")\n    print(f\"   Transit depth: {transit_depth:,.0f} counts ({transit_depth/np.mean(total_flux)*100:.3f}%)\")\n    \nvisualize_transit_data(signal_data, total_flux, planet_id)","metadata":{"trusted":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2025-07-07T02:17:29.755722Z","iopub.execute_input":"2025-07-07T02:17:29.755937Z","iopub.status.idle":"2025-07-07T02:17:30.703825Z","shell.execute_reply.started":"2025-07-07T02:17:29.755921Z","shell.execute_reply":"2025-07-07T02:17:30.70291Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Global Flux Features\n\nExtracts **19 statistical features** from total stellar brightness across the entire observation.\n\n- **Basic stats**: Mean, std, min/max, range, skewness, kurtosis, coefficient of variation\n- **Percentiles**: 9 robust statistical measures (1st to 99th percentiles)  \n- **Transit depth**: Direct measurement of brightness drop","metadata":{}},{"cell_type":"code","source":"%%time\ndef extract_global_flux_features(total_flux):\n    features = {}\n    \n    # Basic statistics\n    features['global_flux_mean'] = np.mean(total_flux)\n    features['global_flux_std'] = np.std(total_flux)\n    features['global_flux_min'] = np.min(total_flux)\n    features['global_flux_max'] = np.max(total_flux)\n    features['global_flux_range'] = features['global_flux_max'] - features['global_flux_min']\n    features['global_flux_skew'] = skew(total_flux)\n    features['global_flux_kurtosis'] = kurtosis(total_flux)\n    features['global_flux_cv'] = features['global_flux_std'] / features['global_flux_mean']\n    \n    # Percentiles\n    for p in [1, 5, 10, 25, 50, 75, 90, 95, 99]:\n        features[f'global_flux_p{p}'] = np.percentile(total_flux, p)\n    \n    # Transit depth features\n    features['global_flux_depth'] = features['global_flux_mean'] - features['global_flux_min']\n    features['global_flux_depth_ratio'] = features['global_flux_depth'] / features['global_flux_mean']\n    \n    return features\n\n# Demo: Global Flux Features\nglobal_features = extract_global_flux_features(total_flux)\nprint(f\"Generated {len(global_features)} global flux features:\")\n\n# Print ALL features organized by category\nprint(\"\\nBasic Statistics:\")\nbasic_stats = ['global_flux_mean', 'global_flux_std', 'global_flux_min', 'global_flux_max', \n               'global_flux_range', 'global_flux_skew', 'global_flux_kurtosis', 'global_flux_cv']\nfor feature in basic_stats:\n    print(f\"  {feature}: {global_features[feature]:.4f}\")\n\nprint(\"\\nPercentiles:\")\nfor p in [1, 5, 10, 25, 50, 75, 90, 95, 99]:\n    feature = f'global_flux_p{p}'\n    print(f\"  {feature}: {global_features[feature]:.4f}\")\n\nprint(\"\\nTransit Depth:\")\ntransit_features = ['global_flux_depth', 'global_flux_depth_ratio']\nfor feature in transit_features:\n    print(f\"  {feature}: {global_features[feature]:.6f}\")\n\nprint(f\"\\nKey transit indicator - depth ratio: {global_features['global_flux_depth_ratio']:.6f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T02:17:30.704744Z","iopub.execute_input":"2025-07-07T02:17:30.704997Z","iopub.status.idle":"2025-07-07T02:17:30.744417Z","shell.execute_reply.started":"2025-07-07T02:17:30.704977Z","shell.execute_reply":"2025-07-07T02:17:30.743764Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Rolling Statistics Features\n\nExtracts **54 multi-timescale time series features** by analyzing brightness patterns using sliding windows of different sizes.\n\n**Time Series Processing**:\n- **Sliding windows**: 50→5000 frames (5 seconds to 8.3 minutes) slide across the 3.75-hour brightness time series\n- **Per window**: Calculates local mean, std, min, max at each time point in the sequence\n- **Temporal analysis**: Finds extremes, ranges, and variability patterns within each rolling window size\n\n**What it captures across time**:\n- **50-frame windows**: Instrument noise and rapid variations\n- **1000-frame windows**: Transit ingress/egress segments  \n- **5000-frame windows**: Full observation trends and baseline drift","metadata":{}},{"cell_type":"code","source":"%%time\ndef extract_rolling_statistics_features(total_flux):\n    features = {}\n    window_sizes = [50, 100, 500, 1000, 2000, 5000]\n    \n    for window in window_sizes:\n        if window < len(total_flux):\n            # Calculate rolling statistics\n            rolling_mean = pd.Series(total_flux).rolling(window=window, center=True).mean().dropna()\n            rolling_std = pd.Series(total_flux).rolling(window=window, center=True).std().dropna()\n            rolling_min = pd.Series(total_flux).rolling(window=window, center=True).min().dropna()\n            rolling_max = pd.Series(total_flux).rolling(window=window, center=True).max().dropna()\n            \n            if len(rolling_mean) > 0:\n                # Statistics of rolling statistics\n                features[f'rolling{window}_mean_min'] = rolling_mean.min()\n                features[f'rolling{window}_mean_max'] = rolling_mean.max()\n                features[f'rolling{window}_mean_std'] = rolling_mean.std()\n                features[f'rolling{window}_mean_range'] = rolling_mean.max() - rolling_mean.min()\n                \n                features[f'rolling{window}_std_mean'] = rolling_std.mean()\n                features[f'rolling{window}_std_max'] = rolling_std.max()\n                \n                # Extreme values\n                features[f'rolling{window}_deepest_dip'] = rolling_min.min()\n                features[f'rolling{window}_highest_peak'] = rolling_max.max()\n                features[f'rolling{window}_volatility'] = rolling_std.std()\n    \n    return features\n\n# Demo: Rolling Statistics Features\nrolling_features = extract_rolling_statistics_features(total_flux)\nprint(f\"Generated {len(rolling_features)} rolling statistics features:\")\n\n# Print ALL features organized by window size\nfor window in [50, 100, 500, 1000, 2000, 5000]:\n    window_features = {k: v for k, v in rolling_features.items() if f'rolling{window}_' in k}\n    if window_features:\n        print(f\"\\nWindow {window} features ({len(window_features)} total):\")\n        for feature_name, value in window_features.items():\n            print(f\"  {feature_name}: {value:.4f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T02:17:30.745221Z","iopub.execute_input":"2025-07-07T02:17:30.745499Z","iopub.status.idle":"2025-07-07T02:17:30.919826Z","shell.execute_reply.started":"2025-07-07T02:17:30.745475Z","shell.execute_reply":"2025-07-07T02:17:30.918738Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Transit Detection Features\n\nExtracts **20 transit-specific features** by removing baseline trends and analyzing dip patterns.\n\n**Data Processing**: **Detrends** (removes slow baseline drift using 5000-frame rolling median) to isolate transit signals, detects dips below 1σ threshold, measures durations and timing.\n- **5 detrended stats**: Min, std, skewness, negative excursions (2σ, 3σ)\n- **4 duration metrics**: Longest dip, number of periods, total time, average duration  \n- **3 timing features**: Deepest point fraction, first-half flag, middle-third flag\n- **4 transit shape**: Transit depth mean, std, asymmetry, flatness\n- **4 observation structure**: First quarter mean, last quarter mean, middle half mean, middle vs edges\n\n**Fallback values are assigned if** no brightness dips below 1σ threshold are detected.","metadata":{}},{"cell_type":"code","source":"%%time\ndef extract_transit_detection_features(total_flux):\n    features = {}\n    \n    # Detrend the light curve\n    baseline = pd.Series(total_flux).rolling(window=5000, center=True).median().fillna(method='bfill').fillna(method='ffill')\n    detrended = total_flux - baseline\n    \n    # Detrended statistics\n    features['detrended_min'] = np.min(detrended)\n    features['detrended_std'] = np.std(detrended)\n    features['detrended_skew'] = skew(detrended)\n    features['detrended_neg_excursions'] = np.sum(detrended < -2 * np.std(detrended))\n    features['detrended_deep_excursions'] = np.sum(detrended < -3 * np.std(detrended))\n    \n    # Transit period detection\n    threshold = np.mean(total_flux) - 1.0 * np.std(total_flux)\n    below_threshold = total_flux < threshold\n    \n    if np.any(below_threshold):\n        # Duration analysis\n        diff = np.diff(np.concatenate(([False], below_threshold, [False])).astype(int))\n        starts = np.where(diff == 1)[0]\n        ends = np.where(diff == -1)[0]\n        durations = ends - starts\n        \n        features['longest_dip_duration'] = np.max(durations) if len(durations) > 0 else 0\n        features['num_dip_periods'] = len(durations)\n        features['total_dip_time'] = np.sum(durations)\n        features['avg_dip_duration'] = np.mean(durations) if len(durations) > 0 else 0\n        \n        # Timing analysis\n        deepest_idx = np.argmin(total_flux)\n        features['deepest_time_fraction'] = deepest_idx / len(total_flux)\n        features['deepest_in_first_half'] = float(deepest_idx < len(total_flux) / 2)\n        features['deepest_in_middle_third'] = float(len(total_flux) / 3 < deepest_idx < 2 * len(total_flux) / 3)\n        \n        # Transit shape\n        transit_flux = total_flux[below_threshold]\n        features['transit_depth_mean'] = np.mean(transit_flux)\n        features['transit_depth_std'] = np.std(transit_flux)\n        features['transit_assymetry'] = skew(transit_flux)\n        features['transit_flatness'] = kurtosis(transit_flux)\n    else:\n        # No transit detected\n        features.update({\n            'longest_dip_duration': 0, 'num_dip_periods': 0, 'total_dip_time': 0, 'avg_dip_duration': 0,\n            'deepest_time_fraction': 0.5, 'deepest_in_first_half': 0, 'deepest_in_middle_third': 0,\n            'transit_depth_mean': np.mean(total_flux), 'transit_depth_std': 0, \n            'transit_assymetry': 0, 'transit_flatness': 0\n        })\n    \n    # Observation structure\n    first_quarter = total_flux[:len(total_flux)//4]\n    last_quarter = total_flux[-len(total_flux)//4:]\n    middle_half = total_flux[len(total_flux)//4:-len(total_flux)//4]\n    \n    features['first_quarter_mean'] = np.mean(first_quarter)\n    features['last_quarter_mean'] = np.mean(last_quarter)\n    features['middle_half_mean'] = np.mean(middle_half)\n    features['middle_vs_edges'] = features['middle_half_mean'] - (features['first_quarter_mean'] + features['last_quarter_mean']) / 2\n    \n    return features\n\n# Demo: Transit Detection Features\ntransit_features = extract_transit_detection_features(total_flux)\nprint(f\"Generated {len(transit_features)} transit detection features:\")\n\n# Print ALL features organized by category\nprint(\"\\nDetrended Statistics:\")\ndetrended_features = ['detrended_min', 'detrended_std', 'detrended_skew', \n                     'detrended_neg_excursions', 'detrended_deep_excursions']\nfor feature in detrended_features:\n    print(f\"  {feature}: {transit_features[feature]:.4f}\")\n\nprint(\"\\nDuration Metrics:\")\nduration_features = ['longest_dip_duration', 'num_dip_periods', 'total_dip_time', 'avg_dip_duration']\nfor feature in duration_features:\n    print(f\"  {feature}: {transit_features[feature]:.4f}\")\n\nprint(\"\\nTiming Features:\")\ntiming_features = ['deepest_time_fraction', 'deepest_in_first_half', 'deepest_in_middle_third']\nfor feature in timing_features:\n    print(f\"  {feature}: {transit_features[feature]:.4f}\")\n\nprint(\"\\nTransit Shape:\")\nshape_features = ['transit_depth_mean', 'transit_depth_std', 'transit_assymetry', 'transit_flatness']\nfor feature in shape_features:\n    print(f\"  {feature}: {transit_features[feature]:.4f}\")\n\nprint(\"\\nObservation Structure:\")\nstructure_features = ['first_quarter_mean', 'last_quarter_mean', 'middle_half_mean', 'middle_vs_edges']\nfor feature in structure_features:\n    print(f\"  {feature}: {transit_features[feature]:.4f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T02:17:30.920861Z","iopub.execute_input":"2025-07-07T02:17:30.921169Z","iopub.status.idle":"2025-07-07T02:17:31.114374Z","shell.execute_reply.started":"2025-07-07T02:17:30.921144Z","shell.execute_reply":"2025-07-07T02:17:31.113401Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Frequency Domain Features\n\nExtracts **16 frequency-based features** using FFT analysis and autocorrelation to separate signals from noise.\n\n- **FFT transformation**: Converts time-series brightness to frequency domain power spectrum\n- **Low-frequency isolation**: Separates transit-like signals (low freq) from instrument noise (high freq)\n- **Autocorrelation analysis**: Detects repeating patterns and periodic structures\n- **Spectral characterization**: Measures where most power is concentrated and how spread out\n\n**Key outputs**: `fft_low_freq_ratio` (transit signal vs noise), `spectral_centroid` (dominant frequency), `autocorr_num_peaks` (periodic pattern detection), power spectrum statistics.\n\nNote: These features are relatively **slow to generate - 3.5s per record** (GPU may help)","metadata":{}},{"cell_type":"code","source":"%%time\ndef extract_frequency_features(total_flux):\n    features = {}\n    \n    # FFT analysis\n    fft_flux = np.fft.fft(total_flux - np.mean(total_flux))\n    fft_power = np.abs(fft_flux)\n    fft_freqs = np.fft.fftfreq(len(total_flux))\n    \n    # Power spectrum\n    features['fft_peak_power'] = np.max(fft_power[1:len(fft_power)//2])\n    features['fft_total_power'] = np.sum(fft_power[1:len(fft_power)//2])\n    features['fft_mean_power'] = np.mean(fft_power[1:len(fft_power)//2])\n    features['fft_std_power'] = np.std(fft_power[1:len(fft_power)//2])\n    \n    # Low frequency analysis\n    low_freq_mask = np.abs(fft_freqs) < 0.01\n    features['fft_low_freq_power'] = np.sum(fft_power[low_freq_mask])\n    features['fft_low_freq_ratio'] = features['fft_low_freq_power'] / features['fft_total_power']\n    \n    # Spectral properties\n    power_spectrum = fft_power[1:len(fft_power)//2]\n    freqs = np.abs(fft_freqs[1:len(fft_freqs)//2])\n    \n    if np.sum(power_spectrum) > 0:\n        features['spectral_centroid'] = np.sum(freqs * power_spectrum) / np.sum(power_spectrum)\n        features['spectral_bandwidth'] = np.sqrt(np.sum(((freqs - features['spectral_centroid'])**2) * power_spectrum) / np.sum(power_spectrum))\n    else:\n        features['spectral_centroid'] = 0\n        features['spectral_bandwidth'] = 0\n    \n    # Autocorrelation\n    autocorr = np.correlate(total_flux - np.mean(total_flux), total_flux - np.mean(total_flux), mode='full')\n    autocorr = autocorr[autocorr.size // 2:]\n    autocorr = autocorr / autocorr[0]\n    \n    # Lag features\n    for lag in [10, 50, 100, 500, 1000]:\n        if lag < len(autocorr):\n            features[f'autocorr_lag{lag}'] = autocorr[lag]\n    \n    # Peak detection\n    peaks, _ = signal.find_peaks(autocorr[1:1000], height=0.1)\n    features['autocorr_num_peaks'] = len(peaks)\n    features['autocorr_first_peak'] = peaks[0] if len(peaks) > 0 else 0\n    features['autocorr_strongest_peak'] = np.max(autocorr[peaks]) if len(peaks) > 0 else 0\n    \n    return features\n\n# Demo: Frequency Features\nfrequency_features = extract_frequency_features(total_flux)\nprint(f\"Generated {len(frequency_features)} frequency domain features:\")\n\n# Print ALL features organized by category\nprint(\"\\nPower Spectrum:\")\npower_features = ['fft_peak_power', 'fft_total_power', 'fft_mean_power', 'fft_std_power']\nfor feature in power_features:\n    print(f\"  {feature}: {frequency_features[feature]:.4f}\")\n\nprint(\"\\nLow Frequency Analysis:\")\nlow_freq_features = ['fft_low_freq_power', 'fft_low_freq_ratio']\nfor feature in low_freq_features:\n    print(f\"  {feature}: {frequency_features[feature]:.4f}\")\n\nprint(\"\\nSpectral Properties:\")\nspectral_features = ['spectral_centroid', 'spectral_bandwidth']\nfor feature in spectral_features:\n    print(f\"  {feature}: {frequency_features[feature]:.6f}\")\n\nprint(\"\\nAutocorrelation Lags:\")\nfor lag in [10, 50, 100, 500, 1000]:\n    feature = f'autocorr_lag{lag}'\n    if feature in frequency_features:\n        print(f\"  {feature}: {frequency_features[feature]:.6f}\")\n\nprint(\"\\nAutocorrelation Peaks:\")\npeak_features = ['autocorr_num_peaks', 'autocorr_first_peak', 'autocorr_strongest_peak']\nfor feature in peak_features:\n    print(f\"  {feature}: {frequency_features[feature]:.4f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T02:17:31.115386Z","iopub.execute_input":"2025-07-07T02:17:31.115683Z","iopub.status.idle":"2025-07-07T02:17:32.93645Z","shell.execute_reply.started":"2025-07-07T02:17:31.115658Z","shell.execute_reply":"2025-07-07T02:17:32.935766Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Spatial Features\n\nExtracts **30 spatial features** by analyzing how the star's position and shape change during transit across key time points.\n\n- **Centroid tracking**: Calculates precise center-of-mass position of star image at start, quarters, middle, end\n- **Concentration analysis**: Measures how much light is focused in center vs edges (8x8 region within 32x32 image)\n- **Movement calculation**: Tracks total path length and range of star position changes\n- **Shape characterization**: Mean brightness and uniformity statistics per time period\n\n**Key outputs**: `centroid_total_movement` (star drift), `centroid_x/y_range` (position shifts), `concentration_range` (shape changes), per-frame spatial statistics.","metadata":{}},{"cell_type":"code","source":"%%time\ndef extract_spatial_features(signal_data):\n    features = {}\n    \n    # Key time points\n    key_frames = [0, len(signal_data) // 4, len(signal_data) // 2, 3 * len(signal_data) // 4, -1]\n    centroids_x, centroids_y, concentrations = [], [], []\n    \n    for i, frame_idx in enumerate(key_frames):\n        frame = signal_data[frame_idx]\n        \n        # Centroid calculation\n        y_indices, x_indices = np.indices(frame.shape)\n        total_spatial_flux = np.sum(frame)\n        \n        if total_spatial_flux > 0:\n            centroid_x = np.sum(x_indices * frame) / total_spatial_flux\n            centroid_y = np.sum(y_indices * frame) / total_spatial_flux\n        else:\n            centroid_x = centroid_y = 16\n        \n        centroids_x.append(centroid_x)\n        centroids_y.append(centroid_y)\n        \n        # Concentration\n        center_region = frame[12:20, 12:20]\n        concentration = np.sum(center_region) / total_spatial_flux if total_spatial_flux > 0 else 0\n        concentrations.append(concentration)\n        \n        # Frame-specific features\n        features[f'frame{i}_spatial_mean'] = np.mean(frame)\n        features[f'frame{i}_spatial_std'] = np.std(frame)\n        features[f'frame{i}_centroid_x'] = centroid_x\n        features[f'frame{i}_centroid_y'] = centroid_y\n        features[f'frame{i}_concentration'] = concentration\n    \n    # Movement analysis\n    features['centroid_x_range'] = np.max(centroids_x) - np.min(centroids_x)\n    features['centroid_y_range'] = np.max(centroids_y) - np.min(centroids_y)\n    features['centroid_total_movement'] = np.sum(np.sqrt(np.diff(centroids_x)**2 + np.diff(centroids_y)**2))\n    features['concentration_range'] = np.max(concentrations) - np.min(concentrations)\n    features['concentration_std'] = np.std(concentrations)\n    \n    return features\n\n# Demo: Spatial Features\nspatial_features = extract_spatial_features(signal_data)\nprint(f\"Generated {len(spatial_features)} spatial features:\")\n\n# Print ALL features organized by frame and analysis type\nprint(\"\\nPer-Frame Spatial Statistics:\")\nfor i in range(5):\n    print(f\"  Frame {i}:\")\n    print(f\"    spatial_mean: {spatial_features[f'frame{i}_spatial_mean']:.4f}\")\n    print(f\"    spatial_std: {spatial_features[f'frame{i}_spatial_std']:.4f}\")\n    print(f\"    centroid_x: {spatial_features[f'frame{i}_centroid_x']:.4f}\")\n    print(f\"    centroid_y: {spatial_features[f'frame{i}_centroid_y']:.4f}\")\n    print(f\"    concentration: {spatial_features[f'frame{i}_concentration']:.6f}\")\n\nprint(\"\\nMovement Analysis:\")\nmovement_features = ['centroid_x_range', 'centroid_y_range', 'centroid_total_movement', \n                    'concentration_range', 'concentration_std']\nfor feature in movement_features:\n    print(f\"  {feature}: {spatial_features[feature]:.6f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T02:17:32.937333Z","iopub.execute_input":"2025-07-07T02:17:32.937553Z","iopub.status.idle":"2025-07-07T02:17:32.948384Z","shell.execute_reply.started":"2025-07-07T02:17:32.937536Z","shell.execute_reply":"2025-07-07T02:17:32.947718Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Gradient & Change Features\n\nExtracts **36 change-pattern features** by analyzing rates of brightness change and dividing observation into segments.\n\n- **Derivatives**: Calculates first (rate of change) and second (acceleration) derivatives of brightness\n- **Change detection**: Identifies periods of high volatility using rolling windows (10, 50, 100 frames)\n- **Segment analysis**: Divides 3.75-hour observation into 4 quarters, calculates trend and statistics per segment\n- **Transit identification**: Counts segments with brightness significantly below global average\n\n**Key outputs**: `flux_diff1_range` (change magnitude), `flux_diff2_extremes` (rapid acceleration events), `diff_volatility_max` (peak change periods), `segment{i}_trend` (linear trends per quarter), `num_transit_segments` (quarters containing transit).","metadata":{}},{"cell_type":"code","source":"%%time\ndef extract_gradient_features(total_flux):\n    features = {}\n    \n    # Derivatives\n    total_flux = np.array(total_flux, dtype=np.float64)\n    flux_diff1 = np.diff(total_flux)\n    flux_diff2 = np.diff(flux_diff1)\n    \n    # First derivative statistics\n    features['flux_diff1_mean'] = np.mean(flux_diff1)\n    features['flux_diff1_std'] = np.std(flux_diff1)\n    features['flux_diff1_min'] = np.min(flux_diff1)\n    features['flux_diff1_max'] = np.max(flux_diff1)\n    features['flux_diff1_range'] = features['flux_diff1_max'] - features['flux_diff1_min']\n    features['flux_diff1_skew'] = skew(flux_diff1)\n    features['flux_diff1_kurtosis'] = kurtosis(flux_diff1)\n    \n    # Second derivative statistics\n    features['flux_diff2_mean'] = np.mean(flux_diff2)\n    features['flux_diff2_std'] = np.std(flux_diff2)\n    features['flux_diff2_extremes'] = np.sum(np.abs(flux_diff2) > 3 * np.std(flux_diff2))\n    \n    # Change point detection\n    for window in [10, 50, 100]:\n        if window < len(flux_diff1):\n            rolling_diff_std = pd.Series(flux_diff1).rolling(window=window).std()\n            features[f'diff_volatility_w{window}_max'] = rolling_diff_std.max()\n            features[f'diff_volatility_w{window}_mean'] = rolling_diff_std.mean()\n    \n    # Segment analysis\n    segments = 4\n    segment_size = len(total_flux) // segments\n    segment_means = []\n    \n    for i in range(segments):\n        start_idx = i * segment_size\n        end_idx = (i + 1) * segment_size if i < segments - 1 else len(total_flux)\n        segment_flux = total_flux[start_idx:end_idx]\n        \n        if len(segment_flux) > 1:\n            trend = np.polyfit(np.arange(len(segment_flux)), segment_flux, 1)[0]\n            features[f'segment{i}_trend'] = trend\n            features[f'segment{i}_mean'] = np.mean(segment_flux)\n            features[f'segment{i}_std'] = np.std(segment_flux)\n            features[f'segment{i}_range'] = np.max(segment_flux) - np.min(segment_flux)\n            segment_means.append(features[f'segment{i}_mean'])\n    \n    # Cross-segment analysis\n    features['segment_mean_range'] = np.max(segment_means) - np.min(segment_means)\n    features['segment_mean_std'] = np.std(segment_means)\n    \n    # Transit segment detection\n    global_mean = np.mean(total_flux)\n    global_std = np.std(total_flux)\n    transit_segments = sum(1 for mean in segment_means if mean < global_mean - global_std)\n    \n    features['num_transit_segments'] = transit_segments\n    features['transit_segment_fraction'] = transit_segments / segments\n    \n    return features\n\n# Demo: Gradient & Change Features\ngradient_features = extract_gradient_features(total_flux)\nprint(f\"Generated {len(gradient_features)} gradient and change features:\")\n\n# Print ALL features organized by category\nprint(\"\\nFirst Derivative Statistics:\")\ndiff1_features = ['flux_diff1_mean', 'flux_diff1_std', 'flux_diff1_min', 'flux_diff1_max', \n                 'flux_diff1_range', 'flux_diff1_skew', 'flux_diff1_kurtosis']\nfor feature in diff1_features:\n    print(f\"  {feature}: {gradient_features[feature]:.4f}\")\n\nprint(\"\\nSecond Derivative Statistics:\")\ndiff2_features = ['flux_diff2_mean', 'flux_diff2_std', 'flux_diff2_extremes']\nfor feature in diff2_features:\n    print(f\"  {feature}: {gradient_features[feature]:.4f}\")\n\nprint(\"\\nChange Point Detection:\")\nfor window in [10, 50, 100]:\n    max_feature = f'diff_volatility_w{window}_max'\n    mean_feature = f'diff_volatility_w{window}_mean'\n    if max_feature in gradient_features:\n        print(f\"  {max_feature}: {gradient_features[max_feature]:.4f}\")\n        print(f\"  {mean_feature}: {gradient_features[mean_feature]:.4f}\")\n\nprint(\"\\nSegment Analysis (4 quarters):\")\nfor i in range(4):\n    print(f\"  Segment {i}:\")\n    print(f\"    trend: {gradient_features[f'segment{i}_trend']:.4f}\")\n    print(f\"    mean: {gradient_features[f'segment{i}_mean']:.4f}\")\n    print(f\"    std: {gradient_features[f'segment{i}_std']:.4f}\")\n    print(f\"    range: {gradient_features[f'segment{i}_range']:.4f}\")\n\nprint(\"\\nCross-Segment Analysis:\")\ncross_features = ['segment_mean_range', 'segment_mean_std', 'num_transit_segments', 'transit_segment_fraction']\nfor feature in cross_features:\n    print(f\"  {feature}: {gradient_features[feature]:.4f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T02:17:32.949253Z","iopub.execute_input":"2025-07-07T02:17:32.949547Z","iopub.status.idle":"2025-07-07T02:17:33.054021Z","shell.execute_reply.started":"2025-07-07T02:17:32.949522Z","shell.execute_reply":"2025-07-07T02:17:33.053256Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Single function that extracts all these features...\n* If you want to enable / disable a specific set of features - do it here!","metadata":{}},{"cell_type":"code","source":"%%time\ndef extract_enhanced_transit_features(signal_data, verbose=True):\n\n    n_frames = signal_data.shape[0]\n    signal_data = apply_adc_correction(signal_data, instrument='FGS1')\n    total_flux = np.sum(signal_data, axis=(1, 2))\n    \n    if verbose:\n        print(f\"Extracting enhanced transit features from {n_frames} frames...\")\n    \n    # Extract features by category (add or remove feature sets here)\n    features = {}\n    features.update(extract_global_flux_features(total_flux))\n    features.update(extract_rolling_statistics_features(total_flux))\n    features.update(extract_transit_detection_features(total_flux))\n    features.update(extract_frequency_features(total_flux))\n    features.update(extract_spatial_features(signal_data))\n    features.update(extract_gradient_features(total_flux))\n    \n    if verbose:\n        print(f\"Generated {len(features)} enhanced transit features\")\n    \n    return features\n\nall_features = extract_enhanced_transit_features(signal_data, verbose=True)\nprint(f\"\\nFinal feature summary:\")\nprint(f\"  Total features: {len(all_features)}\")\nprint(f\"  Feature categories: 6\")\nprint(f\"  Ready for LGBM training!\")\n\n# Show final key features\nprint(f\"\\nKey Enhanced Transit Features:\")\nprint(f\"  Global flux depth: {all_features['global_flux_depth']:.2f}\")\nprint(f\"  Global flux depth ratio: {all_features['global_flux_depth_ratio']:.4f}\")\nprint(f\"  Longest dip duration: {all_features['longest_dip_duration']} frames\")\nprint(f\"  Deepest time fraction: {all_features['deepest_time_fraction']:.3f}\")\nprint(f\"  Rolling1000 deepest dip: {all_features['rolling1000_deepest_dip']:.2f}\")\nprint(f\"  Number of transit segments: {all_features['num_transit_segments']}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T02:17:33.054809Z","iopub.execute_input":"2025-07-07T02:17:33.055046Z","iopub.status.idle":"2025-07-07T02:17:35.556608Z","shell.execute_reply.started":"2025-07-07T02:17:33.055029Z","shell.execute_reply":"2025-07-07T02:17:35.555759Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Do feature extraction on all training data\n* Takes about 2 hours for train data\n* Handles **multiple signals per planet by creating additional records** (planet_id_0, planet_id_0, etc)\n* Function is called in create_submission() to process test data (also about 2 hours)","metadata":{}},{"cell_type":"code","source":"%%time\n\nimport pandas as pd\nimport numpy as np\nimport glob\n\ndef prepare_data_with_enhanced_fgs1(train_df, star_info_df, data_path):\n    \"\"\"Enhanced FGS1 feature extractor — always processes multiple signals per planet with string IDs.\"\"\"\n\n    print(f\"Extracting enhanced FGS1 features for {len(star_info_df)} planets...\")\n    fgs1_features = []\n\n    # Ensure string-based planet IDs throughout\n    star_info_df['planet_id'] = star_info_df['planet_id'].astype(str)\n\n    if train_df is not None:\n        train_df['planet_id'] = train_df['planet_id'].astype(str)\n\n    # Extract features from all available FGS1 signals\n    for i, row in star_info_df.iterrows():\n        base_id = int(float(row['planet_id']))\n        print(f\"\\rProcessing planet {i+1}/{len(star_info_df)} (ID: {base_id})\", end='', flush=True)\n\n        try:\n            signal_paths = sorted(glob.glob(f\"{data_path}{base_id}/FGS1_signal_*.parquet\"))\n\n            for j, path in enumerate(signal_paths):\n                df = pd.read_parquet(path)\n                signal = df.values.reshape(135000, 32, 32)\n                features = extract_enhanced_transit_features(signal, verbose=False)\n                signal_id = f\"{base_id}_{j}\"\n                features['planet_id'] = signal_id\n                fgs1_features.append(features)\n\n        except Exception as e:\n            print(f\"\\n❌ Failed to process planet {base_id}: {e}\")\n            continue\n\n    print(\"\\n✅ Feature extraction complete.\")\n\n    # Build features dataframe\n    features_df = pd.DataFrame(fgs1_features)\n    features_df['planet_id'] = features_df['planet_id'].astype(str)\n    \n    features_df = features_df.set_index('planet_id')\n    print(f\"→ Extracted features for {len(features_df)} entries with {features_df.shape[1]} columns\")\n    \n    # Expand star metadata for each signal\n    expanded_meta = []\n\n    # Normalize IDs for safe comparison\n    star_info_df['planet_id'] = star_info_df['planet_id'].apply(lambda x: str(int(float(x))))\n    \n    for pid in features_df.index:\n        base_id = str(int(float(pid.split(\"_\")[0])))\n        row = star_info_df[star_info_df['planet_id'] == base_id].copy()\n        if not row.empty:\n            row['planet_id'] = pid\n            expanded_meta.append(row)\n\n    star_info_df = pd.concat(expanded_meta, ignore_index=True)\n\n    # Merge metadata and features\n    full_df = star_info_df.set_index('planet_id').join(features_df, how='left')\n    X = full_df.select_dtypes(include=[np.number]).fillna(0).astype(np.float32)\n\n    if train_df is not None:\n        targets_df = train_df.set_index('planet_id')\n        extended_targets = []\n        for pid in X.index:\n            base_id = pid.split(\"_\")[0]\n            if base_id in targets_df.index:\n                y_row = targets_df.loc[base_id].copy()\n                y_row.name = pid\n                extended_targets.append(y_row)\n        y = pd.DataFrame(extended_targets).astype(np.float32)\n\n        print(f\"✅ Final shapes — X: {X.shape}, y: {y.shape}\")\n        return X, y\n\n    else:\n        print(f\"✅ Final shape — X_test: {X.shape}\")\n        return X\n                \nX, y = prepare_data_with_enhanced_fgs1(train_df, train_star_info_df, \"/kaggle/input/ariel-data-challenge-2025/train/\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T03:07:31.077866Z","iopub.execute_input":"2025-07-07T03:07:31.078556Z","iopub.status.idle":"2025-07-07T03:08:03.23859Z","shell.execute_reply.started":"2025-07-07T03:07:31.078529Z","shell.execute_reply":"2025-07-07T03:08:03.237748Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Train models for each target!\n* Trains *283 separate LGBM models per CV fold*.\n* With 5 CV folds - **that's 1415 LGBM models total!**\n* Collects out-of-fold (OOF) predictions on validation splits that we use later for uncertainties","metadata":{}},{"cell_type":"code","source":"%%time\n\ndef train_validate_multioutput_cv(X, y, n_splits=5):\n    print(f\"Training with CV: {X.shape[1]} features → {y.shape[1]} wavelengths\")\n    \n    kf = KFold(n_splits=n_splits, shuffle=True, random_state=42)\n\n    all_models = [[] for _ in range(y.shape[1])]\n    all_fold_metrics = []\n\n    oof_preds = np.zeros_like(y.values, dtype=float)\n    oof_sigmas = np.zeros_like(y.values, dtype=float)\n    residuals_all = [[] for _ in range(y.shape[1])]\n\n    for fold, (train_idx, val_idx) in enumerate(kf.split(X)):\n        print(f\"\\n🔁 Fold {fold+1}/{n_splits}\")\n        X_train, X_val = X.iloc[train_idx], X.iloc[val_idx]\n        y_train, y_val = y.iloc[train_idx], y.iloc[val_idx]\n\n        for i in tqdm(range(y.shape[1]), desc=f\"Training wavelength models for Fold {fold+1}\"):\n            model = lgb.LGBMRegressor(\n                n_estimators=500,\n                learning_rate=0.05,\n                num_leaves=127,\n                max_depth=8,\n                min_child_samples=20,\n                subsample=0.8,\n                colsample_bytree=0.8,\n                reg_alpha=0.1,\n                reg_lambda=0.1,\n                random_state=42,\n                verbose=-1\n            )\n\n            y_single = y_train.iloc[:, i]\n            model.fit(X_train, y_single)\n            all_models[i].append(model)\n\n            # Predict on validation set (OOF)\n            pred_val = model.predict(X_val)\n            oof_preds[val_idx, i] = pred_val\n\n            # Store residuals for uncertainty estimation\n            pred_train = model.predict(X_train)\n            residuals = y_single - pred_train\n            residuals_all[i].extend(residuals.tolist())\n\n            # Temporary sigma just for oof output\n            sigma = np.std(residuals)\n            oof_sigmas[val_idx, i] = max(sigma, 1e-6)\n\n    # Compute final stddevs from full OOF residuals\n    target_uncertainties = [max(np.std(res), 1e-6) for res in residuals_all]\n\n    # Compute metrics\n    r2_scores = [r2_score(y.values[:, i], oof_preds[:, i]) for i in range(y.shape[1])]\n    rmses = [mean_squared_error(y.values[:, i], oof_preds[:, i], squared=False) for i in range(y.shape[1])]\n\n    print(\"\\n📊 CV Performance Summary:\")\n    print(f\"→ Mean R²:   {np.mean(r2_scores):.4f}\")\n    print(f\"→ Mean RMSE: {np.mean(rmses):.4f}\")\n\n    return {\n        'models': all_models,\n        'feature_columns': X.columns.tolist(),\n        'target_columns': y.columns.tolist(),\n        'oof_predictions': oof_preds,\n        'oof_uncertainties': oof_sigmas,\n        'y_true': y.values,\n        'cv_metrics': {\n            'r2_per_target': r2_scores,\n            'rmse_per_target': rmses,\n            'mean_r2': np.mean(r2_scores),\n            'mean_rmse': np.mean(rmses)\n        },\n        'target_uncertainties': target_uncertainties\n    }\n\ntrained_models_cv = train_validate_multioutput_cv(X, y)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T02:24:00.394094Z","iopub.execute_input":"2025-07-07T02:24:00.3944Z","iopub.status.idle":"2025-07-07T02:25:31.735109Z","shell.execute_reply.started":"2025-07-07T02:24:00.394372Z","shell.execute_reply":"2025-07-07T02:25:31.734268Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Feature Importance\n* Saves out to a file for later review\n* It looks like only a few of the features are doing the \"heavy lifting\"\n* LGBM does a great job ignoring unimportant features\n* Some features might become useful with further adjustment and refinement","metadata":{}},{"cell_type":"code","source":"def analyze_feature_importance(trained_models, X, save_file=True):\n    \"\"\"Analyze feature importance across all trained models (supports CV version)\"\"\"\n    \n    print(f\"\\n🔍 Feature Importance Analysis:\")\n    \n    all_models_nested = trained_models['models']  # list of lists: [n_targets][n_folds]\n    all_models_flat = [m for models_per_target in all_models_nested for m in models_per_target]\n    \n    # Calculate average feature importance across all models\n    feature_importance = pd.DataFrame({\n        'feature': trained_models['feature_columns'],\n        'importance': np.mean([model.feature_importances_ for model in all_models_flat], axis=0)\n    }).sort_values('importance', ascending=False)\n    \n    # Display top features\n    print(\"Most important features:\")\n    for i, row in feature_importance.head(25).iterrows():\n        print(f\"  {i+1:2d}. {row['feature']:<35} {row['importance']:.4f}\")\n    \n    # Feature categories analysis\n    categories = {}\n    for _, row in feature_importance.iterrows():\n        category = row['feature'].split('_')[0] if '_' in row['feature'] else 'other'\n        categories.setdefault(category, []).append(row['importance'])\n    \n    print(f\"\\nFeature Category Performance:\")\n    for category, importances in sorted(categories.items(), key=lambda x: np.sum(x[1]), reverse=True):\n        total_contrib = np.sum(importances) / feature_importance['importance'].sum() * 100\n        print(f\"  {category:<15}: {len(importances):>3d} features, {total_contrib:>5.1f}% contribution\")\n    \n    # Quick stats\n    zero_features = (feature_importance['importance'] == 0).sum()\n    print(f\"\\nQuick Stats:\")\n    print(f\"  Zero importance features: {zero_features}\")\n    print(f\"  Top feature: {feature_importance.iloc[0]['feature']}\")\n    \n    # Save to file\n    if save_file:\n        feature_importance.to_csv('feature_importance.csv', index=False)\n        print(f\"\\n💾 Saved feature importance to: feature_importance.csv\")\n    \n    return feature_importance\n    \nfeature_importance_df = analyze_feature_importance(trained_models_cv, X)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T02:25:31.736168Z","iopub.execute_input":"2025-07-07T02:25:31.736432Z","iopub.status.idle":"2025-07-07T02:25:31.880291Z","shell.execute_reply.started":"2025-07-07T02:25:31.736404Z","shell.execute_reply":"2025-07-07T02:25:31.879578Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Perform Inference\n* Predicts targets using corresponding trained models\n* Estimates uncertainty for each prediction using training residuals\n* Can override uncertainty fixed value / be sure to update in create_submission()","metadata":{}},{"cell_type":"code","source":"def predict_with_uncertainty(models, X, fixed_uncertainty=None, target_uncertainties=None):\n    \"\"\"\n    Make predictions with uncertainty estimates using trained models.\n\n    Args:\n        models: List of trained models (flat or nested by CV)\n        X: Features to predict on\n        fixed_uncertainty: Use this fixed value for all targets\n        target_uncertainties: Optional list of std values per target (from OOF residuals)\n\n    Returns:\n        tuple: (predictions, uncertainties) as numpy arrays\n    \"\"\"\n    predictions = []\n    uncertainties = []\n\n    is_cv = isinstance(models[0], list)\n    n_targets = len(models)\n\n    for i in range(n_targets):\n        model_group = models[i] if is_cv else [models[i]]\n        preds = []\n\n        for model in model_group:\n            preds.append(model.predict(X))\n\n        # Average predictions over folds\n        pred_avg = np.mean(preds, axis=0)\n        predictions.append(pred_avg)\n\n        # Use appropriate uncertainty\n        if fixed_uncertainty is not None:\n            unc = fixed_uncertainty\n        elif target_uncertainties is not None:\n            unc = target_uncertainties[i]\n        else:\n            unc = 0.01  # fallback\n\n        unc_array = np.full_like(pred_avg, max(unc, 1e-6))\n        uncertainties.append(unc_array)\n\n    y_pred = np.column_stack(predictions)\n    sigma_pred = np.column_stack(uncertainties)\n\n    # Clean predictions\n    pred_df = pd.DataFrame(y_pred).apply(pd.to_numeric, errors='coerce').fillna(0).clip(lower=0)\n    sigma_df = pd.DataFrame(sigma_pred).apply(pd.to_numeric, errors='coerce').fillna(1e-6).clip(lower=1e-15)\n\n    return pred_df.values, sigma_df.values\n    ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T02:48:23.978369Z","iopub.execute_input":"2025-07-07T02:48:23.978737Z","iopub.status.idle":"2025-07-07T02:48:23.986825Z","shell.execute_reply.started":"2025-07-07T02:48:23.978709Z","shell.execute_reply":"2025-07-07T02:48:23.985965Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Estimate LB Score\n* Metric code imported from https://www.kaggle.com/code/metric/ariel-gaussian-log-likelihood\n* Does Out-of-Fold prediction - shouldn't be leakage\n* Earlier bug has been corrected - hopefully accurate now!","metadata":{}},{"cell_type":"code","source":"sys.path.append('/kaggle/usr/lib/ariel-gaussian-log-likelihood')\nfrom metric import score\n\ndef evaluate_with_gll_metric(y_true, y_pred, sigma_pred, naive_mean, naive_sigma):\n    \"\"\"\n    Evaluate predictions using the official competition GLL metric.\n\n    Args:\n        y_true: Ground truth DataFrame or array (samples x wavelengths)\n        y_pred: Predictions array (samples x wavelengths)\n        sigma_pred: Uncertainty estimates array (samples x wavelengths)\n        naive_mean: Scalar from training set (mean)\n        naive_sigma: Scalar from training set (std)\n\n    Returns:\n        float: GLL score [0, 1]\n    \"\"\"\n    # Ensure y_true is DataFrame\n    if not isinstance(y_true, pd.DataFrame):\n        y_true = pd.DataFrame(y_true)\n\n    y_true = y_true.reset_index(drop=True)\n    n_samples, n_waves = y_pred.shape\n\n    # Ensure no negative or zero sigma\n    sigma_pred = np.clip(sigma_pred, 1e-15, None)\n\n    # Build solution DataFrame\n    solution_df = y_true.copy()\n    solution_df['row_id'] = np.arange(n_samples)\n\n    # Build submission DataFrame in exact format\n    submission_df = pd.DataFrame()\n    for i in range(n_waves):\n        submission_df[f'wavelength_{i}'] = y_pred[:, i]\n    for i in range(n_waves):\n        submission_df[f'wavelength_{i}_std'] = sigma_pred[:, i]\n    submission_df['row_id'] = np.arange(n_samples)\n\n    # Ensure non-negative values (required)\n    submission_df = submission_df.clip(lower=1e-15)\n\n    # Call official scorer\n    gll_score = score(\n        solution=solution_df,\n        submission=submission_df,\n        row_id_column_name='row_id',\n        naive_mean=naive_mean,\n        naive_sigma=naive_sigma,\n        fsg_sigma_true=1e-6,\n        airs_sigma_true=1e-5,\n        fgs_weight=0.4\n    )\n\n    return gll_score\n\n# use OOF predictions and uncertainties collected during training\ny_true = trained_models_cv['y_true']\ny_pred = trained_models_cv['oof_predictions']\nsigma_pred = trained_models_cv['oof_uncertainties']\n\n# Calculate baseline stats using full data\nnaive_mean = y_true.mean()\nnaive_sigma = y_true.std()\n\n# Evaluate with GLL metric\ngll_score = evaluate_with_gll_metric(y_true, y_pred, sigma_pred, naive_mean, naive_sigma)\nprint(f\"🎯 OOF Gaussian Log Likelihood Score: {gll_score:.6f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T02:48:26.740622Z","iopub.execute_input":"2025-07-07T02:48:26.741453Z","iopub.status.idle":"2025-07-07T02:48:27.057865Z","shell.execute_reply.started":"2025-07-07T02:48:26.741424Z","shell.execute_reply":"2025-07-07T02:48:27.056997Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Optimize Sigmas for LB Score Per Wavelength\n* Re-implements competition metric using numpy (and without dataframes) for performance\n* Start with our training residual-based sigma predictions - and then determine optimal adjustments for each wave-length using minimize_scalar\n* Finally - rescore using official metric code\n* Note: Sigma optimization hasn't produced the gains it has for other notebooks.  Unclear if this is a difference in implementation - or the residual error approach is already optimal.\n","metadata":{}},{"cell_type":"code","source":"def fast_gll_score_numpy(y_true, y_pred, sigma_pred, naive_mean, naive_sigma,\n                         fsg_sigma_true=1e-6, airs_sigma_true=1e-5, fgs_weight=0.4):\n    \"\"\"\n    Fast NumPy-based GLL score computation (no DataFrames, optimized for speed).\n    \"\"\"\n    sigma_pred = np.clip(sigma_pred, 1e-15, None)\n    n_samples, n_waves = sigma_pred.shape\n\n    sigma_true = np.append([fsg_sigma_true], np.full(n_waves - 1, airs_sigma_true))\n    sigma_true = np.tile(sigma_true, (n_samples, 1))\n\n    weights = np.append([fgs_weight], np.ones(n_waves - 1))\n    weights = np.tile(weights, (n_samples, 1))\n\n    gll_pred = norm.logpdf(y_true, loc=y_pred, scale=sigma_pred)\n    gll_true = norm.logpdf(y_true, loc=y_true, scale=sigma_true)\n    gll_naive = norm.logpdf(y_true, loc=naive_mean, scale=naive_sigma)\n\n    ind_scores = (gll_pred - gll_naive) / (gll_true - gll_naive + 1e-9)\n    final_score = np.average(ind_scores, weights=weights)\n    return float(np.clip(final_score, 0.0, 1.0))\n\n\ndef optimize_sigma_per_wavelength(y_true, y_pred, sigma_pred, naive_mean, naive_sigma,\n                                       fsg_sigma_true=1e-6, airs_sigma_true=1e-5, fgs_weight=0.4):\n    \"\"\"\n    Optimize sigma scaling per wavelength using fast GLL scoring.\n\n    Returns:\n        best_scales: np.ndarray of optimal scalers (shape: n_wavelengths,)\n    \"\"\"\n    n_waves = sigma_pred.shape[1]\n    best_scales = np.ones(n_waves)\n\n    print(f\"⚡ Optimizing {n_waves} sigma scalers with fast NumPy GLL metric...\\n\")\n\n    for i in tqdm(range(n_waves), desc=\"Optimizing wavelengths\", unit=\"λ\"):\n        def objective(scale):\n            sigma_scaled = sigma_pred.copy()\n            sigma_scaled[:, i] *= scale\n            return -fast_gll_score_numpy(\n                y_true, y_pred, sigma_scaled,\n                naive_mean, naive_sigma,\n                fsg_sigma_true=fsg_sigma_true,\n                airs_sigma_true=airs_sigma_true,\n                fgs_weight=fgs_weight\n            )\n\n        result = minimize_scalar(objective, bounds=(0.25, 4.0), method='bounded')\n        best_scales[i] = result.x\n\n    print(\"✅ Optimization complete.\")\n    return best_scales\n    \n\n# Optimize per-dimension scale factors\nbest_scales = optimize_sigma_per_wavelength(\n    y_true,\n    y_pred,\n    sigma_pred,\n    naive_mean,\n    naive_sigma\n)\n\n# Final GLL score - sigma_pred modified by best_scales\nfinal_score = evaluate_with_gll_metric(\n    y_true, y_pred, sigma_pred * best_scales,\n    naive_mean, naive_sigma\n)\nprint(f\"🎯 Calibrated GLL Score: {final_score:.6f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T02:25:32.23833Z","iopub.execute_input":"2025-07-07T02:25:32.238595Z","iopub.status.idle":"2025-07-07T02:25:43.357344Z","shell.execute_reply.started":"2025-07-07T02:25:32.238568Z","shell.execute_reply":"2025-07-07T02:25:43.356417Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Create our submission\n* Applies our optimized per-wavelength sigma scaling to improve LB score\n* Multiple signals result in multiple predictions (that are then averaged together)","metadata":{"execution":{"iopub.status.busy":"2025-07-06T22:00:01.050691Z","iopub.execute_input":"2025-07-06T22:00:01.052348Z","iopub.status.idle":"2025-07-06T22:00:01.059613Z","shell.execute_reply.started":"2025-07-06T22:00:01.052304Z","shell.execute_reply":"2025-07-06T22:00:01.057909Z"}}},{"cell_type":"code","source":"def create_submission(trained_models, test_star_info_df, output_path='submission.csv', sigma_scalers=None):\n    \"\"\"\n    Generate a submission file from trained models and test data, averaging predictions across multiple signals per planet.\n    \"\"\"\n    # Prepare test data using your custom feature engineering\n    X_test = prepare_data_with_enhanced_fgs1(\n        train_df=None,\n        star_info_df=test_star_info_df,\n        data_path=\"/kaggle/input/ariel-data-challenge-2025/test/\"\n    )\n\n    # Align test data columns with training features\n    X_test_aligned = X_test.reindex(columns=trained_models['feature_columns'], fill_value=0)\n\n    # Predict with uncertainty\n    y_pred, sigma_pred = predict_with_uncertainty(\n        trained_models['models'],\n        X_test_aligned,\n        target_uncertainties=trained_models['target_uncertainties']\n    )\n\n    # Apply per-wavelength scaling to sigma_pred, if provided\n    if sigma_scalers is not None:\n        sigma_pred = sigma_pred * sigma_scalers\n\n    # Parse base planet IDs from signal IDs like \"12345_0\"\n    base_ids = [idx.split(\"_\")[0] for idx in X_test_aligned.index]\n\n    # Convert predictions to DataFrames with signal IDs\n    wl_cols = trained_models['target_columns']\n    sigma_cols = [f\"sigma_{i+1}\" for i in range(y_pred.shape[1])]\n    \n    pred_df = pd.DataFrame(y_pred, index=base_ids, columns=wl_cols)\n    sigma_df = pd.DataFrame(sigma_pred, index=base_ids, columns=sigma_cols)\n\n    # Average over signals per base planet\n    pred_mean = pred_df.groupby(pred_df.index).mean()\n    sigma_mean = sigma_df.groupby(sigma_df.index).mean()\n\n    # Build submission DataFrame\n    submission_df = pd.concat([pred_mean, sigma_mean], axis=1).reset_index()\n    submission_df = submission_df.rename(columns={'index': 'planet_id'})\n    submission_df['planet_id'] = submission_df['planet_id'].astype(int)\n\n    # Save to file\n    submission_df.to_csv(output_path, index=False, float_format='%.5f')\n    print(f\"✅ Submission saved to: {output_path}\")\n\n    return submission_df\n    \nsubmission = create_submission(\n    trained_models_cv,\n    test_star_info_df,\n    output_path='submission.csv',\n    sigma_scalers=best_scales \n)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T03:05:11.286552Z","iopub.execute_input":"2025-07-07T03:05:11.286912Z","iopub.status.idle":"2025-07-07T03:05:21.430957Z","shell.execute_reply.started":"2025-07-07T03:05:11.286887Z","shell.execute_reply":"2025-07-07T03:05:21.430232Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-07T03:05:25.137087Z","iopub.execute_input":"2025-07-07T03:05:25.137636Z","iopub.status.idle":"2025-07-07T03:05:25.153797Z","shell.execute_reply.started":"2025-07-07T03:05:25.137611Z","shell.execute_reply":"2025-07-07T03:05:25.152895Z"}},"outputs":[],"execution_count":null}]}