Implementation of elliptical profiles in lenstool

71 views
Skip to first unread message

Richard Massey

unread,
Aug 18, 2015, 5:52:10 AM8/18/15
to lens...@googlegroups.com
Hi Eric, JP, et al.

I’ve been trying to understand how lenstool implements elliptical mass profiles. If I read the source code right, there is a bug… that, amongst other things, might explain the incorrect magnification of the supernova in A2744.


Gravitational lensing maps the source plane onto the lens plane through a coordinate transformation whose Jacobian is the magnification. Lenstool then manufactures elliptical lenses from circular lenses (which have nice maths) by applying an extra transformation, in e_pot.c. Lenstool rotates the coordinate system so the lens's major axis points along x, then applies a shear via the mapping x->x/(1+e), y->y/(1-e). However, this shear changes the area of a circle: it has Jacobian (1+e^2).

The consequence of this extra Jacobian is probably be that inferred lens masses are wrong by a factor (1+e^2), as their mass normalisation is lowered to compensate. But it could be possible that the consequences become more sinister. If the mass is well-constrained by other means, perhaps the magnification is wrong. None of this will matter for (roughly circular) relaxed clusters, or the quasar lenses that were often modelled in the 1990s. But A2744 is highly elongated (e=0.61 for clump #2), so even an e^2 term could be significant.

As a fix, I would propose that the mapping in e_pot.c be changed to x->x(1+e^2)/(1+e), y->y(1+e^2)/(1-e). Unless this has already been done somewhere in the code that I haven’t spotted.


What do you think? Please tell me I’m wrong.
Thanks,
Richard


Richard Massey
Institute for Computational Cosmology,
Durham University,
South Road, Durham DH1 3LE, UK
Tel: +44 191 334 3652
Mob: +44 7740 648080
Web: http://www.dur.ac.uk/r.j.massey/



Joseph E. Coleman

unread,
Aug 18, 2015, 2:12:36 PM8/18/15
to lens...@googlegroups.com
Hi all,

I have a followup question regarding the topic of elliptical mass
distributions. In the lenstool documentation, the ellipticity is
defined as ϵ = (A^2 - B^2) / (A^2 + B^2), however, in the paper by
Hjorth and Kneib (97), they use the definition ϵ = (A - B) / (A + B).
Is this an issue?

Typically an ellipse with semimajor and semiminor axes A and B
respectively would be defined as

X^2 / A^2 + Y^2 / B^2 = 1.

However, using ellipticity ϵ, the resulting ellipse is not one described
by A and B as it's semiaxes. For example, 1 + ϵ = 2 A^2 / (A^2 + B^2)
for the lenstool definition, and 1 + ϵ = 2 A / (A + B) for the Hjorth
and Kneib definition. One ends up with ellipses

X^2 / (4 A^4 / (A^2+B^2)^2) + Y^2 / (4 B^4 / (A^2+B^2)^2) and
X^2 / (4 A^2 / (A+B)^2) + Y^2 / (4 B^2 / (A+B)^2) respectively for each
of the definitions of ϵ.

In each case, the "effective" semimajor and semiminor axes are different
than an ellipse defined by A and B. The area changes.

I don't know enough to know if this is an issue or not.

Thanks,
Joe Coleman

Kneib Jean-Paul

unread,
Aug 23, 2015, 3:25:31 AM8/23/15
to lens...@googlegroups.com
One cannot have the magnification matrix wrong with a correct deflection field without making the critical and caustic lines showing strange behaviour. So I am certain that the implementation of the magnification matrix is correctly done as both critical and caustic lines are looking good (caustic lines in particular are cuspy as expected).

It is however possibly an issue that for large ellipticity the mass distribution is not anymore looking like an ellipse but more like a peanut shape - the relation between the deflection and magnification is still correct but the mass representation is not anymore physical. And this is certainly a strong limitation.

Note however that the above only apply to approximate elliptical transformation of the NFW or gNFW profile but not for the PIEMD (profile 81) which remains elliptical even for large ellipticities. I believe (but this need to be confirmed) that some of the A2744 models were using PIEMDs and thus should not suffer from the large ellipticities issues …

Best,

Jp
> --
> --
> _____________________________________________________________________
> LENSTOOL: http://projets.lam.fr/projects/lenstool/wiki
> To post to this group, send email to lens...@googlegroups.com
> To unsubscribe from this group, send email to lenstool-u...@googlegroups.com
> For more options, visit this group at http://groups.google.com/group/lenstool?
> ---
> You received this message because you are subscribed to the Google Groups "lenstool" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to lenstool+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

Eric Jullo

unread,
Aug 23, 2015, 4:19:04 PM8/23/15
to lens...@googlegroups.com, Kneib Jean-Paul
Hi Richard

Also, I should add that you'd rather look at files e_grad.c and e_grad2.c for the calculation of the gradient and laplacian of the potential, than at file e_pot.c, which is not finished.

For the PIEMD, all the calculations are detailled in Elliasdottir et al 2008 and Kassiola &Kovner 1993.

Also, you can have a look at Golse&Kneib 2002 for the limitations of the elliptical expansion with the NFW profile.

Cheers
Eric


Eric Jullo

unread,
Aug 23, 2015, 4:41:21 PM8/23/15
to lens...@googlegroups.com, Joseph E. Coleman
Hi Joseph

Indeed ellipticity definitions are many. In Lenstool, we usually use a2-b2/a2+b2 are proposed in Golse&Kneib 2002. This is the emass parameter or ellipticity keyword in the parameter file.

However, in Kassiola&Kovner 1993 and Eliasdottir et al 2007, the PIEMD is defined with epot=a-b/a+b, therefore during optimisation, Lenstool always converts emass into epot (see o_set_lens.c and update_epot.c).

Best
Eric

Richard Massey

unread,
Aug 24, 2015, 12:09:38 PM8/24/15
to lens...@googlegroups.com
Thank you all for the replies,

Lest I offended you, let me be clear that I'm not questioning lenstool... except for the multiple streams of uncommented codes only some which are used :). I'm questioning Kassiola & Kovner (1993)'s original treatment of a PIEMD.

The coordinate transformation used to create an elliptical potential is KK93's equation 2.3.6, which is implemented in Lenstool v6.8 at line 65 of e_gcpx.c (inside routine ci05, which is called from e_grad_pot in e_grad.c).
r_em^2 = x^2/(1+e)^2 + y^2/(1-e)^2
This has a Jacobian different than unity (it is 1-e^2), which might not affect the deflection angles, but could rescale the magnification in a monotonic way (so it wouldn't be noticeable) that is only really incorrect for extremely elliptical objects.

Perhaps the correction for this Jacobian arises as the prefactor (1-e^2) in KK93's equations 2.3.8 and 4.1.2. I had originally thought this came from the conversion between definitions of ellipticity, as Joe said, but it might be more clever than that. KK93 don't exactly spell things out, and they don't do the maths the way I would.

So I'm still a bit unclear. Maybe there's nothing wrong. Either way, thank you for your patience!
Richard


PS: I'm also not questioning the peanut-shaped mass distributions for a PIEP; that's a different issue.
PPS: If I can figure all this out, I'll be adding some mass distributions to the code. Is there a good time to call/skype later in the week, to check that I'm editing all the right routines?

Richard Massey

unread,
Oct 26, 2015, 2:50:07 PM10/26/15
to lens...@googlegroups.com
Hello,

I'm trying to add a new profile to lenstool. The challenge is that this profile has (two) more parameters than a PIEMD, so I have to be able to read in, optimise, and write out those extra values. For now, I'm happy if only standalone mass clumps have these parameters, rather than everything in a potfile. I've got the code running successfully on a best.par file that does no optimisation. However, I can't get the optimisation to see these new parameters.

I'm guessing I need to edit o_global.c, but I can't figure out how (some comments would have helped :) ). Specific questions:
* Some variable names start with "g_": does this mean "global", "grid", or something else (I'm specifically thinking of g_pot)? If it's global, it might matter to me; if it's grid, it doesn't.
* One of my new parameters is an angle that should wrap at 2pi - how do I tell the optimiser not to get stuck in a local minimum at 2pi when a very small positive value would be a better fit?




What I've changed is:

include/dimension.h - Increase NPAMAX from 35 to 37
include/structure.h - Expand pot structure to include two new variables, and add variable ID numbers (35 & 36)
include/fonction.h - Add skewlens function to declarataion list

src/r_potentiel.c - Read in from .par file baseline values of new params
src/r_limit.c - Read in from .par file limits within which to optimise new parameters
src/e_pot.c - Call new subroutine to define potential
src/e_grad.c - Call new subroutine in e_grad_pot() to define potential
src/e_grad2.c - Call new subroutine in e_grad2_pot_ptr() to define potential
src/rotation.c - The new subroutine is added in this file. It could have been anywhere, but hey.
src/o_global.c - ?????
src/o_print_res.c - Write the two new parameters to best.par in writePotentiel(). I haven't changed writeLimit(), which might be a mistake.

Is there anything I've missed, such as routines that pass parameters between sections of the code?

Many thanks,
Richard

Eric Jullo

unread,
Oct 28, 2015, 4:12:13 PM10/28/15
to lens...@googlegroups.com
Hi Richard,

From your email and the list of files that you modified, it looks like you’ve already added the new potential. However, it’s not clear to me whether you really need new parameters in Lenstool, or you can just occupy already existing ones?

Adding a new parameter is a bit difficult, and it’s better if you can avoid this solution. You need to add your parameter in structure.h and make sure to update o_rescale.c and readBayesModel.c, which are called by the MCMC sampler through the bayesapp.c:UserBuild() function. Functions in these files loop over the existing parameters to check whether they are optimised or not.

Cheers
Eric

Richard Massey

unread,
Oct 29, 2015, 1:39:36 PM10/29/15
to lens...@googlegroups.com
Thank you Eric,

I think I've got Lenstool working with extra parameters. With your pointers, it isn't actually too hard - it's just knowing where to look in the first place. For future reference on the mailing list, therefore, I've appended my change log to the bottom of this email.

A couple of curiosities along the way:
1) the ID numbers of new parameters must be between CX and PMASS. You can't add them at the end. The "block" variable, which lists all the parameters to be optimised, includes a load of cosmology parameters, etc. that all get ignored. They should probably be removed, because they are really optimised from "cblock".
2) o_set_lmin.c and o_set_lmax.c are redundant. They do the same as o_set_lens.c, but they act on a variable with a different name.

