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 view(::Array{Float32, 3}, :, :, 1) with eltype Float32:
  56.8072    61.1838    75.8541   …  108.278     92.9383     69.1366
  56.8215    61.6501    80.3612      105.818     96.134      72.3947
  70.1948    83.862    102.748        99.1963    93.6627     75.8874
  93.5305   111.751    125.28         85.5538    89.4707     86.9946
 112.496    128.569    135.987        75.4355    85.6198     97.0798
 116.538    133.269    138.971    …   71.2871    80.273      94.81
 116.572    138.438    141.482        70.5923    78.4015     91.4552
 125.75     148.993    144.755        62.4825    75.7289     94.9209
 127.903    147.538    137.664        46.8252    66.6549     94.3554
 107.914    123.38     117.424        24.5269    49.5657     79.2146
  79.7204    91.204     89.7127   …   -2.02958   23.5113     54.7515
  52.0153    60.2521    59.3639      -27.2743    -0.723793   29.8547
  23.7429    34.0696    36.5107      -39.6659   -17.0888      5.19127
   ⋮                              ⋱                           ⋮
 -12.5636    -8.73611  -20.8456      -12.9603   -27.1412    -26.3102
 -15.9506   -17.0823   -29.9456   …   11.1526   -10.3675    -20.4711
 -16.5631   -20.3669   -31.3976       28.812      5.27584   -11.4729
  -8.66104  -10.8845   -20.6668       36.7089    14.0519     -2.74534
   9.31844    9.54618   -2.95713      38.6585    20.0181      9.09072
  33.9152    30.9992    14.1011       41.4053    31.0849     29.7679
  55.3008    46.4585    30.0591   …   55.5677    51.2278     54.1446
  63.1386    50.8445    37.5707       77.2633    70.8254     68.6741
  54.1671    43.8095    36.9622       92.5483    74.4283     63.445
  44.3789    40.0526    41.9321      105.955     76.4456     55.1704
  48.8498    54.1686    62.0053      117.894     83.9559     57.2095
  56.4118    65.2036    76.6114   …  116.459     89.0266     63.5465

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, Array{Float32, 3}, Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Int64}, true}}:
  56.807228
  56.821537
  70.19481
  93.53048
 112.49566
 116.5383
 116.57192
 125.75002
 127.90324
 107.91376
  79.72035
  52.01527
  23.742939
   ⋮
 -26.31018
 -20.471123
 -11.472891
  -2.7453413
   9.090725
  29.767918
  54.14463
  68.674126
  63.445015
  55.17041
  57.209473
  63.5465

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: 242 samples with 1 evaluation.
 Range (minmax):  17.492 ms22.310 ms   GC (min … max):  0.00% … 17.65%
 Time  (median):     21.815 ms               GC (median):    17.86%
 Time  (mean ± σ):   20.642 ms ±  1.924 ms   GC (mean ± σ):  13.25% ±  8.23%

                                          ▅▁     
  ▅▇█▇▅▃▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▃▅███▇▅▂ ▃
  17.5 ms         Histogram: frequency by time        22.3 ms <

 Memory estimate: 62.98 MiB, allocs estimate: 804.

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

@benchmark Lϕ * f setup=(Lϕ=cache(LenseFlow(ϕ),f))
BenchmarkTools.Trial: 93 samples with 1 evaluation.
 Range (minmax):  31.513 ms37.617 ms   GC (min … max):  0.00% … 10.93%
 Time  (median):     36.040 ms               GC (median):    11.10%
 Time  (mean ± σ):   35.096 ms ±  1.879 ms   GC (mean ± σ):   8.50% ±  4.89%

    ▂                                 █▂▁             
  ▃▇█▆▁▆▁▄▃▁▁▃▃▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▄▇█████▇▄▃▃▁▁▁▃▁▃▃ ▁
  31.5 ms         Histogram: frequency by time        37.5 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.