healpix workspace renamed

Author: xzackli

docs/src/ksz.md

@@ -69,11 +69,12 @@ To generate Healpix maps, you'll need the [Healpix.jl](https://github.com/ziotom
 
 ```@example ksz
 using Healpix
+using XGPaint
 
 nside = 2048
 m_hp = HealpixMap{Float64,RingOrder}(nside)
-max_radius = deg2rad(5.0)  # maximum radius to consider for the profile
-w = HealpixProfileWorkspace(nside, max_radius)
+res = Healpix.Resolution(nside)
+w = HealpixRingProfileWorkspace{Float64}(res)
 
 @time paint!(m_hp, w, model_interp, halo_mass, redshift, ra, dec, proj_v_over_c)
 # Healpix.saveToFITS(m_hp, "!y.fits", typechar="D")  # to save, uncomment

docs/src/sz.md

@@ -79,11 +79,12 @@ To generate Healpix maps, you'll need the [Healpix.jl](https://github.com/ziotom
 
 ```@example tsz
 using Healpix
+using XGPaint
 
 nside = 4096
 m_hp = HealpixMap{Float64,RingOrder}(nside)
-max_radius = deg2rad(5.0)  # maximum radius to consider for the profile
-w = HealpixProfileWorkspace(nside, max_radius)
+res = Healpix.Resolution(nside)
+w = HealpixRingProfileWorkspace{Float64}(res)
 @time paint!(m_hp, w, y_model_interp, halo_mass, redshift, ra, dec)
 # Healpix.saveToFITS(m_hp, "!y.fits", typechar="D")  # to save, uncomment
 

src/XGPaint.jl

@@ -51,6 +51,6 @@ export Radio_Sehgal2009, CIB_Planck2013, CIB_Scarfy, CO_CROWNED, LRG_Yuan23
 export paint!, generate_sources, process_sources, profile_grid, profile_paint!
 export profileworkspace, paint_szp!, profile_grid_szp, profile_paint_szp!, paint_rsz!, profile_grid_rsz, profile_paint_rsz!
 export build_interpolator, Battaglia16ThermalSZProfile, RSZPerturbativeProfile, build_interpolator_szp, build_interpolator_rsz
-export SZPackRSZProfile, nu_to_X, X_to_nu, BattagliaTauProfile, HealpixProfileWorkspace
+export SZPackRSZProfile, nu_to_X, X_to_nu, BattagliaTauProfile, HealpixRingProfileWorkspace
 
 end # module

src/healpixworkspace.jl

@@ -1,14 +1,14 @@
 
-# Ring workspace for fast disc queries
-struct RingWorkspace{T<:AbstractFloat}
+# Ring workspace for fast disc queries; this is designed for CPU for now
+struct HealpixRingProfileWorkspace{T<:AbstractFloat}
     res::Healpix.Resolution
     ring_thetas::Vector{T}     # Colatitude of each ring
     ring_shifts::Vector{T}     # Shift offset for each ring (0.0 or 0.5)
     ring_num_pixels::Vector{Int} # Number of pixels in each ring
     ring_first_pixels::Vector{Int} # First pixel index for each ring
 end
 
-function RingWorkspace{T}(res::Healpix.Resolution) where T<:AbstractFloat
+function HealpixRingProfileWorkspace{T}(res::Healpix.Resolution) where T<:AbstractFloat
     num_rings = res.nsideTimesFour - 1
 
     ring_thetas = Vector{T}(undef, num_rings)
@@ -24,10 +24,10 @@ function RingWorkspace{T}(res::Healpix.Resolution) where T<:AbstractFloat
         ring_first_pixels[ring_idx] = ringinfo.firstPixIdx
     end
 
-    RingWorkspace{T}(res, ring_thetas, ring_shifts, ring_num_pixels, ring_first_pixels)
+    HealpixRingProfileWorkspace{T}(res, ring_thetas, ring_shifts, ring_num_pixels, ring_first_pixels)
 end
 
-RingWorkspace(res::Healpix.Resolution) = RingWorkspace{Float64}(res)
+HealpixRingProfileWorkspace(res::Healpix.Resolution) = HealpixRingProfileWorkspace{Float64}(res)
 
 function get_relevant_rings(res::Healpix.Resolution, θ_center, radius)
     z_top = cos(max(0, θ_center - radius))
