{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# It is an illustration of test sample denoising.\n\nFour days before the end of the competition we realized that the nonstationary noise test data was generated by taking a large single chunk of real detector noise data, adding some white noise, a random selection of a subset of time_ids, sorting them, assigning randomly selected time_stamps, selection of the frequency range, and adding the signal. The frequency also appeared to be shifted by an integer number of Hz. Therefore, noise is identical in multiple test files within overlapping frequency ranges. **This kernel illustrates how to perform sample denoising by matching noisy regions and subtructing one from another.** As a result a noise-free sample remains.\n\nUnfortunetly, not all examples overlap, overlap in nearly all cases is only partial, and the difference does not provide a solid evidence of the wave source. Two days before the end of the competitions we realized that nonstationary noise samples are generated based on O3a + O3b dataset https://www.gw-openscience.org/data/ . Even if some white noise is added, we were able to relayably map test samples to O3 data (but in this case instead of correlation we used L2 distance with torch.cdist). The next step could be denoising of all nonstationary noise samples and feeding them to a regular model: since the residual noise is weak, the confidence of the model could be substantially improved. Running this denoising could be possible within 1-2 days, but since we discovered integer frequency shifts in the last day of the competition, we did not have time to perform this procedure.","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os,gc,random\nimport pickle\nfrom tqdm.auto import tqdm\nfrom collections import OrderedDict\nimport h5py\nimport torch\nimport matplotlib.pyplot as plt\nfrom numpy.lib.stride_tricks import sliding_window_view","metadata":{"execution":{"iopub.status.busy":"2023-01-06T03:04:12.554072Z","iopub.execute_input":"2023-01-06T03:04:12.554506Z","iopub.status.idle":"2023-01-06T03:04:14.213681Z","shell.execute_reply.started":"2023-01-06T03:04:12.554425Z","shell.execute_reply":"2023-01-06T03:04:14.212773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"PATH = '../input/g2net-detecting-continuous-gravitational-waves/test/'\nOUT = 'leak.pickle'\nTARGET_OUT = 'leak.csv'\ndf = pd.read_csv('../input/g2net2022-stationery-list/test_stationery.csv')","metadata":{"execution":{"iopub.status.busy":"2023-01-06T04:02:24.823263Z","iopub.execute_input":"2023-01-06T04:02:24.823951Z","iopub.status.idle":"2023-01-06T04:02:24.84276Z","shell.execute_reply.started":"2023-01-06T04:02:24.823915Z","shell.execute_reply":"2023-01-06T04:02:24.841676Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def extract_data_from_hdf5(path):\n    data = {}\n    with h5py.File(path, \"r\") as f:\n        ID_key = list(f.keys())[0]\n        # Retrieve the frequency data\n        data['freq'] = np.array(f[ID_key]['frequency_Hz'])\n        # Retrieve the Livingston decector data\n        data['L1_SFTs_amplitudes'] = np.array(f[ID_key]['L1']['SFTs'])\n        data['L1_ts'] = np.array(f[ID_key]['L1']['timestamps_GPS'])\n        # Retrieve the Hanford decector data\n        data['H1_SFTs_amplitudes'] = np.array(f[ID_key]['H1']['SFTs'])\n        data['H1_ts'] = np.array(f[ID_key]['H1']['timestamps_GPS'])\n    return data\n\ndef get_correlation(x1,x2):\n    correlation = (torch.einsum('ki,kj->ij', x1, x2)/x1.shape[0] - \\\n                   torch.einsum('i,j->ij', x1.mean(0), x2.mean(0)))/ \\\n                   torch.einsum('i,j->ij', x1.std(0), x2.std(0))\n    return correlation","metadata":{"execution":{"iopub.status.busy":"2023-01-06T03:05:37.286023Z","iopub.execute_input":"2023-01-06T03:05:37.286426Z","iopub.status.idle":"2023-01-06T03:05:37.296136Z","shell.execute_reply.started":"2023-01-06T03:05:37.286394Z","shell.execute_reply":"2023-01-06T03:05:37.294785Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clean_data,denoised_target = {},{}\nTH = 0.95\nfor index, row in tqdm(df.iterrows(), total=len(df)):\n    idx, freq_origin = row[['id','freq']]\n    if idx != 'a2756b16d': continue # remove this line to iterate over all files\n    data_src = extract_data_from_hdf5(os.path.join(PATH, idx+'.hdf5')) \n    buf = {\n        'H1': torch.zeros(data_src['H1_SFTs_amplitudes'].shape),\n        'L1': torch.zeros(data_src['L1_SFTs_amplitudes'].shape),\n        'H1_mask': torch.zeros(data_src['H1_SFTs_amplitudes'].shape, dtype=torch.uint8),\n        'L1_mask': torch.zeros(data_src['L1_SFTs_amplitudes'].shape, dtype=torch.uint8),\n        'denoise_pair_H': [],\n        'denoise_pair_L': [],\n        'denoise_pos_imp': [],\n        'denoise_neg_imp': [],\n        'coverage_H1': 0,\n        'coverage_L1': 0,\n        'empty' : True\n    }\n    for offset in np.arange(-3, 4):\n        freq = freq_origin\n        freq += offset\n        neigbors = df.loc[np.abs(df.freq.values - freq) < 0.2, ['id','freq']].set_index('id').to_dict()['freq']\n        if len(neigbors) == 0: continue\n        \n        for key in neigbors.keys():\n            if key == idx:\n                continue\n            dfreq = round(1800*(freq - neigbors[key]))\n            if 360 - abs(dfreq) < 40: continue\n            f1s,f2s = max(-dfreq,0), min(len(data_src['freq']) - dfreq,len(data_src['freq']))\n            f1t,f2t = max(dfreq,0), min(len(data_src['freq']) + dfreq,len(data_src['freq']))\n            data_tgt = extract_data_from_hdf5(os.path.join(PATH, key+'.hdf5'))\n            \n            src_H = torch.from_numpy(np.abs(data_src['H1_SFTs_amplitudes']*1e22))\n            tgt_H = torch.from_numpy(np.abs(data_tgt['H1_SFTs_amplitudes']*1e22))\n            correlation_H = get_correlation(src_H[f1s:f2s],tgt_H[f1t:f2t])\n            src_L = torch.from_numpy(np.abs(data_src['L1_SFTs_amplitudes']*1e22))\n            tgt_L = torch.from_numpy(np.abs(data_tgt['L1_SFTs_amplitudes']*1e22))\n            correlation_L = get_correlation(src_L[f1s:f2s],tgt_L[f1t:f2t])\n            \n            if correlation_H.max() > TH:\n                values,indices = correlation_H.max(-1)\n                dif_abs = (src_H[f1s:f2s, values > TH] - tgt_H[f1t:f2t,indices[values > TH]]).abs()\n                min_val = torch.min(buf['H1'][f1s:f2s,values > TH], dif_abs)\n                buf['H1'][f1s:f2s,values > TH] = torch.where(buf['H1_mask'][f1s:f2s,values > TH] == 0,\n                        dif_abs, min_val)\n                buf['H1_mask'][f1s:f2s,values > TH] += 1\n                buf['empty'] = False\n                buf['denoise_pair_H'].append(key)\n\n                sample1 = src_H[f1s:f2s, values > TH].mean(1)\n                sample2 = tgt_H[f1t:f2t, indices[values > TH]].mean(1)\n                diff_raw = sample1 - sample2\n                diff_raw_roll10 = sliding_window_view(diff_raw, 10, axis=0).mean(axis=-1)\n                neg_peak_raw = np.abs(diff_raw_roll10[diff_raw_roll10 < -1e-4]).sum()\n                pos_peak_raw = np.abs(diff_raw_roll10[diff_raw_roll10 > 1e-4]).sum()\n                buf['denoise_pos_imp'] = pos_peak_raw\n                buf['denoise_neg_imp'] = neg_peak_raw\n                negative_peak = neg_peak_raw > 1.05e-3\n                positive_peak = pos_peak_raw > 1.05e-3\n                if negative_peak:\n                    if not key in denoised_target.keys():\n                        denoised_target[key] = [1]\n                    else:\n                        denoised_target[key].append(1)\n                else:\n                    if not key in denoised_target.keys():\n                        denoised_target[key] = [0]\n                    else:\n                        denoised_target[key].append(0)\n                if positive_peak:\n                    if not idx in denoised_target.keys():\n                        denoised_target[idx] = [1]\n                    else:\n                        denoised_target[idx].append(1)\n                else:\n                    if not idx in denoised_target.keys():\n                        denoised_target[idx] = [0]\n                    else:\n                        denoised_target[idx].append(0)\n\n            if correlation_L.max() > TH:\n                values,indices = correlation_L.max(-1)\n                dif_abs = (src_L[f1s:f2s, values > TH] - tgt_L[f1t:f2t,indices[values > TH]]).abs()\n                min_val = torch.min(buf['L1'][f1s:f2s,values > TH], dif_abs)\n                buf['L1'][f1s:f2s,values > TH] = torch.where(buf['L1_mask'][f1s:f2s,values > TH] == 0,\n                        dif_abs, min_val)\n                buf['L1_mask'][f1s:f2s,values > TH] += 1\n                buf['empty'] = False\n                buf['denoise_pair_L'].append(key)\n\n    if not buf['empty']:\n        del buf['empty']\n        buf['coverage_H1'] = (buf['H1_mask'].amax(1) > 0).sum() / 360\n        buf['coverage_L1'] = (buf['L1_mask'].amax(1) > 0).sum() / 360\n        clean_data[idx] = buf\n\ndenoised_target = {k: max(v) for k, v in denoised_target.items()}\nfor k, v in denoised_target.items():\n    if k in clean_data.keys():\n        clean_data[k]['label'] = v\n\ndenoised_target = pd.DataFrame.from_dict(denoised_target, orient='index').reset_index()\ndenoised_target.columns = ['id', 'target']\ndenoised_target.to_csv(TARGET_OUT, index=False)\n\nwith open(OUT, 'wb') as f:\n    pickle.dump(clean_data, f, protocol=pickle.HIGHEST_PROTOCOL)","metadata":{"execution":{"iopub.status.busy":"2023-01-06T04:05:31.502914Z","iopub.execute_input":"2023-01-06T04:05:31.50332Z","iopub.status.idle":"2023-01-06T04:06:09.352885Z","shell.execute_reply.started":"2023-01-06T04:05:31.503289Z","shell.execute_reply":"2023-01-06T04:06:09.351877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#the signal is present in both files\ndenoised_target.head()","metadata":{"execution":{"iopub.status.busy":"2023-01-06T04:06:12.507989Z","iopub.execute_input":"2023-01-06T04:06:12.508772Z","iopub.status.idle":"2023-01-06T04:06:12.519766Z","shell.execute_reply.started":"2023-01-06T04:06:12.508735Z","shell.execute_reply":"2023-01-06T04:06:12.51887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(16, 16))\nplt.imshow(clean_data[list(clean_data.keys())[0]]['H1'][:,:512])","metadata":{"execution":{"iopub.status.busy":"2023-01-06T04:06:09.354922Z","iopub.execute_input":"2023-01-06T04:06:09.355452Z","iopub.status.idle":"2023-01-06T04:06:09.825827Z","shell.execute_reply.started":"2023-01-06T04:06:09.355409Z","shell.execute_reply":"2023-01-06T04:06:09.824678Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}