radius definitions for various SBProfiles about galsim HOT 31 CLOSED

This is a good point to raise, and is related to a discussion we had in #33 about the way to define the characteristic angular scale for the optical PSF component. Gary's suggested \lambda / D parameter requires a light rewrite of the SBAiry interface to make the C++ layer consistent with the Python layer. Totally agree that this is not urgent, but consistency across classes is something we probably want - not everyone reads documentation as closely as we might hope!

from galsim.

So here's what I'm thinking (would appreciate comments, @barnabytprowe , @rmjarvis , @TallJimbo , @gbernstein ):

We might want multiple ways to define a Gaussian, with a sigma or a FWHM. We might want multiple ways to define an exponential, such as using a disk scale length (as is currently done) or a half-light radius (which is how it's done for Sersic but not Exponential). I don't anticipate making multiple radius definitions for Moffat, deVauc, Sersic, or Airy unless somebody requests it. So there are two basic questions here:

1. at what level should we do this? Should we be mucking around with C++ SBGaussian etc., or is this a job to be sorted out with the python base classes Gaussian etc.? I'm inclined towards the latter.

2. if we want to have multiple ways of defining profiles, then it means multiple keywords. And then we have to handle it properly if someone does something weird like trying to create a Gaussian using both a sigma and a FWHM with inconsistent numerical values. What's the best way to do that? Should it throw an exception if someone tries to do that? I'm inclined to think that's the right solution.

While we're at it, I am inclined to make the keywords more descriptive, like halfLightRadius and scaleLength instead of re and r0, as Mike requested in #85. (Would do this for all of the base classes, but again, just at the python level.)

I don't intend to handle the lambda/D issue for Airy since I'm not sure what is the best way to handle it, unless someone chimes in on that point.

from galsim.

Multiple keywords seem the most appropriate choice to me as well, with an exception raised when multiple radii are passed. We will have to choose one radius to be the default for when the argument is passed positionally, and I'd recommend we make that consistent by using half-light radius.

I think we should focus on doing this at the Python level - there's no way to do this at the pure-C++ level, really - but we can also try to add it to the Python wrappers for the SBProfile classes. That would require a little effort, but I don't think it would be that hard (and we should start with trying this, because then we could make use of it in the high-level pure-Python classes). I'm happy to do this myself, as it may requiring a little extra Boost.Python special sauce.

from galsim.

I agree about choosing a default for when it is passed positionally, and that this should be the half-light radius (except for Gaussian, for which I thought the two choices would be sigma and FWHM, with sigma being the default).

I hadn't thought about doing this at the level of the python wrappers for SBProfile classes... possibly because I was planning on doing this myself, and working at the level of the pure-Python classes is easier for a newbie like me. I don't entirely understand the advantage of doing this at the level of the wrappers. It seems to me that it's just shifting around the location of the bits of code to handle the different choices of keywords. Do you think it's better to keep the pure-python classes written more cleanly/simply (which motivates having these details handled by the wrappers), or is there some other advantage that I am missing?

I agree that if it's done at the level of python wrappers for SBProfile classes, then you're a better person to do it than I am!

from galsim.

My only motivation for doing it at the SBProfile wrapper level is to keep the GalSim and SBProfile interfaces from diverging in a place where they're otherwise very similar. But I think you're right, too - it's nice to keep all the code that specific to a particular profile in one place (i.e. with the SBProfile object) to extent that's reasonable.

I'll give it a try. I would appreciate input on exactly what the keyword names should be.

from galsim.

Another good reason, so please go for it (am assigning to you now).

from galsim.

I've added multiple-radius support to SBGaussian (half_light_radius, sigma, fwhm) and SBExponential (half_light_radius, scale). I've also renamed the parameters on SBSersic and SBDeVaucouleur to half_light_radius.

What do we want to do about SBMoffat and SBAiry? I assume we might want a FWHM parameter for both of these, and I think that's straightforward for SBMoffat, but I understand we may want to do SBAiry differently. Do we also want half-light radius for these?

from galsim.

I don't know of any simple way to figure out the half-light radius for an Airy (with variable obscuration) or for a Moffat (with arbitrary beta & truncation). You'd need to do the numerical integral & solution for half-light radius when creating each object if it were desirable to specify by half-light radius. On the other hand, for a given telescope, you will know lambda/D, not the Airy half-light radius, when you need to construct the object. So while the consistency of specifying everything with half-light radii would be nice, it would be a bother and not obviously useful.

On Apr 12, 2012, at 5:10 PM, Jim Bosch wrote:

I've added multiple-radius support to SBGaussian (half_light_radius, sigma, fwhm) and SBExponential (half_light_radius, scale). I've also renamed the parameters on SBSersic and SBDeVaucouleur to half_light_radius.

What do we want to do about SBMoffat and SBAiry? I assume we might want a FWHM parameter for both of these, and I think that's straightforward for SBMoffat, but I understand we may want to do SBAiry differently. Do we also want half-light radius for these?

---

Reply to this email directly or view it on GitHub:
#65 (comment)

from galsim.

I agree with Gary. It's nice to have a half-light radius for Sersic profiles, which are often specified in such terms, but I'm not used to thinking about an Airy or Moffat in terms of a hlr.

from galsim.

It looks like we can leave the SBAiry class as it is, then - the input argument is called D_, but the documentation suggests that it's actually D/lambda.

SBMoffat is actually already implemented in terms of a half-light radius, and it does the numerical integration that Gary mentioned. But it looks like it might need to do that numerical integration to normalize the flux anyhow, and I don't completely follow how the truncationFWHM plays into the Fourier-space evaluation parameters.

Gary, I'd be more comfortable if you could look into adding a separate constructor from scale radius or FWHM to SBMoffat in C++ (feel free to work on this branch). I can then wrap that with the existing half-light radius constructor to give it the same multi-radius Python signature the other profiles have. I'm worried I'll mess up the Fourier-space evaluation if I do it all myself.

from galsim.

Jim, I just pushed to branch #65 after changing the SBMoffat constructor. Since C++ constructors are distinguished only by argument signature, and all the sizes would be doubles, I created an extra constructor argument that tells it what kind of size you are specifying. This argument will take one of SBMoffat's enum values:
SBMoffat::FWHM
SBMoffat::HALF_LIGHT_RADIUS
SBMoffat::SCALE_RADIUS
with the 2nd one being the default, so that omitting the argument matches previous behavior.

I have verified that this compiles and links with the existing example code, and that SBDraw.py does what it should for the default case. But I can't do testing of the other options because the python wrapping and other steps needed for a simple test program are beyond my current capabilities. So it'll have to be up to you to make sure it works properly and let me know if there are problems.

On Apr 23, 2012, at 10:52 AM, Jim Bosch wrote:

It looks like we can leave the SBAiry class as it is, then - the input argument is called D_, but the documentation suggests that it's actually D/lambda.

SBMoffat is actually already implemented in terms of a half-light radius, and it does the numerical integration that Gary mentioned. But it looks like it might need to do that numerical integration to normalize the flux anyhow, and I don't completely follow how the truncationFWHM plays into the Fourier-space evaluation parameters.

Gary, I'd be more comfortable if you could look into adding a separate constructor from scale radius or FWHM to SBMoffat in C++ (feel free to work on this branch). I can then wrap that with the existing half-light radius constructor to give it the same multi-radius Python signature the other profiles have. I'm worried I'll mess up the Fourier-space evaluation if I do it all myself.

---

Reply to this email directly or view it on GitHub:
#65 (comment)

from galsim.

Hi Jim - it occurs to me that for this milestone, when implementing multi-object functionality, we might want the ability to specify radii in different ways. I see based on the above commits that you've already done quite a lot here, including unit tests. Is there significant work left to do on this branch? If you don't see yourself having time for this in the next week, is it something that you could imagine handing off to someone else so we can get it merged into master early next week-ish?

from galsim.

All that remains are "real" unit tests that verify that the magic numbers I've used to scale the various radii are correct. I would be happy to hand it off to someone else; I think it's just a matter of making images that are large enough that the half-light radii and FWHM values can be computed somewhat accurately numerically for comparison.

from galsim.

Per discussion in #103, I will take over the unit tests (unless you're back and are itching to do this, Jim..).

from galsim.

But, as a sanity check: rather than numerically integrating light profiles to compare half-light radii and FWHM, could I not simply make two SBProfiles using different keywords, with your magic numbers determining the transformation between the input parameters, and check to make sure that the resulting images are the same?

from galsim.

I am back, but still with very limited great3 time, so I'm happy to pass this off to you.

Your proposal for unit tests sounds fine if you write them without looking at my magic numbers and instead get them from some other source; I'd also like to make sure I didn't do the (admittedly simple) math wrong.

from galsim.

That's probably not sufficient, since if Jim's magic numbers were wrong, you still could get the images to match, since your test would potentially make the same error that Jim (hypothetically) did. e.g. the FWHM isn't guaranteed to be the actual FWHM of the resulting image. So I agree with Jim that making a very high-resolution image where you can explicitly test the FWHM, half-light radius, etc. is a good test to have.

from galsim.

Good point. So as a specific example, I could do my_gaussian = galsim.Exponential(scale = 1.0) and make a high-resolution image. Then I will integrate numerically and make sure that the scale radius is really 1.0. I will also check what is the HLR, and will make another exponential with my_gaussian_2 = galsim.Exponential(half_light_radius = blah), draw it, and make sure it is the same as the first image.

I will make different tests for the different profiles depending on which keywords were defined.

To do the integration, I was planning to cheat and just take a 1d slice through the image centroid to make a vector of values, then integrate via direct summation for total enclosed flux within r is sum(2 pi r I(r) dr). My assumption is that this should be good enough for a high-resolution image.

from galsim.

Cheating is always dangerous for a test...

My rule of thumb is that tests should always do things the dumb slow way with absolutely no optimizations, but that are as obvious as possible to check that they are correct by inspection. So I'd recommend doing the 2-d integral.

from galsim.

Fair enough. But still via summation?

If I can't do it via summation, then I'm not sure how - I'm not familiar with python integration packages (and since what Jim did was at the level of the python wrapping, I can't do the tests in C++ and use the integ package). It would probably take more time than I have this week, given the other issues I'm committed to, to sort that out.

