GalSim-developers/GalSim — pull request #242 by rmandelb (2012-09-07 01:51:24+00:00)
Title: #216 - lensing engine

This pull request is for the "lensing engine", which has two classes that can be used to generate shears:

(1) ```NFWHalo``` - to simulate weak lensing by a dark matter halo following a spherical NFW profile, useful for cluster lensing studies
(2) ```PowerSpectrum``` - to draw random shears that obey some specified power spectrum

The classes work in similar ways: you initialize the class with basic information about the NFW halo or the power spectrum. Then you use methods like ```getShear``` to get the shears at various positions.

(1) was contributed by @pmelchior ; for (2), @joezuntz provided python codes to carry out the core functionality, and I made the interface to GalSim plus unit tests etc. These are not complete in the way we'd want them for the final GalSim release, but enough machinery is in place that I think we should get it onto the master branch and have people play with it before doing further improvements. As usual, others' reactions might guide us as to which improvements are most important!

Looking at (1), the ```NFWHalo```, I see two possible improvements for the near term:
* Should perhaps update the python docstring to specify the halo mass definition - @pmelchior , if you are up for this, then you could update it on this branch while people are reviewing the code.
* I suspect that Mike's interface to config might require ```getShear``` to take single positions rather than full arrays (caching various calculations so subsequent calls are quicker). @rmjarvis , if that's correct, then please let me know and I could probably make that change over the weekend.

Looking at (2), there are several limitations that should eventually be fixed: it only works on grids (there was some code to get shears at non-gridded positions, but I simply didn't have time to test it, so I took it out for now and will add it later); there's a single, arbitrary convention for the normalization of the shears with respect to the input power spectrum (more on that later); we cannot get convergence, only shear; there is a limited set of power spectra and ways to specify them.

And of course, a limitation of both (1) and (2) is that they aren't plugged into the config scheme on #148, but Mike is planning to remedy that once they are merged into master.

So, here are the remaining to-do items (Peter might want to add some for the ```NFWHalo``` as these are mostly for ```PowerSpectrum```)... and the first 3 of them are likely worth discussing here in case we want to make simple changes now, and the last 3 are definitely a substantial amount of work for later:

* We have to decide if the normalization convention I've chosen for the power spectra, as in the dox, makes sense. It's one of only several possible conventions. We could eventually decide to support several conventions if there's demand for it (i.e., choosing via keywords) but at minimum we might discuss now if the one I've adopted is acceptable for the near future. Another issue related to conventions is the sign for g2... I found that if I feed these outputs into the GREAT10 code, for example, then I had to flip the signs of our g2 to get sensible results.

* Currently there is no intelligent unit handling. The dox say that the units for the power spectrum and the grids should be consistent, meaning that for example if you give it P(k) where k has units of inverse radians, then the galaxy positions or grid positions should be in radians. We leave it to the user to handle all this. But we could imagine doing something smarter -- e.g., having the user specify units for power spectrum and for grids, and checking for consistency or even converting if it knows all the conversion factors (e.g., deg vs. radians can be converted trivially; pixels vs. radians cannot be converted without specifying the pixel size).

* Right now, the lensing engine includes a few simple power law power spectra, or the user can specify their own that takes a single argument, k. In future, I can think of a few things we will want: (a) to have a cosmological shear power spectrum for some cosmology and dN/dz; (b) to allow the power spectrum function to take additional kwargs; (c) to write a function that reads in a power spectrum tabulated at some k values from a file, and automatically interpolates it to get P(k) for arbitrary k.

* I want the shear fields from a power spectrum to also include the possibility of getting the convergence.

* I'd like to eliminate the scipy dependence of ```NFWHalo``` if it's not too painful to do so, and also, we should have further discussion on just how much cosmology code we want to include with GalSim.

* Shears at non-gridded positions. Per Mike's request, I will make sure it's possible to do this one position at a time from a saved shear field (in addition to doing it by supplying all positions at once).

Once this is merged, I'll open an issue for these improvements.