From c799b61acf458a3df1f5a7fa5fea869968e60ee2 Mon Sep 17 00:00:00 2001 From: Claudio Galelli Date: Fri, 8 Nov 2024 14:43:52 +0100 Subject: [PATCH] Delete old recipes/fit-psd-lightcurve/fit-psd-lightcurve.ipynb --- .../fit-psd-lightcurve.ipynb | 561 ------------------ 1 file changed, 561 deletions(-) delete mode 100644 recipes/fit-psd-lightcurve/fit-psd-lightcurve.ipynb diff --git a/recipes/fit-psd-lightcurve/fit-psd-lightcurve.ipynb b/recipes/fit-psd-lightcurve/fit-psd-lightcurve.ipynb deleted file mode 100644 index 4d5be98..0000000 --- a/recipes/fit-psd-lightcurve/fit-psd-lightcurve.ipynb +++ /dev/null @@ -1,561 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "b4ca33d6-129c-4a59-89d7-f57e39cacfc0", - "metadata": {}, - "source": [ - "# Generating synthetic lightcurves and fitting the power spectral density of a lightcurve #\n", - "\n", - "This notebook presents the advanced Emmanoulopoulos algorithm for the simulation of synthetic lightcurves. The original paper describing the algorithm is linked [here](https://arxiv.org/pdf/1305.0304.pdf). The version implemented here is compatible with the Gammapy implementation of the Timmer-Koenig algorithm.\n", - "The Timmer-Koenig algorithm generates synthetic lightcurve from a chosen power spectral density (PSD) shape. However, it can only generate time series with a gaussian probability density function (PDF). This is adequate for high-statistics astrophysical domains such as the optical or X-rays, but can be in issue when trying to reproduce curves in the gamma-ray domain, where photon counts are lower and statistics are generally Poissonian. The Emmanoulopoulos algorithm tries to solve this issue, combining a requested PSD and PDF in the simulation. It provides accurate synthetic lightcurves in a range of spectral indexes between -1 and -2 for power-law or similar PSDs.\n", - "\n", - "Together with the simulation algorithm the notebook shows a function to compute the PSD envelope for a lightcurve using either the Timmer-Koenig or the Emmanoulopoulos algorithm. This envelope is then used to fit the PSD fot he observed lightcurve, by passing through a tailored chi-squared-like cost function. This complex fitting is necessary to account for the fact that the periodogram of the observed lightcurve is only a possible realization of the PSD model, moreover convoluted with Poissonian noise and instrumental responses. This can lead to biases or deformation due to random fluctuation of the realization if extracted with a simple curve fit of the periodogram.\n", - "\n", - "The results are satisfactory for power-law or broken-power-law PSDs in a physical interval of spectral indexes, between -1 and -2. Using the Emmanoulopoulos algorithm shows slightly better PSD reconstruction over the Timmer-Koenig after folding of the instrumental effects, at the cost of a longer computation time. In general, using the Timmer-Koenig algorithm is enough if the user is interested only in the frequency domain, while using the Emmanoulopoulos becomes indespensable to model correctly also the time domain.\n", - "\n", - "The functions for the Timmer-Koenig and Emmanoulopoulos algorithms, the envelope, x2-like cost function for fitting, and some helper analytical functions are implemented in the helper package `gammapy_SyLC`." - ] - }, - { - "cell_type": "markdown", - "id": "b693b292-fa5c-4afa-b013-464389e4092a", - "metadata": {}, - "source": [ - "## Imports ##\n", - "\n", - "The first step is importing some usual packages, needed Astropy utilities, scipy tools for PDFs and minimization, and Gammapy functions and classes for the observational part." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a0819ee7", - "metadata": { - "execution": { - "iopub.execute_input": "2024-08-30T15:54:42.755926Z", - "iopub.status.busy": "2024-08-30T15:54:42.755771Z", - "iopub.status.idle": "2024-08-30T15:54:42.762731Z", - "shell.execute_reply": "2024-08-30T15:54:42.762255Z" - } - }, - "outputs": [], - "source": [ - "import gammapy\n", - "\n", - "print(f\"Gammapy version : {gammapy.__version__}\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "4678c090-734e-430d-a81e-ff77cf5ef8a0", - "metadata": { - "execution": { - "iopub.execute_input": "2024-08-30T15:54:42.764787Z", - "iopub.status.busy": "2024-08-30T15:54:42.764421Z", - "iopub.status.idle": "2024-08-30T15:54:44.001370Z", - "shell.execute_reply": "2024-08-30T15:54:44.000908Z" - } - }, - "outputs": [], - "source": [ - "import numpy as np\n", - "from pathlib import Path\n", - "import matplotlib.pyplot as plt\n", - "import inspect\n", - "\n", - "import astropy.units as u\n", - "from astropy.coordinates import SkyCoord, AltAz\n", - "\n", - "from regions import PointSkyRegion\n", - "\n", - "from gammapy.estimators import LightCurveEstimator, FluxPoints\n", - "from gammapy.makers import SpectrumDatasetMaker\n", - "from gammapy.data import Observation, observatory_locations, FixedPointingInfo\n", - "from gammapy.datasets import Datasets, SpectrumDataset\n", - "from gammapy.irf import load_irf_dict_from_file\n", - "from gammapy.maps import MapAxis, RegionGeom, TimeMapAxis, RegionNDMap\n", - "from gammapy.modeling.models import (\n", - " SkyModel,\n", - " PowerLawSpectralModel,\n", - " LightCurveTemplateTemporalModel,\n", - ")\n", - "from gammapy.estimators.utils import compute_lightcurve_fvar\n", - "from gammapy.utils.random import get_random_state\n", - "\n", - "from scipy.optimize import minimize\n", - "from scipy.signal import periodogram\n", - "from scipy.stats import lognorm" - ] - }, - { - "cell_type": "markdown", - "id": "3a7b5a24", - "metadata": {}, - "source": [ - "Additionally, we import the necessary functions for lightcurve simulation and fitting from the helper package [gammapy_SyLC](https://github.com/cgalelli/gammapy_SyLC). The package can be obtained from github and installed by running:\n", - "\n", - "`git clone https://github.com/cgalelli/gammapy_SyLC.git` \n", - "`cd gammapy_SyLC` \n", - "`python -m pip install -e .` " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3d340431", - "metadata": { - "execution": { - "iopub.execute_input": "2024-08-30T15:54:44.003851Z", - "iopub.status.busy": "2024-08-30T15:54:44.003449Z", - "iopub.status.idle": "2024-08-30T15:54:44.006653Z", - "shell.execute_reply": "2024-08-30T15:54:44.006227Z" - } - }, - "outputs": [], - "source": [ - "from gammapy_SyLC import (\n", - " TimmerKonig_lightcurve_simulator,\n", - " Emmanoulopoulos_lightcurve_simulator,\n", - " lightcurve_psd_envelope,\n", - " minimize_x2_fit,\n", - " pl,\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "6ff44b6c-f786-4f7d-bf7c-c27cac5a78dc", - "metadata": {}, - "source": [ - "## Reference Lightcurve ##\n", - "\n", - "As a reference, the notebook uses the H.E.S.S. dataset for the PKS2155 AGN flare of 2006. Data properties such as mean and standard deviation fo the norm, number of points, sampling frequency, are taken from this flare. The synthetic lightcurve will be oversampled by a factor 10." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "330ee851-fcc3-47ad-8b1a-e5bd06d78869", - "metadata": { - "execution": { - "iopub.execute_input": "2024-08-30T15:54:44.008704Z", - "iopub.status.busy": "2024-08-30T15:54:44.008242Z", - "iopub.status.idle": "2024-08-30T15:54:44.283846Z", - "shell.execute_reply": "2024-08-30T15:54:44.283284Z" - } - }, - "outputs": [], - "source": [ - "lc_path = Path(\"$GAMMAPY_DATA/estimators/\")\n", - "lc_filename = \"pks2155_hess_lc/pks2155_hess_lc.fits\"\n", - "\n", - "lc = FluxPoints.read(lc_path / lc_filename, format=\"lightcurve\")\n", - "odata = lc.norm.data.flatten()\n", - "omean = odata.mean()\n", - "ostd = odata.std()\n", - "npoints = len(lc.norm.data) * 10\n", - "times = lc.geom.axes[\"time\"].edges\n", - "tref = lc.geom.axes[\"time\"].reference_time\n", - "smax = np.diff(times).max() / 10\n", - "lc.plot()\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "69638a7d", - "metadata": {}, - "source": [ - "## Simulation ##\n", - "\n", - "\n", - "As a first step, the parameters for the functions used as models in the simulations are setup here:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8b70ce94-ba71-43ef-80be-70b15e1a1d6f", - "metadata": { - "execution": { - "iopub.execute_input": "2024-08-30T15:54:44.285948Z", - "iopub.status.busy": "2024-08-30T15:54:44.285610Z", - "iopub.status.idle": "2024-08-30T15:54:44.288140Z", - "shell.execute_reply": "2024-08-30T15:54:44.287745Z" - } - }, - "outputs": [], - "source": [ - "ln_params = {\"s\": 0.5, \"loc\": 1.5, \"scale\": 1}\n", - "pl_params = {\"index\": -1.4}" - ] - }, - { - "cell_type": "markdown", - "id": "78dd6fa2-c725-42b8-a854-5eefaaec5be2", - "metadata": {}, - "source": [ - "\n", - "Both the TK and EMM algorithms are called with the same power-law PSD. The EMM algorithm uses a lognormal PDF. The difference between TK and EMM algorithms is shown in the leftmost and rightmost plot, where the gaussian vs lognormal shape is evident. The middle plot shows the perfect compatibility in the periodogram. Seed is fixed for reproducibility. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9fa4aa20-d5e4-4352-8a7c-6348d5c882cc", - "metadata": { - "execution": { - "iopub.execute_input": "2024-08-30T15:54:44.289970Z", - "iopub.status.busy": "2024-08-30T15:54:44.289808Z", - "iopub.status.idle": "2024-08-30T15:54:44.811522Z", - "shell.execute_reply": "2024-08-30T15:54:44.810985Z" - } - }, - "outputs": [], - "source": [ - "%%time\n", - "seed = 532019\n", - "# seed = \"random-seed\"\n", - "lctk2, taxis2 = Emmanoulopoulos_lightcurve_simulator(\n", - " lognorm.pdf,\n", - " pl,\n", - " npoints,\n", - " smax,\n", - " pdf_params=ln_params,\n", - " psd_params=pl_params,\n", - " mean=omean,\n", - " std=ostd,\n", - " random_state=seed,\n", - ")\n", - "lctk, taxis = TimmerKonig_lightcurve_simulator(\n", - " pl,\n", - " npoints,\n", - " smax,\n", - " power_spectrum_params=pl_params,\n", - " mean=omean,\n", - " std=ostd,\n", - " random_state=seed,\n", - ")\n", - "freqstk, pgramtk = periodogram(lctk, 1 / smax.value)\n", - "freqstk2, pgramtk2 = periodogram(lctk2, 1 / smax.value)\n", - "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 3))\n", - "ax1.plot(taxis, lctk)\n", - "ax2.loglog(freqstk[1:], pgramtk[1:])\n", - "ax3.hist(lctk)\n", - "ax1.plot(taxis2, lctk2)\n", - "ax2.loglog(freqstk2[1:], pgramtk2[1:])\n", - "ax3.hist(lctk2)\n", - "coeff = np.polyfit(np.log(freqstk[1:]), np.log(pgramtk[1:]), 1)\n", - "coeff2 = np.polyfit(np.log(freqstk2[1:]), np.log(pgramtk2[1:]), 1)\n", - "\n", - "print(coeff, coeff2)" - ] - }, - { - "cell_type": "markdown", - "id": "592cd9e0-55d1-4bbf-93f4-e873b011d43d", - "metadata": {}, - "source": [ - "## Gammapy setup and simulation ##\n", - "\n", - "Setup of geometry for the Gammapy simulation. Generic setup for pointing, energy binning, and IRFs. For realistic simulations, choose IRFs that are consistent with the instrument and observational conditions." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "79fc7316-39ef-490d-9cf9-d714afe9249e", - "metadata": { - "execution": { - "iopub.execute_input": "2024-08-30T15:54:44.813763Z", - "iopub.status.busy": "2024-08-30T15:54:44.813563Z", - "iopub.status.idle": "2024-08-30T15:54:44.889926Z", - "shell.execute_reply": "2024-08-30T15:54:44.889386Z" - } - }, - "outputs": [], - "source": [ - "TimeMapAxis.time_format = \"iso\"\n", - "\n", - "path = Path(\"$GAMMAPY_DATA/cta-caldb\")\n", - "irf_filename = \"Prod5-South-20deg-AverageAz-14MSTs37SSTs.180000s-v0.1.fits.gz\"\n", - "\n", - "irfs = load_irf_dict_from_file(path / irf_filename)\n", - "\n", - "energy_axis = MapAxis.from_energy_bounds(\n", - " energy_min=0.1 * u.TeV, energy_max=100 * u.TeV, nbin=1\n", - ")\n", - "\n", - "energy_axis_true = MapAxis.from_edges(\n", - " np.logspace(-1.2, 2.0, 31), unit=\"TeV\", name=\"energy_true\", interp=\"log\"\n", - ")\n", - "\n", - "time_axis = MapAxis.from_nodes(taxis, name=\"time\", interp=\"lin\")\n", - "\n", - "geom = RegionGeom.create(\n", - " \"galactic;circle(107.65, -40.17, 5)\", axes=[energy_axis]\n", - ")\n", - "\n", - "pointing_position = SkyCoord(107.65, -40.17, unit=\"deg\", frame=\"galactic\")\n", - "pointing = FixedPointingInfo(\n", - " fixed_icrs=SkyCoord(107.65, -40.17, unit=\"deg\", frame=\"galactic\").icrs,\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "8d9aec9b-9f38-4f9b-903e-e5189c7edf6c", - "metadata": {}, - "source": [ - "The time series generated via EMM is taken as a LightCurveTemplateTemporalModel" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c900eb8e-87ce-4771-8ccc-db627c7a541f", - "metadata": { - "execution": { - "iopub.execute_input": "2024-08-30T15:54:44.891959Z", - "iopub.status.busy": "2024-08-30T15:54:44.891781Z", - "iopub.status.idle": "2024-08-30T15:54:44.914387Z", - "shell.execute_reply": "2024-08-30T15:54:44.913901Z" - } - }, - "outputs": [], - "source": [ - "gti_t0 = tref\n", - "\n", - "spectral_model = PowerLawSpectralModel(\n", - " amplitude=1e-10 * u.TeV**-1 * u.cm**-2 * u.s**-1\n", - ")\n", - "\n", - "m = RegionNDMap.create(\n", - " region=PointSkyRegion(center=pointing_position),\n", - " axes=[time_axis],\n", - " unit=\"cm-2s-1TeV-1\",\n", - ")\n", - "\n", - "m.quantity = lctk2\n", - "\n", - "temporal_model = LightCurveTemplateTemporalModel(m, t_ref=gti_t0)\n", - "\n", - "model_simu = SkyModel(\n", - " spectral_model=spectral_model,\n", - " temporal_model=temporal_model,\n", - " name=\"model-simu\",\n", - ")" - ] - }, - { - "cell_type": "markdown", - "id": "7ccc6acb-5054-4879-bba9-ae0bd14e6b02", - "metadata": {}, - "source": [ - "Observation timing setup and simulation fo the datasets. The \"observational\" sampling is taken to be much sparser than the synthetic lightcurve, to avoid aliasing." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "bd889da2-54be-461f-b9a4-a34959e994a7", - "metadata": { - "execution": { - "iopub.execute_input": "2024-08-30T15:54:44.916270Z", - "iopub.status.busy": "2024-08-30T15:54:44.915958Z", - "iopub.status.idle": "2024-08-30T15:54:52.495344Z", - "shell.execute_reply": "2024-08-30T15:54:52.494850Z" - } - }, - "outputs": [], - "source": [ - "lvtm = 10 * u.min\n", - "tstart = gti_t0 + np.arange(npoints / 10) * lvtm\n", - "altaz = pointing_position.transform_to(\n", - " AltAz(obstime=tstart, location=observatory_locations[\"cta_south\"])\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f142057c-852f-4605-88d3-b27d63417803", - "metadata": { - "execution": { - "iopub.execute_input": "2024-08-30T15:54:52.497403Z", - "iopub.status.busy": "2024-08-30T15:54:52.497238Z", - "iopub.status.idle": "2024-08-30T15:55:01.493621Z", - "shell.execute_reply": "2024-08-30T15:55:01.493099Z" - } - }, - "outputs": [], - "source": [ - "datasets = Datasets()\n", - "\n", - "empty = SpectrumDataset.create(\n", - " geom=geom, energy_axis_true=energy_axis_true, name=\"empty\"\n", - ")\n", - "\n", - "maker = SpectrumDatasetMaker(selection=[\"exposure\", \"background\", \"edisp\"])\n", - "\n", - "for idx in range(len(tstart)):\n", - " obs = Observation.create(\n", - " pointing=pointing,\n", - " livetime=lvtm,\n", - " tstart=tstart[idx],\n", - " irfs=irfs,\n", - " reference_time=gti_t0,\n", - " obs_id=idx,\n", - " location=observatory_locations[\"cta_south\"],\n", - " )\n", - " empty_i = empty.copy(name=f\"dataset-{idx}\")\n", - " dataset = maker.run(empty_i, obs)\n", - " dataset.models = model_simu\n", - " dataset.fake()\n", - " datasets.append(dataset)\n", - "\n", - "\n", - "spectral_model = PowerLawSpectralModel(\n", - " amplitude=7e-11 * u.TeV**-1 * u.cm**-2 * u.s**-1\n", - ")\n", - "model_fit = SkyModel(spectral_model=spectral_model, name=\"model-fit\")\n", - "datasets.models = model_fit" - ] - }, - { - "cell_type": "markdown", - "id": "725f0466-21c4-4da2-a5da-21486ef0a43d", - "metadata": {}, - "source": [ - "Lightcurve estimator setup and run." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6150efce-4ced-4237-bf53-12fbd5fe9a30", - "metadata": { - "execution": { - "iopub.execute_input": "2024-08-30T15:55:01.495821Z", - "iopub.status.busy": "2024-08-30T15:55:01.495543Z", - "iopub.status.idle": "2024-08-30T15:55:05.450069Z", - "shell.execute_reply": "2024-08-30T15:55:05.449504Z" - } - }, - "outputs": [], - "source": [ - "lc_maker_1d = LightCurveEstimator(\n", - " energy_edges=[0.1, 100] * u.TeV,\n", - " source=\"model-fit\",\n", - " selection_optional=[\"ul\"],\n", - ")\n", - "\n", - "lc_1d = lc_maker_1d.run(datasets)\n", - "lc_1d.plot();" - ] - }, - { - "cell_type": "markdown", - "id": "b27cd811-a990-48c1-93af-baee1198432c", - "metadata": {}, - "source": [ - "Assessment of the properties of the \"observed\" lightcurve in the time and frequency domain." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e2ee2f7b-c8ee-4251-b276-75d4cd247b66", - "metadata": { - "execution": { - "iopub.execute_input": "2024-08-30T15:55:05.452193Z", - "iopub.status.busy": "2024-08-30T15:55:05.451790Z", - "iopub.status.idle": "2024-08-30T15:55:05.675600Z", - "shell.execute_reply": "2024-08-30T15:55:05.675149Z" - } - }, - "outputs": [], - "source": [ - "data = lc_1d.norm.data.flatten()\n", - "dmean = data.mean()\n", - "dstd = data.std()\n", - "dnpoints = len(data)\n", - "dtimes = lc_1d.geom.axes[\"time\"].edges\n", - "dsmax = np.diff(dtimes).max()\n", - "ffreqs, pgram = periodogram(data, 1 / dsmax.value)\n", - "coeff = np.polyfit(np.log(ffreqs[1:]), np.log(pgram[1:]), 1)\n", - "print(coeff[0])\n", - "plt.loglog(ffreqs[1:], pgram[1:])" - ] - }, - { - "cell_type": "markdown", - "id": "759e1ddd-98dc-4854-8ba7-c1e58facbb87", - "metadata": {}, - "source": [ - "### Fitting ###\n", - "\n", - "The minimize_x2_fit function calls the internal x2_fit as a cost function for the scipy minimizer, providing a fit of the spectral index for the \"observed\" lightcurve, in this case assuming a power-law PSD. Since here only the PSD is fit, and the fluxes are relatively high, the TK algorithm can be used to save time." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b08789fe-2a0a-4056-8b38-1aaea9946059", - "metadata": { "tags": [ - "nbsphinx-thumbnail" - ] - }, - "outputs": [], - "source": [ - "%%time\n", - "index, index_err = minimize_x2_fit(pgram, pl, {\"index\":-2}, dsmax, simulator=\"TK\", nsims=10000, mean=dmean, std=dstd,method=\"Powell\")\n", - "envelopes, freqs = lightcurve_psd_envelope(\n", - " pl,\n", - " dnpoints,\n", - " dsmax,\n", - " psd_params={\"index\": index},\n", - " simulator=\"TK\",\n", - " nsims=10000,\n", - " mean=dmean,\n", - " std=dstd,\n", - ")\n", - "plt.violinplot(envelopes, freqs, widths=np.diff(freqs).min(), showmedians=True)\n", - "plt.plot(freqs, pgram[1:], linewidth=0.7, marker=\"d\")\n", - "plt.yscale(\"log\")\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "67eae614", - "metadata": {}, - "source": [ - "Recipe by [Claudio Galelli](https://github.com/cgalelli/)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.11.9" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -}