-
Notifications
You must be signed in to change notification settings - Fork 101
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Harmonic oscillation reduction #285
Closed
sarajamal57
wants to merge
11
commits into
cesium-ml:main
from
sarajamal57:harmonic-oscillation-reduction
Closed
Changes from all commits
Commits
Show all changes
11 commits
Select commit
Hold shift + click to select a range
f699f09
rescale Lomb-Scargle model coputed on normalized data
sarajamal57 7f4270f
add the full psd vector in the model outputs
sarajamal57 99191da
psd_freqs dict
sarajamal57 62e5abe
reframe
sarajamal57 535ab8b
reframe
sarajamal57 00c58e3
Update lomb_scargle.py and delete new file
acrellin 5008a4b
Clean up docstrings
acrellin 054c639
Apply black formatting
acrellin 8914967
run black
profjsb 4903cb5
add test for opt_normalize
profjsb eb22bd6
add docstring info about what opt_normalize does
profjsb File filter
Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -4,11 +4,19 @@ | |
|
||
|
||
def lomb_scargle_model( | ||
time, signal, error, sys_err=0.05, nharm=8, nfreq=3, tone_control=5.0 | ||
time, | ||
signal, | ||
error, | ||
sys_err=0.05, | ||
nharm=8, | ||
nfreq=3, | ||
tone_control=5.0, | ||
opt_normalize=False, | ||
): | ||
"""Simultaneous fit of a sum of sinusoids by weighted least squares: | ||
y(t) = Sum_k Ck*t^k + Sum_i Sum_j A_ij sin(2*pi*j*fi*(t-t0)+phi_j), | ||
i=[1,nfreq], j=[1,nharm] | ||
"""Simultaneous fit of a sum of sinusoids by weighted least squares. | ||
|
||
y(t) = Sum_k Ck*t^k + Sum_i Sum_j A_ij sin(2*pi*j*fi*(t-t0)+phi_j), | ||
i=[1,nfreq], j=[1,nharm] | ||
|
||
Parameters | ||
---------- | ||
|
@@ -21,11 +29,22 @@ def lomb_scargle_model( | |
error : array_like | ||
Array containing measurement error values. | ||
|
||
nharm : int | ||
Number of harmonics to fit for each frequency. | ||
sys_err : float, optional | ||
Defaults to 0.05 | ||
|
||
nharm : int, optional | ||
Number of harmonics to fit for each frequency. Defaults to 8. | ||
|
||
nfreq : int, optional | ||
Number of frequencies to fit. Defaults to 3. | ||
|
||
nfreq : int | ||
Number of frequencies to fit. | ||
tone_control : float, optional | ||
Defaults to 5.0 | ||
|
||
opt_normalize : boolean, optional | ||
Normalize the timeseries before fitting? This can | ||
help with instabilities seen at large values of `nharm` | ||
Defaults to False. | ||
|
||
Returns | ||
------- | ||
|
@@ -36,13 +55,18 @@ def lomb_scargle_model( | |
fit_lomb_scargle(...) | ||
|
||
""" | ||
|
||
dy0 = np.sqrt(error**2 + sys_err**2) | ||
|
||
wt = 1.0 / dy0**2 | ||
time = time.copy() - min(time) # speeds up lomb_scargle code to have min(time)==0 | ||
signal = signal.copy() | ||
|
||
# Normalization | ||
mmean = np.nanmean(signal) | ||
mscale = np.nanstd(signal) | ||
if opt_normalize: | ||
signal = (signal - mmean) / mscale | ||
error = error / mscale | ||
|
||
dy0 = np.sqrt(error**2 + sys_err**2) | ||
wt = 1.0 / dy0**2 | ||
chi0 = np.dot(signal**2, wt) | ||
|
||
# TODO parametrize? | ||
|
@@ -71,6 +95,7 @@ def lomb_scargle_model( | |
detrend_order=1, | ||
) | ||
model_dict["trend"] = fit["trend_coef"][1] | ||
# | ||
else: | ||
fit = fit_lomb_scargle( | ||
time, | ||
|
@@ -84,22 +109,75 @@ def lomb_scargle_model( | |
nharm=nharm, | ||
detrend_order=0, | ||
) | ||
model_dict["freq_fits"].append(fit) | ||
signal -= fit["model"] | ||
model_dict["freq_fits"][-1]["resid"] = signal.copy() | ||
|
||
# remove current estim_freq for next run [**NORM] | ||
norm_residual = signal - fit["model"] | ||
signal = norm_residual | ||
|
||
# store current estim_freq results, [**RESCALED for first fund freq] | ||
current_fit = ( | ||
rescale_lomb_model(fit, mmean, mscale) | ||
if (opt_normalize & (i == 0)) | ||
else fit | ||
) | ||
model_dict["freq_fits"].append(current_fit) | ||
model_dict["freq_fits"][-1]["resid"] = norm_residual | ||
|
||
if opt_normalize & (i == 0): | ||
model_dict["freq_fits"][-1]["resid"] *= mscale | ||
|
||
if i == 0: | ||
model_dict["varrat"] = np.dot(signal**2, wt) / chi0 | ||
model_dict["varrat"] = ( | ||
np.dot(norm_residual**2, wt) / chi0 | ||
) # no need to rescale | ||
|
||
model_dict["nfreq"] = nfreq | ||
model_dict["nharm"] = nharm | ||
model_dict["chi2"] = fit["chi2"] | ||
model_dict["chi2"] = current_fit["chi2"] | ||
model_dict["f0"] = f0 | ||
model_dict["df"] = df | ||
model_dict["numf"] = numf | ||
|
||
return model_dict | ||
|
||
|
||
def rescale_lomb_model(norm_model, mmean, mscale): | ||
"""Rescale Lomb-Scargle model if compiled on normalized data. | ||
|
||
Parameters | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. docstring needs formatting |
||
---------- | ||
norm_model : dict | ||
Lomb-Scargle model computed on normalized data | ||
|
||
mmean : float64 | ||
mean of input_signal | ||
|
||
mscale : float64 | ||
standard-deviation of input_signal | ||
|
||
Returns | ||
------- | ||
rescaled_model: dict | ||
acrellin marked this conversation as resolved.
Show resolved
Hide resolved
|
||
rescaled Lomb-Scargle model (in adequacy with the input_signal) | ||
|
||
""" | ||
rescaled_model = norm_model.copy() | ||
# unchanged : 'psd', 'chi0', 'freq', 'chi2', 'trace', 'nu0', 'nu', 'npars', | ||
# 'rel_phase', 'rel_phase_error', 'time0', 'signif' | ||
rescaled_model["s0"] = rescaled_model["s0"] / mscale**2 | ||
rescaled_model["lambda"] = rescaled_model["lambda"] / mscale**2 | ||
rescaled_model["model"] = rescaled_model["model"] * mscale + mmean | ||
rescaled_model["model_error"] = rescaled_model["model_error"] * mscale | ||
rescaled_model["trend"] = rescaled_model["trend"] * mscale + mmean | ||
rescaled_model["trend_error"] = rescaled_model["trend_error"] * mscale | ||
rescaled_model["trend_coef"] = rescaled_model["trend_coef"] * mscale + mmean | ||
rescaled_model["amplitude"] = rescaled_model["amplitude"] * mscale | ||
rescaled_model["amplitude_error"] = rescaled_model["amplitude_error"] * mscale | ||
rescaled_model["y_offset"] = rescaled_model["y_offset"] * mscale | ||
|
||
return rescaled_model | ||
|
||
|
||
def lprob2sigma(lprob): | ||
"""Translate a log_e(probability) to units of Gaussian sigmas.""" | ||
if lprob > -36.0: | ||
|
@@ -131,7 +209,8 @@ def fit_lomb_scargle( | |
lomb_scargle_model in order to produce a fit with multiple distinct | ||
frequencies. | ||
|
||
Inputs: | ||
Parameters | ||
---------- | ||
time : array_like | ||
Array containing time values. | ||
|
||
|
@@ -304,6 +383,12 @@ def fit_lomb_scargle( | |
out_dict["nu"] = ntime - npars | ||
out_dict["npars"] = npars | ||
|
||
allfreqs = [ | ||
f0 + df * ii + (ifreq / freq_zoom - 1 / 2.0) * df for ii in range(len(psd)) | ||
] | ||
out_dict["freqs_vector"] = np.asarray(allfreqs) | ||
out_dict["psd_vector"] = psd | ||
|
||
A0, B0 = soln[0:nharm], soln[nharm:] | ||
hat_hat /= np.outer( | ||
np.hstack(((1.0 + ii) ** 2, (1.0 + ii) ** 2)), | ||
|
@@ -439,3 +524,16 @@ def get_lomb_trend(lomb_model): | |
def get_lomb_y_offset(lomb_model): | ||
"""Get the y-intercept of a fitted Lomb-Scargle model.""" | ||
return lomb_model["freq_fits"][0]["y_offset"] | ||
|
||
|
||
def get_lomb_psd(lomb_model): | ||
""" | ||
Get the power spectrum of a fitted Lomb-Scargle model per fitted frequency. | ||
""" | ||
dict_psd = {} | ||
for i in range(lomb_model["nfreq"]): | ||
dict_psd[i] = { | ||
"freqs_vector": lomb_model["freq_fits"][i]["freqs_vector"], | ||
"psd_vector": lomb_model["freq_fits"][i]["psd_vector"], | ||
} | ||
return dict_psd |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What does opt_normalize do?
We also need tests for these new flags.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The normalization resolves the numerical overflow observed for instance in issue #284
Modified/Added functions
lomb_scargle_model()
rescale_lomb_model()
Procedure
If
opt_normalize==False
, this refers to the original Lomb-Scargle periodogram function.Otherwise:
Step 1 - the time-serie (TS) is normalized, and the (mean, std scale) are stored
Step 2 - Lomb-Scargle periodogram & the best fitted model are computed on the normalized (signal, error)
Step 3 - Results are scaled back to the original/raw format using the stored (mean, scale)
Step 1
In function
lomb_scargle_model()
:Step 2
The lomb-scargle (LS) periodogram is estimated on the normalized time-serie (signal, associated error).
Step 3
The best fitted Lomb-Scargle model is rescaled back.
In function
lomb_scargle_model()
:The first component (
i==0
) of varfreq_fits
requires to re-inject the scale