Lensing a flat-sky map

using CMBLensing, PythonPlot

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

(;ds, 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 view(::Matrix{Float32}, :, :) with eltype Float32:
   4.48888   17.2694    34.8424   …  -23.5769     -20.686     -7.06294
  22.7025    30.3171    41.4717      -36.9175     -20.2562     7.50488
  57.3455    57.8147    56.0213      -28.7438       1.41138   40.1413
  82.2528    83.1207    72.5825       -9.23557     25.246     64.437
  82.1103    87.4416    79.9426       12.6929      37.9921    65.3094
  61.2477    71.8372    75.5525   …   30.5254      40.0434    50.295
  37.2227    49.5737    60.2065       39.9043      35.3133    32.5213
  24.684     29.3388    37.6891       54.6839      36.5755    26.7114
  14.2699     7.31913    9.04513      78.8774      49.053     27.6696
  -9.82414  -21.5851   -22.9233       84.6395      45.6952    11.6015
 -44.0145   -53.3747   -48.3444   …   64.2089      20.403    -19.5261
 -72.3802   -76.9141   -59.4429       39.9176      -7.57407  -47.9719
 -86.0061   -82.7494   -55.3906       26.8705     -28.0548   -67.3992
   ⋮                              ⋱                            ⋮
  63.4774    63.6027    52.8885      116.72        80.7807    63.8097
  74.8238    67.5884    46.1743   …  123.343       92.7058    77.968
  83.0827    70.0984    43.2644      118.822       99.3417    89.9547
  77.3378    65.1908    39.0946       94.0322      84.0634    81.6987
  63.8614    57.1189    30.716        55.6022      48.4997    57.4085
  48.5046    47.7565    24.1952       22.4734      14.3792    30.1207
  35.9928    39.5574    21.4412   …    5.68852     -3.08851   13.5288
  26.9482    27.4882    17.5991       -0.211002    -7.16303    8.59088
  13.1345    12.1954    16.1932       -0.0904846   -5.79985    6.16924
   4.98869    5.4613    21.0252        3.54408      1.68818    7.07243
   3.76636    8.03388   26.4022        7.4056       5.90948    7.29354
   1.9529    12.5011    31.2693   …   -1.13697     -3.31901   -1.12265

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 LambertMap{SubArray{Float32, 2, Matrix{Float32}, Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}}, true}}:
   4.488884
  22.70248
  57.345524
  82.25276
  82.110344
  61.24771
  37.222744
  24.683971
  14.269928
  -9.824142
 -44.014496
 -72.38022
 -86.00611
   ⋮
  63.80973
  77.96801
  89.95467
  81.6987
  57.408516
  30.120705
  13.528763
   8.590881
   6.169239
   7.0724335
   7.293541
  -1.1226501

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 precompute!!(LenseFlow(ϕ),f)
BenchmarkTools.Trial: 271 samples with 1 evaluation.
 Range (minmax):  14.959 ms22.800 ms   GC (min … max):  0.00% … 32.62%
 Time  (median):     19.052 ms               GC (median):    19.60%
 Time  (mean ± σ):   18.473 ms ±  1.480 ms   GC (mean ± σ):  17.06% ±  7.51%

                                                    ▁▃▄█▄     
  ▃▃▆▅▆▃▂▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂▄▅█████▇▆▄ ▃
  15 ms           Histogram: frequency by time        19.5 ms <

 Memory estimate: 62.99 MiB, allocs estimate: 804.

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

@benchmark Lϕ * f setup=(Lϕ=precompute!!(LenseFlow(ϕ),f))
BenchmarkTools.Trial: 99 samples with 1 evaluation.
 Range (minmax):  30.385 ms35.487 ms   GC (min … max): 0.00% … 10.79%
 Time  (median):     30.752 ms               GC (median):    0.00%
 Time  (mean ± σ):   31.776 ms ±  1.772 ms   GC (mean ± σ):  3.19% ±  4.92%

    █                                                          
  ▃▇██▄▃▃▄▅▃▃▃▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▃▃▅▆▃▃▃▃▃▃ ▁
  30.4 ms         Histogram: frequency by time        35.2 ms <

 Memory estimate: 16.13 MiB, allocs estimate: 441.

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.