Field Basics

using CMBLensing

Base Fields

The basic building blocks of CMBLensing.jl are CMB "fields", like temperature, Q or U polarization, or the lensing potential $\phi$. These types are all encompassed by the abstract type Field, with some concrete examples including FlatMap for a flat-sky map projection, or FlatQUMap for Q/U polarization, etc...

Flat fields are just thin wrappers around Julia arrays, e.g.

Ix = rand(2,2)
2×2 Matrix{Float64}:
 0.459912  0.966169
 0.696713  0.516243
f = FlatMap(Ix)
4-element 2×2-pixel 1.0′-resolution LambertMap{Array{Float64, 2}}:
 0.4599124117435962
 0.6967131274166374
 0.9661685855933742
 0.5162432729683115

When displayed, you can see the pixels in the 2x2 map have been splayed out into a length-4 array. This is intentional, as even though the maps themselves are two-dimensional, it is extremely useful conceptually to think of fields as vectors (which they are, in fact, as they form an abstract vector space). This tie to vector spaces is deeply rooted in CMBLensing, to the extent that Field objects are a subtype of Julia's own AbstractVector type,

f isa AbstractVector
true

The data itself, however, is still stored as the original 2x2 matrix, and can be accessed as follows,

f.Ix
2×2 Matrix{Float64}:
 0.459912  0.966169
 0.696713  0.516243

But since Fields are vectors, they can be tranposed,

f'
1×4 adjoint(::LambertMap{Array{Float64, 2}}) with eltype Float64:
 0.459912  0.696713  0.966169  0.516243

inner products can be computed,

f' * f
1.8969174610630186

and they can be added with each other as well as multiplied by scalars,

2*f+f
4-element 2×2-pixel 1.0′-resolution LambertMap{Array{Float64, 2}}:
 1.3797372352307886
 2.0901393822499124
 2.8985057567801222
 1.5487298189049346

Diagonal operators

Vector spaces have linear operators which act on the vectors. Linear operators correpsond to matrices, thus for a map with $N$ total pixels, a general linear operator would be an $N$-by-$N$ matrix, which for even modest map sizes becomes far too large to actually store. Thus, an important class of linear operators are ones which are diagonal, since these can actually be stored. CMBLensing uses Julia's builtin Diagonal to represent these. Diagonal(f) takes a vector f and puts it on the diagonal of the matrix:

Diagonal(f)
4×4 Diagonal{Float64, LambertMap{Array{Float64, 2}}}:
 0.459912   ⋅         ⋅         ⋅ 
  ⋅        0.696713   ⋅         ⋅ 
  ⋅         ⋅        0.966169   ⋅ 
  ⋅         ⋅         ⋅        0.516243

Multiplying this operator by the original map is then a matrix-vector product:

Diagonal(f) * f
4-element 2×2-pixel 1.0′-resolution LambertMap{Array{Float64, 2}}:
 0.21151942647581118
 0.4854091819146716
 0.9334817357875012
 0.2665071168850346

Note that this is also equal to the the pointwise multiplication of f with itself:

f .* f
4-element 2×2-pixel 1.0′-resolution LambertMap{Array{Float64, 2}}:
 0.21151942647581118
 0.4854091819146716
 0.9334817357875012
 0.2665071168850346

Field Tuples

You can put Fields together into tuples. For example,

a = FlatMap(rand(2,2))
b = FlatMap(rand(2,2));
FieldTuple(a,b)
8-element Field-2-Tuple{LambertMap{Array{Float64, 2}}, LambertMap{Array{Float64, 2}}}:
 0.09914786391952712
 0.18598410611302874
 0.6349936166558303
 0.6142728589097468
 0.9058994959088245
 0.5481505048566739
 0.6767661193126735
 0.6113430607196322

The components can also have names:

ft = FieldTuple(a=a, b=b)
8-element Field-(a,b)-Tuple{LambertMap{Array{Float64, 2}}, LambertMap{Array{Float64, 2}}}:
 0.09914786391952712
 0.18598410611302874
 0.6349936166558303
 0.6142728589097468
 0.9058994959088245
 0.5481505048566739
 0.6767661193126735
 0.6113430607196322

which can be accessed later:

ft.a
4-element 2×2-pixel 1.0′-resolution LambertMap{Array{Float64, 2}}:
 0.09914786391952712
 0.18598410611302874
 0.6349936166558303
 0.6142728589097468

FieldTuples have all of the same behavior of individual fields. Indeed, spin fields like QU or IQU are simply special FieldTuples:

fqu = FlatQUMap(a,b)
fqu isa FieldTuple
false

Field Vectors

in progress

Basis Conversion

All fields are tagged as to which basis they are stored in. You can convert them to other bases by calling the basis type on them:

f
4-element 2×2-pixel 1.0′-resolution LambertMap{Array{Float64, 2}}:
 0.4599124117435962
 0.6967131274166374
 0.9661685855933742
 0.5162432729683115
