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 Array{Float32,2}:
-3.31253 14.0475 27.3787 … -10.557 -6.519 -10.19
2.20245 16.5875 29.0216 -2.42106 1.42281 -1.18434
5.9539 16.229 28.4277 6.82347 11.9188 8.74619
1.69503 9.10588 17.3534 10.1779 12.8088 5.35868
-9.1377 -3.95205 -4.11731 2.73597 -0.7896 -9.82864
-28.9114 -18.8228 -18.7714 … -16.3268 -29.1551 -36.0737
-50.257 -26.8411 -15.1703 -38.8677 -59.2921 -65.0812
-51.5699 -21.0664 -3.42119 -49.0822 -68.9358 -71.6908
-33.0832 -11.5336 1.08528 -32.1746 -46.8265 -48.7597
-13.9335 -6.189 -3.18058 2.42665 -8.96469 -16.5865
-3.7652 -7.91744 -15.7189 … 33.3868 20.7694 5.03484
0.415972 -13.427 -27.8734 56.0851 42.5448 19.2422
3.94791 -18.5734 -39.8846 71.606 59.2228 32.6523
⋮ ⋱ ⋮
137.327 125.434 103.363 152.084 153.853 144.135
119.191 116.355 93.5048 … 95.2818 103.603 108.325
100.466 104.721 85.4876 47.8358 63.3314 80.1434
80.9991 87.9602 78.0988 16.6928 38.6817 61.1498
66.9741 71.3036 66.4138 -7.2442 21.0998 50.0837
57.9981 52.2749 49.5042 -23.9316 9.67734 45.752
50.7451 34.3516 31.2096 … -29.6499 10.064 47.2559
46.5442 29.0273 26.6146 -23.4841 16.9568 48.5907
38.9765 34.2579 34.1723 -11.9524 14.7834 35.2082
20.9971 28.952 35.4361 -2.19332 4.01476 12.1652
1.18856 16.298 28.5187 -2.91687 -3.03068 -4.95314
-5.94884 10.9649 25.0863 … -9.9566 -7.01671 -10.9484
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}}:
-3.3125343
2.2024498
5.9539013
1.695034
-9.137705
-28.911417
-50.257034
-51.56987
-33.083244
-13.933515
-3.7651958
0.41597176
3.9479141
⋮
144.13489
108.32519
80.14345
61.149803
50.083736
45.75199
47.255924
48.5907
35.208202
12.165173
-4.9531403
-10.948444
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): 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 cache(LenseFlow(ϕ),f)
BenchmarkTools.Trial:
memory estimate: 92.23 MiB
allocs estimate: 1875
--------------
minimum time: 45.408 ms (6.29% GC)
median time: 48.340 ms (5.91% GC)
mean time: 51.973 ms (12.15% GC)
maximum time: 215.740 ms (77.86% GC)
--------------
samples: 97
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.17 MiB
allocs estimate: 615
--------------
minimum time: 28.763 ms (0.00% GC)
median time: 30.916 ms (0.00% GC)
mean time: 31.270 ms (1.27% GC)
maximum time: 34.722 ms (7.91% GC)
--------------
samples: 62
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.