[SSO] Performances

[JulienPeloton]

SHG1G2

Code in profiling.py:

import io
import requests
import pandas as pd
import numpy as np

from fink_utils.sso.spins import estimate_sso_params

r = requests.post(
    "https://api.ztf.fink-portal.org/api/v1/sso",
    json={"n_or_d": "223", "withEphem": True, "output-format": "json"},
)

pdf = pd.read_json(io.BytesIO(r.content))

shg1g2 = estimate_sso_params(
    pdf["i:magpsf_red"].values,
    pdf["i:sigmapsf"].values,
    np.deg2rad(pdf["Phase"].values),
    pdf["i:fid"].values,
    np.deg2rad(pdf["i:ra"].values),
    np.deg2rad(pdf["i:dec"].values),
    model="SHG1G2",
    normalise_to_V=False,
    remap=True,
)

Run with:

kernprof -l profiling.py
python -m line_profiler -rmt profiling.py.lprof

# main
  0.56 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:175 - func_hg1g2
  0.57 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:205 - func_shg1g2
  0.63 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:1533 - build_eqs_for_spins
  0.66 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:2332 - fit_spin
  0.66 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:1757 - estimate_sso_params

versus

# branch issue/264/perf
  0.01 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:171 - func_hg1g2
  0.02 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:198 - func_shg1g2
  0.05 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:1549 - build_eqs_for_spins
  0.08 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:2368 - fit_spin
  0.08 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:1784 - estimate_sso_params

Test suite has been improved to further validate code change.

SOCCA

I use the following script:

import io
import requests
import pandas as pd
import numpy as np
from astropy.coordinates import SkyCoord
import astropy.units as u

from fink_utils.sso.utils import compute_light_travel_correction
from fink_utils.sso.cleaning import dxy_cleaning, iterative_cleaning

from asteroid_spinprops.ssolib import modelfit

ZTF_PIX_SIZE_ARCSEC = 1.0
FILTERING = False

SSOID = "2000 HR66"
r = requests.post(
    "https://api.ztf.fink-portal.org/api/v1/sso",
    json={"n_or_d": SSOID, "withEphem": True, "output-format": "json"},
)

pdf = pd.read_json(io.BytesIO(r.content))

jd_lt = compute_light_travel_correction(
    pdf["i:jd"],
    pdf["Dobs"],
)

# dx, dy from deviation to ephemerides
position = SkyCoord(ra=pdf["i:ra"], dec=pdf["i:dec"], unit="deg")
ephemeride = SkyCoord(ra=pdf["RA"], dec=pdf["DEC"], unit="deg")

dra, ddec = position.spherical_offsets_to(ephemeride)

dx = dra.to(u.arcsec) * ZTF_PIX_SIZE_ARCSEC
dy = ddec.to(u.arcsec) * ZTF_PIX_SIZE_ARCSEC


pdf = pd.DataFrame({
    "cmred": pdf["i:magpsf_red"],
    "csigmapsf": pdf["i:sigmapsf"],
    "Phase": pdf["Phase"],
    "cfid": pdf["i:fid"],
    "ra": pdf["i:ra"],
    "dec": pdf["i:dec"],
    "cjd": jd_lt,
    "i:raephem": pdf["RA"],
    "i:decephem": pdf["DEC"],
    "ra_s": pdf["RA_h"],
    "dec_s": pdf["DEC_h"],
    "cdx": dx,
    "cdy": dy,
    "Dhelio": pdf["Dhelio"],
})
pdf = pdf.sort_values("cjd")

if FILTERING:
    pdf["dxy"] = np.sqrt(pdf["cdx"] ** 2 + pdf["cdy"] ** 2)
    pdf, _ = dxy_cleaning(
        pdf,
        pdf["dxy"],
        pdf["cmred"],
        threshold=0.95,
    )

    pdf, _ = iterative_cleaning(
        pdf,
        pdf["cmred"],
        pdf["csigmapsf"],
        pdf["Phase"],
        pdf["cfid"],
        pdf["ra"],
        pdf["dec"],
    )

# Wrap columns inplace
pdf_transposed = pd.DataFrame({
    colname: [pdf[colname].to_numpy()] for colname in pdf.columns
})

base_kwargs = dict(
    use_angles=True,
    use_filter_dependent=True,
    use_phase=True,
    use_shape=True,
)

current_kwargs = base_kwargs.copy()

outdic = modelfit.get_fit_params(
    data=pdf_transposed,
    flavor="SOCCA",
    shg1g2_constrained=True,
    period_blind=True,
    pole_blind=False,
    period_in=None,
    period_quality_flag=True,
    terminator=True,
    time_me=True,
    remap=True,
    remap_kwargs=current_kwargs,
)

Run with:

kernprof -l profiling_socca.py
python -m line_profiler -rmt profiling.py.lprof

and here are profiling results (I checked parameter values are the same):