Best,
Richard

CHANGES TO LENSTOOL TO ADD NEW PARAMETERS TO A POTENTIAL

include/dimension.h - Increase NPAMAX from 35 to 37
include/structure.h - Expand pot structure to include new variables... these ***must be between CX and PMASS***
include/fonction.h - Add new functions that define the new potential to declarataion list
src/r_potentiel.c - Read in from .par file baseline values of skew and skewtheta
src/r_limit.c - Read in from .par file limits within which to optimise skew and skewtheta
src/e_pot.c - Add code to create the new potential
src/e_grad.c - Add code to create the new potential
src/e_grad2.c - Add code to create the new potential
src/o_print_res.c - Write the two new parameters to best.par in writePotentiel() and bestopt.par in writeLimit().
src/o_rescale.c - List new variables in all three LT<->Bayesys lookup tables
src/o_set_lens.c - List variables in all three LT<->Bayesys lookup tables
src/o_set_lmin.c - List variables in both LT<->Bayesys lookup tables (these routines seem superfluous & o_set_lens could have been reused)
src/o_set_lmax.c - List variables in both LT<->Bayesys lookup tables (these routines seem superfluous & o_set_lens could have been reused)
src/bayesapp.c - Define human readable names for the two new variables in getParName(), convert degrees to radians for angle param in lt2bayes()
src/o_print.c - Change output to map.res... but I don't know what that is, or whether it is needed???
src/main.c - Print a quick notice that this is a forked version of LT

