{"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":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-10-17T17:33:08.515822Z","iopub.execute_input":"2022-10-17T17:33:08.516162Z","iopub.status.idle":"2022-10-17T17:33:09.852035Z","shell.execute_reply.started":"2022-10-17T17:33:08.516129Z","shell.execute_reply":"2022-10-17T17:33:09.847919Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The contents of this notebook were adapted from [the PyFstat tutorials](https://github.com/PyFstat/PyFstat/tree/master/examples/tutorials).","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import sys","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:33:09.853823Z","iopub.execute_input":"2022-10-17T17:33:09.854317Z","iopub.status.idle":"2022-10-17T17:33:09.860421Z","shell.execute_reply.started":"2022-10-17T17:33:09.854251Z","shell.execute_reply":"2022-10-17T17:33:09.859152Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install git+https://github.com/PyFstat/PyFstat@python37","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-17T17:33:09.862314Z","iopub.execute_input":"2022-10-17T17:33:09.862661Z","iopub.status.idle":"2022-10-17T17:33:23.91715Z","shell.execute_reply.started":"2022-10-17T17:33:09.86263Z","shell.execute_reply":"2022-10-17T17:33:23.916128Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from typing import TYPE_CHECKING, Iterable, Optional\nimport logging","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-10-17T17:33:23.919613Z","iopub.execute_input":"2022-10-17T17:33:23.9201Z","iopub.status.idle":"2022-10-17T17:33:23.926424Z","shell.execute_reply.started":"2022-10-17T17:33:23.920048Z","shell.execute_reply":"2022-10-17T17:33:23.925004Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Utils","metadata":{}},{"cell_type":"code","source":"\"\"\"\nUtils\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\n#plt.rcParams[\"font.family\"] = \"Arial\"\n#plt.rcParams[\"font.size\"] = 20\n\n\ndef plot_real_imag_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 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,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-17T17:33:23.928304Z","iopub.execute_input":"2022-10-17T17:33:23.92872Z","iopub.status.idle":"2022-10-17T17:33:23.953051Z","shell.execute_reply.started":"2022-10-17T17:33:23.928677Z","shell.execute_reply":"2022-10-17T17:33:23.952128Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def set_up_logger(\n    outdir: Optional[str] = None,\n    label: Optional[str] = \"pyfstat\",\n    log_level: str = \"INFO\",  # FIXME: Requires Python 3.8 Literal[\"CRITICAL\", \"ERROR\", \"WARNING\", \"INFO\", \"DEBUG\"] = \"INFO\",\n    streams: Optional[Iterable[\"io.TextIOWrapper\"]] = (sys.stdout,),\n    append: bool = True,\n) -> logging.Logger:\n#     \"\"\"Add file and stream handlers to the `pyfstat` logger.\n#     Handler names generated from ``streams`` and ``outdir, label``\n#     must be unique and no duplicated handler will be attached by\n#     this function.\n#     Parameters\n#     ----------\n#     outdir:\n#         Path to outdir directory. If ``None``, no file handler will be added.\n#     label:\n#         Label for the file output handler, i.e.\n#         the log file will be called `label.log`.\n#         Required, in conjunction with ``outdir``, to add a file handler.\n#         Ignored otherwise.\n#     log_level:\n#         Level of logging. This level is imposed on the logger itself and\n#         *every single handler* attached to it.\n#     streams:\n#         Stream to which logging messages will be passed using a\n#         StreamHandler object. By default, log to ``sys.stdout``.\n#         Other common streams include e.g. ``sys.stderr``.\n#     append:\n#         If ``False``, removes all handlers from the `pyfstat` logger\n#         before adding new ones. This removal is not propagated to\n#         handlers on the `root` logger.\n#     Returns\n#     -------\n#     obj:\n#         Configured instance of the ``logging.Logger`` class.\n#     \"\"\"\n    logger = logging.getLogger(\"pyfstat\")\n    logger.setLevel(log_level)\n\n    if not append:\n        for handler in logger.handlers:\n            logger.removeHandler(handler)\n    else:\n        for handler in logger.handlers:\n            handler.setLevel(log_level)\n\n    stream_names = [\n        handler.stream.name\n        for handler in logger.handlers\n        if type(handler) == logging.StreamHandler\n    ]\n    file_names = [\n        handler.baseFilename\n        for handler in logger.handlers\n        if type(handler) == logging.FileHandler\n    ]\n\n    common_formatter = logging.Formatter(\n        \"%(asctime)s.%(msecs)03d %(name)s %(levelname)-8s: %(message)s\",\n        datefmt=\"%y-%m-%d %H:%M:%S\",  # intended to match LALSuite's format\n    )\n\n    for stream in streams or []:\n        if stream.name in stream_names:\n            continue\n        stream_handler = logging.StreamHandler(stream)\n        stream_handler.setFormatter(common_formatter)\n        stream_handler.setLevel(log_level)\n        logger.addHandler(stream_handler)\n\n    if label and outdir:\n        os.makedirs(outdir, exist_ok=True)\n        log_file = os.path.join(outdir, f\"{label}.log\")\n\n        if log_file not in file_names:\n\n            file_handler = logging.FileHandler(log_file)\n            file_handler.setFormatter(common_formatter)\n            file_handler.setLevel(log_level)\n            logger.addHandler(file_handler)\n\n    return logger","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-17T17:35:06.480122Z","iopub.execute_input":"2022-10-17T17:35:06.480633Z","iopub.status.idle":"2022-10-17T17:35:06.495583Z","shell.execute_reply.started":"2022-10-17T17:35:06.48058Z","shell.execute_reply":"2022-10-17T17:35:06.494458Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pyfstat\nfrom pyfstat.utils import get_sft_as_arrays\n\n# Local module to simplify plotting\n#import tutorial_utils\n\nlogger = set_up_logger(label=\"0_generating_noise\", log_level=\"INFO\")\n\n%matplotlib inline","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-17T17:35:06.49678Z","iopub.execute_input":"2022-10-17T17:35:06.497608Z","iopub.status.idle":"2022-10-17T17:35:06.51804Z","shell.execute_reply.started":"2022-10-17T17:35:06.497569Z","shell.execute_reply":"2022-10-17T17:35:06.516708Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The contents of this notebook were adapted from [the PyFstat tutorials](https://github.com/PyFstat/PyFstat/tree/master/examples/tutorials).","metadata":{}},{"cell_type":"markdown","source":"# Part 1 CW Noise","metadata":{}},{"cell_type":"markdown","source":"### Using pyfstat.Writer\nThe most basic example is to generate Gaussian noise as measured by a single detector. \nThis operation can be performed using pyfstat.Writer as follows:","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()\n","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:35:25.405878Z","iopub.execute_input":"2022-10-17T17:35:25.406341Z","iopub.status.idle":"2022-10-17T17:35:26.987912Z","shell.execute_reply.started":"2022-10-17T17:35:25.406302Z","shell.execute_reply":"2022-10-17T17:35:26.986368Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"SFT data is saved at the path specified in writer.sftfilepath. This binary format can be opened as Numpy arrays using pyfstat.helper.get_sft_as_arrays as follows:","metadata":{}},{"cell_type":"code","source":"# Read SFT data into numpy arrays and plot real and imaginary parts\nfrequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:35:26.991055Z","iopub.execute_input":"2022-10-17T17:35:26.991685Z","iopub.status.idle":"2022-10-17T17:35:27.016804Z","shell.execute_reply.started":"2022-10-17T17:35:26.991638Z","shell.execute_reply":"2022-10-17T17:35:27.015729Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"SFT data could contain different detectors which may be operating during different times. Thus, timestamps and fourier_data are dictionaries whose keys correspond to detector names (H1 for LIGO Hanford, L1 for LIGO Livingston). timestamps labels the time at which the data was taken using GPS seconds, while fourier_data contains the Fourier amplitude of such data segment. frequency, which is common for all the detectors, is a 1D array labeling the frequency bins.","metadata":{}},{"cell_type":"code","source":"plot_real_imag_spectrograms(\n    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:35:27.575546Z","iopub.execute_input":"2022-10-17T17:35:27.576513Z","iopub.status.idle":"2022-10-17T17:35:29.200205Z","shell.execute_reply.started":"2022-10-17T17:35:27.576449Z","shell.execute_reply":"2022-10-17T17:35:29.198973Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Since we are generating zero-mean white Gaussian noise, there is a simple relation between the single-sided amplitude spectral density $\\sqrt{S_{\\text{n}}}$  and the variance of the distribution given by\n$$\\sigma^2 = \\frac{1}{4} T_{\\text{SFT}} S_{\\text{n}}$$\n \nwhere the two 1/2 factors are due to the use of a single-sided ASD and the fact that this standard deviation applies to both the real and imaginary parts of the Fourier transform.","metadata":{}},{"cell_type":"code","source":"theoretical_stdev = np.sqrt(0.25 * writer_kwargs[\"Tsft\"]) * writer_kwargs[\"sqrtSX\"]\nplot_real_imag_histogram(fourier_data[\"H1\"], theoretical_stdev);","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:35:29.202151Z","iopub.execute_input":"2022-10-17T17:35:29.203174Z","iopub.status.idle":"2022-10-17T17:35:29.910636Z","shell.execute_reply.started":"2022-10-17T17:35:29.203126Z","shell.execute_reply":"2022-10-17T17:35:29.909349Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Non-stationary noise\nReal data, on the other hand, is hardly ever stationary over long periods of time. This is equivalent to having a time-varying amplitude spectral density sqrtSX, which can be easily implemented by running several instances of pyfstat.Writer and concatenating the resulting file paths using ;.","metadata":{}},{"cell_type":"code","source":"segment_lengths = [5 * 86400, 3 * 86400, 4 * 86400]\nsegment_sqrtSX = [4e-23, 1e-23, 3e-23]\n\nsft_path = []\n\n# Setup Writer\nwriter_kwargs = {\n    \"outdir\": \"PyFstat_example_data\",\n    \"tstart\": 1238166018,\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\",\n    \"SFTWindowBeta\": 0.01,\n}\n\nfor segment in range(len(segment_lengths)):\n    writer_kwargs[\"label\"] = f\"segment_{segment}\"\n    writer_kwargs[\"duration\"] = segment_lengths[segment]\n    writer_kwargs[\"sqrtSX\"] = segment_sqrtSX[segment]\n\n    if segment > 0:\n        writer_kwargs[\"tstart\"] += writer_kwargs[\"Tsft\"] + segment_lengths[segment - 1]\n\n    writer = pyfstat.Writer(**writer_kwargs)\n    writer.make_data()\n\n    sft_path.append(writer.sftfilepath)\n\nsft_path = \";\".join(sft_path)  # Concatenate different files using ;\nfrequency, timestamps, fourier_data = get_sft_as_arrays(sft_path)","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:35:29.912131Z","iopub.execute_input":"2022-10-17T17:35:29.912512Z","iopub.status.idle":"2022-10-17T17:35:30.001005Z","shell.execute_reply.started":"2022-10-17T17:35:29.912475Z","shell.execute_reply":"2022-10-17T17:35:30.000077Z"},"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-17T17:35:30.003212Z","iopub.execute_input":"2022-10-17T17:35:30.003618Z","iopub.status.idle":"2022-10-17T17:35:32.809847Z","shell.execute_reply.started":"2022-10-17T17:35:30.003581Z","shell.execute_reply":"2022-10-17T17:35:32.808705Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Gaps\nMoreover, interferometric detectors are not taking science-quality data throughout the full observing run, meaning a real datastream will contain \"gaps\" during which no data is present. These gaps can be either scheduled downtime periods to conduct maintenance in the detectors or the result of an environmental perturbation driving the detector away from its operating point.\n\nFunctionality to simulate gaps is implemented via the timestamps keyword in pyfstat.Writer. For each GPS timestamp, pyfstat.Writer will produce an SFT of data starting at such time and spanning Tsft seconds. Mind that real data can contain arbitrarily long / short gaps, meaning the time span between the end of an SFT and the beginning of the next one does not have to correspond to a multiple of Tsft.","metadata":{}},{"cell_type":"code","source":"timestamps = {\"H1\": 1238166018 + 1800 * np.array([0, 2, 4, 6])}\n\n# Setup Writer\nwriter_kwargs = {\n    \"label\": \"single_detector_gaps\",\n    \"outdir\": \"PyFstat_example_data\",\n    \"timestamps\": timestamps,\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()\nfrequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:35:32.811213Z","iopub.execute_input":"2022-10-17T17:35:32.8116Z","iopub.status.idle":"2022-10-17T17:35:34.32778Z","shell.execute_reply.started":"2022-10-17T17:35:32.811565Z","shell.execute_reply":"2022-10-17T17:35:34.326642Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Special attention should be payed to the plotting function on use, as some functions like pcolormesh may distort the plotting grid not to be squared.**","metadata":{}},{"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-17T17:35:34.329564Z","iopub.execute_input":"2022-10-17T17:35:34.329914Z","iopub.status.idle":"2022-10-17T17:35:35.076541Z","shell.execute_reply.started":"2022-10-17T17:35:34.329878Z","shell.execute_reply":"2022-10-17T17:35:35.075346Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Narrow instrumental artifacts\nAnother characteristic of real data that affects CW searches are persistent narrow instrumental artifacts, also known as lines, which appear as strong, monochromatic, features in the detector data. Some of these lines have a known origin (e.g. couplings to the power lines oscillating at 60 Hz, vibrational modes of the mirror suspensions, LEDs blinking) but others are poorly understood.\n\nThis kind of artifacts can be simulated using pyfstat.LineWriter, which allows to specify the frequency F0, initial phase phi and amplitude h0 of the narrow instrumental artifact.","metadata":{}},{"cell_type":"code","source":"writer_kwargs = {\n    \"label\": \"single_detector_spectral_line\",\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    \"phi\": 1.3,  # Initial phase of the spectral line\n    \"Band\": 1.0,  # Frequency band-width around F0 [Hz]                \"h0\": 1e-24,              # Amplitude of the spectral line\n    \"sqrtSX\": 1e-23,  # Single-sided Amplitude Spectral Density of the noise\n    \"Tsft\": 1800,  # Fourier transform time duration\n    \"SFTWindowType\": \"tukey\",\n    \"SFTWindowBeta\": 0.01,\n}\n\nwriter = pyfstat.LineWriter(**writer_kwargs)\nwriter.make_data()\n\nfrequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:35:35.080402Z","iopub.execute_input":"2022-10-17T17:35:35.081163Z","iopub.status.idle":"2022-10-17T17:35:35.123844Z","shell.execute_reply.started":"2022-10-17T17:35:35.081115Z","shell.execute_reply":"2022-10-17T17:35:35.122681Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_amplitude_phase_spectrograms(\n    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:35:35.125431Z","iopub.execute_input":"2022-10-17T17:35:35.125778Z","iopub.status.idle":"2022-10-17T17:35:35.939752Z","shell.execute_reply.started":"2022-10-17T17:35:35.125747Z","shell.execute_reply":"2022-10-17T17:35:35.938881Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Multiple detectors\nFinally, the generation of data for multiple detectors is supported by specifying the detector names as a comma-separated string in the detectors key word. The list of supported detectors includes, amongst others, the advanced LIGO detectors (H1 for LIGO Hanford, L1 for LIGO Livingston) and the advanced Virgo detector V1.\n\nsqrtSX can be specified as a single, common value for all detectors or as a comma-separated string containing one value for each of the specified detectors. Likewise, timestamps could be specified as a list of common timestamps for all detectors or as a dictionary of lists, using detector names as keys.","metadata":{}},{"cell_type":"code","source":"writer_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,L1,V1\",  # 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,1e-24,1e-25\",  # 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\nwriter.make_data()\n\nfrequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:35:35.941103Z","iopub.execute_input":"2022-10-17T17:35:35.941842Z","iopub.status.idle":"2022-10-17T17:35:37.731723Z","shell.execute_reply.started":"2022-10-17T17:35:35.941804Z","shell.execute_reply":"2022-10-17T17:35:37.73052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for ind, ifo in enumerate(timestamps.keys()):\n    theoretical_stdev = np.sqrt(0.25 * writer_kwargs[\"Tsft\"]) * float(\n        writer_kwargs[\"sqrtSX\"].split(\",\")[ind]\n    )\n\n    fig, ax = plot_real_imag_spectrograms(\n        timestamps[ifo], frequency, fourier_data[ifo]\n    )\n    fig.suptitle(ifo)\n\n    plot_real_imag_histogram(fourier_data[ifo], theoretical_stdev);","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:35:37.733702Z","iopub.execute_input":"2022-10-17T17:35:37.734496Z","iopub.status.idle":"2022-10-17T17:35:44.943614Z","shell.execute_reply.started":"2022-10-17T17:35:37.734453Z","shell.execute_reply":"2022-10-17T17:35:44.942362Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Part 2 Signals\nGenerating signals\n### Example on how to generate continuous gravitational-wave signals.","metadata":{}},{"cell_type":"code","source":"logger = pyfstat.set_up_logger(label=\"1_generating_signals\", log_level=\"INFO\")","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:35:44.944955Z","iopub.execute_input":"2022-10-17T17:35:44.94534Z","iopub.status.idle":"2022-10-17T17:35:44.952309Z","shell.execute_reply.started":"2022-10-17T17:35:44.945306Z","shell.execute_reply":"2022-10-17T17:35:44.950804Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Introduction\nContinuous gravitational-wave signals (CWs) are long-lasting forms of gravitational radiation. They are usually characterized using two sets of parameters, namely the amplitude parameters and the Doppler parameters, respectively referred to as $\\mathcal{A}$ and $\\lambda$ .\n\nFor the typical case of a rapidly-spinning neutron star, the amplitude parameter contain the nominal CW amplitude \nh0, the (cosine of the) source inclination angle with respect to the line of sight $\\cos\\iota$ , the polarization angle $\\psi$ and the initial phase of the wave $\\psi$. Depending on the emission mechanism, h0\n can be further described using further physical quantities such as the source's frequency, ellipticity, or distance to the detector.\n\nDoppler parameters describe the evolution of the gravitational-wave frequency both due to physical processes undergoing at the source and the motion of the interferometric detector with respect to the Solar system barycenter (i.e. the Sun). In this tutorial, we will limit ourselves to gravitational __wave frequency f0__\n, __spindown f1__\n and __sky position__  $\\hat{n}$ \n, which we will parametrize using the right ascension $\\alpha$ and declination  angles $\\delta$.\n\nA detailed explanation of these parameters can be found in this technical document.\nhttps://dcc.ligo.org/LIGO-T0900149/public\npyfstat.Writer allows these variables to be given as inputs to produce SFTs containing a CW signal and (optionally) Gaussian noise. Signal parameters are both available as attributes in the pyfstat.Writer instance as as a .cff file in outdir. Alternatively, as exemplified in PyFstat_example_injecting_into_noise_sfts.py, the noiseSFTs can be used to provide a pre-generated set of SFTs as background noise.\n> https://github.com/PyFstat/PyFstat/blob/ec86602bb2f93238492a7242ad90995f6654eab7/examples/other_examples","metadata":{}},{"cell_type":"code","source":"writer_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\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    \"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)","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:35:44.95411Z","iopub.execute_input":"2022-10-17T17:35:44.95457Z","iopub.status.idle":"2022-10-17T17:35:49.154812Z","shell.execute_reply.started":"2022-10-17T17:35:44.954535Z","shell.execute_reply":"2022-10-17T17:35:49.153444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_real_imag_spectrograms(\n    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n)\nfig, ax = plot_amplitude_phase_spectrograms(\n    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:35:49.156933Z","iopub.execute_input":"2022-10-17T17:35:49.157483Z","iopub.status.idle":"2022-10-17T17:35:58.829395Z","shell.execute_reply.started":"2022-10-17T17:35:49.15744Z","shell.execute_reply":"2022-10-17T17:35:58.82833Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"These plots show the typical features of CW signals.\n\nShort-scale amplitude modulations are due to the rotation of the Earth, as the detector has different sensitivities depending on its orientation.\n\nThe frequency is modulated in a quasi-periodic fashion due to the translation of the interferometric detector around the Sun. The downwards trend is due to the presence of a negative source spindown due to the emission of energy as gravitational waves.\n\nThe instantaneous frequency of a CW at the detector can be expressed as $$f(t) = \\left[f_0 + f_1 \\left(t - t_{\\text{ref}}\\right)\\right]\\left[1 + \\frac{\\vec{v} \\cdot \\hat{n}}{c}\\right],$$\nwhere $\\vec{v}/c$\n is the velocity of the detector normalized by the speed of light in vaccum and $t_{\\text{ref}}$\n is a fiducial reference time at which $f_0$\n and $f_1$\n are specified.\n\nDetector velocities can be retrieved from SFT files using pyfstat.DetectorStates:","metadata":{}},{"cell_type":"code","source":"states = pyfstat.DetectorStates().get_multi_detector_states_from_sfts(\n    writer.sftfilepath, central_frequency=writer.F0, time_offset=0\n)\n\nts = np.array([data.tGPS.gpsSeconds for data in states.data[0].data])\nvelocities = np.vstack([data.vDetector for data in states.data[0].data]).T\n\nn = np.array(\n    [\n        [\n            np.cos(writer.Alpha) * np.cos(writer.Delta),  # Cartesian X\n            np.sin(writer.Alpha) * np.cos(writer.Delta),  # Cartesian Y\n            np.sin(writer.Delta),  # Cartesian Z\n        ]\n    ]\n)\n\nf_inst = ((writer.F0 + (ts - writer.tref) * writer.F1) * (1 + np.dot(n, velocities)))[0]\n\nfor a in ax:\n    for off in [-2, 2]:\n        a.plot(\n            (ts - ts[0]) / 1800,\n            f_inst + off / writer.Tsft,\n            color=\"white\",\n            ls=\"--\",\n            label=\"Instantaneous frequency\" if off == 2 else \"\",\n        )\n    a.legend()\nfig","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:35:58.830752Z","iopub.execute_input":"2022-10-17T17:35:58.831088Z","iopub.status.idle":"2022-10-17T17:36:01.49661Z","shell.execute_reply.started":"2022-10-17T17:35:58.831057Z","shell.execute_reply":"2022-10-17T17:36:01.495522Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"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\n $$\\mathcal{D} = \\frac{\\sqrt{S_{\\text{n}}}}{h_0}\\;.$$\n\nCurrent searches for CW signals from unknown sources are able to detect a sensitivity depth ranging between $10\\;\\text{Hz}^{-1/2}$\n and $10\\;\\text{Hz}^{-1/2}$\n \n__Visual features, however, tend to disappear around a depth of__  $20\\;\\text{Hz}^{-1/2}$ \n","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)","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:36:01.502734Z","iopub.execute_input":"2022-10-17T17:36:01.503531Z","iopub.status.idle":"2022-10-17T17:36:05.718441Z","shell.execute_reply.started":"2022-10-17T17:36:01.503486Z","shell.execute_reply":"2022-10-17T17:36:05.717207Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_real_imag_spectrograms(\n    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n)\nplot_amplitude_phase_spectrograms(\n    timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:36:05.720414Z","iopub.execute_input":"2022-10-17T17:36:05.72081Z","iopub.status.idle":"2022-10-17T17:36:15.28859Z","shell.execute_reply.started":"2022-10-17T17:36:05.720771Z","shell.execute_reply":"2022-10-17T17:36:15.287583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Alternatively, CW signals can be characterized in terms of their (squared) signal-to-noise ratio (SNR) $\\rho^2$\n. As opposed to sensitivity depth, which focuses on instantaneous features of the CW, $\\rho^2$\n is an integrated quantity along the full duration of the observing run. The dependency of the optimal SNR (i.e. assuming the Doppler parameters match those of the signal) depends in a non-trivial way on the amplitude parameters and sky position of the source (due to the anisotropy of the detector response function) [Eq. (77)]; the average optimal SNR for a uniform distribution of sources across the sky with an isotropic polarization angle, however, can be readily expressed as\n $$\\langle \\rho^2 \\rangle_{\\vec{n}, \\psi} = \n\\frac{1}{20}\n\\frac{h_0^2 \\left(\\cos^4\\iota + 6 \\cos^2\\iota + 1\\right)}\n{S_{\\textrm{n}}}\nT_{\\textrm{obs}},$$\n \nwhere $T_{\\textrm{obs}}$\n represents full duration of the data stream.\n\nThe optimal SNR for a specific template can be computed using pyfstat.SignalToNoiseRatio assuming Gaussian noise with a given sqrtSX value or from a specific set of SFTs.","metadata":{}},{"cell_type":"code","source":"snr = pyfstat.SignalToNoiseRatio.from_sfts(F0=writer.F0, sftfilepath=writer.sftfilepath)\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)\nprint(f\"SNR: {np.sqrt(squared_snr)}\")","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:36:15.289974Z","iopub.execute_input":"2022-10-17T17:36:15.290573Z","iopub.status.idle":"2022-10-17T17:36:17.443931Z","shell.execute_reply.started":"2022-10-17T17:36:15.290537Z","shell.execute_reply":"2022-10-17T17:36:17.442532Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Generating an ensemble of signals\nCharacterizing a CW search or follow-up method usually involves analyzing its response to an ensemble of signals from a population of interest. For the case of all-sky searches, this population is usually consistent with a uniform distribution of sources across the sky with isotropic orientation.\n\nThe pyfstat.InjectionParametersGenerator provides methods to draw parameters from generic distributions in a suitable format, ready to be fed into pyfstat.Writer. Uniform sampling across the sky is baked in the child class pyfstat.AllSkyInjectionParametersGenerator, and isotropic priors on amplitude parameters are available in pyfstat.injection_parameters.isotropic_amplitude_priors.\n\nAll-sky searches tend to perform injection campaigns at fixed values of h0 (or, equivalently,$\\mathcal{D}$ ).","metadata":{}},{"cell_type":"code","source":"# Generate 10 signals with parameters drawn from a specific population\nnum_signals = 10\n\nwriter_kwargs = {\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\nsignal_parameters_generator = pyfstat.AllSkyInjectionParametersGenerator(\n    priors={\n        \"F0\": {\"uniform\": {\"low\": 100.0, \"high\": 100.1}},\n        \"F1\": -1e-10,\n        \"F2\": 0,\n        \"h0\": writer_kwargs[\"sqrtSX\"] / 10,  # Fix amplitude at depth 10.\n        **pyfstat.injection_parameters.isotropic_amplitude_priors,\n        \"tref\": writer_kwargs[\"tstart\"],\n    },\n)\n\nfor ind in range(num_signals):\n\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()","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:36:17.447361Z","iopub.execute_input":"2022-10-17T17:36:17.447788Z","iopub.status.idle":"2022-10-17T17:36:49.74375Z","shell.execute_reply.started":"2022-10-17T17:36:17.447752Z","shell.execute_reply":"2022-10-17T17:36:49.742399Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Part 3\n#### Questions and hypothesis\n1. Are the signals we are looking for located in the center of the frequency band? Not necessarily. You can assume that samples labeled with a 1 **fully contain a signal within the 0.2 Hz band**. Now, that signal may be in the center, it may be in the upper half, it may even cross the from top left corner to the bottom right corner, you name it.\n1.  How do I generate SFTs with signals not centered around [F0+Band/2,F0-Band/2]. The make_sfts doesn't seem to support the non-centered signals? The Writer class is thought to generate CW signals with a broad-enough frequency band so that our analysis codes down the line are able to estimate noise properties from the data (for example, it makes sure there are enough extra bins to compute running medians around the band of interest). I though one could bypass those extra bins easily **but there are no options to do so in PyFstat.**\n1.  troubles to read the .sft files written at the end of tutorial, Once the data is generated as SFT files, you should use the get_sft_as_arrays function (from pyfstat.utils) to read the SFTs you just created into Numpy arrays (see last line of the first cell of Tutorial 1).In your case, if your SFTs are in PyFstat_example_data_ensemble/Signal_3/H-5760_H1_1800SFT_Signal_3-1238166018-10368000.sft, you could read them as example below.\n1. I have a question about the fact that competitors need to generate additional data for training: why do we have to generate our own samples, i.e. why did not you generate all samples by yourself and then upload them as the full training data? Is it because you except each competitor to have its own training dataset and also except that each competitor tries to generate data using personal approach (e.g. based on some parameter estimations using the available training instances)? limited space.\n1. I am wondering about the SNR of the real data. What range of SNRs do physicist's actually expect?Is it right to assume that a SFT given out of blue might have high SNR but might not have the CW signal since signal can be from say electrical interference (60Hz)? **Since we haven't detected any of these signals yet, we have no idea about how strong they actually are. You can check some of the LIGO papers (look for the CW tag) to get a sense of how sensitive our searches are (i.e. how strong could this signals be while remaining undetected); I wouldn't quote any specific numbers, however, since they are highly dependent on search setup and detector configuration.** Minor comment on notation: when we talk about SNR, we refer to the SNR of a specific model. That is, given a set of parameters characterizing a signal, I can compute the SNR associated to that specific signal in the data. Data itself doesn't have an \"SNR\" in our language.\n1. Hi.Is there a fixed window length for fourier transform？**Yes, each timestamp labels a period of time of 1800s (you can check that by noting the frequency resolution of 1/1800 Hz), but timestamps need not to be consecutive (i.e. there may be times at which no data is collected).** Fourier transforms are taken over continuous segments. For every timestamp, you can assume there is a continuous span of time on which we applied the Fourier transform. Whenever we hit a gap, we stop and wait until we have again enough data to compute another Fourier transform.\n1. I am unclear about the \"target\" part while generating the data. The kenrel provided to generate the data does not seem to the information about the label (target 1, 0, or -1). Or am I missing something? Thanks in advance.**0/1 labels refer to whether we included a simulated signal to that sample (1) or if it consists only on noise.**","metadata":{}},{"cell_type":"markdown","source":"#### SFTs with signals not centered around \nWhat I would recommend is to generate broad-band SFTs and the slice out the frequency band of interest once they are read as a NumPy array. This can be done by 1) running Writer to generate noise-only SFTs (so that Band and F0 can be used to specify a frequency band and 2) running Writer using the previously generated SFTs as noiseSFTs without specifying\n\nFor example, suppose I want a signal within the [150., 150.2]Hz band.\n\nFirst, I'd create a set of noise SFTs covering that band","metadata":{}},{"cell_type":"code","source":"import pyfstat\nimport numpy as np\n\n# 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()","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:36:49.746118Z","iopub.execute_input":"2022-10-17T17:36:49.746695Z","iopub.status.idle":"2022-10-17T17:36:49.957999Z","shell.execute_reply.started":"2022-10-17T17:36:49.746639Z","shell.execute_reply":"2022-10-17T17:36:49.955868Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Then, I'd inject the signal","metadata":{}},{"cell_type":"code","source":"# for that signal to fit (including a few extra frequency bins!).\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()","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:36:49.960252Z","iopub.execute_input":"2022-10-17T17:36:49.961553Z","iopub.status.idle":"2022-10-17T17:36:50.354445Z","shell.execute_reply.started":"2022-10-17T17:36:49.961464Z","shell.execute_reply":"2022-10-17T17:36:50.353189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Finally, I read as a Numpy array and slice out the relevant frequency band:","metadata":{}},{"cell_type":"code","source":"# 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\nprint(\"*\" * 20)\nprint(\"*\" * 20)\nprint(f\"These should be 150. (got {freqs[0]}) \"\n      f\"and 150.2 (got {freqs[-1]}).\")","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:36:50.355968Z","iopub.execute_input":"2022-10-17T17:36:50.356343Z","iopub.status.idle":"2022-10-17T17:36:50.842805Z","shell.execute_reply.started":"2022-10-17T17:36:50.356308Z","shell.execute_reply":"2022-10-17T17:36:50.841547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Read generated Sft files\nAs 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!).**\nnoiseSFTs is used to add noise or a signal into a specific file of SFTs. In that case, one needs to know which window\nfunction was used, as failing to do so may significantly bias the SNR of a signal. You don't need this to read your data, butit's worth to keep it in mind what window you are using whenever you generate more.","metadata":{}},{"cell_type":"code","source":"from pyfstat.utils import get_sft_as_arrays","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:36:50.844646Z","iopub.execute_input":"2022-10-17T17:36:50.845126Z","iopub.status.idle":"2022-10-17T17:36:50.851057Z","shell.execute_reply.started":"2022-10-17T17:36:50.84508Z","shell.execute_reply":"2022-10-17T17:36:50.849741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"##### Load a saved file : signal 9 ","metadata":{}},{"cell_type":"code","source":"frequency, timestamps, sfts = get_sft_as_arrays(\"PyFstat_example_data_ensemble/Signal_9/H-17520_H1_1800SFT_Signal_9-1238166018-31536000.sft\")","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:36:50.852798Z","iopub.execute_input":"2022-10-17T17:36:50.853929Z","iopub.status.idle":"2022-10-17T17:36:51.375745Z","shell.execute_reply.started":"2022-10-17T17:36:50.853856Z","shell.execute_reply":"2022-10-17T17:36:51.374433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"timestamps","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:36:51.377766Z","iopub.execute_input":"2022-10-17T17:36:51.378143Z","iopub.status.idle":"2022-10-17T17:36:51.387451Z","shell.execute_reply.started":"2022-10-17T17:36:51.378109Z","shell.execute_reply":"2022-10-17T17:36:51.385618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"frequency","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:36:51.389087Z","iopub.execute_input":"2022-10-17T17:36:51.389653Z","iopub.status.idle":"2022-10-17T17:36:51.402972Z","shell.execute_reply.started":"2022-10-17T17:36:51.389607Z","shell.execute_reply":"2022-10-17T17:36:51.401772Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sfts","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:36:51.404435Z","iopub.execute_input":"2022-10-17T17:36:51.404925Z","iopub.status.idle":"2022-10-17T17:36:51.418668Z","shell.execute_reply.started":"2022-10-17T17:36:51.40489Z","shell.execute_reply":"2022-10-17T17:36:51.417059Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_real_imag_spectrograms(\n    timestamps[\"H1\"], frequency, sfts[\"H1\"]\n)\nplot_amplitude_phase_spectrograms(\n    timestamps[\"H1\"], frequency, sfts[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:36:51.42022Z","iopub.execute_input":"2022-10-17T17:36:51.420705Z","iopub.status.idle":"2022-10-17T17:36:59.64923Z","shell.execute_reply.started":"2022-10-17T17:36:51.420667Z","shell.execute_reply":"2022-10-17T17:36:59.648155Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"H-17520_H1_1800SFT_Signal_3-1238166018-31536000.sft","metadata":{}},{"cell_type":"markdown","source":"##### Load a saved file : signal 3 ","metadata":{}},{"cell_type":"code","source":"frequency, timestamps, sfts = get_sft_as_arrays(\"PyFstat_example_data_ensemble/Signal_3/H-17520_H1_1800SFT_Signal_3-1238166018-31536000.sft\")","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:36:59.650816Z","iopub.execute_input":"2022-10-17T17:36:59.651135Z","iopub.status.idle":"2022-10-17T17:37:00.230615Z","shell.execute_reply.started":"2022-10-17T17:36:59.651104Z","shell.execute_reply":"2022-10-17T17:37:00.229446Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_real_imag_spectrograms(\n    timestamps[\"H1\"], frequency, sfts[\"H1\"]\n)\nplot_amplitude_phase_spectrograms(\n    timestamps[\"H1\"], frequency, sfts[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:37:00.232504Z","iopub.execute_input":"2022-10-17T17:37:00.232827Z","iopub.status.idle":"2022-10-17T17:37:08.217553Z","shell.execute_reply.started":"2022-10-17T17:37:00.232797Z","shell.execute_reply":"2022-10-17T17:37:08.216544Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Part 3 PyFstat examples","metadata":{}},{"cell_type":"markdown","source":"### Compute a spectrogram","metadata":{}},{"cell_type":"code","source":"# not github-action compatible\n# plt.rcParams[\"font.family\"] = \"serif\"\n# plt.rcParams[\"font.size\"] = 18\n# plt.rcParams[\"text.usetex\"] = True\n\n# workaround deprecation warning\n# see https://github.com/matplotlib/matplotlib/issues/21723\nplt.rcParams[\"axes.grid\"] = False\n\nlabel = \"PyFstat_example_spectrogram\"\noutdir = os.path.join(\"PyFstat_example_data\", label)\nlogger = pyfstat.set_up_logger(label=label, outdir=outdir)\n\ndepth = 5\n\ndata_parameters = {\n    \"sqrtSX\": 1e-23,\n    \"tstart\": 1000000000,\n    \"duration\": 2 * 365 * 86400,\n    \"detectors\": \"H1\",\n    \"Tsft\": 1800,\n}\n\nsignal_parameters = {\n    \"F0\": 100.0,\n    \"F1\": 0,\n    \"F2\": 0,\n    \"Alpha\": 0.0,\n    \"Delta\": 0.5,\n    \"tp\": data_parameters[\"tstart\"],\n    \"asini\": 25.0,\n    \"period\": 50 * 86400,\n    \"tref\": data_parameters[\"tstart\"],\n    \"h0\": data_parameters[\"sqrtSX\"] / depth,\n    \"cosi\": 1.0,\n}\n\n# making data\ndata = pyfstat.BinaryModulatedWriter(\n    label=label, outdir=outdir, **data_parameters, **signal_parameters\n)\ndata.make_data()\n\nlogger.info(\"Loading SFT data and computing normalized power...\")\nfreqs, times, sft_data = pyfstat.utils.get_sft_as_arrays(data.sftfilepath)\nsft_power = sft_data[\"H1\"].real ** 2 + sft_data[\"H1\"].imag ** 2\nnormalized_power = (\n    2 * sft_power / (data_parameters[\"Tsft\"] * data_parameters[\"sqrtSX\"] ** 2)\n)\n\nplotfile = os.path.join(outdir, label + \".png\")\nlogger.info(f\"Plotting to file: {plotfile}\")\nfig, ax = plt.subplots(figsize=(0.8 * 16, 0.8 * 9))\nax.set(xlabel=\"Time [days]\", ylabel=\"Frequency [Hz]\", ylim=(99.98, 100.02))\nc = ax.pcolormesh(\n    (times[\"H1\"] - times[\"H1\"][0]) / 86400,\n    freqs,\n    normalized_power,\n    cmap=\"inferno_r\",\n    shading=\"nearest\",\n)\nfig.colorbar(c, label=\"Normalized Power\")\nplt.tight_layout()\nfig.savefig(plotfile)","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:37:49.431368Z","iopub.execute_input":"2022-10-17T17:37:49.432471Z","iopub.status.idle":"2022-10-17T17:38:07.121875Z","shell.execute_reply.started":"2022-10-17T17:37:49.432422Z","shell.execute_reply":"2022-10-17T17:38:07.1204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"def get_predict_fstat_parameters_from_dict(signal_parameters, transientWindowType=None):\n#     \"\"\"Extract a subset of parameters as needed for predicting F-stats.\n#     Given a dictionary with arbitrary signal parameters,\n#     this extracts only those ones required by `helper_functions.predict_fstat()`:\n#     Freq, Alpha, Delta, h0, cosi, psi.\n#     Also preserves transient parameters, if included in the input dict.\n#     Parameters\n#     ----------\n#     signal_parameters: dict\n#         Dictionary containing at least those signal parameters required by\n#         helper_functions.predict_fstat.\n#         This dictionary's keys must follow\n#         the PyFstat convention (e.g. F0 instead of Freq).\n#     transientWindowType: str\n#         Transient window type to store in the output dict.\n#         Currently required because the typical input dicts\n#         produced by various PyFstat functions\n#         tend not to store this property.\n#         If there is a key with this name already, its value will be overwritten.\n#     Returns\n#     -------\n#     predict_fstat_params: dict\n#         The dictionary of selected parameters.\n#     \"\"\"\n    required_keys = [\"F0\", \"Alpha\", \"Delta\", \"h0\", \"cosi\", \"psi\"]\n    transient_keys = {\n        \"transientWindowType\": \"transientWindowType\",\n        \"transient_tstart\": \"transientStartTime\",\n        \"transient_duration\": \"transientTau\",\n    }\n    predict_fstat_params = {key: signal_parameters[key] for key in required_keys}\n    for key in transient_keys:\n        if key in signal_parameters:\n            predict_fstat_params[transient_keys[key]] = signal_parameters[key]\n    if transientWindowType is not None:\n        predict_fstat_params[\"transientWindowType\"] = transientWindowType\n    return predict_fstat_params","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:41:36.703671Z","iopub.execute_input":"2022-10-17T17:41:36.704155Z","iopub.status.idle":"2022-10-17T17:41:36.714429Z","shell.execute_reply.started":"2022-10-17T17:41:36.704113Z","shell.execute_reply":"2022-10-17T17:41:36.712979Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"from pyfstat import (\n    AllSkyInjectionParametersGenerator,\n    InjectionParametersGenerator,\n    Writer,\n    set_up_logger,\n)\n\nlabel = \"PyFstat_example_InjectionParametersGenerator\"\noutdir = os.path.join(\"PyFstat_example_data\", label)\nlogger = set_up_logger(label=label, outdir=outdir)\n\n# Properties of the GW data\ngw_data = {\n    \"sqrtSX\": 1e-23,\n    \"tstart\": 1000000000,\n    \"duration\": 86400,\n    \"detectors\": \"H1,L1\",\n    \"Band\": 1,\n    \"Tsft\": 1800,\n}\n\nlogger.info(\"Drawing random signal parameters...\")\n\n# Draw random signal phase parameters.\n# The AllSkyInjectionParametersGenerator covers [Alpha,Delta] priors automatically.\n# The rest can be a mix of nontrivial prior distributions and fixed values.\nphase_params_generator = AllSkyInjectionParametersGenerator(\n    priors={\n        \"F0\": {\"uniform\": {\"low\": 29.0, \"high\": 31.0}},\n        \"F1\": -1e-10,\n        \"F2\": 0,\n    },\n    seed=23,\n)\nphase_parameters = phase_params_generator.draw()\nphase_parameters[\"tref\"] = gw_data[\"tstart\"]\n\n# Draw random signal amplitude parameters.\n# Here we use the plain InjectionParametersGenerator class.\namplitude_params_generator = InjectionParametersGenerator(\n    priors={\n        \"h0\": {\"normal\": {\"loc\": 1e-24, \"scale\": 1e-24}},\n        \"cosi\": {\"uniform\": {\"low\": 0.0, \"high\": 1.0}},\n        \"phi\": {\"uniform\": {\"low\": 0.0, \"high\": 2 * np.pi}},\n        \"psi\": {\"uniform\": {\"low\": 0.0, \"high\": np.pi}},\n    },\n    seed=42,\n)\namplitude_parameters = amplitude_params_generator.draw()\n\n# Now we can pass the parameter dictionaries to the Writer class and make SFTs.\ndata = Writer(\n    label=label,\n    outdir=outdir,\n    **gw_data,\n    **phase_parameters,\n    **amplitude_parameters,\n)\ndata.make_data()\n\n# Now we draw many phase parameters and check the sky distribution\nNdraws = 10000\nphase_parameters = [phase_params_generator.draw() for n in range(Ndraws)]\nAlphas = np.array([p[\"Alpha\"] for p in phase_parameters])\nDeltas = np.array([p[\"Delta\"] for p in phase_parameters])\nplotfile = os.path.join(outdir, label + \"_allsky.png\")\nlogger.info(f\"Plotting sky distribution of {Ndraws} points to file: {plotfile}\")\nplt.subplot(111, projection=\"aitoff\")\nplt.plot(Alphas - np.pi, Deltas, \".\", markersize=1)\nplt.savefig(plotfile, dpi=300)\nplt.close()\nplotfile = os.path.join(outdir, label + \"_alpha_hist.png\")\nlogger.info(f\"Plotting Alpha distribution of {Ndraws} points to file: {plotfile}\")\nplt.hist(Alphas, 50)\nplt.xlabel(\"Alpha\")\nplt.ylabel(\"draws\")\nplt.savefig(plotfile, dpi=100)\nplt.close()\nplotfile = os.path.join(outdir, label + \"_delta_hist.png\")\nlogger.info(f\"Plotting Delta distribution of {Ndraws} points to file: {plotfile}\")\nplt.hist(Deltas, 50)\nplt.xlabel(\"Delta\")\nplt.ylabel(\"draws\")\nplt.savefig(plotfile, dpi=100)\nplt.close()\nplotfile = os.path.join(outdir, label + \"_sindelta_hist.png\")\nlogger.info(f\"Plotting sin(Delta) distribution of {Ndraws} points to file: {plotfile}\")\nplt.hist(np.sin(Deltas), 50)\nplt.xlabel(\"sin(Delta)\")\nplt.ylabel(\"draws\")\nplt.savefig(plotfile, dpi=100)\nplt.close()\n","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:44:46.268223Z","iopub.execute_input":"2022-10-17T17:44:46.268648Z","iopub.status.idle":"2022-10-17T17:44:49.456061Z","shell.execute_reply.started":"2022-10-17T17:44:46.268615Z","shell.execute_reply":"2022-10-17T17:44:49.454143Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"### Compute the cumulative coherent F-statistic of a signal candidate.","metadata":{}},{"cell_type":"code","source":"label = \"PyFstat_example_twoF_cumulative\"\noutdir = os.path.join(\"PyFstat_example_data\", label)\nlogger = pyfstat.set_up_logger(label=label, outdir=outdir)\n\n# Properties of the GW data\ngw_data = {\n    \"sqrtSX\": 1e-23,\n    \"tstart\": 1000000000,\n    \"duration\": 100 * 86400,\n    \"detectors\": \"H1,L1\",\n    \"Band\": 4,\n    \"Tsft\": 1800,\n}\n\n# Properties of the signal\ndepth = 100\nphase_parameters = {\n    \"F0\": 30.0,\n    \"F1\": -1e-10,\n    \"F2\": 0,\n    \"Alpha\": np.radians(83.6292),\n    \"Delta\": np.radians(22.0144),\n    \"tref\": gw_data[\"tstart\"],\n    \"asini\": 10,\n    \"period\": 10 * 3600 * 24,\n    \"tp\": gw_data[\"tstart\"] + gw_data[\"duration\"] / 2.0,\n    \"ecc\": 0,\n    \"argp\": 0,\n}\namplitude_parameters = {\n    \"h0\": gw_data[\"sqrtSX\"] / depth,\n    \"cosi\": 1,\n    \"phi\": np.pi,\n    \"psi\": np.pi / 8,\n}\n\nPFS_input = get_predict_fstat_parameters_from_dict(\n    {**phase_parameters, **amplitude_parameters}\n)\n\n# Let me grab tref here, since it won't really be needed in phase_parameters\ntref = phase_parameters.pop(\"tref\")\ndata = pyfstat.BinaryModulatedWriter(\n    label=label,\n    outdir=outdir,\n    tref=tref,\n    **gw_data,\n    **phase_parameters,\n    **amplitude_parameters,\n)\ndata.make_data()\n\n# The predicted twoF, given by lalapps_predictFstat can be accessed by\ntwoF = data.predict_fstat()\nlogger.info(\"Predicted twoF value: {}\\n\".format(twoF))\n\n# Create a search object for each of the possible SFT combinations\n# (H1 only, L1 only, H1 + L1).\nifo_constraints = [\"L1\", \"H1\", None]\ncompute_fstat_per_ifo = [\n    pyfstat.ComputeFstat(\n        sftfilepattern=os.path.join(\n            data.outdir,\n            (f\"{ifo_constraint[0]}*.sft\" if ifo_constraint is not None else \"*.sft\"),\n        ),\n        tref=data.tref,\n        binary=phase_parameters.get(\"asini\", 0),\n        minCoverFreq=-0.5,\n        maxCoverFreq=-0.5,\n    )\n    for ifo_constraint in ifo_constraints\n]\n\nfor ind, compute_f_stat in enumerate(compute_fstat_per_ifo):\n    compute_f_stat.plot_twoF_cumulative(\n        label=label + (f\"_{ifo_constraints[ind]}\" if ind < 2 else \"_H1L1\"),\n        outdir=outdir,\n        savefig=True,\n        CFS_input=phase_parameters,\n        PFS_input=PFS_input,\n        custom_ax_kwargs={\n            \"title\": \"How does 2F accumulate over time?\",\n            \"label\": \"Cumulative 2F\"\n            + (f\" {ifo_constraints[ind]}\" if ind < 2 else \" H1 + L1\"),\n        },\n    )\n","metadata":{"execution":{"iopub.status.busy":"2022-10-17T17:41:41.69763Z","iopub.execute_input":"2022-10-17T17:41:41.698092Z","iopub.status.idle":"2022-10-17T17:43:49.808491Z","shell.execute_reply.started":"2022-10-17T17:41:41.698041Z","shell.execute_reply":"2022-10-17T17:43:49.807178Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}