possible overflow? (cesium lomb_scargle model)

[sarajamal57]

(cesium.features.lomb_scargle.lomb_scargle_model.py)

When varying the number of harmonics (nharm) in function lomb_scargle_model(), the computed model can showcase some instabilities for high nharmonics, not in a regular pattern though (see attached file).
Must specify that these instabilities (when changing nharm) do not appear for all light-curves tested, but still noticed for some data.
In the following, one identified case is reported.

Possible source of error
In real astronomical TS, the flux/mags values would refer to negative values or a different numerical precision that could possibly be the source of underflow/overflow within the optimization part (refine_psd, get_eigs) in _lomb_scargle.h and _eigs.h

Possible solution
normalize the signal before fitting the lomb_scargle (temporary solution?)
As a result, the computed model would be a normalized version of the usual model (=computed on initial mag/flux measurements).

Proposition
within cesium.features.lomb_scargle.lomb_scargle_model.py :

• compute internally the normalization of the lc = signal entry in functions fit_lomb_scargle() and lomb_scargle_model()
• compute the lomb_scargle model using fit_lomb_scargle() and _lomb_scargle.h. Output: normalized fitted model
• compute adjustments (*scale and +mean) to provide the model on initial flux/mags

Attached files
displays of initial light_curve superimposed with the estimated trend from cesium_fit for a varying nharm.
The cesium_period is specified for each run.
When computing the cesium model on normalized data, no instabilities are detected.

snippet:

from cesium.features.lomb_scargle import lomb_scargle_model
lc ## cesium TS object from MACHO, object_name='6.6692.9', red band
times = lc.time.copy(); mags = lc.measurement.copy(); errors = lc.error.copy(); 
sys_err=0.00; nfreq=1; tone_control=5.0
for nharm in range(1,21):
    model_cesium = lomb_scargle_model(times-min(times), mags, errors,
                                         sys_err=sys_err, nharm=nharm, nfreq=nfreq, tone_control=tone_control)
    period_cesium  = 1/model_cesium['freq_fits'][0]['freq']
    trend_cesium   = model_cesium['freq_fits'][0]['trend']