cpu - compare ex2d_patch with specter

[sbailey]

I pushed a branch "compare_ex2d_patch" that has a notebook that has the seeds of a comparison between gpu_specter.extract.cpu.ex2d_patch and specter.extract.ex2d.ex2d_patch. Please expand, debug, turn into units tests, etc.

In this case we don't necessarily expect exactly agreement because there could be reasonable differences in exactly how regularization or tiny eigenvalues are handled, but they should probably agree more than they currently do.

Assigning to @dmargala to explore. Let me know if you have trouble getting the pieces up and running either on your laptop or at NERSC.

[dmargala]

Thanks @sbailey, this looks much more elegant than the frankenstein test setup that I started working on to compare ex2d_patch.

[dmargala]

I added a test for ex2d_patch based on the notebook in the branch. The test comparing the output with specter passes when the ndecorr=True is passed to specter.extract.ex2d_patch. Does that seem like it should be right?

There are a few other differences but that seems to be the main one.

The next biggest gain occurs when I drop the weak flux=0 prior (by commenting out this line)

[sbailey]

Nice.

"ndecorr" means "noise decorrelation of neighboring fibers", which comes at the expense of correlating the signal, which we usually don't want. i.e. you can chose to decorrelate the signal or the noise, but not both, and we'd rather have correlated noise than correlated signal.

So we will want gpu_specter to implement the equivalent of what specter ex2d does when ndecorr=False (i.e. the opposite of what it does now). If we keep the option to do either, let's pick a clearer name since I always read "ndecorr" as meaning "no decorrellation" which is actually the opposite of what it does (noise decorrelation = do decorrelate the noise...).

Regarding regularization: we need some kind of regularization for 3 cases:

• if the PSF spot for some fiber,wavelength covers entirely masked pixels or is entirely off the edge of the CCD, it becomes unconstrained and the matrices become singular. A tiny amount of regularization fixes that and results in flux=0 error=huge which is what we want.
• if the PSF spot only touches one or two unmasked pixels, then it is mostly unconstrained and could have a crazy value.
• Imperfect PSF models can lead to ringing in the extracted spectra

Please start with implementing the equivalent of ndecorr=False, and then we can revisit exactly how and how much we do regularization. gpu_specter is probably too simplistic, and original specter is probably too messy.

for reference Bolton & Schlegel 2009 equations 17-19 is the math part about the (non)decorrelation of the neighboring noise, that gets implemented within specter.extract.ex2d.resolution_from_icov.

[dmargala]

Thanks for the explanation.

I'll start working on implementing the equivalent of the ndecorr=False case.

I'm still a little confused about the default behavior in specter. Can you confirm that the ndecorr (and decorr!) branching logic for the default case is correct in these two places:

• https://github.com/desihub/specter/blob/8b1ecb12dc04206809e34505a4426cf146d5814c/py/specter/extract/ex2d.py#L402-L405
• https://github.com/desihub/specter/blob/8b1ecb12dc04206809e34505a4426cf146d5814c/py/specter/extract/ex2d.py#L542-L553

[sbailey]

I also find the specter implementation to be confusing with name switches and whether the option is supposed to be True/False/None/Other. I think the specter branching logic is correct, albeit confusing.

ndecorr = noise decorrelation = default for resolution_from_icov = same as what gpu_specter currently does. In that case, ndecorr=True results in decorr=None because the code has switched to meaning decorrelate the signal (decorr) rather than decorrelate the noise (ndecorr). (!)

ndecorr = False = don't decorrelate the noise = do decorrelate the signal (good) = call resolution_from_icov with decorr = a list of how many wavelengths to include in each block (which is more generic than it needs to be, since every block is always the same size...).

I'll try to put together a simple example to document this more, but I'll submit this comment now in case it is enough for you to keep going.

[dmargala]

Thanks, that helps.

And just to be sure, I confirmed that the default behavior starting from desi_extract_spectra goes down the the ndecorr=False:

• https://github.com/desihub/desispec/blob/8e7ebe22ec6e1987860f1b66dd90953e11853984/py/desispec/scripts/extract.py#L411
• https://github.com/desihub/desispec/blob/8e7ebe22ec6e1987860f1b66dd90953e11853984/py/desispec/scripts/extract.py#L60

[dmargala]

I've implemented B&S equations 17-19 in ex2d_patch as the new default for ex2d_patch. I did add an option to switch to the other mode. Let me know if you have any thoughts about how that is implemented.

I also updated the test to compare the updated ex2d_patch to specter's ex2d_patch with ndecorr=False.

The outputs are now all passing an np.allclose comparison.