Lensing a flat-sky map
using CMBLensing, PyPlotFirst 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 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.5465If 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.5465For 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 cache(LenseFlow(ϕ),f)BenchmarkTools.Trial: 242 samples with 1 evaluation.
Range (min … max): 17.492 ms … 22.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 (min … max): 31.513 ms … 37.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.