Skip to content

v0.5.0
Compare
Choose a tag to compare
@github-actions github-actions released this 25 Jan 23:42
· 442 commits to master since this release
v0.5.0
10cd697

CMBLensing v0.5.0

Diff since v0.4.1

Closed issues:

  • Make Nside a tuple by default (#40)

Merged pull requests:

Description

The main thing this release updates is to change the low-level struct representing a Field, making it cleaner, more efficient, and more extendable. Its not meant to break or change any high level code. Things like FlatMap(rand(10,10), θpix=3), f.Ix, ∇*f, etc... all should work the same. Similarly, code in "high-level" functions like MAP_marg, MAP_joint, sample_joint, the LenseFlow velocities, etc... is mostly or completely unchanged. However, some "low-level" code contains breaking changes, hence the version number is increased as a breaking release.

Some cool things new things:

  • Compatible with Julia 1.6-beta1

  • GPU polarization calculations are sped up by ~30% (by virtue of Q/U and E/B being stored in a single contiguous array, thus requiring less GPU kernel calls to operate on it).

  • Precomilation is no longer triggered every time you change Nside or θpix. So after the first (slow) load_sim(Nside=128), subsequent load_sim(Nside=129), load_sim(Nside=130), etc... will all be fast. This also opens up the possibility of adding some precompiling directives so the compilation price is only paid when precompiling CMBLensing, rather than every time you load it (this is for a future PR, but wouldn't have made sense with the old system where you would have needed to precompile every single Nside and θpix one might imagine using).

  • Automatic differentiation is much more robust, correct, and works in more cases. It should now always be possible to "drop down" to the underlying arrays and work with those, opening up the full range of everything that is differentiable with normal arrays in Julia. E.g.:

    julia> f = FlatMap(rand(10,10));

julia> gradient(f -> prod(sum(sin.(f.arr) .+ 1, dims=1)), f)[1]
100-element 10×10-pixel 1.0′-resolution FlatMap{Array{Float64,2},ProjLambert{Float64}}:
2.580277728774914e10
2.6926717106096638e10
⋮
1.9684716048961834e10
2.21171314125531e10

  • You no longer need to using CUDA before using CMBLensing, you can load CUDA later and it works fine.

  • Fields can hold anything, not just Reals, including e.g. Unitful quantities. The entire system isn't (yet) setup to handle units, but several things work:

    julia> using Unitful

julia> f = FlatMap(rand(10,10) .* u"μK")
100-element 10×10-pixel 1.0′-resolution FlatMap{Array{Quantity{Float64,𝚯,Unitful.FreeUnits{(μK,),𝚯,nothing}},2},ProjLambert{Float64}}:
0.5814487830649309 μK
0.02371155697423344 μK
         ⋮
0.15031100142211362 μK
0.5543493616436406 μK

julia> f*f
100-element 10×10-pixel 1.0′-resolution FlatMap{Array{Quantity{Float64,𝚯^2,Unitful.FreeUnits{(μK^2,),𝚯^2,nothing}},2},ProjLambert{Float64}}:
0.33808268732768904 μK^2
0.33619069438247673 μK^2
         ⋮
0.022593397148518643 μK^2
0.3073032147547118 μK^2

julia> f+f
100-element 10×10-pixel 1.0′-resolution FlatMap{Array{Quantity{Float64,𝚯,Unitful.FreeUnits{(μK,),𝚯,nothing}},2},ProjLambert{Float64}}:
1.1628975661298617 μK
0.04742311394846688 μK
         ⋮
0.30062200284422724 μK
1.1086987232872811 μK

    Note that there is no overhead for unitful arrays compared to simple float arrays (since Julia stores the dimension information in the type).

  • You no longer need fieldinfo(f) to access field metadata like Ny, ℓx, etc... You can just directly access f.Ny, f.ℓx, etc... fieldinfo is left for backwards compatibilty and is just a no-op, and will be removed in a future version.

Under the hood:

  • The FlatMap (and similar for FlatFourier) type used to look like

    struct FlatMap{P<:Flat,T<:Real,M<:AbstractRank2or3Array{T}} <: Field{Map,S0,P,T}
    Ix :: M
end

    where Flat held some info about the pixelization (ie "metadata" about the field):

    abstract type Flat{Nside,θpix,∂mode<:∂modes,D} <: Pix end

    and Ix was an Ny×Nx(×Nbatch)-sized array.

    Polarization fields were FieldTuples, e.g. a FlatQUMap was FieldTuple of two FlatS0s, etc...

  • The new basic structure to hold any kind of field is:

    struct BaseField{B, M, T, A<:AbstractArray{T}} <: Field{B, T}
    arr :: A
    metadata :: M
end

    where B is a Basis, and metadata is any arbitrary metadata we want to carry around with the fields. For flat maps, the array is a Ny×Nx(xNpol×Nbatch)-sized array, where Npol is 1, 2, or 3, for I, P, or IP polarization. This means FlatS2 and FlatS02 are no longer FieldTuples of FlatS0s, they're their own single struct and the I/Q/U data is stored in one contiguous array.

    The fundamental difference basically is that before, all the field "metadata" had to go in the type, whereas now some if can go in the f.metadata field. This makes it easier to add new arbitrary metadata, and makes it so we don't trigger recompilation every time the metadata changes, only when the type of the metadata changes (e.g. changing Nside=10 to Nside=11 does not change the type).

    Broadcasting and all the usual stuff is implemented efficietly, such that there continues to be no overhead for operating on Fields instead of the underlying arrays.

Assets 2
Loading