from galsim.

I'd definitely go for direct summation (with a sufficiently hi-res image). Much easier, much more transparent, and for a sufficiently hi-res images I think it's safe to approximate those trapeziums as rectangles...

from galsim.

Yes. Definitely direct summation, and then just figure out how accurate you expect the value to be and make sure your tests only require agreement to that precision.

from galsim.

Okay, one more thing for this inexperienced python-er: If I use summation to create a table of r vs. I(<r), and want to interpolate to find, say, the HLR, then I'll need some way to interpolate in python. Recommendations?

from galsim.

I think for half-light radius, I would test this in the opposite direction. You know what the half-light radius is supposed to be. So do I(<r_HLR) and see if it is approximately equal to 0.5 * I_tot.

from galsim.

The same for FWHM. You have a value to test. So find I(FWHM/2) and see if it is 0.5 * I_peak.

from galsim.

Yes, good point.

from galsim.

Okay... I just pushed a first attempt at this for Gaussian only, since I thought it would be useful to have feedback before implementing very similar tests for Moffat, Sersic, Exponential. Here's the deal:

• The absolute clunkiest / simplest way to evaluate a 2d image I(<r) and I(r) is only accurate at high resolution, but
• The actual resolution that is needed is not totally obvious to me, and depends sensitively on the type of profile.

