GalSim-developers/GalSim — pull request #332 by rmandelb (2012-12-26 03:04:54+00:00)
Title: #331, #316, #310: a 3-in-1 pull request

I know it's a holiday and nobody will be able to look at this soon, but I wanted to write this up while it's fresh in my mind, so - first of all, sorry for the holiday e-mail spam.

This pull request is meant to address issues #331, #316, and #310, with all edits being done on branch #331. A summary of changes:

For #316: outputs of shear estimation from ```EstimateShearHSM``` are returned as e1/e2 or g1/g2 (depending on the method), not as a ```Shear``` object. The reason for this is that shear estimators are not constrained to lie in the physical range |e|<1, so the methods associated with the ```Shear``` class will not necessarily work. (Note: the results of moments estimation are constrained to |e|<1 so those are still stored as a ```Shear```.) This change required minor edits in the C++ and python, and the unit tests and demos were updated accordingly.

For #331, there were a few changes:
* Mike had made a version of PSFCorr.cpp that used the ```Image``` class and TMV, and I put that in place of the previous PSFCorr.cpp that relied on other structures and ways of handling memory. I also updated the .h file documentation to be consistent with the slight changes in PSFCorr.cpp.

* There was code in PSFCorr.cpp to handle masks for the objects being measured, but there was no way to actually access it from python. I enabled the masks by making them required in C++ and optional in python (where the python code creates one for you if you don't provide it, before calling the C++ routines). In addition to the obvious changes to the interface to these routines, this involved tracking down a somewhat subtle bug deep in the guts of the C++. Since the masks are optional in python there was no need to change unit tests or demos, though I did add some new unit tests to make sure the code behaves properly when given a mask. Note: the C++ also had a place to handle masks in the PSF model used for shear estimation, but I did not make a way to supply a mask for the PSF model, because I don't see it being regularly useful. In typical data we might have to mask out some pixels by a galaxy due to cosmic rays or neighboring objects, but usually the PSF model does not have such defects. If someone strongly disagrees I can enable this code too; since I just did this for the galaxy masks, it wouldn't be too hard.

* Optimization: Mike had pointed out that a lot of time was being spent in ```std::exp``` for calculating the weights that go into adaptive moment calculation. I implemented a simple scheme to avoid calculating the weight all the time, simply by checking whether the ```rho2``` parameter that goes into ```exp(-rho2/2)``` is less than some value. If it is, then the calculation proceeds, otherwise, the weight is set to zero without calling ```exp```. Based on some timing and accuracy tests (more details to follow in a comment after this, for those who are interested and/or masochists) I selected a cutoff value of 25 (i.e., 5-sigma) as being optimal. For typical image sizes, ratio of PSF to pixel scale, and ratio of galaxy to PSF size, this optimization speeds up moments measurement by factors of 3-4, or PSF correction by factors of 2-3. For images that have even more white space around the edges than is typical, the speedup can be even more dramatic. I think that this factor of ~3 is enough that it's worth getting this onto master, rather than waiting for other optimizations (tabulating exp?), but opinions are welcome. Also, for those who want to play with this, the optimization can easily be shut off or changed by commenting out / editing one line of code (a new define statement in PSFCorr.cpp).

For #310: I updated the dox for the draw method (and the quick reference guide) to explain issues with draw with respect to flux normalization. There's no particular reason to group that change with this other stuff, but I was bored, so I did it.

There are some scripts in devutils/ that I used for the timing tests of the optimization, and tests of the mask functionality. I figured they might be useful later so I left them in there.

As mentioned above, there are new unit tests for the mask functionality that I enabled in the moments/shear estimation. However, suggestions for more are welcome, as I'm not completely satisfied that what I have is thorough enough. (It does include some tests that would reveal things like the subtle bug I fixed today.)