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)

┌ Warning: `vendor()` is deprecated, use `BLAS.get_config()` and inspect the output instead
│   caller = npyinitialize() at numpy.jl:67
└ @ PyCall /home/cosmo/.julia/packages/PyCall/L0fLP/src/numpy.jl:67

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 view(::Array{Float32, 3}, :, :, 1) with eltype Float32:
   41.8353     54.0098     50.4641   …    38.43       34.3761     34.1809
   81.7465     92.5089     83.8645        65.6474     65.5049     71.2286
  100.451     106.702      94.0928        78.2358     80.342      88.5211
  101.145     104.178      89.7162        77.7212     83.9819     91.8985
   96.2317     94.399      82.577         70.7831     84.4593     92.9933
   88.3665     83.2739     80.441    …    61.5532     79.7542     90.1173
   83.936      81.4221     92.6486        51.0435     72.1622     85.9434
   90.7489     93.4388    116.749         42.5304     67.644      87.3587
   95.9975    103.625     129.782         38.9614     65.8996     88.503
   97.6787    108.879     135.209         36.8716     64.6788     87.1346
  105.046     122.181     150.381    …    27.136      60.4673     88.7536
  109.676     133.492     158.806          6.03852    41.8113     80.6451
  101.034     129.116     149.835        -11.7109     18.4085     61.0714
    ⋮                                ⋱                             ⋮
  -41.5532    -46.4189    -47.6034       -67.8714    -52.2147    -40.2831
  -67.8366    -71.1633    -74.7697   …   -84.8917    -78.9973    -69.2222
  -95.8941    -93.8205    -98.6414      -100.355    -107.633    -102.299
 -109.13     -101.932    -102.562        -96.4713   -111.967    -114.778
  -97.0635    -88.7144    -85.4534       -63.6396    -84.8965    -98.0221
  -69.738     -67.6184    -64.774        -20.5698    -42.7033    -64.4003
  -39.8173    -44.0026    -45.5737   …     3.72352    -8.426     -28.8865
  -23.6978    -23.9288    -25.5096         5.14941     1.38932   -16.2738
  -27.651     -18.5649    -15.5099        -2.68699   -11.9077    -26.6581
  -37.5459    -27.0872    -23.7102        -8.6441    -25.2734    -38.0606
  -32.0033    -25.6989    -26.448         -2.52193   -18.7722    -29.9463
   -2.98109     5.98887     3.38267  …    13.674       3.09321    -4.1559

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}}:
   41.83526
   81.74647
  100.45109
  101.14453
   96.23173
   88.36646
   83.93596
   90.74887
   95.99751
   97.67871
  105.04609
  109.67584
  101.03371
    ⋮
  -40.283092
  -69.222244
 -102.299255
 -114.77829
  -98.02206
  -64.40033
  -28.886452
  -16.273758
  -26.658096
  -38.060562
  -29.94632
   -4.155899

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: 104 samples with 1 evaluation.
 Range (min … max):  44.724 ms … 59.184 ms  ┊ GC (min … max): 10.13% … 8.38%
 Time  (median):     48.006 ms              ┊ GC (median):    10.22%
 Time  (mean ± σ):   48.489 ms ±  2.476 ms  ┊ GC (mean ± σ):  11.44% ± 2.99%

       ▃    ▄   ▁█ ▃▄   ▄                                      
  ▄▄▆▄▆█▇▄▆▆█▇▇▇██▆██▇▆▆█▄▄▆▇▆▆▆▁▄▁▇▆▁▁▄▁▁▄▄▄▄▁▆▄▄▁▄▁▁▁▁▁▁▁▁▄ ▄
  44.7 ms         Histogram: frequency by time        55.7 ms <

 Memory estimate: 92.40 MiB, allocs estimate: 1623.

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

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

BenchmarkTools.Trial: 53 samples with 1 evaluation.
 Range (min … max):  42.106 ms … 57.010 ms  ┊ GC (min … max): 0.00% … 8.55%
 Time  (median):     44.333 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   45.750 ms ±  3.345 ms  ┊ GC (mean ± σ):  3.56% ± 4.85%

  ▁▁▄██▁▁  █▁  ▁▁      ▁        ▁   ▁▁                         
  ███████▆▆██▁▁██▆▁▆▁▁▁█▁▁▁▁▆▆▆▆█▆▁▆██▆▆▆▁▁▁▁▁▁▁▁▁▁▆▁▁▁▁▁▁▁▁▆ ▁
  42.1 ms         Histogram: frequency by time        54.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.