{"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":"### ℹ️ **Info**\n* **forked original great work kernels**\n    * https://www.kaggle.com/code/sergeifironov/ariel-only-correlation\n\n* **2024/09/08 My Changed**\n    * scipy minimize() param & other params update\n* **2024/09/22 My Changed**\n    * improve LB .517 -> .522","metadata":{}},{"cell_type":"markdown","source":"---\n---","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport seaborn as sns\nimport scipy.stats\nfrom tqdm import tqdm\n\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import r2_score, mean_squared_error\nimport itertools\nfrom scipy.optimize import minimize\nfrom functools import partial\nimport random, os\nfrom astropy.stats import sigma_clip","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-10-17T18:32:06.725163Z","iopub.execute_input":"2024-10-17T18:32:06.725576Z","iopub.status.idle":"2024-10-17T18:32:11.238383Z","shell.execute_reply.started":"2024-10-17T18:32:06.725526Z","shell.execute_reply":"2024-10-17T18:32:11.237115Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/test_adc_info.csv',\n                           index_col='planet_id')\naxis_info = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/axis_info.parquet')","metadata":{"execution":{"iopub.status.busy":"2024-10-17T18:32:28.790364Z","iopub.execute_input":"2024-10-17T18:32:28.790954Z","iopub.status.idle":"2024-10-17T18:32:29.062025Z","shell.execute_reply.started":"2024-10-17T18:32:28.790916Z","shell.execute_reply":"2024-10-17T18:32:29.060981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/test_adc_info.csv')\n\nplanet_ids = adc_info['planet_id']\n\ncut_inf, cut_sup = 39, 321\nbinning = 100\nlinear_corr_dict = {\"AIRS-CH0\":(6, 32, 356), \"FGS1\":(6, 32, 32)}\nbinned_dict = {\"AIRS-CH0\":[11250 // binning // 2, 282], \"FGS1\":[135000 // binning // 2]}\nsensor_sizes_dict = {\"AIRS-CH0\":[[11250, 32, 356], [1, 32, cut_sup-cut_inf]], \"FGS1\":[[135000, 32, 32], [1, 32, 32]]}\n\npath_folder = '/kaggle/input/ariel-data-challenge-2024/'\npath_out = '/kaggle/tmp/data_light_raw/' \noutput_dir = '/kaggle/tmp/data_light_raw/'\n\nif not os.path.exists(path_out):\n    os.makedirs(path_out)\n    print(f\"Directory {path_out} created.\")\nelse:\n    print(f\"Directory {path_out} already exists.\")\n    \ntrain_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_adc_info.csv')\ntest_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/test_adc_info.csv')\n\ntrain_labels = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_labels.csv')\nwavelengths = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/wavelengths.csv')\naxis_info = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/axis_info.parquet')\n\ndef ADC_convert(signal, gain, offset):\n    signal = signal.astype(np.float64)\n    signal /= gain\n    signal += offset\n    return signal\n\nimport numpy as np\nfrom astropy.stats import sigma_clip\nfrom scipy.ndimage import convolve\n\ndef correct_bad_pixels(signal, dead, dark, sigma=5, maxiters=5):\n    \"\"\"\n    Masks and corrects dead and hot pixels in the input signal by interpolation.\n\n    Parameters:\n    signal : np.ndarray\n        The input signal frames, assumed of shape (num_frames, height, width).\n    dead : np.ndarray\n        A boolean mask indicating dead pixels (True where pixels are dead), of shape (height, width).\n    dark : np.ndarray\n        The dark frame used for identifying hot pixels, of shape (height, width).\n    sigma : float, optional\n        The sigma level to use for identifying hot pixels (default is 5).\n    maxiters : int, optional\n        The maximum number of iterations for sigma clipping (default is 5).\n\n    Returns:\n    np.ndarray\n        The signal with hot and dead pixels corrected using interpolation.\n    \"\"\"\n    # Create mask for hot pixels\n    hot_mask = sigma_clip(dark, sigma=sigma, maxiters=maxiters).mask  # Boolean mask for hot pixels\n    combined_mask = np.logical_or(hot_mask, dead)  # Combine hot and dead masks\n\n    # Prepare output array\n    corrected_signal = np.copy(signal)\n\n    # Create a kernel for convolution (average 3x3)\n    kernel = np.ones((3, 3)) / 9.0  # Normalized averaging kernel\n\n    # Iterate over frames\n    for i in range(signal.shape[0]):\n        frame = signal[i]\n        mask = combined_mask  # Use the combined mask for the current frame\n        \n        # Create a temporary frame for correction\n        temp_frame = np.where(mask, np.nan, frame)  # Mask the frame\n\n        # Use convolution to fill NaN values with the average of surrounding pixels\n        # Use a temporary array to perform the convolution\n        filled_frame = convolve(np.nan_to_num(temp_frame, nan=0), kernel, mode='constant', cval=0)\n\n        # Update the corrected signal with the filled values where the mask is True\n        corrected_signal[i] = np.where(mask, filled_frame, frame)\n\n    return corrected_signal\n\ndef subtract_read_noise(signal, read_noise):\n    # Broadcast read_noise to match signal shape if needed\n    read_noise = np.expand_dims(read_noise, axis=0)  # Assuming read_noise has shape (rows, columns)\n\n    # Subtract read noise from the signal frames\n    signal = signal - read_noise\n\n    return signal\n\ndef apply_linear_corr(linear_corr, clean_signal):\n    # Step 1: Flip the coefficients along the first axis\n    linear_corr = np.array(linear_corr)\n    linear_corr = linear_corr.reshape(6, 32, clean_signal.shape[2])\n    linear_corr = np.flip(linear_corr, axis=0)\n\n    # Step 2: Loop through all pixels (ignoring the time dimension) and apply correction\n    for x, y in itertools.product(range(clean_signal.shape[1]), range(clean_signal.shape[2])):\n        # Create a polynomial object with the correction coefficients for pixel (x, y)\n        poli = np.poly1d(linear_corr[:, x, y])\n        \n        # Apply the polynomial to each pixel's signal values (along the time axis)\n        clean_signal[:, x, y] = poli(clean_signal[:, x, y])\n    \n    # Step 3: Return the corrected signal\n    return clean_signal\n\n\ndef correct_flat_field(flat, dead, signal):\n    # Transpose the flat field and dead pixel maps to match the signal dimensions\n    flat = np.array(flat)\n    #flat = flat.transpose()\n    #dead = dead.transpose()\n\n    # Mask out the dead pixels from the flat field map\n    #flat = np.ma.masked_where(dead, flat)\n    flat[flat<0.5] = 1\n    # Replicate the flat field map for each time slice in the signal\n    flat = np.tile(flat, (signal.shape[0], 1, 1))\n\n    # Correct the signal by dividing by the flat field map\n    signal = signal / flat\n    return signal\n\ndef clean_dark(signal, dead, dark, dt):\n    \"\"\"\n    Cleans the dark current from the signal.\n    \n    Parameters:\n    - signal: 3D numpy array (time, rows, columns) representing the signal over time.\n    - dead: 2D numpy array (rows, columns) representing dead pixels.\n    - dark: 2D numpy array (rows, columns) representing the dark current map.\n    - dt: 1D numpy array (time,) representing the time scaling factor for each frame.\n    \n    Returns:\n    - Cleaned signal after dark current subtraction.\n    \"\"\"\n\n    # Mask the dark current map where dead pixels are present\n    dark = np.where(dead, 0, dark)\n    \n    # Subtract the scaled dark current from the signal using broadcasting\n    signal -= (dark[np.newaxis, :, :] * dt[:, np.newaxis, np.newaxis]).astype(np.uint16)\n    \n    return signal\n\ndef get_cds(signal):\n    cds = signal[:,1::2,:] - signal[:,::2,:]\n    return cds\n\ndef correct_flat_field(flat, dead, signal):\n    # Transpose the flat field and dead pixel maps to match the signal dimensions\n    #flat = flat.transpose(1, 0)  # Assuming flat is initially (rows, columns)\n    #dead = dead.transpose(1, 0)  # Assuming dead is initially (rows, columns)\n\n    # Mask out the dead pixels from the flat field map\n    flat = np.ma.masked_where(dead, flat)\n\n    # Ensure that flat is a 3D array for broadcasting\n    flat = np.expand_dims(flat, axis=0)  # Shape becomes (1, rows, columns)\n\n    # Perform flat-field correction using broadcasting\n    # Now signal is (time, rows, columns), flat is (1, rows, columns)\n    signal = signal / flat  # Broadcasting happens here\n\n    return signal\n\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\n\ndef new_preproc(dataset, adc_info, sensor, binning):\n    feats = []\n    cut_inf, cut_sup = 39, 321\n    linear_corr_dict = {\"AIRS-CH0\":(6, 32, 356), \"FGS1\":(6, 32, 32)}\n    binned_dict = {\"AIRS-CH0\":[11250 // binning // 2, 282], \"FGS1\":[135000 // binning // 2]}\n    sensor_sizes_dict = {\"AIRS-CH0\":[[11250, 32, 356], [1, 32, cut_sup-cut_inf]], \"FGS1\":[[135000, 32, 32], [1, 32, 32]]}\n\n    for i, planet_id in tqdm(list(enumerate(planet_ids[:5]))):\n        # Load the signal and calibration frames\n        signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/{sensor}_signal.parquet').to_numpy()\n        dark_frame = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/{sensor}_calibration/dark.parquet').to_numpy()\n        dead_frame = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/{sensor}_calibration/dead.parquet').to_numpy()\n        flat_frame = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/{sensor}_calibration/flat.parquet').to_numpy()\n        linear_corr = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/{sensor}_calibration/linear_corr.parquet').values.astype(np.float64).reshape(linear_corr_dict[sensor])\n\n        signal = signal.reshape(sensor_sizes_dict[sensor][0]) \n        gain = adc_info[f'{sensor}_adc_gain'].values[i]\n        offset = adc_info[f'{sensor}_adc_offset'].values[i]\n        signal = ADC_convert(signal, gain, offset)\n\n        \n        if sensor != \"FGS1\":\n            signal = signal[:, :, cut_inf:cut_sup]\n            dt = np.ones(len(signal)) * 0.1\n            dt[1::2] += 4.5\n            linear_corr = linear_corr[:, :, cut_inf:cut_sup]\n            dark_frame = dark_frame[:, cut_inf:cut_sup]\n            dead_frame = dead_frame[:, cut_inf:cut_sup]\n            flat_frame = flat_frame[:, cut_inf:cut_sup]\n        else:\n            dt = np.ones(len(signal)) * 0.1\n            dt[1::2] += 0.1\n        \n        \n        signal = correct_bad_pixels(signal, dead_frame, dark_frame, sigma=5, maxiters=5)\n\n        signal = signal.clip(0)\n\n        linear_corr_signal = apply_linear_corr(linear_corr, signal)\n        \n        signal = clean_dark(linear_corr_signal, dead_frame, dark_frame, dt)\n        \n        signal = correct_flat_field(flat_frame, dead_frame, signal)\n        \n        if sensor == \"FGS1\":\n            signal = signal.reshape((sensor_sizes_dict[sensor][0][0], sensor_sizes_dict[sensor][0][1] * sensor_sizes_dict[sensor][0][2]))\n            \n            \n        mean_signal = np.nanmean(signal, axis=1) # mean over the 32*32(FGS1) or 32(CH0) pixels\n        cds_signal = (mean_signal[1::2] - mean_signal[0::2])\n        \n        binned = np.zeros((binned_dict[sensor]))\n        for j in range(cds_signal.shape[0] // binning):\n            binned[j] = cds_signal[j*binning:j*binning+binning].mean(axis=0)\n                   \n        if sensor == \"FGS1\":\n            binned = binned.reshape((binned.shape[0],1))\n        \n        feats.append(binned)\n        \n    return np.stack(feats)\n    \npre_train = np.concatenate([new_preproc('test', test_adc_info, \"FGS1\", 30*12), new_preproc('test', test_adc_info, \"AIRS-CH0\", 30)], axis=2)    ","metadata":{"execution":{"iopub.status.busy":"2024-10-17T18:32:30.742576Z","iopub.execute_input":"2024-10-17T18:32:30.743006Z","iopub.status.idle":"2024-10-17T18:33:12.923236Z","shell.execute_reply.started":"2024-10-17T18:32:30.742968Z","shell.execute_reply":"2024-10-17T18:33:12.921898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pre_train.shape","metadata":{"execution":{"iopub.status.busy":"2024-10-17T18:33:14.974291Z","iopub.execute_input":"2024-10-17T18:33:14.974665Z","iopub.status.idle":"2024-10-17T18:33:14.983551Z","shell.execute_reply.started":"2024-10-17T18:33:14.974637Z","shell.execute_reply":"2024-10-17T18:33:14.98225Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### fit polynoms for each sample","metadata":{}},{"cell_type":"code","source":"def phase_detector(signal):\n    phase1, phase2 = None, None\n    best_drop = 0\n    for i in range(2//2,187//2):        \n        t1 = signal[i:i+20//2].max() - signal[i:i+20//2].min()\n        if t1 > best_drop:\n            phase1 = i+(20+5)//2\n            best_drop = t1\n            \n        # 25 se 75\n        # i to i+10 me max-min\n    \n    best_drop = 0\n    for i in range(187//2,372//2):\n        t1 = signal[i:i+20//2].max() - signal[i:i+20//2].min()\n        if t1 > best_drop:\n            phase2 = i-5//2\n            best_drop = t1\n    \n    return phase1, phase2\n\ndef try_s(signal, p1, p2, deg, s):\n    out = list(range(p1-30)) + list(range(p2+30,signal.shape[0]))\n    x, y = out, signal[out].tolist()\n    x = x + list(range(p1,p2))\n\n    y = y + (signal[p1:p2] * (1 + s[0])).tolist()\n    z = np.polyfit(x, y, deg)\n    p = np.poly1d(z)\n    q = np.abs(p(x) - y).mean()\n\n    if s < 1e-4:\n        return q + 1e3\n\n    return q\n    \ndef calibrate_signal(signal):\n    p1,p2 = phase_detector(signal)\n\n    best_deg, best_score = 1, 1e12\n    for deg in range(1, 4):\n        f = partial(try_s, signal, p1, p2, deg)\n        r = minimize(f, [0.001], method = 'Nelder-Mead')\n        s = r.x[0]\n\n        out = list(range(p1-30)) + list(range(p2+30,signal.shape[0]))\n        x, y = out, signal[out].tolist()\n        x = x + list(range(p1,p2))\n        y = y + (signal[p1:p2] * (1 + s)).tolist()\n    \n        z = np.polyfit(x, y, deg)\n        p = np.poly1d(z)\n        q = np.abs(p(x) - y).mean()\n        \n        if q < best_score:\n            best_score = q\n            best_deg = deg\n        \n        print(deg, q)\n            \n    z = np.polyfit(x, y, best_deg)\n    p = np.poly1d(z)\n\n    return s, x, y, p(x)\n\ndef calibrate_train(signal):\n    p1,p2 = phase_detector(signal)\n    \n    best_deg, best_score = 1, 1e12\n    for deg in range(1, 4):\n        f = partial(try_s, signal, p1, p2, deg)\n        r = minimize(f, [0.0001], method = 'Nelder-Mead')\n        s = r.x[0]\n\n        out = list(range(p1-30)) + list(range(p2+30,signal.shape[0]))\n        x, y = out, signal[out].tolist()\n        x = x + list(range(p1,p2))\n        y = y + (signal[p1:p2] * (1 + s)).tolist()\n    \n        z = np.polyfit(x, y, deg)\n        p = np.poly1d(z)\n        q = np.abs(p(x) - y).mean()\n        \n        if q < best_score:\n            best_score = q\n            best_deg = deg\n            \n    z = np.polyfit(x, y, best_deg)\n    p = np.poly1d(z)\n    \n    return s, p(np.arange(signal.shape[0])), p1, p2\n\n\ntrain = pre_train.copy()\nall_s = []\nfor i in range(len(test_adc_info)):\n    signal = train[i,:,1:].mean(axis=1)\n    s, p, p1, p2 = calibrate_train(pre_train[i,:,1:].mean(axis=1))\n    all_s.append(s)\n        \n#copy answer 283 times because we predict mean value\ntrain_s = np.repeat(np.array(all_s), 283).reshape((len(all_s), 283))        \ntrain_sigma = np.ones_like(train_s) * 0.000176","metadata":{"execution":{"iopub.status.busy":"2024-10-17T18:33:16.803003Z","iopub.execute_input":"2024-10-17T18:33:16.80349Z","iopub.status.idle":"2024-10-17T18:33:16.885719Z","shell.execute_reply.started":"2024-10-17T18:33:16.803457Z","shell.execute_reply":"2024-10-17T18:33:16.884578Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Probably we can accurately estimate sigma from train","metadata":{}},{"cell_type":"code","source":"n = 0\ns, x, y, y_new = calibrate_signal(pre_train[n,:,1:].mean(axis=1))\nplt.scatter(x,y)\nplt.scatter(x,y_new)","metadata":{"execution":{"iopub.status.busy":"2024-10-17T18:33:23.061577Z","iopub.execute_input":"2024-10-17T18:33:23.062024Z","iopub.status.idle":"2024-10-17T18:33:23.437508Z","shell.execute_reply.started":"2024-10-17T18:33:23.061989Z","shell.execute_reply":"2024-10-17T18:33:23.436346Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I call the orange line \"starline\". This is probably what we would see if the planet weren't in the way.","metadata":{}},{"cell_type":"markdown","source":"### Making submission","metadata":{}},{"cell_type":"code","source":"ss = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/sample_submission.csv')\n\npreds = train_s.clip(0)\nsigmas = train_sigma\nsubmission = pd.DataFrame(np.concatenate([preds,sigmas], axis=1), columns=ss.columns[1:])\nsubmission.index = test_adc_info.index\nsubmission.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2024-10-17T18:33:25.749704Z","iopub.execute_input":"2024-10-17T18:33:25.750619Z","iopub.status.idle":"2024-10-17T18:33:25.793894Z","shell.execute_reply.started":"2024-10-17T18:33:25.750573Z","shell.execute_reply":"2024-10-17T18:33:25.792544Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission","metadata":{"execution":{"iopub.status.busy":"2024-10-17T18:33:26.260963Z","iopub.execute_input":"2024-10-17T18:33:26.261359Z","iopub.status.idle":"2024-10-17T18:33:26.293392Z","shell.execute_reply.started":"2024-10-17T18:33:26.261332Z","shell.execute_reply":"2024-10-17T18:33:26.292269Z"},"trusted":true},"execution_count":null,"outputs":[]}]}