{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.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":70367,"databundleVersionId":9188054,"sourceType":"competition"}],"dockerImageVersionId":30746,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"References\n\n1. [Why Did Calibration Lead to a Lower Public Score When Combining Two Kaggle Notebooks?](https://www.kaggle.com/competitions/ariel-data-challenge-2024/discussion/530472)\n\n2. [Fork of NeurIPS Ariel 2024 - Starter 5be123](https://www.kaggle.com/code/regisvargas/fork-of-neurips-ariel-2024-starter-5be123)\n\n3. [NeurIPS Ariel 2024 - Starter withdifferentparametr](https://www.kaggle.com/code/bingyuniu/neurips-ariel-2024-starter-withdifferentparametr)\n\n4. [[UPDATE]Calibrating and Binning Astronomical Data](https://www.kaggle.com/code/gordonyip/update-calibrating-and-binning-astronomical-data)\n\n5. [[UPDATE]Calibrating and Binning Astronomical Data (copy)](https://www.kaggle.com/code/aaronjday/update-calibrating-and-binning-astronomical-data)\n\n6. [ariel_only_correlation](https://www.kaggle.com/code/sergeifironov/ariel-only-correlation)","metadata":{}},{"cell_type":"markdown","source":"# Initialization\n\nThis competition seems requires strong scientific background and I had lot of confusion during EDA process. Therefore, I just build a simple starter for future coding.\n\n## Load Library","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns\nimport pandas as pd\nimport polars as pl\nimport numpy as np\nimport torch\nfrom scipy import signal\nimport random\nfrom sklearn.metrics import mean_squared_error\n# import torch.nn as nn\n# import torch.optim as optim\n# from torch.utils.data import DataLoader, Dataset\nfrom tqdm import tqdm\nimport pickle\nimport time\nimport os\nimport pickle\nimport seaborn as sns\nimport scipy.stats\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import r2_score, mean_squared_error\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.cross_decomposition import CCA\nfrom scipy.stats import shapiro, anderson\nfrom statsmodels.tsa.stattools import acf","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load Meta-Data\nPATH = \"/kaggle/input/ariel-data-challenge-2024\"\ntrain_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_adc_info.csv', \n                             index_col='planet_id')\ntrain_labels = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_labels.csv',\n                           index_col='planet_id')\nwavelengths = pd.read_csv(f'{PATH}/wavelengths.csv')\naxis_info = pd.read_parquet(os.path.join(PATH,'axis_info.parquet'))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_adc_info['AIRS-CH0_adc_gain'].loc[785834]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Pre-Processing\n## Load Functions","metadata":{}},{"cell_type":"code","source":"import pywt\ndef wavelet_denoising(data, wavelet = 'db10', sigma=None):\n    \"\"\"Denoises the signal using SURE wavelet shrinkage with the specified wavelet.\"\"\"\n    # Criar uma cópia do array para evitar o erro de \"read-only\"\n    data = np.array(data, copy=True)\n    \n    # Decompose the signal using discrete wavelet transform\n    coeffs = pywt.wavedec(data, wavelet)\n    \n    # Estimate noise level if sigma is not provided\n    if sigma is None:\n        # Using the Median Absolute Deviation (MAD) estimator for noise level\n        sigma = np.median(np.abs(coeffs[-1])) / 0.6745\n    # Apply thresholding (SURE or hard/soft)\n    threshold = sigma * np.sqrt(2 * np.log(len(data)))\n    new_coeffs = [pywt.threshold(c, threshold, mode='soft') for c in coeffs]\n    # Reconstruct the signal using the modified coefficients\n    return pywt.waverec(new_coeffs, wavelet)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import itertools\ndef apply_linear_corr(linear_corr,clean_signal):\n    linear_corr = np.flip(linear_corr, axis=0)\n    for x, y in itertools.product(\n                range(clean_signal.shape[1]), range(clean_signal.shape[2])\n            ):\n        poli = np.poly1d(linear_corr[:, x, y])\n        clean_signal[:, x, y] = poli(clean_signal[:, x, y])\n    return clean_signal","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from astropy.stats import sigma_clip\ndef mask_hot_dead(signal, dead, dark):\n    # Apply sigma_clip to the dark frame and convert the mask to boolean\n    hot = sigma_clip(dark, sigma=5, maxiters=5).mask.astype(bool)\n    # Ensure the dead mask is also in boolean format\n    dead = dead.astype(bool)\n    # Create a combined mask using bitwise OR\n    combined_mask = np.tile(hot, (signal.shape[0], 1, 1)) | np.tile(dead, (signal.shape[0], 1, 1))\n    # Calculate the mean of the unmasked values for each frame\n    mean_signal = np.mean(signal[~combined_mask])\n    # Replace masked values with the mean of the unmasked values\n    signal[combined_mask] = mean_signal\n    return signal","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def clean_dark(signal, dark, dt):\n    dark = np.tile(dark, (signal.shape[0], 1, 1))\n    signal -= dark* dt[:, np.newaxis, np.newaxis]\n    return signal","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_cds(signal):\n    cds = signal[1::2,:,:] - signal[::2,:,:]\n    return cds","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def remove_outliers(data):\n    q1 = np.percentile(data, 25)\n    q3 = np.percentile(data, 75)\n    iqr = q3 - q1\n    lower_bound = q1 - 1.5 * iqr\n    upper_bound = q3 + 1.5 * iqr\n    return data[(data >= lower_bound) & (data <= upper_bound)]\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def mean_without_outliers(data):\n    data_no_outliers = remove_outliers(data)\n    if len(data_no_outliers) > 0:\n        return data_no_outliers.mean()\n    else:\n        return np.nan  # Retorna NaN se todos os valores forem outliers","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def phase_detector(signal, window_size=3000, search_range_obscured=(5000, 25000), search_range_unobscured=(30000, 60000)):\n    \"\"\"\n    Detecta as fases obscured e unobscured de um sinal.\n    \n    Parâmetros:\n    - signal: array, o sinal no qual detectar as fases.\n    - window_size: int, tamanho da janela de busca para identificar mudanças.\n    - search_range_obscured: tuple, intervalo de busca para a fase obscured.\n    - search_range_unobscured: tuple, intervalo de busca para a fase unobscured.\n    \n    Retorna:\n    - phase1: índice do final da fase \"obscured\".\n    - phase2: índice do início da fase \"unobscured\".\n    \"\"\"\n    phase1, phase2 = None, None\n    best_drop_obscured = 0\n    best_drop_unobscured = 0\n    \n    # Detecta a maior queda no sinal dentro do intervalo da fase obscured\n    for i in range(search_range_obscured[0], search_range_obscured[1]):\n        # Verifica se a janela não extrapola os limites do array\n        if i + window_size <= len(signal):\n            t1 = signal[i:i + window_size].max() - signal[i:i + window_size].min()\n            if t1 > best_drop_obscured:\n                phase1 = i + window_size\n                best_drop_obscured = t1\n    \n    # Detecta a maior queda no sinal dentro do intervalo da fase unobscured\n    for i in range(search_range_unobscured[0], search_range_unobscured[1]):\n        # Verifica se a janela não extrapola os limites do array\n        if i + window_size <= len(signal):\n            t1 = signal[i:i + window_size].max() - signal[i:i + window_size].min()\n            if t1 > best_drop_unobscured:\n                phase2 = i + window_size\n                best_drop_unobscured = t1\n    return phase1, phase2\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%writefile utils.py\nimport pandas as pd\nimport polars as pl\nimport numpy as np\nfrom tqdm import tqdm\nimport pickle\nPATH = \"/kaggle/input/ariel-data-challenge-2024\"\ndef load_signal_data(planet_id, dataset, instrument, img_size):\n    file_path = f'{PATH}/{dataset}/{planet_id}/{instrument}_signal.parquet'\n    signal = pd.read_parquet(file_path)\n    if instrument==\"AIRS-CH0\":\n        signal = signal.values.astype(np.float64).reshape((signal.shape[0], 32, 356))\n    else:\n        signal = signal.values.astype(np.float64).reshape((signal.shape[0], 32, 32))\n    if dataset == 'train':\n        gain = train_adc_info[instrument+'_adc_gain'].loc[planet_id]\n        offset = train_adc_info[instrument+'_adc_offset'].loc[planet_id]\n        signal = ADC_convert(signal, gain, offset)\n    else:\n        gain = test_adc_info[instrument+'_adc_gain'].loc[planet_id]\n        offset = test_adc_info[instrument+'_adc_offset'].loc[planet_id]\n        signal = ADC_convert(signal, gain, offset)\n    signal = signal.reshape(signal.shape[0], signal.shape[1] * signal.shape[2])\n    mean_signal = signal.mean(axis=1)\n    mean_signal = mean_signal / np.linalg.norm(mean_signal)\n    net_signal = mean_signal[1::2] - mean_signal[0::2]\n    return wavelet_denoising(net_signal)\ndef ADC_convert(signal, gain, offset):\n    signal = signal.astype(np.float64)\n    signal /= gain\n    signal += offset\n    return signal\ndef read_and_preprocess(dataset, planet_ids, instrument = \"AIRS-CH0\"):\n    \"\"\"Read the files for all planet_ids and extract the time series.\n    Parameters\n    dataset: 'train' or 'test'\n    planet_ids: list of planet ids\n    instrument: the instrument of observation, 'AIRS-CH0' or 'FGS1', default to 'AIRS-CH0'\n    Returns\n    dataframe with one row per planet_id and 67500 values per row for FGS1 and 5624 for AIRS-CH0\n    \"\"\"\n    img_size = 1024 if instrument == \"FGS1\" else 32*356\n    column_num = 67500 if instrument == 'FGS1' else 5626\n    raw_train = np.full((len(planet_ids), column_num), np.nan, dtype=np.float32)\n    for i, planet_id in tqdm(list(enumerate(planet_ids))):\n        raw_train[i] = load_signal_data(planet_id, dataset, instrument, img_size)\n    return raw_train\ndef feature_engineering(f_raw, a_raw, adc_info, window_size=50, step_size=15, dataset='train'):\n    \"\"\"Create a dataframe with combined features from the raw data, including sliding window and time-series statistics.\n    \n    Parameters:\n    f_raw: ndarray of shape (n_planets, 67500)\n    a_raw: ndarray of shape (n_planets, 5625)\n    window_size: int, size of the sliding window for time-series statistics\n    step_size: int, step size for the sliding window\n    \n    Return value:\n    df: DataFrame of shape (n_planets, several features), valid_indices: list of indices after removing outliers\n    \"\"\"\n    # Para o f_raw (sinal \"f\")\n    f_obscured_phases = []\n    f_unobscured_phases = []\n    for signal in f_raw:\n        phase1, phase2 = phase_detector(signal)\n        if phase1 == None:\n            phase1 = 23500\n        if phase2 == None:\n            phase2 = 44000\n        # Calcula a média removendo outliers para a parte obscured\n        f_obscured_segment = signal[phase1:phase2]\n        f_obscured_phases.append(mean_without_outliers(f_obscured_segment))\n    \n        # Calcula a média removendo outliers para a parte unobscured\n        f_unobscured_segment1 = signal[:phase1]\n        f_unobscured_segment2 = signal[phase2:]\n        f_unobscured1 = mean_without_outliers(f_unobscured_segment1)\n        f_unobscured2 = mean_without_outliers(f_unobscured_segment2)\n    \n        f_unobscured_phases.append((f_unobscured1 + f_unobscured2) / 2)\n    f_obscured = np.array(f_obscured_phases)\n    f_unobscured = np.array(f_unobscured_phases)\n    f_relative_reduction = (f_unobscured - f_obscured) / f_unobscured\n    f_std_dev = f_raw.std(axis=1)\n    f_signal_to_noise = f_unobscured / f_std_dev\n    # Para o a_raw (sinal \"a\")\n    a_obscured_phases = []\n    a_unobscured_phases = []\n    for signal in a_raw:\n        phase1, phase2 = phase_detector(signal)\n        if phase1 == None:\n            phase1 = 1958\n        if phase2 == None:\n            phase2 = 3666    \n        # Calcula a média removendo outliers para a parte obscured\n        a_obscured_segment = signal[phase1:phase2]\n        a_obscured_phases.append(mean_without_outliers(a_obscured_segment))\n    \n        # Calcula a média removendo outliers para a parte unobscured\n        a_unobscured_segment1 = signal[:phase1]\n        a_unobscured_segment2 = signal[phase2:]\n        a_unobscured1 = mean_without_outliers(a_unobscured_segment1)\n        a_unobscured2 = mean_without_outliers(a_unobscured_segment2)\n    \n        a_unobscured_phases.append((a_unobscured1 + a_unobscured2) / 2)\n    a_obscured = np.array(a_obscured_phases)\n    a_unobscured = np.array(a_unobscured_phases)\n    a_relative_reduction = (a_unobscured - a_obscured) / a_unobscured\n    a_std_dev = a_raw.std(axis=1)\n    a_signal_to_noise = a_unobscured / a_std_dev\n    f_variance = f_raw.var(axis=1)\n    a_variance = a_raw.var(axis=1)\n    \n    f_skewness = pd.DataFrame(f_raw).skew(axis=1).values\n    a_skewness = pd.DataFrame(a_raw).skew(axis=1).values\n    f_kurtosis = pd.DataFrame(f_raw).kurtosis(axis=1).values\n    a_kurtosis = pd.DataFrame(a_raw).kurtosis(axis=1).values\n    \n    f_half_obscured1 = f_raw[:, 20500:23500].mean(axis=1)\n    f_half_obscured2 = f_raw[:, 44000:47000].mean(axis=1)\n    f_half_reduction1 = (f_unobscured - f_half_obscured1) / f_unobscured\n    f_half_reduction2 = (f_unobscured - f_half_obscured2) / f_unobscured\n    a_half_obscured1 = a_raw[:, 1708:1958].mean(axis=1)\n    a_half_obscured2 = a_raw[:, 3666:3916].mean(axis=1)\n    a_half_reduction1 = (a_unobscured - a_half_obscured1) / a_unobscured\n    a_half_reduction2 = (a_unobscured - a_half_obscured2) / a_unobscured\n    \n    # Sliding window features\n    def sliding_window_features(data, window_size, step_size):\n        features = []\n        max_index = data.shape[1]\n        for start in range(0, max_index - window_size + 1, step_size):\n            end = start + window_size\n            window = data[:, start:end]\n            features.append([\n                np.mean(window, axis=1),\n                np.std(window, axis=1),\n                np.min(window, axis=1),\n                np.max(window, axis=1)\n            ])\n        if features:\n            return np.vstack(features).T  # Stack vertically and transpose to get the correct shape\n        else:\n            return np.empty((data.shape[0], 0))  # Return empty array with correct shape\n    \n    f_sliding_features = sliding_window_features(f_raw, window_size, step_size)\n    a_sliding_features = sliding_window_features(a_raw, window_size, step_size)\n    \n    df = pd.DataFrame({\n        'f_relative_reduction': f_relative_reduction,\n        'f_signal_to_noise': f_signal_to_noise,\n        'f_variance': f_variance,\n        'f_skewness': f_skewness,\n        'f_kurtosis': f_kurtosis,\n        'a_relative_reduction': a_relative_reduction,\n        'a_signal_to_noise': a_signal_to_noise,\n        'a_variance': a_variance,\n        'a_skewness': a_skewness,\n        'a_kurtosis': a_kurtosis,\n        'f_half_reduction1': f_half_reduction1,\n        'f_half_reduction2': f_half_reduction2,\n        'a_half_reduction1': a_half_reduction1,\n        'a_half_reduction2': a_half_reduction2,\n        'a_f_mean_reduction': a_relative_reduction,\n        'a_f_mean_reduction_smooth': a_relative_reduction\n    })\n    def linearity_loss(smoothed_signal):\n        n = len(smoothed_signal)\n        x = np.arange(n)\n        p = np.polyfit(x, smoothed_signal, 1)\n        fitted_line = np.polyval(p, x)\n        mse_linear = np.mean((smoothed_signal - fitted_line) ** 2)\n        return mse_linear\n    def normality_loss(residuals):\n        _, p_value = shapiro(residuals)\n        penalidade_normalidade = max(0, 0.05 - p_value)\n        return penalidade_normalidade\n    def autocorrelation_loss(residuals):\n        acf_vals = acf(residuals, fft=True, nlags=10)\n        penalidade_autocorr = np.sum(np.abs(acf_vals[1:]))\n        return penalidade_autocorr\n    def wavelet_denoising_loss(smoothed_signal, original_signal, alpha=0, beta=1):\n        residuals = original_signal - smoothed_signal\n        loss_linearity = linearity_loss(smoothed_signal)\n        loss_normality = normality_loss(residuals)\n        loss_autocorr = autocorrelation_loss(residuals)\n        loss = alpha * loss_linearity + beta * (loss_normality + loss_autocorr)\n        return loss\n    def wavelet_denoising_with_last_value(signal):\n        \n        if len(signal) % 2 != 0:\n            \n            signal_to_smooth = signal.iloc[:-1]\n            smoothed_signal = wavelet_denoising(signal_to_smooth)  \n            smoothed_signal = np.append(smoothed_signal, signal.iloc[-1])\n        else:\n            smoothed_signal = wavelet_denoising(signal)\n        return smoothed_signal\n    def evaluate_ordering(signal, order):\n        ordered_signal = signal[order]\n        \n        smoothed_signal = wavelet_denoising_with_last_value(ordered_signal)\n        \n        min_val = np.min(ordered_signal)\n        max_val = np.max(ordered_signal)\n        identity_line = np.linspace(min_val, max_val, len(ordered_signal))\n        #loss = mean_squared_error(smoothed_signal[np.argsort(order)], identity_line)\n        #loss  = max(abs(ordered_signal-smoothed_signal))\n        #loss += max(abs(smoothed_signal[np.argsort(order)]-identity_line))\n        normalized_original, mean_orig, std_orig = z_score_normalization(ordered_signal)\n        normalized_smoothed, mean_smooth, std_smooth = z_score_normalization(smoothed_signal)\n        loss = wavelet_denoising_loss(normalized_smoothed, normalized_original)\n        #denormalized_signal = inverse_z_score(normalized_smoothed, mean_smooth, std_smooth)\n        return loss\n    def z_score_normalization(signal):\n        mean = np.mean(signal)\n        std = np.std(signal)\n        normalized_signal = (signal - mean) / std\n        return normalized_signal, mean, std\n    def inverse_z_score(normalized_signal, mean, std):\n        return (normalized_signal * std) + mean\n    def simulated_annealing(signal, max_iter=100000, initial_temp=100, cooling_rate=0.9, min_temp=0.1):\n        n = len(signal)\n        current_order = np.arange(n)\n        current_loss = evaluate_ordering(signal, current_order)\n        best_order = np.copy(current_order)\n        best_loss = current_loss\n        temp = initial_temp\n        for i in range(max_iter):\n            # Troca de dois índices aleatórios no order\n            new_order = np.copy(current_order)\n            idx1, idx2 = np.random.randint(0, n, size=2)\n            new_order[idx1], new_order[idx2] = new_order[idx2], new_order[idx1]\n            new_loss = evaluate_ordering(signal, new_order)\n            # Verifica se a temperatura não caiu abaixo do limite\n            if temp > min_temp:\n                # Probabilidade de aceitação de uma solução pior\n                acceptance_prob = np.exp((current_loss - new_loss) / temp)\n            else:\n                acceptance_prob = 0  # Se a temperatura estiver muito baixa, não aceitamos piores soluções\n            # Atualizar a solução atual se for melhor ou com base na probabilidade\n            if new_loss < current_loss or np.random.rand() < acceptance_prob:\n                current_order = new_order\n                current_loss = new_loss\n                if current_loss < best_loss:\n                    best_order = np.copy(current_order)\n                    best_loss = current_loss\n            # Reduzir a temperatura\n            temp = max(temp * cooling_rate, min_temp)  # Impedir que a temperatura fique muito baixa\n        return best_order, best_loss\n    if dataset == 'train':\n        signal_original = df['a_f_mean_reduction']\n        best_order, best_loss = simulated_annealing(signal_original)\n        print(best_loss)\n        signal_best_ordered = signal_original[best_order]\n        signal_smoothed = wavelet_denoising_with_last_value(signal_best_ordered)\n        df['a_f_mean_reduction_smooth'] = signal_smoothed[np.argsort(best_order)]\n        print(\"ppppp\")\n        print(mean_squared_error(signal_original,df['a_f_mean_reduction']))\n    if f_sliding_features.size > 0:\n        f_sliding_df = pd.DataFrame(f_sliding_features, columns=[f'f_slide_{i}' for i in range(f_sliding_features.shape[1])])\n        df = pd.concat([df, f_sliding_df], axis=1)\n    if a_sliding_features.size > 0:\n        a_sliding_df = pd.DataFrame(a_sliding_features, columns=[f'a_slide_{i}' for i in range(a_sliding_features.shape[1])])\n        df = pd.concat([df, a_sliding_df], axis=1)\n    \n    df = pd.concat([df, adc_info.reset_index().iloc[:, 1:6]], axis=1)\n    return df","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%writefile -a utils.py\n\ndef postprocessing(pred_array, index, sigma_pred):\n    \"\"\"Create a submission dataframe from its components\n    \n    Parameters:\n    pred_array: ndarray of shape (n_samples, 283)\n    index: pandas.Index of length n_samples with name 'planet_id'\n    sigma_pred: float\n    \n    Return value:\n    df: DataFrame of shape (n_samples, 566) with planet_id as index\n    \"\"\"\n    return pd.concat([pd.DataFrame(pred_array.clip(0, None), index=index, columns=wavelengths.columns),\n                      pd.DataFrame(sigma_pred, index=index, columns=[f\"sigma_{i}\" for i in range(1, 284)])],\n                     axis=1)\n\nclass ParticipantVisibleError(Exception):\n    pass\n\ndef competition_score(\n        solution: pd.DataFrame,\n        submission: pd.DataFrame,\n        naive_mean: float,\n        naive_sigma: float,\n        sigma_true: float,\n        row_id_column_name='planet_id',\n    ) -> float:\n    '''\n    This is a Gaussian Log Likelihood based metric. For a submission, which contains the predicted mean (x_hat) and variance (x_hat_std),\n    we calculate the Gaussian Log-likelihood (GLL) value to the provided ground truth (x). We treat each pair of x_hat,\n    x_hat_std as a 1D gaussian, meaning there will be 283 1D gaussian distributions, hence 283 values for each test spectrum,\n    the GLL value for one spectrum is the sum of all of them.\n\n    Inputs:\n        - solution: Ground Truth spectra (from test set)\n            - shape: (nsamples, n_wavelengths)\n        - submission: Predicted spectra and errors (from participants)\n            - shape: (nsamples, n_wavelengths*2)\n        naive_mean: (float) mean from the train set.\n        naive_sigma: (float) standard deviation from the train set.\n        sigma_true: (float) essentially sets the scale of the outputs.\n    '''\n\n    del solution[row_id_column_name]\n    del submission[row_id_column_name]\n\n    if submission.min().min() < 0:\n        raise ParticipantVisibleError('Negative values in the submission')\n    for col in submission.columns:\n        if not pd.api.types.is_numeric_dtype(submission[col]):\n            raise ParticipantVisibleError(f'Submission column {col} must be a number')\n\n    n_wavelengths = len(solution.columns)\n    if len(submission.columns) != n_wavelengths*2:\n        raise ParticipantVisibleError('Wrong number of columns in the submission')\n\n    y_pred = submission.iloc[:, :n_wavelengths].values\n    # Set a non-zero minimum sigma pred to prevent division by zero errors.\n    sigma_pred = np.clip(submission.iloc[:, n_wavelengths:].values, a_min=10**-15, a_max=None)\n    y_true = solution.values\n\n    GLL_pred = np.sum(scipy.stats.norm.logpdf(y_true, loc=y_pred, scale=sigma_pred))\n    GLL_true = np.sum(scipy.stats.norm.logpdf(y_true, loc=y_true, scale=sigma_true * np.ones_like(y_true)))\n    GLL_mean = np.sum(scipy.stats.norm.logpdf(y_true, loc=naive_mean * np.ones_like(y_true), scale=naive_sigma * np.ones_like(y_true)))\n\n    submit_score = (GLL_pred - GLL_mean)/(GLL_true - GLL_mean)\n    return float(np.clip(submit_score, 0.0, 1.0))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"exec(open('utils.py', 'r').read())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load Data","metadata":{}},{"cell_type":"code","source":"%%time\nif os.path.exists(\"/kaggle/input/adc24-intro-training/f_raw_train.pickle\"):\n    f_raw_train = np.load('/kaggle/input/adc24-intro-training/f_raw_train.pickle', allow_pickle=True)\nelse:\n    f_raw_train = read_and_preprocess('train', train_labels.index, 'FGS1')\n    with open('f_raw_train.pickle', 'wb') as f:\n        pickle.dump(f_raw_train, f)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ncontrole = 0\nif os.path.exists(\"/kaggle/input/adc24-intro-training/a_raw_train.pickle\"):\n    a_raw_train = np.load('/kaggle/input/adc24-intro-training/a_raw_train.pickle', allow_pickle=True)\nelse:\n    a_raw_train = read_and_preprocess('train', train_labels.index)\n    with open('a_raw_train.pickle', 'wb') as f:\n        pickle.dump(a_raw_train, f)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Feature Engineering","metadata":{}},{"cell_type":"code","source":"len(train_labels['wl_1'])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ntrain = feature_engineering(f_raw_train, a_raw_train, train_adc_info)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_labels.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data Plot","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(6, 2))\nplt.plot(f_raw_train.mean(axis=0))\nfor time_step in [20500, 23500, 44000, 47000]:\n    plt.axvline(time_step, color='gray')\nplt.xlabel('time step')\nplt.title('FGS1: Overall mean')\nplt.show()\n\nplt.figure(figsize=(6, 2))\nplt.plot(a_raw_train.mean(axis=0))\nfor time_step in [20500, 23500, 44000, 47000]:\n    plt.axvline(time_step * 11250 // 135000, color='gray')\nplt.xlabel('time step')\nplt.title('AIRS-CH0: Overall mean')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"color_array = np.array(plt.rcParams['axes.prop_cycle'].by_key()['color'])\nplt.scatter(train.a_relative_reduction, train_labels.wl_1, s=15, alpha=0.5,\n            c=color_array[train_adc_info.star])\nplt.xlabel('relative signal reduction when planet is in front')\nplt.ylabel('target')\nplt.title('Correlation between relative signal reduction and target')\n# plt.gca().set_aspect('equal')\npoints = [plt.Line2D([0], [0], label=f'star {i}', marker='o', markersize=3,\n         markeredgecolor=color_array[i], markerfacecolor=color_array[i], linestyle='') for i in range(2)]\n\nplt.legend(handles=points)\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Suponha que você já tenha os dados de a_f_mean_reduction e wl_1 no DataFrame train_labels\n# Ordenar os valores de wl_1 e a_f_mean_reduction\n#sorted_indices = np.argsort(train_labels['wl_1'].values)\n# Ordenar a_f_mean_reduction e wl_1 com base nos índices ordenados\n#wl_1_sorted = train_labels['wl_1'].values[sorted_indices]\n#a_f_mean_reduction_sorted = train['a_f_mean_reduction'].values[sorted_indices]\n# Plotar os dados ordenados em um gráfico de linha\nplt.plot(train_labels['wl_1'].values, train['a_f_mean_reduction'].values, marker='o', linestyle='-', color='b')\nplt.xlabel('wl_1')\nplt.ylabel('a_f_mean_reduction')\nplt.title('a_f_mean_reduction vs wl_1 (Ordered)')\nplt.grid(True)\n# Exibir o gráfico\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Suponha que você já tenha os dados de a_f_mean_reduction e wl_1 no DataFrame train_labels\n# Ordenar os valores de wl_1 e a_f_mean_reduction\nsorted_indices = np.argsort(train_labels['wl_1'].values)\n# Ordenar a_f_mean_reduction e wl_1 com base nos índices ordenados\nwl_1_sorted = train_labels['wl_1'].values[sorted_indices]\na_f_mean_reduction_sorted = train['a_f_mean_reduction'].values[sorted_indices]\n# Plotar os dados ordenados em um gráfico de linha\nplt.plot(wl_1_sorted, a_f_mean_reduction_sorted, marker='o', linestyle='-', color='b')\nplt.xlabel('wl_1')\nplt.ylabel('a_f_mean_reduction')\nplt.title('a_f_mean_reduction vs wl_1 (Ordered)')\nplt.grid(True)\n# Exibir o gráfico\nplt.show()\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"color_array = np.array(plt.rcParams['axes.prop_cycle'].by_key()['color'])\nplt.scatter(train.a_f_mean_reduction_smooth, train_labels.wl_1, s=15, alpha=0.5,\n            c=color_array[train_adc_info.star])\nplt.xlabel('relative signal reduction when planet is in front')\nplt.ylabel('target')\nplt.title('Correlation between relative signal reduction and target')\n# plt.gca().set_aspect('equal')\npoints = [plt.Line2D([0], [0], label=f'star {i}', marker='o', markersize=3,\n         markeredgecolor=color_array[i], markerfacecolor=color_array[i], linestyle='') for i in range(2)]\n\nplt.legend(handles=points)\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"color_array = np.array(plt.rcParams['axes.prop_cycle'].by_key()['color'])\nplt.scatter(train.a_f_mean_reduction, train.a_f_mean_reduction_smooth, s=15, alpha=0.5,\n            c=color_array[train_adc_info.star])\nplt.xlabel('relative signal reduction when planet is in front')\nplt.ylabel('target')\nplt.title('Correlation between relative signal reduction and target')\n# plt.gca().set_aspect('equal')\npoints = [plt.Line2D([0], [0], label=f'star {i}', marker='o', markersize=3,\n         markeredgecolor=color_array[i], markerfacecolor=color_array[i], linestyle='') for i in range(2)]\n\nplt.legend(handles=points)\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"max(abs(train.a_f_mean_reduction-train.a_f_mean_reduction_smooth))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model\n## Rigde Model","metadata":{}},{"cell_type":"code","source":"train = train.iloc[:,:-1]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.drop('a_f_mean_reduction', axis=1, inplace=True)\ntrain.drop('a_f_mean_reduction_smooth', axis=1, inplace=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scaler_X = StandardScaler()\nscaler_Y = StandardScaler()\ntrain = scaler_X.fit_transform(train) \ntrain_labels = scaler_Y.fit_transform(train_labels)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.DataFrame(train).head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.DataFrame(train_labels).head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = Ridge(alpha=100)\n\noof_pred = cross_val_predict(model, train, train_labels)\n\nprint(f\"# R2 score: {r2_score(train_labels, oof_pred):.4f}\")\nsigma_pred = mean_squared_error(train_labels, oof_pred, squared=False)\nprint(f\"# Root mean squared error: {sigma_pred:.7f}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.fit(train, train_labels)\nwith open('model.pickle', 'wb') as f:\n    pickle.dump(model, f)\nwith open('sigma_pred.pickle', 'wb') as f:\n    pickle.dump(sigma_pred, f)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Inference","metadata":{}},{"cell_type":"code","source":"# Load the data\ntest_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/test_adc_info.csv',\n                           index_col='planet_id')\nsample_submission = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/sample_submission.csv',\n                                index_col='planet_id')\nf_raw_test = read_and_preprocess('test', sample_submission.index, 'FGS1')\na_raw_test = read_and_preprocess('test', sample_submission.index)\ntest = feature_engineering(f_raw_test, a_raw_test, test_adc_info, dataset = 'test')\ntest = test.iloc[: , :-1]\ntest.drop('a_f_mean_reduction', axis=1, inplace=True)\ntest.drop('a_f_mean_reduction_smooth', axis=1, inplace=True)\ntest = scaler_X.fit_transform(test) \n# Load the model\nwith open('model.pickle', 'rb') as f:\n    model = pickle.load(f)\nwith open('sigma_pred.pickle', 'rb') as f:\n    sigma_pred = pickle.load(f)\n# Predict\ntest_pred = model.predict(test)\ntest_pred = scaler_Y.inverse_transform(test_pred)\n# Package into submission file\nsub_df = sub_df = postprocessing(test_pred,\n                        test_adc_info.index,\n                        sigma_pred=np.tile(np.where(test_adc_info[['star']] <= 1, 0.0001555, 0.00085), (1, 283)))\ndisplay(sub_df)\nsub_df.to_csv('submission.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}