DM-55875: Implement Gomes 2025 with the help of Claude/Fable5. - #37
Conversation
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #37 +/- ##
==========================================
+ Coverage 89.13% 91.74% +2.60%
==========================================
Files 9 11 +2
Lines 939 1429 +490
==========================================
+ Hits 837 1311 +474
- Misses 102 118 +16 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
Going to merge this to try to make rubin-env 14. |
cmsaunders
left a comment
cmsaunders
left a comment
There was a problem hiding this comment.
Overall, I think the structure of the empirical_2pcf and GridConvolutionGP classes works well. They are reasonably easy to understand, at least with help from the technotes. How they are used in GPInterpolation is harder to follow, but whether you want to change anything depends on how much you are worried about outside users. I've made some comments about possible changes inline in gp_interp.py, or I still think you could consider a whole new derived class, which would help indicate that the empirical-2pcf method is a radical departure from the earlier methods.
A couple more points:
- I might have written this inline somewhere, but I'm still not sure: would using optimizer=
empirical-2pcfwith solve_method=directgive you the Gomes method or something mathematically identical, or is this its own method? In either case, should it be clarified that it's not recommended? - I have an overall question about the kernel arrays: in a bunch of places, it says that you have a square kernel, with a shape that is an even number of pixels, but the place with zero lag is at pixel N /2, and I think this means at the center of pixel N/2. What is the reason for having this off-center array, rather than having zero lag at the point where the four central pixels meet?
- I'll blame Claude, but the code would be more readable with fewer acronyms. GRF for Gaussian Random Field and PSD for Positive Semi-Definite are possible to figure out, but slow me down when I'm trying to comprehend what's going on in the code.
| measured 2-points correlation function (`Léget et al 2021 | ||
| <https://doi.org/10.1051/0004-6361/202140463>`_), the measured anisotropic 2D | ||
| 2-points correlation function can be used directly as the kernel, following | ||
| `Gomes et al 2025 <https://doi.org/10.3847/1538-3881/ae1a7b>`_. The measured |
There was a problem hiding this comment.
I think this part needs to mention that you have made algorithmic changes, and that this is not a straight adoption of exactly what is done in Gomes+.
|
|
||
| Input is expected to be 2-dimensional, i.e. X.shape = (n_samples, 2). | ||
|
|
||
| :param x_grid: Lag coordinates of the grid columns, zero lag |
There was a problem hiding this comment.
As discussed on Slack, the term "lag coordinates" is new to me. This would be a good place to define it.
| This is equivalent to the singular value clipping used by | ||
| Gomes et al. (2025). | ||
|
|
||
| Input is expected to be 2-dimensional, i.e. X.shape = (n_samples, 2). |
There was a problem hiding this comment.
This doesn't seem to match any of the inputs below.
| hyperparameters are fitted: the kernel argument must be left to | ||
| its default (an error is raised otherwise), and min_sep and | ||
| nbins are ignored (the grid is controlled by max_sep and | ||
| pixel_size). Conversely, an EmpiricalCorrelationKernel has no |
There was a problem hiding this comment.
This sentence starting with "Conversely" doesn't quite work for me. I would suggest something like "The corresponding kernel, EmpiricalCorrelationKernel, has no hyperparameters and can only be used with the "empirical-2pcf" or "none" optimizers; other options will raise an error."
| :param cg_rtol: Relative tolerance of the conjugate gradient solve of | ||
| solve_method="spectral". [default: 1e-7] | ||
| :param cg_maxiter: Maximum number of conjugate gradient iterations of | ||
| solve_method="spectral". [default: 500] |
There was a problem hiding this comment.
This feels like too many new optional arguments that are specific to 'empirical-2pcf'. Looking briefly at the __init__ for .empirical_2pcf(), it looks like these are just set, not used. Could you let the user modify GPInterpolation.optimizer if they want to set these, rather than expose them all here?
| fill_value=0.0, | ||
| ) | ||
|
|
||
| def _eval_lags(self, dx, dy): |
There was a problem hiding this comment.
Can you document what this is doing a little more?
| spreading (one fine-grid pixel wide, i.e. pixel_size / upsample) and | ||
| the PSD projection is done on the padded grid instead of on the | ||
| point-set covariance; both effects are second order and shrink with | ||
| upsample. |
There was a problem hiding this comment.
Can you rephrase this? I'm not sure what the intended meaning is, but "with upsample" doesn't make sense to me.
| at pixel N//2, whose first row and column (the most negative lag) | ||
| have no positive counterpart on the grid; the output has odd side | ||
| length N + 1, zero lag exactly at its center, and satisfies | ||
| out[c + i, c + j] == out[c - i, c - j]. This is the grid-space |
There was a problem hiding this comment.
Does c mean the central pixel?
| m = n * upsample | ||
| # Spectrum with the frequency origin at the center, frequencies | ||
| # running from -n/2 to n/2 - 1 in each axis. | ||
| f_shift = np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(xi))) |
There was a problem hiding this comment.
Why is there an inverse shift of the kernel before the fft?
|
|
||
|
|
||
| @timer | ||
| def test_empirical_2pcf_psd_repair(): |
There was a problem hiding this comment.
Remove? Or, looking at this more, would you still use this if you were using the "empirical-2pcf" optimizer with solve_method="direct"? I guess so.
2 participants
With the help of Claude (Fable5), this PR implementing the method described in Gomes et al. 2025. The idea is to have an empirical kernel that is derived directly from the measured 2D two-point correlation function and so get free of model choice for the kernel and hyper parameters fitting.