Add the first composite type: Sum #157
Conversation
piff/sumpsf.py
Outdated
| self.chisq_thresh = chisq_thresh | ||
| self.max_iter = max_iter | ||
| self.kwargs = { | ||
| # If components is a list, mark the number of components here for I/O purposed. |
There was a problem hiding this comment.
typo: purposed -> purposes
piff/sumpsf.py
Outdated
| else: | ||
| self._num = None | ||
|
|
||
| # else components are not yet built. This will be called again when that's done. |
There was a problem hiding this comment.
Is this comment meant to be here or above the else statement?
| default_init='zero' | ||
|
|
||
| # If any components are degenerate, mark the sum as degenerate. | ||
| self.degenerate_points = any([comp.degenerate_points for comp in self.components]) |
There was a problem hiding this comment.
Do basis interpolation schemes still work when other components in the Sum have non-basis interpolation methods?
Edit: hmm, after looking at single_iteration, I think I can answer my own question that each component still fits separately, so there's no "interaction" between interpolations. But I'll leave it up in case I'm wrong...
There was a problem hiding this comment.
That's correct. The composition is of full PSF objects, so they don't interfere with each other.
The only thing this assignment does is to not do outlier rejection on the first iteration if any of the components use basis interpolation. I'm not sure how important that is, but it seemed sensible.
piff/sumpsf.py
Outdated
|
|
||
| # Draw the current models for each component | ||
| current_models = [] | ||
| for k, comp in enumerate(self.components): |
There was a problem hiding this comment.
technically don't need k and enumerate here
| if kk != k: | ||
| new_image -= current_models[kk][i].image |
There was a problem hiding this comment.
To make sure I understand what happens on the second and later iterations through all the components: on the first iteration, the first model component (c1) isn't affected by later components (say, c2 and c3) because they're initialized to zero flux. But on later iterations, c1 will be fitting to the data with the best-fit c2 and c3 models from the previous iteration subtracted (and c2 will fit with c1 and c3 subtracted, and so on). Correct?
There was a problem hiding this comment.
That's right. We always subtract off the current best model of all other components, and fit the component being worked on to the residual.
I think this should be pretty stable when applying some tweak to a primary model that is pretty good at describing most of the pattern, but has a residual with a very different interpolation pattern. E.g. maybe we add on a PixelGrid with much larger extent, but no interpolation, so it just captures the average large-scale residuals around a standard basis-interpolated model over a smaller footprint. None of my unit tests tried this out, but that's what I mostly have in mind for how the SumPSF would get used in practice.
piff/sumpsf.py
Outdated
| # The fit components in new_stars are what we want, but the data parts are | ||
| # corrupted by subtracting the other components, which we don't want. | ||
| # Just copy over the fits. |
There was a problem hiding this comment.
Something about this wording was very confusing. Possibly-too-wordy suggestion:
new_stars now has the updated fit components, but the data parts are corrupted by subtracting the other components during fitting. Update the stars by copying the undisturbed data from the original stars and the updated fits from new_stars.
| # Note: There is no useful test to make of the component parameters. | ||
| #print('comp0: ',test_star.fit.params[0]) | ||
| #print('comp1: ',test_star.fit.params[1]) |
There was a problem hiding this comment.
What do you mean here?
There was a problem hiding this comment.
Usually I had been testing the fitted parameters to some known true values. But for this test, both models are PixelGrid, so the parameters aren't anything I know a priori. The only thing to test is to use the models to draw images of the total PSF model and see if it matches the data acceptably.
tests/test_sumpsf.py
Outdated
| print('image max diff = ',np.max(np.abs(test_image[b].array - test_star.image[b].array))) | ||
| np.testing.assert_allclose(test_image[b].array, test_star.image[b].array, atol=1.e-8) | ||
|
|
||
| # Repeat with both components being PixelGrid mdoels. |
There was a problem hiding this comment.
typo: mdoels -> models
|
@jmeyers314 @aaronroodman |
|
I'm game. Thoughts on what to use for the truth? Maybe a GalSim atmospheric screen and then use GP for the interpolation of Kolmogorov for the primary component. Would that have a usefully interesting residual to mop up? Or would I need to include some optical component to give it something that wouldn't be well described by Kolmogorov? |
|
I'm not sure. If you keep the simulated exposure time short enough I'm sure an atmsopheric screen true psf would have interesting residuals to mop up. You might be able to get away with something simpler though, e.g., truth = Moffat with quickly spatially varying beta maybe? |
rmjarvis
force-pushed
the
sum
branch
2 times, most recently
from
92cf9c5 to
03b0d10
Compare
|
Hi Josh, I finally got back to this after 8 months of letting your very good suggestion sit unaddressed. I added a test using Atmospheric + Optical screens for the PSF. The model is a Kolmogorov PSF with GP interpolation + a PixelGrid with just a Mean interpolation. If I don't include the PixelGrid part, here are the images of the residuals (made by the new unit test when running from main): When I include the PixelGrid as a second component in the Sum, here are the residuals: This seems pretty good to me and shows a more realistic use case for how the Sum PSF might be used in practice. |
This is finally the first composite type that I've been working up to. Sum is probably a little simpler than Convolve, so I wanted to start with that. But the other two I have in mind are Convolve and Transform, similar to the corresponding GalSIm object types.
A Sum PSF is a model made from adding together multiple other PSF types. Each component can have its own model and interpolation. The parameters are stored in a list, so each component has its own numpy.array of parameters in the star.fit objects, each of which can be whatever length is required.
The initialization defaults to hsm for the first component and zero for the others, but the user can override this as desired. I suspect for most use cases, the first component will describe most of the PSF patter, and the later components will fill in the bits that weren't well modeled earlier (e.g. the tree-ring residuals we've seen in DES), so I think this will be a reasonable default initialization. Certainly we don't want multiple components all using the hsm initialization, since the sum would then be N x the PSF, which wouldn't be a good initial guess for the components. I suspect that getting good initial guesses for various components will be a big part of making composite PSFs work well, so there may be more work needed eventually on this front.