Add fiberflat gradient correction option

[rongpu]

This PR adds an option in desi_autocalib_fiberflat to apply a linear gradient correction to the fiberflats with reference fiberflats from a specific night, e.g., --gradient-ref-night 20000101. Some details:

• It also saves the gradient vectors in the GRADX and GRADY header keywords.
• The gradient fit is performed on the full focal plane (5000 fibers) and independently for each camera. The fit should be fairly robust against outliers.
• It produces an error when the gradient edge-to-edge variation exceeds -+20% (the largest gradient in October and December is ~12%).
• While it also fits an intercept, it is not used and thus the mean fiberflat value is unchanged.

[sbailey]

Thanks for adding this feature, including clear code and nicely documented docstrings.

Overall this looks good, but please try to make a version of poly_fit_2d that does not use statsmodels. Although it is in our NERSC desiconda environment, I checked on its provenance and it comes as a dependency of seaborn, which is a "for your convenience" addition not a core dependency (it would be different if e.g. scipy requires statsmodels, since we require scipy).

FYI it looks like https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.polyfit.html supports 2D fitting by specifying deg as an array instead of an integer, with https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.polyval2d.html, as the matching 2D evaluator. I haven't tried those myself, we and might have to add our own iterative outlier rejection.

For the pipeline, I think we'll also need to add an alternate way to specify the reference night so that we don't end up in a dependency loop (night A needs night B to be processed first as the reference, but night B will crash downstream because it needs tiles from A, thus requiring iterative processing...). But I don't have a specific alternative to suggest yet, so this is just a heads up. If needed, that alternate method could also be added in a different PR when calling this correction is added to the pipeline dataflow. Heads up @akremin .

[rongpu]

@sbailey I switched from using statsmodels to numpy with 4-sigma clipping. The resulting changes in the fitting results are quite negligible.

[sbailey]

Since a single reference night appears to be sufficient for this (#2172 recommended 20240125), to avoid cross-night dependencies while processing I suggest we do the following:

• Put the daily/calibnight/20240125/fiberflatnight*.fits files into $DESI_SPECTRO_CALIB, and update all configurations after 20230801 to use that as the default fiberflat.
• Update desi_autocalib_fiberflat to have a boolean option --solve-gradient (or similar name), which by default would use desispec.calibfinder.CalibFinder to look up those default flats, or optionally allow --gradient-ref-night REFNIGHT to override that.

In the code this would look something like:

if args.solve_gradient or args.gradient_ref_night is not None:
    ref_fiberflats = {}
    if args.gradient_ref_night is None:
        #- find default fiberflats to use as reference
        from desispec.calibfinder import CalibFinder
        for spectro, fiberflat in fiberflats.items():
            cf = CalibFinder([fiberflat.header,])
            ref_fiberflats[spectro] = cf.find('FIBERFLAT')
    else:
        #- use fiberflats from given night as reference
        for spectro in fiberflats.keys():
            ref_filename = findfile('fiberflatnight', night=args.gradient_ref_night, camera="{}{}".format(args.arm, spectro))
            ref_fiberflats[spectro] = read_fiberflat(ref_filename)
    ...

Then the pipeline can optionally call desi_autocalib_fiberflat --solve-gradient ... and it will use the default fiberflats from $DESI_SPECTRO_CALIB without having to wait for another night to be processed first, but the --gradient-ref-night option still allows overriding that for testing with alternate reference nights.

@rongpu @akremin let's discuss further details in person

[rongpu]

The last commit added the --solve-gradient option and desispec.calibfinder.CalibFinder.

The code is not complete yet and the line ref_fiberflats[spectro] = cf.find('FIBERFLAT') will break because CalibFinder does not actually load the reference fiber flats. It does still work when --gradient-ref-night is specified.

[sbailey]

Ah, you found the typo in my "from memory" example. CalibFinder returns the file path, so you need to wrap that in read_fiberflat, e.g.

ref_fiberflats[spectro] = read_fiberflat (cf.find('FIBERFLAT'))

Please try that. For now it will load a very old fiberflat as reference, but it should be able to find something to use without crashing. I'll update $DESI_SPECTRO_CALIB with a copy of 20240125 as the default fiberflat for more recent nights.

[rongpu]

Thanks for clarifying. I fixed the typo. It should work with the --solve-gradient option now.

[sbailey]

Looks good, thanks. I have added the fiberflatnight-*-20240125.fits files to $DESI_SPECTRO_CALIB and updated the config yamls to use them for any time period after 20241001. I tested this on 20241015 using --solve-gradient to have it automatically select the reference to use, and then compared the results using a different reference night (20241101): Plots are median(fiberflat)/median(fiberflat_ref) per fiber.

A future pipeline-focused PR will add the hooks to call this from the standard pipeline for the nights identified in #2172.