Skip to content

Commit 705f7ca

Browse files
JulienPelotonJulienPeloton
authored
[ZTF] SOCCA integration (#598)
* Add initial code to run SOCCA. Force installation of asteroid_spinprops at the CI level. * Ignore python version defined in asteroid_spinprops for the installation. * Change according to asteroid_spinprops version * Rename fields * Rename SSHG1G2 into SOCCA * Rmeove unused functions * Test with other data * Run only the SSOFT * Fix path * Fix tests
1 parent 1e0929f commit 705f7ca

File tree

3 files changed

+97
-60
lines changed

3 files changed

+97
-60
lines changed

‎.github/workflows/run_test_ztf.yml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -44,6 +44,7 @@ jobs:
4444
- name: Run test suites for ZTF
4545
run: |
4646
pip install onnxruntime==1.16.3
47+
pip install asteroid_spinprops --no-deps
4748
rm -f /tmp/forest_*.onnx
4849
./run_tests.sh -s ztf
4950
curl -s https://codecov.io/bash | bash

‎fink_science/data/alerts/test_SOCCA.parquet‎

57.8 KB
Binary file not shown.

‎fink_science/ztf/ssoft/processor.py‎

Lines changed: 96 additions & 60 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
# Copyright 2023-2024 AstroLab Software
1+
# Copyright 2023-2025 AstroLab Software
22
# Author: Julien Peloton
33
#
44
# Licensed under the Apache License, Version 2.0 (the "License");
@@ -29,6 +29,7 @@
2929
from fink_utils.sso.spins import estimate_sso_params
3030
from fink_utils.sso.spins import extract_obliquity
3131
from fink_utils.sso.utils import rockify, extract_array_from_series
32+
from fink_utils.sso.utils import compute_light_travel_correction
3233

3334
from fink_science import __file__
3435
from fink_science.tester import spark_unit_tests
@@ -39,6 +40,8 @@
3940
from astropy.coordinates import SkyCoord
4041
import astropy.units as u
4142

43+
from asteroid_spinprops.ssolib import modelfit
44+
4245
import logging
4346

4447

@@ -169,7 +172,7 @@
169172
"version": {"type": "str", "description": "Version of the SSOFT YYYY.MM"},
170173
}
171174

