cpp “abs” vs. “fabs” coding bug; frozen misalign. polarisation matrix

bin/fba_launch

@@ -31,6 +31,7 @@ from fiberassign.fba_launch_io import (
     mv_temp2final,
     copy_to_svn,
 )
+from fiberassign.utils import get_fba_use_fabs
 from argparse import ArgumentParser
 
 # AR allowed steps in fba_launch
@@ -88,6 +89,23 @@ def main():
         args.rundate = get_latest_rundate(log=log, step="rundate", start=start)
         log.info("{:.1f}s\trundate\tsetting args.rundate={}".format(time() - start, args.rundate))
 
+    # AR fba_use_fabs: if None, set according to what is in the data/cutoff-dates.yaml
+    if args.fba_use_fabs is None:
+        log.info("")
+        log.info("")
+        log.info("{:.1f}s\tfba_use_fabs\tTIMESTAMP={}".format(time() - start, Time.now().isot))
+        args.fba_use_fabs = get_fba_use_fabs(args.rundate)
+        log.info("{:.1f}s\tfba_use_fabs\tsetting args.fba_use_fabs={}".format(time() - start, args.fba_use_fabs))
+        if os.getenv("FIBERASSIGN_USE_FABS") is not None:
+            if os.getenv("FIBERASSIGN_USE_FABS") != str(args.fba_use_fabs):
+                msg = "already defined environment variable FIBERASSIGN_USE_FABS={} will be overwritten by args.fba_use_fabs={}".format(
+                    os.getenv("FIBERASSIGN_USE_FABS"),
+                    args.fba_use_fabs
+                )
+                log.warning("{:.1f}s\tfba_use_fabs\t{}".format(time() - start, msg))
+    os.environ["FIBERASSIGN_USE_FABS"] = str(args.fba_use_fabs)
+    log.info("{:.1f}s\tfba_use_fabs\tsetting $FIBERASSIGN_USE_FABS={}".format(time() - start, os.getenv("FIBERASSIGN_USE_FABS")))
+
     # AR mtltime
     # AR if possible: setting to the latest of:
     # AR - latest timestamp for the program
@@ -739,6 +757,12 @@ if __name__ == "__main__":
         action="store_true",
         default=False,
     )
+    parser.add_argument(
+        "--fba_use_fabs",
+        help="value to determine the cpp behavior, see PR470 (default: based on the rundate)",
+        type=str,
+        default=None
+    )
 
     args = parser.parse_args()
     log = Logger.get()
@@ -766,6 +790,9 @@ if __name__ == "__main__":
     # AR mtltime (2021-09-01): if set to None, now setting it inside main()
     # AR          so that it is recorded in the fiberassign-TILEID.log file
 
+    # AR fba_use_fabs (2025-04-xx): if set to None, setting it inside main()
+    # AR          so that it is recorded in the fiberassign-TILEID.log file
+
     # AR goaltype
     if args.goaltype is None:
         args.goaltype = args.program

doc/changes.rst

@@ -6,8 +6,10 @@ fiberassign change log
 5.7.3 (unreleased)
 ------------------
 
+* Fix C++ bug with use of std::abs and improve floating point reproducibility (PR `#470`_).
 * Remove ``DesiTest`` from setup.py and warn about other deprecated features (PR `#464`_).
 
+.. _`#470`: https://github.com/desihub/fiberassign/pull/470
 .. _`#464`: https://github.com/desihub/fiberassign/pull/464
 
 5.7.2 (2023-10-04)

py/fiberassign/data/cutoff-dates.yaml

@@ -1,3 +1,4 @@
 rundate:
     std_wd: "2021-12-01T19:00:00+00:00" # std_wd not counted in per-petal-std after that rundate
     etc_stuck: "2022-04-28T19:00:00+00:00" # etc fibers are treated as stuck fibers after this date
+    fba_use_fabs: {"2019-09-16T00:00:00+00:00": 0, "2025-05-12T19:00:00+00:00": 1} # to set the env. variable for the c++ code (see PR470)

py/fiberassign/fba_launch_io.py

