clean up shearing conventions

[rmandelb]

As discussed in a few issues (e.g., #71) there are some ambiguities in the definitions of shears. For example, the following code snippet:

x = galsim.Shear(0.1,0)
print x.getG()

results in a shear value of 0.0501256289338. That's because in the first line, the 0.1 is being interpreted as a distortion (a^2-b^2)/(a^2+b^2). This could easily get confusing, so it should be revisited after the 2nd milestone.

We might also revisit some of the naming decisions we've made, e.g. applyShear/createSheared.

[rmjarvis]

To add to the confusion, I'll point out that applyShear explicitly uses the (g1,g2) definition rather than the (e1,e2) definition, which I think is what we want "shear" to mean. But that should be part of the conversation here.

[rmandelb]

Yes, good point, that's a worrisome inconsistency.

[rmjarvis]

When I wrote the code for applyShear, I didn't realize that the Shear class was wrapped, so I thought I was just changing the convention between python and c++ (which is worrisome enough, but I thought manageable for the time being). But since you pointed out that the Shear class is wrapped, it's definitely a problem.

[reikonakajima]

Just a short note: I was looking through Doxygen notes on galsim.psfcorr.EstimateShearHSM (from issue #151), and the otherwise crystal-clear example given has the shear/distortion convention issue that confused me for a bit. But I guess the issue is already raised here.

[rmandelb]

What about the new version that I just pushed (to master)? It's more explicit about what the conventions are.

On May 15, 2012, at 4:42 PM, Reiko Nakajima wrote:

Just a short note: I was looking through Doxygen notes on galsim.psfcorr.EstimateShearHSM (from issue #151), and the otherwise crystal-clear example given has the shear/distortion convention issue that confused me for a bit. But I guess the issue is already raised here.

---

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

---

Rachel Mandelbaum
http://www.astro.princeton.edu/~rmandelb
rmandelb@astro.princeton.edu

[rmandelb]

This issue on shearing conventions is one that I would like to get cleared up as part of this 2-week cleanup period between milestones. So, I'm going to ping some of the people who have been actively coding or discussing recently - @barnabytprowe @rmjarvis @TallJimbo @gbernstein @joergdietrich @pmelchior @joezuntz (apologies to those I might have left off accidentally).

I'll start by summarizing the problem and then proposing a way to address it. Let's discuss in this issue, and once we have a plan we can decide who is implementing it and when. (Or, please feel free to let me know if you think this is not worth dealing with now...!)

• the C++ SBProfile code has a Shear class with several constructors. The constructor that takes 2 doubles as arguments actually takes distortions (e1, e2) where |e| = (1-q^2)/(1+q^2); these are not shears |g| = (1-q)/(1+q) or conformal shears eta (though once you have a Shear object you can use the various getG and getEta methods to pull those bits of information out, or you can reset the value using setG1G2 etc.). This class is wrapped, so you can use it in python, e.g.,

x = galsim.Shear(0.1, 0.0)
print x.getG()

would give a value of 0.0501.

• we also should think about the Ellipse class and its initialization
• the C++ SBProfile member function shear takes e1/e2 as arguments. This method is also wrapped and accessible in python.
• the python methods that shear our GalSim base classes (Gaussian, Moffat, etc.), applyShear and createSheared, use g1/g2 rather than e1/e2.
• the HSMShapeData class has members like observed_shape, where shape means distortion. (I think this is okay in general but maybe not in a situation where we already have shear/distortion-related confusion going on! But it might be worth giving this class member functions to get out shears...)

So as a consequence, if we work with our base classes and do this:

gauss1 = galsim.Gaussian(sigma = 1.)
gauss1s = gauss1.createSheared(0.1, 0.0)

then the equivalent operation in terms of SBProfiles would be

gauss2 = galsim.SBGaussian(sigma = 1.)
gauss2s = gauss2.shear(math.tanh(2*math.atanh(0.1)), 0.0)

In a word, yuck.

Are there any other aspects to the problem that I'm missing?

Given that statement of the problem, I think that we have a few possible approaches:

1. to pick a default quantity to work in terms of (presumably shear or distortion, not conformal shear), then make the Shear constructor take that quantity in its constructor, and make all the C++ and python shearing methods work in terms of it as well. If we choose shear then we have to change the C++, if we choose distortion then we have to change python.

2. to pick a default quantity as in (1), but allow an optional keyword that would change to use the other. i.e., if we choose shear g, then

gauss1 = galsim.Gaussian(sigma = 1.0)
gauss1s = gauss1.createSheared(0.1, 0.0)

would apply a shear g1=0.1, g2=0.0, whereas

gauss1s = gauss1.createSheared(0.1, 0.0, distortion = True)

would correspond to specifying a distortion of (0.1, 0.0). (We would do this at the C++ level as well, presumably.)

3. as for radii, we could imagine requiring that the convention be specified explicitly. In other words, x = galsim.Shear(0.1, 0.0) or gauss1s = gauss1.createSheared(0.1, 0.0) would raise an exception, but x = galsim.Shear(g1 = 0.1, g2 = 0.0) or x = galsim.Shear(e1 = 0.1, e2 = 0.0) is acceptable. Or, we could make the syntax x = galsim.Shear(0.1, 0.0, type = "e") (i.e., require type = "e", "g", or "eta"). We could decide if we want to make this change in both C++ and python, or as for the radius specification, if we want to enforce this in the boost python wrapping.

Are there other options I'm missing?
My preference is for (2), though I am vacillating between that and (3). Also I am thinking mostly as a developer... from the perspective of user-friendliness, I have worries about all the options, really!

[rmjarvis]

I think my preference is to use the word "shear" to always mean g1,g2 and the word "distortion" to always mean e1,e2. Then we can have parallel methods like obj.applyShear(g1,g2) and obj.applyDistortion(e1,e2). Then the user can always use their preferred version with a minimum of fuss.

The only problem with this is that there is already an applyDistortion which takes an Ellipse object, rather than e1,e2. So I would propose changing that to something like applyEllipse or applyTransform. Or probably even remove it entirely, since I don't think it's something we need. At least not until we start using the general AffineTransform stuff eventually. But we can put off the appropriate notation to use that until then.

On the C++ side, I'd have the current Shear class be renamed to something like GenShear, which would be a base class for two derived classes Shear and Distortion (or three if we want one for conformal shear), which mostly just change the meaning of the constructor. Everything else (e.g. the getG1G2 methods) would be done via GenShear. In fact, we could just let Distortion be the base class rather than GenShear, and have Shear be a derived class from that. Either way.

[TallJimbo]

The approach I took in LSST-land was to have different ellipse,
ellipticity, and radius classes for each parameterization, with common base
classes and implicit conversions between them. It's a lot of overhead, but
it makes a pretty nice interface for the user, I think.

You can find the headers here:

http://dev.lsstcorp.org/cgit/LSST/DMS/afw.git/tree/include/lsst/afw/geom/ellipses

I'd do some things differently if I were to do it again, of course, and I'd
be curious what a Mike Jarvis review and rewrite of those classes would
look like, but it might be a useful starting point.
On May 16, 2012 4:12 PM, "Rachel Mandelbaum" <
reply@reply.github.com>
wrote:

This issue on shearing conventions is one that I would like to get cleared
up as part of this 2-week cleanup period between milestones. So, I'm going
to ping some of the people who have been actively coding or discussing
recently - @barnabytprowe @rmjarvis @TallJimbo @gbernstein @joergdietrich
@pmelchior @joezuntz (apologies to those I might have left off
accidentally).

I'll start by summarizing the problem and then proposing a way to address
it. Let's discuss in this issue, and once we have a plan we can decide who
is implementing it and when. (Or, please feel free to let me know if you
think this is not worth dealing with now...!)

• the C++ SBProfile code has a Shear class with several
  constructors. The constructor that takes 2 doubles as arguments actually
  takes distortions (e1, e2) where |e| = (1-q^2)/(1+q^2); these are not
  shears |g| = (1-q)/(1+q) or conformal shears eta (though once you have a
  Shear object you can use the various getG and getEta
  methods to pull those bits of information out, or you can reset the value
  using setG1G2 etc.). This class is wrapped, so you can use it in
  python, e.g.,

x = galsim.Shear(0.1, 0.0)
print x.getG()

would give a value of 0.0501.

• we also should think about the Ellipse class and its initialization
• the C++ SBProfile member function shear takes e1/e2 as arguments.
  This method is also wrapped and accessible in python.
• the python methods that shear our GalSim base classes (Gaussian,
  Moffat, etc.), applyShear and createSheared, use g1/g2
  rather than e1/e2.
• the HSMShapeData class has members like observed_shape,
  where shape means distortion. (I think this is okay in general but maybe
  not in a situation where we already have shear/distortion-related confusion
  going on! But it might be worth giving this class member functions to get
  out shears...)

So as a consequence, if we work with our base classes and do this:

gauss1 = galsim.Gaussian(sigma = 1.)
gauss1s = gauss1.createSheared(0.1, 0.0)

then the equivalent operation in terms of SBProfiles would be

gauss2 = galsim.SBGaussian(sigma = 1.)
gauss2s = gauss2.shear(math.tanh(2*math.atanh(0.1)), 0.0)

In a word, yuck.

Are there any other aspects to the problem that I'm missing?

Given that statement of the problem, I think that we have a few possible
approaches:

1. to pick a default quantity to work in terms of (presumably shear or
   distortion, not conformal shear), then make the Shear constructor
   take that quantity in its constructor, and make all the C++ and python
   shearing methods work in terms of it as well. If we choose shear then we
   have to change the C++, if we choose distortion then we have to change
   python.

2. to pick a default quantity as in (1), but allow an optional keyword
   that would change to use the other. i.e., if we choose shear g, then

gauss1 = galsim.Gaussian(sigma = 1.0)
gauss1s = gauss1.createSheared(0.1, 0.0)

would apply a shear g1=0.1, g2=0.0, whereas

gauss1s = gauss1.createSheared(0.1, 0.0, distortion = True)

would correspond to specifying a distortion of (0.1, 0.0). (We would do
this at the C++ level as well, presumably.)

3. as for radii, we could imagine requiring that the convention be
   specified explicitly. In other words, x = galsim.Shear(0.1, 0.0) or
   gauss1s = gauss1.createSheared(0.1, 0.0) would raise an exception,
   but x = galsim.Shear(g1 = 0.1, g2 = 0.0) or x = galsim.Shear(e1 = 0.1, e2 = 0.0) is acceptable. Or, we could make the syntax x = galsim.Shear(0.1, 0.0, type = "e") (i.e., require type = "e", "g", or
   "eta"). We could decide if we want to make this change in both C++ and
   python, or as for the radius specification, if we want to enforce this in
   the boost python wrapping.

Are there other options I'm missing?
My preference is for (2), though I am vacillating between that and (3).
Also I am thinking mostly as a developer... from the perspective of
user-friendliness, I have worries about all the options, really!

---

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

[pmelchior]

I'm with Mike: "shear" needs to mean (g1, g2). Anything else is just confusing. If applyDistortion is not needed, I'd suggest to remove it, and reserve applyTransform for more general transformations.

[rmandelb]

Okay, first of all: I realize now that I expressed a preference for (2) but I didn't say which quantity I would like as default. I want shear to mean (g1, g2), so it seems like so far we're all in agreement on that point.

I need to think some more about these suggestions -- more tonight / tomorrow...

[rmandelb]

Also, I notice that the documentation on applyDistortion in galsim/base.py is a little incomplete in that it suggests that it's for making an e1/e2 distortion. In reality, the SBProfile distort method that it calls can be used to apply an Ellipse transformation that represents a distortion in the e1/e2 sense and/or a magnification and/or a shift. If we remove applyDistortion from the python level, then we still have applyShift and applyShear but we'd be missing something that can do magnification, unless I'm missing something. So applyDistortion does have some additional functionality here -- it could perhaps be changed to applyTransform as Peter suggests.

[pmelchior]

I actually did not mean to suggest that, but as long as applyDistortion has a purpose, I would keep it. I think applyTransform is a bit opaque a name for it.

[rmandelb]

Really sorry for misrepresenting your opinion, I misread the end of your second sentence.

I do think that having applyDistortion as a name for something that can be used for magnifying or shifting is a bit odd. So I want to keep that functionality but without a name that sounds like it should just be doing an (e1, e2)-type transformation.

[rmjarvis]

How about adding applyDilation for this functionality. (Better IMO than applyMagnification since that is ambiguous about referring to size or flux or both.)

[barnabytprowe]

To enter the fray, I definitely agree that Shear should be (g1, g2) and refer to (a - b)(a + b) ellipticities.

I think it makes sense to also adopt eta = conformal shear as per B&J 2002 and this leaves (e1, e2) up for grabs by the commonly adopted |e| = (a^2 - b^2)/(a^2 + b^2) -type ellipticity definition. This latter I also think should be called "Distortion" (or perhaps "Polarization"), although what exact name we choose is less important than that we document clearly and precisely what we mean by each term.

If we go for Distortion as the name for these (e1, e2) ellipticities, then we need to steer well clear of talking about distortions when we mean the full Affine transformation, so I would back talking about applyTransform type name conventions.

In terms of 1), 2) or 3) in the choices you mention above Rachel, I'm afraid I'm also somewhat undecided between 2 and 3. I think I prefer 3...