My solution was to write the test such that it computes the quantities we want (i.e., I(r=sigma)/I(r=0), I(r=0.5*FWHM)/I(r=0), I(r<hlr)/I_total) for some fairly low resolution image, and then keeps increasing the resolution by factors of 2 until the results converge to within a fraction of a %. Then I compare the converged result with the ideal expectation at the 2nd decimal place. This is potentially dangerous as we could end up iterating a long time, and it results in some repeated code chunks. But, it does work for the Gaussian case, and I cannot see myself doing better in the ~day before we want to merge this in so it can be used for this milestone. Conclusion: either we do similar tests for Moffat, Sersic, Exponential -- or, we merge this in rather dangerously without unit tests and someone writes more sophisticated ones later (I could do a better job if I had more time) -- or, I fill in the rest of the tests this way since they DO work, and we open an issue to do a more sophisticated job later. Thoughts?

from galsim.

I think you should probably just make a new issue for the tests.

I think this approach may not work in a reasonable amount of time for some of the profiles that have extended wings, and I think it's more important to let other code using it move ahead.

(I also feel a little bad about dumping this issue on you, and I'd be happy to take it back after the deadline).

from galsim.

FWIW, the 19 tests in test_SBProfile.py (i.e., all the old tests + these new ones for Gaussian) take a total of 17s to run in my machine. I would definitely be worried about timing issues for deVauc, but we don't have a python base class for it (yet), so it seems like it wouldn't be horrible to put this in now for Moffat, Exponential, and Sersic with two relatively low values of n (e.g., 1.5, 2). At least that gives us some basic assurance that the magic numbers in pysrc/SBProfile.cpp work, and then we can merge it and open a new issue for more sophisticated tests?

On May 8, 2012, at 10:30 PM, Jim Bosch wrote:

I think you should probably just make a new issue for the tests.

I think this approach may not work in a reasonable amount of time for some of the profiles that have extended wings, and I think it's more important to let other code using it move ahead.

(I also feel a little bad about dumping this issue on you, and I'd be happy to take it back after the deadline).

---

Reply to this email directly or view it on GitHub:
#65 (comment)

Rachel Mandelbaum
http://www.astro.princeton.edu/~rmandelb
[email protected]

from galsim.

I think that sounds reasonable, but only if it doesn't take too long! (I'm totally bogged down /recently stuck in #101 and feel the walls closing in :) :) :)

from galsim.

Now that I did it for 3 radius types for Gaussians, it will take 5 minutes to do so for the other types. I'll do it and make sure it's not hideously slow, then make a pull request and open an issue for more sophisticated unit tests.

from galsim.