Timing methods

This notebook demonstrates some of the “quick-look” timing analyses available in phoptic, as well as how to extend the available timing analyses using `stingray <https://docs.stingray.science/en/stable/>`__.

Generating and Reducing Data

Before we can explore the different timing analyses provided in phoptic, we need some light curves to analyse. For this example, I’ll generate and reduce some gappy observations (see the reduction tutorial for more details on data reduction):

import os

os.environ['OMP_NUM_THREADS'] = '1'

from pathlib import Path

import phoptic


out_dir = Path('out')

phoptic.generate_gappy_observations(
    out_directory=out_dir / 'data',  # path to the directory where the generated data will be saved
    circular_aperture=False,  # disable circular aperture shadow
    n_images=500,
    )

reducer = phoptic.Reducer(
    out_directory=out_dir / 'reduced',
    data_directory=out_dir / 'data',
)

reducer.create_catalogs()

phot = phoptic.OptimalPhotometer()
reducer.photometry(phot)

dphot = phoptic.DifferentialPhotometer(
    out_directory=out_dir / 'reduced',
)

unnormed_analyser = dphot.get_relative_light_curve(
    key='1:g',
    target=2,
    comparisons=[1, 3, 6],
    phot_label='optimal',
    prefix='example',
    match_other_cameras=True,
    show_diagnostics=False,
)

We now have an Analyzer instance that provides a number of quick-look timing analyses routines. Under-the-hood, these routines are mostly wrappers around astropy routines.

By default, light curves are not normalised when performing differential photometry. To normalise our light curves to a mean of 1, we can define a new Analyzer instance with norm='mean':

analyser = phoptic.Analyzer(
    out_directory=out_dir / 'reduced',
    light_curves=unnormed_analyser.light_curves,
    norm='mean',
)

phoptic.Analyzer routines

Let’s first take a look at our Analyzer light curves using the plot() method:

analyser.plot()

As we can see, we have three variable, gappy light curves, all with a mean flux of 1. The filter for each light curve is shown in the top right. We can also see that this figure was saved to timing_analysis_tutorial/reduced/plots/source_2_optimal_light_curves.pdf.

We can inspect our light curve data via the light_curves attribute. Within an Analyzer, light curves are stored as `astropy.timeseries.TimeSeries <https://docs.astropy.org/en/stable/timeseries/index.html>`__ instances:

We can also rebin our light curves to a lower time resolution using the rebin() method:

from astropy import units as u

binned_analyser = analyser.rebin(time_bin_size=10 * u.s)

This method takes a time_bin_size parameter, which must be a astropy.units.Quantity, and returns a new Analyzer instance containing the binned light curves. Under-the-hood, this uses astropy’s aggregate_downsample() function, which is slow for long timeseries.

Let’s take a look at the binned light curves:

binned_analyser.plot(save=False)

phoptic handles the error propagation during binning and avoids padding gaps in the original light curves with zeros. In this case, the bin size was too large, and the variability of the light curves has been lost. As such, we’ll carry on using the un-binned light curves.

Let’s now perform a few quick-look timing analyses on these light curves.

Lomb-Scargle Periodograms

Often times, you might want to compute the Lomb-Scargle periodograms of your light curves to check for periododic signals. astropy implements two flavours of the Lomb-Scargle periodogram: `LombScargle <https://docs.astropy.org/en/stable/timeseries/lombscargle.html>`__ and `LombScargleMultiband <https://docs.astropy.org/en/stable/timeseries/lombscarglemb.html>`__. For convenience, phoptic’s Analyzer class implements wrappers for both routines that handle the boilerplate, allowing you to get straight to the analysis.

Lomb-Scargle

Let’s first take a look at the standard Lomb-Scargle periodogram:

lsps = analyser.lomb_scargle(scale='semilogx')

We can see that each light curve has a dedicated periodogram, and the corresponding filter is shown in the top right of each panel. In this case, all three periodograms show a sharp peak between 0.1 and 0.2 Hz.

By default, the autofrequency() method of astropy’s LombScargle class is used to construct the frequency grid of the periodogram. Alternatively, you may pass your own frequency grid. To avoid confusion with units, custom frequency grids must be passed as an astropy `Quantity <https://docs.astropy.org/en/stable/units/quantity.html>`__. For example, let’s recompute these periodograms in the 0.1-0.2 Hz frequency range. However, we do not want to overwrite our previous plot, so we’ll also pass save=False to prevent the new plot from being saved:

