Generalize the definition of Map and Alm

src/alm.jl

@@ -1,24 +1,54 @@
 # 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(lmax, mmax) = Alm{ComplexF64}(lmax, mmax)
+Alm(lmax, mmax, arr::AA) where {T <: Number, AA <: AbstractArray{T,1}} =
+    Alm{T, AA}(lmax, mmax, arr)
+
 ################################################################################
 
 """

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

src/map.jl

@@ -27,70 +27,118 @@ 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
+
+    # 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)
 """
-mutable struct Map{T, O <: Order} <: GenericMap{T}
-    pixels::Array{T,1}
+mutable struct Map{T, O <: Order, AA <: AbstractArray{T, 1}} <: GenericMap{T}
+    pixels::AA
     resolution::Resolution
 
     """
         Map{T, O <: Order}(nside) -> Map{T, O}
 
     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

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,16 +54,16 @@ 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)
 
     interpolate(m, θ, ϕ, pixbuf, weightbuf)
 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.