@@ -683,6 +683,7 @@ def print_config_infos(
     Notes:
         20210928 : add DUST_DIR, as assign.merge_results_tile() requires it to populate EBV=0 values.
         20211109 : add SKYHEALPIXS_DIR
+        20250418 : add FIBERASSIGN_USE_FABS
     """
     # AR machine
     log.info(
@@ -711,6 +712,13 @@ def print_config_infos(
                 time() - start, step, required_env_var, os.getenv(required_env_var)
             )
         )
+    for env_var in ["FIBERASSIGN_USE_FABS"]:
+        log.info(
+            "{:.1f}s\t{}\t{}={}".format(
+                time() - start, step, env_var, os.getenv(env_var)
+            )
+        )
+
 
 def get_ledger_paths(
     survey,
@@ -1975,6 +1983,7 @@ def launch_onetile_fa(
             gfafn,
         ]
     opts += ["--lookup_sky_source", args.lookup_sky_source,]
+    opts += ["--fba_use_fabs", str(args.fba_use_fabs)]
     log.info(
         "{:.1f}s\t{}\ttileid={:06d}: running raw fiber assignment (run_assign_full) with opts={}".format(
             time() - start, step, tileid, " ; ".join(opts)
@@ -2090,6 +2099,7 @@ def update_fiberassign_header(
         20211119 : added lookup_sky_source keyword
         20211227 : use newly defined args.goaltype (instead of args.program previously)
         20221031 : keywords too2, too3, etc could be automatically added (see get_desitarget_paths())
+        20250418 : add fba_fabs to track the FIBERASSIGN_USE_FABS quantity (--args.fba_use_fabs)
     """
     # AR sanity check on faflavor
     if faflavor != "{}{}".format(hdr_survey, hdr_faprgrm):
@@ -2169,6 +2179,8 @@ def update_fiberassign_header(
     )
     # AR lookup_sky_source
     fd["PRIMARY"].write_key("LKSKYSRC", args.lookup_sky_source)
+    # AR fba_use_fabs
+    fd["PRIMARY"].write_key("USE_FABS", str(args.fba_use_fabs))
     fd.close()
 
 

py/fiberassign/hardware.py

@@ -458,7 +458,7 @@ def load_hardware_args(focalplane=None, rundate=None, add_margins={}):
     return aa, time_lo, time_hi
 
 def radec2xy(hw, tile_ra, tile_dec, tile_obstime, tile_obstheta, tile_obsha,
-             ra, dec, use_cs5, threads=0):
+             ra, dec, use_cs5, threads=0, use_hardcoded_polmis_rotmat=True):
     '''
     For the tile pointed at (tilera, tiledec), project the (ra, dec)
     value into X/Y mm.
@@ -485,22 +485,38 @@ def radec2xy(hw, tile_ra, tile_dec, tile_obstime, tile_obstheta, tile_obsha,
     #y = np.array([y for x,y in xy])
     from astropy.time import Time
     from desimeter.fiberassign import fiberassign_radec2xy_cs5, fiberassign_radec2xy_flat
+    from inspect import getfullargspec
+
+    # Check for the precomputed polar misalignment option in desispec.
+    precomputed_matrix = 'use_hardcoded_polmis_rotmat' in getfullargspec(fiberassign_radec2xy_cs5).args
+
     # Note that MJD is only used for precession, so no need for
     # high precision.
     t = Time(tile_obstime, format='isot')
     mjd = t.mjd
 
     # Don't pass adc[12]: Let desimeter use its pm-alike routines
-    if use_cs5:
-        x, y = fiberassign_radec2xy_cs5(ra, dec, tile_ra, tile_dec, mjd,
-                                        tile_obsha, tile_obstheta)
+    if precomputed_matrix:
+        if use_cs5:
+            x, y = fiberassign_radec2xy_cs5(ra, dec, tile_ra, tile_dec, mjd,
+                                            tile_obsha, tile_obstheta, use_hardcoded_polmis_rotmat=use_hardcoded_polmis_rotmat)
+        else:
+            x, y = fiberassign_radec2xy_flat(ra, dec, tile_ra, tile_dec, mjd,
+                                             tile_obsha, tile_obstheta, use_hardcoded_polmis_rotmat=use_hardcoded_polmis_rotmat)
     else:
-        x, y = fiberassign_radec2xy_flat(ra, dec, tile_ra, tile_dec, mjd,
-                                         tile_obsha, tile_obstheta)
+        log = Logger.get()
+        log.warning('fiberassign_radec2xy_(cs5|flat) without use_hardcoded_polmis_rotmat option is deprecated')
+
+        if use_cs5:
+            x, y = fiberassign_radec2xy_cs5(ra, dec, tile_ra, tile_dec, mjd,
+                                            tile_obsha, tile_obstheta)
+        else:
+            x, y = fiberassign_radec2xy_flat(ra, dec, tile_ra, tile_dec, mjd,
+                                             tile_obsha, tile_obstheta)
     return x,y
 
 def xy2radec(hw, tile_ra, tile_dec, tile_obstime, tile_obstheta, tile_obsha,
-             x, y, use_cs5, threads=0):
+             x, y, use_cs5, threads=0, use_hardcoded_polmis_rotmat=True):
     '''
     For the tile pointed at (tilera, tiledec), compute the RA,Dec
     pointing of the specified X/Y location in millimeters.
