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̃);

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 Matrix{Float32}:
  -2.49814   -30.0118   -43.6573    …   19.8136    17.9947     16.5335
   1.04933   -26.3951   -45.8329         5.1517    10.1516     14.0571
   0.382877  -24.5437   -49.9398        -3.77188    4.66585     8.95565
  -5.46502   -31.7157   -61.9552        -2.1017     4.47664     4.33123
 -17.8687    -47.3206   -77.371         -7.64374   -1.72768    -3.54133
 -25.1634    -54.1203   -79.446     …  -13.5888    -5.57747    -7.71585
 -23.3903    -44.8229   -61.0596        -3.74521   -0.863227   -7.25933
 -15.4169    -26.2807   -38.6727        13.4585     5.63699    -6.89403
  -5.72618   -11.5179   -23.0099        22.5172     8.6737     -2.38327
  11.9088      1.43139  -12.2472        32.4663    21.7151     16.6796
  30.056      11.2163    -6.96114   …   44.5476    43.0212     42.5209
  33.4909      8.28001  -12.0563        40.4607    47.0118     50.167
  22.9502     -1.19586  -22.1942        17.6934    26.4536     34.5812
   ⋮                                ⋱                           ⋮
 132.433     109.636     76.9897        84.3316   106.871     130.217
 120.938     106.796     78.5164    …   51.3455    76.5028    107.008
 111.424     104.259     81.0615        33.7821    60.107      93.5488
 112.579     108.457     86.3929        31.421     56.0891     90.6019
 109.626     108.075     83.8451        38.563     51.4494     81.7154
  85.1565     86.4289    61.9635        47.6084    42.8823     59.2365
  51.5797     49.5673    29.4711    …   56.6562    39.9719     40.6048
  21.1516     13.7437    -0.181247      60.5717    34.2714     24.2754
  -4.76942   -14.4064   -23.9993        51.7697    19.6946      3.79754
 -19.4234    -30.1529   -37.743         34.1821     8.82407    -7.83921
 -18.531     -33.6154   -40.534         24.5933    10.5523     -2.99186
 -10.013     -32.4222   -41.5184    …   23.9008    16.8637      8.3991

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}}:
  -2.4981384
   1.0493317
   0.38287735
  -5.465025
 -17.86869
 -25.16341
 -23.390274
 -15.416857
  -5.726185
  11.9087925
  30.056025
  33.490944
  22.950226
   ⋮
 130.21707
 107.007866
  93.54883
  90.60188
  81.71538
  59.236465
  40.604824
  24.275356
   3.797536
  -7.8392086
  -2.9918566
   8.399103

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]);

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:

• 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"]);

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.35 MiB
  allocs estimate:  1653
  --------------
  minimum time:     51.047 ms (5.76% GC)
  median time:      53.815 ms (5.73% GC)
  mean time:        55.404 ms (6.78% GC)
  maximum time:     83.476 ms (3.56% GC)
  --------------
  samples:          91
  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.14 MiB
  allocs estimate:  441
  --------------
  minimum time:     35.937 ms (0.00% GC)
  median time:      38.916 ms (0.00% GC)
  mean time:        39.678 ms (2.28% GC)
  maximum time:     47.413 ms (0.00% GC)
  --------------
  samples:          51
  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.