(Jim, while I found it very interesting and instructive to look at I can't help but think that the approach you took in LSST would be a lot of overhead here, but I'm willing to be outvoted on that! Oh, I think the docs for ReducedShear.h actually refers to the ConformalShear while I look at it...!)

[gbernstein]

My opinion on the shear conventions is that you should have a consistent default that people can enter as a pair of numbers, but a means to use another and have it be clear that you're doing so (Rachel's #2 basically). Having required classes for angles, radii, shears, etc makes the usage more ornate and may deter new users from just jumping in and trying the code to do their simple tasks first.

Jim, I wrote classes for spherical coordinates that work like the scheme you describe, that I'm very happy with - all transforms between any two coordinate systems are done transparently by assigning coordinates from one derived class to another. So it's elegant. But I also think it's overkill for this job (there are only 3 systems of shear that are likely to be used) and erects another level of complication for future people who might want to learn to use & alter the code. And I say this as someone who usually errs on the side of wasting too much time trying to be elegant !

[reikonakajima]

Rachel: the update on the Doxygen text (as seen in the diff) is now clear.
(Minor problem: I ran Doxygen on the master node, but I don't see the new text reflected on the html pages, though I can tell that it has changed from before. If someone could tell me if I'm doing something wrong, please let me know.)

From a "new user's" perspective: I think Gary's note on "having a consistent default that people can enter as a pair of numbers" is the way to go, and it would make sense that that be the (g1, g2) shear. It would be great to have an option of specifying (e1,e1), (eta1,eta2) and/or shifts and dilations, via Mike's functions-with-explicit-function-names. Although the few examples I've seen so far have been, really, amazingly crystal-clear.

