{"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":196522559,"sourceType":"kernelVersion"},{"sourceId":198267231,"sourceType":"kernelVersion"},{"sourceId":203621335,"sourceType":"kernelVersion"}],"dockerImageVersionId":30761,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"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\nfrom functools import partial\nfrom astropy.stats import sigma_clip\n\nfrom scipy.ndimage import gaussian_filter\nfrom scipy.signal import savgol_filter\nimport scipy\nfrom multiprocessing import Pool\nfrom astropy.convolution import convolve, Box1DKernel\ndef smooth_data(data, window_size):\n    return savgol_filter(data, window_size, 3)  ","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-10-22T13:57:13.547881Z","iopub.execute_input":"2024-10-22T13:57:13.548444Z","iopub.status.idle":"2024-10-22T13:57:14.97239Z","shell.execute_reply.started":"2024-10-22T13:57:13.548392Z","shell.execute_reply":"2024-10-22T13:57:14.970351Z"},"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')\n\nvert_flux = np.load('/kaggle/input/exoplanets-signal-processing-v5/vert_flux.npy')\nwl_symmetry_correction = vert_flux[:,:].mean(0)\nreverse_corr = wl_symmetry_correction[::-1]/wl_symmetry_correction","metadata":{"execution":{"iopub.status.busy":"2024-10-22T13:52:05.99838Z","iopub.execute_input":"2024-10-22T13:52:05.999089Z","iopub.status.idle":"2024-10-22T13:52:06.229768Z","shell.execute_reply.started":"2024-10-22T13:52:05.999045Z","shell.execute_reply":"2024-10-22T13:52:06.228303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Calibration","metadata":{}},{"cell_type":"code","source":"%%writefile process.py\nimport 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\nfrom functools import partial\nfrom astropy.stats import sigma_clip\n\nfrom scipy.ndimage import gaussian_filter\nfrom scipy.signal import savgol_filter\nimport scipy\nfrom multiprocessing import Pool\nfrom astropy.convolution import convolve, Box1DKernel\ndef smooth_data(data, window_size):\n    return savgol_filter(data, window_size, 3)  \ndataset = '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')\n\n\nvert_flux = np.load('/kaggle/input/exoplanets-signal-processing-v5/vert_flux.npy')\nwl_symmetry_correction = vert_flux[:,:].mean(0)\nreverse_corr = wl_symmetry_correction[::-1]/wl_symmetry_correction\n\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\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 wavelength_interpolate(signal, distance = 1, fill_boundaries = False):\n    average = (signal[:,:,:-distance*2]+signal[:,:,distance*2:])/2 #fix dead sensors using interpolation across wavelength\n    if fill_boundaries:\n        return np.concatenate([signal[:,:,1:distance+1],average,signal[:,:,-distance-1:-1]], axis = 2)\n   \n    else:\n        return np.concatenate([signal[:,:,0:distance],average,signal[:,:,-distance:]], axis = 2)\n           \ndef flux_dist_interpolate(signal, distance = 1, fill_boundaries = False):\n    average = (signal[:,:-distance*2]+signal[:,distance*2:])/2 #fix dead sensors using interpolation across wavelength\n    if fill_boundaries:\n        return np.concatenate([signal[:,1:distance+1],average,signal[:,-distance-1:-1]], axis = 1)\n   \n    else:\n        return np.concatenate([signal[:,0:distance],average,signal[:,-distance:]], axis = 1)        \n\ndef process_planet(i, planet_id):\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] \n        signal = signal.reshape(sensor_sizes_dict[sensor][0][0],144) \n\n    if sensor != \"FGS1\":\n        signal = signal[:,10:22,:] \n\n    cds_signal = (signal[1::2] - signal[0::2])\n    if sensor != \"FGS1\":\n        wl_flux_dist = np.nanmean(cds_signal[:,:,:]/cds_signal[:,:,:].sum(1)[:,None,:],0)\n        fix = flux_dist_interpolate(wl_flux_dist,distance = 1, fill_boundaries = False)\n        fixed_wl_flux_dist = np.where(np.isnan(wl_flux_dist), fix, wl_flux_dist)\n\n        for dist in range(2,10):\n            nan_filt = np.isnan(fixed_wl_flux_dist)\n            if nan_filt.sum() ==0:\n                break\n            fixed_wl_flux_dist = np.where(nan_filt,flux_dist_interpolate(fixed_wl_flux_dist,\n                                                                             distance = dist,fill_boundaries = 0),\n                                              fixed_wl_flux_dist)\n        fixed_wl_corr = fixed_wl_flux_dist[::-1,:]/fixed_wl_flux_dist\n        fixed_wl_corr = np.where(np.isnan(fixed_wl_corr),1,fixed_wl_corr)\n        fixed_wl_flux_dist = np.where(np.isnan(fixed_wl_flux_dist),vert_flux.mean(0),fixed_wl_flux_dist)\n        error_rate = np.isnan(cds_signal)\n        \n        nan_filt = np.isnan(cds_signal)\n        if nan_filt.sum() >0:\n            signal_flux_reverse = cds_signal[:,::-1,:]#*fixed_wl_corr\n            signal_flux_fill = (cds_signal[:,:,:-2]+cds_signal[:,:,2:])/2 #fix dead sensors using interpolation across wavelength\n            signal_flux_fill = np.concatenate([cds_signal[:,:,1,None],signal_flux_fill,cds_signal[:,:,-2,None]], axis = 2)\n            \n            cds_signal = np.where(nan_filt, (signal_flux_reverse+signal_flux_fill)/2, cds_signal) #fix dead sensors using symmetry\n        \n        \n        nan_filt = np.isnan(cds_signal)\n        if nan_filt.sum() >0:\n            signal_flux_reverse = cds_signal[:,::-1,:]#*fixed_wl_corr\n            cds_signal = np.where(nan_filt, signal_flux_reverse, cds_signal) #fix dead sensors using symmetry\n        \n#         nan_filt = np.isnan(cds_signal)\n#         if nan_filt.sum() >0:\n#             signal_flux_fill = (cds_signal[:,:,:-2]+cds_signal[:,:,2:])/2 #fix dead sensors using interpolation across wavelength\n#             signal_flux_fill = np.concatenate([cds_signal[:,:,0,None],signal_flux_fill,cds_signal[:,:,-1,None]], axis = 2)\n#             cds_signal = np.where(np.isnan(cds_signal), signal_flux_fill, cds_signal)\n#         nan_filt = np.isnan(cds_signal)\n        nan_filt = np.isnan(cds_signal)\n        if nan_filt.sum() >0:\n            for dist in range(1,10):\n                nan_filt = np.isnan(cds_signal)\n                if nan_filt.sum() ==0:\n                    break\n                cds_signal = np.where(nan_filt,wavelength_interpolate(cds_signal, distance = dist,fill_boundaries = 0),cds_signal)\n        \n        \n        nan_filt = np.isnan(cds_signal)\n        if nan_filt.sum() >0:\n            signal_flux_fill = (cds_signal[:,:,:-2]+cds_signal[:,:,2:])/2 #fix dead sensors using interpolation across wavelength\n            signal_flux_fill = np.concatenate([cds_signal[:,:,1,None],signal_flux_fill,cds_signal[:,:,-2,None]], axis = 2)\n            cds_signal = np.where(nan_filt, signal_flux_fill, cds_signal)\n        \n        \n        error_rate_scaled = (error_rate.mean(0)*fixed_wl_flux_dist).sum(0)\n        #smoothed = gaussian_filter(cds_signal,1,axes = (0,2)) #na's\n        #smoothed = savgol_filter(cds_signal,20,3,axis = 0)  #**** updates ****, smooth across the spectrometer's vertical\n    cds_signal = np.nanmean(cds_signal, axis=1)\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    if sensor != \"FGS1\":   \n        return cds_signal,error_rate_scaled\n    else:\n        return cds_signal\n\ndataset = f'{dataset}'\nsensor = \"FGS1\"\nbinning = 30*12\n\ncut_inf, cut_sup = 39, 321\nsensor_sizes_dict = {\"AIRS-CH0\":[[11250, 32, 356], [1, 32, cut_sup-cut_inf]], \"FGS1\":[[135000, 32, 32], [1, 32, 32]]}\nbinned_dict = {\"AIRS-CH0\":[11250 // binning // 2, 282], \"FGS1\":[135000 // binning // 2]}\nlinear_corr_dict = {\"AIRS-CH0\":(6, 32, 356), \"FGS1\":(6, 32, 32)}\nplanet_ids = adc_info.index\n\nfeats = []\nerrors = []\nwith Pool(4) as p:\n    feats = p.starmap(process_planet, [(i, planet_id) for i, planet_id in list(enumerate(planet_ids))])  \ncolor =  np.stack(feats)\n    \ndataset = f'{dataset}'\nsensor =  \"AIRS-CH0\"\nbinning = 30\n\n\ncut_inf, cut_sup = 39, 321\nsensor_sizes_dict = {\"AIRS-CH0\":[[11250, 32, 356], [1, 32, cut_sup-cut_inf]], \"FGS1\":[[135000, 32, 32], [1, 32, 32]]}\nbinned_dict = {\"AIRS-CH0\":[11250 // binning // 2, 282], \"FGS1\":[135000 // binning // 2]}\nlinear_corr_dict = {\"AIRS-CH0\":(6, 32, 356), \"FGS1\":(6, 32, 32)}\nplanet_ids = adc_info.index\n\nfeats = []\nerrors = []\nwith Pool(4) as p:\n    feats = p.starmap(process_planet, [(i, planet_id) for i, planet_id in list(enumerate(planet_ids))])  \nbins = [t[0] for t in feats]\nerror_rate = [t[1] for t in feats]\n#vert_flux = [t[2] for t in feats]\n    \nspectra = np.stack(bins)\nerror_rate = np.stack(error_rate)\nnp.save('error_rate.npy',error_rate)\nnp.save('data_train.npy', spectra)\nnp.save('data_train_color.npy', color)\n","metadata":{"execution":{"iopub.status.busy":"2024-10-22T13:52:08.340881Z","iopub.execute_input":"2024-10-22T13:52:08.341357Z","iopub.status.idle":"2024-10-22T13:52:22.64165Z","shell.execute_reply.started":"2024-10-22T13:52:08.341314Z","shell.execute_reply":"2024-10-22T13:52:22.639791Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!python process.py","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"color = np.load('data_train_color.npy')\n\nbinning = 12\nbinned = np.zeros((len(color),5625))\n\nfor j in range(color.shape[1] // binning):\n    binned[:,j] = color[:,j*binning:j*binning+binning].mean(axis=1) ","metadata":{"execution":{"iopub.status.busy":"2024-10-22T19:55:35.020399Z","iopub.execute_input":"2024-10-22T19:55:35.021134Z","iopub.status.idle":"2024-10-22T19:55:35.329495Z","shell.execute_reply.started":"2024-10-22T19:55:35.02108Z","shell.execute_reply":"2024-10-22T19:55:35.328251Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pre_train = np.concatenate([binned[:,:,None],\n                            np.load('data_train.npy')], axis=2)\npre_train = np.concatenate([pre_train[:,:,0,None],pre_train[:,:,::-1][:,:,:-1]], axis = 2)\n\nimport gc\ndel binned\ndel color\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2024-10-22T13:55:49.679545Z","iopub.execute_input":"2024-10-22T13:55:49.680041Z","iopub.status.idle":"2024-10-22T13:55:49.890594Z","shell.execute_reply.started":"2024-10-22T13:55:49.679991Z","shell.execute_reply":"2024-10-22T13:55:49.888799Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy.optimize import minimize\nfrom scipy import optimize\nfrom functools import partial\n\ndef smooth_data(data, window_size):\n    return savgol_filter(data, window_size, 3)  # window size 51, polynomial order 3\n\ndef objective(signal, deg, s, enhanced = False,region1 = None, transit = None, region2 = None):\n    \n    if enhanced:\n        x = region1 + transit + region2\n        y = signal[region1].tolist() + (signal[transit] * (1 + s)).tolist() + signal[region2].tolist()\n        \n    else:\n        delta = 240\n        x = list(range(0,p1-delta))+list(range(p1+delta,p2 - delta))+list(range(p2+delta,len(signal)))\n        y = signal[:p1-delta].tolist() + (signal[p1+delta:p2 - delta] * (1 + s)).tolist() + signal[p2+delta:].tolist()\n    \n    #print(len(x), len(y))\n    z = np.polyfit(x, y, deg=deg)\n    p = np.poly1d(z)\n    q = np.abs(p(x) - y).mean()\n    \n    return q\n\ndef objective2(signal, deg, s, enhanced = False,p1 = None, p2= None):\n    \n    if enhanced:\n        x = region1 + transit + region2\n        y = signal[region1].tolist() + (signal[transit] * (1 + s)).tolist() + signal[region2].tolist()\n        \n    else:\n        delta = 240\n        x = list(range(0,p1-delta))+list(range(p1+delta,p2 - delta))+list(range(p2+delta,len(signal)))\n        y = signal[:p1-delta].tolist() + (signal[p1+delta:p2 - delta] * (1 + s)).tolist() + signal[p2+delta:].tolist()\n    \n    #print(len(x), len(y))\n    z = np.polyfit(x, y, deg=deg)\n    p = np.poly1d(z)\n    q = np.abs(p(x) - y).mean()\n    \n    return q\n\n\ndef best_fit2(deg, s, signal, enhanced = False,region1 = None, transit = None, region2 = None):\n    if enhanced:\n        x = region1 + transit + region2\n        y = signal[region1].tolist() + (signal[transit] * (1 + s)).tolist() + signal[region2].tolist()\n        \n    else:\n        delta = 2\n        x = list(range(0,p1-delta))+list(range(p1+delta,p2 - delta))+list(range(p2+delta,len(signal)))\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=deg)\n    p = np.poly1d(z)\n    q = np.abs(p(x) - y).mean()\n    \n    return p(x),np.array(y)\n\ndef best_fit3(signal,deg, s, enhanced=False,region1 = None, transit = None, region2 = None, p1 = None, p2 = None):\n    if enhanced:\n        x = region1 + transit + region2\n        y = signal[region1].tolist() + (signal[transit] * (1 + s)).tolist() + signal[region2].tolist()\n        \n    else:\n        delta = 240\n        x = list(range(0,p1-delta))+list(range(p1+delta,p2 - delta))+list(range(p2+delta,len(signal)))\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=deg)\n    \n    return z\n\n\ndef phase_detector(signal):\n    window = 120\n    box_kernel = Box1DKernel(window)\n    #smoothed_signal = convolve(signal, box_kernel, boundary = 'extend')\n    smoothed_signal = smooth_data(signal, window)\n    difference = np.gradient(smoothed_signal)\n    #cut out boundaries\n    cutoff = int((window-1)/2)\n    \n    difference = convolve(difference, box_kernel, boundary = 'extend')\n    #print(difference)\n    #plt.plot(difference)\n    p1 = difference[cutoff:-cutoff].argmin()+cutoff #prevent errors with endpoints\n    p2 = difference[cutoff:-cutoff].argmax()+cutoff\n    return p1, p2\ndef phase_detection2(signal,p1,p2,pred, deg = 3):\n    z = best_fit3(signal,deg, pred, p1 = p1, p2 = p2)\n    p = np.poly1d(z)\n    px = p(np.array(range(len(signal))))\n    limit = 8*30\n    noise_var = (\n        np.abs(signal[list(range(p1-limit))+list(range(p2+limit,len(signal)))]-px[list(range(p1-limit))+list(range(p2+limit,len(signal)))]).mean()\n    )*1\n    #noise_var = np.abs(signal[list(range(p1-8))+list(range(p2+8,len(signal)))]-px[list(range(p1-8))+list(range(p2+8,len(signal)))]).mean()*3\n    for step in range(1,limit+1):\n        margin = np.abs(px[p1-step]-signal[p1-step])\n        if (margin<noise_var):\n            break\n    d1=step-1\n\n    for step in range(1,limit+1):\n        margin = np.abs(px[p1+step]-signal[p1+step]*(1+pred))\n        if (margin<noise_var):\n            break\n    d2 = step-1\n\n    for step in range(1,limit+1):\n        margin = np.abs(px[p2-step]-signal[p2-step]*(1+pred))\n        if (margin<noise_var):\n            break\n    d3=step-1\n\n    for step in range(1,limit+1):\n        margin = np.abs(px[p2+step]-signal[p2+step])\n        if (margin<noise_var):\n            break\n    d4 = step-1\n    return d1,d2,d3,d4\n","metadata":{"execution":{"iopub.status.busy":"2024-10-22T13:56:09.338367Z","iopub.execute_input":"2024-10-22T13:56:09.338885Z","iopub.status.idle":"2024-10-22T13:56:09.379496Z","shell.execute_reply.started":"2024-10-22T13:56:09.338841Z","shell.execute_reply":"2024-10-22T13:56:09.377945Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def fit_mean_curve(i):\n    signal = pre_train[i,:,:240].mean(axis=1)#phase detection on unscaled is more reliable, since higher wavelength less noise\n    p1,p2 = phase_detector(signal)\n    signal = smooth_data(signal, 121)\n\n\n    limit = 240\n    #     noise_var = np.std(signal[list(range(p1-limit))+list(range(p2+limit,len(signal)))])\n    #     signal = gaussian_filter(signal,sigma = noise_var)*2\n    r = minimize(\n                    partial(objective2,signal,4, p1 = p1, p2 = p2),\n                    [0.0001],\n                    method= 'Nelder-Mead',tol = 0.00001\n                      )\n    pred = r.final_simplex[0].mean()\n\n\n    d1,d2,d3,d4 = phase_detection2(signal,p1,p2,pred,4)\n    pad = 0\n    d1+=pad\n    d2+=pad\n    d3+=pad\n    d4+=pad\n    #phase.append([d1,d2,d3,d4])\n    region1 = list(range(p1-d1))\n    region1a =  list(range(p1-d1,p1+d2+1))\n    transit = list(range(p1+d2+1,p2-d3))\n    region2a = list(range(p2-d3,p2+d4+1))\n    region2 = list(range(p2+d4+1,len(signal)))\n    b1,b2,b3,b4 = p1-d1, p1+d2, p2-d3, p2+d4\n    #cp.append([b1,b2,b3,b4])\n\n    signal = 1000*(pre_train[i,:,:240]/pre_train[i,:,:240].max(0)).mean(axis = 1)#pre_train[i,:,1:].mean(axis=1)\n    orig_signal = signal.copy()\n\n\n    window = 31\n    box_kernel = Box1DKernel(window)\n    #signal = convolve(signal, box_kernel, boundary = 'extend')\n    signal[region1] = convolve(signal[region1], box_kernel, boundary = 'extend')\n    signal[transit] = convolve(signal[transit], box_kernel, boundary = 'extend')\n    signal[region2] = convolve(signal[region2], box_kernel, boundary = 'extend')\n\n    \n\n#     cutoff = int((window-1)/2)\n#     region1 = region1[:-cutoff]\n#     transit = transit[cutoff:-cutoff]\n#     region2 = region2[cutoff:]\n    \n\n#     noise_var = np.std(signal[region1 +transit+ region2])\n#     signal[region1] = gaussian_filter(signal[region1],sigma = noise_var)\n#     signal[transit] = gaussian_filter(signal[transit],sigma = noise_var)\n#     signal[region2] = gaussian_filter(signal[region2],sigma = noise_var)\n\n#     signal[region1] = savgol_filter(signal[region1],30,3)\n#     signal[transit] = savgol_filter(signal[transit],30,3)\n#     signal[region2] = savgol_filter(signal[region2],30,3)\n    #     signal[region1] = convolve(signal[region1], box_kernel, boundary = 'extend')\n    #     signal[transit] = convolve(signal[transit], box_kernel, boundary = 'extend')\n    #     signal[region2] = convolve(signal[region2], box_kernel, boundary = 'extend')\n\n    ss = []\n    ls = []\n    r2 = []\n    mae = []\n    mse = []\n    orig_mae = []\n    linearity = []\n    dl_corr = []\n    for j in range(6):\n        r = minimize(\n                    partial(objective,signal,j,enhanced=True,\n                            region1 = region1, \n                            transit = transit, \n                            region2 = region2),\n                    [pred],\n                    method= 'Nelder-Mead', tol = 1e-5\n                      )\n        s = r.final_simplex[0].mean()\n        ls.append(r.fun)\n        ss.append(s)\n        yhat,y = best_fit2(j, s, signal,enhanced=True,region1 = region1, \n                            transit = transit, \n                            region2 = region2)\n        sse = np.abs(yhat-y)\n        r2.append(1-(sse**2).sum()/np.square((y-y.mean())).sum())\n        mse.append((sse**2).mean())\n        x,y = best_fit2(j, s, orig_signal,enhanced=True,region1 = region1, \n                            transit = transit, \n                            region2 = region2)\n        #z = best_fit3(j, s, enhanced=False)\n        #x = np.poly1d(z)(region1 + transit + region2)\n        #y = np.concatenate([orig_signal[region1],orig_signal[transit]*(1+s),orig_signal[region2]],axis = 0)\n        orig_mae.append(np.abs((x-y)/x).mean())\n\n        x = np.array(range(len(signal)))\n        z = best_fit3(signal,j, s, enhanced=True,region1 = region1, \n                            transit = transit, \n                            region2 = region2)\n        y = np.poly1d(z)(x)\n        z = np.polyfit(x, y, deg=1)\n        p = np.poly1d(z)\n        q = np.abs((p(x) - y)/p(x)).mean()\n        linearity.append(q)\n\n        z = best_fit3(signal,4, pred, enhanced = True,region1 = region1, \n                            transit = transit, \n                            region2 = region2)\n        p = np.poly1d(z)\n        px = p(np.array(range(len(signal))))\n        q_est = signal[transit]/px[transit]\n        z = np.polyfit(transit, q_est, deg=2)\n        p = np.poly1d(z)\n        z2 = np.polyfit([transit[0],transit[-1]], [p(transit[0]),p(transit[-1])], deg=1)\n        lin_corr = np.poly1d(z2)\n\n        corrected_ld = p(transit)-lin_corr(transit)+ min([p(transit[0]),p(transit[-1])])\n        transit_mid = (len(transit)+1)/2\n        if (len(transit)+1)%2:\n            mid_value = (corrected_ld[int(transit_mid)]+corrected_ld[int(transit_mid+.6)])/2\n        else:\n            mid_value = corrected_ld[int(transit_mid)]\n        dl_corr.append(corrected_ld.mean()/mid_value)\n\n    return ss+ls+r2+mse+orig_mae+linearity+dl_corr+[p1,p2,d1,d2,d3,d4]","metadata":{"execution":{"iopub.status.busy":"2024-10-22T13:56:32.329113Z","iopub.execute_input":"2024-10-22T13:56:32.329589Z","iopub.status.idle":"2024-10-22T13:56:32.363523Z","shell.execute_reply.started":"2024-10-22T13:56:32.329539Z","shell.execute_reply":"2024-10-22T13:56:32.361528Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"results = []\n\nwith Pool(4) as p:\n    results = p.map(fit_mean_curve, [i for i in range(len(pre_train))])\n\ndata2 = pd.DataFrame(results, columns = [f\"s{i}\" for i in range(6)]+\\\n                    [f\"mae{i}\" for i in range(6)]+[f\"r2_{i}\" for i in range(6)]+\\\n                    [f\"mse{i}\" for i in range(6)]+\\\n                   [f\"orig_mae{i}\" for i in range(6)]+[f\"linearity{i}\" for i in range(6)]+\\\n                   [f\"dl{i}\" for i in range(6)] + ['p1','p2','d1','d2','d3','d4'])","metadata":{"execution":{"iopub.status.busy":"2024-10-22T13:57:44.661483Z","iopub.execute_input":"2024-10-22T13:57:44.662192Z","iopub.status.idle":"2024-10-22T13:57:45.437689Z","shell.execute_reply.started":"2024-10-22T13:57:44.662131Z","shell.execute_reply":"2024-10-22T13:57:45.435499Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data2['45'] = data2['mse4']/data2['mse5']-1\ndata2['24'] = data2['mse2']/data2['mse4']-1\ndata2['34'] = data2['mse3']/data2['mse4']-1\ndata2['23'] = data2['mse2']/data2['mse3']-1\ndata2['12'] = data2['mse1']/data2['mse2']-1","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"s_points = data2[['p1','p2']].values/30\nphase = data2[['d1','d2','d3','d4']].values/30\nunbinned_phase = phase.copy()","metadata":{"execution":{"iopub.status.busy":"2024-10-22T14:03:08.085823Z","iopub.execute_input":"2024-10-22T14:03:08.087178Z","iopub.status.idle":"2024-10-22T14:03:08.096923Z","shell.execute_reply.started":"2024-10-22T14:03:08.087103Z","shell.execute_reply":"2024-10-22T14:03:08.095383Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"binning = 30\nbinned_pre_train = np.zeros((len(pre_train),5625//binning,283))\nfor k in range(5625 // binning):\n    binned_pre_train[:,k] = np.mean(pre_train[:,k*binning:k*binning+binning],1)\n    ","metadata":{"execution":{"iopub.status.busy":"2024-10-22T13:58:48.100705Z","iopub.execute_input":"2024-10-22T13:58:48.101201Z","iopub.status.idle":"2024-10-22T13:58:48.124194Z","shell.execute_reply.started":"2024-10-22T13:58:48.101147Z","shell.execute_reply":"2024-10-22T13:58:48.12279Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pre_train = binned_pre_train\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2024-10-22T13:59:30.68727Z","iopub.execute_input":"2024-10-22T13:59:30.687712Z","iopub.status.idle":"2024-10-22T13:59:30.941464Z","shell.execute_reply.started":"2024-10-22T13:59:30.687671Z","shell.execute_reply":"2024-10-22T13:59:30.940028Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# pre_train = np.load('/kaggle/input/exoplanets-signal-processing/data_train.npy')\n# pre_train = np.concatenate([pre_train[:,:,0,None],pre_train[:,:,::-1][:,:,:-1]], axis = 2)\n# pre_train = pre_train[:200]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Modelization","metadata":{}},{"cell_type":"code","source":"def objective(signal,deg, s, enhanced = False):\n    \n    if enhanced:\n        x = region1 + transit + region2\n        y = signal[region1].tolist() + (signal[transit] * (1 + s)).tolist() + signal[region2].tolist()\n        \n    else:\n        delta = 2\n        x = list(range(0,p1-delta))+list(range(p1+delta,p2 - delta))+list(range(p2+delta,len(signal)))\n        y = signal[:p1-delta].tolist() + (signal[p1+delta:p2 - delta] * (1 + s)).tolist() + signal[p2+delta:].tolist()\n    \n    #print(len(x), len(y))\n    z = np.polyfit(x, y, deg=deg)\n    p = np.poly1d(z)\n    q = np.abs(p(x) - y).mean()\n    \n    return q\n# def phase_detector_old(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[1:]) + 1 #prevent errors with endpoints\n#     phase2 = np.argmax(first_derivative2[:-1]) + MIN\n\n#     return phase1, phase2\ndef phase_detector(signal):\n    \n    window = 4\n    box_kernel = Box1DKernel(window)\n    #smoothed_signal = convolve(signal, box_kernel, boundary = 'extend')\n    #smoothed_signal = smooth_data(signal, window)\n    difference = np.gradient(signal)\n    #cut out boundaries\n    cutoff = int((window)/2)\n    \n    difference = convolve(difference, box_kernel, boundary = 'extend')\n    #print(difference)\n    #plt.plot(difference)\n    p1 = difference[cutoff:-cutoff].argmin()+cutoff #prevent errors with endpoints\n    p2 = difference[cutoff:-cutoff].argmax()+cutoff\n    return p1, p2\n\n    return phase1, phase2\n\ndef best_fit(deg, s):\n    i = deg\n    delta = 2\n    x = list(range(0,p1-delta))+list(range(p1+delta,p2 - delta))+list(range(p2+delta,len(signal)))\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    return x,p(x)\n\ndef best_fit2(deg, s, signal, enhanced):\n    if enhanced:\n        x = region1 + transit + region2\n        y = signal[region1].tolist() + (signal[transit] * (1 + s)).tolist() + signal[region2].tolist()\n        \n    else:\n        delta = 2\n        x = list(range(0,p1-delta))+list(range(p1+delta,p2 - delta))+list(range(p2+delta,len(signal)))\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=deg)\n    p = np.poly1d(z)\n    q = np.abs(p(x) - y).mean()\n    \n    return p(x),np.array(y)\n\ndef best_fit3(deg, s, enhanced=False):\n    if enhanced:\n        x = region1 + transit + region2\n        y = signal[region1].tolist() + (signal[transit] * (1 + s)).tolist() + signal[region2].tolist()\n        \n    else:\n        delta = 2\n        x = list(range(0,p1-delta))+list(range(p1+delta,p2 - delta))+list(range(p2+delta,len(signal)))\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=deg)\n    \n    return z\n\ndef phase_detection2(signal,p1,p2,pred, deg = 3):\n    z = best_fit3(deg, pred)\n    p = np.poly1d(z)\n    px = p(np.array(range(187)))\n    noise_var = np.sqrt(\n        np.square(signal[list(range(p1-8))+list(range(p2+8,len(signal)))]-px[list(range(p1-8))+list(range(p2+8,len(signal)))]).mean()\n    )*3\n    #noise_var = np.abs(signal[list(range(p1-8))+list(range(p2+8,len(signal)))]-px[list(range(p1-8))+list(range(p2+8,len(signal)))]).mean()*3\n    for step in [1,2,3,4,5,6,7,8]:\n        margin = np.abs(px[p1-step]-signal[p1-step])\n        if (margin<noise_var):\n            break\n    d1=step-1\n\n    for step in [1,2,3,4,5,6,7,8]:\n        margin = np.abs(px[p1+step]-signal[p1+step]*(1+pred))\n        if (margin<noise_var):\n            break\n    d2 = step-1\n\n    for step in [1,2,3,4,5,6,7,8]:\n        margin = np.abs(px[p2-step]-signal[p2-step]*(1+pred))\n        if (margin<noise_var):\n            break\n    d3=step-1\n\n    for step in [1,2,3,4,5,6,7,8]:\n        margin = np.abs(px[p2+step]-signal[p2+step])\n        if (margin<noise_var):\n            break\n    d4 = step-1\n    return d1,d2,d3,d4\n","metadata":{"execution":{"iopub.status.busy":"2024-10-22T14:03:15.064226Z","iopub.execute_input":"2024-10-22T14:03:15.065774Z","iopub.status.idle":"2024-10-22T14:03:15.103369Z","shell.execute_reply.started":"2024-10-22T14:03:15.065719Z","shell.execute_reply":"2024-10-22T14:03:15.101476Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"results = []\nphase = []\nfor i in range(len(pre_train)):\n\n    signal = pre_train[i,:,:240].mean(axis=1)#phase detection on unscaled is more reliable #phase_detector(signal)\n#     p1,p2 = s_points[i]\n#     p1 = int(np.round(p1))\n#     p2 = int(np.round(p2))\n    p1,p2 = phase_detector(signal)\n#    p1,p2 = phase_detector(signal)\n    #if p2-p1 <20: assert False\n    \n    noise_var = np.std(signal[list(range(p1-8))+list(range(p2+8,len(signal)))])\n    signal = gaussian_filter(signal,sigma = noise_var)\n    r = minimize(\n                    partial(objective,signal,3),\n                    [0.0001],\n                    method= 'Nelder-Mead',tol = 0.00001\n                      )\n    pred = r.final_simplex[0].mean()\n    \n    signal = pre_train[i,:,:240].mean(axis=1)\n    d1,d2,d3,d4 = phase_detection2(signal,p1,p2,pred)\n\n    phase.append([d1,d2,d3,d4])\n\n    region1 = list(range(p1-d1))\n    region1a =  list(range(p1-d1,p1+d2+1))\n    transit = list(range(p1+d2+1,p2-d3))\n    region2a = list(range(p2-d3,p2+d4+1))\n    region2 = list(range(p2+d4+1,len(signal)))\n    \n    signal = 1000*(pre_train[i,:,:240]/pre_train[i,:,:240].max(0)).mean(axis = 1)#pre_train[i,:,1:].mean(axis=1)\n    orig_signal = signal.copy()\n    #signal = 1000*(pre_train[i,:,:]/np.concatenate([pre_train[i,:p1-8,:],pre_train[i,p2+8:,:]],axis = 0).mean(0)).mean(axis = 1)\n    #signal = 1000*(pre_train[i,:,:]/(0.5*pre_train[i,p1-8,:]+0.5*pre_train[i,p2+8,:])).mean(axis = 1)\n    noise_var = np.std(signal[region1 +transit+ region2])\n    #signal = gaussian_filter(signal,sigma = noise_var)\n    signal[region1] = gaussian_filter(signal[region1],sigma = noise_var)\n    signal[transit] = gaussian_filter(signal[transit],sigma = noise_var)\n    signal[region2] = gaussian_filter(signal[region2],sigma = noise_var)\n        \n    ss = []\n    ls = []\n    r2 = []\n    mae = []\n    mse = []\n    orig_mae = []\n    for j in range(6):\n        r = minimize(\n                    partial(objective,signal,j,enhanced=True),\n                    [pred],\n                    method= 'Nelder-Mead', tol = 1e-5\n                      )\n        s = r.final_simplex[0].mean()\n        ls.append(r.fun)\n        ss.append(s)\n        yhat,y = best_fit2(j, s, signal,enhanced=True)\n        sse = np.abs(yhat-y)\n        r2.append(1-(sse**2).sum()/np.square((y-y.mean())).sum())\n        mse.append((sse**2).mean())\n        x,y = best_fit2(j, s, orig_signal,enhanced=True)\n        orig_mae.append(np.abs((x-y)/x).mean())\n    results.append(ss+ls+r2+mse+orig_mae)\n    \ndata = pd.DataFrame(results, columns = [f\"s{i}\" for i in range(6)]+[f\"mae{i}\" for i in range(6)]+[f\"r2_{i}\" for i in range(6)]+[f\"mse{i}\" for i in range(6)]+[f\"orig_mae{i}\" for i in range(6)])\nphase = np.array(phase)\n\ndata['45'] = data['mse4']/data['mse5']-1\ndata['24'] = data['mse2']/data['mse4']-1\ndata['34'] = data['mse3']/data['mse4']-1\ndata['23'] = data['mse2']/data['mse3']-1\ndata['12'] = data['mse1']/data['mse2']-1","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#smoothing by wl\npoly = {}\ndef fit_wl(wl, smoothing = 30):\n    split = 0.1\n    \n    orig_signal = 1000*(pre_train[index2,:,max(wl-smoothing,0):wl+smoothing+1].mean(1)/pre_train[index2,:,max(wl-smoothing,0):wl+smoothing+1].mean(1).max(0))#pre_train[index2,:,1:].mean(axis=1)\n    noise_var = np.std(orig_signal[region1 + transit+ region2])*filter_scale\n    signal = orig_signal.copy()#*wl_scale\n    signal[region1] = gaussian_filter(signal[region1],sigma = noise_var)\n    signal[transit] = gaussian_filter(signal[transit],sigma = noise_var)\n    signal[region2] = gaussian_filter(signal[region2],sigma = noise_var)\n\n    r = minimize(\n                partial(objective,signal,j, enhanced= True),\n                [last_prediction],\n                method= 'Nelder-Mead', tol = 1e-5\n                  )\n    s = r.final_simplex[0].mean()\n    ls_wl = r.fun\n    ss_wl = max(s,0.00001)\n    x,y = best_fit2(j,s,signal, enhanced= True)\n    mae_scaled = (np.abs(x-y)/x).mean()\n    \n    x,y = best_fit2(j,s,orig_signal, enhanced= True)\n    orig_mae = (np.abs(x-y)/x).mean()\n    #x,y = best_fit2(j,s,signal, enhanced= True)\n    x,y = best_fit2(j,ss[-1],signal, enhanced= True)\n    mean_mae = np.abs(x-y).mean()\n    mean_mae_scaled = (np.abs(x-y)/x).mean()\n    ss_wl_30 = ss_wl\n    orig_mae30 = orig_mae\n    mean_mae30 = mean_mae\n    mae30 = ls_wl\n#                 if s/orig_mae/10 < 1.6:\n#                     #print(wl,s, orig_mae,smoothing, s/orig_mae/10)\n#                     break\n    previous_result = [ss_wl,ls_wl,orig_mae,mean_mae,smoothing,mae_scaled,mean_mae_scaled,ss_wl_30,orig_mae30,mean_mae30,mae30]\n    if orig_mae*1000 >split:\n        return previous_result\n    \n    for smoothing in [25,20,15,10,5,2]:\n\n        orig_signal = 1000*(pre_train[index2,:,max(wl-smoothing,0):wl+smoothing+1].mean(1)/pre_train[index2,:,max(wl-smoothing,0):wl+smoothing+1].mean(1).max(0))#pre_train[index2,:,1:].mean(axis=1)\n        noise_var = np.std(orig_signal[region1 + transit+ region2])*filter_scale\n        signal = orig_signal.copy()#*wl_scale\n        signal[region1] = gaussian_filter(signal[region1],sigma = noise_var)\n        signal[transit] = gaussian_filter(signal[transit],sigma = noise_var)\n        signal[region2] = gaussian_filter(signal[region2],sigma = noise_var)\n\n        r = minimize(\n                    partial(objective,signal,j, enhanced= True),\n                    [last_prediction],\n                    method= 'Nelder-Mead', tol = 1e-5\n                      )\n        s = r.final_simplex[0].mean()\n        ls_wl = r.fun\n        ss_wl = max(s,0.00001)\n        x,y = best_fit2(j,s,signal, enhanced= True)\n        mae_scaled = (np.abs(x-y)/x).mean()\n\n        x,y = best_fit2(j,s,orig_signal, enhanced= True)\n        orig_mae = (np.abs(x-y)/x).mean()\n        #x,y = best_fit2(j,s,signal, enhanced= True)\n        x,y = best_fit2(j,ss[-1],signal, enhanced= True)\n        mean_mae = np.abs(x-y).mean()\n        mean_mae_scaled = (np.abs(x-y)/x).mean()\n        if orig_mae*1000 >split:\n            break\n        previous_result = [ss_wl,ls_wl,orig_mae,mean_mae,smoothing,mae_scaled,mean_mae_scaled,ss_wl_30,orig_mae30,mean_mae30,mae30]\n    if orig_mae > split + 0.03:\n        return previous_result\n    return [ss_wl,ls_wl,orig_mae,mean_mae,smoothing,mae_scaled,mean_mae_scaled,ss_wl_30,orig_mae30,mean_mae30,mae30]\n\nfor j in [2,3,4]:\n    smoothing = 30\n    smooth_window =20 #for wavelength curve\n    filter_scale = 1\n\n    poly[j] = {}\n    poly[j]['wl'] = []\n    poly[j]['mae'] = []\n    poly[j]['orig_mae'] = []\n    poly[j]['smoothed_wl'] = []\n    poly[j]['mean_mae']= []\n    poly[j]['smoothing']= []\n    poly[j]['wl_mae']= []\n    poly[j]['10_90']= []\n    poly[j]['mae_scaled'] = []\n    poly[j]['mean_mae_scaled'] = []\n    poly[j]['ss_wl_30'] = []\n    poly[j]['orig_mae30'] = []\n    poly[j]['mean_mae30'] = []\n    poly[j]['mae30'] = []\n    for index2 in tqdm(range(len(data))):\n\n        pred = data[f's{j}'].iloc[index2]#sub.iloc[index2,1]\n\n#         p1,p2 = s_points[index2]\n#         p1 = int(np.round(p1))\n#         p2 = int(np.round(p2))\n        signal = pre_train[index2,:,:].mean(axis=1)\n        p1,p2 = phase_detector(signal)\n\n        signal = 1000*(pre_train[index2,:,:]/pre_train[index2,:,:].max(0)).mean(axis = 1)#pre_train[index2,:,1:].mean(axis=1)\n\n\n        d1,d2,d3,d4 = phase[index2]\n\n#         d1 = int(np.round(d1))\n#         d2 = int(np.round(d2))\n#         d3 = int(np.round(d3))\n#         d4 = int(np.round(d4))\n\n        region1 = list(range(p1-d1))\n        region1a =  list(range(p1-d1,p1+d2+1))\n        transit = list(range(p1+d2+1,p2-d3))\n        region2a = list(range(p2-d3,p2+d4+1))\n        region2 = list(range(p2+d4+1,len(signal)))\n\n        noise_var = np.std(signal[region1 + transit+ region2])\n        signal[region1] = gaussian_filter(signal[region1],sigma = noise_var)\n        signal[transit] = gaussian_filter(signal[transit],sigma = noise_var)\n        signal[region2] = gaussian_filter(signal[region2],sigma = noise_var)\n\n\n        ss = []\n        ls = []\n        tls = []\n\n        #fit mean curve\n        r = minimize(\n                    partial(objective,signal,j,enhanced= True),\n                    [pred],\n                    method= 'Nelder-Mead',tol = 1e-5\n                      )\n        s = r.final_simplex[0].mean()\n        ls.append(r.fun)\n        ss.append(s)\n        z = best_fit3(j,s,enhanced=True)#np.polyfit(transit, signal[transit], deg=2)\n        p = np.poly1d(z)\n        wl_scale = p(range(187)).max()/p(range(187))\n\n        pred = s\n        last_prediction = s\n\n        results = []\n        with Pool(4) as p:\n            results = p.starmap(fit_wl, [(wl, smoothing) for wl in range(283)])  \n\n        results = np.array(results)\n        wl_mae = np.abs(results[:,0]-results[:,0].mean()).mean()\n        q90_q10 = np.quantile(results[:,0],0.9)- np.quantile(results[:,0],0.1)\n        #last check for extremities\n#         if q90_q10>0.0008:\n#             results = []\n#             smoothing = 10\n#             with Pool(4) as p:\n#                 results = p.starmap(fit_wl, [(wl, smoothing) for wl in range(283)])  \n#         elif q90_q10 >0.0005:\n#             results = []\n#             smoothing = 20\n#             with Pool(4) as p:\n#                 results = p.starmap(fit_wl, [(wl, smoothing) for wl in range(283)])  \n\n        results = np.array(results)\n        poly[j]['wl'].append(results[:,0])\n        poly[j]['mae'].append(results[:,1])\n        poly[j]['orig_mae'].append(results[:,2])\n        poly[j]['mean_mae'].append(results[:,3])\n        poly[j]['smoothing'].append(results[:,4])\n        poly[j]['mae_scaled'].append(results[:,5])\n        poly[j]['mean_mae_scaled'].append(results[:,6])\n        poly[j]['smoothed_wl'].append(smooth_data(results[:,0],smooth_window))\n        poly[j]['wl_mae'].append(wl_mae)\n        poly[j]['10_90'].append(q90_q10)\n        poly[j]['ss_wl_30'].append(results[:,7])\n        poly[j]['orig_mae30'].append(results[:,8])\n        poly[j]['mean_mae30'].append(results[:,9])\n        poly[j]['mae30'].append(results[:,10])\n    for key in poly[j]:\n        print(key)\n        poly[j][key] = np.stack(poly[j][key])\npoly2 = poly","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"limit =200\nmed2b = np.median(poly2[2]['smoothed_wl'][:,:limit],1)\nfilt = np.abs(med2b-data2['s2'].values)>1e-4\nmed2b[filt] = data2['s2'].values[filt]\nmed3b = np.median(poly2[3]['smoothed_wl'][:,:limit],1)\nfilt = np.abs(med3b-data2['s3'].values)>1e-4\nmed3b[filt] = data2['s3'].values[filt]\nmed4b = np.median(poly2[4]['smoothed_wl'][:,:limit],1)\nfilt = np.abs(med4b-data2['s4'].values)>1e-4\nmed4b[filt] = data2['s4'].values[filt]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"poly = {}","metadata":{"execution":{"iopub.status.busy":"2024-10-22T14:05:11.275464Z","iopub.execute_input":"2024-10-22T14:05:11.276814Z","iopub.status.idle":"2024-10-22T14:05:11.283389Z","shell.execute_reply.started":"2024-10-22T14:05:11.276757Z","shell.execute_reply":"2024-10-22T14:05:11.281341Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def fit_wl(wl, smoothing = 30):\n    orig_signal = 1000*(pre_train[index2,:,max(wl-smoothing,0):wl+smoothing+1].mean(1)/pre_train[index2,:,max(wl-smoothing,0):wl+smoothing+1].mean(1).max(0))#pre_train[index2,:,1:].mean(axis=1)\n    noise_var = np.std(orig_signal[region1 + transit+ region2])*filter_scale\n    signal = orig_signal.copy()#*wl_scale\n    signal[region1] = gaussian_filter(signal[region1],sigma = noise_var)\n    signal[transit] = gaussian_filter(signal[transit],sigma = noise_var)\n    signal[region2] = gaussian_filter(signal[region2],sigma = noise_var)\n\n    r = minimize(\n                partial(objective,signal,j, enhanced= True),\n                [last_prediction],\n                method= 'Nelder-Mead', tol = 1e-5\n                  )\n    s = r.final_simplex[0].mean()\n    ls_wl = r.fun\n    ss_wl = max(s,0.00001)\n    x,y = best_fit2(j,s,orig_signal, enhanced= True)\n    orig_mae = (np.abs(x-y)/x).mean()\n    #x,y = best_fit2(j,s,signal, enhanced= True)\n    x,y = best_fit2(j,ss[-1],signal, enhanced= True)\n    mean_mae = np.abs(x-y).mean()\n#                 if s/orig_mae/10 < 1.6:\n#                     #print(wl,s, orig_mae,smoothing, s/orig_mae/10)\n#                     break\n    return [ss_wl,ls_wl,orig_mae,mean_mae,smoothing]\n\n","metadata":{"execution":{"iopub.status.busy":"2024-10-22T14:03:25.033842Z","iopub.execute_input":"2024-10-22T14:03:25.034443Z","iopub.status.idle":"2024-10-22T14:03:25.049075Z","shell.execute_reply.started":"2024-10-22T14:03:25.034391Z","shell.execute_reply":"2024-10-22T14:03:25.047563Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for j in [2,3,4]:\n\n    smoothing = 30\n    smoothing_scale = np.round(np.linspace(10,smoothing,200)).astype(int)\n    smooth_window =20 #for wavelength curve\n    filter_scale = 1\n\n    poly[j] = {}\n    poly[j]['wl'] = []\n    poly[j]['mae'] = []\n    poly[j]['orig_mae'] = []\n    poly[j]['smoothed_wl'] = []\n    poly[j]['mean_mae']= []\n    poly[j]['smoothing']= []\n    poly[j]['wl_mae']= []\n    poly[j]['10_90']= []\n    for index2 in tqdm(range(len(data))):\n\n        pred = data[f's{j}'].iloc[index2]#sub.iloc[index2,1]\n\n#         p1,p2 = s_points[index2]\n#         p1 = int(np.round(p1))\n#         p2 = int(np.round(p2))\n        signal = pre_train[index2,:,:].mean(axis=1)\n        p1,p2 = phase_detector(signal)\n        \n        signal = 1000*(pre_train[index2,:,:]/pre_train[index2,:,:].max(0)).mean(axis = 1)#pre_train[index2,:,1:].mean(axis=1)\n        \n    \n        d1,d2,d3,d4 = phase[index2]\n        \n#         d1 = int(np.round(d1))\n#         d2 = int(np.round(d2))\n#         d3 = int(np.round(d3))\n#         d4 = int(np.round(d4))\n\n        region1 = list(range(p1-d1))\n        region1a =  list(range(p1-d1,p1+d2+1))\n        transit = list(range(p1+d2+1,p2-d3))\n        region2a = list(range(p2-d3,p2+d4+1))\n        region2 = list(range(p2+d4+1,len(signal)))\n\n        noise_var = np.std(signal[region1 + transit+ region2])\n        signal[region1] = gaussian_filter(signal[region1],sigma = noise_var)\n        signal[transit] = gaussian_filter(signal[transit],sigma = noise_var)\n        signal[region2] = gaussian_filter(signal[region2],sigma = noise_var)\n\n\n        ss = []\n        ls = []\n        tls = []\n\n        #fit mean curve\n        r = minimize(\n                    partial(objective,signal,j,enhanced= True),\n                    [pred],\n                    method= 'Nelder-Mead',tol = 1e-5\n                      )\n        s = r.final_simplex[0].mean()\n        ls.append(r.fun)\n        ss.append(s)\n        z = best_fit3(j,s,enhanced=True)#np.polyfit(transit, signal[transit], deg=2)\n        p = np.poly1d(z)\n        wl_scale = p(range(187)).max()/p(range(187))\n        \n        pred = s\n        last_prediction = s\n\n        results = []\n        smoothing = 30\n        with Pool(4) as p:\n            results = p.starmap(fit_wl, [(wl, smoothing) for wl in range(283)])  \n\n        results = np.array(results)\n        wl_mae = np.abs(results[:,0]-results[:,0].mean()).mean()\n        q90_q10 = np.quantile(results[:,0],0.9)- np.quantile(results[:,0],0.1)\n        #last check for extremities\n#         if q90_q10>0.0008:\n#             results = []\n#             smoothing = 10\n#             with Pool(4) as p:\n#                 results = p.starmap(fit_wl, [(wl, smoothing) for wl in range(283)])  \n#         elif q90_q10 >0.0005:\n#             results = []\n#             smoothing = 20\n#             with Pool(4) as p:\n#                 results = p.starmap(fit_wl, [(wl, smoothing) for wl in range(283)])  \n        if (q90_q10>np.exp(-7.5))|(wl_mae>np.exp(-9)):\n            results = []\n            smoothing = 10\n            with Pool(4) as p:\n                results = p.starmap(fit_wl, [(wl, smoothing) for wl in range(283)])  \n    \n        results = np.array(results)\n        poly[j]['wl'].append(results[:,0])\n        poly[j]['mae'].append(results[:,1])\n        poly[j]['orig_mae'].append(results[:,2])\n        poly[j]['mean_mae'].append(results[:,3])\n        poly[j]['smoothing'].append(results[:,4])\n        poly[j]['smoothed_wl'].append(smooth_data(results[:,0],smooth_window))\n        poly[j]['wl_mae'].append(wl_mae)\n        poly[j]['10_90'].append(q90_q10)\n    for key in poly[j]:\n        print(key)\n        poly[j][key] = np.stack(poly[j][key])","metadata":{"execution":{"iopub.status.busy":"2024-10-22T14:06:30.999417Z","iopub.execute_input":"2024-10-22T14:06:30.999895Z","iopub.status.idle":"2024-10-22T14:06:33.619796Z","shell.execute_reply.started":"2024-10-22T14:06:30.999854Z","shell.execute_reply":"2024-10-22T14:06:33.617493Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"improve_fit = {}\nimprove_fit['24'] = poly[2]['mae']/poly[4]['mae']-1\nimprove_fit['34'] = poly[3]['mae']/poly[4]['mae']-1\nimprove_fit['23'] = poly[2]['mae']/poly[3]['mae']-1\n# med2 = np.median(poly[2]['smoothed_wl'],1)\n# med3 = np.median(poly[3]['smoothed_wl'],1)\n# med4 = np.median(poly[4]['smoothed_wl'],1)\nlimit =200\nmed2 = np.median(poly[2]['smoothed_wl'][:,:limit],1)\nfilt = np.abs(med2-data['s2'].values)>1e-4\nmed2[filt] = data['s2'].values[filt]\nmed3 = np.median(poly[3]['smoothed_wl'][:,:limit],1)\nfilt = np.abs(med3-data['s3'].values)>1e-4\nmed3[filt] = data['s3'].values[filt]\nmed4 = np.median(poly[4]['smoothed_wl'][:,:limit],1)\nfilt = np.abs(med4-data['s4'].values)>1e-4\nmed4[filt] = data['s4'].values[filt]","metadata":{"execution":{"iopub.status.busy":"2024-10-22T14:06:36.585031Z","iopub.execute_input":"2024-10-22T14:06:36.585487Z","iopub.status.idle":"2024-10-22T14:06:36.596037Z","shell.execute_reply.started":"2024-10-22T14:06:36.585446Z","shell.execute_reply":"2024-10-22T14:06:36.594288Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"weights = [0,2,1,0,2,1]\npoly_choice = (poly[2]['smoothed_wl']*weights[0]+poly[3]['smoothed_wl']*weights[1]+poly[4]['wl']*weights[2]\\\n               +poly2[2]['smoothed_wl']*weights[3] + poly2[3]['smoothed_wl']*weights[4] +poly2[4]['smoothed_wl']*weights[5]\n              )/sum(weights)\n# for i in range(len(poly_choice)):\n#     poly_choice[i] = smooth_data(poly_choice[i],20)\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"error_rate = np.load(\"error_rate.npy\")\nerror_rate = np.concatenate([np.zeros([len(data),1]),error_rate[:,::-1]], axis = 1)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#no constant values at all\nj = 3\nmae_diff = np.abs(poly[j]['mean_mae']-poly[j]['mae'])/poly[j]['orig_mae']\n#mae_diff = np.abs(poly2[j]['mean_mae30']-poly2[j]['mae30'])/poly2[j]['orig_mae30']\n#mae_diff2 = np.abs(poly[2]['mean_mae']-poly[2]['mae'])/poly[2]['orig_mae']\n\nlimit =283\n\nfilt = (mae_diff <0)\nmean_pred = np.mean([\n    np.where(data2['24'].values<1, med2,med4),\n    np.where(data2['24'].values<1, med2b,med4b),\n    #np.where(data['24'].values<1, med2,med4),\n((med4)),\n    med4b,\n   np.where(data['24'].values<1, data['s2'].values/3+data['s3'].values/3+data['s4'].values/3,data['s4'].values),\n    np.where(data['24'].values<1, data['s2'].values/3+data['s3'].values/3+data['s4'].values/3,data['s4'].values),\nnp.where(data2['24'].values<1, data2['s2'].values/3+data2['s3'].values/3+data2['s4'].values/3,data2['s4'].values),\nnp.where(data2['24'].values<1, data2['s2'].values/3+data2['s3'].values/3+data2['s4'].values/3,data2['s4'].values),\n#np.where(data6['24'].values<.3, data6['s2'].values,data6['s4'].values),\n    #np.where(data3['24'].values<.2, data3['s2'].values,data3['s4'].values),\n#np.where(data5['24'].values<1, data5['s2'].values/3+data5['s3'].values/3+data5['s4'].values/3,data5['s4'].values)/2,   \n\n],0)\n\ncurve_choice = poly[j]['smoothed_wl']\ncurve_choice = poly_choice\nmean_lc = np.repeat(mean_pred,283).reshape((len(data),283))\nmean_shifted = curve_choice-curve_choice.mean(1)[:,None]+mean_lc\n#mean_shifted = poly_choice-poly_choice.mean(1)[:,None]+mean_lc\norig_mae30 = poly2[j]['orig_mae30']\nmae = np.abs((poly2[j]['ss_wl_30']-poly2[j]['ss_wl_30'].mean(1)[:,None])).mean(1)\n\n\ntarget_scale = poly[j]['orig_mae'].mean(1)*.7\neg = curve_choice.copy()\neg_diff = eg-eg.mean(1)[:,None]\neg_scaled = eg_diff/eg_diff.max(1)[:,None]*target_scale[:,None]\n\npred = np.where(filt,mean_lc + eg_scaled,\n                np.concatenate(\n                    [\n                        mean_shifted\n                       #curve_choice\n                     #poly_choice\n                        ,\n                        (mean_lc + eg_scaled)[:,limit:]\n                    ],axis = 1)\n               )\nbase_std = orig_mae30*.3+poly[j]['wl_mae'][:,None]*.6\n\n\nstd = np.where(filt,base_std*.9,\n               base_std*.9,\n\n              )\n\npercentage_filter = (mae_diff<55).mean(1) > 0.75\nstd[percentage_filter] = mae[percentage_filter][:,None]*0.65 +np.repeat(data['orig_mae4'].values[percentage_filter,None]*.35,283, axis = 1)\ntarget_scale = (std[percentage_filter])/2\neg = curve_choice[percentage_filter]\neg_diff = eg-eg.mean(1)[:,None]\neg_scaled = eg_diff/eg_diff.max(1)[:,None]*target_scale\npred[percentage_filter] = mean_lc[percentage_filter] + eg_scaled\nprint(filt.sum()/673/282, sum(percentage_filter))\n\n\npercentage_filter = (mae_diff<25).mean(1) > 0.75\nstd[percentage_filter] = mae[percentage_filter][:,None]*0.55 +np.repeat(data['orig_mae4'].values[percentage_filter,None]*.35,283, axis = 1)\ntarget_scale = (std[percentage_filter])/2\neg =curve_choice[percentage_filter]\neg_diff = eg-eg.mean(1)[:,None]\neg_scaled = eg_diff/eg_diff.max(1)[:,None]*target_scale\npred[percentage_filter] = mean_lc[percentage_filter]  +eg_scaled\nprint(filt.sum()/673/282, sum(percentage_filter))\n\npercentage_filter =( (mae_diff<15).mean(1) > 0.75 )\nstd[percentage_filter] = mae[percentage_filter][:,None]*0.45+np.repeat(data['orig_mae4'].values[percentage_filter,None]*.35,283, axis = 1)\ntarget_scale = (std[percentage_filter])/3\neg = curve_choice[percentage_filter]\neg_diff = eg-eg.mean(1)[:,None]\neg_scaled = eg_diff/eg_diff.max(1)[:,None]*target_scale\npred[percentage_filter] = mean_lc[percentage_filter] +eg_scaled\nprint(filt.sum()/673/282, sum(percentage_filter))\n\nstd = std*1.1\n#std = std*(1+error_rate*2)#np.where((error_rate>0.1), std*1.5, std)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"all_s = pred + 0.000004\n#sigma = deg4_mae/1000#np.ones_like(all_s) * 0.000130\nsigma = std","metadata":{"execution":{"iopub.status.busy":"2024-10-22T14:06:47.853426Z","iopub.execute_input":"2024-10-22T14:06:47.854299Z","iopub.status.idle":"2024-10-22T14:06:47.861574Z","shell.execute_reply.started":"2024-10-22T14:06:47.854238Z","shell.execute_reply":"2024-10-22T14:06:47.860111Z"},"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\npred = all_s.clip(0) \nsubmission = pd.DataFrame(np.concatenate([pred,sigma], axis=1), columns=ss.columns[1:])\nsubmission.index = adc_info.index\nsubmission.to_csv('submission.csv')\nsubmission\n","metadata":{"execution":{"iopub.status.busy":"2024-10-22T14:06:49.743308Z","iopub.execute_input":"2024-10-22T14:06:49.743849Z","iopub.status.idle":"2024-10-22T14:06:49.820295Z","shell.execute_reply.started":"2024-10-22T14:06:49.743805Z","shell.execute_reply":"2024-10-22T14:06:49.818803Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}