DEV Community: pupper emeritus

GSoC Week 11-13

pupper emeritus — Thu, 31 Aug 2023 20:38:33 +0000

I am reaching the final stages of my project well ahead of the schedule.

Both the Crosspectrum and Powerspectrum classes, the Lomb Scargle Fourier Transform and Basic Documentation have been completed. Tests are almost done too

I have double checked with Matteo that the results are consistent using multiple examples and they are close to expected results.

The following are the documentation notebooks on how to use the classes.

GSoC Week 9-10

pupper emeritus — Sun, 13 Aug 2023 10:23:47 +0000

Big progress on the front of algorithm working. Turns out there wasn't that much of a problem in the algorithm. I just had to subtract the mean from the data before taking the fourier transform. The Lomb Scargle seems to work on data that has mean subtracted from it. Furthermore they dont seem to work that well or at all in full spectrum.
The LSFT also is highly sensitive to the time intervals that are input to it. I have found some more clues as to how to make it even better. I will expound upon this further within 2 weeks since my final exams are going on and I have limited time.

Minimum Working Code Example

def data_func(time, freq=1.2324235252):
    return  2 * np.sin(2 * np.pi * time * freq)

t0 = 0
t1 = 100
dt = 0.1

np.random.seed(43)
time_uniform = np.arange(t0, t1, dt)
time_nonuniform = np.sort(np.random.uniform(t0, t1, time_uniform.size))

npft = np.fft.fft(data_func(time_uniform))
npfreqs = np.fft.fftfreq(npft.size, dt)
npft = npft[npfreqs>=0]
npfreqs = npfreqs[npfreqs>=0]
lsfreqs = np.linspace(np.min(npfreqs), np.max(npfreqs), npfreqs.size * 8)
lsfreqs = lsfreqs[lsfreqs>=0]
np.random.seed(43)
lsft_slow_arr = lsft_slow(data_func(time_nonuniform), time_nonuniform, lsfreqs,sign=-1, fullspec=False)
lsft_fast_arr = lsft_fast(data_func(time_nonuniform), time_nonuniform, lsfreqs,sign=-1, fullspec=False,oversampling=10)
plt.plot(time_nonuniform,lsft_slow_inv(lsft_slow_arr,freqs=time_nonuniform, t=lsfreqs).real ,alpha=0.5,label="Slow")
plt.plot(time_nonuniform,lsft_fast_inv(lsft_fast_arr,freqs=time_nonuniform, t=lsfreqs).real,alpha=0.5,label="Fast")
plt.plot(time_nonuniform,data_func(time_nonuniform),label="Original Data")
plt.legend()
plt.xlim(0,10)
plt.ylim(-3,3)

GSoC Week 5-8

pupper emeritus — Sun, 30 Jul 2023 07:35:19 +0000

Brief

I have worked on creating unit tests for the Lomb Scargle Cross Spectrum class, cross verifying the algorithm by comparing with the papers and fixed typos in docstrings. Apologies for the delay. Had my exams.

Details

I have noticed a few faults in both the fast and the slow algorithms. I have gone back to the drawing board and tried to address those issues by following the papers as closely as possible. All the changes are visible in the following draft pull request.
https://github.com/StingraySoftware/stingray/pull/737

After the fixing, the fast and the slow algorithm have started giving very similar outputs. I am starting to suspect that the time lags might be broken in the algorithms themselves. It is starting to get a little suspicious when different methods are giving very similar results and they are still not what that is expected. Last time around the fast and slow algorithms have given different results. After cross verifying with the papers, The results from both have fast and slow algorithms converged.

To keep the project sailing along while I wait for confirmation that this is an issue with the implementation or the algorithm , I have decided to work on writing unit tests for the various classes and methods. Furthermore I also worked on fixing the docstrings.

GSoC Week 4 Update

pupper emeritus — Mon, 03 Jul 2023 19:17:26 +0000

Brief

This week I successfully finished implementing the fast algorithm. Now my LombScargleCrossspectrum and LombScarglePowerspectrum are that much closer to completion. Only things left to sort out/implement are time lags and phase lag functions and checking the phase of the output.

Details

Testing on the following synthetic data has been conducted to compare the outputs with the existing cross spectrum and power spectrum for evenly spaced data first then checking the outputs of the lomb scargle variants on unevenly sampled data

rand = np.random.default_rng(42)
n = 1000
t = np.linspace(0, 10, n)
y = np.sin(2 * np.pi * 3.0 * t) + 0.1 * rand.standard_normal(n)
y2 = np.sin(2 * np.pi * 3.0 * t) + 0.1 * rand.standard_normal(n)
y -= np.min(y)
y2 -= np.min(y2)

The Cross spectra for evenly sampled data

The time lags for evenly sampled data

As it is evident the time lags need work.

The power spectra for evenly sampled data

One quirk is that the power spectrum class is returning the power spectrum with a negative sign. This is a known bug. The values otherwise are within margin of error.