[TallJimbo]

I'm happy to concede that a lot of new classes may be overkill for this purpose.

Barney, thanks for catching the documentation bug in the LSST code; you're absolutely right.

[rmandelb]

In comparison with the cases of radii and angles, both of which we've dealt with recently:

a) We dealt with different radius specifications by preserving the SBProfile C++ almost completely, modifying the wrappers in pysrc/ and the base classes instead.

b) We dealt with angles by making an Angle class. My feeling at the time was that this made things super clear for developers but I do have a lingering worry that users might find it cumbersome.

We can't take the analogy with those cases too far, and here we already have a Shear class defined for us with lots of methods for handling conversions and so on. My gut reaction is to be leery of introducing several new classes, so I think I still prefer my case (2) above, because

• it gives a consistent definition of shearing by default, so that when we intend to shear by (g1, g2) as we often will, the function calls will be very simple, i.e. applyShear(g1, g2) without additional keywords etc.
• it allows for deviations from the default definition for those who prefer other definitions, by simply using a keyword parameter

However, I'm undecided as to how to implement this. That is, we could imagine completely preserving the C++ as is and introducing all changes at the level of the wrappers and the various python base classes and their methods, just as for radii. But - since there is a Shear class in C++ and shearing is so central to what we are trying to do, I kind of feel that here we do want consistency between the C++ and the python in terms of shearing conventions, which would dictate changes in the Shear constructor, the shear method of SBProfiles, distort, and possibly other things (basically I would want to read through all the .h files carefully to find other places this might crop up!).