Richard Massey

unread,
Oct 30, 2015, 10:05:55 AM10/30/15
to lens...@googlegroups.com
Eric,

I spoke too soon. Bayesys is currently only exploring regions of parameter space (for my new parameters) where their value is positive. This is presumably determined by Common->MassInf.

You previously encountered this for positions (which can be positive or negative), and solved it by adding the minimum value as an offset so that Bayesys only ever sees positive values. But that’s commented out! It seems you had a more elegant solution. I can’t figure out what you did do - please could you give me a hint?


In a related note, it also looks from the source code that it’s possible to search for x and y with a circular top hat prior, rather than a square. Does this actually work? (It’s not mentioned in the manual, and optimisation option 4 is actually used in the manual for a different case). If this can be done, it’s actually what I’d prefer for a pair of my new parameters.

Thanks,
Richard

David Harvey

unread,
Oct 30, 2015, 10:25:18 AM10/30/15
to lens...@googlegroups.com
Hi Richard,

As far as I remember, the parameters are always passed to bayesys between 0 and 1 and then o_rescale rescales the numbers back to the correct value.
So its in o_rescale.c you want to look. i think.
David

Eric Jullo

unread,
Oct 30, 2015, 12:58:21 PM10/30/15
to lens...@googlegroups.com
Hi

