{"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":"gpu","dataSources":[{"sourceId":101849,"databundleVersionId":13093295,"sourceType":"competition"},{"sourceId":9629432,"sourceType":"datasetVersion","datasetId":5846888},{"sourceId":12685517,"sourceType":"datasetVersion","datasetId":8008661},{"sourceId":247163208,"sourceType":"kernelVersion"}],"dockerImageVersionId":31040,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"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\n\nimport itertools\nfrom astropy.stats import sigma_clip\nimport joblib","metadata":{"trusted":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2025-08-06T07:18:23.21388Z","iopub.execute_input":"2025-08-06T07:18:23.214123Z","iopub.status.idle":"2025-08-06T07:18:29.085737Z","shell.execute_reply.started":"2025-08-06T07:18:23.214103Z","shell.execute_reply":"2025-08-06T07:18:29.085051Z"},"jupyter":{"source_hidden":true}},"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-08-06T07:18:29.087015Z","iopub.execute_input":"2025-08-06T07:18:29.087516Z","iopub.status.idle":"2025-08-06T07:18:29.319594Z","shell.execute_reply.started":"2025-08-06T07:18:29.087494Z","shell.execute_reply":"2025-08-06T07:18:29.318887Z"},"jupyter":{"source_hidden":true}},"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-08-06T07:18:29.32029Z","iopub.execute_input":"2025-08-06T07:18:29.320513Z","iopub.status.idle":"2025-08-06T07:18:29.324974Z","shell.execute_reply.started":"2025-08-06T07:18:29.320496Z","shell.execute_reply":"2025-08-06T07:18:29.324285Z"},"jupyter":{"source_hidden":true}},"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-08-06T07:18:29.326633Z","iopub.execute_input":"2025-08-06T07:18:29.326838Z","iopub.status.idle":"2025-08-06T07:18:30.786244Z","shell.execute_reply.started":"2025-08-06T07:18:29.326813Z","shell.execute_reply":"2025-08-06T07:18:30.785647Z"},"jupyter":{"source_hidden":true}},"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    # 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_visualized =  apply_adc_correction(signal_data,'FGS1')\ntotal_flux = np.sum(signal_data_visualized, axis=(1, 2))  \nprint(signal_data_visualized.shape, total_flux.shape)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:18:30.786872Z","iopub.execute_input":"2025-08-06T07:18:30.787131Z","iopub.status.idle":"2025-08-06T07:18:31.174941Z","shell.execute_reply.started":"2025-08-06T07:18:30.787112Z","shell.execute_reply":"2025-08-06T07:18:31.174063Z"},"jupyter":{"source_hidden":true}},"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_visualized, total_flux, planet_id)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:18:31.175725Z","iopub.execute_input":"2025-08-06T07:18:31.176004Z","iopub.status.idle":"2025-08-06T07:18:32.057974Z","shell.execute_reply.started":"2025-08-06T07:18:31.175977Z","shell.execute_reply":"2025-08-06T07:18:32.057162Z"},"jupyter":{"source_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Calibration - Tanh","metadata":{}},{"cell_type":"code","source":"!pip install -q --no-index --find-links=/kaggle/input/ariel-2024-pqdm pqdm\nimport pandas as pd\nimport numpy as np\nimport pandas.api.types\nimport scipy.stats\n\nfrom tqdm import tqdm\nfrom pqdm.threads import pqdm\nimport itertools\nimport pickle\n\nfrom scipy.optimize import minimize\nfrom sklearn.metrics import mean_squared_error\n\nimport plotly.express as px\n\nfrom astropy.stats import sigma_clip\nfrom scipy.signal import savgol_filter\nclass Config:\n    DATA_PATH = '/kaggle/input/ariel-data-challenge-2025'\n    DATASET = \"train\" # non -ML thì chỉ cần dùng cho test là được, còn bài này cần extraction hết\n\n    SCALE = 0.93960\n    SIGMA = 0.0009\n    \n    CUT_INF = 0\n    CUT_SUP = 356\n    \n    SENSOR_CONFIG = {\n        \"AIRS-CH0\": {\n            \"raw_shape\": [11250, 32, 356],\n            \"calibrated_shape\": [1, 32, CUT_SUP - CUT_INF],\n            \"linear_corr_shape\": (6, 32, 356),\n            \"dt_pattern\": (0.1, 4.5), \n            \"binning\": 30\n        },\n        \"FGS1\": {\n            \"raw_shape\": [135000, 32, 32],\n            \"calibrated_shape\": [1, 32, 32],\n            \"linear_corr_shape\": (6, 32, 32),\n            \"dt_pattern\": (0.1, 0.1),\n            \"binning\": 30 * 12\n        }\n    }\n    \n    MODEL_PHASE_DETECTION_SLICE = slice(30, 140)\n    MODEL_OPTIMIZATION_DELTA = 7\n    MODEL_POLYNOMIAL_DEGREE = 3\n    @classmethod\n    def set_dataset(cls, new_dataset: str):\n        cls.DATASET = new_dataset\n    \n    N_JOBS = 4\nclass SignalProcessor:\n    def __init__(self, config):\n        self.cfg = config\n        self.adc_info = pd.read_csv(f\"{self.cfg.DATA_PATH}/adc_info.csv\")\n        self.planet_ids = pd.read_csv(f'{self.cfg.DATA_PATH}/{self.cfg.DATASET}_star_info.csv', index_col='planet_id').index.astype(int)\n\n        \n    def set_dataset(self, new_dataset: str):\n        \"\"\"\n        Cập nhật dataset (train/test) và reload metadata tương ứng.\n        \"\"\"\n        self.cfg.set_dataset(new_dataset)\n        \n    def _apply_linear_corr(self, linear_corr, signal):\n        linear_corr_flipped = np.flip(linear_corr, axis=0)\n        corrected_signal = signal.copy()\n        \n        for x, y in itertools.product(range(signal.shape[1]), range(signal.shape[2])):\n            poly = np.poly1d(linear_corr_flipped[:, x, y])\n            corrected_signal[:, x, y] = poly(corrected_signal[:, x, y])\n        #print(f\"corrected_signal\\n{corrected_signal.shape}\")\n        return corrected_signal\n\n    def _calibrate_single_signal(self, planet_id, sensor):\n        sensor_cfg = self.cfg.SENSOR_CONFIG[sensor]\n\n        signal = pd.read_parquet(f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_signal_0.parquet\").to_numpy()\n        dark = pd.read_parquet(f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/dark.parquet\").to_numpy()\n        dead = pd.read_parquet(f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/dead.parquet\").to_numpy()\n        flat = pd.read_parquet(f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/flat.parquet\").to_numpy()\n        linear_corr = pd.read_parquet(f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/linear_corr.parquet\").values.astype(np.float64).reshape(sensor_cfg[\"linear_corr_shape\"])\n\n        signal = signal.reshape(sensor_cfg[\"raw_shape\"])\n        \n        #print(f\"\\nsignal.shape:\\n{signal.shape}\")\n        \n        gain = self.adc_info[f\"{sensor}_adc_gain\"].iloc[0]\n        offset = self.adc_info[f\"{sensor}_adc_offset\"].iloc[0]\n        signal = signal / gain + offset\n\n        hot = sigma_clip(dark, sigma=5, maxiters=5).mask\n\n        if sensor == \"AIRS-CH0\":\n            signal = signal[:, :, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            linear_corr = linear_corr[:, :, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            dark = dark[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            dead = dead[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            flat = flat[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            hot = hot[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n        \n        base_dt, increment = sensor_cfg[\"dt_pattern\"]\n        dt = np.ones(len(signal)) * base_dt\n        dt[1::2] += increment\n        \n        signal = signal.clip(0)\n        signal = self._apply_linear_corr(linear_corr, signal)\n        signal -= dark * dt[:, np.newaxis, np.newaxis]\n        \n        flat = flat.reshape(sensor_cfg[\"calibrated_shape\"])\n        flat[dead.reshape(sensor_cfg[\"calibrated_shape\"])] = np.nan\n        flat[hot.reshape(sensor_cfg[\"calibrated_shape\"])] = np.nan\n        \n        signal = signal / flat\n        #print(signal.shape)\n        return signal\n\n\n    def _process_planet_sensor(self, args):\n        planet_id, sensor = args['planet_id'], args['sensor']\n        calibrated = self._calibrate_single_signal(planet_id, sensor)\n        #print(f\"calibrated:{calibrated.shape}\")\n        return calibrated\n\n    def process_all_data(self):\n        args_fgs1 = [dict(planet_id=planet_id, sensor=\"FGS1\") for planet_id in self.planet_ids]\n        preprocessed_fgs1 = pqdm(args_fgs1, self._process_planet_sensor, n_jobs=self.cfg.N_JOBS)\n\n        args_airs_ch0 = [dict(planet_id=planet_id, sensor=\"AIRS-CH0\") for planet_id in self.planet_ids]\n        preprocessed_airs_ch0 = pqdm(args_airs_ch0, self._process_planet_sensor, n_jobs=self.cfg.N_JOBS)\n        \n        return preprocessed_fgs1[0],preprocessed_airs_ch0[0]\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:18:33.740483Z","iopub.execute_input":"2025-08-06T07:18:33.740762Z","iopub.status.idle":"2025-08-06T07:18:37.858443Z","shell.execute_reply.started":"2025-08-06T07:18:33.740736Z","shell.execute_reply":"2025-08-06T07:18:37.857707Z"},"jupyter":{"source_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Extraction function","metadata":{}},{"cell_type":"code","source":"%%time\n# =============================\n## Global Flux Features\n# =============================\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# =============================\n## Rolling Statistics Features\n# =============================\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# =============================\n## Transit Detection Features\n# =============================\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# =============================\n## Frequency Domain Features\n# =============================\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# =============================\n## Spatial Features\n# =============================\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.nansum(frame)\n        \n        if total_spatial_flux > 0:\n            centroid_x = np.nansum(x_indices * frame) / total_spatial_flux\n            centroid_y = np.nansum(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.nansum(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.nanmean(frame)\n        features[f'frame{i}_spatial_std'] = np.nanstd(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    centroids_x_movements = np.array(centroids_x)\n    centroids_y_movements = np.array(centroids_y)\n    movement = np.sqrt(np.diff(centroids_x_movements)**2 + np.diff(centroids_y_movements)**2)\n    # Movement analysis\n    features['centroid_x_range'] = np.nanmax(centroids_x) - np.nanmin(centroids_x)\n    features['centroid_y_range'] = np.nanmax(centroids_y) - np.nanmin(centroids_y)\n    features['centroid_total_movement'] = np.nansum(movement)  \n    features['concentration_range'] = np.nanmax(concentrations) - np.nanmin(concentrations)\n    features['concentration_std'] = np.nanstd(concentrations)\n    \n    return features\n\n# =============================\n## Gradient & Change Features\n# =============================\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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:18:37.859728Z","iopub.execute_input":"2025-08-06T07:18:37.860866Z","iopub.status.idle":"2025-08-06T07:18:37.882113Z","shell.execute_reply.started":"2025-08-06T07:18:37.86083Z","shell.execute_reply":"2025-08-06T07:18:37.881386Z"},"jupyter":{"source_hidden":true}},"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, planet_id,verbose=True):\n\n    n_frames = signal_data.shape[0]\n    \n    args_fgs1 = dict(planet_id=planet_id,sensor=\"FGS1\")\n    signal_data_fgs1 = signal_processor._process_planet_sensor(args_fgs1)\n    \n    # args_airs_ch0 = dict(planet_id=planet_id,sensor=\"AIRS-CH0\")\n    # signal_data_airs = signal_processor._process_planet_sensor(args_airs_ch0)\n    \n    total_flux = np.nansum(signal_data_fgs1, 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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:18:45.461124Z","iopub.execute_input":"2025-08-06T07:18:45.461382Z","iopub.status.idle":"2025-08-06T07:18:45.467066Z","shell.execute_reply.started":"2025-08-06T07:18:45.461363Z","shell.execute_reply":"2025-08-06T07:18:45.466317Z"},"jupyter":{"source_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\"\"\"config = Config()\nsignal_processor = SignalProcessor(config)\nsignal_data_fgs1,signal_data_airs = signal_processor.process_all_data()\ntotal_flux = np.nansum(signal_data_fgs1, axis=(1, 2)) # mask hot, dead pixel so that contains nan value\nprint(signal_data_fgs1.shape, total_flux.shape)\nall_features = extract_enhanced_transit_features(signal_data_fgs1, verbose=True)\n\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-08-06T07:18:46.809166Z","iopub.execute_input":"2025-08-06T07:18:46.809454Z","iopub.status.idle":"2025-08-06T07:18:46.81477Z","shell.execute_reply.started":"2025-08-06T07:18:46.809433Z","shell.execute_reply":"2025-08-06T07:18:46.814188Z"},"jupyter":{"source_hidden":true}},"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\nfrom pprint import pprint\n\ndef prepare_data_with_enhanced_fgs1(signal_processor,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        signal_paths = sorted(glob.glob(f\"{data_path}{base_id}/FGS1_signal_*.parquet\"))\n        #print(signal_paths)\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,base_id,verbose=False)\n            signal_id = f\"{base_id}_{j}\"\n            features['planet_id'] = signal_id\n            fgs1_features.append(features)\n\n        #print(i,fgs1_features)\n        #return 0\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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:18:48.490415Z","iopub.execute_input":"2025-08-06T07:18:48.490691Z","iopub.status.idle":"2025-08-06T07:18:48.499775Z","shell.execute_reply.started":"2025-08-06T07:18:48.490669Z","shell.execute_reply":"2025-08-06T07:18:48.498988Z"},"jupyter":{"source_hidden":true}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Dataset (Feature) & Model","metadata":{}},{"cell_type":"markdown","source":"Cho LGBM của bạn, mỗi mẫu (row) đầu vào sẽ gồm:\n\n1. **Feature matrix X**\n\n   * **175 engineered features** từ mỗi tín hiệu FGS1, gồm:\n\n     * 19 Global Flux Statistics\n     * 54 Rolling Statistics\n     * 20 Transit Detection\n     * 16 Frequency Domain\n     * 30 Spatial\n     * 36 Gradient & Change\n   * **8 star\\_info features**:\n\n     * Rs (stellar radius)\n     * Ms (stellar mass)\n     * Ts (stellar temperature)\n     * Mp (planet mass)\n     * e (eccentricity)\n     * P (orbital period)\n     * sma (semi-major axis)\n     * i (inclination)\n\n   Tổng input dimension cho mỗi row: **175 + 8 = 183**.\n\n2. **Target vector y**\n\n   * Mỗi model LGBM tương ứng một bước sóng (1 trong 283), nên\n\n     * với model thứ k (k=1…283), y là giá trị thực đo (flux hoặc residual) tại bước sóng k cho mỗi tín hiệu.\n   * Kết quả sau predict sẽ được lấy trung bình (average) qua tất cả các signals của cùng một planet.\n\n3. **Lưu ý về cấu trúc dữ liệu**\n\n   * Vì có nhiều signals per planet, bạn có nhiều hơn 1 row X–y cho cùng một planet.\n   * Dùng K-Fold CV: tách X, y thành các fold, train mỗi model trên 5 fold và ghi lại OOF/preds.\n\n---\n\nVí dụ khởi tạo X, y cho model bước sóng thứ k:\n\n```python\n# giả sử `features_df` đã chứa 175 feature + 8 star_info,\n# và `targets_df` có shape (n_rows, 283) chứa flux thực đo cho 283 bước sóng\n\nX = features_df.values              # shape (n_rows, 183)\ny = targets_df.iloc[:, k].values    # shape (n_rows,)\n```\n\nSau đó,train LGBMRegressor trên (X, y) rồi lặp cho k từ 0 đến 282.\n","metadata":{}},{"cell_type":"markdown","source":"## LGBM\nTrain 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_lgbm_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    }","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:18:51.965329Z","iopub.execute_input":"2025-08-06T07:18:51.966057Z","iopub.status.idle":"2025-08-06T07:18:51.974053Z","shell.execute_reply.started":"2025-08-06T07:18:51.966035Z","shell.execute_reply":"2025-08-06T07:18:51.973406Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## CatBoost","metadata":{}},{"cell_type":"code","source":"%%capture\n!pip install -q catboost","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:18:56.859634Z","iopub.execute_input":"2025-08-06T07:18:56.859908Z","iopub.status.idle":"2025-08-06T07:18:59.888034Z","shell.execute_reply.started":"2025-08-06T07:18:56.859885Z","shell.execute_reply":"2025-08-06T07:18:59.887039Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom catboost import CatBoostRegressor, Pool\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import mean_squared_error\ndef train_validate_catboost_multioutput_cv(X, y, n_splits=5, random_state=42):\n    \"\"\"\n    Args:\n        X: ndarray of shape (n_samples, n_features) # (,175+8)\n        y: ndarray of shape (n_samples,)\n        n_splits: int, number of folds\n        random_state: int\n\n    Returns:\n        oof_preds: ndarray of shape (n_samples,)\n        models: list of trained models\n        scores: list of RMSE per fold\n    \"\"\"\n    print(f\"Training with CV: {X.shape[1]} features → {y.shape[1]} wavelengths\")\n    kf = KFold(n_splits=n_splits, shuffle=True, random_state=random_state)\n    \n    oof_preds = np.zeros(len(X))\n    models = []\n    scores = []\n\n    for fold, (train_idx, val_idx) in enumerate(kf.split(X)):\n        X_train, X_val = X[train_idx], X[val_idx]\n        y_train, y_val = y[train_idx], y[val_idx]\n\n        model = CatBoostRegressor(\n            iterations=1000,\n            learning_rate=0.05, # mặc định của CatBoost là 0.03\n            depth=6,\n            loss_function='RMSE',\n            eval_metric='RMSE',\n            verbose=0,\n            random_seed=random_state\n        )\n\n        model.fit(X_train,\n                  y_train,\n                  eval_set=(X_val, y_val),\n                 # early_stopping_rounds=50\n                 )\n\n        val_preds = model.predict(X_val)\n        oof_preds[val_idx] = val_preds\n\n        rmse = mean_squared_error(y_val, val_preds, squared=False)\n        print(f\"Fold {fold+1} RMSE: {rmse:.5f}\")\n\n        scores.append(rmse)\n        models.append(model)\n\n    avg_rmse = np.mean(scores)\n    print(f\"Overall OOF RMSE: {avg_rmse:.5f}\")\n\n    return oof_preds, models, scores\n    \n# X: (1210, 183)\n# y: (1210, 283)\n\n#X_train = X\n#y_train = y.iloc[:, 0].values   # bước sóng 0\n\n#oof_preds, models, scores = train_catboost_for_wavelength(X_train, y_train)\n\"\"\"Fold 1 RMSE: 0.00151\nFold 2 RMSE: 0.00138\nFold 3 RMSE: 0.00092\nFold 4 RMSE: 0.00161\nFold 5 RMSE: 0.00083\nOverall OOF RMSE: 0.00125\"\"\"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:18:59.889921Z","iopub.execute_input":"2025-08-06T07:18:59.890158Z","iopub.status.idle":"2025-08-06T07:19:00.325433Z","shell.execute_reply.started":"2025-08-06T07:18:59.890134Z","shell.execute_reply":"2025-08-06T07:19:00.324657Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## XGBoost","metadata":{}},{"cell_type":"code","source":"from xgboost import XGBRegressor\nfrom sklearn.metrics import r2_score, mean_squared_error\nfrom sklearn.model_selection import KFold\nimport numpy as np\nfrom tqdm import tqdm\n\ndef train_validate_xgb_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\"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 = XGBRegressor(\n                n_estimators=500, # số lượng cây (DT) huấn luyện\n                learning_rate=0.05,\n                max_depth=8,\n                min_child_weight=20, # số lượng tối thiểu của tổng trọng số\n                subsample=0.8, # ratio data lấy mẫu để huấn luyện\n                colsample_bytree=0.8, # ratio features đc sử dụng cho mỗi cây\n                reg_alpha=0.1, # lasso để giảm complexity\n                reg_lambda=0.1, # rigde làm mượt trọng số\n                random_state=42,\n                verbosity=0\n            )\n            \n            # Lấy target cho step sóng i\n            y_single = y_train.iloc[:, i]\n            model.fit(X_train, y_single)\n            all_models[i].append(model)\n\n            # Dự đoán trên tập validation (OOF)\n            pred_val = model.predict(X_val)\n            oof_preds[val_idx, i] = pred_val\n\n            # Tính residuals trên tập train để ước lượng độ lệch\n            pred_train = model.predict(X_train)\n            residuals = y_single - pred_train\n            residuals_all[i].extend(residuals.tolist())\n\n            # Sigma tạm thời dùng độ lệch chuẩn của residuals\n            sigma = np.std(residuals)\n            oof_sigmas[val_idx,i] = max(sigma, 1e-6)\n            \n    target_uncertainties = [max(np.std(res),1e-6) for res in residuals_all] \n\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            ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:19:00.326309Z","iopub.execute_input":"2025-08-06T07:19:00.327388Z","iopub.status.idle":"2025-08-06T07:19:00.337572Z","shell.execute_reply.started":"2025-08-06T07:19:00.327369Z","shell.execute_reply":"2025-08-06T07:19:00.336743Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"!pip install -q optuna\n\nimport optuna\nfrom xgboost import XGBRegressor\nfrom sklearn.model_selection import cross_val_score, KFold\nimport random\n\ndef objective(trial):\n    # Hyperparameter space\n    params = {\n        'n_estimators': trial.suggest_int('n_estimators', 100, 1000),\n        'learning_rate': trial.suggest_float('learning_rate', 1e-3, 0.3, log=True),\n        'max_depth': trial.suggest_int('max_depth', 3, 10),\n        'min_child_weight': trial.suggest_int('min_child_weight', 1, 30),\n        'subsample': trial.suggest_float('subsample', 0.5, 1.0),\n        'colsample_bytree': trial.suggest_float('colsample_bytree', 0.5, 1.0),\n        'reg_alpha': trial.suggest_float('reg_alpha', 0, 1.0),\n        'reg_lambda': trial.suggest_float('reg_lambda', 0, 1.0),\n        'random_state': 42,\n        'verbosity': 0,\n    }\n    \"\"\"model = XGBRegressor(**params)\n    \n        # Dùng wavelength đầu tiên (ví dụ y[:, 0]) để tune\n        scores = cross_val_score(model, X, y.iloc[:, 0], \n                                 cv=KFold(n_splits=3, shuffle=True, random_state=42),\n                                 scoring='neg_root_mean_squared_error')\n        \n        return -scores.mean()  # minimize RMSE\"\"\"\n    \n    # Chọn ngẫu nhiên k bước sóng khác nhau\n    k = 50\n    total_wavelengths = y.shape[1]\n    random.seed(42)  # cố định seed để tái lập được\n    selected_indices = random.sample(range(total_wavelengths), k)\n\n    rmse_list = []\n\n    for i in selected_indices:\n        y_target = y.iloc[:, i]\n        model = XGBRegressor(**params)\n        scores = cross_val_score(model, X, y_target,\n                                 cv=KFold(n_splits=3, shuffle=True, random_state=42),\n                                 scoring='neg_root_mean_squared_error')\n        rmse_list.append(-scores.mean())  # convert to positive RMSE\n\n    return np.mean(rmse_list)  # minimize average RMSE\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:46:35.686944Z","iopub.execute_input":"2025-08-06T07:46:35.687506Z","iopub.status.idle":"2025-08-06T07:46:38.773328Z","shell.execute_reply.started":"2025-08-06T07:46:35.687481Z","shell.execute_reply":"2025-08-06T07:46:38.772274Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Tabnet","metadata":{}},{"cell_type":"code","source":"!pip install --no-deps /kaggle/input/support-neurips/pytorch_tabnet-4.1.0-py3-none-any.whl\nfrom pytorch_tabnet.tab_model import TabNetClassifier, TabNetRegressor\nprint(\"TabNet đã được import thành công!\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:29:10.197821Z","iopub.execute_input":"2025-08-06T07:29:10.198377Z","iopub.status.idle":"2025-08-06T07:29:11.70086Z","shell.execute_reply.started":"2025-08-06T07:29:10.198351Z","shell.execute_reply":"2025-08-06T07:29:11.70009Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 3.2 TabNet-based Regressor\n\nWe experiment with TabNet (Arik & Pfister, 2021) as an alternative to tree-based models for multi-output flux prediction.\n\n**Data Preparation**  \nThe input matrix includes 183 features per signal (175 FGS1 + 8 planetary info). Targets are 283 wavelength channels. Features are standardized before training.\n\n**Training Configuration**  \nWe apply K-Fold cross-validation (k=5) and train TabNet models with the following parameters:\n\n- `n_d`/`n_a`: 32\n- `n_steps`: 5\n- `mask_type`: sparsemax\n- `max_epochs`: 200, early stopping with patience=20\n\n**Evaluation Metrics**  \nWe report `R²` and `RMSE` for each wavelength channel. TabNet performance is compared against XGBoost and LGBM as baselines.\n\n**Observations**  \n- TabNet performs competitively on wavelength groups with more variance.\n- However, it is more computationally expensive.\n- It offers model interpretability through attention masks (optional).\n\n","metadata":{}},{"cell_type":"code","source":"from pytorch_tabnet.tab_model import TabNetRegressor\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import r2_score, mean_squared_error\nfrom sklearn.preprocessing import StandardScaler\nimport numpy as np\nimport torch\n\ndef train_validate_tabnet_multioutput_cv(X_df, y_df, batch_size,n_splits=5,max_epochs=200, seed=42):\n    np.random.seed(seed)\n    torch.manual_seed(seed)\n\n#    X_np = X_df.values.astype(np.float32)\n#    y_np = y_df.values.astype(np.float32) \n\n    # scaling features\n    scaler=StandardScaler()\n    X = scaler.fit_transform(X_df.values) # type np.array\n    y = y_df.values # type np.array\n\n    n_targets = y.shape[1] # dimension: 283 \n    kf = KFold(n_splits=n_splits,shuffle=True,random_state=seed)\n\n    # train loop\n    all_models = [[None]*n_splits for _ in range(n_targets)]\n    oof_preds = np.zeros_like(y, dtype=float)\n    oof_sigmas = np.zeros_like(y, dtype=float)\n    residuals_all = [[] for _ in range(n_targets)]\n\n    for fold, (train_idx,val_idx) in enumerate(kf.split(X)):\n        X_train, X_val = X[train_idx],X[val_idx]\n        y_train, y_val = y[train_idx],y[val_idx]\n        #print(X_train.shape)\n        \n        #train model_i cho mỗi wave_length_i vì TabNetRegressor chỉ generate single-output. \n        for i in range(n_targets):\n            model = TabNetRegressor(\n                n_d=16, #decision\n                n_a=16, #attention\n                n_steps=5, #step acculative feature transformer (attention + decision) => 1 bước chọn 1 feature khác\n                gamma=1.5,\n                lambda_sparse=1e-5, # L1 penalty lên attention mask ? chưa hiểu lắm\n                optimizer_fn=torch.optim.Adam,\n                optimizer_params=dict(lr=2e-2), # cần xem xét\n                mask_type='sparsemax',\n                seed=seed,\n                verbose=0\n            )\n            model.fit(\n                X_train, y_train[:, i].reshape(-1,1),\n                eval_set=[(X_val, y_val[:, i].reshape(-1,1))],\n                eval_metric=['rmse'],\n                max_epochs=max_epochs,\n                patience=10, # default\n                batch_size=batch_size,\n                virtual_batch_size=batch_size//4,\n                #num_workers=0\n            )\n            # lưu model\n            all_models[i][fold] = model\n            \n            # OOF predict\n            pred_val = model.predict(X_val).flatten()\n            oof_preds[val_idx, i] = pred_val\n\n            # train residual để ước lượng uncertainty\n            pred_train = model.predict(X_train).flatten()\n            res = y_train[:, i] - pred_train\n            residuals_all[i].extend(res.tolist())\n            oof_sigmas[val_idx, i] = max(np.std(res), 1e-6)\n\n    # 4. Tính target_uncertainties\n    target_uncertainties = [max(np.std(res), 1e-6) for res in residuals_all]\n    \n    # 5. Tính metrics\n    r2_scores = [r2_score(y[:, i], oof_preds[:, i]) for i in range(n_targets)]\n    rmses = [mean_squared_error(y[:, i], oof_preds[:, i], squared=False) for i in range(n_targets)]\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    # 6. Trả về cấu trúc giống XGB\n    return {\n        'models': all_models,\n        'feature_columns': X_df.columns.tolist(),\n        'target_columns': y_df.columns.tolist(),\n        'oof_predictions': oof_preds,\n        'oof_uncertainties': oof_sigmas,\n        'y_true': y,\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    }","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:29:25.346822Z","iopub.execute_input":"2025-08-06T07:29:25.34802Z","iopub.status.idle":"2025-08-06T07:29:25.359379Z","shell.execute_reply.started":"2025-08-06T07:29:25.347991Z","shell.execute_reply":"2025-08-06T07:29:25.358587Z"}},"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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:39:16.466156Z","iopub.execute_input":"2025-08-06T07:39:16.466734Z","iopub.status.idle":"2025-08-06T07:39:16.474133Z","shell.execute_reply.started":"2025-08-06T07:39:16.466709Z","shell.execute_reply":"2025-08-06T07:39:16.473323Z"}},"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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:39:17.685624Z","iopub.execute_input":"2025-08-06T07:39:17.686335Z","iopub.status.idle":"2025-08-06T07:39:17.692748Z","shell.execute_reply.started":"2025-08-06T07:39:17.686302Z","shell.execute_reply":"2025-08-06T07:39:17.692015Z"}},"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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:39:17.913856Z","iopub.execute_input":"2025-08-06T07:39:17.914118Z","iopub.status.idle":"2025-08-06T07:39:17.927324Z","shell.execute_reply.started":"2025-08-06T07:39:17.914098Z","shell.execute_reply":"2025-08-06T07:39:17.92671Z"}},"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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:39:19.288016Z","iopub.execute_input":"2025-08-06T07:39:19.288301Z","iopub.status.idle":"2025-08-06T07:39:19.296597Z","shell.execute_reply.started":"2025-08-06T07:39:19.288251Z","shell.execute_reply":"2025-08-06T07:39:19.295847Z"}},"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(signal_processor,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        signal_processor,\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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:39:19.517138Z","iopub.execute_input":"2025-08-06T07:39:19.517649Z","iopub.status.idle":"2025-08-06T07:39:19.524381Z","shell.execute_reply.started":"2025-08-06T07:39:19.517626Z","shell.execute_reply":"2025-08-06T07:39:19.523666Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Running","metadata":{}},{"cell_type":"code","source":"config = Config()\nsignal_processor = SignalProcessor(config)\nsignal_processor.set_dataset('train')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:39:20.902488Z","iopub.execute_input":"2025-08-06T07:39:20.903085Z","iopub.status.idle":"2025-08-06T07:39:20.915072Z","shell.execute_reply.started":"2025-08-06T07:39:20.90306Z","shell.execute_reply":"2025-08-06T07:39:20.914296Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def _is_prepare_data(is_prepare_data: bool = True):\n    if is_prepare_data: \n        X, y = prepare_data_with_enhanced_fgs1(signal_processor,train_df, train_star_info_df, \"/kaggle/input/ariel-data-challenge-2025/train/\")\n        # Lưu file X,Y để load tiết kiệm thời gian\n        joblib.dump(X, \"X_enhanced_fgs1.pkl\")\n        joblib.dump(y, \"y_enhanced_fgs1.pkl\")\n    else:\n        # Đọc lại sau này\n        X = joblib.load(\"/kaggle/input/support-neurips/X_enhanced_fgs1.pkl\")\n        y = joblib.load(\"/kaggle/input/support-neurips/y_enhanced_fgs1.pkl\")\n    return X,y\nX,y = _is_prepare_data(False)\nprint(f\"✅ Final shapes — X: {X.shape}, y: {y.shape}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:46:39.76868Z","iopub.execute_input":"2025-08-06T07:46:39.769486Z","iopub.status.idle":"2025-08-06T07:46:39.783472Z","shell.execute_reply.started":"2025-08-06T07:46:39.769452Z","shell.execute_reply":"2025-08-06T07:46:39.78266Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"study = optuna.create_study(direction='minimize')\nstudy.optimize(objective, n_trials=30)  # thử 30 cấu hình\nprint(\"✅ Best trial:\")\nprint(study.best_trial)\n\nprint(\"\\n📌 Best params:\")\nfor key, value in study.best_trial.params.items():\n    print(f\"{key}: {value}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:47:09.346834Z","iopub.execute_input":"2025-08-06T07:47:09.347113Z","iopub.status.idle":"2025-08-06T07:47:17.404675Z","shell.execute_reply.started":"2025-08-06T07:47:09.347091Z","shell.execute_reply":"2025-08-06T07:47:17.404116Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## EDA","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom sklearn.preprocessing import StandardScaler\n\n# 1. Load your dataframes\n# features_df: (n_rows, 183), targets_df: (n_rows, 283)\n# Assume they've been loaded already.\nfeatures_df = X\ntargets_df = y\ndesc_all = features_df.describe()\nprint(\"=== Descriptive stats for all features ===\")\ndisplay(desc_all.head())\n\ncols = ['Rs','Ms','Ts','Mp','e','P']\n\nfig, axes = plt.subplots(2, 3, figsize=(12, 6))\nfor ax, col in zip(axes.flatten(), cols):\n    ax.boxplot(features_df[col].dropna())\n    ax.set_title(f\"{col}\")\n    ax.set_ylabel(col)\nplt.suptitle(\"Boxplots of star_info features\")\nplt.tight_layout(rect=[0, 0.03, 1, 0.95])\nplt.show()\n\n# 2. Quick EDA: histograms of a few features\n\nfig, axes = plt.subplots(2, 3, figsize=(12, 6))\nfor ax, col in zip(axes.flat, features_df.columns[:6]):\n    ax.hist(features_df[col], bins=30)\n    ax.set_title(col)\nplt.suptitle(\"Before scaling\"); plt.tight_layout(); plt.show()\n\n\n# 3. Scale all features\nscaler = StandardScaler()\nX_scaled = scaler.fit_transform(features_df)\n\n# 4. Plot same features after scaling\nfig, axes = plt.subplots(2, 3, figsize=(12, 6))\nfor ax, idx in zip(axes.flat, range(6)):\n    ax.hist(X_scaled[:, idx], bins=30)\n    ax.set_title(features_df.columns[idx])\nplt.suptitle(\"After StandardScaler\"); plt.tight_layout(); plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:47:33.763581Z","iopub.execute_input":"2025-08-06T07:47:33.763902Z","iopub.status.idle":"2025-08-06T07:47:36.537867Z","shell.execute_reply.started":"2025-08-06T07:47:33.763875Z","shell.execute_reply":"2025-08-06T07:47:36.537069Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%capture\n#trained_models_cv = train_validate_lgbm_multioutput_cv(X.iloc[:10], y.iloc[:10])\n\n#trained_models_cv = train_validate_xgb_multioutput_cv(X, y)\n\ntrained_models_cv = train_validate_tabnet_multioutput_cv(X,\n                                                         y,\n                                                         batch_size=32,\n                                                         n_splits=5,\n                                                         max_epochs=200, \n                                                         seed=42)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:32:24.512945Z","iopub.execute_input":"2025-08-06T07:32:24.513957Z","iopub.status.idle":"2025-08-06T07:38:02.595637Z","shell.execute_reply.started":"2025-08-06T07:32:24.513915Z","shell.execute_reply":"2025-08-06T07:38:02.594953Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# In ra tóm tắt\nprint(\"▶️ Mean R²:\", trained_models_cv['cv_metrics']['mean_r2'])\nprint(\"▶️ Mean RMSE:\", trained_models_cv['cv_metrics']['mean_rmse'])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:38:02.597061Z","iopub.execute_input":"2025-08-06T07:38:02.597357Z","iopub.status.idle":"2025-08-06T07:38:02.60143Z","shell.execute_reply.started":"2025-08-06T07:38:02.59733Z","shell.execute_reply":"2025-08-06T07:38:02.600698Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"feature_importance_df = analyze_feature_importance(trained_models_cv, X)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:39:24.718998Z","iopub.execute_input":"2025-08-06T07:39:24.7193Z","iopub.status.idle":"2025-08-06T07:39:24.746291Z","shell.execute_reply.started":"2025-08-06T07:39:24.719251Z","shell.execute_reply":"2025-08-06T07:39:24.745631Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 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}\")\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-08-06T07:39:39.030532Z","iopub.execute_input":"2025-08-06T07:39:39.030793Z","iopub.status.idle":"2025-08-06T07:39:43.972558Z","shell.execute_reply.started":"2025-08-06T07:39:39.030775Z","shell.execute_reply":"2025-08-06T07:39:43.97182Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"signal_processor.set_dataset('test')\nsubmission = create_submission(\n    signal_processor,\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-08-06T07:38:02.632563Z","iopub.status.idle":"2025-08-06T07:38:02.632894Z","shell.execute_reply.started":"2025-08-06T07:38:02.632733Z","shell.execute_reply":"2025-08-06T07:38:02.632748Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-06T07:38:02.634531Z","iopub.status.idle":"2025-08-06T07:38:02.63486Z","shell.execute_reply.started":"2025-08-06T07:38:02.634672Z","shell.execute_reply":"2025-08-06T07:38:02.634689Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}