172-
COLUMNS_SSHG1G2 = {
175+
COLUMNS_SOCCA = {
173176
"G1_1": {
174177
"type": "double",
175178
"description": "G1 phase parameter for the ZTF filter band g",
@@ -408,32 +411,6 @@
408411
}
409412

410413

411-
def estimate_period():
412-
"""TBD
413-
414-
Should be done only for SSHG1G2
415-
416-
sb_method: str
417-
Specify the single-band lomb scargle implementation to use.
418-
See https://docs.astropy.org/en/stable/api/astropy.timeseries.LombScargleMultiband.html#astropy.timeseries.LombScargleMultiband.autopower
419-
If nifty-ls is installed, one can also specify fastnifty. Although
420-
in this case it does not work yet for Nterms_* higher than 1.
421-
"""
422-
pass
423-
424-
425-
def extract_ssoft_parameters_sshg1g2():
426-
"""TBD
427-
428-
Notes
429-
-----
430-
For the SSHG1G2 model, the strategy is the following:
431-
1. Compute parameters as if it was SHG2G1 model (incl. period estimation)
432-
2. Using previously computed parameters, compute parameters from SSHG1G2
433-
"""
434-
pass
435-
436-
437414
@pandas_udf(MapType(StringType(), FloatType()), PandasUDFType.SCALAR)
438415
@profile
439416
def extract_ssoft_parameters(
@@ -457,7 +434,7 @@ def extract_ssoft_parameters(
457434
Notes
458435
-----
459436
Only works for HG, HG1G2, and SHG1G2. Rotation period
460-
is not estimated here. For SSHG1G2, see <TBD>
437+
is not estimated here. For SOCCA, see <TBD>
461438
462439
Parameters
463440
----------
@@ -483,7 +460,7 @@ def extract_ssoft_parameters(
483460
Use only the former on the Spark Cluster (local installation of ephemcc),
484461
otherwise use `rest` to call the ssodnet web service.
485462
model: str
486-
Model name. Available: HG, HG1G2, SHG1G2
463+
Model name. Available: HG, HG1G2, SHG1G2, SOCCA
487464
488465
489466
Returns
@@ -502,7 +479,7 @@ def extract_ssoft_parameters(
502479
[30, 1, 1, 1, 2 * np.pi, np.pi / 2],
503480
),
504481
},
505-
"SSHG1G2": {
482+
"SOCCA": {
506483
"p0": [15.0, 0.15, 0.15, np.pi, 0.0, 5.0, 1.05, 1.05, 0.0],
507484
"bounds": (
508485
[-3, 0, 0, 0, -np.pi / 2, 2.2 / 24.0, 1, 1, -np.pi / 2],
@@ -525,35 +502,81 @@ def extract_ssoft_parameters(
525502
extract_array_from_series(dobs, index, float)
526503
* extract_array_from_series(dhelio, index, float)
527504
)
528-
pdf = pd.DataFrame({
529-
"i:ssnamenr": [ssname] * len(raobs.to_numpy()[index]),
530-
"i:magpsf": extract_array_from_series(magpsf, index, float),
531-
"i:sigmapsf": extract_array_from_series(sigmapsf, index, float),
532-
"i:jd": extract_array_from_series(jd, index, float),
533-
"i:fid": extract_array_from_series(fid, index, int),
534-
"i:ra": extract_array_from_series(raobs, index, float),
535-
"i:dec": extract_array_from_series(decobs, index, float),
536-
"i:raephem": extract_array_from_series(raephem, index, float),
537-
"i:decephem": extract_array_from_series(decephem, index, float),
538-
"i:magpsf_red": magpsf_red,
539-
"Phase": extract_array_from_series(phase, index, float),
540-
"Dobs": extract_array_from_series(dobs, index, float),
541-
})
542-
pdf = pdf.sort_values("i:jd")
543-
544-
outdic = estimate_sso_params(
545-
pdf["i:magpsf_red"].to_numpy(),
546-
pdf["i:sigmapsf"].to_numpy(),
547-
np.deg2rad(pdf["Phase"].to_numpy()),
548-
pdf["i:fid"].to_numpy(),
549-
np.deg2rad(pdf["i:ra"].to_numpy()),
550-
np.deg2rad(pdf["i:dec"].to_numpy()),
551-
jd=pdf["i:jd"].to_numpy(),
552-
p0=MODELS[model_name]["p0"],
553-
bounds=MODELS[model_name]["bounds"],
554-
model=model_name,
555-
normalise_to_V=False,
556-
)
505+
if model_name == "SOCCA":
506+
jd_lt = compute_light_travel_correction(
507+
extract_array_from_series(jd, index, float),
508+
extract_array_from_series(dobs, index, float),
509+
)
510+
pdf = pd.DataFrame({
511+
"cmred": magpsf_red,
512+
"csigmapsf": extract_array_from_series(sigmapsf, index, float),
513+
"Phase": extract_array_from_series(phase, index, float),
514+
"cfid": extract_array_from_series(fid, index, int),
515+
"ra": extract_array_from_series(raobs, index, float),
516+
"dec": extract_array_from_series(decobs, index, float),
517+
"cjd": jd_lt,
518+
"i:raephem": extract_array_from_series(raephem, index, float),
519+
"i:decephem": extract_array_from_series(decephem, index, float),
520+
})
521+
pdf = pdf.sort_values("cjd")
522+
523+
# Wrap columns inplace
524+
pdf_transposed = pd.DataFrame({
525+
colname: [pdf[colname].to_numpy()] for colname in pdf.columns
526+
})
527+
528+
# parameter estimation
529+
outdic = modelfit.get_fit_params(
530+
pdf_transposed,
531+
flavor=model_name,
532+
shg1g2_constrained=True,
533+
period_blind=True,
534+
pole_blind=True,
535+
alt_spin=False,
536+
period_in=None,
537+
terminator=False,
538+
)
539+
540+
# replace names inplace for the remaning computation
541+
pdf = pdf.rename(
542+
columns={
543+
"ra": "i:ra",
544+
"dec": "i:dec",
545+
"cfid": "i:fid",
546+
"cjd": "i:jd", # FIXME: this is lighttime corrected
547+
}
548+
)
549+
else:
550+
pdf = pd.DataFrame({
551+
"i:ssnamenr": [ssname] * len(raobs.to_numpy()[index]),
552+
"i:magpsf": extract_array_from_series(magpsf, index, float),
553+
"i:sigmapsf": extract_array_from_series(sigmapsf, index, float),
554+
"i:jd": extract_array_from_series(jd, index, float),
555+
"i:fid": extract_array_from_series(fid, index, int),
556+
"i:ra": extract_array_from_series(raobs, index, float),
557+
"i:dec": extract_array_from_series(decobs, index, float),
558+
"i:raephem": extract_array_from_series(raephem, index, float),
559+
"i:decephem": extract_array_from_series(decephem, index, float),
560+
"i:magpsf_red": magpsf_red,
561+
"Phase": extract_array_from_series(phase, index, float),
562+
"Dobs": extract_array_from_series(dobs, index, float),
563+
})
564+
565+
pdf = pdf.sort_values("i:jd")
566+
567+
outdic = estimate_sso_params(
568+
pdf["i:magpsf_red"].to_numpy(),
569+
pdf["i:sigmapsf"].to_numpy(),
570+
np.deg2rad(pdf["Phase"].to_numpy()),
571+
pdf["i:fid"].to_numpy(),
572+
np.deg2rad(pdf["i:ra"].to_numpy()),
573+
np.deg2rad(pdf["i:dec"].to_numpy()),
574+
jd=pdf["i:jd"].to_numpy(),
575+
p0=MODELS[model_name]["p0"],
576+
bounds=MODELS[model_name]["bounds"],
577+
model=model_name,
578+
normalise_to_V=False,
579+
)
557580

558581
# Add astrometry
559582
fink_coord = SkyCoord(
@@ -674,6 +697,19 @@ def build_the_ssoft(
674697
>>> col_ssoft_shg1g2 = sorted(ssoft_shg1g2.columns)
675698
>>> expected_cols = sorted({**COLUMNS, **COLUMNS_SHG1G2}.keys())
676699
>>> assert col_ssoft_shg1g2 == expected_cols, (col_ssoft_shg1g2, expected_cols)
700+
701+
>>> ssoft_socca = build_the_ssoft(
702+
... aggregated_filename=aggregated_filename,
703+
... nparts=1,
704+
... nmin=50,
705+
... frac=None,
706+
... model='SOCCA',
707+
... version=None,
708+
... ephem_method="rest",
709+
... sb_method="fastnifty")
710+
<BLANKLINE>
711+
>>> assert len(ssoft_socca) == 2, ssoft_socca
712+
>>> assert "period" in ssoft_socca.columns, ssoft_socca.columns
677713
"""
678714
spark = SparkSession.builder.getOrCreate()
679715
spark.sparkContext.setLogLevel("WARN")

0 commit comments

Comments
 (0)