The Lomb Scargle cross spectrum and power spectrum when data is unevenly sampled

t = np.sort(rand.random(n))*10

The cross spectrum

The power spectrum

The time lags

They are off here too. Which will be fixed in the coming week.

For exhaustive testing code refer

GSoC Week 2-3 Update

pupper emeritus — Sun, 25 Jun 2023 16:15:57 +0000

Brief

These weeks I refactored the LombScargleCrossspectrum and LombScarglePowerspectrum classes to accommodate the fast algorithm which went smoothly.
However when it comes to the fast algorithm. I had tunnel vision and unconsciously made the lsft_fast function compute the power spectrum instead of the fourier transform. Right now I am working towards isolating the algorithm to compute the fourier transform using the Press and Rybicki optimizations(https://ui.adsabs.harvard.edu/abs/1989ApJ...338..277P/abstract).

Challenges Faced

Integrating the optimization to the existing slow algorithm is giving me a bit of trouble. I'm still figuring out how to add the optimizations. If this is done, I can move onto making the time lag, phase lag functions and then onto testing and documentation.

Details

Added the following parameters to both the classes in order to accommodate choice between the fast and slow algorithm.

method : str
        The method to be used by the Lomb-Scargle Fourier Transformation function. `fast`
        and `slow` are the allowed values. Default is `fast`. fast uses the optimized Press
        and Rybicki O(n*log(n))

    oversampling : float, optional, default: 5
        Interpolation Oversampling Factor (for the fast algorithm)

For full code refer https://github.com/StingraySoftware/stingray/pull/737

Most important part of the process is the Lomb Scargle Fourier Transform.
The wrapper class is trivial, they only wrap the fast and slow lomb scargle fourier transform functions.

Results using the slow algorithm

On synthetic data

rand = np.random.default_rng(42)
    n = 100
    t = np.sort(rand.random(n)) * n
    y = np.cos(2 * np.pi * 5 * t) + 0.01 * rand.standard_normal(n)
    y -= np.min(y)
    lc1 = Lightcurve(t, y, err_dist="poisson")
    y2 = np.cos(2 * np.pi * 5.0 * (t)) + 0.01 * rand.standard_normal(n)
    y2 -= np.min(y2)
    lc2 = Lightcurve(t, y2, err_dist="poisson")

On real data

The lightcurve

https://heasarc.gsfc.nasa.gov/FTP/nicer/data/obs/2018_03//1200120106/xti/event_cl/ni1200120106_0mpu7_cl.evt.gz

GSoC Week 1 Progress Update

pupper emeritus — Thu, 08 Jun 2023 07:29:01 +0000

Brief

Week 1 has been quite eventful with me creating a Lomb Scargle Fourier Transform function and get a working class for Lomb Scargle Cross Spectrum and Power Spectrum.

I inherited my LS cross spectrum from the regular cross spectrum class. Had to rewrite the constructor, initial_checks , make_crossspectrum, _make_auxil_pds and _initialize_empty. And also wrote a new _ls_cross method which just returns the frequencies and cross spectra for given light curves and it is an internal function only to be used by the class.

The original slow implementation has been completed. I am still working on the fast version.

As this is not a project that can be completed in a bunch of small PRs, I will push to a single PR which will be merged after completion of the project. The following draft PR is the one to which I will be pushing to.

https://github.com/StingraySoftware/stingray/pull/737/

Details

The following are the APIs for the classes

Cross Spectrum

Make a cross spectrum from an unevenly sampled light curve.
You can also make an empty :class:`Crossspectrum` object to populate with your
own Fourier-transformed data (this can sometimes be useful when making
binned power spectra).

Parameters
----------
data1: :class:`stingray.Lightcurve` or :class:`stingray.events.EventList`, optional, default ``None``
    The dataset for the first channel/band of interest.

data2: :class:`stingray.Lightcurve` or :class:`stingray.events.EventList`, optional, default ``None``
    The dataset for the second, or "reference", band.

norm: {``frac``, ``abs``, ``leahy``, ``none``}, default ``none``
    The normalization of the (real part of the) cross spectrum.

power_type: string, optional, default ``real``
    Parameter to choose among complete, real part and magnitude of the cross spectrum.

fullspec: boolean, optional, default ``False``
    If False, keep only the positive frequencies, or if True, keep all of them .

Other Parameters
----------------
dt: float
    The time resolution of the light curve. Only needed when constructing
    light curves in the case where ``data1``, ``data2`` are
    :class:`EventList` objects

skip_checks: bool
    Skip initial checks, for speed or other reasons (you need to trust your
    inputs!)

min_freq : float
    Minimum frequency to take the Lomb-Scargle Fourier Transform

max_freq: float
    Maximum frequency to take the Lomb-Scargle Fourier Transform

df : float
    The time resolution of the light curve. Only needed where ``data1``, ``data2`` are

method : str
    The method to be used by the Lomb-Scargle Fourier Transformation function. `fast`
    and `slow` are the alloowed values. Default is `fast`. fast uses the optimized Press
    and Rybicki O(n*log(n))

Attributes
----------
freq: numpy.ndarray
    The array of mid-bin frequencies that the Fourier transform samples

power: numpy.ndarray
    The array of cross spectra (complex numbers)

power_err: numpy.ndarray
    The uncertainties of ``power``.
    An approximation for each bin given by ``power_err= power/sqrt(m)``.
    Where ``m`` is the number of power averaged in each bin (by frequency
    binning, or averaging more than one spectra). Note that for a single
    realization (``m=1``) the error is equal to the power.

df: float
    The frequency resolution

m: int
    The number of averaged cross-spectra amplitudes in each bin.

n: int
    The number of data points/time bins in one segment of the light
    curves.

k: array of int
    The rebinning scheme if the object has been rebinned otherwise is set to 1.

nphots1: float
    The total number of photons in light curve 1

nphots2: float
    The total number of photons in light curve 2

Power Spectrum

Make a :class:`LombScarglePowerspectrum` (also called periodogram) from a unevenly sampled (binned)
light curve. Periodograms can be normalized by either Leahy normalization,
fractional rms normalization, absolute rms normalization, or not at all.

You can also make an empty :class:`LombScarglePowerspectrum` object to populate with
your own fourier-transformed data (this can sometimes be useful when making
binned power spectra).

Parameters
----------
data: :class:`stingray.lightcurve.Lightcurve` or :class:`stingray.events.EventList` object, optional, default ``None``
    The light curve data to be Fourier-transformed.

norm: {"leahy" | "frac" | "abs" | "none" }, optional, default "frac"
    The normaliation of the power spectrum to be used. Options are
    "leahy", "frac", "abs" and "none", default is "frac".

Other Parameters
----------------
dt: float
    The time resolution of the light curve. Only needed when constructing
    light curves in the case where ``data`` is a
    :class:`EventList` object

skip_checks: bool
    Skip initial checks, for speed or other reasons (you need to trust your
    inputs!).

min_freq : float
    Minimum frequency to take the Lomb-Scargle Fourier Transform

max_freq: float
    Maximum frequency to take the Lomb-Scargle Fourier Transform

df : float
    The time resolution of the light curve. Only needed where ``data`` is a :class`stingray.Eventlist` object

method : str
    The method to be used by the Lomb-Scargle Fourier Transformation function. `fast`
    and `slow` are the alloowed values. Default is `fast`. fast uses the optimized Press
    and Rybicki O(n*log(n))

Attributes
----------
norm: {"leahy" | "frac" | "abs" | "none" }
    The normalization of the power spectrum.

freq: numpy.ndarray
    The array of mid-bin frequencies that the Fourier transform samples.

power: numpy.ndarray
    The array of normalized squared absolute values of Fourier
    amplitudes.

power_err: numpy.ndarray
    The uncertainties of ``power``.
    An approximation for each bin given by ``power_err= power/sqrt(m)``.
    Where ``m`` is the number of power averaged in each bin (by frequency
    binning, or averaging power spectra of segments of a light curve).
    Note that for a single realization (``m=1``) the error is equal to the
    power.

df: float
    The frequency resolution.

m: int
    The number of averaged powers in each bin.

n: int
    The number of data points in the light curve.

nphots: float
    The total number of photons in the light curve.

Community Bonding Period

pupper emeritus — Thu, 01 Jun 2023 03:11:02 +0000

Brief

Hello everyone this is Sri Guru Datta P.(pupperemeritus). I am extremely excited to get started with my GSoC journey. It has been a dream come true for me.

Most of the progress in these weeks has been done on the fronts of

Re-adjusting the proposed schedule : Re-evaluating and taking a closer look at requirements
Gathering more information required to implement : Been provided reference implementation in MATLAB by Dr. Jeff Scargle himself. And validating the plan and research done during the period of creating the project proposal.
Had weekly meets to break the ice and communicate more directly through weekly meets with mentors.

Environment setup

I had already completed my environment setup earlier in order to make pull requests to stingray. I just followed the documentation.
https://docs.stingray.science/install.html

Codebase familiarity

I already had a bit of familiarity with powerspectrum, crossspectrum and bexvar modules as I had worked on various pull requests on them. My mentor advised me to look into multitaper, which I did. It gave me a few ideas on how I could implement this project.

Coding

I had already started a few days early and went ahead and implemented the Lomb-Scargle Fourier Transform which will be the backbone of the power spectrum and cross spectrum. 2 days into the coding period and I also have started coding up the wrapper classes. It turns out to be a very tangled process in order to get every small nook and cranny of the class to work. But I have seemed to get the hang of it.

GSoC Blog Zero

pupper emeritus — Sat, 06 May 2023 21:49:25 +0000

This is a test post to check the status of the blog and RSS feed. I am very excited to contribute to Stingray @ Open Astronomy this summer.