Lensing a flat-sky map
using CMBLensing, PythonPlotFirst 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̃);
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.12265If 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.1226501For 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:
\[ 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): likePowerLens, 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 precompute!!(LenseFlow(ϕ),f)BenchmarkTools.Trial: 271 samples with 1 evaluation.
Range (min … max): 14.959 ms … 22.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 (min … max): 30.385 ms … 35.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.