GalSim-developers/GalSim — pull request #232 by rmjarvis (2012-08-08 05:46:46+00:00)
Title: #200 Various optimizations

This pull request is basically just some speed improvements in the code. It includes the topics from issues #200 and #211, plus a couple of other improvements I noticed as well.

The upshot is that MultiObjectDemo script 4 is much faster for the OpticalPSF and Kolmogorov PSFs. The current output on my laptop is:
```
OpticalPSF: Total time = 1.511389 (fft: 0.661542, phot: 0.830254)
Kolmogorov * Airy: Total time = 1.155592 (fft: 0.206877, phot: 0.928913)
```
The output using the current master branch is:
```
OpticalPSF: Total time = 10.606787 (fft: 6.806380, phot: 3.780898)
Kolmogorov * Airy: Total time = 7.297645 (fft: 1.144141, phot: 6.134237)
```

Code changes:
- Added functions calculateStepK and calculateMaxK to the SBInterpolatedImage class. These will calculate the appropriate values of stepK and maxK based on the actual image provided rather than just basing them on the dimensions of the image. These are not done automatically, but if you think the provided image might be larger than was really needed (as is the case with OpticalPSF, since it doesn't know a priori how big it will need, so it goes bigger than is usually required), you can let the code update the stepK and/or maxK values, so subsequent draw operations can be faster.
- Use these functions for both OpticalPSF and AtmosphericPSF. It was more important for the former, but it helped with both.
- Tabulate the xvalues for Lanczos, since that is faster than doing the multiple sin functions. Also switched the uvalue table to include the actual values when fluxConserve=true rather than an intermediate value from which the final answer is calculated.
- Added a new SBKolmogorov class which implements the Kolmogorov spectrum in C++. Some advantages over the AtmosphericPSF class include
1. it does the FFT in 1-d using a Hankel transform rather than a 2-d fft,
2. direct calculation of the stepk and maxk using the normal definitions, so no need for the calls to calculateStepK or calculateMaxK.
3. drawing with fft (the normal way) uses the exact Kolmogorov spectrum rather than the fft of an fft, which had introduced some modest errors.
4. drawing with photon shooting is faster, since it knows the profile is radially symmetric.
- Added a Kolmogorov class in python that wraps SBKolmogorov. I didn't replace AtmosphericPSF, since the functionality isn't identical. I think the new one is more accurate, since it avoids one fft. But the fact is that images created by the two classes disagree at the few percent level. So for now, I left AtmosphericPSF in case people want to run some comparisons to confirm that the difference is acceptable. And looking forward, I thought the existing code might provide a starting point for a more realistic multi-layer model of the atmosphere as per issue #169.
- I calculated a more precise conversion between lambda/r0 and fwhm for a Kolmogorov spectrum using numerical integration in the C++ code.
- I also calculated the conversion between lambda/r0 and the half-light radius. So that is now implemented for Kolmogorov. (I didn't back port it to AtmosphericPSF, but the same number should apply.)
- Changed ProbabilityTree from being a trinary tree to a (faster) binary tree.
- Added a shortcut entry into the binary tree that speeds up the binary search in most cases.
- Added a specialization for Airy when obscuration == 0.0. (Probably a common case.)
- Made Table inherit from FluxDensity, which lets a Table be used directly as a radial function for OneDimensionalDeviate. (A feature used by the new SBKolmogorov class.)
- Removed some extraneous floor functions. The floor function is actually fairly slow, so it is worth removing when possible.

Test suite changes:
- Changed the spline test from using OpticalPSF to using in image based on sheared and shifted Gaussians. Barney already had to update this test when making a change to the OpticalPSF class, and I had to again. Since it seems plausible that this might happen again, I switched the spline test to use a simpler image based on 3 sheared gaussians added together. This can (probably) remain fixed even with future code improvements.
- Added tests of new Kolmogorov class analogous to other GSObject tests
- Removed the `sbprofile_` part of the name of many tests in test_SBProfile.py. They just seemed superfluous given the name of the file they are in.
- Changed tests of OpticalyPSF and AtmosphericPSF to match new stepk and maxk values.
- Changed some tests not in test_SBProfile.py to use GSObjects rather than direct SBProfiles. I think it makes sense to have SBProfile tests in test_SBProfile.py. But everywhere else, I think we should be using the GSObject versions, since they have a cleaner interface.

Finally a change in the script, which acts as a proposal for a change to OpticalPSF:
- In the code for script 4, I specify OpticalPSF to use bilinear interpolation rather than Lanczos interpolation. This makes a big difference for the OpticalPSF photon shooting times, since the Lanczos interpolant has a lot of negative flux, which is hard for the shooting method to deal with. So, should we keep this? Should we make bilinear the default action in OpticalPSF if the interpolant is not specified? Is there a compelling reason to keep Lanczos? If so (or even if not), should we even allow the interpolant to be specified by the user? Wouldn't it be better to choose the "correct" interpolant and just use that, rather than have interpolantxy be an option?