{"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":"# preface \nI was inspired by @AmbrosM idea https://www.kaggle.com/competitions/ariel-data-challenge-2024/discussion/530152#2969648\n\nThe idea of ​​this approach is to predict the average answer for each test sample. Instead of building models, we will search for the correct answer experimentally. We will select for each spectrum such a multiplier that the transit part multiplied by it will \"line up\" with the other points. The line can be either a straight line or a polynomial up to the 5th degree. For selection, we will use the Nelder-Mead method. The found multiplication factor minus one is our answer.\n\nThere are some changes in data preparation here.\n* dt for dark calibration changed in favor https://www.kaggle.com/code/gordonyip/update-calibrating-and-binning-astronomical-data/comments#2964759\n* signal clipped to zero due this https://www.kaggle.com/competitions/ariel-data-challenge-2024/discussion/530247#2970709\n* I masked hot and dead pixels with NaN in flat and averaging through spatial dimension.","metadata":{}},{"cell_type":"markdown","source":"### **upd for version 2**: it seems I forgot to remove **5x sigma** copied from somewhere. I wonder how this will affect the score.","metadata":{}},{"cell_type":"markdown","source":"### preprocess the data","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-08-27T21:12:59.527415Z","iopub.execute_input":"2024-08-27T21:12:59.527861Z","iopub.status.idle":"2024-08-27T21:12:59.534303Z","shell.execute_reply.started":"2024-08-27T21:12:59.527818Z","shell.execute_reply":"2024-08-27T21:12:59.533041Z"},"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-08-27T21:12:59.5364Z","iopub.execute_input":"2024-08-27T21:12:59.536891Z","iopub.status.idle":"2024-08-27T21:12:59.565766Z","shell.execute_reply.started":"2024-08-27T21:12:59.53685Z","shell.execute_reply":"2024-08-27T21:12:59.564668Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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 = 39, 321\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, 282], \"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] #11250 * 32 * 282\n            #dt = axis_info['AIRS-CH0-integration_time'].dropna().values\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        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        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([preproc('test', test_adc_info, \"FGS1\", 25*12), preproc('test', test_adc_info, \"AIRS-CH0\", 25)], axis=2)","metadata":{"execution":{"iopub.status.busy":"2024-08-27T21:12:59.567898Z","iopub.execute_input":"2024-08-27T21:12:59.568281Z","iopub.status.idle":"2024-08-27T21:13:09.429149Z","shell.execute_reply.started":"2024-08-27T21:12:59.568251Z","shell.execute_reply":"2024-08-27T21:13:09.427848Z"},"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(50//2,150//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    best_drop = 0\n    for i in range(200//2,250//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, 7):\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, 6):\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    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-08-28T06:13:21.523111Z","iopub.execute_input":"2024-08-28T06:13:21.523501Z","iopub.status.idle":"2024-08-28T06:13:21.862027Z","shell.execute_reply.started":"2024-08-28T06:13:21.523472Z","shell.execute_reply":"2024-08-28T06:13:21.860458Z"},"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-08-27T21:13:09.49489Z","iopub.execute_input":"2024-08-27T21:13:09.495295Z","iopub.status.idle":"2024-08-27T21:13:09.820027Z","shell.execute_reply.started":"2024-08-27T21:13:09.49526Z","shell.execute_reply":"2024-08-27T21:13:09.818647Z"},"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.0001)\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-08-27T21:13:09.821627Z","iopub.execute_input":"2024-08-27T21:13:09.822379Z","iopub.status.idle":"2024-08-27T21:13:09.848749Z","shell.execute_reply.started":"2024-08-27T21:13:09.822337Z","shell.execute_reply":"2024-08-27T21:13:09.847702Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission","metadata":{"execution":{"iopub.status.busy":"2024-08-27T21:13:09.850243Z","iopub.execute_input":"2024-08-27T21:13:09.850898Z","iopub.status.idle":"2024-08-27T21:13:09.874214Z","shell.execute_reply.started":"2024-08-27T21:13:09.850858Z","shell.execute_reply":"2024-08-27T21:13:09.873047Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}