{"metadata":{"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"},{"sourceId":191155259,"sourceType":"kernelVersion"}],"dockerImageVersionId":30746,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true},"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"},"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":"markdown","source":"# NeurIPS - Ariel Data Challenge 2024\n\n## Kaggle competitions\n\nCompleted by Timur Garipov for the Master's Research Seminar at Novosibirsk State University","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":0.009657,"end_time":"2024-08-03T12:42:14.074227","exception":false,"start_time":"2024-08-03T12:42:14.06457","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"# Import","metadata":{}},{"cell_type":"code","source":"import os\nimport pickle\nimport numpy as np\nimport scipy.stats\nimport pandas as pd\nimport polars as pl\nimport seaborn as sns\nfrom tqdm import tqdm\nimport matplotlib.pyplot as plt\n\n\nfrom sklearn.linear_model import Ridge\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.metrics import r2_score, mean_squared_error\n\nRANDOM_SEED = 154\nPATH = \"/kaggle/input/ariel-data-challenge-2024\"","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":[],"execution":{"iopub.status.busy":"2024-09-28T15:57:20.030912Z","iopub.execute_input":"2024-09-28T15:57:20.031749Z","iopub.status.idle":"2024-09-28T15:57:22.900281Z","shell.execute_reply.started":"2024-09-28T15:57:20.03171Z","shell.execute_reply":"2024-09-28T15:57:22.899253Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Utils","metadata":{}},{"cell_type":"code","source":"%%writefile utils.py\n\n\ndef load_signal_data(planet_id, dataset, instrument, img_size):\n    signal = pl.read_parquet(f'{PATH}/{dataset}/{planet_id}/{instrument}_signal.parquet')\n    mean_signal = signal.cast(pl.Int32).sum_horizontal().cast(pl.Float32).to_numpy() / img_size\n    net_signal = mean_signal[1::2] - mean_signal[0::2]\n    return net_signal\n\ndef f_read_and_preprocess(dataset, planet_ids, instrument=\"AIRS-CH0\"):\n    img_size = 1024 if instrument == \"FGS1\" else 32*356\n    column_num = 67500 if instrument == 'FGS1' else 5625\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\n\ndef feature_engineering(f_raw, a_raw, adc_info, window_size=50, step_size=15):\n    f_obscured = f_raw[:, 23500:44000].mean(axis=1)\n    f_unobscured = (f_raw[:, :20500].mean(axis=1) + f_raw[:, 47000:].mean(axis=1)) / 2\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    a_obscured = a_raw[:, 1958:3666].mean(axis=1)\n    a_unobscured = (a_raw[:, :1708].mean(axis=1) + a_raw[:, 3916:].mean(axis=1)) / 2\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    # 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    print(f'f_sliding_features.shape: {f_sliding_features.shape}')\n    print(f'a_sliding_features.shape: {a_sliding_features.shape}')\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    })\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    \n    return df\n\n\n\ndef postprocessing(pred_array, index, sigma_pred):\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    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    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":{"execution":{"iopub.status.busy":"2024-09-28T15:57:45.181765Z","iopub.execute_input":"2024-09-28T15:57:45.182314Z","iopub.status.idle":"2024-09-28T15:57:45.194271Z","shell.execute_reply.started":"2024-09-28T15:57:45.182278Z","shell.execute_reply":"2024-09-28T15:57:45.193188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"exec(open('utils.py', 'r').read())","metadata":{"execution":{"iopub.status.busy":"2024-09-28T15:57:45.863425Z","iopub.execute_input":"2024-09-28T15:57:45.863785Z","iopub.status.idle":"2024-09-28T15:57:45.871412Z","shell.execute_reply.started":"2024-09-28T15:57:45.863754Z","shell.execute_reply":"2024-09-28T15:57:45.870419Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data","metadata":{"papermill":{"duration":0.008428,"end_time":"2024-08-03T12:42:17.180297","exception":false,"start_time":"2024-08-03T12:42:17.171869","status":"completed"},"tags":[]}},{"cell_type":"code","source":"train_adc_info = pd.read_csv(f'{PATH}/train_adc_info.csv', index_col='planet_id')\ntrain_labels = pd.read_csv(f'{PATH}/train_labels.csv', index_col='planet_id')\nwavelengths = pd.read_csv(f'{PATH}/wavelengths.csv')","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":[],"execution":{"iopub.status.busy":"2024-09-28T15:57:46.64009Z","iopub.execute_input":"2024-09-28T15:57:46.640943Z","iopub.status.idle":"2024-09-28T15:57:46.776955Z","shell.execute_reply.started":"2024-09-28T15:57:46.640907Z","shell.execute_reply":"2024-09-28T15:57:46.775892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nf_raw_train = f_read_and_preprocess('train', train_labels.index, 'FGS1')\na_raw_train = f_read_and_preprocess('train', train_labels.index)\n\nwith open('f_raw_train.pickle', 'wb') as f:\n    pickle.dump(f_raw_train, f)\n    \nwith open('a_raw_train.pickle', 'wb') as f:\n    pickle.dump(a_raw_train, f)","metadata":{"execution":{"iopub.status.busy":"2024-09-28T15:57:46.98093Z","iopub.execute_input":"2024-09-28T15:57:46.981291Z","iopub.status.idle":"2024-09-28T16:25:33.38541Z","shell.execute_reply.started":"2024-09-28T15:57:46.981259Z","shell.execute_reply":"2024-09-28T16:25:33.384379Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Feature engineering","metadata":{}},{"cell_type":"code","source":"%%time\ntrain = feature_engineering(f_raw_train, a_raw_train, train_adc_info)\ntrain = train.drop(columns=['star'])","metadata":{"execution":{"iopub.status.busy":"2024-09-28T16:25:33.387461Z","iopub.execute_input":"2024-09-28T16:25:33.38805Z","iopub.status.idle":"2024-09-28T16:25:38.161823Z","shell.execute_reply.started":"2024-09-28T16:25:33.388014Z","shell.execute_reply":"2024-09-28T16:25:38.160841Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Modelling","metadata":{"papermill":{"duration":0.076791,"end_time":"2024-08-03T13:00:09.056736","exception":false,"start_time":"2024-08-03T13:00:08.979945","status":"completed"},"tags":[]}},{"cell_type":"code","source":"from sklearn.linear_model import LinearRegression\nfrom sklearn.linear_model import Ridge\n\nlr_model = LinearRegression()\nlr_model.fit(train, train_labels)\nlr_pred = lr_model.predict(train)\n\noof_pred = lr_pred\nsigma_pred = mean_squared_error(train_labels, oof_pred, squared=False)\nprint(f\"R2 score: {r2_score(train_labels, oof_pred):.3f}\")\nprint(f\"Root mean squared error: {sigma_pred:.6f}\")","metadata":{"execution":{"iopub.status.busy":"2024-09-28T16:32:24.784754Z","iopub.execute_input":"2024-09-28T16:32:24.785558Z","iopub.status.idle":"2024-09-28T16:32:29.501756Z","shell.execute_reply.started":"2024-09-28T16:32:24.785522Z","shell.execute_reply":"2024-09-28T16:32:29.500481Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"oof_df = postprocessing(oof_pred, train_adc_info.index, sigma_pred)\ngll_score = competition_score(train_labels.copy().reset_index(),\n                              oof_df.copy().reset_index(),\n                              naive_mean=train_labels.values.mean(),\n                              naive_sigma=train_labels.values.std(),\n                              sigma_true=0.000003)\nprint(f\"Estimated competition score: {gll_score:.4f}\")","metadata":{"_kg_hide-input":true,"papermill":{"duration":0.029243,"end_time":"2024-08-03T12:42:17.162838","exception":false,"start_time":"2024-08-03T12:42:17.133595","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-09-28T16:32:29.50447Z","iopub.execute_input":"2024-09-28T16:32:29.505404Z","iopub.status.idle":"2024-09-28T16:32:29.600471Z","shell.execute_reply.started":"2024-09-28T16:32:29.505358Z","shell.execute_reply":"2024-09-28T16:32:29.599543Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with open('lr_model.pickle', 'wb') as f:\n    pickle.dump(lr_model, f)\nwith open('sigma_pred.pickle', 'wb') as f:\n    pickle.dump(sigma_pred, f)","metadata":{"execution":{"iopub.status.busy":"2024-09-26T19:26:27.84804Z","iopub.execute_input":"2024-09-26T19:26:27.84845Z","iopub.status.idle":"2024-09-26T19:26:27.949972Z","shell.execute_reply.started":"2024-09-26T19:26:27.848416Z","shell.execute_reply":"2024-09-26T19:26:27.94919Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission","metadata":{"papermill":{"duration":0.077479,"end_time":"2024-08-03T13:00:10.305108","exception":false,"start_time":"2024-08-03T13:00:10.227629","status":"completed"},"tags":[]}},{"cell_type":"code","source":"test_adc_info = pd.read_csv(f'{PATH}/test_adc_info.csv', index_col='planet_id')\nsample_submission = pd.read_csv(f'{PATH}/sample_submission.csv', index_col='planet_id')\n\nf_raw_test = f_read_and_preprocess('test', sample_submission.index, 'FGS1')\na_raw_test = f_read_and_preprocess('test', sample_submission.index)\n\ntest = feature_engineering(f_raw_test, a_raw_test, test_adc_info)\ntest = test.iloc[: , :-1]\n\nwith open('lr_model.pickle', 'rb') as f:\n    lr_model = pickle.load(f)\nwith open('sigma_pred.pickle', 'rb') as f:\n    sigma_pred = pickle.load(f)\n\nlr_pred_test = lr_model.predict(test)\n\ntest_pred = lr_pred_test\n\nsub_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)))\nsub_df.to_csv('submission.csv')","metadata":{"papermill":{"duration":1.822768,"end_time":"2024-08-03T13:00:12.204995","exception":false,"start_time":"2024-08-03T13:00:10.382227","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-09-26T19:27:20.310865Z","iopub.execute_input":"2024-09-26T19:27:20.311583Z","iopub.status.idle":"2024-09-26T19:27:23.857745Z","shell.execute_reply.started":"2024-09-26T19:27:20.311546Z","shell.execute_reply":"2024-09-26T19:27:23.856337Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}