From 0c06686c2ff58552564b9e6b3ee0de36ee5c5020 Mon Sep 17 00:00:00 2001 From: Maurizio Tomasi Date: Tue, 28 Apr 2020 09:36:48 +0200 Subject: [PATCH 1/4] Generalize the definition of Map and Alm --- src/alm.jl | 51 ++++++++++---- src/conformables.jl | 43 +++++++----- src/map.jl | 97 +++++++++++++++++++-------- src/map_pixelfunc.jl | 32 +++++---- src/polarizedmap.jl | 67 ++++++++++++------- src/projections.jl | 30 ++++----- src/sphtfunc.jl | 135 ++++++++++++++++++++++++-------------- test/test_mapio.jl | 4 +- test/test_polarizedmap.jl | 2 +- 9 files changed, 303 insertions(+), 158 deletions(-) diff --git a/src/alm.jl b/src/alm.jl index f596599..3eaa36d 100644 --- a/src/alm.jl +++ b/src/alm.jl @@ -1,24 +1,53 @@ # Definition of the composite type Alm -"An array of a_ℓm numbers." -mutable struct Alm{T <: Number} - alm::Array{T, 1} +"""An array of harmonic coefficients (a_ℓm). + +The type `T` is used for the value of each harmonic coefficient, and +it must be a `Number` (one should however only use complex types for +this). The type `AA` is used to store the array of coefficients; a +typical choice is `Vector`. + +A `Alm` type contains the following fields: + +- `alm`: the array of harmonic coefficients +- `lmax`: the maximum value for ``ℓ`` +- `mmax`: the maximum value for ``m`` +- ``tval`: maximum number of ``m`` coefficients for the maximum ``ℓ`` + + +""" +mutable struct Alm{T <: Number, AA <: AbstractArray{T, 1}} + alm::AA lmax::Int mmax::Int tval::Int - Alm{T}(lmax, mmax) where {T <: Number} = new(zeros(T, numberOfAlms(lmax, mmax)), - lmax, - mmax, - 2lmax + 1) - - function Alm{T}(lmax, mmax, arr::Array{T, 1}) where {T <: Number} + Alm{T, AA}( + lmax, + mmax, + ) where {T <: Number, AA <: AbstractArray{T, 1}} = + new{T, AA}( + zeros(T, numberOfAlms(lmax, mmax)), + lmax, + mmax, + 2lmax + 1, + ) + + function Alm{T, AA}( + lmax, + mmax, + arr::AA, + ) where {T <: Number, AA <: AbstractArray{T, 1}} (numberOfAlms(lmax, mmax) == length(arr)) || throw(DomainError()) - new{T}(arr, lmax, mmax, 2lmax + 1) + new{T, AA}(arr, lmax, mmax, 2lmax + 1) end end - + +Alm{T}(lmax, mmax) where {T <: Number} = Alm{T, Array{T, 1}}(lmax, mmax) +Alm{T}(lmax, mmax, arr) where {T <: Number} = + Alm{T, Array{T, 1}}(lmax, mmax, collect(arr)) + ################################################################################ """ diff --git a/src/conformables.jl b/src/conformables.jl index e0606b3..4cd8477 100644 --- a/src/conformables.jl +++ b/src/conformables.jl @@ -1,32 +1,45 @@ ################################################################################ -conformables(map1::Map{T, RingOrder}, map2::Map{S, RingOrder}) where {T, S} = +conformables( + map1::Map{T, RingOrder, AA1}, + map2::Map{S, RingOrder, AA2}, +) where {T, S, AA1, AA2} = map1.resolution.nside == map2.resolution.nside -conformables(map1::Map{T, NestedOrder}, map2::Map{S, NestedOrder}) where {T, S} = +conformables( + map1::Map{T, NestedOrder, AA1}, + map2::Map{S, NestedOrder, AA2}, +) where {T, S, AA1, AA2} = map1.resolution.nside == map2.resolution.nside -conformables(map1::Map{T, O1}, - map2::Map{S, O2}) where {T, S, O1 <: Order, O2 <: Order} = false +conformables(map1::Map{T, O1, AA1}, + map2::Map{S, O2, AA2}) where {T, S, O1, O2, AA1, AA2} = false ################################################################################ -conformables(map1::PolarizedMap{T, RingOrder}, map2::PolarizedMap{S, RingOrder}) where {T, S} = +conformables( + map1::PolarizedMap{T, RingOrder, AA1}, + map2::PolarizedMap{S, RingOrder, AA2}, +) where {T, S, AA1, AA2} = map1.i.resolution.nside == map2.i.resolution.nside -conformables(map1::PolarizedMap{T, NestedOrder}, map2::PolarizedMap{S, NestedOrder}) where {T, S} = +conformables( + map1::PolarizedMap{T, NestedOrder, AA1}, + map2::PolarizedMap{S, NestedOrder, AA2}, +) where {T, S, AA1, AA2} = map1.i.resolution.nside == map2.i.resolution.nside -conformables(map1::PolarizedMap{T, O1}, - map2::PolarizedMap{S, O2}) where {T, S, O1 <: Order, O2 <: Order} = false +conformables(map1::PolarizedMap{T, O1, AA1}, + map2::PolarizedMap{S, O2, AA2}) where {T, S, O1, O2, AA1, AA2} = false """ - conformables{T, S, O1 <: Order, O2 <: Order}(map1::Map{T, O1}, - map2::Map{S, O2}) -> Bool - conformables{T, S, O1 <: Order, O2 <: Order}(map1::PolarizedMap{T, O1}, - map2::PolarizedMap{S, O2}) -> Bool - -Determine if two Healpix maps are "conformables", i.e., if their -shape and ordering are the same. + conformables{T, S, O1, O2}(map1::Map{T, O1, AA1}, + map2::Map{S, O2, AA2}) -> Bool + conformables{T, S, O1, O2}(map1::PolarizedMap{T, O1, AA1}, + map2::PolarizedMap{S, O2, AA2}) -> Bool + +Determine if two Healpix maps are "conformables", i.e., if their shape +and ordering are the same. The array types `AA1` and `AA2` are not considered +in testing conformability. """ conformables diff --git a/src/map.jl b/src/map.jl index f830b26..27ac975 100644 --- a/src/map.jl +++ b/src/map.jl @@ -27,20 +27,50 @@ don't care about the ordering of a map.""" abstract type GenericMap{T} <: AbstractArray{T, 1} end """ - struct Map{T, O <: Order} <: GenericMap{T} + struct Map{T, O <: Order, AA <: AbstractArray{T, 1}} <: GenericMap{T} -A Healpix map. The type `T` is used for the value of the pixels in -a map, and it can be anything (even a string!). The type `O` is used -to specify the ordering of the pixels, and it can either be -`RingOrder` or `NestedOrder`. +A Healpix map. The type `T` is used for the value of the pixels in a +map, and it can be anything (even a string!). The type `O` is used to +specify the ordering of the pixels, and it can either be `RingOrder` +or `NestedOrder`. The type `AA` is used to store the array of pixels; +typical types are `Vector`, `CUArray`, `SharedArray`, etc. A `Map` type contains the following fields: - `pixels`: array of pixels - `resolution`: instance of a `Resolution` object +You can construct a map using one of the following forms: + +- `Map{T, O, AA}(arr)` and `Map{T, O, AA}(nside::Number)` will use + `AA` as basetype + +- `Map{T, O}(arr)` and `Map{T, O}(nside::Number)` will use `Array{T, + 1}` as basetype + +# Examples + +The following example creates a map with `NSIDE=32` in `RING` order, +containing integer values starting from 1: + + mymap = Healpix.Map{Int64, Healpix.RingOrder}(1:Healpix.nside2npix(32)) + +The call to `collect` is required to convert the range in an array. + +This example creates a map in `NESTED` order, with `NSIDE=64`, filled +with zeroes: + + mymap = Healpix.Map{Float64, Healpix.NestedOrder}(64) + +Finally, the following examples show how to use `SharedArray`: + + using SharedArrays + mymap = Healpix.Map{Float64, Healpix.NestedOrder, SharedArray{Float64, 1}}(64) + + pixels = SharedArray{Int64, 1}(1:12 |> collect) + mymap = Healpix.Map{Int64, Healpix.RingOrder, SharedArray{Int64, 1}}(m) """ -mutable struct Map{T, O <: Order} <: GenericMap{T} +mutable struct Map{T, O <: Order, AA <: AbstractArray{T, 1}} <: GenericMap{T} pixels::Array{T,1} resolution::Resolution @@ -49,48 +79,63 @@ mutable struct Map{T, O <: Order} <: GenericMap{T} Create an empty map with the specified NSIDE. """ - Map{T, O}(nside::Number) where {T, O <: Order} = new(zeros(T, nside2npix(nside)), - Resolution(nside)) + function Map{T, O, AA}(nside::Number) where {T, O <: Order, AA <: AbstractArray{T, 1}} + new(zeros(T, nside2npix(nside)), Resolution(nside)) + end """ - Create a map with the specified array of pixels. + Initialize a map from a generic array """ - function Map{T, O}(arr::Array{T,1}) where {T, O <: Order} + function Map{T, O, AA}(arr::AA) where {T, O <: Order, AA <: AbstractArray{T, 1}} nside = npix2nside(length(arr)) new(arr, Resolution(nside)) end end +"""Convenience function that uses `Array{T, 1}` as the type used to +hold the vector of pixels. +""" +function Map{T, O}(nside::Number) where {T, O <: Order} + Map{T, O, Array{T, 1}}(nside) +end + +"""Convenience function that uses `Array{T, 1}` as the type used to +hold the vector of pixels. +""" +function Map{T, O}(arr) where {T, O <: Order} + Map{T, O, Array{T, 1}}(collect(arr)) +end + import Base: +, -, *, / -+(a::Map{T,O}, b::Map{T,O}) where {T <: Number, O} = Map{T, O}(a.pixels .+ b.pixels) --(a::Map{T,O}, b::Map{T,O}) where {T <: Number, O} = Map{T, O}(a.pixels .- b.pixels) -*(a::Map{T,O}, b::Map{T,O}) where {T <: Number, O} = Map{T, O}(a.pixels .* b.pixels) -/(a::Map{T,O}, b::Map{T,O}) where {T <: Number, O} = Map{T, O}(a.pixels ./ b.pixels) ++(a::Map{T,O,AA}, b::Map{T,O,AA}) where {T <: Number, O, AA} = Map{T, O, AA}(a.pixels .+ b.pixels) +-(a::Map{T,O,AA}, b::Map{T,O,AA}) where {T <: Number, O, AA} = Map{T, O, AA}(a.pixels .- b.pixels) +*(a::Map{T,O,AA}, b::Map{T,O,AA}) where {T <: Number, O, AA} = Map{T, O, AA}(a.pixels .* b.pixels) +/(a::Map{T,O,AA}, b::Map{T,O,AA}) where {T <: Number, O, AA} = Map{T, O, AA}(a.pixels ./ b.pixels) -+(a::Map{T,O}, b::Number) where {T <: Number, O} = Map{T, O}(a.pixels .+ b) --(a::Map{T,O}, b::Number) where {T <: Number, O} = a + (-b) -*(a::Map{T,O}, b::Number) where {T <: Number, O} = Map{T, O}(a.pixels .* b) -/(a::Map{T,O}, b::Number) where {T <: Number, O} = Map{T, O}(a.pixels ./ b) ++(a::Map{T,O,AA}, b::Number) where {T <: Number, O, AA} = Map{T, O, AA}(a.pixels .+ b) +-(a::Map{T,O,AA}, b::Number) where {T <: Number, O, AA} = a + (-b) +*(a::Map{T,O,AA}, b::Number) where {T <: Number, O, AA} = Map{T, O, AA}(a.pixels .* b) +/(a::Map{T,O,AA}, b::Number) where {T <: Number, O, AA} = Map{T, O, AA}(a.pixels ./ b) -+(a::Number, b::Map{T,O}) where {T <: Number, O} = b + a --(a::Number, b::Map{T,O}) where {T <: Number, O} = b + (-a) -*(a::Number, b::Map{T,O}) where {T <: Number, O} = b * a -/(a::Number, b::Map{T,O}) where {T <: Number, O} = Map{T, O}(a ./ b.pixels) ++(a::Number, b::Map{T,O,AA}) where {T <: Number, O, AA} = b + a +-(a::Number, b::Map{T,O,AA}) where {T <: Number, O, AA} = b + (-a) +*(a::Number, b::Map{T,O,AA}) where {T <: Number, O, AA} = b * a +/(a::Number, b::Map{T,O,AA}) where {T <: Number, O, AA} = Map{T, O, AA}(a ./ b.pixels) ################################################################################ # Iterator interface -Base.size(m::Map{T, O}) where {T, O} = (m.resolution.numOfPixels,) +Base.size(m::Map{T, O, AA}) where {T, O, AA} = (m.resolution.numOfPixels,) -Base.IndexStyle(::Type{<:Map{T, O}}) where {T, O} = IndexLinear() +Base.IndexStyle(::Type{<:Map{T, O, AA}}) where {T, O, AA} = IndexLinear() -function getindex(m::Map{T, O}, i::Integer) where {T, O} +function getindex(m::Map{T, O, AA}, i::Integer) where {T, O, AA} 1 ≤ i ≤ m.resolution.numOfPixels || throw(BoundsError(m, i)) m.pixels[i] end -function setindex!(m::Map{T, O}, val, i::Integer) where {T, O} +function setindex!(m::Map{T, O, AA}, val, i::Integer) where {T, O, AA} 1 ≤ i ≤ m.resolution.numOfPixels || throw(BoundsError(m, i)) m.pixels[i] = val end diff --git a/src/map_pixelfunc.jl b/src/map_pixelfunc.jl index fe5be4e..474d1ca 100644 --- a/src/map_pixelfunc.jl +++ b/src/map_pixelfunc.jl @@ -1,15 +1,17 @@ ################################################################################ -ang2pix(map::Map{T, RingOrder}, theta, phi) where {T} = ang2pixRing(map.resolution, theta, phi) -ang2pix(map::Map{T, NestedOrder}, theta, phi) where {T} = ang2pixNest(map.resolution, theta, phi) -ang2pix(map::PolarizedMap{T, RingOrder}, theta, phi) where {T} = +ang2pix(map::Map{T, RingOrder, AA}, theta, phi) where {T, AA} = + ang2pixRing(map.resolution, theta, phi) +ang2pix(map::Map{T, NestedOrder, AA}, theta, phi) where {T, AA} = + ang2pixNest(map.resolution, theta, phi) +ang2pix(map::PolarizedMap{T, RingOrder, AA}, theta, phi) where {T, AA} = ang2pixRing(map.i.resolution, theta, phi) -ang2pix(map::PolarizedMap{T, NestedOrder}, theta, phi) where {T} = +ang2pix(map::PolarizedMap{T, NestedOrder, AA}, theta, phi) where {T, AA} = ang2pixNest(map.i.resolution, theta, phi) @doc raw""" - ang2pix{T, O <: Order}(map::Map{T, O}, theta, phi) - ang2pix{T, O <: Order}(map::PolarizedMap{T, O}, theta, phi) + ang2pix{T, O, AA}(map::Map{T, O}, theta, phi) + ang2pix{T, O, AA}(map::PolarizedMap{T, O}, theta, phi) Convert the direction specified by the colatitude `theta` (∈ [0, π]) and the longitude `phi` (∈ [0, 2π]) into the index of the pixel in the @@ -19,11 +21,13 @@ ang2pix ################################################################################ -pix2ang(map::Map{T, RingOrder}, ipix) where {T} = pix2angRing(map.resolution, ipix) -pix2ang(map::Map{T, NestedOrder}, ipix) where {T} = pix2angNest(map.resolution, ipix) -pix2ang(map::PolarizedMap{T, RingOrder}, ipix) where {T} = +pix2ang(map::Map{T, RingOrder, AA}, ipix) where {T, AA} = + pix2angRing(map.resolution, ipix) +pix2ang(map::Map{T, NestedOrder, AA}, ipix) where {T, AA} = + pix2angNest(map.resolution, ipix) +pix2ang(map::PolarizedMap{T, RingOrder, AA}, ipix) where {T, AA} = pix2angRing(map.i.resolution, ipix) -pix2ang(map::PolarizedMap{T, NestedOrder}, ipix) where {T} = +pix2ang(map::PolarizedMap{T, NestedOrder, AA}, ipix) where {T, AA} = pix2angNest(map.i.resolution, ipix) @doc raw""" @@ -39,7 +43,7 @@ pix2ang ################################################################################ # Interpolation -function interpolate(m::Map{T, RingOrder}, θ, ϕ, pixbuf, weightbuf) where {T} +function interpolate(m::Map{T, RingOrder, AA}, θ, ϕ, pixbuf, weightbuf) where {T, AA} getinterpolRing(m.resolution, θ, ϕ, pixbuf, weightbuf) result = zero(weightbuf[1]) @@ -50,7 +54,7 @@ function interpolate(m::Map{T, RingOrder}, θ, ϕ, pixbuf, weightbuf) where {T} result end -function interpolate(m::Map{T, RingOrder}, θ, ϕ) where {T} +function interpolate(m::Map{T, RingOrder, AA}, θ, ϕ) where {T, AA} pixbuf = Array{Int}(undef, 4) weightbuf = Array{Float64}(undef, 4) @@ -58,8 +62,8 @@ function interpolate(m::Map{T, RingOrder}, θ, ϕ) where {T} end """ - interpolate(m::Map{T, RingOrder}, θ, ϕ) -> Value - interpolate(m::Map{T, RingOrder}, θ, ϕ, pixbuf, weightbuf) -> Value + interpolate(m::Map{T, RingOrder, AA}, θ, ϕ) -> Value + interpolate(m::Map{T, RingOrder, AA}, θ, ϕ, pixbuf, weightbuf) -> Value Return an interpolated value of the map along the specified direction. diff --git a/src/polarizedmap.jl b/src/polarizedmap.jl index 8bf7467..b13adfe 100644 --- a/src/polarizedmap.jl +++ b/src/polarizedmap.jl @@ -11,7 +11,7 @@ when you don't care about the ordering of a map. abstract type GenericPolarizedMap{T} end @doc raw""" - mutable struct PolarizedMap{T, O <: Healpix.Order} + mutable struct PolarizedMap{T, O <: Healpix.Order, AA <: AbstractArray{T, 1}} A polarized I/Q/U map. It contains three Healpix maps with the same NSIDE: @@ -22,21 +22,21 @@ A polarized I/Q/U map. It contains three Healpix maps with the same NSIDE: You can create an instance of this type using the function [`PolarizedMap{T,O}`](@ref), which comes in three flavours: -- `PolarizedMap(i::Map{T,O}, q::Map{T,O}, u::Map{T,O})` +- `PolarizedMap(i::Map{T,O,AA}, q::Map{T,O,AA}, u::Map{T,O,AA})` - `PolarizedMap{T,O}(i::AbstractVector{T}, q::AbstractVector{T}, u::AbstractVector{T})` - `PolarizedMap{T,O}(nside::Number)` """ -mutable struct PolarizedMap{T, O <: Order} <: GenericPolarizedMap{T} - i::Map{T,O} - q::Map{T,O} - u::Map{T,O} +mutable struct PolarizedMap{T, O <: Order, AA <: AbstractArray{T, 1}} <: GenericPolarizedMap{T} + i::Map{T,O,AA} + q::Map{T,O,AA} + u::Map{T,O,AA} - function PolarizedMap{T, O}( - i::Map{T, O}, - q::Map{T, O}, - u::Map{T, O}, - ) where {T, O <: Order} + function PolarizedMap{T, O, AA}( + i::Map{T, O, AA}, + q::Map{T, O, AA}, + u::Map{T, O, AA}, + ) where {T, O <: Order, AA <: AbstractArray{T, 1}} ((length(i) != length(q)) || (length(i) != length(q))) && throw( ArgumentError("The three I/Q/U maps must have the same resolution"), ) @@ -44,24 +44,41 @@ mutable struct PolarizedMap{T, O <: Order} <: GenericPolarizedMap{T} new(i, q, u) end - function PolarizedMap{T, O}( - i::AbstractVector{T}, - q::AbstractVector{T}, - u::AbstractVector{T}, - ) where {T, O <: Order} + function PolarizedMap{T, O, AA}( + i::AA, + q::AA, + u::AA, + ) where {T, O <: Order, AA <: AbstractArray{T, 1}} + ((length(i) != length(q)) || (length(i) != length(q))) && throw( + ArgumentError("The three I/Q/U vectors must have the same resolution"), + ) - PolarizedMap{T, O}( - Map{T, O}(i), - Map{T, O}(q), - Map{T, O}(u), + new( + Map{T, O, AA}(i), + Map{T, O, AA}(q), + Map{T, O, AA}(u), ) end - function PolarizedMap{T, O}(nside::Number) where {T, O <: Order} - PolarizedMap{T, O}( - Map{T, O}(nside), - Map{T, O}(nside), - Map{T, O}(nside), + function PolarizedMap{T, O, AA}( + nside::Number + ) where {T, O <: Order, AA <: AbstractArray{T, 1}} + new( + Map{T, O, AA}(nside), + Map{T, O, AA}(nside), + Map{T, O, AA}(nside), ) end end + +function PolarizedMap{T, O}(nside::Number) where {T, O <: Order} + PolarizedMap{T, O, Array{T, 1}}(nside) +end + +function PolarizedMap{T, O}( + i::Array{T, 1}, + q::Array{T, 1}, + u::Array{T, 1}, +) where {T, O <: Order} + PolarizedMap{T, O, Array{T, 1}}(i, q, u) +end diff --git a/src/projections.jl b/src/projections.jl index f3523fd..d787dc2 100644 --- a/src/projections.jl +++ b/src/projections.jl @@ -12,7 +12,7 @@ import RecipesBase const UNSEEN = -1.6375e+30 """ - project(invprojfn, m::Map{T, O}, bmpwidth, bmpheight; kwargs...) where {T <: Number, O <: Order} + project(invprojfn, m::Map{T, O, AA}, bmpwidth, bmpheight; kwargs...) Return a 2D bitmap (array) containing a cartographic projection of the map and a 2D bitmap containing a boolean mask. The size of the bitmap @@ -32,8 +32,8 @@ Return a `Array{Union{Missing, Float32}}` containing the intensity of each pixel. Pixels falling outside the projection are marked as NaN, and unseen pixels are marked as `missing`. """ -function project(invprojfn, m::Map{T,O}, bmpwidth, bmpheight, - projparams = Dict()) where {T <: Number,O <: Healpix.Order} +function project(invprojfn, m::Map{T, O, AA}, bmpwidth, bmpheight, + projparams = Dict()) where {T <: Number, O, AA} center = get(projparams, :center, (0, 0)) unseen = get(projparams, :unseen, UNSEEN) @@ -259,18 +259,18 @@ end ################################################################################ """ - equirectangular(m::Map{T,O}; kwargs...) where {T <: Number, O <: Order} + equirectangular(m::Map{T,O,AA}; kwargs...) where {T <: Number, O, AA} High-level wrapper around `project` for equirectangular projections. """ -function equirectangular(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order} +function equirectangular(m::Map{T, O, AA}, projparams = Dict()) where {T <: Number, O, AA} width = get(projparams, :width, 720) height = get(projparams, :height, width) project(equiprojinv, m, width, height, projparams) end """ - mollweide(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order} + mollweide(m::Map{T, O, AA}, projparams = Dict()) where {T <: Number, O, AA} High-level wrapper around `project` for Mollweide projections. @@ -280,7 +280,7 @@ The following parameters can be set in the `projparams` dictionary: - `height`: height of the image, in pixels; if not specified, it will be assumed to be equal to `width` """ -function mollweide(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order} +function mollweide(m::Map{T, O, AA}, projparams = Dict()) where {T <: Number, O, AA} width = get(projparams, :width, 720) height = get(projparams, :height, width ÷ 2) project(mollweideprojinv, m, width, height, projparams) @@ -299,7 +299,7 @@ The following parameters can be set in the `projparams` dictionary: - `center`: position of the pixel in the middle of the left globe (*latitude* and longitude). """ -function orthographic(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order} +function orthographic(m::Map{T, O, AA}, projparams = Dict()) where {T <: Number, O, AA} width = get(projparams, :width, 720) height = get(projparams, :height, width) ϕ0, λ0 = get(projparams, :center, (0, 0)) @@ -309,7 +309,7 @@ function orthographic(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: end """ - orthographic2(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order} + orthographic2(m::Map{T, O, AA}, projparams = Dict()) where {T <: Number, O, AA} High-level wrapper around `project` for stereo orthographic projections. @@ -321,7 +321,7 @@ The following parameters can be set in the `projparams` dictionary: - `center`: position of the pixel in the middle of the left globe (*latitude* and longitude). Default is (0, 0). """ -function orthographic2(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order} +function orthographic2(m::Map{T, O, AA}, projparams = Dict()) where {T <: Number, O, AA} width = get(projparams, :width, 720) height = get(projparams, :height, width ÷ 2) ϕ0, λ0 = get(projparams, :center, (0, 0)) @@ -331,7 +331,7 @@ ortho2inv(x, y, ϕ0, λ0) end """ - gnomonic(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order} + gnomonic(m::Map{T, O, AA}, projparams = Dict()) where {T <: Number, O, AA} High-level wrapper around `project` for gnomonic projections. @@ -355,7 +355,7 @@ The following parameters can be set in the `projparams` dictionary: plot(m, gnomonic, Dict(:fov_rad = deg2rad(1.5), :center = (0, 0, deg2rad(45)))) ```` """ -function gnomonic(m::Map{T,O}, projparams = Dict()) where {T <: Number,O <: Order} +function gnomonic(m::Map{T, O, AA}, projparams = Dict()) where {T <: Number, O, AA} width = get(projparams, :width, 720) height = get(projparams, :height, width) ϕ0, λ0, ψ0 = get(projparams, :center, (0, 0, 0)) @@ -367,9 +367,9 @@ end ################################################################################ -RecipesBase.@recipe function plot(m::Map{T,O}, +RecipesBase.@recipe function plot(m::Map{T, O, AA}, projection = mollweide, - projparams = Dict()) where {T <: Number,O <: Order} + projparams = Dict()) where {T <: Number, O, AA} img, mask, anymasked = projection(m, projparams) @@ -425,7 +425,7 @@ RecipesBase.@recipe function plot(m::Map{T,O}, end @doc raw""" - plot(m::Map{T,O}, projection = mollweide, projparams = Dict()) + plot(m::Map{T, O, AA}, projection = mollweide, projparams = Dict()) Draw a representation of the map, using some specific projection. The parameter `projection` must be a function returning the diff --git a/src/sphtfunc.jl b/src/sphtfunc.jl index adff12c..80d6e78 100644 --- a/src/sphtfunc.jl +++ b/src/sphtfunc.jl @@ -50,8 +50,8 @@ end """ - map2alm!(map::Map{Float64, RingOrder}, alm::Alm{ComplexF64}; niter::Integer=3) - map2alm!(map::PolarizedMap{Float64, RingOrder}, alm::Array{Alm{ComplexF64},1}; + map2alm!(map::Map{Float64, RingOrder, Array{Float64, 1}}, alm::Alm{ComplexF64, Array{ComplexF64, 1}}; niter::Integer=3) + map2alm!(map::PolarizedMap{Float64, RingOrder, Array{Float64, 1}}, alm::Array{Alm{ComplexF64, Array{ComplexF64, 1}},1}; niter::Integer=3) This function performs a spherical harmonic transform on the map and places the results @@ -61,14 +61,18 @@ done in-place. # Arguments - `map::Union{Map{Float64, RingOrder}, PolarizedMap{Float64, RingOrder}`: the map that to perform the spherical harmonic transform on. -- `alm::Alm{ComplexF64}`: the spherical harmonic coefficients to be written to. +- `alm::Alm{ComplexF64, Array{ComplexF64, 1}}`: the spherical harmonic coefficients to be written to. # Keywords - `niter::Integer`: number of iterations of SHTs to perform, to enhance accuracy """ function map2alm! end -function map2alm!(map::Map{Float64, RingOrder}, alm::Alm{ComplexF64}; niter::Integer=3) +function map2alm!( + map::Map{Float64, RingOrder, Array{Float64, 1}}, + alm::Alm{ComplexF64, Array{ComplexF64, 1}}; + niter::Integer = 3, +) geom_info = Libsharp.make_healpix_geom_info(map.resolution.nside, 1) alm_info = Libsharp.make_triangular_alm_info(alm.lmax, alm.mmax, 1) Libsharp.sharp_execute!( @@ -79,8 +83,11 @@ function map2alm!(map::Map{Float64, RingOrder}, alm::Alm{ComplexF64}; niter::Int end end -function map2alm!(map::PolarizedMap{Float64, RingOrder}, alm::Array{Alm{ComplexF64},1}; - niter::Integer=3) +function map2alm!( + map::PolarizedMap{Float64, RingOrder, Array{Float64, 1}}, + alm::Array{Alm{ComplexF64, Array{ComplexF64, 1}},1}; + niter::Integer = 3, +) geom_info = Libsharp.make_healpix_geom_info(map.i.resolution.nside, 1) alm_info = Libsharp.make_triangular_alm_info(alm[1].lmax, alm[1].mmax, 1) @@ -103,31 +110,36 @@ end """ - map2alm(map::Map{Float64, RingOrder}; + map2alm(map::Map{Float64, RingOrder, Array{Float64, 1}}; lmax=nothing, mmax=nothing, niter::Integer=3) - map2alm(m::Map{T, RingOrder}; lmax=nothing, mmax=nothing, + map2alm(m::Map{T, RingOrder, Array{T, 1}}; lmax=nothing, mmax=nothing, niter::Integer=3) where T <: Real -Compute the spherical harmonic coefficients of a map. To enhance precision, more iterations -of the transforms can be performed by passing a nonzero `niter`. The underlying SHT library -libsharp performs all calculations in Cdouble, so all inputs are converted to types based on -Float64. +Compute the spherical harmonic coefficients of a map. To enhance +precision, more iterations of the transforms can be performed by +passing a nonzero `niter`. The underlying SHT library libsharp +performs all calculations in Cdouble, so all inputs are converted to +types based on Float64. # Arguments -- `map::Union{Map{T, RingOrder}, PolarizedMap{T, RingOrder}`: the map that to - perform the spherical harmonic transform on. +- `map::Union{Map{T, RingOrder, AA}, PolarizedMap{T, RingOrder, AA}`: + the map that to perform the spherical harmonic transform on. # Keywords -- `lmax::Integer`: the maximum ℓ coefficient, will default to 3*nside-1 if not specified. +- `lmax::Integer`: the maximum ℓ coefficient, will default to + 3*nside-1 if not specified. - `mmax::Integer`: the maximum m coefficient -- `niter::Integer`: number of SHT iterations, to enhance precision. Defaults to 3 +- `niter::Integer`: number of SHT iterations, to enhance precision. + Defaults to 3 # Returns -- `Alm{ComplexF64}`: the spherical harmonic coefficients corresponding to the map +- `Alm{ComplexF64, Array{ComplexF64, 1}}`: the spherical harmonic + coefficients corresponding to the map + """ function map2alm end -function map2alm(map::Map{Float64, RingOrder}; +function map2alm(map::Map{Float64, RingOrder, Array{Float64, 1}}; lmax=nothing, mmax=nothing, niter::Integer=3) nside = map.resolution.nside @@ -140,7 +152,7 @@ function map2alm(map::Map{Float64, RingOrder}; return alm end -function map2alm(map::PolarizedMap{Float64, RingOrder}; +function map2alm(map::PolarizedMap{Float64, RingOrder, Array{Float64, 1}}; lmax=nothing, mmax=nothing, niter::Integer=3) nside = map.i.resolution.nside @@ -156,15 +168,23 @@ function map2alm(map::PolarizedMap{Float64, RingOrder}; end # convert maps to Float64 -function map2alm(map::Map{T, RingOrder}; lmax=nothing, mmax=nothing, - niter::Integer=3) where T <: Real +function map2alm( + map::Map{T, RingOrder, Array{T, 1}}; + lmax = nothing, + mmax = nothing, + niter::Integer = 3 +) where T <: Real map_float = Map{Float64, RingOrder}(convert(Array{Float64,1}, map.pixels)) return map2alm(map_float, lmax=lmax, mmax=mmax, niter=niter) end # convert PolarizedMap to Float64 -function map2alm(map::PolarizedMap{T, RingOrder}; lmax=nothing, mmax=nothing, - niter::Integer=3) where T <: Real +function map2alm( + map::PolarizedMap{T, RingOrder, Array{T, 1}}; + lmax = nothing, + mmax = nothing, + niter::Integer = 3, +) where T <: Real m_i = convert(Array{Float64,1}, map.i) m_q = convert(Array{Float64,1}, map.q) m_u = convert(Array{Float64,1}, map.u) @@ -174,23 +194,28 @@ end """ - alm2map!(alm::Alm{ComplexF64}, map::Map{Float64, RingOrder}) - alm2map!(alm::Array{Alm{ComplexF64},1}, map::PolarizedMap{Float64, RingOrder}) + alm2map!(alm::Alm{ComplexF64, Array{ComplexF64, 1}}, map::Map{Float64, RingOrder, Array{Float64, 1}}) + alm2map!(alm::Array{Alm{ComplexF64, Array{ComplexF64, 1}},1}, map::PolarizedMap{Float64, RingOrder, Array{Float64, 1}}) -This function performs a spherical harmonic transform on the map and places the results -in the passed `alm` object. This function requires types derived from Float64, since it is -done in-place. +This function performs a spherical harmonic transform on the map and +places the results in the passed `alm` object. This function requires +types derived from Float64, since it is done in-place. # Arguments -- `alm::Alm{ComplexF64}`: the spherical harmonic coefficients to perform the spherical - harmonic transform on. -- `map::Union{Map{Float64, RingOrder}, PolarizedMap{Float64, RingOrder}`: the map to be - written to. + +- `alm::Alm{ComplexF64, Array{ComplexF64, 1}}`: the spherical harmonic + coefficients to perform the spherical harmonic transform on. +- `map::Union{Map{Float64, RingOrder, Array{Float64, 1}}, + PolarizedMap{Float64, RingOrder, Array{Float64, 1}}`: the map to be written to. + """ function alm2map! end # in-place alm2map for spin-0 -function alm2map!(alm::Alm{ComplexF64}, map::Map{Float64, RingOrder}) +function alm2map!( + alm::Alm{ComplexF64, Array{ComplexF64, 1}}, + map::Map{Float64, RingOrder, Array{Float64, 1}}, +) geom_info = Libsharp.make_healpix_geom_info(map.resolution.nside, 1) alm_info = Libsharp.make_triangular_alm_info(alm.lmax, alm.mmax, 1) Libsharp.sharp_execute!( @@ -199,7 +224,10 @@ function alm2map!(alm::Alm{ComplexF64}, map::Map{Float64, RingOrder}) end # in-place alm2map for TEB to IQU -function alm2map!(alm::Array{Alm{ComplexF64},1}, map::PolarizedMap{Float64, RingOrder}) +function alm2map!( + alm::Array{Alm{ComplexF64, Array{ComplexF64, 1}},1}, + map::PolarizedMap{Float64, RingOrder, Array{Float64, 1}}, +) geom_info = Libsharp.make_healpix_geom_info(map.i.resolution.nside, 1) alm_info = Libsharp.make_triangular_alm_info(alm[1].lmax, alm[1].mmax, 1) @@ -216,30 +244,36 @@ end """ - alm2map(alm::Alm{ComplexF64}, nside::Integer) - alm2map(alm::Alm{T}, nside::Integer) where T - alm2map(alm::Array{Alm{ComplexF64},1}, nside::Integer) - alm2map(alms::Array{Alm{T},1}, nside::Integer) where T + alm2map(alm::Alm{ComplexF64, Array{Float64, 1}}, nside::Integer) + alm2map(alm::Alm{T, Array{T, 1}}, nside::Integer) where T + alm2map(alm::Array{Alm{ComplexF64, Array{ComplexF64, 1}},1}, nside::Integer) + alm2map(alms::Array{Alm{T, Array{T, 1}},1}, nside::Integer) where T -Compute a map from spherical harmonic coefficients. The underlying SHT library libsharp -performs all calculations in Cdouble, so all inputs are converted to types based on Float64. +Compute a map from spherical harmonic coefficients. The underlying SHT +library libsharp performs all calculations in Cdouble, so all inputs +are converted to types based on Float64. # Arguments -- `alm`: the spherical harmonic coefficients to transform. If of type `Alm{T}`, we assume a - spin-0 spherical harmonic transform. If an array of `Alm` is passed, we assume that - the components correspond to T, E, and B coefficients. + +- `alm`: the spherical harmonic coefficients to transform. If of type + `Alm{T}`, we assume a spin-0 spherical harmonic transform. If an + array of `Alm` is passed, we assume that the components correspond + to T, E, and B coefficients. # Keywords - `nside::Integer`: Healpix resolution parameter # Returns -- `Map{Float64, RingOrder}` or `PolarizedMap{Float64, RingOrder}` depending on if the input - alm is of type `Alm{T}` or `Array{Alm{T}}` respectively. + +- `Map{Float64, RingOrder, Array{Float64, 1}}` or + `PolarizedMap{Float64, RingOrder, Array{Float64, 1}}` depending on + if the input alm is of type `Alm{T, Array{T, 1}}` or `Array{Alm{T, + Array{T, 1}}}` respectively. """ function alm2map end # create a new set of spin-0 maps and project the coefficients to the map -function alm2map(alm::Alm{ComplexF64}, nside::Integer) +function alm2map(alm::Alm{ComplexF64, Array{ComplexF64, 1}}, nside::Integer) npix = nside2npix(nside) map = Map{Float64, RingOrder}(zeros(Float64, npix)) alm2map!(alm, map) @@ -247,7 +281,7 @@ function alm2map(alm::Alm{ComplexF64}, nside::Integer) end # create a new set of IQU maps and project the coefficients to the map -function alm2map(alm::Array{Alm{ComplexF64},1}, nside::Integer) +function alm2map(alm::Array{Alm{ComplexF64, Array{ComplexF64, 1}},1}, nside::Integer) npix = nside2npix(nside) map = PolarizedMap{Float64, RingOrder}( zeros(Float64, npix), zeros(Float64, npix), zeros(Float64, npix)) @@ -257,13 +291,16 @@ end # convert to ComplexF64 Alm for spin-0 if passed some other type function alm2map(alm::Alm{T}, nside::Integer) where T - alm_float = Alm{ComplexF64}( + alm_float = Alm{ComplexF64, Array{ComplexF64, 1}}( alm.lmax, alm.mmax, convert(Array{ComplexF64,1}, alm.alm)) return alm2map(alm_float, nside) end # convert to ComplexF64 Alm for TEB if passed some other type -function alm2map(alms::Array{Alm{Complex{T}},1}, nside::Integer) where T <: Real +function alm2map( + alms::Array{Alm{Complex{T}, Array{Complex{T}, 1}},1}, + nside::Integer, +) where T <: Real lmax = alms[1].lmax mmax = alms[1].mmax alm_t = Alm{ComplexF64}(lmax, mmax, diff --git a/test/test_mapio.jl b/test/test_mapio.jl index 86427ad..16d2131 100644 --- a/test/test_mapio.jl +++ b/test/test_mapio.jl @@ -1,12 +1,12 @@ # Map loading m = Healpix.readMapFromFITS("float_map.fits", 1, Float32) -@test typeof(m) == Healpix.Map{Float32,Healpix.RingOrder} +@test typeof(m) == Healpix.Map{Float32, Healpix.RingOrder, Array{Float32, 1}} @test m.resolution.nside == 4 @test m.pixels == [Float32(x) for x in 0:(12 * 4^2 - 1)] m = Healpix.readMapFromFITS("int_map.fits", 1, Int8) -@test typeof(m) == Healpix.Map{Int8,Healpix.RingOrder} +@test typeof(m) == Healpix.Map{Int8, Healpix.RingOrder, Array{Int8, 1}} @test m.resolution.nside == 1 @test m.pixels == [Int8(x) for x in 0:11] diff --git a/test/test_polarizedmap.jl b/test/test_polarizedmap.jl index 2da564f..89c5955 100644 --- a/test/test_polarizedmap.jl +++ b/test/test_polarizedmap.jl @@ -16,7 +16,7 @@ polmap = Healpix.PolarizedMap{Int8, Healpix.RingOrder}(128) @test length(polmap.i) == length(polmap.q) == length(polmap.u) @test length(polmap.i) == Healpix.nside2npix(128) -@test_throws ArgumentError("The three I/Q/U maps must have the same resolution") Healpix.PolarizedMap{Int, Healpix.RingOrder}( +@test_throws ArgumentError("The three I/Q/U vectors must have the same resolution") Healpix.PolarizedMap{Int, Healpix.RingOrder}( pixels_nside1, pixels_nside2, pixels_nside4, From 9f9a33df756d4075126222d7fbb9d177069d8e40 Mon Sep 17 00:00:00 2001 From: Maurizio Tomasi Date: Fri, 1 May 2020 08:04:31 +0200 Subject: [PATCH 2/4] Fix the definition of Map and map2alm This commit removes several trailing whitespaces as well. --- src/map.jl | 2 +- src/sphtfunc.jl | 54 ++++++++++++++++++++++++------------------------- 2 files changed, 28 insertions(+), 28 deletions(-) diff --git a/src/map.jl b/src/map.jl index 27ac975..ea9ffac 100644 --- a/src/map.jl +++ b/src/map.jl @@ -71,7 +71,7 @@ Finally, the following examples show how to use `SharedArray`: mymap = Healpix.Map{Int64, Healpix.RingOrder, SharedArray{Int64, 1}}(m) """ mutable struct Map{T, O <: Order, AA <: AbstractArray{T, 1}} <: GenericMap{T} - pixels::Array{T,1} + pixels::AA resolution::Resolution """ diff --git a/src/sphtfunc.jl b/src/sphtfunc.jl index 80d6e78..ce372a6 100644 --- a/src/sphtfunc.jl +++ b/src/sphtfunc.jl @@ -2,13 +2,13 @@ """ iterate_map2alm!( - maps::Array{Array{Float64,1},1}, alms::Array{Array{Complex{Float64},1},1}, - geom_info::Libsharp.GeomInfo, alm_info::Libsharp.AlmInfo, + maps::Array{Array{Float64,1},1}, alms::Array{Array{Complex{Float64},1},1}, + geom_info::Libsharp.GeomInfo, alm_info::Libsharp.AlmInfo, niter::Integer, spin::Integer) -This is an internal function for implementing iterative map2alm, used when the parameter -`niter` is greater than 0. It synthesizes a map from the alms, subtracts it from the map to -form a residual map, and then adds the harmonic coefficients of the residual map to the +This is an internal function for implementing iterative map2alm, used when the parameter +`niter` is greater than 0. It synthesizes a map from the alms, subtracts it from the map to +form a residual map, and then adds the harmonic coefficients of the residual map to the alms. It repeats this `niter` times. It performs this in-place on arrays of Float64 and Complex{Float64}. @@ -21,8 +21,8 @@ Complex{Float64}. - `spin::Integer`: spin of the field, 0 or 2 """ function iterate_map2alm!( - maps::Array{Array{Float64,1},1}, alms::Array{Array{Complex{Float64},1},1}, - geom_info::Libsharp.GeomInfo, alm_info::Libsharp.AlmInfo, + maps::Array{Array{Float64,1},1}, alms::Array{Array{Complex{Float64},1},1}, + geom_info::Libsharp.GeomInfo, alm_info::Libsharp.AlmInfo, niter::Integer, spin::Integer) ncomp = size(maps, 1) @@ -51,7 +51,7 @@ end """ map2alm!(map::Map{Float64, RingOrder, Array{Float64, 1}}, alm::Alm{ComplexF64, Array{ComplexF64, 1}}; niter::Integer=3) - map2alm!(map::PolarizedMap{Float64, RingOrder, Array{Float64, 1}}, alm::Array{Alm{ComplexF64, Array{ComplexF64, 1}},1}; + map2alm!(map::PolarizedMap{Float64, RingOrder, Array{Float64, 1}}, alm::Array{Alm{ComplexF64, Array{ComplexF64, 1}},1}; niter::Integer=3) This function performs a spherical harmonic transform on the map and places the results @@ -59,7 +59,7 @@ in the passed `alm` object. This function requires types derived from Float64, s done in-place. # Arguments -- `map::Union{Map{Float64, RingOrder}, PolarizedMap{Float64, RingOrder}`: the map that to +- `map::Union{Map{Float64, RingOrder}, PolarizedMap{Float64, RingOrder}`: the map that to perform the spherical harmonic transform on. - `alm::Alm{ComplexF64, Array{ComplexF64, 1}}`: the spherical harmonic coefficients to be written to. @@ -85,7 +85,7 @@ end function map2alm!( map::PolarizedMap{Float64, RingOrder, Array{Float64, 1}}, - alm::Array{Alm{ComplexF64, Array{ComplexF64, 1}},1}; + alm::Array{Alm{ComplexF64, Array{ComplexF64, 1}},1}; niter::Integer = 3, ) geom_info = Libsharp.make_healpix_geom_info(map.i.resolution.nside, 1) @@ -110,10 +110,10 @@ end """ - map2alm(map::Map{Float64, RingOrder, Array{Float64, 1}}; + map2alm(map::Map{Float64, RingOrder, AA}; lmax=nothing, mmax=nothing, niter::Integer=3) - map2alm(m::Map{T, RingOrder, Array{T, 1}}; lmax=nothing, mmax=nothing, - niter::Integer=3) where T <: Real + map2alm(m::Map{T, RingOrder, AA}; lmax=nothing, mmax=nothing, + niter::Integer=3) where {T <: Real, AA <: AbstractArray{T, 1} } Compute the spherical harmonic coefficients of a map. To enhance precision, more iterations of the transforms can be performed by @@ -139,7 +139,7 @@ types based on Float64. """ function map2alm end -function map2alm(map::Map{Float64, RingOrder, Array{Float64, 1}}; +function map2alm(map::Map{Float64, RingOrder, Array{Float64, 1}}; lmax=nothing, mmax=nothing, niter::Integer=3) nside = map.resolution.nside @@ -169,22 +169,22 @@ end # convert maps to Float64 function map2alm( - map::Map{T, RingOrder, Array{T, 1}}; + map::Map{T, RingOrder, AA}; lmax = nothing, - mmax = nothing, + mmax = nothing, niter::Integer = 3 -) where T <: Real +) where { T <: Real, AA <: AbstractArray{T, 1} } map_float = Map{Float64, RingOrder}(convert(Array{Float64,1}, map.pixels)) return map2alm(map_float, lmax=lmax, mmax=mmax, niter=niter) end # convert PolarizedMap to Float64 function map2alm( - map::PolarizedMap{T, RingOrder, Array{T, 1}}; + map::PolarizedMap{T, RingOrder, AA}; lmax = nothing, - mmax = nothing, + mmax = nothing, niter::Integer = 3, -) where T <: Real +) where { T <: Real, AA <: AbstractArray{T, 1} } m_i = convert(Array{Float64,1}, map.i) m_q = convert(Array{Float64,1}, map.q) m_u = convert(Array{Float64,1}, map.u) @@ -216,7 +216,7 @@ function alm2map!( alm::Alm{ComplexF64, Array{ComplexF64, 1}}, map::Map{Float64, RingOrder, Array{Float64, 1}}, ) - geom_info = Libsharp.make_healpix_geom_info(map.resolution.nside, 1) + geom_info = Libsharp.make_healpix_geom_info(map.resolution.nside, 1) alm_info = Libsharp.make_triangular_alm_info(alm.lmax, alm.mmax, 1) Libsharp.sharp_execute!( Libsharp.SHARP_ALM2MAP, 0, [alm.alm], [map.pixels], @@ -228,16 +228,16 @@ function alm2map!( alm::Array{Alm{ComplexF64, Array{ComplexF64, 1}},1}, map::PolarizedMap{Float64, RingOrder, Array{Float64, 1}}, ) - geom_info = Libsharp.make_healpix_geom_info(map.i.resolution.nside, 1) + geom_info = Libsharp.make_healpix_geom_info(map.i.resolution.nside, 1) alm_info = Libsharp.make_triangular_alm_info(alm[1].lmax, alm[1].mmax, 1) Libsharp.sharp_execute!( - Libsharp.SHARP_ALM2MAP, 0, + Libsharp.SHARP_ALM2MAP, 0, [alm[1].alm], [map.i.pixels], geom_info, alm_info, Libsharp.SHARP_DP) Libsharp.sharp_execute!( - Libsharp.SHARP_ALM2MAP, 2, + Libsharp.SHARP_ALM2MAP, 2, [alm[2].alm, alm[3].alm], [map.q.pixels, map.u.pixels], geom_info, alm_info, Libsharp.SHARP_DP) end @@ -303,11 +303,11 @@ function alm2map( ) where T <: Real lmax = alms[1].lmax mmax = alms[1].mmax - alm_t = Alm{ComplexF64}(lmax, mmax, + alm_t = Alm{ComplexF64}(lmax, mmax, convert(Array{ComplexF64,1}, alms[1].alm)) - alm_e = Alm{ComplexF64}(lmax, mmax, + alm_e = Alm{ComplexF64}(lmax, mmax, convert(Array{ComplexF64,1}, alms[2].alm)) - alm_b = Alm{ComplexF64}(lmax, mmax, + alm_b = Alm{ComplexF64}(lmax, mmax, convert(Array{ComplexF64,1}, alms[3].alm)) return alm2map([alm_t, alm_e, alm_b], nside) end From 3f24051e67feaf1f16d99d395d80413129ad859b Mon Sep 17 00:00:00 2001 From: Maurizio Tomasi Date: Fri, 1 May 2020 08:05:14 +0200 Subject: [PATCH 3/4] Improve the example for Map --- src/map.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/map.jl b/src/map.jl index ea9ffac..c2bd698 100644 --- a/src/map.jl +++ b/src/map.jl @@ -65,8 +65,11 @@ with zeroes: Finally, the following examples show how to use `SharedArray`: using SharedArrays + + # Create a map with all pixels set to zero mymap = Healpix.Map{Float64, Healpix.NestedOrder, SharedArray{Float64, 1}}(64) + # Create a map and initialize pixel values with a SharedArray pixels = SharedArray{Int64, 1}(1:12 |> collect) mymap = Healpix.Map{Int64, Healpix.RingOrder, SharedArray{Int64, 1}}(m) """ From 7052d412cc4ce5f89a31c707882ad9e2648090f6 Mon Sep 17 00:00:00 2001 From: Maurizio Tomasi Date: Sat, 9 May 2020 07:26:03 +0200 Subject: [PATCH 4/4] Provide simpler constructors for Alm --- src/alm.jl | 5 +++-- src/sphtfunc.jl | 23 ++++++++++------------- test/test_alm.jl | 8 +++++++- test/test_sphtfunc.jl | 16 ++++++++-------- 4 files changed, 28 insertions(+), 24 deletions(-) diff --git a/src/alm.jl b/src/alm.jl index 3eaa36d..6dd867b 100644 --- a/src/alm.jl +++ b/src/alm.jl @@ -45,8 +45,9 @@ mutable struct Alm{T <: Number, AA <: AbstractArray{T, 1}} end Alm{T}(lmax, mmax) where {T <: Number} = Alm{T, Array{T, 1}}(lmax, mmax) -Alm{T}(lmax, mmax, arr) where {T <: Number} = - Alm{T, Array{T, 1}}(lmax, mmax, collect(arr)) +Alm(lmax, mmax) = Alm{ComplexF64}(lmax, mmax) +Alm(lmax, mmax, arr::AA) where {T <: Number, AA <: AbstractArray{T,1}} = + Alm{T, AA}(lmax, mmax, arr) ################################################################################ diff --git a/src/sphtfunc.jl b/src/sphtfunc.jl index ce372a6..ac1da4f 100644 --- a/src/sphtfunc.jl +++ b/src/sphtfunc.jl @@ -146,7 +146,7 @@ function map2alm(map::Map{Float64, RingOrder, Array{Float64, 1}}; lmax = isnothing(lmax) ? 3 * nside - 1 : lmax mmax = isnothing(mmax) ? lmax : mmax nalms = numberOfAlms(lmax, mmax) - alm = Alm{ComplexF64}(lmax, mmax, zeros(ComplexF64, nalms)) + alm = Alm(lmax, mmax, zeros(ComplexF64, nalms)) map2alm!(map, alm; niter=niter) return alm @@ -159,9 +159,9 @@ function map2alm(map::PolarizedMap{Float64, RingOrder, Array{Float64, 1}}; lmax = isnothing(lmax) ? 3 * nside - 1 : lmax mmax = isnothing(mmax) ? lmax : mmax nalms = numberOfAlms(lmax, mmax) - alms = [Alm{ComplexF64}(lmax, mmax, zeros(ComplexF64, nalms)), - Alm{ComplexF64}(lmax, mmax, zeros(ComplexF64, nalms)), - Alm{ComplexF64}(lmax, mmax, zeros(ComplexF64, nalms))] + alms = [Alm(lmax, mmax, zeros(ComplexF64, nalms)), + Alm(lmax, mmax, zeros(ComplexF64, nalms)), + Alm(lmax, mmax, zeros(ComplexF64, nalms))] map2alm!(map, alms; niter=niter) return alms @@ -256,9 +256,9 @@ are converted to types based on Float64. # Arguments - `alm`: the spherical harmonic coefficients to transform. If of type - `Alm{T}`, we assume a spin-0 spherical harmonic transform. If an - array of `Alm` is passed, we assume that the components correspond - to T, E, and B coefficients. + `Alm{T, Array{T, 1}}`, we assume a spin-0 spherical harmonic + transform. If an array of `Alm` is passed, we assume that the + components correspond to T, E, and B coefficients. # Keywords - `nside::Integer`: Healpix resolution parameter @@ -303,11 +303,8 @@ function alm2map( ) where T <: Real lmax = alms[1].lmax mmax = alms[1].mmax - alm_t = Alm{ComplexF64}(lmax, mmax, - convert(Array{ComplexF64,1}, alms[1].alm)) - alm_e = Alm{ComplexF64}(lmax, mmax, - convert(Array{ComplexF64,1}, alms[2].alm)) - alm_b = Alm{ComplexF64}(lmax, mmax, - convert(Array{ComplexF64,1}, alms[3].alm)) + alm_t = Alm(lmax, mmax, convert(Array{ComplexF64,1}, alms[1].alm)) + alm_e = Alm(lmax, mmax, convert(Array{ComplexF64,1}, alms[2].alm)) + alm_b = Alm(lmax, mmax, convert(Array{ComplexF64,1}, alms[3].alm)) return alm2map([alm_t, alm_e, alm_b], nside) end diff --git a/test/test_alm.jl b/test/test_alm.jl index 33f0344..ceb3249 100644 --- a/test/test_alm.jl +++ b/test/test_alm.jl @@ -6,6 +6,12 @@ @test_throws DomainError(-1, "`mmax` is not positive or zero") Healpix.numberOfAlms(4, -1) @test_throws DomainError((5, 7), "`lmax` and `mmax` are inconsistent") Healpix.numberOfAlms(5, 7) +alm = Healpix.Alm(10, 8) +@test Healpix.almIndex(alm, 4, 2) == 24 +@test Healpix.almIndex(alm, 5, 2) == 25 +@test Healpix.almIndex(alm, 5, 3) == 33 +@test Healpix.almIndex(alm, [4, 6, 5], [3, 4, 5]) == [32, 41, 46] + alm = Healpix.Alm{ComplexF32}(10, 8) @test Healpix.almIndex(alm, 4, 2) == 24 @test Healpix.almIndex(alm, 5, 2) == 25 @@ -43,5 +49,5 @@ alm = Healpix.readAlmFromFITS("alm.fits", ComplexF64) @test alm[28] ≈ (-6.698490836781e-01 + 4.661675665246e+00im) atol = eps ## test alm2cl -testalm = Healpix.Alm{Complex{Float64}}(2,2,ComplexF64.(1:6)) +testalm = Healpix.Alm(2,2,ComplexF64.(1:6)) @test isapprox(Healpix.alm2cl(testalm), [1., 12., 26.2]) diff --git a/test/test_sphtfunc.jl b/test/test_sphtfunc.jl index 9a6d88a..28975de 100644 --- a/test/test_sphtfunc.jl +++ b/test/test_sphtfunc.jl @@ -43,7 +43,7 @@ reference_Δ = [ nside = 2 lmax = 2 nalms = Healpix.numberOfAlms(lmax, lmax) -alm = Healpix.Alm{ComplexF64}(lmax, lmax, 2 .* ones(ComplexF64, nalms)) +alm = Healpix.Alm(lmax, lmax, 2 .* ones(ComplexF64, nalms)) map = Healpix.alm2map(alm, nside) test_map_spin0 = [ @@ -60,7 +60,7 @@ test_map_spin0 = [ @test isapprox(map.pixels, test_map_spin0) ## test type convert -alm_bf = Healpix.Alm{BigFloat}(lmax, lmax, 2 .* ones(BigFloat, nalms)) +alm_bf = Healpix.Alm(lmax, lmax, 2 .* ones(BigFloat, nalms)) map_bf = Healpix.alm2map(alm, nside) @test isapprox(map_bf.pixels, test_map_spin0) @@ -113,9 +113,9 @@ reference_Δ = [ nside = 2 lmax = 2 nalms = Healpix.numberOfAlms(lmax, lmax) -alm_t = Healpix.Alm{ComplexF64}(lmax, lmax, 2 .* ones(ComplexF64, nalms)) -alm_e = Healpix.Alm{ComplexF64}(lmax, lmax, 2 .* ones(ComplexF64, nalms)) -alm_b = Healpix.Alm{ComplexF64}(lmax, lmax, 2 .* ones(ComplexF64, nalms)) +alm_t = Healpix.Alm(lmax, lmax, 2 .* ones(ComplexF64, nalms)) +alm_e = Healpix.Alm(lmax, lmax, 2 .* ones(ComplexF64, nalms)) +alm_b = Healpix.Alm(lmax, lmax, 2 .* ones(ComplexF64, nalms)) maps = Healpix.alm2map([alm_t, alm_e, alm_b], nside) @@ -148,9 +148,9 @@ test_map_spin2_u = [ @test isapprox(maps.u , test_map_spin2_u) ## test type conversion -alm_t = Healpix.Alm{Complex{Float16}}(lmax, lmax, 2 .* ones(Complex{Float16}, nalms)) -alm_e = Healpix.Alm{Complex{Float16}}(lmax, lmax, 2 .* ones(Complex{Float16}, nalms)) -alm_b = Healpix.Alm{Complex{Float16}}(lmax, lmax, 2 .* ones(Complex{Float16}, nalms)) +alm_t = Healpix.Alm(lmax, lmax, 2 .* ones(Complex{Float16}, nalms)) +alm_e = Healpix.Alm(lmax, lmax, 2 .* ones(Complex{Float16}, nalms)) +alm_b = Healpix.Alm(lmax, lmax, 2 .* ones(Complex{Float16}, nalms)) maps_float = Healpix.alm2map([alm_t, alm_e, alm_b], nside) @test isapprox(maps_float.i , test_map_spin0)