{"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":"Hello fellow Kagglers,\n\ninspiration: [Generating continuous gravitational-wave signals\n](https://www.kaggle.com/code/rodrigotenorio/generating-continuous-gravitational-wave-signals)\n\nThis notebook demonstrates the generation of 7,000 training noise samples.\n\n[Training Notebook](https://www.kaggle.com/code/markwijkhuizen/g2net-efficientnetv2-s-generated-data-tf)\n\n[Inference Notebook](https://www.kaggle.com/code/markwijkhuizen/g2net-efficientnetv2-s-generated-data-tf-inference)\n\nV2: 7K -> 5K samples","metadata":{}},{"cell_type":"code","source":"# Kaggle notebooks run on Python 3.7, which was dropped by PyFstat a few relases back.\n# Please, use the following command to install PyFstat on a Kaggle notebook.\n# This will install an up-to-date version of PyFstat with Python 3.7 support.\n# Do use the latest version of PyFstat if you use your own Python >= 3.8 installation.\n!pip install -q git+https://github.com/PyFstat/PyFstat@python37","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport sys\nimport cv2\nimport h5py\nimport random\nimport math\nimport glob\nimport time\nimport multiprocessing\n\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport matplotlib as mpl\nimport tensorflow as tf\n\nimport pyfstat\n\nfrom scipy import stats\nfrom tqdm.notebook import tqdm\nfrom multiprocessing import cpu_count\n\n%matplotlib inline","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:32:07.888402Z","iopub.execute_input":"2022-11-07T19:32:07.888773Z","iopub.status.idle":"2022-11-07T19:32:18.003889Z","shell.execute_reply.started":"2022-11-07T19:32:07.888737Z","shell.execute_reply":"2022-11-07T19:32:18.002814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Config","metadata":{}},{"cell_type":"code","source":"# MatplotLib Global Settings\nmpl.rcParams.update(mpl.rcParamsDefault)\nmpl.rcParams['xtick.labelsize'] = 16\nmpl.rcParams['ytick.labelsize'] = 16\nmpl.rcParams['axes.labelsize'] = 18\nmpl.rcParams['axes.titlesize'] = 24","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:32:18.004999Z","iopub.execute_input":"2022-11-07T19:32:18.005586Z","iopub.status.idle":"2022-11-07T19:32:18.014533Z","shell.execute_reply.started":"2022-11-07T19:32:18.005554Z","shell.execute_reply":"2022-11-07T19:32:18.012415Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_labels = pd.read_csv('/kaggle/input/g2net-detecting-continuous-gravitational-waves/train_labels.csv')\n    \ndisplay(train_labels.info())\n\ndisplay(train_labels.head())","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:32:18.017516Z","iopub.execute_input":"2022-11-07T19:32:18.017988Z","iopub.status.idle":"2022-11-07T19:32:18.084668Z","shell.execute_reply.started":"2022-11-07T19:32:18.017949Z","shell.execute_reply":"2022-11-07T19:32:18.08335Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"IS_INTERACTIVE = os.environ['KAGGLE_KERNEL_RUN_TYPE'] == 'Interactive'\n\nTRAIN_DIR = '/kaggle/input/g2net-detecting-continuous-gravitational-waves/train'\nN_SAMPLES = len(train_labels)\n\nTARGET_HEIGHT = 360\nTARGET_HEIGHT_COMPRESSED = 360\nTARGET_HEIGHT_COMPRESS_RATIO = TARGET_HEIGHT // TARGET_HEIGHT_COMPRESSED\n\nTARGET_WIDTH = 4096\nTARGET_WIDTH_COMPRESSED = 256\nTARGET_WIDTH_COMPRESS_RATIO = TARGET_WIDTH // TARGET_WIDTH_COMPRESSED\n\nNUM_SIGNALS = 100 if IS_INTERACTIVE else 5000","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:32:18.086409Z","iopub.execute_input":"2022-11-07T19:32:18.087572Z","iopub.status.idle":"2022-11-07T19:32:18.094498Z","shell.execute_reply.started":"2022-11-07T19:32:18.087532Z","shell.execute_reply":"2022-11-07T19:32:18.093138Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Label Samples","metadata":{}},{"cell_type":"code","source":"N = 8 if IS_INTERACTIVE else 32\ndf = train_labels.loc[train_labels['target'] == 1].sample(N, random_state=42)\nfor row_idx, row in tqdm(df.iterrows(), total=N):\n    train_id = row['id']\n    target = row['target']\n    if target is not 1:\n        continue        \n\n    with h5py.File(f'{TRAIN_DIR}/{train_id}.hdf5', 'r') as file:\n        SFT_H = np.array(file[train_id]['H1']['SFTs'], dtype=np.complex128)\n        SFT_L = np.array(file[train_id]['L1']['SFTs'], dtype=np.complex128)\n            \n    # Normalize Signal\n    signal_norm = {\n        'H1': SFT_H.real ** 2 + SFT_H.imag ** 2,\n        'L1': SFT_L.real ** 2 + SFT_L.imag ** 2,\n    }\n    \n    fig, ax = plt.subplots(1, 2, figsize=(15, 8))\n    \n    fig.suptitle(f\"Train ID: {train_id}, TARGET: {target}, Row IDX: {row_idx}\")\n        \n    for d_idx, (detector, signal) in enumerate(signal_norm.items()):\n        ax[d_idx].set(\n            xlabel=\"Timestamps [GPS]\",\n            ylabel=\"Frequency [Hz]\",\n            title=f\"{detector} - Normalised Power\",\n        )\n        \n        # Save as PNG\n        patch_uint8 = signal[:, :4096].reshape(\n                TARGET_HEIGHT_COMPRESSED, TARGET_HEIGHT_COMPRESS_RATIO, TARGET_WIDTH_COMPRESSED, TARGET_WIDTH_COMPRESS_RATIO\n            ).mean(axis=(1,3))\n        patch_uint8 = (patch_uint8 - patch_uint8.min())\n        patch_uint8 = patch_uint8 * (255 / patch_uint8.max())\n        patch_uint8 = patch_uint8.astype(np.uint8)\n        \n        c = ax[d_idx].pcolormesh(\n            np.arange(TARGET_WIDTH_COMPRESSED),\n            np.arange(TARGET_HEIGHT_COMPRESSED),\n            patch_uint8,\n            cmap=\"inferno_r\",\n            shading=\"nearest\",\n        )\n        \n    plt.colorbar(c, ax=ax[d_idx], label=\"Normalized Power\")\n        \n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:32:18.096456Z","iopub.execute_input":"2022-11-07T19:32:18.09681Z","iopub.status.idle":"2022-11-07T19:32:31.228548Z","shell.execute_reply.started":"2022-11-07T19:32:18.096781Z","shell.execute_reply":"2022-11-07T19:32:31.227332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Gaps Distribution","metadata":{}},{"cell_type":"code","source":"time_delta_recordings = []\n\nfor row_idx, row in tqdm(train_labels.iterrows(), total=N_SAMPLES):\n    train_id = row['id']\n    c = 0\n    with h5py.File(f'{TRAIN_DIR}/{train_id}.hdf5', 'r') as file:\n        timestamps_h = np.array(file[train_id]['H1']['timestamps_GPS'])\n        for td in (timestamps_h[1:] - timestamps_h[:-1]):\n            if td // 1800 > 0:\n                time_delta_recordings.append(td // 1800)\n            \n        timestamps_l = np.array(file[train_id]['L1']['timestamps_GPS'])\n        for td in (timestamps_l[1:] - timestamps_l[:-1]):\n            if td // 1800 > 0:\n                time_delta_recordings.append(td // 1800)\n            \n                \ntime_delta_recordings = np.array(time_delta_recordings, dtype=np.int32)","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:32:31.230009Z","iopub.execute_input":"2022-11-07T19:32:31.230372Z","iopub.status.idle":"2022-11-07T19:32:49.908143Z","shell.execute_reply.started":"2022-11-07T19:32:31.23034Z","shell.execute_reply":"2022-11-07T19:32:49.906521Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display(pd.Series(time_delta_recordings).describe().to_frame('Ratio'))\n\ndisplay(pd.Series(time_delta_recordings).value_counts().to_frame('Count').head())\n\n\nplt.figure(figsize=(15,8))\nplt.title('Time Delta Measurements Distribution')\npd.Series(time_delta_recordings).plot(kind='hist', bins=100)\nplt.xlim(0, plt.xlim()[1])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:32:49.909961Z","iopub.execute_input":"2022-11-07T19:32:49.91042Z","iopub.status.idle":"2022-11-07T19:32:51.002542Z","shell.execute_reply.started":"2022-11-07T19:32:49.910379Z","shell.execute_reply":"2022-11-07T19:32:51.001323Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_timestamp_idxs():\n    while True:\n        return np.cumsum(np.random.choice(time_delta_recordings, TARGET_WIDTH))","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:32:51.003909Z","iopub.execute_input":"2022-11-07T19:32:51.004251Z","iopub.status.idle":"2022-11-07T19:32:51.010037Z","shell.execute_reply.started":"2022-11-07T19:32:51.00422Z","shell.execute_reply":"2022-11-07T19:32:51.008812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Target Directories","metadata":{}},{"cell_type":"code","source":"!rm -rf train_samples\n!mkdir -p train_samples/x\n!ls -l train_samples","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:32:51.015465Z","iopub.execute_input":"2022-11-07T19:32:51.016245Z","iopub.status.idle":"2022-11-07T19:32:54.337802Z","shell.execute_reply.started":"2022-11-07T19:32:51.016198Z","shell.execute_reply":"2022-11-07T19:32:54.336471Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plot Signals","metadata":{}},{"cell_type":"code","source":"def plot_signal(frequency, timestamps, amplitudes, ind=None, snr=None, target=None):\n    fig, ax = plt.subplots(2, 2, figsize=(20, 15))\n    \n    if None not in [ind, snr, target]:\n        fig.suptitle(f\"Signal {ind} - SNR: {snr:.2f}, TARGET: {target}\")\n        \n    for d_ind, detector in enumerate(amplitudes.keys()):\n        ax[d_ind][0].set(xlabel=\"Timestamps [GPS]\",\n                         ylabel=\"Frequency [Hz]\",\n                         title=f\"{detector} - Real part\")\n        ax[d_ind][1].set(xlabel=\"Timestamps [GPS]\",\n                         ylabel=\"Frequency [Hz]\",\n                         title=f\"{detector} - Imaginary part\")\n\n        c0 = ax[d_ind][0].pcolormesh(\n                (timestamps[detector] - timestamps[detector][0]) / 86400,\n                frequency,\n                amplitudes[detector].real,\n                cmap=\"inferno_r\",\n                shading=\"nearest\",\n            )\n        c1 = ax[d_ind][1].pcolormesh(\n            (timestamps[detector] - timestamps[detector][0]) / 86400,\n            frequency,\n            amplitudes[detector].imag,\n            cmap=\"inferno_r\",\n            shading=\"nearest\",\n        )\n\n        fig.colorbar(c0, ax=ax[d_ind][0])\n        fig.colorbar(c1, ax=ax[d_ind][1])\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:32:54.340036Z","iopub.execute_input":"2022-11-07T19:32:54.343596Z","iopub.status.idle":"2022-11-07T19:32:54.35663Z","shell.execute_reply.started":"2022-11-07T19:32:54.343548Z","shell.execute_reply":"2022-11-07T19:32:54.354864Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plot Normalized Signal","metadata":{}},{"cell_type":"code","source":"def plot_normalized_signal(signal_norm, frequency, timestamps, ind=None, snr=None, target=None):\n    fig, ax = plt.subplots(2, 1, figsize=(15, 15))\n    \n    if None not in [ind, snr, target]:\n        fig.suptitle(f\"Signal {ind} - SNR: {snr:.2f}, TARGET: {target}\", size=24)\n        \n    for d_ind, (detector, signal) in enumerate(signal_norm.items()):\n        ax[d_ind].set(\n            xlabel=\"Timestamps [GPS]\",\n            ylabel=\"Frequency [Hz]\",\n            title=f\"{detector} - Normalised Power\",\n        )\n        \n        c = ax[d_ind].pcolormesh(\n            (timestamps[detector] - timestamps[detector][0]) / 86400,\n            frequency,\n            signal,\n            cmap=\"inferno_r\",\n            shading=\"nearest\",\n        )\n        \n        plt.colorbar(c, ax=ax[d_ind], label=\"Normalized Power\")\n        \n        yticks = np.arange(frequency[0], frequency[-1] + 1e-3, 0.02)\n        ax[d_ind].set_yticks(yticks)\n        ax[d_ind].set_yticklabels([f'{t:.2f}' for t in yticks])\n\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:32:54.358492Z","iopub.execute_input":"2022-11-07T19:32:54.359011Z","iopub.status.idle":"2022-11-07T19:32:54.379659Z","shell.execute_reply.started":"2022-11-07T19:32:54.358963Z","shell.execute_reply":"2022-11-07T19:32:54.378381Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Signal Config","metadata":{}},{"cell_type":"code","source":"# https://pyfstat.readthedocs.io/en/latest/pyfstat.html\n\ndef get_signal_config():\n    t_start = 1238166018\n\n    # These parameters describe background noise and data format\n    writer_kwargs = {\n            'sqrtSX': 1e-23, # Single-sided Amplitude Spectral Density of the noise\n            'Tsft': 1800, # Fourier transform time duration\n            \"SFTWindowType\": \"tukey\",  # Window function to compute short Fourier transforms\n            \"SFTWindowBeta\": 0.01,  # Parameter associated to the window function\n            'timestamps': {\n                'H1': t_start + 1800 * get_timestamp_idxs(),\n                'L1': t_start + 1800 * get_timestamp_idxs(),\n            }\n       }\n\n    # This class allows us to sample signal parameters from a specific population.\n    # Implicitly, sky positions are drawn uniformly across the celestial sphere.\n    # PyFstat also implements a convenient set of priors to sample a population\n    # of isotropically oriented neutron stars.\n    signal_params = {\n        'F0': np.random.uniform(50, 500),\n        'Band': 0.3, # Frequency band-width around F0 [Hz]\n    }\n    \n    return writer_kwargs, signal_params","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:38:39.343965Z","iopub.status.idle":"2022-11-07T19:38:39.344994Z","shell.execute_reply.started":"2022-11-07T19:38:39.344425Z","shell.execute_reply":"2022-11-07T19:38:39.344473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Train Data Generation","metadata":{}},{"cell_type":"code","source":"def generate_cw(ind):\n    np.random.seed(ind)\n    for i in range(10):\n        try:\n            writer_kwargs, signal_params = get_signal_config()\n            # Draw signal parameters.\n            # Noise can be drawn by setting `params[\"h0\"] = 0\n            # Noise\n            if True:\n                target = 0\n            else:\n                target = 1\n\n            writer_kwargs['outdir'] = f'PyFstat_example_data_ensemble/Signal_{ind}'\n            writer_kwargs['label'] = 'Signal'\n\n            writer = pyfstat.Writer(**writer_kwargs, **signal_params)\n\n            writer.make_data()\n\n            snr = 0\n            \n            # Data can be read as a numpy array using PyFstat\n            frequency, timestamps, amplitudes = pyfstat.utils.get_sft_as_arrays(\n                writer.sftfilepath\n            )\n\n            # Cast to complex 128\n            for detector, signal in amplitudes.items():\n                amplitudes[detector] = signal.astype(np.complex128)\n\n            # Normalize Signal\n            signal_norm = {\n                'H1': amplitudes['H1'].real ** 2 + amplitudes['H1'].imag ** 2,\n                'L1': amplitudes['L1'].real ** 2 + amplitudes['L1'].imag ** 2,\n            }\n\n            if ind < 10:\n                display(pd.Series(signal_params).to_frame('Signal Params').T)\n                plot_normalized_signal(signal_norm, frequency, timestamps, ind, snr, target)\n\n            # Save Samples\n            x_dir = f'train_samples/x/{ind}'\n            !mkdir $x_dir\n\n            if False and amplitudes['H1'].shape != (TARGET_SIZE, TARGET_SIZE * N_PATCHES_PER_SAMPLE):\n                print(f'ind {ind} incorrect dimension {amplitudes[\"H1\"].shape}, retrying')\n                sys.stdout.flush()\n                continue\n\n            for detector, signal in amplitudes.items():\n                amplitudes[detector] = signal.reshape(-1, TARGET_WIDTH_COMPRESSED, TARGET_WIDTH_COMPRESS_RATIO).mean(axis=2)\n\n            if ind < 10:\n                fig, ax = plt.subplots(1, 2, figsize=(15, 8))\n                fig.suptitle(f\"Train ID: {train_id}, TARGET: {target}, Row IDX: {row_idx}\")\n\n            for d_idx, (detector, signal) in enumerate(signal_norm.items()):\n                # Save as PNG\n                offset_y = np.random.randint(0, signal.shape[0] - TARGET_HEIGHT)\n                patch_uint8 = signal[offset_y:offset_y + TARGET_HEIGHT, :4096]\n                patch_uint8 = patch_uint8[:, :4096].reshape(\n                        TARGET_HEIGHT_COMPRESSED, TARGET_HEIGHT_COMPRESS_RATIO, TARGET_WIDTH_COMPRESSED, TARGET_WIDTH_COMPRESS_RATIO\n                    ).mean(axis=(1,3))\n                patch_uint8 = (patch_uint8 - patch_uint8.min())\n                patch_uint8 = patch_uint8 * (255 / patch_uint8.max())\n                patch_uint8 = patch_uint8.astype(np.uint8)\n                cv2.imwrite(f'train_samples/x/{ind}/{detector}.png', patch_uint8, [cv2.IMWRITE_PNG_COMPRESSION, 1])\n\n                if ind < 10:\n                    ax[d_idx].set(\n                        xlabel=\"X\",\n                        ylabel=\"Frequency [Hz]\",\n                        title=f\"{detector} - Normalised Power\",\n                    )\n                    \n                    c = ax[d_idx].pcolormesh(\n                        np.arange(TARGET_WIDTH_COMPRESSED),\n                        np.arange(TARGET_HEIGHT_COMPRESSED),\n                        patch_uint8,\n                        cmap=\"inferno_r\",\n                        shading=\"nearest\",\n                    )\n                    plt.colorbar(c, ax=ax[d_idx], label=\"Normalized Power\")\n\n            if ind < 10:\n                plt.show()\n\n            if IS_INTERACTIVE or (ind % 100) == 0:\n                print(f'succesfuly created ind {ind}, shape: {amplitudes[\"H1\"].shape}, snr: {int(snr)}')\n                sys.stdout.flush()\n\n            # Remove Temp Files\n            for fp in glob.glob(f'PyFstat_example_data_ensemble/Signal_{ind}/*'):\n                os.remove(fp)\n\n            return ind, target, snr\n        except Exception as e:\n            if IS_INTERACTIVE:\n                print(f'ind {ind} retrying | {e}')\n                sys.stdout.flush()","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:40:13.176395Z","iopub.execute_input":"2022-11-07T19:40:13.177094Z","iopub.status.idle":"2022-11-07T19:40:13.272565Z","shell.execute_reply.started":"2022-11-07T19:40:13.177037Z","shell.execute_reply":"2022-11-07T19:40:13.271183Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Generate Noise Samples in Parallel\nNUMBER_OF_PROCESSES = multiprocessing.cpu_count()\n\npool = multiprocessing.Pool(processes=NUMBER_OF_PROCESSES)\ntry:\n    map_jobs = [(i,) for i in range(0, NUM_SIGNALS)]\n    result = pool.starmap(generate_cw, map_jobs)\nfinally:\n    pool.close()\n    pool.join()","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:40:13.918526Z","iopub.execute_input":"2022-11-07T19:40:13.919112Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"result_sorted = sorted(result, key=lambda t: t[0])\n_, TARGETS, SNRS = np.array(result_sorted).T\nTARGETS = TARGETS.astype(np.int8)\nSNRS = SNRS.astype(np.float32)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Generation Test","metadata":{}},{"cell_type":"code","source":"N_SAMPLES_CREATED = len(glob.glob('train_samples/x/*/*.png'))\nprint(f'N_SAMPLES_CREATED: {N_SAMPLES_CREATED}')","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:32:57.262773Z","iopub.status.idle":"2022-11-07T19:32:57.263205Z","shell.execute_reply.started":"2022-11-07T19:32:57.263001Z","shell.execute_reply":"2022-11-07T19:32:57.263021Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Tar Training Samples\n!tar -c --remove-files -f train_samples.tar -C train_samples x\n\n# Save targets and signal2noise ratios\nnp.save('TARGETS.npy', TARGETS)\nnp.save('SNRS.npy', SNRS)","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:32:57.264511Z","iopub.status.idle":"2022-11-07T19:32:57.264904Z","shell.execute_reply.started":"2022-11-07T19:32:57.264712Z","shell.execute_reply":"2022-11-07T19:32:57.26473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Target Distribution","metadata":{}},{"cell_type":"code","source":"# Sanity check\ndisplay(pd.Series(TARGETS).value_counts(normalize=True).to_frame('Ratio'))","metadata":{"execution":{"iopub.status.busy":"2022-11-07T19:32:57.265997Z","iopub.status.idle":"2022-11-07T19:32:57.266388Z","shell.execute_reply.started":"2022-11-07T19:32:57.266189Z","shell.execute_reply":"2022-11-07T19:32:57.266206Z"},"trusted":true},"execution_count":null,"outputs":[]}]}