# main
  0.00 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:146 - func_hg12
  0.64 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:1536 - build_eqs_for_spins
  1.71 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:207 - func_shg1g2
  9.19 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:1610 - build_eqs_for_spin_shape
 10.77 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:2338 - fit_spin
 10.77 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:1762 - estimate_sso_params
 11.62 seconds - /home/peloton/codes/asteroid-spinprops/asteroid_spinprops/ssolib/periodest.py:34 - get_period_estimate
 11.67 seconds - /home/peloton/codes/asteroid-spinprops/asteroid_spinprops/ssolib/periodest.py:329 - perform_residual_resampling
 26.98 seconds - /home/peloton/codes/asteroid-spinprops/asteroid_spinprops/ssolib/modelfit.py:18 - get_fit_params

# branch issue/264/perf -- speed up from func spins only
  0.02 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:1547 - build_eqs_for_spins
  0.09 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:198 - func_shg1g2
  0.17 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:171 - func_hg1g2
  2.59 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:1630 - build_eqs_for_spin_shape
  3.49 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:2359 - fit_spin
  3.49 seconds - /home/peloton/codes/fink-utils/fink_utils/sso/spins.py:1789 - estimate_sso_params
 11.59 seconds - /home/peloton/codes/asteroid-spinprops/asteroid_spinprops/ssolib/periodest.py:34 - get_period_estimate
 11.65 seconds - /home/peloton/codes/asteroid-spinprops/asteroid_spinprops/ssolib/periodest.py:329 - perform_residual_resampling
 18.57 seconds - /home/peloton/codes/asteroid-spinprops/asteroid_spinprops/ssolib/modelfit.py:20 - get_fit_params

Almost factor of 2 speed-up.

Cleaning procedure

The results above do not use the cleaning procedure from Odysseas. If I trigger it, I get the following results:

Cleaning	rms	f_obs	Period_class	get_fit_params time
True	0.079	10.871	0	33 sec
False	0.091	9.869	0	17 sec

Values are close, but get_fit_params struggles a lot to converge if I get supposedly cleaned values! Note that I am using ZTF only, so maybe the cleaning procedure decimates too much the lightcurve, hence the fits struggles (for 2000 HR66 used here, we have initially 21 measurements in g and 46 measurements in r). I tried on Wallenquist (173 in g, 228 in r): the timings are the same with or without cleaning, but the rms is smaller with cleaning (ouf !):

Cleaning	rms	f_obs	Period_class	get_fit_params time
True	0.076	8.715	0	39 sec
False	0.080	8.715	0	38 sec

So we should keep an eye on undersampled lightcurves which might be artificially degraded by cleaning procedures (which are relative and not absolute).

Further speed-up

At this point, one of the limiting function was model.autopower which was called 50+ times. One of the parameter of that function is samples_per_peak (default is 5). If I modify this parameter, I get (results using 2000 HR66 and cleaning procedures):

samples_per_peak	rms	f_obs	Period_class	get_fit_params time
5 (default)	0.079	10.871	0	33 sec
3	0.087	9.869	0	24 sec
2	0.079	10.872	1	12 sec

Note that the case for 5 and 2 samples per peak give the same period, but Period_class is different (and there is a factor 2 in terms of performance!). Weird?

Note that if I use no cleaning + samples_per_peak = 2, I get the same result than with cleaning, which might indicates that samples_per_peak = 5 is overkill. And indeed, using a more densely sampled lightcurve (Wallenquist, 173 in g, 228 in r):

samples_per_peak	rms	f_obs	Period_class	get_fit_params time
5 (default)	0.076	8.715	0	41 sec
2	0.076	8.715	0	18 sec

Same results, but 2x faster with 2 samples per peak.

[OdysseasXenos]

The speedup is great for SOCCA and I think it will shine most when we use sHG1G2 for cleaning data - which used to be one of the bottlenecks when dealing with ATLAS.

Regarding cleaning data in the case of ZTF I think it does more harm than good. We have a relatively low number of points per lightcurve (100-300 vs the 1,000-3,000 of ATLAS) and in any of my analyses until this point I don't filter any outliers from ZTF (as is already done for all other models in the ssoft).

On the 5 points per peak, it comes from VanderPlas 2018 but the choice here is the lower limit of a suggestion :

To ensure that our grid sufficiently samples each peak, it is
prudent to oversample by some factor—say, n_0 samples per peak

and:

So what is a good choice for n_0? Values ranging from n_0=5
(Schwarzenberg-Czerny 1996; VanderPlas & Ivezic 2015) to
n_0=10 (Debosscher et al. 2007; Richards et al. 2012) seem to
be common in the literature; for periodograms computed in this
paper, we use n_0=5.

hence the choice for SOCCA.

Given that the samples per peak argument is directly proportional to the number of periodograms evaluatedI'm not surprised to see a significant speed up when it's =2 instead of =5 especially when computing the bootstrap score of the period estimate.

I am unsure of how it will affect the quality of the period evaluations at scale, but I plan on taking a sub-sample of SOCCA solutions using ZTF data that are very near the low period limit of the period parameter space (~2 hours or so) and reavulating them with a different period parameter space (say 0.1 - 2 hours). This is will be detrimental in terms of the number of periodogram evaluations so the only way to remedy that will be by reducing the oversampling of the periodogram peaks. So that can serve as a benchmark for the quality of solutions when the peak sampling is lower.