@@ -53,8 +53,11 @@ Returns two ranges to handle φ wraparound at the 0/2π boundary:
 - range1: Main contiguous range of pixels  
 - range2: Additional range when disc crosses φ=0/2π seam (empty if no crossing)
 """
-function get_ring_disc_ranges(workspace::RingWorkspace, ring_idx::Int, 
+function get_ring_disc_ranges(workspace::HealpixRingProfileWorkspace{T}, ring_idx::Int, 
                               center_theta::T, center_phi::T, radius::T) where T
+    # Validate that phi is in the valid range [0, 2π)
+    @assert 0 <= center_phi < T(2π) "center_phi must be in the range [0, 2π), got $(center_phi)"
+    
     nr = workspace.ring_num_pixels[ring_idx]
     ring_theta = workspace.ring_thetas[ring_idx]
     shift = workspace.ring_shifts[ring_idx]

src/profiles.jl

@@ -350,9 +350,9 @@ end
 #     end
 # end
 
-function profile_paint_generic!(m::HealpixMap{T, RingOrder}, workspace::RingWorkspace{T}, 
+function profile_paint_generic!(m::HealpixMap{T, RingOrder}, workspace::HealpixRingProfileWorkspace{T}, 
         model, Mh, z, α₀, δ₀, θmax, normalization=1) where T
-    ϕ₀ = α₀  
+    ϕ₀ = mod(T(α₀), T(2π))  # Normalize RA to [0, 2π)
     θ₀ = T(π)/2 - δ₀
     x₀, y₀, z₀ = ang2vec(θ₀, ϕ₀)
     θmin = compute_θmin(model)
@@ -391,7 +391,7 @@ function profile_paint_generic!(m::HealpixMap{T, RingOrder}, workspace::RingWork
 end
 
 # fall back to generic profile painter if no specialized painter is defined for the model
-function profile_paint!(m::HealpixMap{T, RingOrder}, w::RingWorkspace{T}, model, 
+function profile_paint!(m::HealpixMap{T, RingOrder}, w::HealpixRingProfileWorkspace{T}, model, 
                         Mh, z, α₀, δ₀, θmax, normalization=1) where T
     profile_paint_generic!(m, w, model, Mh, z, α₀, δ₀, θmax, normalization)
 end

src/profiles_rsz.jl

@@ -133,7 +133,7 @@ end
 
 # like the usual paint, but use the sign of the null as the sign of the perturbation
 function profile_paint!(m::HealpixMap{T, RingOrder}, 
-        workspace::RingWorkspace{T}, model::RSZPerturbativeProfile, 
+        workspace::HealpixRingProfileWorkspace{T}, model::RSZPerturbativeProfile, 
         Mh, z, α₀, δ₀, θmax) where T
     X_0 = calc_null(model, Mh, z)
     profile_paint_generic!(m, workspace, model, Mh, z, α₀, δ₀, θmax, sign(model.X - X_0))

src/profiles_szp.jl

@@ -99,7 +99,7 @@ function profile_paint!(m::Enmap{T, 2, Matrix{T}, CarClenshawCurtis{T}},
 end
 
 function profile_paint!(m::HealpixMap{T, RingOrder}, 
-        workspace::RingWorkspace{T}, model::SZPackRSZProfile, 
+        workspace::HealpixRingProfileWorkspace{T}, model::SZPackRSZProfile, 
         Mh, z, α₀, δ₀, θmax) where T
     rsz_factor_I_over_y = compute_rsz_factor_I_over_y(model, Mh, z)
     profile_paint_generic!(m, workspace, model.y_model_interp, Mh, z, α₀, δ₀, θmax, rsz_factor_I_over_y)

test/runtests.jl

@@ -10,7 +10,7 @@ using PhysicalConstants.CODATA2018
 
 # Include HEALPix ring workspace tests
 include("test_healpix_rings.jl")
-include("test_ringworkspace.jl")
+include("test_HealpixRingProfileWorkspace.jl")
 
 # relative background evolutions differ by 1e-3 between Julia and Python 2
 rtol = 1e-3

test/test_healpix_rings.jl

@@ -4,10 +4,10 @@ using XGPaint
 
 @testset "Ring-Slice HEALPix Disc Query Tests" begin
 
-    @testset "RingWorkspace Construction" begin
+    @testset "HealpixRingProfileWorkspace Construction" begin
         nside = 32
         res = Healpix.Resolution(nside)
-        workspace = XGPaint.RingWorkspace(res)
+        workspace = XGPaint.HealpixRingProfileWorkspace(res)
 
         @test workspace.res === res
         @test length(workspace.ring_thetas) == res.nsideTimesFour - 1
@@ -17,7 +17,7 @@ using XGPaint
         @test eltype(workspace.ring_shifts) == Float64
 
         # Test type parameter
-        workspace_f32 = XGPaint.RingWorkspace{Float32}(res)
+        workspace_f32 = XGPaint.HealpixRingProfileWorkspace{Float32}(res)
         @test eltype(workspace_f32.ring_thetas) == Float32
         @test eltype(workspace_f32.ring_shifts) == Float32
     end
@@ -44,7 +44,7 @@ using XGPaint
     @testset "get_ring_disc_ranges Basic Cases" begin
         nside = 32
         res = Healpix.Resolution(nside)
-        workspace = XGPaint.RingWorkspace(res)
+        workspace = XGPaint.HealpixRingProfileWorkspace(res)
 
         # Test empty disc (radius = 0)
         range1, range2 = XGPaint.get_ring_disc_ranges(workspace, 10, π/4, 0.0, 0.0)
@@ -61,9 +61,9 @@ using XGPaint
     @testset "Systematic Validation vs HEALPix Reference" begin
         nside = 32  # Smaller for faster testing
         res = Healpix.Resolution(nside)
-        workspace = XGPaint.RingWorkspace(res)
+        workspace = XGPaint.HealpixRingProfileWorkspace(res)
 
-        function query_disc_rings!(pixel_list::Vector{Int}, workspace::XGPaint.RingWorkspace, 
+        function query_disc_rings!(pixel_list::Vector{Int}, workspace::XGPaint.HealpixRingProfileWorkspace, 
                                   θ_center, ϕ_center, radius)
             empty!(pixel_list)
             ring_min, ring_max = XGPaint.get_relevant_rings(workspace.res, θ_center, radius)
@@ -105,9 +105,9 @@ using XGPaint
     @testset "Random Points Validation" begin
         nside = 32
         res = Healpix.Resolution(nside)
-        workspace = XGPaint.RingWorkspace(res)
+        workspace = XGPaint.HealpixRingProfileWorkspace(res)
 
-        function query_disc_rings!(pixel_list::Vector{Int}, workspace::XGPaint.RingWorkspace, 
+        function query_disc_rings!(pixel_list::Vector{Int}, workspace::XGPaint.HealpixRingProfileWorkspace, 
                                   θ_center, ϕ_center, radius)
             empty!(pixel_list)
             ring_min, ring_max = XGPaint.get_relevant_rings(workspace.res, θ_center, radius)
@@ -150,7 +150,7 @@ using XGPaint
     @testset "Edge Cases" begin
         nside = 32
         res = Healpix.Resolution(nside)
-        workspace = XGPaint.RingWorkspace(res)
+        workspace = XGPaint.HealpixRingProfileWorkspace(res)
 
         # Test exact north pole
         range1, range2 = XGPaint.get_ring_disc_ranges(workspace, 1, 0.0, 0.0, deg2rad(5.0))

test/test_ringworkspace.jl

@@ -1,5 +1,5 @@
 using XGPaint
-using XGPaint: RingWorkspace
+using XGPaint: HealpixRingProfileWorkspace
 using Healpix
 using Test
 
@@ -42,16 +42,16 @@ end
 (m::TestModel)(θ, Mh, z) = m.value
 XGPaint.compute_θmin(::TestModel{T}) where T = eps(T)
 
-@testset "RingWorkspace Profile Painting Tests" begin
+@testset "HealpixRingProfileWorkspace Profile Painting Tests" begin
 
-    @testset "RingWorkspace vs Brute Force" begin
+    @testset "HealpixRingProfileWorkspace vs Brute Force" begin
         # Setup
         nside = 64
         res = Healpix.Resolution(nside)
-        workspace = RingWorkspace(res)
+        workspace = HealpixRingProfileWorkspace(res)
 
         # Create test maps
-        map_ringworkspace = HealpixMap{Float64, RingOrder}(zeros(Float64, nside2npix(nside)))
+        map_HealpixRingProfileWorkspace = HealpixMap{Float64, RingOrder}(zeros(Float64, nside2npix(nside)))
         map_bruteforce = HealpixMap{Float64, RingOrder}(zeros(Float64, nside2npix(nside)))
 
         # Test model
@@ -69,26 +69,26 @@ XGPaint.compute_θmin(::TestModel{T}) where T = eps(T)
 
         for (i, test_case) in enumerate(test_cases)
             # Reset maps
-            fill!(map_ringworkspace.pixels, 0.0)
+            fill!(map_HealpixRingProfileWorkspace.pixels, 0.0)
             fill!(map_bruteforce.pixels, 0.0)
 
             # Paint with both methods
-            XGPaint.profile_paint_generic!(map_ringworkspace, workspace, model, 
+            XGPaint.profile_paint_generic!(map_HealpixRingProfileWorkspace, workspace, model, 
                                          test_case.Mh, test_case.z, test_case.α₀, test_case.δ₀, test_case.θmax)
 
             profile_paint_bruteforce!(map_bruteforce, model, 
                                     test_case.Mh, test_case.z, test_case.α₀, test_case.δ₀, test_case.θmax)
 
             # Compare results - should be exactly equal for constant model
-            @test map_ringworkspace.pixels ≈ map_bruteforce.pixels
-            @test count(x -> x != 0, map_ringworkspace.pixels) == count(x -> x != 0, map_bruteforce.pixels)
+            @test map_HealpixRingProfileWorkspace.pixels ≈ map_bruteforce.pixels
+            @test count(x -> x != 0, map_HealpixRingProfileWorkspace.pixels) == count(x -> x != 0, map_bruteforce.pixels)
         end
     end
 
-    @testset "RingWorkspace Edge Cases" begin
+    @testset "HealpixRingProfileWorkspace Edge Cases" begin
         nside = 64
         res = Healpix.Resolution(nside)
-        workspace = RingWorkspace(res)
+        workspace = HealpixRingProfileWorkspace(res)
 
         map_test = HealpixMap{Float64, RingOrder}(zeros(Float64, nside2npix(nside)))
         model = TestModel(1.0)
@@ -110,14 +110,14 @@ XGPaint.compute_θmin(::TestModel{T}) where T = eps(T)
         end
     end
 
-    @testset "RingWorkspace vs Brute Force with Real Y Profile" begin
+    @testset "HealpixRingProfileWorkspace vs Brute Force with Real Y Profile" begin
         # Setup with real y profile model
         nside = 32  # Smaller for faster testing with real profile
         res = Healpix.Resolution(nside)
-        workspace = RingWorkspace(res)
+        workspace = HealpixRingProfileWorkspace(res)
 
         # Create test maps
-        map_ringworkspace = HealpixMap{Float64, RingOrder}(zeros(Float64, nside2npix(nside)))
+        map_HealpixRingProfileWorkspace = HealpixMap{Float64, RingOrder}(zeros(Float64, nside2npix(nside)))
         map_bruteforce = HealpixMap{Float64, RingOrder}(zeros(Float64, nside2npix(nside)))
 
         # Real y profile model
@@ -132,26 +132,26 @@ XGPaint.compute_θmin(::TestModel{T}) where T = eps(T)
 
         for (i, test_case) in enumerate(test_cases)
             # Reset maps
-            fill!(map_ringworkspace.pixels, 0.0)
+            fill!(map_HealpixRingProfileWorkspace.pixels, 0.0)
             fill!(map_bruteforce.pixels, 0.0)
 
             # Paint with both methods
-            XGPaint.profile_paint_generic!(map_ringworkspace, workspace, y_model, 
+            XGPaint.profile_paint_generic!(map_HealpixRingProfileWorkspace, workspace, y_model, 
                                          test_case.Mh, test_case.z, test_case.α₀, test_case.δ₀, test_case.θmax)
 
             profile_paint_bruteforce!(map_bruteforce, y_model, 
                                     test_case.Mh, test_case.z, test_case.α₀, test_case.δ₀, test_case.θmax)
 
             # Compare results - should be very close (allowing for floating point differences)
-            @test map_ringworkspace.pixels ≈ map_bruteforce.pixels rtol=1e-12
+            @test map_HealpixRingProfileWorkspace.pixels ≈ map_bruteforce.pixels rtol=1e-12
 
             # Check that we actually painted something non-zero
-            @test count(x -> x != 0, map_ringworkspace.pixels) > 0
+            @test count(x -> x != 0, map_HealpixRingProfileWorkspace.pixels) > 0
             @test count(x -> x != 0, map_bruteforce.pixels) > 0
 
             # Check same number of non-zero pixels
-            @test count(x -> x != 0, map_ringworkspace.pixels) == count(x -> x != 0, map_bruteforce.pixels)
+            @test count(x -> x != 0, map_HealpixRingProfileWorkspace.pixels) == count(x -> x != 0, map_bruteforce.pixels)
         end
     end
 
-end  # RingWorkspace Profile Painting Tests
+end  # HealpixRingProfileWorkspace Profile Painting Tests