Multiplane lensing only agrees with lenstronomy to 1e-2

[ConnorStoneAstro]

When running multiplane lensing on two SIEs aligned in two separate planes there is a difference in the caustic results and lenstronomy results of order 2e-3. This can be reproduced by running the multiplane_test.py code in the test folder.

When running a single SIE plane the results agree to high precision:

When adding a second plane there must be something magnifying the errors because suddenly the difference goes up to 2e-3:

There is also clear structure in the errors. Note that these images are created just by taking the log of the abs difference between the two alpha maps (in the x-axis).

For this issue we need to determine the source of these differences. Possible explanations could be: problem with our cosmology calculator, some unit conversion is different between packages, lenstronomy sub-pixel integration, numerical precision (float32 for pytorch vs float64 for lenstronomy), numerical stability when evaluating SIE is different between packages, something else?

[adam-coogan]

Some notes based on tests on the param-refactor branch:

Difference is around 5e-3

Using the default lenstronomy cosmology gives worse agreement. Since our cosmology doesn't have radiation, should instead use astropy.cosmology.FlatLambdaCDM with Tcmb0=0. With that, the comoving distance agrees fairly well with our FlatLambdaCDM cosmology.

Lenstronomy doesn't appear to be doing any fancy sub-pixel integrations -- see ray_shooting_partial here.

Using float32 vs float64 makes no difference.