{"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":"code","source":"import h5py\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-01-05T07:21:44.889463Z","iopub.execute_input":"2023-01-05T07:21:44.890802Z","iopub.status.idle":"2023-01-05T07:21:44.897265Z","shell.execute_reply.started":"2023-01-05T07:21:44.890745Z","shell.execute_reply":"2023-01-05T07:21:44.895655Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In this notebook shows how my team generated noise online.","metadata":{}},{"cell_type":"markdown","source":"# Load origin sample","metadata":{}},{"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'])*1e22\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'])*1e22\n        data['H1_ts'] = np.array(f[ID_key]['H1']['timestamps_GPS'])\n    return data\n\n\ndef to_spectrogram(sfts):\n    return sfts.real**2 + sfts.imag**2","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:23:31.393432Z","iopub.execute_input":"2023-01-05T07:23:31.393803Z","iopub.status.idle":"2023-01-05T07:23:31.403401Z","shell.execute_reply.started":"2023-01-05T07:23:31.393773Z","shell.execute_reply":"2023-01-05T07:23:31.402062Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_data = extract_data_from_hdf5('../input/g2net-detecting-continuous-gravitational-waves/test/00054c878.hdf5')\nspec_h1 = to_spectrogram(sample_data['H1_SFTs_amplitudes'])","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:26:16.629453Z","iopub.execute_input":"2023-01-05T07:26:16.629909Z","iopub.status.idle":"2023-01-05T07:26:17.173553Z","shell.execute_reply.started":"2023-01-05T07:26:16.629872Z","shell.execute_reply":"2023-01-05T07:26:17.17222Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.imshow(spec_h1[:, :512])","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:26:17.175515Z","iopub.execute_input":"2023-01-05T07:26:17.175873Z","iopub.status.idle":"2023-01-05T07:26:17.412354Z","shell.execute_reply.started":"2023-01-05T07:26:17.175841Z","shell.execute_reply":"2023-01-05T07:26:17.411261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"spec_h1.mean(), spec_h1.std() # mean, std","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:30:52.048531Z","iopub.execute_input":"2023-01-05T07:30:52.049324Z","iopub.status.idle":"2023-01-05T07:30:52.06223Z","shell.execute_reply.started":"2023-01-05T07:30:52.049288Z","shell.execute_reply":"2023-01-05T07:30:52.06059Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Get statistics from the sample\nCalculate mean amptitude for each data point.","metadata":{}},{"cell_type":"code","source":"mean_amp_time = spec_h1.mean(0)","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:26:45.954006Z","iopub.execute_input":"2023-01-05T07:26:45.954461Z","iopub.status.idle":"2023-01-05T07:26:45.960351Z","shell.execute_reply.started":"2023-01-05T07:26:45.954422Z","shell.execute_reply":"2023-01-05T07:26:45.958945Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(20, 3))\nplt.plot(mean_amp_time)","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:28:39.42933Z","iopub.execute_input":"2023-01-05T07:28:39.429717Z","iopub.status.idle":"2023-01-05T07:28:39.688843Z","shell.execute_reply.started":"2023-01-05T07:28:39.429687Z","shell.execute_reply":"2023-01-05T07:28:39.687742Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Reconstruct from statistics","metadata":{}},{"cell_type":"code","source":"def reconstruct_from_stat_complex(mean_arr):\n    real = np.random.normal(size=(360, len(mean_arr)))\n    imag = np.random.normal(size=(360, len(mean_arr)))\n    for t, mean in enumerate(mean_arr):\n        factor = mean / (real[:, t]**2+imag[:, t]**2).mean()\n        real[:, t] *= np.sqrt(factor)\n        imag[:, t] *= np.sqrt(factor)\n    return (real + imag * 1j).astype(np.complex64)","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:33:22.777499Z","iopub.execute_input":"2023-01-05T07:33:22.77795Z","iopub.status.idle":"2023-01-05T07:33:22.786325Z","shell.execute_reply.started":"2023-01-05T07:33:22.777906Z","shell.execute_reply":"2023-01-05T07:33:22.784723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"new_sfts_h1 = reconstruct_from_stat_complex(mean_amp_time)\nnew_spec_h1 = to_spectrogram(new_sfts_h1)","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:33:23.698448Z","iopub.execute_input":"2023-01-05T07:33:23.698819Z","iopub.status.idle":"2023-01-05T07:33:23.961964Z","shell.execute_reply.started":"2023-01-05T07:33:23.698789Z","shell.execute_reply":"2023-01-05T07:33:23.960652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.imshow(new_spec_h1[:, :512])","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:33:24.027824Z","iopub.execute_input":"2023-01-05T07:33:24.028489Z","iopub.status.idle":"2023-01-05T07:33:24.330551Z","shell.execute_reply.started":"2023-01-05T07:33:24.028438Z","shell.execute_reply":"2023-01-05T07:33:24.329401Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"new_spec_h1.mean(), new_spec_h1.std() # mean, std","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:33:25.122286Z","iopub.execute_input":"2023-01-05T07:33:25.122705Z","iopub.status.idle":"2023-01-05T07:33:25.133671Z","shell.execute_reply.started":"2023-01-05T07:33:25.122669Z","shell.execute_reply":"2023-01-05T07:33:25.132571Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Reconstructed spectrogram has almost identical mean/std to the original spectrogram.","metadata":{}},{"cell_type":"markdown","source":"# Some other images\n## Stationary noise","metadata":{}},{"cell_type":"code","source":"def resize_image(img, img_size=360):\n    f, t = img.shape\n    resize_f = t // img_size\n    t2 = img_size * resize_f\n    img = np.mean(\n        img[:, :t2].reshape(f, t2//resize_f, resize_f), axis=2)\n    return img\n    \n\ndef generate_noise(path):\n    sample_data = extract_data_from_hdf5(path)\n    spec_h1 = to_spectrogram(sample_data['H1_SFTs_amplitudes'])\n    print(f'origin statistics: {spec_h1.mean():.5f} {spec_h1.std():.5f}')\n    new_sfts_h1 = reconstruct_from_stat_complex(spec_h1.mean(0))\n    new_spec_h1 = to_spectrogram(new_sfts_h1)\n    print(f'generated statistics: {new_spec_h1.mean():.5f} {new_spec_h1.std():.5f}')\n    plt.figure(figsize=(12, 6))\n    plt.subplot(1, 2, 1)\n    plt.imshow(resize_image(spec_h1))\n    plt.title('origin')\n    plt.subplot(1, 2, 2)\n    plt.imshow(resize_image(new_spec_h1))\n    plt.title('generated')\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:42:40.551367Z","iopub.execute_input":"2023-01-05T07:42:40.552454Z","iopub.status.idle":"2023-01-05T07:42:40.563053Z","shell.execute_reply.started":"2023-01-05T07:42:40.5524Z","shell.execute_reply":"2023-01-05T07:42:40.561825Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"generate_noise('../input/g2net-detecting-continuous-gravitational-waves/test/00529384a.hdf5')","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:42:41.64418Z","iopub.execute_input":"2023-01-05T07:42:41.644593Z","iopub.status.idle":"2023-01-05T07:42:42.754084Z","shell.execute_reply.started":"2023-01-05T07:42:41.64456Z","shell.execute_reply":"2023-01-05T07:42:42.753262Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"generate_noise('../input/g2net-detecting-continuous-gravitational-waves/test/0074f67f0.hdf5')","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:43:21.269896Z","iopub.execute_input":"2023-01-05T07:43:21.270322Z","iopub.status.idle":"2023-01-05T07:43:22.448097Z","shell.execute_reply.started":"2023-01-05T07:43:21.270285Z","shell.execute_reply":"2023-01-05T07:43:22.447055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Nonstationary noise","metadata":{}},{"cell_type":"code","source":"generate_noise('../input/g2net-detecting-continuous-gravitational-waves/test/55c067665.hdf5')","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:44:17.722453Z","iopub.execute_input":"2023-01-05T07:44:17.722861Z","iopub.status.idle":"2023-01-05T07:44:18.962717Z","shell.execute_reply.started":"2023-01-05T07:44:17.722827Z","shell.execute_reply":"2023-01-05T07:44:18.960034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"generate_noise('../input/g2net-detecting-continuous-gravitational-waves/test/55e54564a.hdf5')","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:44:50.427545Z","iopub.execute_input":"2023-01-05T07:44:50.427974Z","iopub.status.idle":"2023-01-05T07:44:51.615655Z","shell.execute_reply.started":"2023-01-05T07:44:50.42794Z","shell.execute_reply":"2023-01-05T07:44:51.614537Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**They look almost idenitcal!**","metadata":{}},{"cell_type":"markdown","source":"# Appendix: artifact injection\nIn some samples, you can find artifacts. There were many types of artifacts: random walk, point, line etc. We used a simple amplitude based criteria to extract artifact from a spectrogram, and copy-and-pasted it to the generated noise.","metadata":{}},{"cell_type":"code","source":"def extract_artifact(spec, n_sigma=8):\n    if np.iscomplexobj(spec):\n        spec = spec.real**2 + spec.imag**2\n    spec_std = spec.std()\n    spec_min = spec.min()\n    amp_map = (spec - spec_min) / spec_std\n    artifact_map = amp_map > n_sigma\n    return artifact_map\n\n\ndef generate_noise2(path, artifact_sigma=8):\n    sample_data = extract_data_from_hdf5(path)\n    spec_h1 = to_spectrogram(sample_data['H1_SFTs_amplitudes'])\n    artifact_map = extract_artifact(spec_h1, 8)\n    print(f'origin statistics: {spec_h1.mean():.5f} {spec_h1.std():.5f}')\n    new_sfts_h1 = reconstruct_from_stat_complex(spec_h1.mean(0).clip(0, 10))\n    new_spec_h1 = to_spectrogram(new_sfts_h1)\n    new_spec_h1[np.where(artifact_map)] = spec_h1[np.where(artifact_map)] # inject\n    print(f'generated statistics: {new_spec_h1.mean():.5f} {new_spec_h1.std():.5f}')\n    plt.figure(figsize=(12, 6))\n    plt.subplot(1, 2, 1)\n    plt.imshow(resize_image(spec_h1).clip(0, 20)) # clip values to see background\n    plt.title('origin')\n    plt.subplot(1, 2, 2)\n    plt.imshow(resize_image(new_spec_h1).clip(0, 20))\n    plt.title('generated')\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:56:33.400855Z","iopub.execute_input":"2023-01-05T07:56:33.401539Z","iopub.status.idle":"2023-01-05T07:56:33.412909Z","shell.execute_reply.started":"2023-01-05T07:56:33.401489Z","shell.execute_reply":"2023-01-05T07:56:33.411619Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"generate_noise2('../input/g2net-detecting-continuous-gravitational-waves/test/575f47724.hdf5', 8)","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:56:34.020153Z","iopub.execute_input":"2023-01-05T07:56:34.020947Z","iopub.status.idle":"2023-01-05T07:56:35.201388Z","shell.execute_reply.started":"2023-01-05T07:56:34.020902Z","shell.execute_reply":"2023-01-05T07:56:35.198407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"generate_noise2('../input/g2net-detecting-continuous-gravitational-waves/test/56b090eaf.hdf5', 6)","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:56:35.203236Z","iopub.execute_input":"2023-01-05T07:56:35.203699Z","iopub.status.idle":"2023-01-05T07:56:36.539355Z","shell.execute_reply.started":"2023-01-05T07:56:35.203667Z","shell.execute_reply":"2023-01-05T07:56:36.538326Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"generate_noise2('../input/g2net-detecting-continuous-gravitational-waves/test/2cab32f98.hdf5', 6)","metadata":{"execution":{"iopub.status.busy":"2023-01-05T07:59:03.725417Z","iopub.execute_input":"2023-01-05T07:59:03.725838Z","iopub.status.idle":"2023-01-05T07:59:04.886786Z","shell.execute_reply.started":"2023-01-05T07:59:03.725806Z","shell.execute_reply":"2023-01-05T07:59:04.885936Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Although, it is not a perfet way to extract artifacts, I believe this method simulates the impact of artifacts good enough.","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}