GalSim-developers/GalSim — pull request #504 by jmeyers314 (2014-02-01 01:29:25+00:00)
Title: #467

Hi all. I think now is a good time to get a review for the chromaticity code I've been adding to GalSim. Here's a rough summary of what's new:

New files

- galsim/bandpass.py
- Creates a `Bandpass` filter from congruent arrays of wave, throughput. (I chose bandpass since filter is a reserved word in python)
- `Bandpass` object is callable, returns throughput as fn of wave using `galsim.LookupTable` interpolation.

- galsim/sed.py
- Create a spectral energy distribution (SED) object from congruent arrays of wavelength and either flambda, fnu, or fphoton.
- `SED` object is callable, returns photon density as fn of wavelength.
- You can add and multiply `SED`s.
- You can set normalization at a specific wavelength or integrated against a `Bandpass`.

- galsim/chromatic.py
- `ChromaticObject`
- abstract base class for wavelength dependent sb profiles.
- features a generic .draw() method, usually overridden by ChromaticConvolution.draw() in practice.
- syntax for .draw() now requires a `Bandpass` over which to integrate.
- .draw() can optionally accept a number of "image integrators" (midpoint, trapezoidal, Simpsons, etc.) and associated parameters too.
- derived classes must implement an `evaluateAtWavelength` method that takes a wavelength argument and returns a `GSObject`
- `Chromatic`
- Chromaticize a `GSObject` by attaching it to an `SED`, which dynamically adjusts its flux.
- `ChromaticSum`
- Add `ChromaticObject`s (and `GSObject`s) together, by drawing them one at a time.
- `ChromaticConvolution`:
- Convolve `ChromaticObject`s (and `GSObject`s) together
- This one has some fancy-ish optimization in order to minimize the number of required convolutions by integrating over wavelength first when possible.
- `ChromaticShiftAndDilate`:
- Implement wavelength-dependent shifts and dilations of a fiducial PSF.
- This one is good for differential chromatic refraction (DCR), wavelength-dependent seeing, and a wavelength-dependent diffraction limit.
- `ChromaticAtmosphere`:
- Subclass of `ChromaticShiftAndDilate` specifically for DCR and chromatic seeing.

- galsim/dcr.py
- calculate DCR

Changed files

For better or worse, my emacs is configured to remove trailing whitespace on saving, so there are a lot of lines marked in the diff of the files below that are only whitespace. Sorry about that. I'll explicitly mention the line numbers that have actually changed...

-galsim/base.py
- 119-121: added a couple attributes to `GSObject` so it will duck type as a `ChromaticObject`
- 128-141: commented out `__iadd__` and `_isub__`, since the new `Add` as a function paradigm can occasionally transform a `GSObject` into a `ChromaticObject`, which I think can't be implemented with `__iadd__`.
- 587: added `evaluateAtWavelength` so a `GSObject` can duck type into a `ChromaticObject`.

- galsim/compound.py
- 63: Rename `Add` class to `Sum`
- 36-61: Make function `Add`, which returns either a `Sum` object or `ChromaticSum` argument depending on its arguments.
- 128-156: Similarly `Convolve` -> `Convolution` and `Convolve` becomes a function.

-galsim/integ.py
- added functions to integrate functions that return `galsim.Image`s over wavelength.
- `midpoint_int_image`, `trapezoidal_int_image`, `simpsons_int_image` are straightforward generalizations of the 1d equivalents.
- globally_adaptive_GK_int_image uses some ideas out of a recent adaptive quadrature review article I found, but doesn't actually seem to be much more efficient than the above algorithms, at least for the tests I've run. I kept it in for now, though.

New tests

- tests/test_chromatic.py
- `test_draw_add_commutivity()`
- `ChromaticConvolution` assumes that drawing and adding commute in galsim, so check it.
- `test_ChromaticConvolution_InterpolatedImage()`
- `ChromaticConvolution` also assumes that it can store intermediate results into an `InterpolatedImage` without accuracy penalty. Check this!
- `test_chromatic_add()`
- make sure the `__add__` operator is working as expected.
- `test_dcr_moments()`
- Test that moments of images drawn using different SEDs as input to `ChromaticAtmosphere` reproduce expected analytic results.
- `test_chromatic_seeing_moments()`
- Same as `test_dcr_moments` but for chromatic seeing.
- `test_monochromatic_filter`
- Check that `ChromaticObject.draw()` drawn through a very narrow filter reproduces equivalent `GSObject.draw()`
- `test_chromatic_flux`
- Check that the total drawn flux matches the expected normalization of a `ChromaticObject`.

- tests/test_sed.py
- `test_SED_add`
- Check that `SED`s add as expected

- tests/compare_image_integrators.py
- I wanted to compare the performance of different image integrators. They all seem to be about the same, converging somewhere around 100-1000 image evaluations, depending on your desired accuracy. (I'm sure this is SED/bandpass dependent, but those are the numbers for this particular setup).

- examples/demo12.py
- Three demonstrations in here. Use `Chromatic` to:
1) Draw a g-band dropout galaxy through LSST filters
2) Draw a bulge+disk galaxy with color gradient
3) Draw galaxy convolved with a chromatic atmospheric PSF

Finally, I have a few more ideas for galsim chromaticity that I haven't attempted yet. The first one is to cache the intermediate `InterpolatedImage`s produced by `ChromaticConvolution` in a dictionary. For my own studies, I know this can lead to orders of magnitude improvement in speed for things like ring tests, where you might generate many galaxy images all with the same SED and PSF. You might then fit all of these images using a slightly different (but unchanging) PSF. So what should be 2 InterpolatedImage creations can easily become hundreds or more.

The tricky part is figuring out how to index the dictionary. My idea is to use hashlib to give GSObjects and SEDs an MD5 hash. The hash would need to be updated upon transformations like `applyShear` or `setRedshift`. These hashes could then be used as keys in a dictionary that belonging to `ChromaticConvolution`. I'm curious if anybody has any ideas that are less invasive though.

The other unsolved problem I have is finding good default parameters for image integrators. N=250 midpoint rule seems to work reasonably for the LSST filters I've been working with, but it's not obvious to me that this will generally be a good choice. Parameterizing some kind of error estimate might be better. Or maybe just evaluating at precisely the wavelengths used to create a `Bandpass` object would also work, leaving it up to the user to sample his or her bandpass appropriately.

-Josh