f′ = Fourier(f)
4-element 2×2-pixel 1.0′-resolution LambertFourier{Array{ComplexF64, 2}}:
  2.6390373977219195 + 0.0im
 0.21312459695202146 + 0.0im
 -0.3257863194014521 + 0.0im
 -0.6867260282981038 + 0.0im

Basis conversion is usually done automatically for you. E.g. here f′ is automatically converted to a FlatMap before addition:

f + f′
4-element 2×2-pixel 1.0′-resolution LambertMap{Array{Float64, 2}}:
 0.9198248234871925
 1.3934262548332748
 1.9323371711867483
 1.032486545936623

A key feature of Diagonal operators is they convert the field they are acting on to the right basis before multiplication:

Diagonal(f) * f′
4-element 2×2-pixel 1.0′-resolution LambertMap{Array{Float64, 2}}:
 0.2115194264758112
 0.4854091819146717
 0.9334817357875012
 0.2665071168850346

A FlatMap times a FlatFourier doesn't have a natural linear algebra meaning so its an error:

f * f′
MethodError: no method matching *(::LambertMap{Array{Float64, 2}}, ::LambertFourier{Array{ComplexF64, 2}})



Stacktrace:

 [1] top-level scope

   @ In[21]:1

 [2] eval

   @ ./boot.jl:360 [inlined]

 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)

   @ Base ./loading.jl:1094
┌ Error: Error showing method candidates, aborted
│   exception = (ErrorException("invoke: argument type error"), Union{Ptr{Nothing}, Base.InterpreterIP}[Ptr{Nothing} @0x00007f5a34887859, Ptr{Nothing} @0x00007f5a3490dc63, Ptr{Nothing} @0x00007f5a1fcaa456, Ptr{Nothing} @0x00007f5a1fcaa46f, Ptr{Nothing} @0x00007f5a348fc719, Ptr{Nothing} @0x00007f5a263c8ac6, Ptr{Nothing} @0x00007f5a1fcaa359, Ptr{Nothing} @0x00007f5a1fcaa3e5, Ptr{Nothing} @0x00007f5a348fc719, Ptr{Nothing} @0x00007f5a26707914, Ptr{Nothing} @0x00007f5a1fc38603, Ptr{Nothing} @0x00007f5a1fc38865, Ptr{Nothing} @0x00007f5a348fc719, Ptr{Nothing} @0x00007f5a25ff4b07, Ptr{Nothing} @0x00007f5a25ff4baa, Ptr{Nothing} @0x00007f5a348fc719, Ptr{Nothing} @0x00007f5a2660f0ea, Ptr{Nothing} @0x00007f5a1fca9b1f, Ptr{Nothing} @0x00007f5a1fca9b81, Ptr{Nothing} @0x00007f5a348fc719, Ptr{Nothing} @0x00007f5a1fca8e57, Ptr{Nothing} @0x00007f5a1fca96a7, Ptr{Nothing} @0x00007f5a1fca970c, Ptr{Nothing} @0x00007f5a348fc719, Ptr{Nothing} @0x00007f5a1fca2262, Ptr{Nothing} @0x00007f5a1fca60e8, Ptr{Nothing} @0x00007f5a1fca7bbc, Ptr{Nothing} @0x00007f5a1fca7cc7, Ptr{Nothing} @0x00007f5a1fca7d33, Ptr{Nothing} @0x00007f5a348fc719, Ptr{Nothing} @0x00007f5a3490d526, Ptr{Nothing} @0x00007f5a1fca148e, Ptr{Nothing} @0x00007f5a1fca1513, Ptr{Nothing} @0x00007f5a348fc719, Ptr{Nothing} @0x00007f5a3490f808, Ptr{Nothing} @0x00007f5a1fca1282, Ptr{Nothing} @0x00007f5a348fc719, Ptr{Nothing} @0x00007f5a3490f808, Ptr{Nothing} @0x00007f5a1fca1118, Ptr{Nothing} @0x00007f5a348fc719, Ptr{Nothing} @0x00007f5a1fc96e24, Ptr{Nothing} @0x00007f5a1fc970d1, Ptr{Nothing} @0x00007f5a1fc97135, Ptr{Nothing} @0x00007f5a348fc719, Ptr{Nothing} @0x00007f5a1fc2e299, Ptr{Nothing} @0x00007f5a348fc719, Ptr{Nothing} @0x00007f5a3490d526, Ptr{Nothing} @0x00007f5a1fbf58be, Ptr{Nothing} @0x00007f5a1fbf5aee, Ptr{Nothing} @0x00007f5a1fbf5b0c, Ptr{Nothing} @0x00007f5a348fc719, Ptr{Nothing} @0x00007f5a3491d84a, Ptr{Nothing} @0x0000000000000000])
└ @ Base errorshow.jl:329

Properties and indices

FlatMap and FlatFourier can be indexed directly like arrays. If given 1D indices, this is the index into the vector representation:

