Skip to content

Commit

Permalink
Documentation, Algorithm Fixes, API Change
Browse files Browse the repository at this point in the history
  • Loading branch information
@pupperemeritus
pupperemeritus committed Aug 23, 2023
1 parent 0942450 commit a415e9b
Show file tree
Hide file tree
Showing 2 changed files with 156 additions and 86 deletions.
150 changes: 90 additions & 60 deletions stingray/fourier.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1932,11 +1932,12 @@ def lsft_fast(
t: npt.ArrayLike,
freqs: npt.ArrayLike,
sign: Optional[int] = 1,
fullspec: Optional[bool] = False,
oversampling: Optional[int] = 5,
):
) -> npt.NDArray:
"""
Calculates the Lomb-Scargle Fourier transform of a light curve.
Only considers non-negative frequencies.
Subtracts mean from data as it is required for the working of the algorithm.
Parameters
----------
Expand All @@ -1963,46 +1964,42 @@ def lsft_fast(
ft_res : numpy.ndarray
An array of Fourier transformed data.
"""
if sign not in [1, -1]:
raise ValueError("sign must be 1 or -1")

y_ = copy.deepcopy(y) - np.mean(y)
freqs = freqs[freqs >= 0]
# Constants initialization
sum_xx = np.sum(y)
num_xt = len(y)
sum_xx = np.sum(y_)
num_xt = len(y_)
num_ww = len(freqs)
const1 = np.sqrt(0.5) * np.sqrt(num_ww)
const2 = const1

# Arrays initialization
ft_real = ft_imag = np.zeros((num_ww))
f0, df, N = freqs[0], freqs[1] - freqs[0], len(freqs)

# Sum (y_i * cos(wt - wtau))
Sh, Ch = trig_sum(t, y, df, N, f0, oversampling=oversampling)
Sh, Ch = trig_sum(t, y_, df, N, f0, oversampling=oversampling)

# Summation of (cos(wt - wtau))^2 and (sin(wt - wtau))^2
_, scos2x = trig_sum(
t,
np.ones_like(y) / len(y),
df,
N,
f0,
freq_factor=2,
oversampling=oversampling,
)
# cos^2(x) = (1 + cos(2x))/2
# sin^2(x) = (1 - cos(2x))/2
C2, S2 = (1 + scos2x) / 2, (1 - scos2x) / 2
# Angular Frequency
w = freqs * 2 * np.pi

# Summation of cos(2wt) and sin(2wt)
csum, ssum = trig_sum(
t, np.ones_like(y) / len(y), df, N, f0, freq_factor=2, oversampling=oversampling
t, np.ones_like(y_) / len(y_), df, N, f0, freq_factor=2, oversampling=oversampling
)

wtau = 0.5 * np.arctan2(ssum, csum)

S2, C2 = trig_sum(
t,
np.ones_like(y_),
df,
N,
f0,
freq_factor=2,
oversampling=oversampling,
)

const = np.sqrt(0.5) * np.sqrt(num_ww)

coswtau = np.cos(wtau)
sinwtau = np.sin(wtau)

Expand All @@ -2015,32 +2012,29 @@ def lsft_fast(
scos2 = 0.5 * (N + C2 * cos2wtau + S2 * sin2wtau)
ssin2 = 0.5 * (N - C2 * cos2wtau - S2 * sin2wtau)

# Calculating the real and imaginary components of the Fourier transform
ft_real = const1 * sumr / np.sqrt(scos2)
ft_imag = const2 * sumi / np.sqrt(ssin2)
ft_real = const * sumr / np.sqrt(scos2)
ft_imag = const * np.sign(sign) * sumi / np.sqrt(ssin2)

# Looping through all the frequencies
ft_real[0] = sum_xx / np.sqrt(num_xt)
ft_imag[0] = 0
ft_real[freqs == 0] = sum_xx / np.sqrt(num_xt)
ft_imag[freqs == 0] = 0

# Resultant fourier transform for current angular frequency
ft_res = np.complex128(ft_real + (1j * ft_imag)) * np.exp(-1j * w * t[0])
phase = wtau - w * t[0]

if sign == -1:
ft_res = np.conjugate(ft_res)
ft_res = np.complex128(ft_real + (1j * ft_imag)) * np.exp(-1j * phase)

return ft_res if fullspec else ft_res[freqs > 0]
return ft_res


def lsft_slow(
y: npt.ArrayLike,
t: npt.ArrayLike,
freqs: npt.ArrayLike,
sign: Optional[int] = 1,
fullspec: Optional[bool] = False,
):
) -> npt.NDArray:
"""
Calculates the Lomb-Scargle Fourier transform of a light curve.
Only considers non-negative frequencies.
Subtracts mean from data as it is required for the working of the algorithm.
Parameters
----------
Expand All @@ -2064,43 +2058,79 @@ def lsft_slow(
ft_res : numpy.ndarray
An array of Fourier transformed data.
"""
if sign not in (1, -1):
raise ValueError("sign must be 1 or -1")