from astropy import units as u
import numpy as np

frequency = np.linspace(0.1, 0.2, 100) * u.Hz  # units of Hz

lsps = analyser.lomb_scargle(frequency=frequency, save=False, scale='semilogx')

Now we can see that all three bands suggest a signal with a frequency of between 0.12 and 0.14 Hz.

Like many of the timing analysis methods of Analyzer, lomb_scargle() returns a dictionary of results that may be used to conduct further analyses. In this case, the keys list the filters and the values list the corresponding `astropy.timeseries.LombScargle <https://docs.astropy.org/en/stable/api/astropy.timeseries.LombScargle.html#astropy.timeseries.LombScargle.from_timeseries>`__ instances:

lsps

Let’s use the results returned by the lomb_scargle() method to precisely infer the peak frequencies:

import numpy as np

for fltr, lsp in lsps.items():
    power = lsp.power(frequency=frequency)
    peak_index = np.argmax(power)
    print(f'{fltr} Peak frequency = {frequency[peak_index]}')

Now we can see that all three bands show a periodicity at \(\sim\) 0.135 Hz which, as we can see at the top of this notebook, is the true frequency of variability for our source of interest.

Multiband Lomb-Scargle

In some cases, it may be preferable to compute a single Lomb-Scargle periodogram for all light curves; this can be done using astropy’s `LombScargleMultiband <https://docs.astropy.org/en/stable/timeseries/lombscarglemb.html>`__:

periodograms = analyser.multiband_lomb_scargle(scale='semilogx')

As with lomb_scargle(), custom frequency grids can also be passed to multiband_lomb_scargle():

periodograms = analyser.multiband_lomb_scargle(frequency=frequency, save=False, scale='semilogx')

Phase Folding/Binning

Now that we have identified a candidate signal, we may want fold our light curves to see what it looks like. astropy’s TimeSeries class implements a fold() method, and phoptic’s Analyzer class provides a convenience wrapper for this method with some added functionality (binning). To fold a light curve, we need a period on which to fold; to avoid confusion with units, phoptic requires that the period being passed is a Quantity:

f_signal = 0.135 * u.Hz  # periodogram peak frequency
p_signal = (1 / f_signal).to(u.s)  # fold period in seconds

folded_lcs = analyser.fold(
    period=p_signal,
    save=False,
)

The fold() method of phoptic.Analyzer returns a TimeSeries whose "time" column quantifies the folded phase:

folded_lcs.pprint()

Additionally, the fold() method can be used to phase bin your light curves by passing the desired number of bins to to the nbins parameter:

folded_lcs = analyser.fold(
    period=p_signal,
    nbins=10,
    save=False,
    )

Phase binning more clearly reveals the pulse shape compared to simply folding. From the above plots, it’s clear that there is a lag in the pulsations between the different bands. Currently, phoptic does not support quantify lags natively. However, phoptic’s Analyzer class can be used to export your light curves to `stingray <https://docs.stingray.science/en/stable/>`__, which provides a much more complete timing analysis framework - including ways to compute lags. Let’s take a look at this now.

stingray

stingray is a feature-rich spectral timing Python package developed primarily for use with X-ray data. stingray implements a number of generic timing analyis routines that can be applied to arbitrary time series. However, since stingray was developed with X-ray data in mind, it makes certain assumptions about the data. These assumptions can cause issues for non-X-ray data, as we’ll see soon.

To export your light curves to stingray, call the export_light_curves_to_stingray() method of phoptic.Analyzer:

stingray_lcs = analyser.export_light_curves_to_stingray()

We can see that a warning about error distributions has been triggered. We’ll take a look at this in more detail later. For now, let’s ignore it.

export_light_curves_to_stingray() returns a dictionary of {filter: `stingray.Lightcurve <https://docs.stingray.science/en/stable/notebooks/Lightcurve/Lightcurve%20tutorial.html>`__} pairs:

stingray_lcs

Let’s take a closer look at one of these light curves using the plot() method of stingray.Lightcurve:

ax = stingray_lcs['1:g'].plot()

As we can see, the light curve times are converted into seconds, with \(t=0\) corresponding to the t_ref attribute of the Analyzer instance. Additionally, phoptic has automatically inferred GTIs for each light curve, which are required for compatibility with many of stingray’s methods. The shaded regions in the above plot represent the time between GTIs. Let’s use these stingray.Lightcurve instances to do some more advanced analyses.

Lags

Cross-correlation

