Implementation of the power spectrum inference with GPs

[mlefkir]

This is the Python-JAX implementation of a method to infer the power spectral density of irregular time series using Gaussian process regression. The method is described in a forthcoming paper and in the Julia package Pioran.jl, it relies on approximating a bending power-law model in a sum of scalable kernels implemented in tinygp.

Relevant Issue(s)/PR(s)

Provide an overview of the implemented solution or the fix and elaborate on the modifications.

Is there a new dependency introduced by your contribution? If so, please specify.

Any other comments?

[pep8speaks]

Hello @mlefkir! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

• In the file stingray/modeling/gpmodeling.py:

Line 98:80: E501 line too long (85 > 79 characters) Line 157:80: E501 line too long (95 > 79 characters) Line 185:80: E501 line too long (92 > 79 characters) Line 186:80: E501 line too long (93 > 79 characters) Line 187:80: E501 line too long (83 > 79 characters) Line 205:80: E501 line too long (86 > 79 characters) Line 216:80: E501 line too long (97 > 79 characters) Line 263:80: E501 line too long (83 > 79 characters) Line 266:80: E501 line too long (92 > 79 characters) Line 267:80: E501 line too long (99 > 79 characters) Line 270:80: E501 line too long (89 > 79 characters) Line 325:80: E501 line too long (81 > 79 characters) Line 328:80: E501 line too long (86 > 79 characters) Line 352:80: E501 line too long (94 > 79 characters) Line 390:80: E501 line too long (91 > 79 characters) Line 391:80: E501 line too long (88 > 79 characters) Line 394:80: E501 line too long (94 > 79 characters) Line 396:80: E501 line too long (87 > 79 characters) Line 397:80: E501 line too long (86 > 79 characters) Line 412:80: E501 line too long (87 > 79 characters) Line 413:80: E501 line too long (86 > 79 characters) Line 426:80: E501 line too long (94 > 79 characters) Line 454:80: E501 line too long (92 > 79 characters) Line 484:80: E501 line too long (81 > 79 characters) Line 487:80: E501 line too long (86 > 79 characters) Line 498:80: E501 line too long (95 > 79 characters) Line 539:80: E501 line too long (92 > 79 characters) Line 542:80: E501 line too long (91 > 79 characters) Line 590:80: E501 line too long (86 > 79 characters) Line 594:80: E501 line too long (80 > 79 characters) Line 606:80: E501 line too long (94 > 79 characters) Line 609:80: E501 line too long (96 > 79 characters) Line 614:80: E501 line too long (89 > 79 characters) Line 702:80: E501 line too long (89 > 79 characters) Line 995:80: E501 line too long (92 > 79 characters) Line 1026:80: E501 line too long (83 > 79 characters) Line 1044:80: E501 line too long (81 > 79 characters) Line 1045:80: E501 line too long (89 > 79 characters) Line 1046:80: E501 line too long (87 > 79 characters) Line 1194:80: E501 line too long (81 > 79 characters) Line 1195:80: E501 line too long (89 > 79 characters) Line 1247:80: E501 line too long (94 > 79 characters) Line 1317:80: E501 line too long (91 > 79 characters) Line 1321:80: E501 line too long (92 > 79 characters)

• In the file stingray/modeling/tests/test_gpmodeling.py:

Line 118:80: E501 line too long (81 > 79 characters) Line 145:80: E501 line too long (90 > 79 characters) Line 336:80: E501 line too long (94 > 79 characters) Line 342:80: E501 line too long (80 > 79 characters) Line 347:80: E501 line too long (83 > 79 characters) Line 350:80: E501 line too long (91 > 79 characters) Line 362:80: E501 line too long (98 > 79 characters) Line 369:80: E501 line too long (91 > 79 characters) Line 384:80: E501 line too long (91 > 79 characters) Line 392:80: E501 line too long (90 > 79 characters) Line 400:80: E501 line too long (91 > 79 characters) Line 404:80: E501 line too long (83 > 79 characters) Line 411:80: E501 line too long (89 > 79 characters) Line 412:80: E501 line too long (96 > 79 characters)

Comment last updated at 2024-10-03 11:36:06 UTC

[matteobachetti]

@mlefkir thanks for your contribution to Stingray, and sorry for my late reply! I'm not an expert in this kind of methods (@Gaurav17Joshi implemented the current GP model and might have some comments). Would you mind showing the usage of this method through an example?

From the point of view of the requirements for PRs to Stingray:

1. most formatting issues can be solved by running black on the new files. See https://docs.stingray.science/en/stable/contributing.html#coding-style-and-conventions. Don't worry about PEP8 still complaining after black has run ;)
2. A changelog entry is highly appreciated, it helps keeping track of the evolution of Stingray. https://docs.stingray.science/en/stable/contributing.html#updating-and-maintaining-the-changelog

[Gaurav17Joshi]

Yes, @mlefkir, do you have any use case example for this method (Something like a .ipynb notebook). Also @dhuppenkothen, can you also have a look into the usefulness of this method for Stingray, and whether it should go into the same file as the gpmodeling part?

[mlefkir]

@Gaurav17Joshi @matteobachetti I made two examples available here Examples, one uses nested sampling with the jaxns sampler already called in Stingray and the other one uses NUTS with NumPyro. I will update the changelog at some point.

[matteobachetti]

@mlefkir I'm playing with your PR. Really sorry for the slow progress, but my knowledge of these methods is pretty poor and it takes me a lot of time to just understand how it works, and... my agenda is pretty full 😅 . Just to guide me in the review: does your model also work with the spectral shapes that @Gaurav17Joshi set up? Is it just working with unevenly sampled data with error bars? What is the difference between your implementation and the existing one? Also, feel free to advertise your paper when it is ready!

[mlefkir]

@matteobachetti This model is an improvement on available models. Currently, available models have a fixed low and high-frequency slope. For instance, the DRW/Exponential kernel has a Lorentzian power spectrum with a low-frequency slope 0 and a high-frequency slope of -2. This method allows modelling spectral shapes with flexible bends frequencies and slopes which can be between 0 and 4 using a sum of basis functions as shown in the figure below:

This method is designed for Gaussian process regression so it can be used for any Gaussian time series with (or without) error bars. While the algorithm in tinygp is fast there are limitations on the number of points in terms of computational cost so I would use the method only for irregularly sampled data or data with gaps with less than 10,000 points. This is a simple implementation of the method which should be easier to maintain than the Python library I wrote a few months ago, I intend to archive it at some point.

[mlefkir]

@matteobachetti If you are interested, the preprint of the paper describing the method is available here: https://arxiv.org/pdf/2501.05886

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 24.21053% with 144 lines in your changes missing coverage. Please review.

Project coverage is 94.71%. Comparing base (8896659) to head (30aab7f).

Files with missing lines	Patch %	Lines
stingray/modeling/gpmodeling.py	24.21%	144 Missing :warning:

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #832      +/-   ##
==========================================
- Coverage   96.03%   94.71%   -1.33%     
==========================================
  Files          48       48              
  Lines        9770     9950     +180     
==========================================
+ Hits         9383     9424      +41     
- Misses        387      526     +139     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[matteobachetti]

@mlefkir do you think you can dedicate some time to this in the next weeks/months? It would be a great addition to Stingray

[mlefkir]

@mlefkir do you think you can dedicate some time to this in the next weeks/months? It would be a great addition to Stingray

it is on my list of things to do once I finished writing my thesis. I aim to complete the PR before the end of December.