# Constants initialization
const1 = np.sqrt(0.5) * np.sqrt(len(y))
const2 = const1
y_ = copy.deepcopy(y) - np.mean(y)
freqs = freqs[freqs >= 0]

ft_real = np.zeros_like(freqs)
ft_imag = np.zeros_like(freqs)
ft_res = np.zeros_like(freqs, dtype=np.complex128)

for i, freq in enumerate(freqs):
wrun = freq * 2 * np.pi
if i == 0:
ft_real[i] = sum(y) / np.sqrt(len(y))
ft_imag[i] = 0
num_y = len(y_)
num_freqs = len(freqs)
sum_y = np.sum(y_)
const1 = np.sqrt(0.5) * np.sqrt(num_y)
const2 = const1 * np.sign(sign)
ft_real = ft_imag = np.zeros(num_freqs)
ft_res = np.zeros(num_freqs, dtype=np.complex128)

for i in range(num_freqs):
wrun = freqs[i] * 2 * np.pi
if wrun == 0:
ft_real = sum_y / np.sqrt(num_y)
ft_imag = 0
phase_this = 0
else:
csum = np.sum(np.cos(2.0 * wrun * t))
ssum = np.sum(np.sin(2.0 * wrun * t))

watan = np.arctan2(ssum, csum)
wtau = 0.5 * watan

cos_wt_wtau = np.cos(wrun * t - wtau)
sin_wt_wtau = np.sin(wrun * t - wtau)
sumr = np.sum(y_ * np.cos(wrun * t - wtau))
sumi = np.sum(y_ * np.sin(wrun * t - wtau))

sumr = np.sum(y * cos_wt_wtau)
sumi = np.sum(y * sin_wt_wtau)
scos2 = np.sum(np.power(np.cos(wrun * t - wtau), 2))
ssin2 = np.sum(np.power(np.sin(wrun * t - wtau), 2))

scos2 = np.sum(np.power(cos_wt_wtau, 2))
ssin2 = np.sum(np.power(sin_wt_wtau, 2))
ft_real = const1 * sumr / np.sqrt(scos2)
ft_imag = const2 * sumi / np.sqrt(ssin2)
phase_this = wtau - wrun * t[0]
ft_res[i] = np.complex128(ft_real + (1j * ft_imag)) * np.exp(-1j * phase_this)
return ft_res

ft_real[i] = const1 * sumr / np.sqrt(scos2)
ft_imag[i] = const2 * sumi / np.sqrt(ssin2)

ft_res[i] = np.complex128(ft_real[i] + (1j * ft_imag[i])) * np.exp(-1j * wrun * t[0])
if sign == -1:
ft_res = np.conjugate(ft_res)
def impose_symmetry_lsft(
ft_res: npt.ArrayLike,
sum_y: float,
len_y: int,
freqs: npt.ArrayLike,
) -> npt.ArrayLike:
"""
Impose symmetry on the input fourier transform.
return ft_res if fullspec else ft_res[freqs > 0]
Parameters
----------
ft_res : np.array
The Fourier transform of the signal.
sum_y : float
The sum of the values of the signal.
len_y : int
The length of the signal.
freqs : np.array
An array of frequencies at which the transform is sampled.
Returns
-------
lsft_res : np.array
The Fourier transform of the signal with symmetry imposed.
freqs_new : np.array
The new frequencies
"""
zero_present = np.amy(freqs == 0)
if not zero_present:
ft_res_new = np.concatenate([np.conjugate(np.flip(ft_res)), 0, ft_res])
freqs_new = np.concatenate([-np.flip(freqs), 0, freqs])
else:
ft_res_new = np.concatenate([np.conjugate(np.flip(ft_res))[1:], ft_res])
freqs_new = np.concatenate([-np.flip(freqs)[1:], freqs])
return ft_res_new, freqs_new
92 changes: 66 additions & 26 deletions stingray/lombscargle.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
from .crossspectrum import Crossspectrum
from .events import EventList
from .exceptions import StingrayError
from .fourier import lsft_fast, lsft_slow
from .fourier import lsft_fast, lsft_slow, impose_symmetry_lsft
from .lightcurve import Lightcurve
from .utils import simon