Having recognised a lag between our light curves, there are several ways we can quantify it. For example, we can use the cross-correlations between light curves. This is a relatively simple approach, but it has the benefit of working with arbitrary time series. We can compute cross-correlations in stingray using `CrossCorrelation <https://docs.stingray.science/en/stable/notebooks/CrossCorrelation/cross_correlation_notebook.html>`__:

from stingray.crosscorrelation import CrossCorrelation

cross_corr = CrossCorrelation(
    lc1=stingray_lcs['1:g'],
    lc2=stingray_lcs['3:i'],
    )

ax = cross_corr.plot(
    labels=['Lag [s]', 'Correlation'],
    )

We can see that the cross-correlation spectrum peaks at a non-zero lag. We can get the exact time lag from the time_shift attribute of the CrossCorrelation instance:

print(f'Time lag: {cross_corr.time_shift:.2f}s')

We find a negative lag of about 4 s, meaning \(i\) is lagging \(g\) (as expected).

When we generated this data, we were told that the \(i\)-band variability has a phase lag of 3.142 (i.e., \(1 \pi\)) radians, while the \(g\)-band variability has a phase lag of 0 radians. We can convert this time lag into a phase lag by dividing it by the period of the variability:

print(f'Phase lag: {2 * np.pi * cross_corr.time_shift / p_signal.to_value(u.s):.2f} radians')

The phase lag is a little higher than expected, but this is due to the time resolution of our light curves.

Cross-correlating generally works well, provided the light curves are dominated by the signal of interest. In many cases, however, light curves show variability on many time-scales, and we may want to isolate the lags of signals at specific frequencies. In such cases, it may be preferable to compute the cross-spectrum.

Cross-spectra

A Cautionary Tale

To compute cross-spectra, we can use stingray’s AveragedCrossspectrum class:

from stingray import AveragedCrossspectrum
from matplotlib import pyplot as plt

segment_size = 30 * u.s

avg_cs = AveragedCrossspectrum.from_lightcurve(
    lc1=stingray_lcs['1:g'],
    lc2=stingray_lcs['3:i'],
    norm='frac',
    segment_size=segment_size.to_value(u.s),  # stingray assumes all time values are in seconds
)

avg_cs_amplitude = np.abs(avg_cs.power)

fig, ax = plt.subplots(tight_layout=True)

ax.step(
    avg_cs.freq,  # convert freq from cyc/d to Hz
    avg_cs_amplitude,
    where='mid',
    c='k',
)

ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Power [(frac. rms)$^2$ / Hz]')

plt.show()

As we can see, the cross spectrum amplitude looks similar to the Lomb-Scargle periodograms we saw earlier, which is what we would expect. So far so good. Let’s try to compute the coherence between these two light curves and check it’s in the expected range (\(0 \leq \text{coherence} \leq 1\)):

coh, coh_err = avg_cs.coherence()

print(np.logical_and(coh >= 0, coh <= 1).all())

Uh oh, we have unphysical coherence values! Let’s take a look at some of them:

for i in range(5):
    print(f'Coherence: {coh[i]}, error: {coh_err[i]}')

Our coherence estimates are in the thousands, and our error estimates in the negative thousands. Something is very wrong…

The issue is due to stingray assuming Poisson statistics. Since we know the phase lag between these two light curves, let’s see what happens if we try to compute it from the averaged cross spectrum:

plags, plags_err = avg_cs.phase_lag()

fig, ax = plt.subplots(tight_layout=True)

ax.errorbar(
    avg_cs.freq,
    plags,
    plags_err,
    fmt='kx',
)

ax.axvline(
    0.135,
    c='red',
    label='Lag frequency'
    )
ax.axhline(
    np.pi,
    c='blue',
    label='Expected lag'
    )
ax.legend()

ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Phase lag [rad]')

plt.show()

We get the correct phase lag at the correct frequency. However, we don’t seem to have any error bars, and so it’s impossible to quantify the significance of this lag. Let’s see what our error bars are:

print(plags_err)

Our errors are all NaNs! This is another consequence of our light curves not being Poissonian. For cross spectral analyses, stingray may therefore not be suitable. For power spectra, however, stingray can be used without issue provided a suitable normalisation is specified (e.g., norm="frac").

That concludes the timing methods tutorial for phoptic! The timing methods of Analyzer demonstrated here are intended as “quick-look” analyses, while more complete analyses are made easy thanks to stingray. Note, however, that stingray may not always work as intended when applied to relative light curves.