Implement fast centroid fit to sky camera spots

[dkirkby]

The centroid fit implemented by @HenriCoquinot in PR #13 significantly improves the sky background estimates but requires 8-9s (on perlmutter) to fit each new SkyCam exposure, compared with ~120ms without the centroid fit. The nominal SkyCam cadence is 60s but the GFA cadence is 5s, so processing a new SkyCam exposure disrupts the processing of GFA images. This PR is to re-implement the centroid fit to run in less than 1s while obtaining essentially the same results.

[dkirkby]

The current implementation consists of:

• foreach camera (2)
  • foreach spot (10)
    • initialize (dx,dy) = (0,0)
    • for each fit iteration (10)
      • calculate shifted spot profile
      • calculate chisq and its gradients, then perform a simple gradient descent step to update (dx,dy)
  • calculate mean (dx,dy) over all (unmasked) spots
  • calculate shifted spot profiles for all fibers
  • fit for (flux,bg) using final mean-shifted spot profiles

I suspect most of the time is spent calculated shifted profiles 2 x (10 x 10 + 10) = 220 times, but I need to confirm this.

[dkirkby]

My estimate (270) of the number of times shifted profiles are calculated was low because shifts are also used to estimate partial derivatives and there is some branching logic so the number of calls varies each time. It looks like the average number of calls is around 731.

Each call to desietc.util.shifted_profile takes 10.7ms (also on perlmutter) so this accounts for 7.8s of the total 8-9s, confirming that this is the bottleneck.

[dkirkby]

The basic outline of a faster implementation is to:

• precompute each profile on a uniform grid of (dx,dy) shifts in the SkyCamera constructor
• for each new exposure
  • compute the chisquare of all observed spots for each (dx,dy) on the uniform grid (using the precomputed shifted profiles)
  • estimate the (dx,dy) where the chisquare is maximum
  • estimate the flux and bg level at this (dx,dy)

Note that chisquare values are computed assuming a global (dx,dy) shift for all spots, instead of first optimizing chisquare separately for each fiber, then using their average for a final fit.

With this approach, the best fit model at each (dx,dy) on the grid can be computed in one step with just linear algebra using the existing desietc.util.fit_spots method.

[dkirkby]

The SkyCamera.setraw() method now has a new option fast_centroids that enables the new code. I am currently testing it on a full night's data (20240222) to compare the timing and results with the original (slower) implementation.

Specifically, there are 3 modes to compare (w/o temp corrections), depending on the options passed to SkyCamera.setraw():

• original mode, currently in use, with refit=True, fit_centroids=False
• slow centroid fits with refit=False, fit_centroids=True, fast_centroids=False
• fast centroid fits with refit=False, fit_centroids=True, fast_centroids=True

The fast implementation takes ~1.0s per pair of SkyCam exposures, based on initial tests, compared with 8-9s for the original implementation.

[dkirkby]

Sky levels from each pair of the 421 pairs of SkyCam exposures taken during 20240222 calculated with the slow vs fast implementation of centroid fits:

There is a small (-0.1%) systematic shift to lower fluxes between the methods, and some outliers with up to 1% deviations at fluxes below ~40, but the overall agreement looks good.

To go further, will need to analyze more nights and compare with offline sky estimates...