Expand All @@ -33,7 +33,7 @@ class LombScargleCrossspectrum(Crossspectrum):
norm: {``frac``, ``abs``, ``leahy``, ``none``}, string, optional, default ``none``
The normalization of the cross spectrum.
power_type: {``real``, ``absolute``, ``all`}, string, optional, default ``real``
power_type: {``real``, ``absolute``, ``all`}, string, optional, default ``all``
Parameter to choose among complete, real part and magnitude of the cross spectrum.
fullspec: boolean, optional, default ``False``
Expand Down Expand Up @@ -115,7 +115,7 @@ def __init__(
data1: Optional[Union[EventList, Lightcurve]] = None,
data2: Optional[Union[EventList, Lightcurve]] = None,
norm: Optional[str] = "none",
power_type: Optional[str] = "real",
power_type: Optional[str] = "all",
dt: Optional[float] = None,
fullspec: Optional[bool] = False,
skip_checks: bool = False,
Expand Down Expand Up @@ -453,7 +453,9 @@ def _ls_cross(self, lc1, lc2, freq=None, fullspec=False, method="fast", oversamp
fit_mean=False,
center_data=False,
normalization="psd",
).autofrequency(minimum_frequency=self.min_freq, maximum_frequency=self.max_freq),
).autofrequency(
minimum_frequency=max(self.min_freq, 0), maximum_frequency=self.max_freq
),
)[0]
freqs2 = (
LombScargle(
Expand All @@ -462,34 +464,43 @@ def _ls_cross(self, lc1, lc2, freq=None, fullspec=False, method="fast", oversamp
fit_mean=False,
center_data=False,
normalization="psd",
).autofrequency(minimum_frequency=self.min_freq, maximum_frequency=self.max_freq),
).autofrequency(
minimum_frequency=max(self.min_freq, 0), maximum_frequency=self.max_freq
),
)[0]
if max(freqs2) > max(freq):
freq = freqs2

if method == "slow":
lsft1 = lsft_slow(lc1.counts, lc1.time, freq, sign=-1, fullspec=fullspec)
lsft2 = lsft_slow(lc2.counts, lc2.time, freq, sign=1, fullspec=fullspec)
lsft1 = lsft_slow(lc1.counts, lc1.time, freq)
lsft2 = lsft_slow(lc2.counts, lc2.time, freq)
elif method == "fast":
lsft1 = lsft_fast(
lc1.counts,
lc1.time,
freq,
sign=-1,
fullspec=fullspec,
oversampling=oversampling,
)
lsft2 = lsft_fast(
lc2.counts,
lc2.time,
freq,
fullspec=fullspec,
oversampling=oversampling,
)
cross = np.multiply(lsft1, lsft2)
if not fullspec:
freq = freq[freq > 0]
cross = cross[freq > 0]
if fullspec:
lsft1, _ = impose_symmetry_lsft(
lsft1,
np.sum((lc1.counts)),
lc1.n,
freq,
)
lsft2, freq = impose_symmetry_lsft(
lsft2,
np.sum(lc2.counts),
lc2.n,
freq,
)
cross = np.multiply(lsft1, np.conjugate(lsft2))
return freq, cross

def _initialize_empty(self):
Expand All @@ -515,22 +526,51 @@ def _initialize_empty(self):
self.variance2 = None
return

def phase_lag(self):
"""Not applicable for unevenly sampled data"""
raise AttributeError(
"Object has no attribute named 'phase_lag' ! Not applicable for unevenly sampled data"
)

def time_lag(self):
r"""Calculate the fourier time lag of the cross spectrum.
The time lag is calculated by taking the phase lag :math:`\phi` and
"""Not applicable for unevenly sampled data"""
raise AttributeError(
"Object has no attribute named 'time_lag' ! Not applicable for unevenly sampled data"
)

..math::
def classical_significances(self, threshold=1, trial_correction=False):
"""Not applicable for unevenly sampled data"""
raise AttributeError(
"Object has no attribute named 'classical_significances' ! Not applicable for unevenly sampled data"
)

\tau = \frac{\phi}{\two pi \nu}
def from_time_array():
"""Not applicable for unevenly sampled data"""
raise AttributeError(
"Object has no attribute named 'from_time_array' ! Not applicable for unevenly sampled data"
)

where :math:`\nu` is the center of the frequency bins.
"""
if self.__class__ == LombScargleCrossspectrum:
ph_lag = self.phase_lag()
def from_events():
"""Not applicable for unevenly sampled data"""
raise AttributeError(
"Object has no attribute named 'from_events' ! Not applicable for unevenly sampled data"
)

return ph_lag / (2 * np.pi * self.freq)
else:
raise AttributeError("Object has no attribute named 'time_lag' !")
def from_lightcurve():
"""Not applicable for unevenly sampled data"""
raise AttributeError(
"Object has no attribute named 'from_lightcurve' ! Not applicable for unevenly sampled data"
)

def from_lc_iterable():
"""Not applicable for unevenly sampled data"""
raise AttributeError(
"Object has no attribute named 'from_lc_iterable' ! Not applicable for unevenly sampled data"
)

def _initialize_from_any_input():
"""Not required for unevenly sampled data"""
raise AttributeError("Object has no attribute named '_initialize_from_any_input' !")


class LombScarglePowerspectrum(LombScargleCrossspectrum):
Expand All @@ -552,7 +592,7 @@ class LombScarglePowerspectrum(LombScargleCrossspectrum):
norm: {``frac``, ``abs``, ``leahy``, ``none``}, string, optional, default ``none``
The normalization of the cross spectrum.
power_type: {``real``, ``absolute``, ``all`}, string, optional, default ``real``
power_type: {``real``, ``absolute``, ``all`}, string, optional, default ``all``
Parameter to choose among complete, real part and magnitude of the cross spectrum.
fullspec: boolean, optional, default ``False``
Expand Down

0 comments on commit a415e9b

Please sign in to comment.