Delay calculation of SED until it is actually needed. #1245
Conversation
| flux_ratio = self._flux_ratio(wave) | ||
| return self._jac, self._offset, flux_ratio | ||
|
|
||
| def _shoot(self, photons, rng): |
There was a problem hiding this comment.
Does this get the Poisson statistics right? If the SED is red, then we want more red photons not weightier red photons, no?
There was a problem hiding this comment.
Hm. Good point I'm pretty sure this was just the normal ChromaticTransformation._shoot function stripped of the parts that weren't relevant. So I think that one also has the wrong noise properties. It's quite possible others do as well.
I'll see if I can think of a good unit test for this that we can run on all the chromatic objects to see which ones break the expected Poisson statistics. It might not be possible to get it perfect for non-separable chromatic profiles. But we can at least get it right for the separable ones.
In any case, the right implementation for this one is, as you say, to not rescale the fluxes, but just assign wavelengths the way WavelengthSampler does.
There was a problem hiding this comment.
So, I looked into this a bit, and this function only gets used for PSFs (i.e. dimensionless chromatic objects). The base separable galaxy object with an SED does use the WavelengthSampler already. But then the PSFs are applied as photon operators, which act on the photons as they come in. So if one of them has a wavelength dependent flux scaling (not sure how that would arise physically, but we allow it), then it scales the flux values of the photons as requested. I think this is probably right. It's just a weird edge case that probably doesn't happen much in practice.
There was a problem hiding this comment.
Sounds good. If this is only for PSFs (sounds right) then yeah, I don't see how this could ever even reasonably get triggered. Even chromatic PSFs ought to have flux=1 at every wavelength I think. Thanks for taking a look.
| self.wave_list = obj.wave_list | ||
|
|
||
| @property | ||
| def sed(self): |
There was a problem hiding this comment.
Shouldn't deconvolved_sed = 1/original_sed like the original? Seems like there may be a missing round-trip sort of unit test here.
| self.wave_list = obj.wave_list | ||
|
|
||
| @property | ||
| def sed(self): |
There was a problem hiding this comment.
Similar to above, isn't sqrt_sed = sqrt(original_sed) ?
There was a problem hiding this comment.
Good catch on these two. This was a sloppy copy/paste error. Clearly we're missing a unit test that checks this.
| @@ -435,7 +435,7 @@ def interpolated(self): return False | |||
| @property | |||
| def deinterpolated(self): return self | |||
| @property | |||
| def SED(self): | |||
| def sed(self): | |||
There was a problem hiding this comment.
I think we need a deprecated uppercase SED here.
There was a problem hiding this comment.
I guess so. I don't think users ever probably duck type GSObjects as a ChromaticObject. That's more of a convenience for us in implementation. But it doesn't hurt to add.
|
Actually, I have one more thing while we're here. I think the following really ought to work, but it doesn't: sed = galsim.SED(
galsim.LookupTable(
[0, 621, 622, 623, np.inf],
# [0, 621, 622, 623, 10000],
[0, 0, 1, 0, 0],
interpolant='linear'
),
wave_type='nm',
flux_type='fphotons'
)
bp = galsim.Bandpass("LSST_r.dat", 'nm')
print(sed.calculateFlux(bp)) # gives nanIf I swap the np.inf for 10000 though, it does seem to work. Something about |
…eir calculations.
rmjarvis
added
chromatic
I was looking through the profile Jim made recently of an imSim run, which inspired the
combine_wave_listspeed up in #1243. This PR continues in that vein with various optimizations, mostly related to the SEDs.combine_wave_listtwice whenever we do those calculations. Once internal to the SED object, and once with a separatewave_listattribute in the chromatic object. But AFAICT, there is no possible way for these two wave lists to get out of sync. The only difference between them is that the one in the ChromaticObject can sometimes have an initial 0 wave and a final np.inf wave. Other than that, they were always identical. Furthermore, those bounding values were never useful for anything. So I switched thewave_listin the ChromaticObject to be a property that just returnsself.SED.wave_list. This didn't require any other adjustments in the code or tests, so I'm pretty confident it's functionally equivalent and saves half the calls tocombine_wave_list.<frozen importlib._bootstrap>:389(parent)in the profile, and I finally figured out that these are whenimportstatements happen inside a function rather than at the top of a module. I guess it's when python checks to see if it needs to import anything and realizes it already has everything imported. It's not a huge amount of time, but it's wasteful. So I went through the whole code base and moved as many imports as possible out of functions to module scope. Some of these required a little care to avoid circular import patterns. But since it was a lot of lines changed with no real content, I did this directly on main. If anyone wants to look, it's commit 5388b0e. But I don't think there is anything really substantive there. Anyway, this also led to a non-trivial speed up for faint objects that don't require a lot of work, so tiny bits of overhead are noticeable.obj.SED. So I decided while I was working on all this to change it to lowercase. (The capital one works still, but it deprecated.) This matches our prevailing style in GalSim to use lowercase attributes, even for abbreviations that are otherwise usually capitalized. (e.g. wcs, psf)The relevant timing script is
devel/time_faint_bd_gals.py. On my laptop, the time to draw 1000 galaxies went from 1.48 seconds before these changes down to 0.47 seconds after. So a factor of 3 faster, which feels pretty good to me. The remaining tall (ish) poles are initializing the knots in c++, multiplying the sed by a scalar to get the final flux right, making the non-chromatic Transformation object, and initializing the SED objects. Here is the full profile now, fwiw(The
__get__call in the second line is our lazy_property decorator, so it's just farming out to various other things and then writing the result as an attribute.) Nothing jumps out as particularly egregious now.