As for applyDistortion, I'd be in favor of replacing with applyDilation (as Mike suggested) and/or applyMagnification. The rest of the applyDistortion functionality is already there in the applyShift and applyShear. Actually, there is actually room for both applyDilation and applyMagnification, as we'd probably want to reserve the word magnification for the lensing use of the term, i.e. at fixed surface brightness which doesn't preserve flux, whereas dilation could be at fixed flux. Hmmm.

[rmjarvis]

I mostly agree with the growing consensus, except that I still think that

obj.applyDistortion(e1,e2)

is clearer than

obj.applyShear(e1,e2,distortion=True)

As for the C++ side, the easiest implementation would be to just change all "distort" to "transform" and "shear" to "distortion". Then if we want, we can add a Shear class that derives from Distortion that would have the regular shear definition for its constructor.

[rmandelb]

Mike, I agree with that in the context of this issue, but in a more general context I think applyDistortion could be confused with instrumental effects. If we have a separate method for these two different conventions for specifying a shape (which I would be amenable to) then I think we need a different name.

Likewise, I don't want to change shear to distortion in C++ for that reason.

On May 17, 2012, at 12:16 AM, Mike Jarvis wrote:

I mostly agree with the growing consensus, except that I still think that

obj.applyDistortion(e1,e2)

