This package is a quick adaptation of the RainbowFit function from the original : https://github.com/light-curve/light-curve-python Rainbow here is adapted to fit rising TDE. Temperature evolution is supposed to be constant. Bolometric flux is a sigmoid.
First you should install rust from here : https://www.rust-lang.org/tools/install
Then :
sudo apt-get install cmake
sudo apt-get install libgsl-dev
cd light-curve-python/light-curve
pip install .
Minimum supported Python version is 3.7. We provide binary wheels via PyPi for a number of platforms and architectures, both for CPython and PyPy. We also provide binary wheels for stable CPython ABI, so the package is guaranteed to work with all future CPython3 versions.
Arch \ OS | Linux glibc | Linux musl | macOS | Windows light-curve#186 |
---|---|---|---|---|
x86-64 | wheel (MKL) | wheel (MKL) | wheel | wheel (no Ceres, no GSL) |
i686 | src | src | — | not tested |
aarch64 | wheel | wheel | src light-curve#5 | not tested |
ppc64le | wheel | not tested (no Rust toolchain) | — | — |
- "wheel": binary wheel is available on pypi.org, local building is not required for the platform, the only pre-requirement is a recent
pip
version. For Linux x86-64 we provide binary wheels built with Intel MKL for better periodogram performance, which is not a default build option. For Windows x86-64 we provide wheel with no Ceres and no GSL support, which is not a default build option. - "src": the package is confirmed to be built and pass unit tests locally, but testing and package building is not supported by CI. It is required to have the GNU scientific library (GSL) v2.1+ and the Rust toolchain v1.62+ to install it via
pip install
. - "not tested": building from the source code is not tested, please report us building status via issue/PR/email.
We build aarch64 macOS 12.0+ Python 3.8+ wheels locally and submit them running this command in light-curve
directory:
rm -rf ./wheelhouse
curl https://sh.rustup.rs -sSf | sh -s -- --default-toolchain stable -y
brew install gsl
CIBW_BUILD='cp3*-macosx_arm64' CIBW_ENVIRONMENT="MATURIN_PEP517_ARGS='--locked --no-default-features --features ceres-source,fftw-source,gsl'" CIBW_BEFORE_ALL='' python3 -mcibuildwheel --platform macos
CIBW_BUILD='cp310-macosx_arm64' CIBW_ENVIRONMENT="MATURIN_PEP517_ARGS='--locked --no-default-features --features ceres-source,fftw-source,gsl,abi3'" CIBW_BEFORE_ALL='' python3 -mcibuildwheel --platform macos
twine upload wheelhouse/*.whl
Most of the classes implement various feature evaluators useful for light-curve based astrophysical source classification and characterisation.
import light_curve as lc
import numpy as np
# Time values can be non-evenly separated but must be an ascending array
n = 101
t = np.linspace(0.0, 1.0, n)
perfect_m = 1e3 * t + 1e2
err = np.sqrt(perfect_m)
m = perfect_m + np.random.normal(0, err)
# Half-amplitude of magnitude
amplitude = lc.Amplitude()
# Fraction of points beyond standard deviations from mean
beyond_std = lc.BeyondNStd(nstd=1)
# Slope, its error and reduced chi^2 of linear fit
linear_fit = lc.LinearFit()
# Feature extractor, it will evaluate all features in more efficient way
extractor = lc.Extractor(amplitude, beyond_std, linear_fit)
# Array with all 5 extracted features
result = extractor(t, m, err, sorted=True, check=False)
print('\n'.join(f"{name} = {value:.2f}" for name, value in zip(extractor.names, result)))
# Run in parallel for multiple light curves:
results = amplitude.many(
[(t[:i], m[:i], err[:i]) for i in range(n // 2, n)],
n_jobs=-1,
sorted=True,
check=False,
)
print("Amplitude of amplitude is {:.2f}".format(np.ptp(results)))
If you're confident in your inputs you could use sorted = True
(t
is in ascending order)
and check = False
(no NaNs in inputs, no infs in t
or m
) for better performance.
Note that if your inputs are not valid and are not validated by
sorted=None
and check=True
(default values) then all kind of bad things could happen.
Print feature classes list
import light_curve as lc
print([x for x in dir(lc) if hasattr(getattr(lc, x), "names")])
Read feature docs
import light_curve as lc
help(lc.BazinFit)
See the complete list of available feature evaluators and documentation in light-curve-feature
Rust crate docs.
Feature name | Description | Min data points | Features number |
---|---|---|---|
Amplitude | Half amplitude of magnitude: |
1 | 1 |
AndersonDarlingNormal | Unbiased Anderson–Darling normality test statistic:
where |
4 | 1 |
BazinFit | Five fit parameters and goodness of fit (reduced |
6 | 1 |
BeyondNStd | Fraction of observations beyond |
2 | 1 |
Cusum | A range of cumulative sums:
|
2 | 1 |
Eta | Von Neummann |
2 | 1 |
EtaE | Modernisation of Eta for unevenly time series:
|
2 | 1 |
ExcessVariance | Measure of the variability amplitude:
|
2 | 1 |
InterPercentileRange |
|
1 | 1 |
Kurtosis | Excess kurtosis of magnitude:
|
4 | 1 |
LinearFit | The slope, its error and reduced |
3 | 3 |
LinearTrend | The slope and its error of the light curve in the linear fit of a magnitude light curve without respect to the observation error |
2 | 2 |
MagnitudePercentageRatio | Magnitude percentage ratio:
|
1 | 1 |
MaximumSlope | Maximum slope between two sub-sequential observations:
|
2 | 1 |
Mean | Mean magnitude:
|
1 | 1 |
MeanVariance | Standard deviation to mean ratio:
|
2 | 1 |
Median | Median magnitude | 1 | 1 |
MedianAbsoluteDeviation | Median of the absolute value of the difference between magnitude and its median:
|
1 | 1 |
MedianBufferRangePercentage | 1 | 1 | |
OtsuSplit | Difference of subset means, standard deviation of the lower subset, standard deviation of the upper
subset and lower-to-all observation count ratio for two subsets of magnitudes obtained by Otsu's method split.
Otsu's method is used to perform automatic thresholding. The algorithm returns a single threshold that separate values into two classes. This threshold is determined by minimizing intra-class intensity variance |
2 | 4 |
PercentAmplitude | Maximum deviation of magnitude from its median:
|
1 | 1 |
PercentDifferenceMagnitudePercentile | Ratio of |
1 | 1 |
RainbowFit | Four fit parameters and goodness of fit (reduced |
6 | 1 |
ReducedChi2 | Reduced |
2 | 1 |
Skew | Skewness of magnitude:
|
3 | 1 |
StandardDeviation | Standard deviation of magnitude:
|
2 | 1 |
StetsonK |
Stetson K coefficient described light curve shape:
|
2 | 1 |
VillarFit | Seven fit parameters and goodness of fit (reduced |
8 | 8 |
WeightedMean | Weighted mean magnitude:
|
1 | 1 |
Meta-features can accept other feature extractors and apply them to pre-processed data.
This feature transforms time-series data into the Lomb-Scargle periodogram, providing an estimation of the power spectrum. The peaks argument corresponds to the number of the most significant spectral density peaks to return. For each peak, its period and "signal-to-noise" ratio are returned.
The optional features argument accepts a list of additional feature evaluators, which are applied to the power spectrum: frequency is passed as "time," power spectrum is passed as "magnitude," and no uncertainties are set.
Binning time series to bins with width
Binned time series is defined by
As of v0.8, experimental extractors (see below), support multi-band light-curve inputs.
import numpy as np
from light_curve.light_curve_py import LinearFit
t = np.arange(20, dtype=float)
m = np.arange(20, dtype=float)
sigma = np.full_like(t, 0.1)
bands = np.array(["g"] * 10 + ["r"] * 10)
feature = LinearFit(bands=["g", "r"])
values = feature(t, m, sigma, bands)
print(values)
Rainbow (Russeil+23) is a black-body parametric model for transient light curves.
It uses Bazin function as a model for bolometric flux evolution and a logistic function for the temperature evolutiion.
This example demonstrates the reconstruction of a syntetic light curve with this model.
RainbowFit
requires iminuit
package and Python 3.8+.
import numpy as np
from light_curve.light_curve_py import RainbowFit
def bb_nu(wave_aa, T):
"""Black-body spectral model"""
nu = 3e10 / (wave_aa * 1e-8)
return 2 * 6.626e-27 * nu**3 / 3e10**2 / np.expm1(6.626e-27 * nu / (1.38e-16 * T))
# Effective wavelengths in Angstrom
band_wave_aa = {"g": 4770.0, "r": 6231.0, "i": 7625.0, "z": 9134.0}
# Parameter values
reference_time = 60000.0 # time close to the peak time
# Bolometric flux model parameters
amplitude = 1.0 # bolometric flux semiamplitude, arbitrary (non-spectral) flux/luminosity units
rise_time = 5.0 # exponential growth timescale, days
fall_time = 30.0 # exponential decay timescale, days
# Temperature model parameters
Tmin = 5e3 # temperature on +infinite time, kelvins
delta_T = 10e3 # (Tmin + delta_T) is temperature on -infinite time, kelvins
k_sig = 4.0 # temperature evolution timescale, days
rng = np.random.default_rng(0)
t = np.sort(rng.uniform(reference_time - 3 * rise_time, reference_time + 3 * fall_time, 1000))
band = rng.choice(list(band_wave_aa), size=len(t))
waves = np.array([band_wave_aa[b] for b in band])
# Temperature evolution is a sigmoid function
temp = Tmin + delta_T / (1.0 + np.exp((t - reference_time) / k_sig))
# Bolometric flux evolution is the Bazin function
lum = amplitude * np.exp(-(t - reference_time) / fall_time) / (1.0 + np.exp(-(t - reference_time) / rise_time))
# Spectral flux density for each given pair of time and passband
flux = np.pi * bb_nu(waves, temp) / (5.67e-5 * temp**4) * lum
# S/N = 5 for minimum flux, scale for Poisson noise
flux_err = np.sqrt(flux * np.min(flux) / 5.0)
flux += rng.normal(0.0, flux_err)
feature = RainbowFit.from_angstrom(band_wave_aa, with_baseline=False)
values = feature(t, flux, sigma=flux_err, band=band)
print(dict(zip(feature.names, values)))
print(f"Goodness of fit: {values[-1]}")
Note, that while we don't use precise physical constant values to generate the data, RainbowFit
uses CODATA 2018 values.
From the technical point of view the package consists of two parts: a wrapper for light-curve-feature
Rust crate (light_curve_ext
sub-package) and pure Python sub-package light_curve_py
.
We use the Python implementation of feature extractors to test Rust implementation and to implement new experimental extractors.
Please note, that the Python implementation is much slower for most of the extractors and doesn't provide the same functionality as the Rust implementation.
However, the Python implementation provides some new feature extractors you can find useful.
You can manually use extractors from both implementations:
import numpy as np
from numpy.testing import assert_allclose
from light_curve.light_curve_ext import LinearTrend as RustLinearTrend
from light_curve.light_curve_py import LinearTrend as PythonLinearTrend
rust_fe = RustLinearTrend()
py_fe = PythonLinearTrend()
n = 100
t = np.sort(np.random.normal(size=n))
m = 3.14 * t - 2.16 + np.random.normal(size=n)
assert_allclose(rust_fe(t, m), py_fe(t, m),
err_msg="Python and Rust implementations must provide the same result")
This should print a warning about experimental status of the Python class
You can run all benchmarks from the Python project folder with python3 -mpytest --benchmark-enable tests/test_w_bench.py
, or with slow benchmarks disabled python3 -mpytest -m "not (nobs or multi)" --benchmark-enable tests/test_w_bench.py
.
Here we benchmark the Rust implementation (rust
) versus feets
package and our own Python implementation (lc_py
) for a light curve having n=1000 observations.
The plot shows that the Rust implementation of the package outperforms other ones by a factor of 1.5—50.
This allows to extract a large set of "cheap" features well under one ms for n=1000.
The performance of parametric fits (BazinFit
and VillarFit
) and Periodogram
depend on their parameters, but the typical timescale of feature extraction including these features is 20—50 ms for few hundred observations.
Benchmark results of several features for both the pure-Python and Rust implementations of the "light-curve" package, as a function of the number of observations in a light curve. Both the x-axis and y-axis are on a logarithmic scale.
Processing time per a single light curve for extraction of features subset presented in first benchmark versus the number of CPU cores used. The dataset consists of 10,000 light curves with 1,000 observations in each.
See benchmarks' descriptions in more details in "Performant feature extraction for photometric time series".
Class DmDt
provides dm–dt mapper (based on Mahabal et al. 2011, Soraisam et al. 2020). It is a Python wrapper for light-curve-dmdt
Rust crate.
import numpy as np
from light_curve import DmDt
from numpy.testing import assert_array_equal
dmdt = DmDt.from_borders(min_lgdt=0, max_lgdt=np.log10(3), max_abs_dm=3, lgdt_size=2, dm_size=4, norm=[])
t = np.array([0, 1, 2], dtype=np.float32)
m = np.array([0, 1, 2], dtype=np.float32)
desired = np.array(
[
[0, 0, 2, 0],
[0, 0, 0, 1],
]
)
actual = dmdt.points(t, m)
assert_array_equal(actual, desired)
If you found this project useful for your research please cite Malanchev et al., 2021
@ARTICLE{2021MNRAS.502.5147M,
author = {{Malanchev}, K.~L. and {Pruzhinskaya}, M.~V. and {Korolev}, V.~S. and {Aleo}, P.~D. and {Kornilov}, M.~V. and {Ishida}, E.~E.~O. and {Krushinsky}, V.~V. and {Mondon}, F. and {Sreejith}, S. and {Volnova}, A.~A. and {Belinski}, A.~A. and {Dodin}, A.~V. and {Tatarnikov}, A.~M. and {Zheltoukhov}, S.~G. and {(The SNAD Team)}},
title = "{Anomaly detection in the Zwicky Transient Facility DR3}",
journal = {\mnras},
keywords = {methods: data analysis, astronomical data bases: miscellaneous, stars: variables: general, Astrophysics - Instrumentation and Methods for Astrophysics, Astrophysics - Solar and Stellar Astrophysics},
year = 2021,
month = apr,
volume = {502},
number = {4},
pages = {5147-5175},
doi = {10.1093/mnras/stab316},
archivePrefix = {arXiv},
eprint = {2012.01419},
primaryClass = {astro-ph.IM},
adsurl = {https://ui.adsabs.harvard.edu/abs/2021MNRAS.502.5147M},
adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}