f
4-element 2×2-pixel 1.0′-resolution LambertMap{Array{Float64, 2}}:
 0.4599124117435962
 0.6967131274166374
 0.9661685855933742
 0.5162432729683115
f[1], f[2], f[3], f[4]
(0.4599124117435962, 0.6967131274166374, 0.9661685855933742, 0.5162432729683115)
f[5]
BoundsError: attempt to access 2×2 Matrix{Float64} at index [5]



Stacktrace:

 [1] getindex

   @ ./array.jl:801 [inlined]

 [2] getindex(f::LambertMap{Array{Float64, 2}}, I::Int64)

   @ CMBLensing ~/CMBLensing/src/base_fields.jl:31

 [3] top-level scope

   @ In[24]:1

 [4] eval

   @ ./boot.jl:360 [inlined]

 [5] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)

   @ Base ./loading.jl:1094

Or with a 2D index, this indexes directly into the 2D map:

f[1,1], f[2,1], f[1,2], f[2,2]
(0.4599124117435962, 0.6967131274166374, 0.9661685855933742, 0.5162432729683115)

Note: there is no overhead to indexing f in this way as compared to working directly on the underlying array.

For other fields which are built on FieldTuples, 1D indexing will instead index the tuple indices:

ft
8-element Field-(a,b)-Tuple{LambertMap{Array{Float64, 2}}, LambertMap{Array{Float64, 2}}}:
 0.09914786391952712
 0.18598410611302874
 0.6349936166558303
 0.6142728589097468
 0.9058994959088245
 0.5481505048566739
 0.6767661193126735
 0.6113430607196322
ft[1]
4-element 2×2-pixel 1.0′-resolution LambertMap{Array{Float64, 2}}:
 0.09914786391952712
 0.18598410611302874
 0.6349936166558303
 0.6142728589097468
ft[2]
4-element 2×2-pixel 1.0′-resolution LambertMap{Array{Float64, 2}}:
 0.9058994959088245
 0.5481505048566739
 0.6767661193126735
 0.6113430607196322
ft[3]
BoundsError: attempt to access NamedTuple{(:a, :b), Tuple{LambertMap{Array{Float64, 2}}, LambertMap{Array{Float64, 2}}}} at index [3]



Stacktrace:

 [1] getindex

   @ ./namedtuple.jl:117 [inlined]

 [2] getindex(f::Field-(a,b)-Tuple{LambertMap{Array{Float64, 2}}, LambertMap{Array{Float64, 2}}}, i::Int64)

   @ CMBLensing ~/CMBLensing/src/field_tuples.jl:31

 [3] top-level scope

   @ In[29]:1

 [4] eval

   @ ./boot.jl:360 [inlined]

 [5] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)

   @ Base ./loading.jl:1094

To get the underlying data arrays, use the object's properties:

f.Ix
2×2 Matrix{Float64}:
 0.459912  0.966169
 0.696713  0.516243

You can always find out what properties are available by typing f.<Tab>. For example, if you typed ft then hit <Tab> you'd get:

ft |> propertynames
(:fs, :a, :b)

For a FieldTuple like the FlatQUMap object, fqu, you can get each individual Q or U field:

fqu.Q
4-element 2×2-pixel 1.0′-resolution LambertMap{SubArray{Float64, 2, Array{Float64, 3}, Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Int64}, true}}:
 0.09914786391952712
 0.18598410611302874
 0.6349936166558303
 0.6142728589097468

Or fqu.Qx which is shorthand for fqu.Q.Ix:

fqu.Q.Ix === fqu.Qx
true

If you convert f to Fourier space, it would have the Il property to get the Fourier coefficients of the $I$ component:

Fourier(f).Il
2×2 Matrix{ComplexF64}:
  2.63904+0.0im  -0.325786+0.0im
 0.213125+0.0im  -0.686726+0.0im

For convenience, you can index fields with brackets [] and any necessary conversions will be done automatically:

f[:Il]
2×2 Matrix{ComplexF64}:
  2.63904+0.0im  -0.325786+0.0im
 0.213125+0.0im  -0.686726+0.0im

This works between any bases. For example. fqu is originally QUMap but we can convert to EBFourier and get the El coefficients:

fqu[:El]
2×2 view(::Array{ComplexF64, 3}, :, :, 1) with eltype ComplexF64:
    -1.5344-0.0im   0.964135+0.0im
 -0.0661155+0.0im  -0.292326+0.0im

The general rule to keep in mind for these two ways of accessing the underlying data is:

  • Properties (i.e. f.Ix) are type-stable and get you the underlying data arrays, even recursively from special FieldTuples like FlatQUMap, etc... If these arrays are modified, they affect the original field.
  • Indices (i.e. f[:Ix]) are not type-stable, and may or may not be one of the underlying data arrays (because a basis conversion may have been performed). They should be used for getting (not setting) data, and in non-performance-critical code.