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}:
  -21.3988    -31.2764    -31.4762   …    13.107      15.9506     0.188393
   21.1677      9.74529     1.80161       37.2107     45.4682    37.0359
   42.5893     29.3561     14.7012        31.6938     42.9514    48.5313
   32.5094     22.7532     10.8725         8.53284    16.7403    29.9706
   -4.25489    -8.91726   -14.2057       -20.5206    -19.6585    -7.93431
  -40.9182    -43.7168    -45.2873   …   -43.0891    -44.3664   -40.1209
  -60.3192    -64.3219    -67.6532       -59.9282    -57.6744   -56.3624
  -70.7922    -75.83      -81.909        -76.5602    -72.2982   -69.9258
  -85.3       -88.599     -93.4104       -96.6156    -97.1364   -91.009
  -98.0427    -99.1436    -99.0743      -113.113    -116.868   -107.507
 -105.799    -102.063     -97.7122   …  -125.024    -126.663   -116.647
 -102.832     -93.7473    -90.4179      -130.636    -126.843   -116.368
  -88.9204    -79.1908    -83.3106      -130.583    -121.213   -106.162
    ⋮                                ⋱                            ⋮
 -125.982    -111.371     -86.6567      -154.923    -144.697   -135.956
 -156.124    -148.6      -131.106    …  -182.796    -172.557   -163.505
 -164.282    -158.233    -149.572       -197.385    -193.996   -179.491
 -158.876    -148.002    -141.678       -186.318    -191.641   -178.598
 -154.2      -137.686    -126.841       -158.498    -171.74    -168.419
 -149.632    -134.013    -124.461       -130.882    -152.612   -158.42
 -141.22     -135.52     -133.462    …  -111.001    -135.92    -144.412
 -129.132    -130.542    -135.033       -103.26     -126.577   -131.655
 -115.05     -121.505    -131.172       -101.428    -121.311   -121.25
 -103.251    -117.214    -129.534        -94.4354   -109.568   -106.441
  -91.7414   -108.93     -117.518        -72.7479    -80.5761   -83.072
  -65.5341    -79.3468    -79.9654   …   -31.9194    -34.1041   -46.6624

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 FlatMap{256×256 map, 3′ pixels, fourier∂, Array{Float32,2}}:
  -21.39875
   21.167694
   42.589256
   32.50937
   -4.2548904
  -40.918236
  -60.319244
  -70.7922
  -85.30002
  -98.042725
 -105.79917
 -102.83204
  -88.92036
    ⋮
 -135.95642
 -163.50516
 -179.49149
 -178.59776
 -168.41861
 -158.4197
 -144.41196
 -131.65543
 -121.24965
 -106.440735
  -83.07201
  -46.662354

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:  61.10 MiB
  allocs estimate:  1479
  --------------
  minimum time:     22.827 ms (0.00% GC)
  median time:      26.768 ms (13.88% GC)
  mean time:        25.731 ms (10.20% GC)
  maximum time:     27.315 ms (14.17% GC)
  --------------
  samples:          195
  evals/sample:     1

Once cached, it's very fast and memory non-intensive to repeatedly apply the operator:

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

BenchmarkTools.Trial: 
  memory estimate:  16.13 MiB
  allocs estimate:  329
  --------------
  minimum time:     30.037 ms (0.00% GC)
  median time:      30.458 ms (0.00% GC)
  mean time:        31.189 ms (2.50% GC)
  maximum time:     34.298 ms (11.08% GC)
  --------------
  samples:          86
  evals/sample:     1