is clearer than

obj.applyShear(e1,e2,distortion=True)

As for the C++ side, the easiest implementation would be to just change all "distort" to "transform" and "shear" to "distortion". Then if we want, we can add a Shear class that derives from Distortion that would have the regular shear definition for its constructor.

---

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

---

Rachel Mandelbaum
http://www.astro.princeton.edu/~rmandelb
rmandelb@astro.princeton.edu

[rmandelb]

I get the sense that there's a reasonable consensus for what we should do, with the main sticking point being names of functions and so on. I usually think of the quantity (a^2-b^2)/(a^2+b^2) as either ellipticity or distortion. Distortion seems like a bad choice for the reasons already mentioned above (too general a term, could mean a few things) and ellipticity seems like a bad choice because there are a gazillion different ellipticity definitions (okay, not a gazillion, but at least 3). I do think it is okay/non-confusing to take the current distort and make it into transform. Any further comments on that point would be welcome.

Unless someone else particularly wants to take this on, then I am going to volunteer myself because I'm interested in the work and I think it's a reasonable fit for my coding skills. Due to timing constraints (moving in 5 weeks!) I can't promise to finish it all next week before starting the next milestone, but I don't think it's a big problem if this issue takes a little bit longer than that.

[gbernstein]

On May 17, 2012, at 12:01 AM, Rachel Mandelbaum wrote:

