Add an option to parametrize the shift when computing ephemerides (#148)

* Add an option to parametrize the shift when computing ephemerides

* PEP8

* change the arguments for sso.period.compute_residuals

* Update the way we call compute_residuals

* PEP8

* Fix typo

* Assertion in case API does not answer correctly

Author: JulienPeloton

fink_utils/sso/periods.py

@@ -96,35 +96,49 @@ def extract_physical_parameters(pdf, flavor):
     return outdic
 
 
-def compute_residuals(pdf, flavor, phyparam):
+def compute_residuals(
+    mag_red, phase, filters, flavor, phyparam, ra=None, dec=None, times=None
+):
     """Compute residuals between data and predictions
 
     Parameters
     ----------
-    pdf: pandas DataFrame
-        DataFrame with SSO data from Fink REST API
+    mag_red: array of float
+        Reduced magnitude (m-d_obs) for each observation
+    phase: array of float
+        Phase angle for each observation, in degree
+    filters: array of int
+        Filter ID for each observation
     flavor: str
         Model flavor: SHG1G2, HG1G2, HG12, or HG
     phyparam: dict
         Dictionary containing reduced chi2, and estimated parameters and
         error on each parameters.
+    ra: array of float, optional
+        Array of RA angles, in degree. Only used for flavor in ["SHG1G2", "SSHG1G2"]
+    dec: array of float, optional
+        Array of Dec angles, in degree. Only used for flavor in ["SHG1G2", "SSHG1G2"]
+    times: array of float, optional
+        Array of times, in JD with light travel correction applied.
+        Only used for flavor="SSHG1G2"
 
     Returns
     -------
     pd.Series
         Series containing `observation - model` in magnitude
     """
-    pdf["preds"] = 0.0
-    for filtnum in pdf["i:fid"].unique():
-        if filtnum == 3:
-            continue
-        cond = pdf["i:fid"] == filtnum
+    model = np.zeros_like(phase, dtype=float)
+    for filtnum in np.unique(filters):
+        # if filtnum == 3:
+        #    continue
+        cond = filters == filtnum
 
         if flavor == "SSHG1G2":
             pha = [
-                np.deg2rad(pdf["Phase"][cond]),
-                np.deg2rad(pdf["i:ra"][cond]),
-                np.deg2rad(pdf["i:dec"][cond]),
+                np.deg2rad(phase[cond]),
+                np.deg2rad(ra[cond]),
+                np.deg2rad(dec[cond]),
+                times[cond],
             ]
             preds = func_sshg1g2(
                 pha,
@@ -140,9 +154,9 @@ def compute_residuals(pdf, flavor, phyparam):
             )
         elif flavor == "SHG1G2":
             pha = [
-                np.deg2rad(pdf["Phase"][cond]),
-                np.deg2rad(pdf["i:ra"][cond]),
-                np.deg2rad(pdf["i:dec"][cond]),
+                np.deg2rad(phase[cond]),
+                np.deg2rad(ra[cond]),
+                np.deg2rad(dec[cond]),
             ]
             preds = func_hg1g2_with_spin(
                 pha,
@@ -155,26 +169,26 @@ def compute_residuals(pdf, flavor, phyparam):
             )
         elif flavor == "HG":
             preds = func_hg(
-                np.deg2rad(pdf["Phase"][cond]),
+                np.deg2rad(phase[cond]),
                 phyparam["H_{}".format(filtnum)],
                 phyparam["G_{}".format(filtnum)],
             )
         elif flavor == "HG12":
             preds = func_hg12(
-                np.deg2rad(pdf["Phase"][cond]),
+                np.deg2rad(phase[cond]),
                 phyparam["H_{}".format(filtnum)],
                 phyparam["G12_{}".format(filtnum)],
             )
         elif flavor == "HG1G2":
             preds = func_hg1g2(
-                np.deg2rad(pdf["Phase"][cond]),
+                np.deg2rad(phase[cond]),
                 phyparam["H_{}".format(filtnum)],
                 phyparam["G1_{}".format(filtnum)],
                 phyparam["G2_{}".format(filtnum)],
             )
-        pdf.loc[cond, "preds"] = preds
+        model[cond] = preds
 
-    return pdf["i:magpsf_red"] - pdf["preds"]
+    return mag_red - model
 
 
 @profile
@@ -278,6 +292,8 @@ def estimate_synodic_period(
                 "https://api.fink-portal.org/api/v1/sso",
                 json={"n_or_d": ssnamenr, "withEphem": True, "output-format": "json"},
             )
+
+            assert r.status_code == 200, r.content
         else:
             _LOG.error("You need to specify either `ssnamenr` or `pdf`.")
 
@@ -287,17 +303,26 @@ def estimate_synodic_period(
         # get the physical parameters with the latest data
         phyparam = extract_physical_parameters(pdf, flavor)
 
-    # Compute the residuals (obs - model)
-    residuals = compute_residuals(pdf, flavor, phyparam)
-
     if lt_correction:
         # Speed of light in AU/day
-        time = compute_light_travel_correction(pdf["i:jd"], pdf["Dobs"])
+        times = compute_light_travel_correction(pdf["i:jd"], pdf["Dobs"])
     else:
-        time = pdf["i:jd"]
+        times = pdf["i:jd"]
+
+    # Compute the residuals (obs - model)
+    residuals = compute_residuals(
+        pdf["i:magpsf_red"],
+        pdf["Phase"],
+        pdf["i:fid"],
+        flavor,
+        phyparam,
+        ra=pdf["i:ra"],
+        dec=pdf["i:dec"],
+        times=times,
+    )
 
     model = LombScargleMultiband(
-        time,
+        times,
         residuals,
         pdf["i:fid"],
         pdf["i:sigmapsf"],
@@ -316,7 +341,7 @@ def estimate_synodic_period(
     # TODO: I need to be convinced...
     best_period = 2 / freq_maxpower
 
-    out = model.model(time.to_numpy(), freq_maxpower)
+    out = model.model(times.to_numpy(), freq_maxpower)
     prediction = np.zeros_like(residuals)
     for index, filt in enumerate(pdf["i:fid"].unique()):
         if filt == 3:

fink_utils/sso/utils.py

@@ -222,6 +222,7 @@ def get_miriade_data(
     method="rest",
     parameters=None,
     timeout=30,
+    shift=15.0,
     uid=None,
 ):
     """Add ephemerides information from Miriade to a Pandas DataFrame with SSO lightcurve
@@ -278,6 +279,7 @@ def get_miriade_data(
                 observer=observer,
                 rplane=rplane,
                 tcoor=tcoor,
+                shift=shift,
                 timeout=timeout,
             )
         elif method == "ephemcc":
@@ -287,6 +289,7 @@ def get_miriade_data(
                 observer=observer,
                 rplane=rplane,
                 tcoor=tcoor,
+                shift=shift,
                 parameters=parameters,
                 uid=uid,
             )
@@ -310,6 +313,7 @@ def get_miriade_data(
                         pdf_sub["i:jd"],
                         observer=observer,
                         rplane="2",
+                        shift=shift,
                         timeout=timeout,
                     )
                 elif method == "ephemcc":
@@ -318,6 +322,7 @@ def get_miriade_data(
                         pdf_sub["i:jd"],
                         observer=observer,
                         rplane="2",
+                        shift=shift,
                         parameters=parameters,
                         uid=uid,
                     )