{"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":"# How to generate continuous gravitational-wave signals\n\nThis notebook presents a quick hands-on introduction on generating your very own set of continuous gravitational-wave signals (CWs).\n\nWe recommend having a look at [the PyFstat tutorials](https://www.kaggle.com/competitions/g2net-detecting-continuous-gravitational-waves/discussion/347052) for a more in-depth introduction to the topic.\n\n1. [Parameters](#1)\n2. [Generating data](#2)\n3. [Generating specific frequency bands](#3)","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 git+https://github.com/PyFstat/PyFstat@python37","metadata":{"execution":{"iopub.status.busy":"2022-10-15T07:18:34.105673Z","iopub.execute_input":"2022-10-15T07:18:34.106871Z","iopub.status.idle":"2022-10-15T07:19:17.181317Z","shell.execute_reply.started":"2022-10-15T07:18:34.106742Z","shell.execute_reply":"2022-10-15T07:19:17.180332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport sys\n\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nimport pyfstat\n\nfrom scipy import stats\n\n%matplotlib inline","metadata":{"execution":{"iopub.status.busy":"2022-10-15T07:20:47.330646Z","iopub.execute_input":"2022-10-15T07:20:47.331597Z","iopub.status.idle":"2022-10-15T07:20:47.339005Z","shell.execute_reply.started":"2022-10-15T07:20:47.331545Z","shell.execute_reply":"2022-10-15T07:20:47.337826Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"1\"></a> <br>\n# 1. Parameters\n\nStandard CW signals can be parameterised in terms of two sets of parameters: the Doppler-modulation parameters $\\lambda$ and the amplitude parameters $\\mathcal{A}$.\nThe former encode how the frequency of a signal modulates due to its intrinsic frequency evolution and the movement of the Earth in the Solar system, \nwhile the latter describes the overall amplitude of a CW depending on the parameters of the source.\n\nFor a CW emmitted by a rapidly-spinning and isolated neutron star (NS), Doppler-modulation parameters include the frequency `F0` and the linear spindown parameter `F1`,\nboth taken at a reference time `tref`, and the sky position in terms of the right ascension `Alpha` and declination `Delta` angles of equatorial cordinates. \nAmplitude parameters, on the other hand, include the average amplitude of a CW signal `h0`, the initial phase of the signal `phi`, the polarization angle `psi`\nand (the cosine of) the inclination angle of the source `cosi`, which gives us the relative orientation of the NS with respect to the detector.\n\nAs described in [the signal tutorial](https://github.com/PyFstat/PyFstat/blob/master/examples/tutorials/1_generating_signals.ipynb),\nthe amplitude of a CW signal is usually expressed in terms of the noises's amplitude using depth $\\mathcal{D}$ or signal-to-noise ratio (SNR) $\\rho$.\nFor our purposes, the former is essentially a quotient\n$$\n\\mathcal{D} = \\frac{\\sqrt{S_{\\mathrm{n}}}}{h_0}\n$$\nwhile the latter is a more involved expression wich also depends on the duration of the dataset at hand and the detector's response function.\nIt is important to note, however, that $\\rho$ and $\\mathcal{D}$ scale reciprocally: \"weak\" signals have a *low* SNR and a *high* depth (since they are \"buried deeper into the noise\" than a strong signal).","metadata":{}},{"cell_type":"markdown","source":"<a id=\"2\"></a> <br>\n# 2. Generating data\n\nA specific sample requires of background noise and optionally a signal. \nIn order to generate noise, one needs to specify a set of `detectors` (`H1` or `L1` in this case), the duration of the sample and the Amplitude Spectral Density of the noise `sqrtSX`. \nCW analyses are simple in this front, as `sqrtSX` is proportional to the (stationary) standard deviation of an underlying zero-mean Gaussian process.\n\nSample duration can be specified in two ways. If the sample contains contiguous data (i.e. the detector was taking science-quality data uninterrupted), one can simply specify the starting time and duration of the sample using `tstart` and `duration`. Data with gaps, on the other hand, can be generated by specifying a specific set of timestamps using the `timestamps` option.\n\nData is saved as a list of Short Fourier Transforms (SFTs). The duration and windowing of these SFTs can also be modified using `Tsft`, `SFTWindowType` and `SFTWindowBeta`.\nMost analyses tune `Tsft` around 1800 seconds order to ensure the power of a putative CW signal stays within a bin.","metadata":{}},{"cell_type":"code","source":"# Generate signals with parameters drawn from a specific population\nnum_signals = 5\n\n# These parameters describe background noise and data format\nwriter_kwargs = {\n                \"tstart\": 1238166018,\n                \"duration\": 4 * 30 * 86400,  \n                \"detectors\": \"H1,L1\",        \n                \"sqrtSX\": 1e-23,          \n                \"Tsft\": 1800,             \n                \"SFTWindowType\": \"tukey\", \n                \"SFTWindowBeta\": 0.01,\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.\nsignal_parameters_generator = pyfstat.AllSkyInjectionParametersGenerator(\n    priors={\n        \"tref\": writer_kwargs[\"tstart\"],\n        \"F0\": {\"uniform\": {\"low\": 100.0, \"high\": 100.1}},\n        \"F1\": lambda: -10**stats.uniform(-12, 4).rvs(),\n        \"F2\": 0,\n        \"h0\": lambda: writer_kwargs[\"sqrtSX\"] / stats.uniform(1, 10).rvs(),\n        **pyfstat.injection_parameters.isotropic_amplitude_priors,\n    },\n)\n\nsnrs = np.zeros(num_signals)\n\nfor ind in range(num_signals):\n\n    # Draw signal parameters.\n    # Noise can be drawn by setting `params[\"h0\"] = 0\n    params = signal_parameters_generator.draw()\n    writer_kwargs[\"outdir\"] = f\"PyFstat_example_data_ensemble/Signal_{ind}\"\n    writer_kwargs[\"label\"] = f\"Signal_{ind}\"\n    \n    writer = pyfstat.Writer(**writer_kwargs, **params)\n    writer.make_data()\n    \n    # SNR can be compute from a set of SFTs for a specific set\n    # of parameters as follows:\n    snr = pyfstat.SignalToNoiseRatio.from_sfts(\n        F0=writer.F0, sftfilepath=writer.sftfilepath\n    )\n    squared_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    )\n    snrs[ind] = np.sqrt(squared_snr)\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    fig, ax = plt.subplots(2, 2, figsize=(16, 10))\n    fig.suptitle(f\"Signal {ind} - SNR: {snrs[ind]:.2f}\")\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(timestamps[detector], frequency,\n                                     amplitudes[detector].real)\n        c1 = ax[d_ind][1].pcolormesh(timestamps[detector], frequency,\n                                     amplitudes[detector].imag)\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-10-15T07:20:48.803501Z","iopub.execute_input":"2022-10-15T07:20:48.803915Z","iopub.status.idle":"2022-10-15T07:23:16.301547Z","shell.execute_reply.started":"2022-10-15T07:20:48.803884Z","shell.execute_reply":"2022-10-15T07:23:16.29994Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"3\"></a> <br>\n# 3. Generating specific frequency bands\n\nThe `pyfstat.Writer` class generates CW signals following the conventions used by LIGO data-analysis codes; more specifically, frequency bands are always made wide enough so that any requested CW signal can fit into it and there are enough extra bins to compute running quantities. This may be a problem, as in some cases data with a very specific frequency band may need to be generated.\n\nA simple way to generate specific output would be to generate a broad SFT and slice the frequency band of interest whenever reading as a numpy array.\nThis can be done as follows (taken from [this conversation with DennisSakva](https://www.kaggle.com/competitions/g2net-detecting-continuous-gravitational-waves/discussion/347052#1986395):\n\n1. Use `Writer` to generate noise-only SFTs, using Band and F0 to specify a broad frequency band [F0 - Band/2, F0 + Band/2]\n2. Use another `Writer` to inject a signal in the previous SFTs via the `noiseSFTs` option. This will return the same band as the original SFTs.\n3. Read into numpy array and select the band of interest\n\nBelow there is an example to generate a signal within the [150., 150.2) Hz band:","metadata":{}},{"cell_type":"code","source":"# Generate SFTs noise-only SFTs covering the band of interest\nnoise_kwargs = {\n    \"tstart\": 1238166018,\n    \"duration\": 4 * 30 * 86400,\n    \"sqrtSX\": 1e-23,\n    \"detectors\": \"H1,L1\",\n    \"Tsft\": 1800,\n    \"F0\": 150.1, # No signals: [F0 - Band/2, F0 + Band/2]\n    \"Band\": 0.5, \n    \"SFTWindowType\": \"tukey\",\n    \"SFTWindowBeta\": 0.001,\n}\n\nnoise_writer = pyfstat.Writer(label=\"custom_band_noise\", **noise_kwargs)\nnoise_writer.make_data()\n\n# Inject signal into noise SFTs. Note the lack of `Band` argument.\nsignal_kwargs = {\n        \"noiseSFTs\": noise_writer.sftfilepath,\n        \"F0\": 150.15,\n        \"F1\": -1e-8,\n        \"Alpha\": 0.3,\n        \"Delta\": 0,\n        \"h0\": 1e-23/10,\n        \"cosi\": 1,\n        \"psi\": 0.2,\n        \"phi\": 0.\n        }\nfor key in [\"SFTWindowType\", \"SFTWindowBeta\"]:\n    signal_kwargs[key] = noise_kwargs[key]\n\nsignal_writer = pyfstat.Writer(label=\"custom_band_signal\", **signal_kwargs)\nsignal_writer.make_data()\n\n# Slice out the band of interest\nfreqs, times, sft_data = pyfstat.utils.get_sft_as_arrays(signal_writer.sftfilepath)\n\nfirst_index = np.argmin(np.abs(freqs - 150.))\nlast_index = np.argmin(np.abs(freqs - 150.2))\n\nfreqs = freqs[first_index:last_index+1]\namplitudes = {key: val[first_index:last_index + 1, :]\n        for key, val in sft_data.items()}\n\nfreqs","metadata":{"execution":{"iopub.status.busy":"2022-10-15T07:33:01.273801Z","iopub.execute_input":"2022-10-15T07:33:01.274222Z","iopub.status.idle":"2022-10-15T07:33:11.107527Z","shell.execute_reply.started":"2022-10-15T07:33:01.274187Z","shell.execute_reply":"2022-10-15T07:33:11.106493Z"},"trusted":true},"execution_count":null,"outputs":[]}]}