{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.10.13"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"},{"sourceId":194955113,"sourceType":"kernelVersion"}],"dockerImageVersionId":30746,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false},"papermill":{"default_parameters":{},"duration":1087.424982,"end_time":"2024-08-03T13:00:18.164826","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2024-08-03T12:42:10.739844","version":"2.5.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd\nimport polars as pl\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport seaborn as sns\nimport scipy.stats\nfrom tqdm import tqdm\nimport pickle\n\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import r2_score, mean_squared_error\n\nimport os\nimport itertools\nfrom scipy.interpolate import splrep, BSpline\nimport multiprocessing\n\nimport numba as nb","metadata":{"_kg_hide-input":true,"papermill":{"duration":3.023083,"end_time":"2024-08-03T12:42:17.107326","exception":false,"start_time":"2024-08-03T12:42:14.084243","status":"completed"},"tags":[],"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:45:11.571442Z","iopub.execute_input":"2024-11-01T16:45:11.571846Z","iopub.status.idle":"2024-11-01T16:45:17.038585Z","shell.execute_reply.started":"2024-11-01T16:45:11.571813Z","shell.execute_reply":"2024-11-01T16:45:17.037024Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_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('/kaggle/input/ariel-data-challenge-2024/wavelengths.csv')\naxis_info = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/axis_info.parquet')\nstar = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_adc_info.csv').star.to_numpy()","metadata":{"papermill":{"duration":0.188106,"end_time":"2024-08-03T12:42:17.37713","exception":false,"start_time":"2024-08-03T12:42:17.189024","status":"completed"},"tags":[],"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:45:17.179204Z","iopub.execute_input":"2024-11-01T16:45:17.180474Z","iopub.status.idle":"2024-11-01T16:45:17.564318Z","shell.execute_reply.started":"2024-11-01T16:45:17.180427Z","shell.execute_reply":"2024-11-01T16:45:17.562702Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Data processing funcs","metadata":{}},{"cell_type":"code","source":"path_folder = '/kaggle/input/ariel-data-challenge-2024/'\ncut_inf, cut_sup = 39, 321","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:45:17.565862Z","iopub.execute_input":"2024-11-01T16:45:17.566264Z","iopub.status.idle":"2024-11-01T16:45:17.5725Z","shell.execute_reply.started":"2024-11-01T16:45:17.566231Z","shell.execute_reply":"2024-11-01T16:45:17.571056Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile ADC_convert.py\ndef ADC_convert(signal, gain, offset):\n    signal = signal.astype(np.float64)\n    signal /= gain\n    signal += offset\n    return signal","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:45:17.574437Z","iopub.execute_input":"2024-11-01T16:45:17.575016Z","iopub.status.idle":"2024-11-01T16:45:17.588139Z","shell.execute_reply.started":"2024-11-01T16:45:17.574976Z","shell.execute_reply":"2024-11-01T16:45:17.586688Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile apply_linear_corr.py\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":{"iopub.status.busy":"2024-11-01T16:45:17.589639Z","iopub.execute_input":"2024-11-01T16:45:17.590045Z","iopub.status.idle":"2024-11-01T16:45:17.600107Z","shell.execute_reply.started":"2024-11-01T16:45:17.590002Z","shell.execute_reply":"2024-11-01T16:45:17.59868Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile clean_dark.py\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","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:45:17.60182Z","iopub.execute_input":"2024-11-01T16:45:17.602365Z","iopub.status.idle":"2024-11-01T16:45:17.611547Z","shell.execute_reply.started":"2024-11-01T16:45:17.602277Z","shell.execute_reply":"2024-11-01T16:45:17.610257Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile correct_flat_field.py\ndef correct_flat_field(flat, signal):\n    flat = np.tile(flat, (signal.shape[0], 1, 1))\n    signal = signal / flat\n    return signal","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:45:17.61296Z","iopub.execute_input":"2024-11-01T16:45:17.61344Z","iopub.status.idle":"2024-11-01T16:45:17.623951Z","shell.execute_reply.started":"2024-11-01T16:45:17.613395Z","shell.execute_reply":"2024-11-01T16:45:17.622719Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"exec(open('ADC_convert.py', 'r').read())\nexec(open('apply_linear_corr.py', 'r').read())\nexec(open('clean_dark.py', 'r').read())\nexec(open('correct_flat_field.py', 'r').read())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:45:17.629588Z","iopub.execute_input":"2024-11-01T16:45:17.629981Z","iopub.status.idle":"2024-11-01T16:45:17.638102Z","shell.execute_reply.started":"2024-11-01T16:45:17.62995Z","shell.execute_reply":"2024-11-01T16:45:17.636742Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile a_read_and_preprocess.py\ndef a_read_and_preprocess(dataset, adc_info, planet_ids):\n    a_raw_train = np.full((len(planet_ids), 5625, 282), np.nan, dtype=np.float64)\n    for i, planet_id in tqdm(list(enumerate(planet_ids))):\n        signal = pl.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/AIRS-CH0_signal.parquet')\n        signal = np.reshape(signal.to_numpy(), (-1, 32, 356)).astype(np.float64)\n        \n        gain = adc_info['AIRS-CH0_adc_gain'].loc[planet_id]\n        offset = adc_info['AIRS-CH0_adc_offset'].loc[planet_id]\n        signal = ADC_convert(signal, gain, offset)\n        signal = np.clip(signal, a_min = 0, a_max = None)\n        dt_airs = axis_info['AIRS-CH0-integration_time'].dropna().values\n        signal = signal[:, :, cut_inf:cut_sup]\n        \n        flat = pd.read_parquet(os.path.join(path_folder,f'{dataset}/{planet_id}/AIRS-CH0_calibration/flat.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n        dark = pd.read_parquet(os.path.join(path_folder,f'{dataset}/{planet_id}/AIRS-CH0_calibration/dark.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n        dead = pd.read_parquet(os.path.join(path_folder,f'{dataset}/{planet_id}/AIRS-CH0_calibration/dead.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n        linear_corr = pd.read_parquet(os.path.join(path_folder,f'{dataset}/{planet_id}/AIRS-CH0_calibration/linear_corr.parquet')).values.astype(np.float64).reshape((6, 32, 356))[:, :, cut_inf:cut_sup]\n\n        signal = apply_linear_corr(linear_corr,signal)\n        \n        signal = clean_dark(signal, dark, dt_airs)\n        \n        \n        signal = correct_flat_field(flat, signal)\n        \n        \n        dead_norm_repair = signal.copy()\n        dead_norm_repair = np.concatenate([dead_norm_repair[:,:,1:2], dead_norm_repair, dead_norm_repair[:,:,-2:-1]], axis = -1)\n        dead_norm_repair = (dead_norm_repair[:,:,:-2]+dead_norm_repair[:,:,2:])/2\n        signal = np.where(dead, dead_norm_repair , signal)\n        \n\n        signal = np.mean((signal[1::2] - signal[0::2]), axis = 1)\n        a_raw_train[i] = signal\n    return a_raw_train","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:45:17.64059Z","iopub.execute_input":"2024-11-01T16:45:17.641562Z","iopub.status.idle":"2024-11-01T16:45:17.651598Z","shell.execute_reply.started":"2024-11-01T16:45:17.641523Z","shell.execute_reply":"2024-11-01T16:45:17.650051Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\nexec(open('a_read_and_preprocess.py', 'r').read())\na_raw_train = a_read_and_preprocess('train', train_adc_info, train_labels.index)\nwith open('a_raw_train.pickle', 'wb') as f:\n    pickle.dump(a_raw_train, f)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:45:17.65341Z","iopub.execute_input":"2024-11-01T16:45:17.653885Z","iopub.status.idle":"2024-11-01T16:46:09.642944Z","shell.execute_reply.started":"2024-11-01T16:45:17.653841Z","shell.execute_reply":"2024-11-01T16:46:09.64157Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Feature engineering","metadata":{}},{"cell_type":"markdown","source":"I based my Numba implementation on [np](https://github.com/numpy/numpy/blob/v2.1.0/numpy/lib/_polynomial_impl.py#L702-L779) and on [this Numba implementation](https://gist.github.com/kadereub/9eae9cff356bb62cdbd672931e8e5ec4) (basically modified to follow the NumPy conventions).","metadata":{}},{"cell_type":"code","source":"@nb.njit\ndef vander(x, deg):\n    v = np.zeros((len(x),deg+1))\n    for n in range(0, deg + 1):\n        v[:, n] = x**n\n    return v\n    \n@nb.jit(nopython=True)\ndef polyfit_numba(x, y, deg):\n    lhs = vander(x, deg)\n    p = np.linalg.lstsq(lhs, y)[0]\n    return p[::-1]\n\n@nb.jit(nopython=True)\ndef polyval_numba(p, x):\n    y = x*0\n    for pv in p:\n        y = y * x + pv\n    return y","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:46:09.644739Z","iopub.execute_input":"2024-11-01T16:46:09.645212Z","iopub.status.idle":"2024-11-01T16:46:09.650763Z","shell.execute_reply.started":"2024-11-01T16:46:09.645168Z","shell.execute_reply":"2024-11-01T16:46:09.649405Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile get_fit_numba.py\n\ndef get_fit_numba(x, y, deg = 2):\n    fit_pol = polyfit_numba(x, y, deg)\n    y_fitted = polyval_numba(fit_pol, x)\n    rms = np.sum((y-y_fitted)**2)\n    return y_fitted, rms, fit_pol","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:46:09.652449Z","iopub.execute_input":"2024-11-01T16:46:09.652836Z","iopub.status.idle":"2024-11-01T16:46:09.665611Z","shell.execute_reply.started":"2024-11-01T16:46:09.652805Z","shell.execute_reply":"2024-11-01T16:46:09.664247Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile get_rms_list_numba.py\n\ndef get_rms_list_numba(points_1, points_2, signal):\n    '''\n    For a signal, a list of indices on the signal points_1 and another list of indices points_2,\n    for every combination of p1 in points_1 and p2 in points_2 such that p1<p2, it calculate \n    the total RMSE between the signal and a combination of 2nd order polynomial fitted on the signal\n    up to p1, a 2nd polynomial fitted on the signal from p2 to the end, and a line that connects the \n    end of the first polynomial with the start of the second one.\n    '''\n    rms_list = []\n    p1_list = []\n    p2_list = []\n    rms_1_list = []\n    eval_1_list = []\n    for p1 in points_1:\n        x1 = np.asarray(range(0, p1+1))*1.0\n        y1 = signal[:p1+1]\n        y_fitted_1, rms_1, fit_pol_1 = get_fit_numba(x1,y1)\n        rms_1_list.append(rms_1)\n        eval_1_list.append(polyval_numba(fit_pol_1, x1[-1]))\n    rms_3_list = []\n    eval_3_list = []\n    for p2 in points_2:\n        x3 = np.asarray(range(p2, len(signal)))*1.0\n        y3 = signal[p2:len(signal)]\n        y_fitted_3, rms_3, fit_pol_3 = get_fit_numba(x3,y3)\n        rms_3_list.append(rms_3)\n        eval_3_list.append(polyval_numba(fit_pol_3, x3[0]))\n    for idx1, p1 in enumerate(points_1):\n        for idx2, p2 in enumerate(points_2):\n            if p1>=5 and p2<len(signal)-5 and p1<p2-5:\n                x1 = np.asarray(range(0, p1+1))*1.0\n                x2 = np.asarray(range(p1, p2+1))*1.0\n                y2 = signal[p1:p2+1]\n                x3 = np.asarray(range(p2, len(signal)))*1.0\n                \n                fit_2 = polyfit_numba(np.asarray([x1[-1], x3[0]]),\n                                 np.asarray([eval_1_list[idx1],\n                                             eval_3_list[idx2]]), deg = 1)\n                y_fitted_2 = polyval_numba(fit_2, x2)\n                rms_2 = np.sum((y2-y_fitted_2)**2)\n                \n                rms_1 = rms_1_list[idx1]\n                rms_3 = rms_3_list[idx2]\n                rms_list.append(rms_1+rms_2+rms_3)\n                p1_list.append(p1)\n                p2_list.append(p2)\n    return rms_list, np.asarray(p1_list), np.asarray(p2_list)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:46:09.667635Z","iopub.execute_input":"2024-11-01T16:46:09.668216Z","iopub.status.idle":"2024-11-01T16:46:09.680206Z","shell.execute_reply.started":"2024-11-01T16:46:09.668142Z","shell.execute_reply":"2024-11-01T16:46:09.678633Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"exec(open('get_fit_numba.py', 'r').read())\nexec(open('get_rms_list_numba.py', 'r').read())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:46:09.681872Z","iopub.execute_input":"2024-11-01T16:46:09.68244Z","iopub.status.idle":"2024-11-01T16:46:09.698738Z","shell.execute_reply.started":"2024-11-01T16:46:09.682404Z","shell.execute_reply":"2024-11-01T16:46:09.696202Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def get_best_points_by_pol_fitting(signal, fit_stride, initial_points):\n    '''\n    Apply get_rms_list_numba on the first half and on the second half of the original signal in order \n    to find the ingress and egress.\n    '''\n    x1,x2,x3,x4,x5,x6,x7,x8 = initial_points\n    points_1 = np.asarray(range(x1, x2, fit_stride))\n    points_2 = np.asarray(range(x3, x4, fit_stride))\n    points_3 = np.asarray(range(x5, x6, fit_stride))\n    points_4 = np.asarray(range(x7, x8, fit_stride))\n\n    window_mean_mid = int(len(signal)/2)\n    rms_list_first, p1_list_first, p2_list_first = get_rms_list_numba(\n        points_1, points_2, signal[:window_mean_mid])\n    rms_list_second, p1_list_second, p2_list_second = get_rms_list_numba(\n        points_3-window_mean_mid, points_4-window_mean_mid, signal[window_mean_mid:])\n    p1_list_second, p2_list_second = p1_list_second+window_mean_mid, p2_list_second+window_mean_mid\n\n    idx_min = np.argmin(rms_list_first)\n    p1 = p1_list_first[idx_min]\n    p2 = p2_list_first[idx_min]\n    idx_min = np.argmin(rms_list_second)\n    p3 = p1_list_second[idx_min]\n    p4 = p2_list_second[idx_min]\n    return p1,p2,p3,p4","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:46:09.700398Z","iopub.execute_input":"2024-11-01T16:46:09.700831Z","iopub.status.idle":"2024-11-01T16:46:09.713102Z","shell.execute_reply.started":"2024-11-01T16:46:09.7008Z","shell.execute_reply":"2024-11-01T16:46:09.711695Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile pol_fit_get_points.py\n\ndef pol_fit_get_points(item):\n    '''\n    Apply get_best_points_by_pol_fitting three times consecutively on the signal, first with a stride of 100,\n    then with a stride of 20 on the best section found in the first iteration, and finally with a stride of 1\n    on the best section found in the second iteration. This is a kind of modified binary search.  \n    '''\n    signal, idx = item\n    cumsum_window = 3\n    cumsum = np.cumsum(np.mean(signal, axis = 1))\n    window_mean = (cumsum[cumsum_window:]-cumsum[:-cumsum_window])/cumsum_window\n\n    window_mean_mid = int(len(window_mean)/2)\n    initial_points = [0,window_mean_mid, 0,window_mean_mid,\n                      window_mean_mid,len(window_mean), window_mean_mid,len(window_mean)]\n    fit_strides = 100\n    p1,p2,p3,p4 = get_best_points_by_pol_fitting(\n        window_mean, fit_stride = fit_strides, initial_points = initial_points)\n\n\n    initial_points = [p1-fit_strides, p1+fit_strides, p2-fit_strides, p2+fit_strides,\n                      p3-fit_strides, p3+fit_strides, p4-fit_strides, p4+fit_strides]\n    fit_strides = 20\n    p1,p2,p3,p4 = get_best_points_by_pol_fitting(\n        window_mean, fit_stride = fit_strides, initial_points = initial_points)\n\n\n    initial_points = [p1-fit_strides, p1+fit_strides, p2-fit_strides, p2+fit_strides,\n                      p3-fit_strides, p3+fit_strides, p4-fit_strides, p4+fit_strides]\n    fit_strides = 1\n    p1,p2,p3,p4 = get_best_points_by_pol_fitting(\n        window_mean, fit_stride = fit_strides, initial_points = initial_points)\n\n    points = np.asarray([p1,p2,p3,p4])+cumsum_window//2\n    if idx%10 == 0:\n        print(idx)\n    return points","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:46:09.714659Z","iopub.execute_input":"2024-11-01T16:46:09.715131Z","iopub.status.idle":"2024-11-01T16:46:09.730978Z","shell.execute_reply.started":"2024-11-01T16:46:09.715091Z","shell.execute_reply":"2024-11-01T16:46:09.729696Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"exec(open('pol_fit_get_points.py', 'r').read())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:46:09.732461Z","iopub.execute_input":"2024-11-01T16:46:09.732837Z","iopub.status.idle":"2024-11-01T16:46:09.74471Z","shell.execute_reply.started":"2024-11-01T16:46:09.732799Z","shell.execute_reply":"2024-11-01T16:46:09.743392Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile pol_fit_get_relreds_freq.py\n\ndef pol_fit_get_relreds_freq(item):\n    '''\n    relred = RELative REDuction, i.e., the absorption magnitude we seek. In this method, I calculate\n    the relred using the points p1,p2,p3,p4 I found in pol_fit_get_points, with p1 and p2 the start/end of the \n    ingress and p3/p4 the start/end of the egress. I fit a 2nd polynom up to p1, between p2 to p3 and from p4 \n    onward, and calculate the difference between the y-values of the end/start of the first/second polynomials and\n    the start/end of the third/second polynomial. This method is noisy and was not used in the final solution, but it is \n    presented here anyway for comparison.\n    '''\n    signal, points = item\n    cumsum_window = 3\n    points = points-cumsum_window//2\n    cumsum = np.cumsum(signal, axis = 0)\n    window_mean = (cumsum[cumsum_window:]-cumsum[:-cumsum_window])/cumsum_window\n    EPS_1 = 0\n    EPS_2 = 0\n\n    p1, p2, p3, p4 = points\n\n    freq_window = 0\n    data_list = []\n    data_list_pols = []\n    for freq_idx in range(282):\n        min_freq = max(0, freq_idx-freq_window)\n        max_freq = min(freq_idx+freq_window+1, 282)\n        signaL_mean_freq = np.mean(window_mean[:, min_freq:max_freq], axis = -1)\n\n        y1 = signaL_mean_freq[EPS_1:p1+1-EPS_2]\n        y2 = signaL_mean_freq[p1:p2+1]\n        y3_1 = signaL_mean_freq[p2+EPS_2:p3+1-EPS_1]\n        y3_2 = signaL_mean_freq[p2+EPS_1:p3+1-EPS_2]\n        y4 = signaL_mean_freq[p3:p4+1-EPS_1]\n        y5 = signaL_mean_freq[p4+EPS_2:len(window_mean)-EPS_1]\n        \n        x1 = np.asarray(range(EPS_1, p1+1-EPS_2))*1.0\n        x3_1 = np.asarray(range(p2+EPS_2, p3+1-EPS_1))*1.0\n        x3_2 = np.asarray(range(p2+EPS_1, p3+1-EPS_2))*1.0\n        x5 = np.asarray(range(p4+EPS_2, len(window_mean)-EPS_1))*1.0\n        y_fitted_1, rms_1, fit_pol_1 = get_fit_numba(x1,y1)\n        y_fitted_3_1, rms_3_1, fit_pol_3_1 = get_fit_numba(x3_1,y3_1)\n        y_fitted_3_2, rms_3_2, fit_pol_3_2 = get_fit_numba(x3_2,y3_2)\n        y_fitted_5, rms_5, fit_pol_5 = get_fit_numba(x5,y5)\n\n        delta_1_y_1 = polyval_numba(fit_pol_1, (p1+p2)/2)\n        delta_1_y_2 = polyval_numba(fit_pol_3_1, (p1+p2)/2)\n\n        delta_2_y_1 = polyval_numba(fit_pol_3_2, (p3+p4)/2)\n        delta_2_y_2 = polyval_numba(fit_pol_5, (p3+p4)/2)\n\n        data = [delta_1_y_1, delta_1_y_2, delta_2_y_1, delta_2_y_2]\n        data_list.append(data)\n        data_list_pols.append([fit_pol_1, fit_pol_3_1, fit_pol_3_2, fit_pol_5])\n    return np.asarray(data_list), data_list_pols","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:46:31.249558Z","iopub.status.idle":"2024-11-01T16:46:31.250047Z","shell.execute_reply.started":"2024-11-01T16:46:31.249818Z","shell.execute_reply":"2024-11-01T16:46:31.249838Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"exec(open('pol_fit_get_relreds_freq.py', 'r').read())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:46:31.252483Z","iopub.status.idle":"2024-11-01T16:46:31.253049Z","shell.execute_reply.started":"2024-11-01T16:46:31.252811Z","shell.execute_reply":"2024-11-01T16:46:31.252834Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile pol_fit_advanced_method.py\ndef pol_fit_advanced_method(item):\n    '''\n    The 'advanced' method for calculating the relreds, inspired by Sergei's work:\n    https://www.kaggle.com/code/sergeifironov/ariel-only-correlation\n    In this approach, we move the middle part of the signal up or down by different absorption constants, then fit \n    a polynom on the side and middle together, and the constant that corresponds to the best fit (lowest RMSE) is the\n    target. Sergei performed the fitting on both sides and the middle together; I preferred to divide it into two parts,\n    fitting separately for the right side with the middle and the left side with the middle. In addition, different\n    from Sergei's use of scipy.optimize.minimize, I wanted direct control and understanding of everything that happened\n    'under the hood,' so I use a custom, modified implementation of binary search.\n    '''\n    signal, points, delta_orig = item\n\n    cumsum_window = 3\n    points = points-cumsum_window//2\n    cumsum = np.cumsum(signal, axis = 0)\n    window_mean = (cumsum[cumsum_window:]-cumsum[:-cumsum_window])/cumsum_window\n    EPS_1 = 0\n    EPS_2 = 0\n    EPS_3 = 10000\n    EPS_4 = 10000\n    DEG = 2 \n    \n    p1, p2, p3, p4 = points\n    freq_window = 0\n\n    baseline_list = []\n    delta_list = []\n    err_list = []\n    fit_pols = []\n    for freq_idx in range(282):\n        min_freq = max(0, freq_idx-freq_window)\n        max_freq = min(freq_idx+freq_window+1, 282)\n        signaL_mean_freq = np.mean(window_mean[:, min_freq:max_freq], axis = -1)\n\n        y1 = signaL_mean_freq[max(p1-EPS_3, 0)+EPS_1:p1+1-EPS_2]\n        y2 = signaL_mean_freq[p1:p2+1]\n        y3 = signaL_mean_freq[p2+EPS_2:min(p2+EPS_4, p3)+1-EPS_1]\n        y4 = signaL_mean_freq[p3:p4+1]\n        y5 = signaL_mean_freq[p4:]\n\n        x1 = np.asarray(range(max(p1-EPS_3, 0)+EPS_1, p1+1-EPS_2))*1.0\n        x3 = np.asarray(range(p2+EPS_2, min(p2+EPS_4, p3)+1-EPS_1))*1.0\n        x5 = np.asarray(range(p4, len(window_mean)))*1.0\n\n        x_concat = np.concatenate([x1,x3])\n        y_concat = np.concatenate([y1,y3])\n        D1 = len(x1)\n        D3 = len(x3)\n        D5 = len(x5)\n\n        mid = 0\n        candidate = y_concat.copy()\n        candidate[D1:D1+D3] = candidate[D1:D1+D3]+mid\n        y_fitted_candidate, rms_candidate, fit_pol_candidate = get_fit_numba(x_concat,candidate, deg = DEG)\n        e_mid = rms_candidate\n        delta = 40\n        for i in range(23):\n            top = mid+delta\n            candidate = y_concat.copy()\n            candidate[D1:D1+D3] = candidate[D1:D1+D3]+top\n            y_fitted_candidate, rms_candidate, fit_pol_candidate = get_fit_numba(x_concat,candidate, deg = DEG)\n            e_top = rms_candidate\n            bott= mid-delta\n            candidate = y_concat.copy()\n            candidate[D1:D1+D3] = candidate[D1:D1+D3]+bott\n            y_fitted_candidate, rms_candidate, fit_pol_candidate = get_fit_numba(x_concat,candidate, deg = DEG)\n            e_bott = rms_candidate\n            if e_top<e_mid:\n                mid = top\n                e_mid = e_top\n            elif e_bott<e_mid:\n                mid = bott\n                e_mid = e_bott\n            delta = delta/2\n        candidate = y_concat.copy()\n        candidate[D1:D1+D3] = candidate[D1:D1+D3]+mid\n        y_fitted_candidate, rms_candidate, fit_pol_candidate = get_fit_numba(x_concat,candidate, deg = DEG)\n        baseline_list.append((y_fitted_candidate[D1]+y_fitted_candidate[D1+1])/2)\n        delta_list.append(mid)\n        err_list.append(rms_candidate)\n        fit_pols.append(fit_pol_candidate)\n\n    baseline_list_2 = []\n    delta_list_2 = []\n    err_list_2 = []\n    fit_pols_2 = []\n\n    for freq_idx in range(282):\n        min_freq = max(0, freq_idx-freq_window)\n        max_freq = min(freq_idx+freq_window+1, 282)\n        signaL_mean_freq = np.mean(window_mean[:, min_freq:max_freq], axis = -1)\n\n        y1 = signaL_mean_freq[:p1+1]\n        y2 = signaL_mean_freq[p1:p2+1]\n        y3 = signaL_mean_freq[max(p3-EPS_4, p2)+EPS_1:p3+1-EPS_2]\n        y4 = signaL_mean_freq[p3:p4+1]\n        y5 = signaL_mean_freq[p4+EPS_2:min(p4+EPS_3+1, len(window_mean))-EPS_1]\n\n        x1 = np.asarray(range(0, p1+1))*1.0\n        x3 = np.asarray(range(max(p3-EPS_4, p2)+EPS_1, p3+1-EPS_2))*1.0\n        x5 = np.asarray(range(p4+EPS_2, min(p4+EPS_3+1, len(window_mean))-EPS_1))*1.0\n\n        x_concat = np.concatenate([x3,x5])\n        y_concat = np.concatenate([y3,y5])\n        D1 = len(x1)\n        D3 = len(x3)\n        D5 = len(x5)\n\n        mid = 0\n        candidate = y_concat.copy()\n        candidate[:D3] = candidate[:D3]+mid\n        y_fitted_candidate, rms_candidate, fit_pol_candidate = get_fit_numba(x_concat,candidate, deg = DEG)\n        e_mid = rms_candidate\n        delta = 40\n        for i in range(23):\n            top = mid+delta\n            candidate = y_concat.copy()\n            candidate[:D3] = candidate[:D3]+top\n            y_fitted_candidate, rms_candidate, fit_pol_candidate = get_fit_numba(x_concat,candidate, deg = DEG)\n            e_top = rms_candidate\n            bott= mid-delta\n            candidate = y_concat.copy()\n            candidate[:D3] = candidate[:D3]+bott\n            y_fitted_candidate, rms_candidate, fit_pol_candidate = get_fit_numba(x_concat,candidate, deg = DEG)\n            e_bott = rms_candidate\n            if e_top<e_mid:\n                mid = top\n                e_mid = e_top\n            elif e_bott<e_mid:\n                mid = bott\n                e_mid = e_bott\n            delta = delta/2\n        candidate = y_concat.copy()\n        candidate[:D3] = candidate[:D3]+mid\n        y_fitted_candidate, rms_candidate, fit_pol_candidate = get_fit_numba(x_concat,candidate, deg = DEG)\n        baseline_list_2.append((y_fitted_candidate[D3]+y_fitted_candidate[D3+1])/2)\n        delta_list_2.append(mid)\n        err_list_2.append(rms_candidate)\n        fit_pols_2.append(fit_pol_candidate)\n    return baseline_list, delta_list, err_list, baseline_list_2, delta_list_2, err_list_2, fit_pols, fit_pols_2","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:46:31.262405Z","iopub.status.idle":"2024-11-01T16:46:31.26286Z","shell.execute_reply.started":"2024-11-01T16:46:31.262642Z","shell.execute_reply":"2024-11-01T16:46:31.26266Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"exec(open('pol_fit_advanced_method.py', 'r').read())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:46:31.265177Z","iopub.status.idle":"2024-11-01T16:46:31.265662Z","shell.execute_reply.started":"2024-11-01T16:46:31.265453Z","shell.execute_reply":"2024-11-01T16:46:31.265472Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile get_range_fit.py\ndef get_range_fit(global_start,global_end, mean_signal):\n    '''\n    Fit 2nd polynomial on different sub-sections. Initially, I used it for various experiments, but in the final submission, I only needed it for\n    the coefficient of the fit on the entire sections.  \n    '''\n    DEG = 2\n    rms_list_1 = []\n    fit_pol_list_1 = []\n    start = global_start\n    for end in range(global_start+100, global_end, 1):\n        x = np.asarray(range(start, end))*1.0\n        y = mean_signal[start:end]\n        y_fitted, rms, fit_pol = get_fit_numba(x,y, deg = DEG)\n        rms_list_1.append(rms)\n        fit_pol_list_1.append(fit_pol)\n    fit_pol_list_1 = np.asarray(fit_pol_list_1)\n\n    rms_list_2 = []\n    fit_pol_list_2 = []\n    end = global_end\n    for start in range(global_start, global_end-100, 1):\n        x = np.asarray(range(start, end))*1.0\n        y = mean_signal[start:end]\n        y_fitted, rms, fit_pol = get_fit_numba(x,y, deg = DEG)\n        rms_list_2.append(rms)\n        fit_pol_list_2.append(fit_pol)\n    fit_pol_list_2 = np.asarray(fit_pol_list_2)\n    \n    return rms_list_1, fit_pol_list_1, rms_list_2, fit_pol_list_2","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"exec(open('get_range_fit.py', 'r').read())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile fit_on_sections.py\ndef fit_on_sections(item):\n    '''\n    Apply the function get_range_fit on the three sections, start middle and end (before between and after ingress/egress).\n    '''\n    signal = item[0]\n    points = item[1]\n\n    mean_signal = np.mean(signal, axis = 1)\n    rms_list_1, fit_pol_list_1, rms_list_2, fit_pol_list_2 = get_range_fit(0,points[0], mean_signal)\n    rms_list_3, fit_pol_list_3, rms_list_4, fit_pol_list_4 = get_range_fit(points[1],points[2], mean_signal)\n    rms_list_5, fit_pol_list_5, rms_list_6, fit_pol_list_6 = get_range_fit(points[3], len(mean_signal), mean_signal)\n    \n    data = [rms_list_1, fit_pol_list_1, rms_list_2, fit_pol_list_2,\n            rms_list_3, fit_pol_list_3, rms_list_4, fit_pol_list_4,\n            rms_list_5, fit_pol_list_5, rms_list_6, fit_pol_list_6]\n    \n    return data","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"exec(open('fit_on_sections.py', 'r').read())","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile feature_engineering.py\n\ndef feature_engineering(a_raw, save = True, points_list_load = None, to_print = True):\n    rel_red_1_list = []\n    rel_red_2_list = []\n    points_list = []\n    if points_list_load is not None:\n        points_list = points_list_load\n    else:\n        items = [[a_raw[i], i] for i in range(len(a_raw))]\n        results = []\n        with multiprocessing.Pool() as pool:\n            for result in pool.map(pol_fit_get_points, items):\n                results.append(result)\n        points_list = np.asarray(results)\n        \n    items = [[a_raw[i], points_list[i]] for i in range(len(a_raw))]\n    results = []\n    with multiprocessing.Pool() as pool:\n        for result in pool.map(pol_fit_get_relreds_freq, items):\n            results.append(result)\n    data_list = [x[0] for x in results]\n    data_list_pols = [x[1] for x in results]\n\n    data_list = np.asarray(data_list)\n    delta_1_y_1 = data_list[:, :, 0]\n    delta_1_y_2 = data_list[:, :, 1]\n    delta_2_y_1 = data_list[:, :, 2]\n    delta_2_y_2 = data_list[:, :, 3]\n    delta1 = delta_1_y_1-delta_1_y_2\n    delta2 = delta_2_y_2-delta_2_y_1\n    rel_red_1_list = (delta1)/delta_1_y_1\n    rel_red_2_list = (delta2)/delta_2_y_2\n        \n    a_relative_reduction = (rel_red_1_list+rel_red_2_list)/2\n    \n    delta_mean = (delta1+delta2)/2\n    items = [[a_raw[i], points_list[i], delta_mean[i]] for i in range(len(a_raw))]\n    advanced_fitting_results = []\n    with multiprocessing.Pool() as pool:\n        for result in pool.map(pol_fit_advanced_method, items):\n            advanced_fitting_results.append(result)\n\n    items = [[a_raw[i], points_list[i]] for i in range(len(a_raw))]\n    results_fit_on_sections = []\n    with multiprocessing.Pool() as pool:\n        for result in pool.map(fit_on_sections, items):\n            results_fit_on_sections.append(result)\n\n    return points_list, data_list_pols, data_list, advanced_fitting_results, results_fit_on_sections","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:46:31.285635Z","iopub.status.idle":"2024-11-01T16:46:31.286209Z","shell.execute_reply.started":"2024-11-01T16:46:31.285917Z","shell.execute_reply":"2024-11-01T16:46:31.285941Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%time\nexec(open('feature_engineering.py', 'r').read())\n\npoints_list_new, data_list_pols, data_list, advanced_fitting_results, results_fit_on_sections  = feature_engineering(\n    a_raw_train, save=False, to_print=True)\npickle.dump(points_list_new, open('points_list.p', 'bw'))\npickle.dump(data_list_pols, open('data_list_pols.p', 'bw'))\npickle.dump(data_list, open('data_list.p', 'bw'))\npickle.dump(advanced_fitting_results, open('advanced_fitting_results.p', 'bw'))\npickle.dump(results_fit_on_sections, open('results_fit_on_sections.p', 'bw'))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2024-11-01T16:46:31.287961Z","iopub.status.idle":"2024-11-01T16:46:31.288424Z","shell.execute_reply.started":"2024-11-01T16:46:31.288184Z","shell.execute_reply":"2024-11-01T16:46:31.288202Z"}},"outputs":[],"execution_count":null}]}