Lensing a flat-sky map

using CMBLensing, PyPlot

First we load a simulated unlensed field, $f$, and lensing potential, $\phi$,

@unpack f,ϕ = load_sim(
    θpix  = 2,       # size of the pixels in arcmin
    Nside = 256,     # number of pixels per side in the map
    T     = Float32, # Float32 or Float64 (former is ~twice as fast)
    pol   = :I       # :I for Intensity, :P for polarization, or :IP for both
);

We can lense the map with LenseFlow,

f̃ = LenseFlow(ϕ) * f;

And flip between lensed and unlensed maps,

animate([f,f̃], fps=1)

The difference between lensed and unlensed,

plot(f-f̃);

png

Loading your own data

CMBLensing flat-sky Field objects like f or ϕ are just thin wrappers around arrays. You can get the underlying data arrays for $I(\mathbf{x})$, $Q(\mathbf{x})$, and $U(\mathbf{x})$ with f[:Ix], f[:Qx], and f[:Ux] respectively, or the Fourier coefficients, $I(\mathbf{l})$, $Q(\mathbf{l})$, and $U(\mathbf{l})$ with f[:Il], f[:Ql], and f[:Ul],

mapdata = f[:Ix]

256×256 Array{Float32,2}:
  -3.31253    14.0475    27.3787   …  -10.557     -6.519    -10.19
   2.20245    16.5875    29.0216       -2.42106    1.42281   -1.18434
   5.9539     16.229     28.4277        6.82347   11.9188     8.74619
   1.69503     9.10588   17.3534       10.1779    12.8088     5.35868
  -9.1377     -3.95205   -4.11731       2.73597   -0.7896    -9.82864
 -28.9114    -18.8228   -18.7714   …  -16.3268   -29.1551   -36.0737
 -50.257     -26.8411   -15.1703      -38.8677   -59.2921   -65.0812
 -51.5699    -21.0664    -3.42119     -49.0822   -68.9358   -71.6908
 -33.0832    -11.5336     1.08528     -32.1746   -46.8265   -48.7597
 -13.9335     -6.189     -3.18058       2.42665   -8.96469  -16.5865
  -3.7652     -7.91744  -15.7189   …   33.3868    20.7694     5.03484
   0.415972  -13.427    -27.8734       56.0851    42.5448    19.2422
   3.94791   -18.5734   -39.8846       71.606     59.2228    32.6523
   ⋮                               ⋱                          ⋮
 137.327     125.434    103.363       152.084    153.853    144.135
 119.191     116.355     93.5048   …   95.2818   103.603    108.325
 100.466     104.721     85.4876       47.8358    63.3314    80.1434
  80.9991     87.9602    78.0988       16.6928    38.6817    61.1498
  66.9741     71.3036    66.4138       -7.2442    21.0998    50.0837
  57.9981     52.2749    49.5042      -23.9316     9.67734   45.752
  50.7451     34.3516    31.2096   …  -29.6499    10.064     47.2559
  46.5442     29.0273    26.6146      -23.4841    16.9568    48.5907
  38.9765     34.2579    34.1723      -11.9524    14.7834    35.2082
  20.9971     28.952     35.4361       -2.19332    4.01476   12.1652
   1.18856    16.298     28.5187       -2.91687   -3.03068   -4.95314
  -5.94884    10.9649    25.0863   …   -9.9566    -7.01671  -10.9484

If you have your own map data in an array you'd like to load into a CMBLensing Field object, you can construct it as follows:

FlatMap(mapdata, θpix=3)

65536-element 256×256-pixel 3.0′-resolution FlatMap{Array{Float32,2},ProjLambert{Float32}}:
  -3.3125343
   2.2024498
   5.9539013
   1.695034
  -9.137705
 -28.911417
 -50.257034
 -51.56987
 -33.083244
 -13.933515
  -3.7651958
   0.41597176
   3.9479141
   ⋮
 144.13489
 108.32519
  80.14345
  61.149803
  50.083736
  45.75199
  47.255924
  48.5907
  35.208202
  12.165173
  -4.9531403
 -10.948444

For more info on Field objects, see Field Basics.

Inverse lensing

You can inverse lense a map with the \ operator (which does A \ b ≡ inv(A) * b):

LenseFlow(ϕ) \ f;

Note that this is true inverse lensing, rather than lensing by the negative deflection (which is often called "anti-lensing"). This means that lensing then inverse lensing a map should get us back the original map. Lets check that this is the case:

Ns = [7 10 20]
plot([f - (LenseFlow(ϕ,N) \ (LenseFlow(ϕ,N) * f)) for N in Ns],
    title=["ODE steps = $N" for N in Ns]);

png

A cool feature of LenseFlow is that inverse lensing is trivially done by running the LenseFlow ODE in reverse. Note that as we crank up the number of ODE steps above, we recover the original map to higher and higher precision.

Other lensing algorithms

We can also lense via:

PowerLens: the standard Taylor series expansion to any order:

\[ f(x+\nabla x) \approx f(x) + (\nabla f)(\nabla \phi) + \frac{1}{2} (\nabla \nabla f) (\nabla \phi)^2 + ... \]

TayLens (Næss&Louis 2013): like PowerLens, but first a nearest-pixel permute step, then a Taylor expansion around the now-smaller residual displacement

plot([(PowerLens(ϕ,2)*f - f̃) (Taylens(ϕ,2)*f - f̃)], 
    title=["PowerLens - LenseFlow" "TayLens - LenseFlow"]);

png

Benchmarking

LenseFlow is highly optimized code since it appears on the inner-most loop of our analysis algorithms. To benchmark LenseFlow, note that there is first a precomputation step, which caches some data in preparation for applying it to a field of a given type. This was done automatically when evaluating LenseFlow(ϕ) * f but we can benchmark it separately since in many cases this only needs to be done once for a given $\phi$, e.g. when Wiener filtering at fixed $\phi$,

using BenchmarkTools

@benchmark cache(LenseFlow(ϕ),f)

BenchmarkTools.Trial: 
  memory estimate:  92.23 MiB
  allocs estimate:  1875
  --------------
  minimum time:     45.408 ms (6.29% GC)
  median time:      48.340 ms (5.91% GC)
  mean time:        51.973 ms (12.15% GC)
  maximum time:     215.740 ms (77.86% GC)
  --------------
  samples:          97
  evals/sample:     1

Once cached, it's faster and less memory intensive to repeatedly apply the operator:

@benchmark Lϕ * f setup=(Lϕ=cache(LenseFlow(ϕ),f))

BenchmarkTools.Trial: 
  memory estimate:  16.17 MiB
  allocs estimate:  615
  --------------
  minimum time:     28.763 ms (0.00% GC)
  median time:      30.916 ms (0.00% GC)
  mean time:        31.270 ms (1.27% GC)
  maximum time:     34.722 ms (7.91% GC)
  --------------
  samples:          62
  evals/sample:     1

Note that this documentation is generated on limited-performance cloud servers. Actual benchmarks are likely much faster locally or on a cluster, and yet (much) faster on GPU.