{"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":"!pip -qq install git+https://github.com/PyFstat/PyFstat@python37","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:11:31.758549Z","iopub.execute_input":"2022-10-09T16:11:31.759068Z","iopub.status.idle":"2022-10-09T16:12:18.730089Z","shell.execute_reply.started":"2022-10-09T16:11:31.758964Z","shell.execute_reply":"2022-10-09T16:12:18.727448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport h5py # import to read hdf5\nimport numpy as np\nimport pandas as pd\nfrom glob import glob\nfrom pathlib import Path\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\n# To create graitational waves\nimport pyfstat\nfrom pyfstat.utils import get_sft_as_arrays","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-10-09T16:47:31.953807Z","iopub.execute_input":"2022-10-09T16:47:31.954427Z","iopub.status.idle":"2022-10-09T16:47:31.96401Z","shell.execute_reply.started":"2022-10-09T16:47:31.954382Z","shell.execute_reply":"2022-10-09T16:47:31.962948Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ROOT_DIR = '../input/g2net-detecting-continuous-gravitational-waves'\nos.path.isdir(ROOT_DIR)","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:12:22.310564Z","iopub.execute_input":"2022-10-09T16:12:22.311314Z","iopub.status.idle":"2022-10-09T16:12:22.322565Z","shell.execute_reply.started":"2022-10-09T16:12:22.311274Z","shell.execute_reply":"2022-10-09T16:12:22.320826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.read_csv(f\"{ROOT_DIR}/train_labels.csv\")\ndf.head()","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:12:22.325346Z","iopub.execute_input":"2022-10-09T16:12:22.325726Z","iopub.status.idle":"2022-10-09T16:12:22.399457Z","shell.execute_reply.started":"2022-10-09T16:12:22.325695Z","shell.execute_reply":"2022-10-09T16:12:22.398003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.target.value_counts()","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:12:22.401041Z","iopub.execute_input":"2022-10-09T16:12:22.401533Z","iopub.status.idle":"2022-10-09T16:12:22.416896Z","shell.execute_reply.started":"2022-10-09T16:12:22.401493Z","shell.execute_reply":"2022-10-09T16:12:22.415806Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> label 1 when a GW signal is present in the data.\n\n> label 0 is when a GW signal is not present in the sample.\n\n> label -1 is for samples where the Physicists are currently unable to determine the status.","metadata":{}},{"cell_type":"markdown","source":"## Read HDF5 File","metadata":{}},{"cell_type":"code","source":"train_files = glob(f\"{ROOT_DIR}/train/*.hdf5\")\ntrain_files[:5]","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:12:22.418481Z","iopub.execute_input":"2022-10-09T16:12:22.41918Z","iopub.status.idle":"2022-10-09T16:12:22.660157Z","shell.execute_reply.started":"2022-10-09T16:12:22.419132Z","shell.execute_reply":"2022-10-09T16:12:22.659235Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATA_IDX = 0","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:12:22.661653Z","iopub.execute_input":"2022-10-09T16:12:22.662535Z","iopub.status.idle":"2022-10-09T16:12:22.668684Z","shell.execute_reply.started":"2022-10-09T16:12:22.662468Z","shell.execute_reply":"2022-10-09T16:12:22.66718Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"file = Path(train_files[DATA_IDX])\nwith h5py.File(file, \"r\") as f:\n    # Print all root level object names (aka keys) \n    # these can be group or dataset names \n    print(\"Top level Keys: %s\" % f.keys())\n    key = list(f.keys())[0]\n    print(\"This is the name of the file \\n\")\n    # get the object type for a_group_key: usually group or dataset\n    print(f\"Type of the data accessed using {key}: {type(f[key])}\")\n    print(\"The type can be either Group or Dataset. Group has more keys while we can get numpy array from Dataset directly.\\n\")\n    # If a_group_key is a group name, \n    # this gets the object names in the group and returns as a list\n    f = f[key]\n    data_keys = list(f)\n    print(\"These are the keys for the underlying group: \", data_keys, \"\\n\")\n    # get type of each data_keys\n    print(type(f[data_keys[0]])) # Group\n    print(type(f[data_keys[1]])) # Group\n    print(type(f[data_keys[2]])) # Dataset\n    print(\"The underlying data type for each key are Group, Group and Dataset. We thus need to go one level deep to access the data from Group. \\n\")\n    # H1\n    h1 = f[data_keys[0]]\n    h1_keys = list(h1)\n    print(\"These are the keys for the Group accessed uisng H1 key: \", h1_keys, \"\\n\")\n    # L1\n    l1 = f[data_keys[1]]\n    l1_keys = list(l1)\n    print(\"These are the keys for the Group accessed uisng L1 key: \", l1_keys, \"\\n\")\n    # frequency_Hz\n    freq_hz = f[data_keys[2]]\n    freq_hz = list(freq_hz)\n    print(\"Since the type accessed using frequency_hz key was Dataset, we can get the array directly. The length of the array is: \", len(freq_hz), \"\\n\")\n    # H1 data\n    h1_stft = h1[\"SFTs\"][()]\n    h1_timestamp = h1[\"timestamps_GPS\"][()]\n    print(\"The H1 SFTs data shape: \", h1_stft.shape)\n    print(\"The H1 timestamp data shape: \", h1_timestamp.shape, \"\\n\")\n    # H2 data\n    l1_stft = l1[\"SFTs\"][()]\n    l1_timestamp = l1[\"timestamps_GPS\"][()]\n    print(\"The L1 SFTs data shape: \", l1_stft.shape)\n    print(\"The L1 timestamp data shape: \", l1_timestamp.shape, \"\\n\")\n    print(\"We have time on the x-axis while frequency on the y-axis. \\n\")\n    \n    print(\"Label: \", df.loc[DATA_IDX].target)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-10-09T16:12:22.671309Z","iopub.execute_input":"2022-10-09T16:12:22.671864Z","iopub.status.idle":"2022-10-09T16:12:23.263834Z","shell.execute_reply.started":"2022-10-09T16:12:22.671812Z","shell.execute_reply":"2022-10-09T16:12:23.262421Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Utility to read hdf5 file\ndef read_data(file: Path):\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 {\n            \"H1\": [h1_stft, h1_timestamp],\n            \"L1\": [l1_stft, l1_timestamp],\n            \"freq_hz\": freq_hz\n        }","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:12:23.265472Z","iopub.execute_input":"2022-10-09T16:12:23.265964Z","iopub.status.idle":"2022-10-09T16:12:23.275313Z","shell.execute_reply.started":"2022-10-09T16:12:23.265896Z","shell.execute_reply":"2022-10-09T16:12:23.274085Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"file = Path(train_files[1])\ndata = read_data(file)\nh1_sft, h1_ts = data[\"H1\"]\nl1_sft, l1_ts = data[\"L1\"]\nfreq_hz = data[\"freq_hz\"]","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:12:23.278838Z","iopub.execute_input":"2022-10-09T16:12:23.279248Z","iopub.status.idle":"2022-10-09T16:12:23.944791Z","shell.execute_reply.started":"2022-10-09T16:12:23.279215Z","shell.execute_reply":"2022-10-09T16:12:23.943471Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Understanding shapes of the data gives a nice insight.","metadata":{}},{"cell_type":"code","source":"print(f\"The shape of H1 SFT is: {h1_sft.shape} while that of it's time stamp is: {h1_ts.shape}\")\nprint(f\"The shape of L1 SFT is: {l1_sft.shape} while that of it's time stamp is: {l1_ts.shape}\")\nprint(f\"The shape of frequency Hz is: {np.array(freq_hz).shape}\")","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:12:23.946146Z","iopub.execute_input":"2022-10-09T16:12:23.947086Z","iopub.status.idle":"2022-10-09T16:12:23.953139Z","shell.execute_reply.started":"2022-10-09T16:12:23.947046Z","shell.execute_reply":"2022-10-09T16:12:23.952165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> The x-axis is the time stamp. Note however that the time stamp may not be continuous. From the data description: \"The SFTs are not always contiguous in time, since the interferometers are not continuously online.\"\n\n> The H1 data corresponds to LIGO Hanford while L1 corresponds to LIGO Livingston.\n\n> The frequency in Hz is measured from 0-360 Hz. Not sure if this is uniform throughtout the dataset.\n\n> The timestamps between the two interferometers (LIGOs) are not the same.","metadata":{}},{"cell_type":"markdown","source":"# H1: LIGO Hanford","metadata":{}},{"cell_type":"markdown","source":"## Time stamp","metadata":{}},{"cell_type":"code","source":"h1_ts","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:12:23.954522Z","iopub.execute_input":"2022-10-09T16:12:23.954995Z","iopub.status.idle":"2022-10-09T16:12:23.966968Z","shell.execute_reply.started":"2022-10-09T16:12:23.954947Z","shell.execute_reply":"2022-10-09T16:12:23.965706Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> The timestamps are in UNIX epoch format. They can be converted to human readable time.","metadata":{}},{"cell_type":"code","source":"from datetime import datetime\nprint(datetime.utcfromtimestamp(h1_ts[0]).strftime('%Y-%m-%d %H:%M:%S'))","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:12:23.968704Z","iopub.execute_input":"2022-10-09T16:12:23.969161Z","iopub.status.idle":"2022-10-09T16:12:23.977157Z","shell.execute_reply.started":"2022-10-09T16:12:23.969118Z","shell.execute_reply":"2022-10-09T16:12:23.97578Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> Thus this particular signal started recording on 27th March, 2009 at 3:07 PM GMT.","metadata":{}},{"cell_type":"markdown","source":"There are gaps in the timestamp but how wild is the gap?","metadata":{}},{"cell_type":"code","source":"gaps = []\na = h1_ts[0]\nfor b in h1_ts[1:]:\n    gap = b-a\n    a = b\n    gaps.append(gap)\n\nplt.figure(figsize=(20, 5))\nplt.plot(h1_ts[1:], gaps, marker='.');","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:47:52.421415Z","iopub.execute_input":"2022-10-09T16:47:52.42246Z","iopub.status.idle":"2022-10-09T16:47:52.665521Z","shell.execute_reply.started":"2022-10-09T16:47:52.422418Z","shell.execute_reply":"2022-10-09T16:47:52.664289Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> Clearly there are a lot of gaps. Each spike is a gap in the timestamp. The height of each spike is the gap in UNIX epoch. The bold horizontal line is the gap that repeats multiple times, in this case it is 1800 UNIX epoch. Basically, the timestamp and thus STFs are not continuous reading (what a bummer?) and the gaps can be wild.","metadata":{}},{"cell_type":"markdown","source":"What about the timestamp of L1?","metadata":{}},{"cell_type":"code","source":"len(l1_ts) == len(h1_ts)","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:12:24.125469Z","iopub.execute_input":"2022-10-09T16:12:24.12929Z","iopub.status.idle":"2022-10-09T16:12:24.145452Z","shell.execute_reply.started":"2022-10-09T16:12:24.129217Z","shell.execute_reply":"2022-10-09T16:12:24.144157Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> The length of both the timestamps are different. That's expected since both the detectors might not be turned off at the same time.","metadata":{}},{"cell_type":"code","source":"l1_gaps = []\na = l1_ts[0]\nfor b in l1_ts[1:]:\n    gap = b-a\n    a = b\n    l1_gaps.append(gap)\n    \nplt.figure(figsize=(20, 5))\nplt.plot(l1_ts[1:], l1_gaps, marker='.');","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:47:57.629196Z","iopub.execute_input":"2022-10-09T16:47:57.630062Z","iopub.status.idle":"2022-10-09T16:47:57.881578Z","shell.execute_reply.started":"2022-10-09T16:47:57.630011Z","shell.execute_reply":"2022-10-09T16:47:57.88024Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> Clearly the gaps are different.","metadata":{}},{"cell_type":"code","source":"h1_ts[0] == l1_ts[0]","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:12:24.268338Z","iopub.execute_input":"2022-10-09T16:12:24.269227Z","iopub.status.idle":"2022-10-09T16:12:24.278533Z","shell.execute_reply.started":"2022-10-09T16:12:24.269178Z","shell.execute_reply":"2022-10-09T16:12:24.276997Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> The detectors are not taking measurements starting at the same time as well.","metadata":{}},{"cell_type":"markdown","source":"### Utilities\n\nThe following code is copied from [`tutorial_utils.py`](https://github.com/PyFstat/PyFstat/blob/master/examples/tutorials/tutorial_utils.py).","metadata":{}},{"cell_type":"code","source":"\"\"\"\nUtils\n=====\n\nUtility functions to simplify tutorials.\n\"\"\"\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib import colors\nfrom scipy import stats\n\nplt.rcParams[\"font.family\"] = \"serif\"\nplt.rcParams[\"font.size\"] = 20\n\n\ndef plot_real_imag_spectrograms(timestamps, frequency, fourier_data):\n    \"uses simple plt.pcolormesh\"\n    fig, axs = plt.subplots(1, 2, figsize=(16, 10))\n\n    for ax in axs:\n        ax.set(xlabel=\"SFT index\", ylabel=\"Frequency [Hz]\")\n\n    time_in_days = (timestamps - timestamps[0]) / 1800\n\n    axs[0].set_title(\"SFT Real part\")\n    c = axs[0].pcolormesh(\n        time_in_days,\n        frequency,\n        fourier_data.real,\n        norm=colors.CenteredNorm(),\n    )\n    fig.colorbar(c, ax=axs[0], orientation=\"horizontal\", label=\"Fourier Amplitude\")\n\n    axs[1].set_title(\"SFT Imaginary part\")\n    c = axs[1].pcolormesh(\n        time_in_days,\n        frequency,\n        fourier_data.imag,\n        norm=colors.CenteredNorm(),\n    )\n\n    fig.colorbar(c, ax=axs[1], orientation=\"horizontal\", label=\"Fourier Amplitude\")\n\n    return fig, axs\n\n\ndef plot_real_imag_spectrograms_with_gaps(timestamps, frequency, fourier_data, Tsft):\n\n    # Fill up gaps with Nans\n    gap_length = timestamps[1:] - (timestamps[:-1] + Tsft)\n\n    gap_data = [fourier_data[:, 0]]\n    gap_timestamps = [timestamps[0]]\n\n    for ind, gap in enumerate(gap_length):\n        if gap > 0:\n            gap_data.append(np.full_like(fourier_data[:, ind], np.nan + 1j * np.nan))\n            gap_timestamps.append(timestamps[ind] + Tsft)\n\n        gap_data.append(fourier_data[:, ind + 1])\n        gap_timestamps.append(timestamps[ind + 1])\n\n    return plot_real_imag_spectrograms(\n        np.hstack(gap_timestamps), frequency, np.vstack(gap_data).T\n    )\n\n\ndef plot_real_imag_histogram(fourier_data, theoretical_stdev=None):\n\n    fig, ax = plt.subplots(figsize=(16, 10))\n    ax.set(xlabel=\"SFT value\", ylabel=\"PDF\", yscale=\"log\")\n\n    ax.hist(\n        fourier_data.real.ravel(),\n        density=True,\n        bins=\"auto\",\n        histtype=\"step\",\n        lw=2,\n        label=\"Real part\",\n    )\n    ax.hist(\n        fourier_data.imag.ravel(),\n        density=True,\n        bins=\"auto\",\n        histtype=\"step\",\n        lw=2,\n        label=\"Imaginary part\",\n    )\n\n    if theoretical_stdev is not None:\n        x = np.linspace(-4 * theoretical_stdev, 4 * theoretical_stdev, 1000)\n        y = stats.norm(scale=theoretical_stdev).pdf(x)\n        ax.plot(x, y, color=\"black\", ls=\"--\", label=\"Gaussian distribution\")\n\n    ax.legend()\n\n    return fig, ax\n\n\ndef plot_amplitude_phase_spectrograms(timestamps, frequency, fourier_data):\n    fig, axs = plt.subplots(1, 2, figsize=(16, 10))\n\n    for ax in axs:\n        ax.set(xlabel=\"SFT index\", ylabel=\"Frequency [Hz]\")\n\n    time_in_days = (timestamps - timestamps[0]) / 1800\n\n    axs[0].set_title(\"SFT absolute value\")\n    c = axs[0].pcolorfast(\n        time_in_days, frequency, np.absolute(fourier_data), norm=colors.Normalize()\n    )\n    fig.colorbar(c, ax=axs[0], orientation=\"horizontal\", label=\"Value\")\n\n    axs[1].set_title(\"SFT phase\")\n    c = axs[1].pcolorfast(\n        time_in_days, frequency, np.angle(fourier_data), norm=colors.CenteredNorm()\n    )\n\n    fig.colorbar(c, ax=axs[1], orientation=\"horizontal\", label=\"Value\")\n\n    return fig, axs","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-10-09T16:48:11.444218Z","iopub.execute_input":"2022-10-09T16:48:11.445172Z","iopub.status.idle":"2022-10-09T16:48:11.47146Z","shell.execute_reply.started":"2022-10-09T16:48:11.445112Z","shell.execute_reply":"2022-10-09T16:48:11.470058Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## STFT\n\nShort time fourier transform computes the fourier transform of the signal at shorter duration of signals given by a window size. The DFT (discrete fourier transform) or FFT (fast fourier transform) are computed by taking the entire duration of the signal into consideration.\n\nNote that the plots generated below are not considering the gaps in the timestamp and considers the data as a single continuous SFT.\n\n","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(2, 2, figsize=(16, 10))\nfig.suptitle(f\"Filename Id: {file.stem}\")\n\nfor d_ind, detector in enumerate([\"H1\", \"L1\"]):\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(data[detector][1], data[\"freq_hz\"], \n                                 data[detector][0].real)\n    c1 = ax[d_ind][1].pcolormesh(data[detector][1], data[\"freq_hz\"], \n                                 data[detector][0].imag)\n\n    fig.colorbar(c0, ax=ax[d_ind][0])\n    fig.colorbar(c1, ax=ax[d_ind][1])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:48:01.662403Z","iopub.execute_input":"2022-10-09T16:48:01.662803Z","iopub.status.idle":"2022-10-09T16:48:06.324917Z","shell.execute_reply.started":"2022-10-09T16:48:01.662771Z","shell.execute_reply":"2022-10-09T16:48:06.323986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Plotting the spectrogram by considering the gaps in SFTs gives a plot like this which is more meaningful. Note that I am providing the \"Tsft\" argument of this function as 1800. This is what the authors of PyFstat have recommened in the tutorials and by looking at the gaps ourself, 1800 is what appears multiple times.","metadata":{}},{"cell_type":"code","source":"plot_real_imag_spectrograms_with_gaps(data[detector][1], data[\"freq_hz\"], data[detector][0], 1800);","metadata":{"execution":{"iopub.status.busy":"2022-10-09T21:37:46.535021Z","iopub.execute_input":"2022-10-09T21:37:46.535536Z","iopub.status.idle":"2022-10-09T21:37:50.934776Z","shell.execute_reply.started":"2022-10-09T21:37:46.535496Z","shell.execute_reply":"2022-10-09T21:37:50.932766Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Understanding By Generating\n\nThe data can better be understood if we generate the data and see how different parameters contribute to the SFTs. This section is based on the awesome tutorials provided by the host of this competition. You can find the tutorials [here](https://github.com/PyFstat/PyFstat/tree/master/examples/tutorials). The following is a distilled information in an attemp to better understand the data.\n\n----\n\n### TLDR;\n\n**The SFTs generated with `h0` parameter equal to 0 have target/label as 0. The SFTs generated with `h0` parameter equal to 1, thus having SNR (signal to noise ratio) > 0 have it's target/label as 1.** The competition hosts were kind enough to clarify this [here](https://www.kaggle.com/competitions/g2net-detecting-continuous-gravitational-waves/discussion/347052#1973327).\n\nThe main piece of code to generate the data is `writer = pyfstat.Writer(**writer_kwargs, **signal_parameters)`. If you don't pass `**signal_parameter` you will get data corresponding to 0 and so on. You can check the `Writer` class [here](https://github.com/PyFstat/PyFstat/blob/b173b3d6e39088fa3f50110033e6c7ec2c1d6f54/pyfstat/make_sfts.py#L21).\n\n----","metadata":{}},{"cell_type":"markdown","source":"The CW wave has the following properties:\n\n* Quasi-monochromatic: \"Wave in which most of the energy is confined to a single wavelength or a very narrow waveband\". However, gravitational continuous waves are not completely monochromatic. ([source](http://www.idigitalphoto.com/dictionary/quasi-monochromatic))\n* Long-standing: The standing waves are stationary (not a travelling wave). The disturbance does not travel in any direction. Standing waves have points of zero amplitude called nodes and points of maximum amplitude called the antinodes. \n* Produced by non-axisymmetric rapidly-spinning neutron stars\n* Produce narrow-banded signals - that's why the y axis (Freq in Hz) in the spectrogram lies in a narrow band of 1-2 Hz.\n* The waves are doppler-modulated due to the motion of the detector in the Solar system.\n\nIt's a common practice to use short time fourier transform to analyze these waves. These are Fourier transforms of short data segments (typically around 30 minutes or **1800 seconds**).","metadata":{}},{"cell_type":"markdown","source":"### Generate Gausian Noise","metadata":{}},{"cell_type":"code","source":"# Setup Writer\nwriter_kwargs = {\n    \"label\": \"single_detector_gaussian_noise\",\n    \"outdir\": \"PyFstat_example_data\",\n    \"tstart\": 1238166018,  # Starting time of the observation [GPS time]\n    \"duration\": 5 * 86400,  # Duration [seconds]\n    \"detectors\": \"H1\",  # Detector to simulate, in this case LIGO Hanford\n    \"F0\": 100.0,  # Central frequency of the band to be generated [Hz]\n    \"Band\": 1.0,  # Frequency band-width around F0 [Hz]\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}\nwriter = pyfstat.Writer(**writer_kwargs)\n\n# Create SFTs\nwriter.make_data()","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:54:17.767342Z","iopub.execute_input":"2022-10-09T16:54:17.767844Z","iopub.status.idle":"2022-10-09T16:54:17.782538Z","shell.execute_reply.started":"2022-10-09T16:54:17.767811Z","shell.execute_reply":"2022-10-09T16:54:17.781362Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Note: \"As a side note, the path where the SFTs are created is stored in the sftfilepath attribute of Writer\n(mind that it gets overwritten every time make_data is called!).\" ([source](https://www.kaggle.com/competitions/g2net-detecting-continuous-gravitational-waves/discussion/347052#1976378))","metadata":{}},{"cell_type":"code","source":"frequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:54:19.851471Z","iopub.execute_input":"2022-10-09T16:54:19.852192Z","iopub.status.idle":"2022-10-09T16:54:19.870526Z","shell.execute_reply.started":"2022-10-09T16:54:19.852151Z","shell.execute_reply":"2022-10-09T16:54:19.868943Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"h1_sft, h1_ts = fourier_data[\"H1\"], timestamps[\"H1\"]\nprint(f\"The shape of H1 SFT is: {h1_sft.shape} while that of it's time stamp is: {h1_ts.shape}\")\nprint(f\"The shape of frequency Hz is: {frequency.shape}\")","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:54:23.781761Z","iopub.execute_input":"2022-10-09T16:54:23.782223Z","iopub.status.idle":"2022-10-09T16:54:23.790555Z","shell.execute_reply.started":"2022-10-09T16:54:23.782187Z","shell.execute_reply":"2022-10-09T16:54:23.788288Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Why the shape of time stamp is (240,)?\n\nThere are three time related parameters:\n```\n\"tstart\": 1238166018,  # Starting time of the observation [GPS time]\n\"duration\": 5 * 86400,  # Duration [seconds]\n\"Tsft\": 1800,  # Fourier transform time duration\n```","metadata":{}},{"cell_type":"code","source":"print(f\"Total duration is : {5*86400}; time segments for SFT computation: {1800}; total timestamps: {(5*86400)/1800}\")","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:54:29.659656Z","iopub.execute_input":"2022-10-09T16:54:29.660383Z","iopub.status.idle":"2022-10-09T16:54:29.666247Z","shell.execute_reply.started":"2022-10-09T16:54:29.660342Z","shell.execute_reply":"2022-10-09T16:54:29.665205Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_real_imag_spectrograms(\n    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:54:32.020999Z","iopub.execute_input":"2022-10-09T16:54:32.021898Z","iopub.status.idle":"2022-10-09T16:54:33.6308Z","shell.execute_reply.started":"2022-10-09T16:54:32.02185Z","shell.execute_reply":"2022-10-09T16:54:33.628989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**The SNR for such a signal is 0 since there are no signal.** :P","metadata":{}},{"cell_type":"markdown","source":"Let's inspect the apmplitude and angle of the generate SFT.","metadata":{}},{"cell_type":"code","source":"fig, ax = plot_amplitude_phase_spectrograms(\n    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-09T16:55:50.481823Z","iopub.execute_input":"2022-10-09T16:55:50.48246Z","iopub.status.idle":"2022-10-09T16:55:51.423058Z","shell.execute_reply.started":"2022-10-09T16:55:50.482416Z","shell.execute_reply":"2022-10-09T16:55:51.421502Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Since it's a gaussian noise, the power (amplitude) and angle in the SFT is random.","metadata":{}},{"cell_type":"markdown","source":"### Generate CW Signal\n\nThe CW waves are primarily dependent on: `amplitude` and `Doppler` parameters.\n\n- the `amplitide` parameter depends on nominal CW amplitude $h0$, the cosine of the source inclination angle with respect to the line of sight $cos(ι)$ (cosi), the polarization angle $\\psi{}$ (psi) and the initial phase of the wave $\\phi{}$ (phi). Depending on the emission mechanism, $h0$ can be further described using further physical quantities such as the source's frequency, ellipticity, or distance to the detector.\n\n- the `doppler` parameter describe the relationship between the source of emmission (and it's motion) and the detectors. The parameter depends on gravitational wave frequency $f0$, spindown $f1$and sky position $\\hat{n}$, which we will parametrize using the right ascension $\\alpha{}$ and declination $\\delta{}$ angles.\n","metadata":{}},{"cell_type":"code","source":"# The same writer params are required to generate gaussian noise.\nwriter_kwargs = {\n    \"label\": \"single_detector_gaussian_noise\",\n    \"outdir\": \"PyFstat_example_data\",\n    \"tstart\": 1238166018,\n    \"duration\": 365 * 86400,\n    \"detectors\": \"H1\",\n    \"sqrtSX\": 1e-23,\n    \"Tsft\": 1800,\n    \"SFTWindowType\": \"tukey\",\n    \"SFTWindowBeta\": 0.01,\n}\n\n# These parameters determine the CW signal that we have to detect.\nsignal_parameters = {\n    \"F0\": 100.0,\n    \"F1\": -1e-9,\n    \"Alpha\": 0.0,\n    \"Delta\": 0.0,\n    \"h0\": 1e-22, # this parameter in particular determines if a signal exist in the generated SFT.\n    \"cosi\": 1,\n    \"psi\": 0.0,\n    \"phi\": 0.0,\n    \"tref\": writer_kwargs[\"tstart\"],\n}\n\nwriter = pyfstat.Writer(**writer_kwargs, **signal_parameters)\nwriter.make_data()\nfrequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)\nplot_real_imag_spectrograms(\n    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-09T21:49:30.950613Z","iopub.execute_input":"2022-10-09T21:49:30.951118Z","iopub.status.idle":"2022-10-09T21:49:43.534337Z","shell.execute_reply.started":"2022-10-09T21:49:30.951079Z","shell.execute_reply":"2022-10-09T21:49:43.533054Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> In this case the signal is really clear and thus SNR should be high.","metadata":{}},{"cell_type":"code","source":"fig, ax = plot_amplitude_phase_spectrograms(\n    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-09T17:02:53.325086Z","iopub.execute_input":"2022-10-09T17:02:53.329277Z","iopub.status.idle":"2022-10-09T17:02:54.814682Z","shell.execute_reply.started":"2022-10-09T17:02:53.32915Z","shell.execute_reply":"2022-10-09T17:02:54.813345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The power and the angle component of the generate SFT is prominent.","metadata":{}},{"cell_type":"code","source":"# SNR can be compute from a set of SFTs for a specific set\n# of parameters as follows:\nsnr = pyfstat.SignalToNoiseRatio.from_sfts(\n    F0=writer.F0, sftfilepath=writer.sftfilepath\n)\nsquared_snr = snr.compute_snr2(\n    Alpha=writer.Alpha, \n    Delta=writer.Delta,\n    psi=writer.psi,\n    phi=writer.phi, \n    h0=writer.h0,\n    cosi=writer.cosi\n)\nsqrt_snr = np.sqrt(squared_snr)\nprint(f\"Signal - SNR: {sqrt_snr}\")","metadata":{"execution":{"iopub.status.busy":"2022-10-09T17:24:43.666839Z","iopub.execute_input":"2022-10-09T17:24:43.667398Z","iopub.status.idle":"2022-10-09T17:24:45.702827Z","shell.execute_reply.started":"2022-10-09T17:24:43.667351Z","shell.execute_reply":"2022-10-09T17:24:45.701943Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> The data of interest (real data) should have such high SNRs (there would be no point of this competition, if signal is easily detected).","metadata":{}},{"cell_type":"markdown","source":"Let's talk about the possible SNR range of CW signals. As per the tutorial shared by competition hosts, \"The presence of visible CW features in the spectrogram can be quantified using the so-called sensitivity depth, which is defined as the ratio between the noise and CW amplitude $D = \\sqrt{Sn}/h0$. Current searches for CW signals from unknown sources are able to detect a sensitivity depth ranging between 10 $Hz^{-0.5}$ and 50 $Hz^{-0.5}$. Visual features, however, tend to disappear around a depth of 20 $Hz^{-0.5}$.\"\n\nLooking at the formula, $\\sqrt{Sn}$ is \"sqrtSX\" which in our case is `1e-23`. And we want to look at visual features where the sensitivity depth is 20 $Hz^{-0.5}$. Thus we will calculate $h0$.","metadata":{}},{"cell_type":"code","source":"# CW signal at a depth of 20\nsignal_parameters[\"h0\"] = writer_kwargs[\"sqrtSX\"] / 20.0\n\nwriter = pyfstat.Writer(**writer_kwargs, **signal_parameters)\n\n# Create SFTs\nwriter.make_data()\n\nfrequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)\n\nplot_real_imag_spectrograms(\n    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n)\n\nplot_amplitude_phase_spectrograms(\n    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-09T20:08:42.524018Z","iopub.execute_input":"2022-10-09T20:08:42.524738Z","iopub.status.idle":"2022-10-09T20:08:52.791337Z","shell.execute_reply.started":"2022-10-09T20:08:42.524688Z","shell.execute_reply":"2022-10-09T20:08:52.789967Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> We really can't see any visual feature in the spectrogram.","metadata":{}},{"cell_type":"code","source":"# SNR can be compute from a set of SFTs for a specific set\n# of parameters as follows:\nsnr = pyfstat.SignalToNoiseRatio.from_sfts(\n    F0=writer.F0, sftfilepath=writer.sftfilepath\n)\nsquared_snr = snr.compute_snr2(\n    Alpha=writer.Alpha,\n    Delta=writer.Delta,\n    psi=writer.psi,\n    phi=writer.phi,\n    h0=writer.h0,\n    cosi=writer.cosi\n)\nsqrt_snr = np.sqrt(squared_snr)\nprint(f\"Signal - SNR: {sqrt_snr}\")","metadata":{"execution":{"iopub.status.busy":"2022-10-09T20:09:16.607648Z","iopub.execute_input":"2022-10-09T20:09:16.608129Z","iopub.status.idle":"2022-10-09T20:09:18.695032Z","shell.execute_reply.started":"2022-10-09T20:09:16.608091Z","shell.execute_reply":"2022-10-09T20:09:18.693757Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"But the SNR > 0 indicates that our target signal (CW waves) is present and we need better system to predict them.","metadata":{}},{"cell_type":"markdown","source":"### So what about the gaps in the timestamp?\n\nThe gaps are becasue the detectors are not turned on 24x7. They need to go into maintainence or something else, and since we have two detectors they might not have readings at the dame time.\n\nLet's see how we can generate SFTs with CW waves with gaps following the tutorials. Note that I am purposefully generating SFT with higher SNR so that it's easy to visualize the CW waves features.","metadata":{}},{"cell_type":"code","source":"# We can provide different timestamps using a dict with values indicating the name of the detector.\n# We are using the timestamp of the H1 detector taken from the first row of the provided csv file.\ntimestamps = {\"H1\": data[detector][1]}\n\n# Setup Writer\nwriter_kwargs = {\n    \"label\": \"single_detector_gaps\",\n    \"outdir\": \"PyFstat_example_data\",\n    \"timestamps\": timestamps,\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}\n\n# These parameters determine the CW signal that we have to detect.\nsignal_parameters = {\n    \"F0\": 100.0,\n    \"F1\": -1e-9,\n    \"Alpha\": 0.0,\n    \"Delta\": 0.0,\n    \"h0\": 1e-22,\n    \"cosi\": 1,\n    \"psi\": 0.0,\n    \"phi\": 0.0,\n}\n\nwriter = pyfstat.Writer(**writer_kwargs, **signal_parameters)\n# writer = pyfstat.Writer(**writer_kwargs)\n\n# Create SFTs\nwriter.make_data()\nfrequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)","metadata":{"execution":{"iopub.status.busy":"2022-10-09T21:51:06.778216Z","iopub.execute_input":"2022-10-09T21:51:06.778851Z","iopub.status.idle":"2022-10-09T21:51:09.998759Z","shell.execute_reply.started":"2022-10-09T21:51:06.7788Z","shell.execute_reply":"2022-10-09T21:51:09.997013Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_real_imag_spectrograms_with_gaps(\n    timestamps[\"H1\"], frequency, fourier_data[\"H1\"], writer_kwargs[\"Tsft\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-09T21:51:19.592316Z","iopub.execute_input":"2022-10-09T21:51:19.592815Z","iopub.status.idle":"2022-10-09T21:51:22.105073Z","shell.execute_reply.started":"2022-10-09T21:51:19.592771Z","shell.execute_reply":"2022-10-09T21:51:22.103952Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# WIP","metadata":{}}]}