@@ -525,15 +541,31 @@ def xy2radec(hw, tile_ra, tile_dec, tile_obstime, tile_obstheta, tile_obsha,
     # ra  = np.array([r for r,d in radec])
     # dec = np.array([d for r,d in radec])
     from desimeter.fiberassign import fiberassign_cs5_xy2radec, fiberassign_flat_xy2radec
+    from inspect import getfullargspec
     from astropy.time import Time
+
+    # Check if we can use the precomputed matrix option.
+    precomputed_matrix = 'use_hardcoded_polmis_rotmat' in getfullargspec(fiberassign_cs5_xy2radec).args
+
     t = Time(tile_obstime, format='isot')
     mjd = t.mjd
-    if use_cs5:
-        ra,dec = fiberassign_cs5_xy2radec(x, y, tile_ra, tile_dec, mjd,
-                                          tile_obsha, tile_obstheta)
+    if precomputed_matrix:
+        if use_cs5:
+            ra,dec = fiberassign_cs5_xy2radec(x, y, tile_ra, tile_dec, mjd,
+                                              tile_obsha, tile_obstheta, use_hardcoded_polmis_rotmat=use_hardcoded_polmis_rotmat)
+        else:
+            ra,dec = fiberassign_flat_xy2radec(x, y, tile_ra, tile_dec, mjd,
+                                               tile_obsha, tile_obstheta, use_hardcoded_polmis_rotmat=use_hardcoded_polmis_rotmat)
     else:
-        ra,dec = fiberassign_flat_xy2radec(x, y, tile_ra, tile_dec, mjd,
-                                           tile_obsha, tile_obstheta)
+        log = Logger.get()
+        log.warning('fiberassign_(cs5|flat)_xy2radec without use_hardcoded_polmis_rotmat option is deprecated')
+
+        if use_cs5:
+            ra,dec = fiberassign_cs5_xy2radec(x, y, tile_ra, tile_dec, mjd,
+                                              tile_obsha, tile_obstheta)
+        else:
+            ra,dec = fiberassign_flat_xy2radec(x, y, tile_ra, tile_dec, mjd,
+                                               tile_obsha, tile_obstheta)
     return ra,dec
 
 def xy2cs5(x, y):

py/fiberassign/scripts/assign.py

@@ -14,7 +14,7 @@
 import argparse
 import re
 
-from ..utils import GlobalTimers, Logger
+from ..utils import GlobalTimers, Logger, get_fba_use_fabs
 
 from ..hardware import load_hardware, get_default_exclusion_margins
 
@@ -197,6 +197,12 @@ def parse_assign(optlist=None):
                         choices=["ls", "gaia"],
                         help="Source for the look-up table for sky positions for stuck fibers:"
                         " 'ls': uses $SKYBRICKS_DIR; 'gaia': uses $SKYHEALPIXS_DIR (default=ls)")
+    parser.add_argument(
+        "--fba_use_fabs",
+        help="value to determine the cpp behavior, see PR470 (default: based on the rundate)",
+        type=str,
+        default=None
+    )
 
     args = None
     if optlist is None:
@@ -207,6 +213,9 @@ def parse_assign(optlist=None):
     if args.sky is None:
         args.sky = list()
 
+    if args.fba_use_fabs is None:
+        args.fba_use_fabs = get_fba_use_fabs(args.rundate)
+
     # If any of the masks are strings, determine the survey type from
     # the first target file to know which bitmask to use
     if isinstance(args.sciencemask, str) or \
@@ -296,6 +305,21 @@ def run_assign_init(args, plate_radec=True):
             for default_main_stdmask().
     """
     log = Logger.get()
+
+    # AR fba_use_fabs
+    if args.fba_use_fabs is not None:
+        if os.getenv("FIBERASSIGN_USE_FABS") is not None:
+            if os.getenv("FIBERASSIGN_USE_FABS") != str(args.fba_use_fabs):
+                msg = "already defined environment variable FIBERASSIGN_USE_FABS={} will be overwritten by args.fba_use_fabs={}".format(
+                    os.getenv("FIBERASSIGN_USE_FABS"),
+                    args.fba_use_fabs
+                )
+                log.warning(msg)
+        os.environ["FIBERASSIGN_USE_FABS"] = str(args.fba_use_fabs)
+
+    for env_var in ["FIBERASSIGN_USE_FABS"]:
+        log.info("run_assign_init(): use {}={}".format(env_var, os.getenv(env_var)))
+
     # Read hardware properties
     hw = load_hardware(rundate=args.rundate, add_margins=args.margins)
 
@@ -378,6 +402,22 @@ def run_assign_full(args, plate_radec=True):
         None
 
     """
