{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","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"},{"sourceId":9331726,"sourceType":"datasetVersion","datasetId":5628160},{"sourceId":9698813,"sourceType":"datasetVersion","datasetId":5930564},{"sourceId":9733481,"sourceType":"datasetVersion","datasetId":5956822},{"sourceId":9751281,"sourceType":"datasetVersion","datasetId":5970158},{"sourceId":196522559,"sourceType":"kernelVersion"}],"dockerImageVersionId":30761,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"This notebook is an update of https://www.kaggle.com/code/sergeifironov/ariel-only-correlation\nfrom Sergei Fironov\n\nUpdates :\n- keep 10:22 pixels from the 32 (the image are well centred)\n- Use the derivative for the determination of the beginning and end of the signal during eclipse (idea from Reza R. Choubeh)\n- 'Simplification' of the code for minimize\n- Degree of polyfit <= 4\n- Predictions of test after training Ridge regression with the modelization results (targets predictions with modelization) and the True targets. ","metadata":{}},{"cell_type":"markdown","source":"# Librairies","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom tqdm import tqdm\nimport joblib\n\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import r2_score, mean_squared_error\nimport itertools\n\nfrom scipy.optimize import minimize\nfrom scipy import optimize\n\nfrom astropy.stats import sigma_clip","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-10-30T13:54:48.976842Z","iopub.execute_input":"2024-10-30T13:54:48.977227Z","iopub.status.idle":"2024-10-30T13:54:51.8121Z","shell.execute_reply.started":"2024-10-30T13:54:48.97717Z","shell.execute_reply":"2024-10-30T13:54:51.810585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nums_parallel = 4","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:54:51.813759Z","iopub.execute_input":"2024-10-30T13:54:51.814304Z","iopub.status.idle":"2024-10-30T13:54:51.820914Z","shell.execute_reply.started":"2024-10-30T13:54:51.81426Z","shell.execute_reply":"2024-10-30T13:54:51.819575Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cut_inf, cut_sup = 39, 321","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:54:51.822591Z","iopub.execute_input":"2024-10-30T13:54:51.823039Z","iopub.status.idle":"2024-10-30T13:54:51.834829Z","shell.execute_reply.started":"2024-10-30T13:54:51.822997Z","shell.execute_reply":"2024-10-30T13:54:51.833647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dataset = 'test'\nadc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/'+f'{dataset}_adc_info.csv',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-30T13:54:51.837803Z","iopub.execute_input":"2024-10-30T13:54:51.838225Z","iopub.status.idle":"2024-10-30T13:54:52.034459Z","shell.execute_reply.started":"2024-10-30T13:54:51.838184Z","shell.execute_reply":"2024-10-30T13:54:52.033258Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Calibration","metadata":{}},{"cell_type":"code","source":"def 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\n\ndef 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\n\ndef preproc(dataset, adc_info, sensor, binning = 15):\n    cut_inf, cut_sup = 0, 356\n    sensor_sizes_dict = {\"AIRS-CH0\":[[11250, 32, 356], [1, 32, cut_sup-cut_inf]], \"FGS1\":[[135000, 32, 32], [1, 32, 32]]}\n    binned_dict = {\"AIRS-CH0\":[11250 // binning // 2, cut_sup-cut_inf], \"FGS1\":[135000 // binning // 2]}\n    linear_corr_dict = {\"AIRS-CH0\":(6, 32, 356), \"FGS1\":(6, 32, 32)}\n    planet_ids = adc_info.index\n    \n    feats = []\n    for i, planet_id in tqdm(list(enumerate(planet_ids))):\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}/' + str(planet_id) + '/' + sensor + '_calibration/dark.parquet', engine='pyarrow').to_numpy()\n        dead_frame = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/' + str(planet_id) + '/' + sensor + '_calibration/dead.parquet', engine='pyarrow').to_numpy()\n        flat_frame = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/' + str(planet_id) + '/' + sensor + '_calibration/flat.parquet', engine='pyarrow').to_numpy()\n        linear_corr = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/' + str(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 = signal / gain + offset\n        \n        hot = sigma_clip(\n            dark_frame, sigma=5, maxiters=5\n        ).mask\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 #@bilzard idea\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            hot = hot[:, cut_inf:cut_sup]\n        else:\n            dt = np.ones(len(signal))*0.1\n            dt[1::2] += 0.1\n            \n        signal = signal.clip(0) #@graySnow idea\n        linear_corr_signal = apply_linear_corr(linear_corr, signal)\n        signal = clean_dark(linear_corr_signal, dark_frame, dt)\n        \n        flat = flat_frame.reshape(sensor_sizes_dict[sensor][1])\n        flat[dead_frame.reshape(sensor_sizes_dict[sensor][1])] = np.nan\n        flat[hot.reshape(sensor_sizes_dict[sensor][1])] = np.nan\n        signal = signal / flat\n        \n        \n        if sensor == \"FGS1\":\n            signal = signal[:,10:22,10:22] # **** updates ****\n            signal = signal.reshape(sensor_sizes_dict[sensor][0][0],144) # # **** updates ****\n\n        if sensor != \"FGS1\":\n            signal = signal[:,10:22,:] # **** updates ****\n\n        mean_signal = np.nanmean(signal, axis=1) \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)","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:54:52.036485Z","iopub.execute_input":"2024-10-30T13:54:52.036996Z","iopub.status.idle":"2024-10-30T13:54:52.075123Z","shell.execute_reply.started":"2024-10-30T13:54:52.036941Z","shell.execute_reply":"2024-10-30T13:54:52.073698Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"f_raw_train = preproc(f'{dataset}', adc_info, \"FGS1\", 30*12)\nf_raw_train.shape","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:54:52.077045Z","iopub.execute_input":"2024-10-30T13:54:52.07756Z","iopub.status.idle":"2024-10-30T13:54:58.659223Z","shell.execute_reply.started":"2024-10-30T13:54:52.077474Z","shell.execute_reply":"2024-10-30T13:54:58.658053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"a_raw_train = preproc(f'{dataset}', adc_info, \"AIRS-CH0\", 30)\na_raw_train.shape","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:54:58.660451Z","iopub.execute_input":"2024-10-30T13:54:58.660838Z","iopub.status.idle":"2024-10-30T13:55:04.910464Z","shell.execute_reply.started":"2024-10-30T13:54:58.660798Z","shell.execute_reply":"2024-10-30T13:55:04.908895Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"window = 37\n\nwv_smooth_a_raw_train = []\nfor i in range(cut_inf, cut_sup):\n    wv_smooth_a_raw_train.append( a_raw_train[:, :, max(i-window, 0):(i+window)].mean(axis=-1) )\n\nwv_smooth_a_raw_train = np.stack(wv_smooth_a_raw_train, axis=-1)\nwv_smooth_a_raw_train.shape","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:04.912673Z","iopub.execute_input":"2024-10-30T13:55:04.913233Z","iopub.status.idle":"2024-10-30T13:55:04.939224Z","shell.execute_reply.started":"2024-10-30T13:55:04.913174Z","shell.execute_reply":"2024-10-30T13:55:04.937728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Fit transit depth","metadata":{"execution":{"iopub.status.busy":"2024-10-03T13:31:39.861441Z","iopub.execute_input":"2024-10-03T13:31:39.861914Z","iopub.status.idle":"2024-10-03T13:31:39.867129Z","shell.execute_reply.started":"2024-10-03T13:31:39.86187Z","shell.execute_reply":"2024-10-03T13:31:39.865794Z"}}},{"cell_type":"code","source":"from concurrent.futures import ProcessPoolExecutor","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:04.940695Z","iopub.execute_input":"2024-10-30T13:55:04.941271Z","iopub.status.idle":"2024-10-30T13:55:04.948224Z","shell.execute_reply.started":"2024-10-30T13:55:04.941192Z","shell.execute_reply":"2024-10-30T13:55:04.946933Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from functools import partial","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:04.950246Z","iopub.execute_input":"2024-10-30T13:55:04.950792Z","iopub.status.idle":"2024-10-30T13:55:04.961819Z","shell.execute_reply.started":"2024-10-30T13:55:04.950736Z","shell.execute_reply":"2024-10-30T13:55:04.960514Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def phase_detector(signal):\n    \n    MIN = np.argmin(signal[30:140])+30\n    signal1 = signal[:MIN ]\n    signal2 = signal[MIN :]\n\n    first_derivative1 = np.gradient(signal1)\n    first_derivative1 /= first_derivative1.max()\n    first_derivative2 = np.gradient(signal2)\n    first_derivative2 /= first_derivative2.max()\n\n    phase1 = np.argmin(first_derivative1)\n    phase2 = np.argmax(first_derivative2) + MIN\n\n    return phase1, phase2\n    \ndef objective(signal, p1, p2, s):\n    \n    best_q = 1e10\n    for i in range(4) :\n        delta = 2\n        x = list(range(signal.shape[0]-delta*4))\n        y = signal[:p1-delta].tolist() + (signal[p1+delta:p2 - delta] * (1 + s)).tolist() + signal[p2+delta:].tolist()\n        \n        z = np.polyfit(x, y, deg=i)\n        p = np.poly1d(z)\n        q = np.abs(p(x) - y).mean()\n    \n    if q < best_q :\n        best_q = q\n    \n    return q","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:04.963415Z","iopub.execute_input":"2024-10-30T13:55:04.963922Z","iopub.status.idle":"2024-10-30T13:55:04.976457Z","shell.execute_reply.started":"2024-10-30T13:55:04.963879Z","shell.execute_reply":"2024-10-30T13:55:04.97503Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def fit_transit_depth(signal, phase_signal):\n    ''' signal: 1d_array\n    '''\n    p1, p2 = phase_detector(phase_signal)\n    f = partial(objective, signal, p1, p2)\n    r = minimize(f, [0.0001], method= 'Nelder-Mead')\n    s = r.x[0]\n    return s","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:04.977923Z","iopub.execute_input":"2024-10-30T13:55:04.978364Z","iopub.status.idle":"2024-10-30T13:55:04.987144Z","shell.execute_reply.started":"2024-10-30T13:55:04.978321Z","shell.execute_reply":"2024-10-30T13:55:04.985875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## A mean depth","metadata":{"execution":{"iopub.status.busy":"2024-10-03T13:40:23.02253Z","iopub.execute_input":"2024-10-03T13:40:23.023691Z","iopub.status.idle":"2024-10-03T13:40:23.028281Z","shell.execute_reply.started":"2024-10-03T13:40:23.023642Z","shell.execute_reply":"2024-10-03T13:40:23.0271Z"}}},{"cell_type":"code","source":"a_raw_train_mean = a_raw_train[:, :, cut_inf:cut_sup].mean(axis=-1)\na_raw_train_mean.shape","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:04.993228Z","iopub.execute_input":"2024-10-30T13:55:04.993713Z","iopub.status.idle":"2024-10-30T13:55:05.002718Z","shell.execute_reply.started":"2024-10-30T13:55:04.99367Z","shell.execute_reply":"2024-10-30T13:55:05.001434Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with ProcessPoolExecutor(nums_parallel) as exe:\n    a_mean_depth = list( tqdm(exe.map(fit_transit_depth, a_raw_train_mean, a_raw_train_mean), total=adc_info.shape[0]) )","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:05.004014Z","iopub.execute_input":"2024-10-30T13:55:05.004431Z","iopub.status.idle":"2024-10-30T13:55:05.154008Z","shell.execute_reply.started":"2024-10-30T13:55:05.004389Z","shell.execute_reply":"2024-10-30T13:55:05.152315Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"a_mean_depth = np.array(a_mean_depth, dtype=np.float32)","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:05.15607Z","iopub.execute_input":"2024-10-30T13:55:05.156487Z","iopub.status.idle":"2024-10-30T13:55:05.163631Z","shell.execute_reply.started":"2024-10-30T13:55:05.156444Z","shell.execute_reply":"2024-10-30T13:55:05.162043Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"a_mean_depth","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:05.165699Z","iopub.execute_input":"2024-10-30T13:55:05.16618Z","iopub.status.idle":"2024-10-30T13:55:05.180771Z","shell.execute_reply.started":"2024-10-30T13:55:05.166136Z","shell.execute_reply":"2024-10-30T13:55:05.17923Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## F depth","metadata":{}},{"cell_type":"code","source":"with ProcessPoolExecutor(nums_parallel) as exe:\n    f_depth = list( tqdm(exe.map(fit_transit_depth, f_raw_train[:, :, 0], a_raw_train_mean), total=adc_info.shape[0]) )","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:05.182353Z","iopub.execute_input":"2024-10-30T13:55:05.182843Z","iopub.status.idle":"2024-10-30T13:55:05.303217Z","shell.execute_reply.started":"2024-10-30T13:55:05.182802Z","shell.execute_reply":"2024-10-30T13:55:05.30159Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"f_depth = np.array(f_depth, dtype=np.float32)","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:05.305445Z","iopub.execute_input":"2024-10-30T13:55:05.305965Z","iopub.status.idle":"2024-10-30T13:55:05.312353Z","shell.execute_reply.started":"2024-10-30T13:55:05.305917Z","shell.execute_reply":"2024-10-30T13:55:05.310975Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"f_depth","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:05.313812Z","iopub.execute_input":"2024-10-30T13:55:05.31416Z","iopub.status.idle":"2024-10-30T13:55:05.328645Z","shell.execute_reply.started":"2024-10-30T13:55:05.314124Z","shell.execute_reply":"2024-10-30T13:55:05.327138Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Every wv depth","metadata":{}},{"cell_type":"code","source":"def fit_every_wv(a_raw, phase_signal):\n    ''' a_raw: T x 282\n    '''\n    res = []\n    for i in range(282):\n        signal = a_raw[:, i]\n        res.append(\n            fit_transit_depth(signal, phase_signal)\n        )\n    return np.array(res, dtype=np.float32).clip(0)","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:05.330219Z","iopub.execute_input":"2024-10-30T13:55:05.33073Z","iopub.status.idle":"2024-10-30T13:55:05.339676Z","shell.execute_reply.started":"2024-10-30T13:55:05.330675Z","shell.execute_reply":"2024-10-30T13:55:05.338389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with ProcessPoolExecutor(nums_parallel) as exe:\n    a_wv_depth = list( tqdm(exe.map(fit_every_wv, wv_smooth_a_raw_train, a_raw_train_mean), total=adc_info.shape[0]) )","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:05.341305Z","iopub.execute_input":"2024-10-30T13:55:05.34176Z","iopub.status.idle":"2024-10-30T13:55:16.175153Z","shell.execute_reply.started":"2024-10-30T13:55:05.34172Z","shell.execute_reply":"2024-10-30T13:55:16.173474Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"a_wv_depth = np.vstack(a_wv_depth)\na_wv_depth.shape","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:16.176933Z","iopub.execute_input":"2024-10-30T13:55:16.177353Z","iopub.status.idle":"2024-10-30T13:55:16.189035Z","shell.execute_reply.started":"2024-10-30T13:55:16.177311Z","shell.execute_reply":"2024-10-30T13:55:16.187579Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Smooth a_wv_depth","metadata":{}},{"cell_type":"code","source":"from scipy.signal import savgol_filter\n\ndef smooth_data(data, window_size):\n    return savgol_filter(data, window_size, 3)  # window size 51, polynomial order 3","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:16.190877Z","iopub.execute_input":"2024-10-30T13:55:16.191315Z","iopub.status.idle":"2024-10-30T13:55:16.280502Z","shell.execute_reply.started":"2024-10-30T13:55:16.191264Z","shell.execute_reply":"2024-10-30T13:55:16.278856Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"window_size = 103\nsmoothed = []\nfor pred in tqdm(a_wv_depth):\n    smooth_pred = smooth_data(pred, window_size=window_size)\n    smoothed.append(smooth_pred)\nsmoothed = np.row_stack(smoothed)\nsmoothed.shape","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:16.282603Z","iopub.execute_input":"2024-10-30T13:55:16.283676Z","iopub.status.idle":"2024-10-30T13:55:16.318503Z","shell.execute_reply.started":"2024-10-30T13:55:16.283627Z","shell.execute_reply":"2024-10-30T13:55:16.317034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"a_wv_depth = smoothed","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:16.320089Z","iopub.execute_input":"2024-10-30T13:55:16.320535Z","iopub.status.idle":"2024-10-30T13:55:16.326115Z","shell.execute_reply.started":"2024-10-30T13:55:16.320467Z","shell.execute_reply":"2024-10-30T13:55:16.324845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 为 FGS1 估算凌日深度","metadata":{}},{"cell_type":"code","source":"window = 64\ntail_a_raw_train = a_raw_train[:, :, -window:].mean(axis=-1)\n    \nwith ProcessPoolExecutor(nums_parallel) as exe:\n    tail_a_depth = list( tqdm(exe.map(fit_transit_depth, tail_a_raw_train, a_raw_train_mean), total=adc_info.shape[0]) )\n\ntail_a_depth = np.array(tail_a_depth, dtype=np.float32)","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:16.327995Z","iopub.execute_input":"2024-10-30T13:55:16.328404Z","iopub.status.idle":"2024-10-30T13:55:16.460088Z","shell.execute_reply.started":"2024-10-30T13:55:16.328362Z","shell.execute_reply":"2024-10-30T13:55:16.458219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"merged_f_depth = tail_a_depth * 0.9 + f_depth * 0.1","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:16.462266Z","iopub.execute_input":"2024-10-30T13:55:16.46273Z","iopub.status.idle":"2024-10-30T13:55:16.472457Z","shell.execute_reply.started":"2024-10-30T13:55:16.462684Z","shell.execute_reply":"2024-10-30T13:55:16.470865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"wv_smooth_pred = np.column_stack([\n    merged_f_depth,\n    a_wv_depth[:, ::-1]\n])\nwv_smooth_pred.shape","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:16.474162Z","iopub.execute_input":"2024-10-30T13:55:16.475386Z","iopub.status.idle":"2024-10-30T13:55:16.488021Z","shell.execute_reply.started":"2024-10-30T13:55:16.475329Z","shell.execute_reply":"2024-10-30T13:55:16.486268Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Estimate sigma by linear regression","metadata":{}},{"cell_type":"code","source":"from sklearn.linear_model import LinearRegression, Ridge","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:16.490062Z","iopub.execute_input":"2024-10-30T13:55:16.490694Z","iopub.status.idle":"2024-10-30T13:55:16.501413Z","shell.execute_reply.started":"2024-10-30T13:55:16.490633Z","shell.execute_reply":"2024-10-30T13:55:16.499981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_wv_smooth_pred = np.load('/kaggle/input/adc2024-train-wv-smooth-pred/train_wv_smooth_pred.npy')\ntrain_wv_smooth_pred.shape","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:16.503112Z","iopub.execute_input":"2024-10-30T13:55:16.503588Z","iopub.status.idle":"2024-10-30T13:55:16.53623Z","shell.execute_reply.started":"2024-10-30T13:55:16.503534Z","shell.execute_reply":"2024-10-30T13:55:16.534773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_label_df = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_labels.csv', index_col='planet_id')\ntrain_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_adc_info.csv', index_col='planet_id')","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:16.537875Z","iopub.execute_input":"2024-10-30T13:55:16.538323Z","iopub.status.idle":"2024-10-30T13:55:16.649714Z","shell.execute_reply.started":"2024-10-30T13:55:16.53828Z","shell.execute_reply":"2024-10-30T13:55:16.648202Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mean_squared_error(train_label_df.values, train_wv_smooth_pred, squared=False)","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:16.651684Z","iopub.execute_input":"2024-10-30T13:55:16.652198Z","iopub.status.idle":"2024-10-30T13:55:16.666765Z","shell.execute_reply.started":"2024-10-30T13:55:16.652141Z","shell.execute_reply":"2024-10-30T13:55:16.664797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sigma_gt = np.abs(train_label_df.values - train_wv_smooth_pred)\nsigma_gt.shape","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:16.668787Z","iopub.execute_input":"2024-10-30T13:55:16.669552Z","iopub.status.idle":"2024-10-30T13:55:16.681902Z","shell.execute_reply.started":"2024-10-30T13:55:16.669482Z","shell.execute_reply":"2024-10-30T13:55:16.68059Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sigma_gt.min(), sigma_gt.mean(), sigma_gt.max(), sigma_gt.std()","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:16.683427Z","iopub.execute_input":"2024-10-30T13:55:16.684Z","iopub.status.idle":"2024-10-30T13:55:16.696967Z","shell.execute_reply.started":"2024-10-30T13:55:16.683957Z","shell.execute_reply":"2024-10-30T13:55:16.695608Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(12, 4))\n_ = plt.hist(sigma_gt.flatten(), bins=500)","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:16.698405Z","iopub.execute_input":"2024-10-30T13:55:16.698826Z","iopub.status.idle":"2024-10-30T13:55:18.047859Z","shell.execute_reply.started":"2024-10-30T13:55:16.698787Z","shell.execute_reply":"2024-10-30T13:55:18.046537Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"#### 行星粒度的 sigma 与各波段 depth 的差异正相关","metadata":{"execution":{"iopub.status.busy":"2024-10-29T07:30:38.314135Z","iopub.execute_input":"2024-10-29T07:30:38.314619Z","iopub.status.idle":"2024-10-29T07:30:38.319941Z","shell.execute_reply.started":"2024-10-29T07:30:38.314578Z","shell.execute_reply":"2024-10-29T07:30:38.318737Z"}}},{"cell_type":"code","source":"x = sigma_gt.mean(axis=-1)\ny = train_wv_smooth_pred.std(axis=-1)\nr = np.corrcoef(x, y)[1, 0]\nplt.scatter(x, y)\nplt.title(f'r={r}')","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:18.049382Z","iopub.execute_input":"2024-10-30T13:55:18.049916Z","iopub.status.idle":"2024-10-30T13:55:18.445213Z","shell.execute_reply.started":"2024-10-30T13:55:18.049866Z","shell.execute_reply":"2024-10-30T13:55:18.443905Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### LR for sigma prediction","metadata":{}},{"cell_type":"code","source":"star0_flag = train_adc_info.star.values == 0\nstar1_flag = train_adc_info.star.values == 1\nstar0_flag.sum(), star1_flag.sum()","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:18.446815Z","iopub.execute_input":"2024-10-30T13:55:18.447214Z","iopub.status.idle":"2024-10-30T13:55:18.456756Z","shell.execute_reply.started":"2024-10-30T13:55:18.447173Z","shell.execute_reply":"2024-10-30T13:55:18.455192Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X = train_wv_smooth_pred.std(axis=-1).reshape(-1, 1)\nX.shape","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:18.45821Z","iopub.execute_input":"2024-10-30T13:55:18.458616Z","iopub.status.idle":"2024-10-30T13:55:18.472869Z","shell.execute_reply.started":"2024-10-30T13:55:18.458576Z","shell.execute_reply":"2024-10-30T13:55:18.47153Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"star0_model = LinearRegression()\nstar0_model.fit(X[star0_flag], sigma_gt[star0_flag])","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:18.474596Z","iopub.execute_input":"2024-10-30T13:55:18.475081Z","iopub.status.idle":"2024-10-30T13:55:18.511726Z","shell.execute_reply.started":"2024-10-30T13:55:18.475039Z","shell.execute_reply":"2024-10-30T13:55:18.510472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"star1_model = LinearRegression()\nstar1_model.fit(X[star1_flag], sigma_gt[star1_flag])","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:18.525121Z","iopub.execute_input":"2024-10-30T13:55:18.52832Z","iopub.status.idle":"2024-10-30T13:55:18.548549Z","shell.execute_reply.started":"2024-10-30T13:55:18.528238Z","shell.execute_reply":"2024-10-30T13:55:18.547467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mix_model = LinearRegression()\nmix_model.fit(X, sigma_gt)","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:18.553693Z","iopub.execute_input":"2024-10-30T13:55:18.556455Z","iopub.status.idle":"2024-10-30T13:55:18.577482Z","shell.execute_reply.started":"2024-10-30T13:55:18.55639Z","shell.execute_reply":"2024-10-30T13:55:18.576405Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### wv average sigma","metadata":{}},{"cell_type":"code","source":"loaded_data = np.load('/kaggle/input/adc2024-train-star-wv-rmse/star_wv_rmse.npz')\navg_wv_rmse = loaded_data['avg_wv_rmse']\nstar0_wv_rmse = loaded_data['star0_wv_rmse']\nstar1_wv_rmse = loaded_data['star1_wv_rmse']","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:18.579268Z","iopub.execute_input":"2024-10-30T13:55:18.580302Z","iopub.status.idle":"2024-10-30T13:55:18.599484Z","shell.execute_reply.started":"2024-10-30T13:55:18.580245Z","shell.execute_reply":"2024-10-30T13:55:18.598366Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"avg_wv_rmse.min(), star0_wv_rmse.min(), star1_wv_rmse.min()","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:18.601218Z","iopub.execute_input":"2024-10-30T13:55:18.602005Z","iopub.status.idle":"2024-10-30T13:55:18.612564Z","shell.execute_reply.started":"2024-10-30T13:55:18.601946Z","shell.execute_reply":"2024-10-30T13:55:18.611354Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"avg_wv_rmse.mean(), star0_wv_rmse.mean(), star1_wv_rmse.mean()","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:18.614813Z","iopub.execute_input":"2024-10-30T13:55:18.615899Z","iopub.status.idle":"2024-10-30T13:55:18.629291Z","shell.execute_reply.started":"2024-10-30T13:55:18.615837Z","shell.execute_reply":"2024-10-30T13:55:18.627902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.plot(star0_wv_rmse, color='red')\nplt.plot(star1_wv_rmse, color='green')\nplt.plot(avg_wv_rmse, color='blue')","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:18.631516Z","iopub.execute_input":"2024-10-30T13:55:18.632477Z","iopub.status.idle":"2024-10-30T13:55:18.944997Z","shell.execute_reply.started":"2024-10-30T13:55:18.63242Z","shell.execute_reply":"2024-10-30T13:55:18.943563Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Sigma prediction","metadata":{}},{"cell_type":"code","source":"pred.shape","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:18.946738Z","iopub.execute_input":"2024-10-30T13:55:18.947276Z","iopub.status.idle":"2024-10-30T13:55:18.956915Z","shell.execute_reply.started":"2024-10-30T13:55:18.947216Z","shell.execute_reply":"2024-10-30T13:55:18.955463Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"wv_smooth_pred = wv_smooth_pred.clip(0)","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:18.958805Z","iopub.execute_input":"2024-10-30T13:55:18.960247Z","iopub.status.idle":"2024-10-30T13:55:18.969158Z","shell.execute_reply.started":"2024-10-30T13:55:18.96018Z","shell.execute_reply":"2024-10-30T13:55:18.96773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sigma = []\nfor star, x in zip(adc_info.star, wv_smooth_pred.std(axis=-1)):\n    if star == 0:\n        out = star0_model.predict(x.reshape(1, -1))\n        out = out.clip(3e-5) * 1.5\n    elif star == 1:\n        out = star1_model.predict(x.reshape(1, -1))\n        out = out.clip(3e-5) * 1.5\n    else:\n        out = mix_model.predict(x.reshape(1, -1))\n        out = out.clip(3e-5) * 1.5\n    \n    sigma.append(out)\n\nsigma = np.row_stack(sigma)\nsigma.shape","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:18.971313Z","iopub.execute_input":"2024-10-30T13:55:18.971814Z","iopub.status.idle":"2024-10-30T13:55:18.988721Z","shell.execute_reply.started":"2024-10-30T13:55:18.971756Z","shell.execute_reply":"2024-10-30T13:55:18.987552Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.plot(sigma[0], color='red')\nplt.plot(star1_wv_rmse * 1.5, color='blue')","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:18.990539Z","iopub.execute_input":"2024-10-30T13:55:18.991421Z","iopub.status.idle":"2024-10-30T13:55:19.280989Z","shell.execute_reply.started":"2024-10-30T13:55:18.991374Z","shell.execute_reply":"2024-10-30T13:55:19.279563Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Submission","metadata":{}},{"cell_type":"code","source":"ss = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/sample_submission.csv')\n\n# sigma = np.ones_like(wv_smooth_pred) * 0.000145\n# sigma = np.tile(\n#     np.where(adc_info[['star']] <= 1, 0.00013, 0.00018),\n#     (1, 283)\n# )\n\npred = wv_smooth_pred.clip(0)  # N x 283\nsubmission = pd.DataFrame(np.concatenate([pred,sigma], axis=1), columns=ss.columns[1:])\nsubmission.index = adc_info.index\nsubmission.to_csv('submission.csv')\nsubmission.head()","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:19.282797Z","iopub.execute_input":"2024-10-30T13:55:19.283351Z","iopub.status.idle":"2024-10-30T13:55:19.347703Z","shell.execute_reply.started":"2024-10-30T13:55:19.28329Z","shell.execute_reply":"2024-10-30T13:55:19.346348Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.plot(submission.iloc[0].values[:283])","metadata":{"execution":{"iopub.status.busy":"2024-10-30T13:55:19.349004Z","iopub.execute_input":"2024-10-30T13:55:19.349403Z","iopub.status.idle":"2024-10-30T13:55:19.622235Z","shell.execute_reply.started":"2024-10-30T13:55:19.349361Z","shell.execute_reply":"2024-10-30T13:55:19.620654Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}