Indeed, all the rescaling is done in o_rescale.c.

I'm not sure I fully understood your way of inserting the new parameter between CX and PMASS.

Why don't you upload your changes on the svn. Do you have the write access?

Eric




LENSTOOL: http://projets.lam.fr/projects/lenstool/wiki
To post to this group, send email to lens...@googlegroups.com
To unsubscribe from this group, send email to lenstool-u...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/lenstool?
---
You received this message because you are subscribed to the Google Groups "lenstool" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lenstool+u...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Richard Massey

unread,
Nov 3, 2015, 10:01:06 AM11/3/15
to lens...@googlegroups.com
Thank you Dave, Eric.
My problem was indeed in o_rescale.c. Now I've properly understood how that routine works, I see I blindly/cut/pasted something stupid. Everything is now running without hitch.

Happy to upload changes in principle if you give me write access. But I adapted David's copy rather than a downloaded version, and I don't want to break anything. I'll see you later this week, so let's figure out the best thing to do then.

Cheers,
Richard
> LENSTOOL: http://projets.lam.fr/projects/lenstool/wiki
> To post to this group, send email to lens...@googlegroups.com
> To unsubscribe from this group, send email to lenstool-u...@googlegroups.com
> For more options, visit this group at http://groups.google.com/group/lenstool?
> ---
> You received this message because you are subscribed to the Google Groups "lenstool" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to lenstool+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
> --
> --
>
> LENSTOOL: http://projets.lam.fr/projects/lenstool/wiki
> To post to this group, send email to lens...@googlegroups.com
> To unsubscribe from this group, send email to lenstool-u...@googlegroups.com
> For more options, visit this group at http://groups.google.com/group/lenstool?
> ---
> You received this message because you are subscribed to the Google Groups "lenstool" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to lenstool+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
> --
> --
> _____________________________________________________________________

Richard Massey

unread,
Nov 3, 2015, 10:04:40 AM11/3/15
to lens...@googlegroups.com
PS: new parameters for a lens profile have to be listed in structure.h between CX and PMASS (rather than at the end of the parameter list) because there are lots of loops like this (in o_rescale and throughout):
for ( ipx = CX; ipx <= PMASS; ipx++ )

On 30 Oct 2015, at 16:58, Eric Jullo wrote:

> LENSTOOL: http://projets.lam.fr/projects/lenstool/wiki
> To post to this group, send email to lens...@googlegroups.com
> To unsubscribe from this group, send email to lenstool-u...@googlegroups.com
> For more options, visit this group at http://groups.google.com/group/lenstool?
> ---
> You received this message because you are subscribed to the Google Groups "lenstool" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to lenstool+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
> --
> --
>
> LENSTOOL: http://projets.lam.fr/projects/lenstool/wiki
> To post to this group, send email to lens...@googlegroups.com
> To unsubscribe from this group, send email to lenstool-u...@googlegroups.com
> For more options, visit this group at http://groups.google.com/group/lenstool?
> ---
> You received this message because you are subscribed to the Google Groups "lenstool" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to lenstool+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
> --
> --
> _____________________________________________________________________

Eric Jullo

unread,
Nov 3, 2015, 11:21:23 AM11/3/15
to lens...@googlegroups.com
This is good news.

Can you register on the readmine of LAM

http://projects.lam.fr

Then, I will add you to the list of developers, and will get write access on the svn.

Cheers
Eric
Reply all
Reply to author
Forward
0 new messages