+    log = Logger.get()
+
+    # AR fba_use_fabs
+    if args.fba_use_fabs is not None:
+        if os.getenv("FIBERASSIGN_USE_FABS") is not None:
+            if os.getenv("FIBERASSIGN_USE_FABS") != str(args.fba_use_fabs):
+                msg = "already defined environment variable FIBERASSIGN_USE_FABS={} will be overwritten by args.fba_use_fabs={}".format(
+                    os.getenv("FIBERASSIGN_USE_FABS"),
+                    args.fba_use_fabs
+                )
+                log.warning(msg)
+        os.environ["FIBERASSIGN_USE_FABS"] = str(args.fba_use_fabs)
+
+    for env_var in ["FIBERASSIGN_USE_FABS"]:
+        log.info("run_assign_full(): use {}={}".format(env_var, os.getenv(env_var)))
+
     gt = GlobalTimers.get()
     gt.start("run_assign_full calculation")
 
@@ -467,6 +507,22 @@ def run_assign_bytile(args):
         None
 
     """
+    log = Logger.get()
+
+    # AR fba_use_fabs
+    if args.fba_use_fabs is not None:
+        if os.getenv("FIBERASSIGN_USE_FABS") is not None:
+            if os.getenv("FIBERASSIGN_USE_FABS") != str(args.fba_use_fabs):
+                msg = "already defined environment variable FIBERASSIGN_USE_FABS={} will be overwritten by args.fba_use_fabs={}".format(
+                    os.getenv("FIBERASSIGN_USE_FABS"),
+                    args.fba_use_fabs
+                )
+                log.warning(msg)
+        os.environ["FIBERASSIGN_USE_FABS"] = str(args.fba_use_fabs)
+
+    for env_var in ["FIBERASSIGN_USE_FABS"]:
+        log.info("run_assign_bytile(): use {}={}".format(env_var, os.getenv(env_var)))
+
     gt = GlobalTimers.get()
     gt.start("run_assign_bytile calculation")
 

py/fiberassign/utils.py

@@ -122,7 +122,51 @@ def get_date_cutoff(datetype, cutoff_case):
     date_cutoff = config[datetype][cutoff_case]
     return date_cutoff
 
-
+
+def get_fba_use_fabs(rundate):
+    """
+    Return the value to which set the $FIBERASSIGN_ALGORITHM_EPOCH environment variable,
+        which drives some way the cpp computation is done.
+
+    Args:
+        rundate: rundate, in the "YYYY-MM-DDThh:mm:ss+00:00" formatting (string)
+
+    Returns:
+        fba_use_fabs: an integer (int)
+
+    Notes:
+        See this PR https://github.com/desihub/fiberassign/pull/470.
+        In the fiberassign running function, then one will set:
+            os.environ["FIBERASSIGN_USE_FABS"] = fba_use_fabs
+        So far:
+        - fba_use_fabs=0 means we reproduce the buggy behavior;
+        - fba_use_fabs=1 means we use the expected behavior.
+    """
+    # AR get the cutoff dates, and corresponding values
+    mydict = get_date_cutoff("rundate", "fba_use_fabs")
+    cutoff_rundates = np.array([key for key in mydict])
+    cutoff_mjds = np.array([get_mjd(_) for _ in cutoff_rundates])
+    values = np.array([mydict[_] for _ in cutoff_rundates])
+    # AR safe, order by increasing mjd
+    ii = cutoff_mjds.argsort()
+    cutoff_rundates, cutoff_mjds, values = cutoff_rundates[ii], cutoff_mjds[ii], values[ii]
+    # AR now pick the earliest date after the input rundate
+    # AR as we have set the first cutoff date as 2019-09-16,
+    # AR    there will always have some index returned here
+    # AR    though still checking if someone queries with an earlier rundate...
+    input_mjd = get_mjd(rundate)
+    if input_mjd < cutoff_mjds[0]:
+        msg = "rundate={} is earlier than DESI! ({})".format(
+            rundate, cutoff_rundates[0]
+        )
+        log.error(msg)
+        raise ValueError(msg)
+
+    i = np.where(cutoff_mjds <= get_mjd(rundate))[0][-1]
+
+    return values[i]
+
+
 def get_svn_version(svn_dir):
     """
     Gets the SVN revision number of an SVN folder.