However, I'm undecided as to how to implement this. That is, we could imagine completely preserving the C++ as is and introducing all changes at the level of the wrappers and the various python base classes and their methods, just as for radii. But - since there is a Shear class in C++ and shearing is so central to what we are trying to do, I kind of feel that here we do want consistency between the C++ and the python in terms of shearing conventions, which would dictate changes in the Shear constructor, the shear method of SBProfiles, distort, and possibly other things (basically I would want to read through all the .h files carefully to find other places this might crop up!).

It would not be difficult to change the C++ Shear class to treat g1/g2 as its default language rather than e1/e2. And I think it's a good idea to have the same convention in all parts of the package. The somewhat harder part would be finding all the code that uses it and making sure it's consistent with the new convention, that's not too bad either though.

If it gets to the point where you want to make this change and no one else has eagerly volunteered to learn how Shear class works and change it, then I could do this work in a few hours.

[rmandelb]

Gary - I'm quite happy to do it. I agree that it's not tricky, just a bit of tediousness in making sure to find all the relevant bits (and for me, there's an additional advantage in that it forces me to go over bits of the code, some of which I'm still learning). I think that the only subtleties left are what to call some of the methods - e.g., shear should be changed to apply a shear (not distortion), but distort might need a more general name like transform (since it does more than just an e1/e2 distortion), and on the python side, we might want an applyShear (same as now) but also applySomething that does an e1/e2 distortion, except that if we call it applyDistortion it could be confused with some kind of instrumental / astrometric distortion, which would be bad.

There's also Mike's idea which I like of having a class for e1/e2 Distortion (or actually, a GenShear base class from which both Shear and Distortion are derived). But again, there's that tricky bit about the name Distortion. Not sure what to do about that.

[rmjarvis]

Given the ambiguity of the word distortion, and nobody jumping in with a name they prefer, I think it's worth taking another look at Rachel's option 3 (or maybe a modification thereof).

We could just have one set of functions for shearing, so no applyDistortion or createDistorted. Likewise, just the one Shear class. Then the arguments to these functions (or Shear's constructor) would always be named arguments defining which version we mean: applyShear(g1=0.03, g2=0.05), applyShear(e1=0.2), applyShear(eta=0.1, theta = 12.6 * galsim.degrees), etc.

If we want, we could have the default meaning be (g1,g2), but I don't mind requiring that they always be specified either.

Then all the different functions that take a shear would just pass the (args, kwargs) to the Shear class constructor for parsing, so that would only need to be done once. And maybe if something takes a shear and also other things (do we have any of those? If not, we might someday), we could require that the argument be an explicit Shear object, so you'd have to write foo(Shear(e=0.3,theta=th), image=im1, ...) rather than intermixing the g1,g2,e1,e2,eta,theta, etc. options among other arguments.

In C++, we'd have to implement this a bit differently of course, and we can discuss what would make the most sense. But I think on the python side, this would be fairly clear in all cases.

[rmandelb]

Then the arguments to these functions (or Shear's constructor) would always be named arguments defining which version we mean<<

So when it comes to the Shear class, I presume you mean for this to be done at the level of the wrappers in pysrc/, similarly to what Jim did with the radii - i.e., requiring that either g1/g2 or e1/e2 or eta1/eta2 or whatever be set such that there is a single uniquely-specified shear g1/g2, converting the inputs to that shear, and then create the Shear object using that g1/g2? (this would obviously requiring changing the Shear class such that the constructor that takes 2 doubles is interpreting them as g1/g2 rather than e1/e2 -- but in principle I think no other changes are needed to the C++ Shear class) And I agree that we should have g1/g2 be the default.

As for the applyShear / createSheared, I think this is fine. Likewise we can make applyDilation etc. And I think we've agreed that distort method / applyDistortion needs to go.

I guess the only remaining question is then the C++ shear method, but similar to the Shear class, we could make it have the rather sensible default be that it takes g1/g2, and do any other magic in the wrappers (in any case, in the python we'll mostly be using our base classes and applyShear etc.).

I like this. Let's leave a final decision for the telecon Tues, but I think it's a reasonable compromise.