{"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":"# G2Net: Winning with External Data\n\n\n## TL;DR\n\nHere I'm giving away a strategy that will likely give a huge advantage if implemented correctly. In short, it might be possible to reverse engineer the injected CW signal with Ligo&Virgo data. \n\n1. Get raw detector data at https://www.gw-openscience.org/archive/O3a_4KHZ_R1/\n2. Generate SFTs using instructions at https://youtu.be/A9iWRcmG0Rs?t=2199 (you'll likely need to google for more details on this step). Use GPS timestamps from `test` set to produce SFTs. \n3. Read generated SFTs using `pyfstat.utils.sft.get_sft_as_arrays` \n4. Match generated SFTs with _real samples_ from `test` set. \n5. Find the difference between SFTs and real samples. `difference > EPS -> label==1; difference <= EPS -> label==0.`\n\nThis is not the solution organizers had in mind, but it's completely legitimate according to the rules: https://www.kaggle.com/competitions/g2net-detecting-continuous-gravitational-waves/rules (see info on external data usage)\n\n## Intuition\n\n### 1. Why use https://www.gw-openscience.org/archive/O3a_4KHZ_R1/ dataset? \n\nSamples come from the approximate period between 1238170136 and 1248567118 (GPS time). You can check at https://www.andrews.edu/~tzs/timeconv/timeconvert.php that it corresponds to [Apr 01, 2019 - Jul 31, 2019] which perfectly falls into `O3a_4KHZ_R1` dataset release. \n\n**Please, upvote if you think this idea is crazy enough to give it a try!**","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"import glob\nimport numpy as np\nimport pandas as pd\nfrom tqdm import tqdm\nfrom pathlib import Path\nimport h5py\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib import colors\nfrom scipy import stats\nimport os\n\nplt.rcParams[\"font.family\"] = \"serif\"\nplt.rcParams[\"font.size\"] = 20\n\n\n\n# Utility to read hdf5 file\ndef read_data(file):\n    file = Path(file)\n    with h5py.File(file, \"r\") as f:\n        filename = file.stem\n        f = f[filename]\n        h1 = f[\"H1\"]\n        l1 = f[\"L1\"]\n        freq_hz = list(f[\"frequency_Hz\"])\n        \n        h1_stft = h1[\"SFTs\"][()]\n        h1_timestamp = h1[\"timestamps_GPS\"][()]\n        # H2 data\n        l1_stft = l1[\"SFTs\"][()]\n        l1_timestamp = l1[\"timestamps_GPS\"][()]\n        \n        return [h1_stft, h1_timestamp],            [l1_stft, l1_timestamp], np.array(freq_hz)\n        \n\n(h1_sfts, h1_ts), (l1_sfts, l1_ts), freq = read_data('/kaggle/input/g2net-detecting-continuous-gravitational-waves/test/00054c878.hdf5')\nprint('~Test dataset period:', min(h1_ts), max(h1_ts))","metadata":{"execution":{"iopub.status.busy":"2022-11-30T11:52:23.055599Z","iopub.execute_input":"2022-11-30T11:52:23.056151Z","iopub.status.idle":"2022-11-30T11:52:24.913688Z","shell.execute_reply.started":"2022-11-30T11:52:23.056039Z","shell.execute_reply":"2022-11-30T11:52:24.912663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### 2. Real vs generated samples in the test set\n\nFirst note that the test dataset contains both real and generated samples. It is actually quite easy to differentiate them:\n* Real samples (~20%) have non-stationary noise and this noise shows consistent pattern accross multiple bands. We will only focus on these samples in this notebook. \n* 80% of test samples use stationary noise and are therefore generated. Normal comptuer vision ML models will be used to tackle this subset. \n\nLet's look at some examples:","metadata":{}},{"cell_type":"code","source":"df_stat_test = []\nnp.random.seed(123)\nfor fname in tqdm(np.random.choice(glob.glob('/kaggle/input/g2net-detecting-continuous-gravitational-waves/test/*'), size=100)):\n    (h1_sfts, h1_ts), (l1_sfts, l1_ts), freq = read_data(fname)\n    stds = np.array([np.std(v) for v in np.array_split(np.absolute(h1_sfts).astype(np.float64), 10,  axis=1)])    \n    \n    (h1_sfts, h1_ts), (l1_sfts, l1_ts), freq = read_data(fname)\n    df_stat_test.append({\n        'fname': fname,\n        'minstd': min(stds),\n        'maxstd': max(stds),\n        'stddiff': max(stds) - min(stds)\n    })\n    \nprint('Test set amplitude STD stats')\ndf_stat_test = pd.DataFrame(df_stat_test)\ndf_stat_test","metadata":{"execution":{"iopub.status.busy":"2022-11-30T11:52:24.915777Z","iopub.execute_input":"2022-11-30T11:52:24.916167Z","iopub.status.idle":"2022-11-30T11:53:08.328697Z","shell.execute_reply.started":"2022-11-30T11:52:24.916124Z","shell.execute_reply":"2022-11-30T11:53:08.327911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15, 5))\nax = df_stat_test.stddiff.plot.hist(bins=100, title='Test data: distribution of variation of noise across time buckets')\nax.axvline(1e-24)","metadata":{"execution":{"iopub.status.busy":"2022-11-30T11:53:08.330111Z","iopub.execute_input":"2022-11-30T11:53:08.330673Z","iopub.status.idle":"2022-11-30T11:53:08.81489Z","shell.execute_reply.started":"2022-11-30T11:53:08.330641Z","shell.execute_reply":"2022-11-30T11:53:08.813451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In the distribution above there is a visibile distinction between a set of samples with <1e-24 and >1e-24 variation of noise STD. \n* `< 1e-24` corresponds to the generated stationary noise\n* `> 1e-24` is nonstationary noise","metadata":{}},{"cell_type":"code","source":"# There are ~80% of samples with stationary noise:\nlen(df_stat_test[df_stat_test.stddiff < 1e-24])","metadata":{"execution":{"iopub.status.busy":"2022-11-30T11:53:08.817756Z","iopub.execute_input":"2022-11-30T11:53:08.818425Z","iopub.status.idle":"2022-11-30T11:53:08.827727Z","shell.execute_reply.started":"2022-11-30T11:53:08.818385Z","shell.execute_reply":"2022-11-30T11:53:08.82642Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# There are only ~20% of samples with non-stationary noise:\nlen(df_stat_test[df_stat_test.stddiff >= 1e-24])","metadata":{"execution":{"iopub.status.busy":"2022-11-30T11:53:08.830101Z","iopub.execute_input":"2022-11-30T11:53:08.830622Z","iopub.status.idle":"2022-11-30T11:53:08.843111Z","shell.execute_reply.started":"2022-11-30T11:53:08.830557Z","shell.execute_reply":"2022-11-30T11:53:08.841749Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_amplitude_spectrogram(timestamps, frequency, fourier_data):\n    \n    ax = plt.subplot(1, 2, 1)\n    ax.set(xlabel=\"SFT index\", ylabel=\"Frequency [Hz]\")\n    time_in_days = (timestamps - timestamps[0]) / 3600 / 24\n    ax.set_title(\"SFT amplitude\")\n    c = ax.pcolorfast(\n        time_in_days, frequency, np.absolute(fourier_data)[:-1, :-1], norm=colors.Normalize()\n    )\n    \n    ax = plt.subplot(1, 2, 2)\n    \n    noise_levels = np.std(np.absolute(fourier_data).astype(np.float64), axis=0)\n    \n    ax.plot(noise_levels)\n    ax.set_title('SFT Amplitude STDDEV')\n    \n","metadata":{"execution":{"iopub.status.busy":"2022-11-30T11:53:08.845086Z","iopub.execute_input":"2022-11-30T11:53:08.845783Z","iopub.status.idle":"2022-11-30T11:53:08.861728Z","shell.execute_reply.started":"2022-11-30T11:53:08.845736Z","shell.execute_reply":"2022-11-30T11:53:08.86034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i, (_, r) in enumerate(df_stat_test.query('stddiff < 1e-24').head(3).iterrows()):\n    (h1_sfts, h1_ts), (l1_sfts, l1_ts), freq = read_data(r.fname)\n    plt.figure(figsize=(15, 7))\n    plt.suptitle(f'Generated test sample: {os.path.basename(r.fname)}')\n    plot_amplitude_spectrogram(\n        h1_ts, np.array(freq), h1_sfts\n    )","metadata":{"execution":{"iopub.status.busy":"2022-11-30T11:53:08.863507Z","iopub.execute_input":"2022-11-30T11:53:08.864476Z","iopub.status.idle":"2022-11-30T11:53:11.521573Z","shell.execute_reply.started":"2022-11-30T11:53:08.864427Z","shell.execute_reply":"2022-11-30T11:53:11.519736Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i, (_, r) in enumerate(df_stat_test.query('stddiff > 1e-24').head(2).iterrows()):\n    (h1_sfts, h1_ts), (l1_sfts, l1_ts), freq = read_data(r.fname)\n    plt.figure(figsize=(15, 7))\n    plt.suptitle(f'Real test sample: {os.path.basename(r.fname)}')\n    plot_amplitude_spectrogram(\n        h1_ts, np.array(freq), h1_sfts\n    )","metadata":{"execution":{"iopub.status.busy":"2022-11-30T11:53:11.523482Z","iopub.execute_input":"2022-11-30T11:53:11.524353Z","iopub.status.idle":"2022-11-30T11:53:13.150164Z","shell.execute_reply.started":"2022-11-30T11:53:11.5243Z","shell.execute_reply":"2022-11-30T11:53:13.149188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Note the difference between std plots of generated and real samples. \n* Generated samples always use stationary noise \n* Real samples use non-stationary noise. Additionally you can see that the noise levels drop down over time. This corresponds to incremental improvements of Ligo&Virgo instruments over time. ","metadata":{}},{"cell_type":"markdown","source":"### 3. How do I know that all CW signals (label=1) were artificially injected into \"real\" test samples?\n\nSo what is the problem here? For now consider only real test samples. For `label=1` there might in theory be 2 possibilities:\n1. Some real CW signal was found by physicists and correctly labelled for the competition dataset. \n2. The signal was artificially injected. \n\nI want to argue here that it's always the latter case for all the 20% of real samples here. Correct me if I'm wrong, but I couldn't find any evidence that Ligo&Virgo collaboration ever found reliable CW signal in `O3a_4KHZ_R1` data. \n\nSo if it's artificially injected, we can find a difference between ground truth SFTs and competition data, right?","metadata":{}},{"cell_type":"markdown","source":"### 4. Why not just use traditional ML models for all test samples?\n\nFirst of all, clearly, having 20% of test data accurately predicted is a very nice bonus, but I want to make an even stronger point here. \n\nFor some reason the trainig set contains only **generated samples** as can be seen below. \n\nThis means that: \n1. Any model trained on only `train` data will experience significant distribution shift when predicting the 20% of \"real\" test samples and will likely get them wrong. \n2. Any model that uses train data augmentations and/or `pyfstat`-generated synthetic data with likely will perform **worse** than the model trained on samples acquired directly from https://www.gw-openscience.org/archive/O3a_4KHZ_R1/\n\nSo anyway using `O3a_4KHZ_R1` dataset sound like a good idea. ","metadata":{}},{"cell_type":"code","source":"df_train = pd.read_csv('/kaggle/input/g2net-detecting-continuous-gravitational-waves/train_labels.csv')\ndf_train","metadata":{"execution":{"iopub.status.busy":"2022-11-30T11:53:13.151508Z","iopub.execute_input":"2022-11-30T11:53:13.152044Z","iopub.status.idle":"2022-11-30T11:53:13.171747Z","shell.execute_reply.started":"2022-11-30T11:53:13.152011Z","shell.execute_reply":"2022-11-30T11:53:13.170962Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_stat = []\nnp.random.seed(123)\nfor _, r in tqdm(df_train.query('target >= 0').sample(100, random_state=42).iterrows(), total=100):\n    fname = f'/kaggle/input/g2net-detecting-continuous-gravitational-waves/train/{r.id}.hdf5'\n    (h1_sfts, h1_ts), (l1_sfts, l1_ts), freq = read_data(fname)\n    stds = np.array([np.std(v) for v in np.array_split(np.absolute(h1_sfts).astype(np.float64), 10,  axis=1)])    \n    \n    (h1_sfts, h1_ts), (l1_sfts, l1_ts), freq = read_data(fname)\n    df_stat.append({\n        'fname': fname,\n        'minstd': min(stds),\n        'maxstd': max(stds),\n        'stddiff': max(stds) - min(stds)\n    })\n    \nprint('Training set amplitude STD stats')\ndf_stat = pd.DataFrame(df_stat)\ndf_stat","metadata":{"execution":{"iopub.status.busy":"2022-11-30T11:53:13.174527Z","iopub.execute_input":"2022-11-30T11:53:13.175483Z","iopub.status.idle":"2022-11-30T11:53:53.999241Z","shell.execute_reply.started":"2022-11-30T11:53:13.175441Z","shell.execute_reply":"2022-11-30T11:53:53.998202Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15, 5))\nax = df_stat.stddiff.plot.hist(bins=100, title='Train data: distribution of variation of noise across time buckets')\nax.axvline(1e-24)","metadata":{"execution":{"iopub.status.busy":"2022-11-30T11:53:54.000549Z","iopub.execute_input":"2022-11-30T11:53:54.001712Z","iopub.status.idle":"2022-11-30T11:53:54.424938Z","shell.execute_reply.started":"2022-11-30T11:53:54.001679Z","shell.execute_reply":"2022-11-30T11:53:54.423744Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('No \"real\" samples in the training set: ', len(df_stat.query('stddiff > 1e-24')))","metadata":{"execution":{"iopub.status.busy":"2022-11-30T11:53:54.426389Z","iopub.execute_input":"2022-11-30T11:53:54.428226Z","iopub.status.idle":"2022-11-30T11:53:54.438791Z","shell.execute_reply.started":"2022-11-30T11:53:54.428178Z","shell.execute_reply":"2022-11-30T11:53:54.437388Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i, (_, r) in enumerate(df_stat.head(3).iterrows()):\n    (h1_sfts, h1_ts), (l1_sfts, l1_ts), freq = read_data(fname)\n    plt.figure(figsize=(15, 7))\n    plt.suptitle(f'Generated train sample: {os.path.basename(r.fname)}')\n    plot_amplitude_spectrogram(\n        h1_ts, np.array(freq), h1_sfts\n    )","metadata":{"execution":{"iopub.status.busy":"2022-11-30T11:53:54.440386Z","iopub.execute_input":"2022-11-30T11:53:54.441064Z","iopub.status.idle":"2022-11-30T11:53:57.094164Z","shell.execute_reply.started":"2022-11-30T11:53:54.441028Z","shell.execute_reply":"2022-11-30T11:53:57.092779Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Detect all real test samples\n\nThe code detects all non-stationary noise in test dataset. \nUse `test.csv` in your notebooks to differentiate samples with real vs generated data. ","metadata":{}},{"cell_type":"code","source":"%%time\n\nimport multiprocessing as mp\n\ndef gen_stddiff(file):\n    file = Path(file)\n    with h5py.File(file, \"r\") as f:\n        filename = file.stem\n        f = f[filename]\n        h1 = f[\"H1\"]        \n        h1_sft = h1[\"SFTs\"][()]        \n        stds = np.array([np.std(v) for v in np.array_split(np.absolute(h1_sft).astype(np.float64), 10,  axis=1)])            \n        return {\n            'id': os.path.basename(file)[:-5],\n            'noise_diff': max(stds) - min(stds),\n        }\n\nwith mp.Pool() as p:\n    df_test = pd.DataFrame(p.map(gen_stddiff, glob.glob('/kaggle/input/g2net-detecting-continuous-gravitational-waves/test/*')))\n    \ndf_test    ","metadata":{"execution":{"iopub.status.busy":"2022-11-30T11:53:57.095723Z","iopub.execute_input":"2022-11-30T11:53:57.096119Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(15, 7))\nax = df_test.noise_diff.plot.hist(bins=1000, title='Noise variation, all test data')\nax.axvline(1e-24)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test['is_generated_noise'] = df_test.noise_diff < 1e-24\nprint('pct of generated samples:', df_test['is_generated_noise'].mean())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test[['id', 'is_generated_noise']].to_csv('test.csv', index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Caveats\n\nThere are certain caveats that I want to mention here. \n1. It's possible, that there is some additional noise injected into \"real\" samples. In this case matching them with ground truth SFTs becomes more challenging. \n2. It's not easy to reproduce SFT files (step 2). Scarce documentation and not too many examples of how to do it. But I'm pretty sure organizers took this path, so it's definitely reproducible. \n\n\n**Please, upvote if you think this idea is crazy enough to give it a try!**","metadata":{}}]}