{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","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":13093295,"sourceType":"competition"}],"dockerImageVersionId":31089,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"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":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-07-27T21:02:03.520505Z","iopub.execute_input":"2025-07-27T21:02:03.520809Z","iopub.status.idle":"2025-07-27T21:02:11.172756Z","shell.execute_reply.started":"2025-07-27T21:02:03.520782Z","shell.execute_reply":"2025-07-27T21:02:11.171578Z"}},"outputs":[],"execution_count":null},{"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,"execution":{"iopub.status.busy":"2025-07-27T21:02:19.052866Z","iopub.execute_input":"2025-07-27T21:02:19.054258Z","iopub.status.idle":"2025-07-27T21:02:19.300527Z","shell.execute_reply.started":"2025-07-27T21:02:19.054214Z","shell.execute_reply":"2025-07-27T21:02:19.299518Z"}},"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-27T21:02:23.423542Z","iopub.execute_input":"2025-07-27T21:02:23.423985Z","iopub.status.idle":"2025-07-27T21:02:23.431398Z","shell.execute_reply.started":"2025-07-27T21:02:23.423937Z","shell.execute_reply":"2025-07-27T21:02:23.430028Z"}},"outputs":[],"execution_count":null},{"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,"execution":{"iopub.status.busy":"2025-07-27T21:02:26.533235Z","iopub.execute_input":"2025-07-27T21:02:26.533743Z","iopub.status.idle":"2025-07-27T21:02:28.812399Z","shell.execute_reply.started":"2025-07-27T21:02:26.533707Z","shell.execute_reply":"2025-07-27T21:02:28.811501Z"}},"outputs":[],"execution_count":null},{"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-27T21:02:31.822875Z","iopub.execute_input":"2025-07-27T21:02:31.82321Z","iopub.status.idle":"2025-07-27T21:02:32.24117Z","shell.execute_reply.started":"2025-07-27T21:02:31.823188Z","shell.execute_reply":"2025-07-27T21:02:32.240136Z"}},"outputs":[],"execution_count":null},{"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,"execution":{"iopub.status.busy":"2025-07-27T21:02:36.083305Z","iopub.execute_input":"2025-07-27T21:02:36.08371Z","iopub.status.idle":"2025-07-27T21:02:37.104846Z","shell.execute_reply.started":"2025-07-27T21:02:36.083681Z","shell.execute_reply":"2025-07-27T21:02:37.103614Z"}},"outputs":[],"execution_count":null},{"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-27T21:02:42.769899Z","iopub.execute_input":"2025-07-27T21:02:42.7703Z","iopub.status.idle":"2025-07-27T21:02:42.813027Z","shell.execute_reply.started":"2025-07-27T21:02:42.770267Z","shell.execute_reply":"2025-07-27T21:02:42.812103Z"}},"outputs":[],"execution_count":null},{"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-27T21:02:47.842122Z","iopub.execute_input":"2025-07-27T21:02:47.842442Z","iopub.status.idle":"2025-07-27T21:02:48.061105Z","shell.execute_reply.started":"2025-07-27T21:02:47.84242Z","shell.execute_reply":"2025-07-27T21:02:48.059854Z"}},"outputs":[],"execution_count":null},{"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-27T21:02:53.593842Z","iopub.execute_input":"2025-07-27T21:02:53.594156Z","iopub.status.idle":"2025-07-27T21:02:53.775879Z","shell.execute_reply.started":"2025-07-27T21:02:53.594135Z","shell.execute_reply":"2025-07-27T21:02:53.774634Z"}},"outputs":[],"execution_count":null},{"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-27T21:02:59.730234Z","iopub.execute_input":"2025-07-27T21:02:59.731114Z","iopub.status.idle":"2025-07-27T21:03:02.070115Z","shell.execute_reply.started":"2025-07-27T21:02:59.731081Z","shell.execute_reply":"2025-07-27T21:03:02.069105Z"}},"outputs":[],"execution_count":null},{"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-27T21:03:07.5461Z","iopub.execute_input":"2025-07-27T21:03:07.546498Z","iopub.status.idle":"2025-07-27T21:03:07.560168Z","shell.execute_reply.started":"2025-07-27T21:03:07.546471Z","shell.execute_reply":"2025-07-27T21:03:07.559241Z"}},"outputs":[],"execution_count":null},{"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-27T21:03:12.697629Z","iopub.execute_input":"2025-07-27T21:03:12.698346Z","iopub.status.idle":"2025-07-27T21:03:12.782637Z","shell.execute_reply.started":"2025-07-27T21:03:12.698317Z","shell.execute_reply":"2025-07-27T21:03:12.781615Z"}},"outputs":[],"execution_count":null},{"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-27T21:03:21.265896Z","iopub.execute_input":"2025-07-27T21:03:21.266266Z","iopub.status.idle":"2025-07-27T21:03:24.486425Z","shell.execute_reply.started":"2025-07-27T21:03:21.266241Z","shell.execute_reply":"2025-07-27T21:03:24.485402Z"}},"outputs":[],"execution_count":null},{"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-27T21:03:29.129865Z","iopub.execute_input":"2025-07-27T21:03:29.130235Z","iopub.status.idle":"2025-07-27T22:34:18.663259Z","shell.execute_reply.started":"2025-07-27T21:03:29.130207Z","shell.execute_reply":"2025-07-27T22:34:18.661789Z"}},"outputs":[],"execution_count":null},{"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-27T22:34:22.655962Z","iopub.execute_input":"2025-07-27T22:34:22.656384Z","iopub.status.idle":"2025-07-27T22:45:56.277604Z","shell.execute_reply.started":"2025-07-27T22:34:22.656352Z","shell.execute_reply":"2025-07-27T22:45:56.276535Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import os\nimport joblib\nimport pickle\n\ndef save_trained_models(model_dict, save_dir='trained_models'):\n    \"\"\"\n    Save trained models and related metadata to disk.\n\n    Parameters:\n    - model_dict: the dictionary returned by train_validate_multioutput_cv\n    - save_dir: directory where models and metadata will be stored\n    \"\"\"\n    os.makedirs(save_dir, exist_ok=True)\n\n    # Save models per target per fold\n    for target_idx, model_list in enumerate(model_dict['models']):\n        for fold_idx, model in enumerate(model_list):\n            model_path = os.path.join(save_dir, f'model_target{target_idx}_fold{fold_idx}.pkl')\n            joblib.dump(model, model_path)\n\n    # Save metadata (features, targets, metrics, uncertainties, etc.)\n    metadata = {\n        'feature_columns': model_dict['feature_columns'],\n        'target_columns': model_dict['target_columns'],\n        'cv_metrics': model_dict['cv_metrics'],\n        'target_uncertainties': model_dict['target_uncertainties']\n    }\n\n    with open(os.path.join(save_dir, 'trained_model_info.pkl'), 'wb') as f:\n        pickle.dump(metadata, f)\n\n    print(f\"✅ All models and metadata saved to: {save_dir}\")\nsave_trained_models(trained_models_cv, save_dir='trained_models')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-27T22:46:12.92168Z","iopub.execute_input":"2025-07-27T22:46:12.922565Z","iopub.status.idle":"2025-07-27T22:46:23.185769Z","shell.execute_reply.started":"2025-07-27T22:46:12.92252Z","shell.execute_reply":"2025-07-27T22:46:23.184751Z"}},"outputs":[],"execution_count":null},{"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-27T22:46:26.759019Z","iopub.execute_input":"2025-07-27T22:46:26.759415Z","iopub.status.idle":"2025-07-27T22:46:26.947304Z","shell.execute_reply.started":"2025-07-27T22:46:26.75939Z","shell.execute_reply":"2025-07-27T22:46:26.946126Z"}},"outputs":[],"execution_count":null},{"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-27T22:46:43.912046Z","iopub.execute_input":"2025-07-27T22:46:43.912741Z","iopub.status.idle":"2025-07-27T22:46:43.923689Z","shell.execute_reply.started":"2025-07-27T22:46:43.912702Z","shell.execute_reply":"2025-07-27T22:46:43.921946Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\n# Dummy GLL score function – replace with official one later\ndef score(solution, submission, row_id_column_name, naive_mean, naive_sigma,\n          fsg_sigma_true, airs_sigma_true, fgs_weight):\n    \"\"\"\n    Dummy version of GLL scoring function.\n\n    Returns a value between 0 and 1 (higher is better).\n    \"\"\"\n    y_true = solution.drop(columns=[row_id_column_name])\n    y_pred = submission[[col for col in submission.columns if not col.endswith('_std') and col != row_id_column_name]]\n    sigma = submission[[col for col in submission.columns if col.endswith('_std')]]\n\n    # Gaussian log likelihood (simplified)\n    log_likelihoods = -0.5 * np.log(2 * np.pi * sigma.values ** 2) - ((y_true.values - y_pred.values) ** 2) / (2 * sigma.values ** 2)\n    gll = np.mean(log_likelihoods)\n\n    # Normalize score between 0 and 1 (dummy scale for example)\n    return np.clip((gll + 10) / 10, 0, 1)  # Simulated range\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-27T22:52:49.621267Z","iopub.execute_input":"2025-07-27T22:52:49.622067Z","iopub.status.idle":"2025-07-27T22:52:49.632947Z","shell.execute_reply.started":"2025-07-27T22:52:49.622024Z","shell.execute_reply":"2025-07-27T22:52:49.631756Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def evaluate_with_gll_metric(y_true, y_pred, sigma_pred, naive_mean, naive_sigma):\n    \"\"\"\n    Evaluate predictions using GLL-like metric (dummy version).\n    \"\"\"\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\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 positive values\n    submission_df = submission_df.clip(lower=1e-15)\n\n    # Dummy GLL 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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-27T22:53:05.846601Z","iopub.execute_input":"2025-07-27T22:53:05.846994Z","iopub.status.idle":"2025-07-27T22:53:05.856582Z","shell.execute_reply.started":"2025-07-27T22:53:05.846968Z","shell.execute_reply":"2025-07-27T22:53:05.855037Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Dummy test data\nnp.random.seed(42)\nsamples, wavelengths = 100, 5\ny_true = pd.DataFrame(np.random.rand(samples, wavelengths))\ny_pred = np.random.rand(samples, wavelengths)\nsigma_pred = np.abs(np.random.rand(samples, wavelengths)) + 1e-3\n\nnaive_mean = y_true.mean().mean()\nnaive_sigma = y_true.stack().std()\n\n# Evaluate\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}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-27T22:53:20.111433Z","iopub.execute_input":"2025-07-27T22:53:20.112297Z","iopub.status.idle":"2025-07-27T22:53:20.150927Z","shell.execute_reply.started":"2025-07-27T22:53:20.112266Z","shell.execute_reply":"2025-07-27T22:53:20.149805Z"}},"outputs":[],"execution_count":null},{"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-27T22:53:27.32765Z","iopub.execute_input":"2025-07-27T22:53:27.328056Z","iopub.status.idle":"2025-07-27T22:53:27.441963Z","shell.execute_reply.started":"2025-07-27T22:53:27.328029Z","shell.execute_reply":"2025-07-27T22:53:27.440893Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\ndef 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\n\n# ✅ Define best_scales (1.0 = no scaling; adjust if needed)\nbest_scales = np.ones(len(trained_models_cv['target_columns']))\n\n\n# ✅ Create the submission\nsubmission = create_submission(\n    trained_models_cv,\n    test_star_info_df,\n    output_path='submission.csv',\n    sigma_scalers=best_scales\n)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-27T22:53:54.176437Z","iopub.execute_input":"2025-07-27T22:53:54.17718Z","iopub.status.idle":"2025-07-27T22:54:07.580274Z","shell.execute_reply.started":"2025-07-27T22:53:54.177152Z","shell.execute_reply":"2025-07-27T22:54:07.579167Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-27T22:54:28.479207Z","iopub.execute_input":"2025-07-27T22:54:28.479568Z","iopub.status.idle":"2025-07-27T22:54:28.501997Z","shell.execute_reply.started":"2025-07-27T22:54:28.479543Z","shell.execute_reply":"2025-07-27T22:54:28.500742Z"}},"outputs":[],"execution_count":null}]}