From 969295ab0cd90f2ff39161eb92e378d926ce9f37 Mon Sep 17 00:00:00 2001 From: xzackli Date: Thu, 9 Nov 2023 16:25:14 -0800 Subject: [PATCH 1/4] profile szpack table reading --- src/profiles_szp.jl | 22 +++++++--------------- 1 file changed, 7 insertions(+), 15 deletions(-) diff --git a/src/profiles_szp.jl b/src/profiles_szp.jl index a3bdd42..a535dd6 100644 --- a/src/profiles_szp.jl +++ b/src/profiles_szp.jl @@ -10,21 +10,13 @@ using DelimitedFiles using Interpolations -# write a function that makes szpack_interp (i.e. set_up_szpack_table(file_name)) and pass it to SZpack -#function setup_szpack_table(file_name) -# table = readdlm("/home/lkuhn/CITA-2023/Notes/szpack_interp.dat") -# nu_vector = LinRange(log(35.6888844460172*1e9),log(5353.33266690298*1e9),3000) -# temp_vector = LinRange(1.0e-3,30.0,100) -# szpack_interp = scale(Interpolations.interpolate(table, BSpline(Cubic(Line(OnGrid())))), (temp_vector), (nu_vector)) -# -# return szpack_interp - - -table = readdlm("/home/lkuhn/CITA-2023/Notes/szpack_interp.dat") -nu_vector = LinRange(log(35.6888844460172*1e9),log(5353.33266690298*1e9),3000) -temp_vector = LinRange(1.0e-3,30.0,100) -szpack_interp = scale(Interpolations.interpolate(table, BSpline(Cubic(Line(OnGrid())))), (temp_vector), (nu_vector)) - +function read_szpack_table(filename) + table = readdlm(filename) + nu_vector = LinRange(log(35.6888844460172*1e9),log(5353.33266690298*1e9),3000) + temp_vector = LinRange(1.0e-3,30.0,100) + szpack_interp = scale(Interpolations.interpolate(table, BSpline(Cubic(Line(OnGrid())))), (temp_vector), (nu_vector)) + return szpack_interp +end function SZpack(𝕡, M_200, z, r, τ=0.01) """ From d43bdd8849c87702b17ab40dc6d0df992a8fe232 Mon Sep 17 00:00:00 2001 From: xzackli Date: Thu, 9 Nov 2023 16:40:05 -0800 Subject: [PATCH 2/4] collect deps in main file --- src/XGPaint.jl | 6 +++++- src/profiles.jl | 9 --------- src/profiles_rsz.jl | 10 ---------- src/profiles_szp.jl | 11 ----------- 4 files changed, 5 insertions(+), 31 deletions(-) diff --git a/src/XGPaint.jl b/src/XGPaint.jl index d45ceeb..cad1ef8 100644 --- a/src/XGPaint.jl +++ b/src/XGPaint.jl @@ -18,9 +18,13 @@ using JLD2, FileIO import ThreadsX import Distributions import DataInterpolations # used for uneven spaced interpolators - using DelimitedFiles +import PhysicalConstants.CODATA2018 as constants +const M_sun = 1.98847e30u"kg" +const T_cmb = 2.725 * u"K" +const P_e_factor = constants.σ_e / (constants.m_e * constants.c_0^2) + include("./util.jl") include("./model.jl") include("./profiles.jl") diff --git a/src/profiles.jl b/src/profiles.jl index 17137a3..0508a33 100644 --- a/src/profiles.jl +++ b/src/profiles.jl @@ -1,13 +1,4 @@ -import PhysicalConstants.CODATA2018 as constants -using Unitful -const M_sun = 1.98847e30u"kg" -const T_cmb = 2.725 * u"K" -const P_e_factor = constants.σ_e / (constants.m_e * constants.c_0^2) -using Cosmology -using QuadGK - - # RECTANGULAR WORKSPACES diff --git a/src/profiles_rsz.jl b/src/profiles_rsz.jl index c9a5898..79caa77 100644 --- a/src/profiles_rsz.jl +++ b/src/profiles_rsz.jl @@ -1,13 +1,4 @@ -import PhysicalConstants.CODATA2018 as constants -using Unitful -const M_sun = 1.98847e30u"kg" -const P_e_factor = constants.σ_e / (constants.m_e * constants.c_0^2) -const T_cmb = 2.725 * u"K" -using Cosmology -using QuadGK - - struct Battaglia16RelativisticSZProfile{T,C} <: AbstractGNFW{T} f_b::T # Omega_b / Omega_c = 0.0486 / 0.2589 cosmo::C @@ -89,7 +80,6 @@ function rSZ(𝕡, M_200, z, r) I = (X^3/(ℯ^X-1)) * (2*(2π)^4*(constants.k_B*T_cmb)^3)/((constants.h*constants.c_0)^2) * n T = I/abs((2 * constants.h^2 * ω^4 * ℯ^X)/(constants.k_B * constants.c_0^2 * T_cmb * (ℯ^X - 1)^2)) - #return T return abs(T) end diff --git a/src/profiles_szp.jl b/src/profiles_szp.jl index a535dd6..9e4004a 100644 --- a/src/profiles_szp.jl +++ b/src/profiles_szp.jl @@ -1,15 +1,4 @@ -import PhysicalConstants.CODATA2018 as constants -using Unitful -const M_sun = 1.98847e30u"kg" -const P_e_factor = constants.σ_e / (constants.m_e * constants.c_0^2) -const T_cmb = 2.725 * u"K" -using Cosmology -using QuadGK -using DelimitedFiles -using Interpolations - - function read_szpack_table(filename) table = readdlm(filename) nu_vector = LinRange(log(35.6888844460172*1e9),log(5353.33266690298*1e9),3000) From 8de243c2d159b48865987c64a8977756271bdfe3 Mon Sep 17 00:00:00 2001 From: Martine Lokken Date: Fri, 24 Mar 2023 14:02:57 -0400 Subject: [PATCH 3/4] Incorporated break model; added healpix example attempt --- examples/make_map_raw_car_martine.jl | 204 +++++++++++++++++++++++++++ src/profiles.jl | 3 - 2 files changed, 204 insertions(+), 3 deletions(-) create mode 100644 examples/make_map_raw_car_martine.jl diff --git a/examples/make_map_raw_car_martine.jl b/examples/make_map_raw_car_martine.jl new file mode 100644 index 0000000..d97b176 --- /dev/null +++ b/examples/make_map_raw_car_martine.jl @@ -0,0 +1,204 @@ +# modified from Zack Li + +using Pixell, WCS, XGPaint +using Cosmology +using Interpolations +import XGPaint: AbstractProfile +using HDF5 +import JSON +using JLD2 +using ThreadsX +using FileIO +using Healpix + +print("Threads: ", Threads.nthreads(), "\n") +modeltype::String = "break" +healpix = true + + +if healpix + nside = 1024 + m = HealpixMap{Float64,RingOrder}(nside) +else + shape, wcs = fullsky_geometry(0.5 * Pixell.arcminute) + # for a box: + # box = [30 -30; # RA + # -25 25] * Pixell.degree # DEC + # shape, wcs = geometry(CarClenshawCurtis{Float64}, box, 0.5 * Pixell.arcminute) + m = Enmap(zeros(shape), wcs) + # precomputed sky angles + α_map, δ_map = posmap(shape, wcs) + psa = (sin_α=sin.(α_map), cos_α=cos.(α_map), sin_δ=sin.(δ_map), cos_δ=cos.(δ_map)) +end + +## +fid = h5open("/mnt/raid-cita/mlokken/buzzard/catalogs/halos/buzzard_halos.hdf5", "r") +# fid = h5open("/mnt/scratch-lustre/mlokken/pkpatch/halos_hdf5.h5", "r") +ra, dec = deg2rad.(fid["ra"]), deg2rad.(fid["dec"]) +redshift = collect(fid["z"]) +halo_mass = collect(fid["m200c"]) +# choose the cutoff for the pressure profiles here +cutoff = 4 # the standard +apply_beam = true +## +perm = sortperm(dec, alg=ThreadsX.MergeSort) +ra = ra[perm] +dec = dec[perm] +redshift = redshift[perm] +halo_mass = halo_mass[perm] + + + +## +print("Precomputing the model profile grid.\n") + +# set up a profile to paint +if modeltype=="break" + p = XGPaint.BreakModel(Omega_c=0.24, Omega_b=0.046, h=0.7, alpha_break=1.5) # Buzzard cosmology +elseif modeltype=="battaglia" + p = XGPaint.BattagliaProfile(Omega_c=0.24, Omega_b=0.046, h=0.7) # Buzzard cosmology +end + +# beam stuff +N_logθ = 512 +rft = RadialFourierTransform(n=N_logθ, pad=256) +beamstr = "" + +if apply_beam + ells = rft.l + fwhm = deg2rad(1.6/60.) # ACT beam + σ² = fwhm^2 / (8log(2)) + lbeam = @. exp(-ells * (ells+1) * σ² / 2) + beamstr = "_1p6arcmin" +end +model_file::String = "cached_$(modeltype)_$(beamstr).jld2" +if isfile(model_file) + print("Found cached Battaglia profile model. Loading from disk.\n") + model = load(model_file) + prof_logθs, prof_redshift, prof_logMs, prof_y = model["prof_logθs"], + model["prof_redshift"], model["prof_logMs"], model["prof_y"] +else + print("Didn't find a cached profile model. Computing and saving.\n") + logθ_min, logθ_max = log(minimum(rft.r)), log(maximum(rft.r)) + @time prof_logθs, prof_redshift, prof_logMs, prof_y = profile_grid(p; + N_logθ=length(rft.r), logθ_min=logθ_min, logθ_max=logθ_max) + + if apply_beam + # now apply the beam + @time XGPaint.transform_profile_grid!(prof_y, rft, lbeam) + prof_y = prof_y[begin+rft.pad:end-rft.pad, :, :] + prof_logθs = prof_logθs[begin+rft.pad:end-rft.pad] + XGPaint.cleanup_negatives!(prof_y) + end + jldsave(model_file; prof_logθs, prof_redshift, prof_logMs, prof_y) + +end + +if ~healpix # do this for CAR + itp = Interpolations.interpolate(log.(prof_y), BSpline(Cubic(Line(OnGrid())))) + sitp = scale(itp, prof_logθs, prof_redshift, prof_logMs); +end + + +## +function paint_map!(m::Enmap, p::XGPaint.AbstractProfile, psa, sitp, masses, + redshifts, αs, δs, irange; mult=4) + for i in irange + α₀ = αs[i] + δ₀ = δs[i] + mh = masses[i] + z = redshifts[i] + θmax = XGPaint.θmax(p, mh * XGPaint.M_sun, z, mult=mult) + profile_paint!(m, α₀, δ₀, psa, sitp, z, mh, θmax) + end +end + +function paint_map!(m::HealpixMap, p::XGPaint.AbstractProfile{T}, masses, + redshifts, αs, δs, irange; mult=4) where {T} + for i in irange + α₀ = αs[i] + δ₀ = δs[i] + mh = masses[i] + z = redshifts[i] + θmax = XGPaint.θmax(p, mh * XGPaint.M_sun, z, mult=mult) + θmin = 0.0 # don't know what this should be # error + beam_interp = XGPaint.realspacegaussbeam(T,fwhm)[0] #error + w = XGPaint.HealpixPaintingWorkspace(nside, θmin, θmax, beam_interp) + profile_paint!(m, α₀, δ₀, w, z, mh, θmax) + end +end + +function chunked_paint!(m::Enmap, p::XGPaint.AbstractProfile, psa, sitp, masses, + redshifts, αs, δs; mult=4) + m .= 0.0 + + N_sources = length(masses) + chunksize = ceil(Int, N_sources / (2Threads.nthreads())) + chunks = chunk(N_sources, chunksize); + + Threads.@threads for i in 1:Threads.nthreads() + chunk_i = 2i + i1, i2 = chunks[chunk_i] + paint_map!(m, p, psa, sitp, masses, redshifts, αs, δs, i1:i2, mult=mult) + end + + Threads.@threads for i in 1:Threads.nthreads() + chunk_i = 2i - 1 + i1, i2 = chunks[chunk_i] + paint_map!(m, p, psa, sitp, masses, redshifts, αs, δs, i1:i2, mult=mult) + end +end + +function chunked_paint!(m::HealpixMap, p::XGPaint.AbstractProfile, masses, + redshifts, αs, δs; mult=4) + m .= 0.0 + + N_sources = length(masses) + chunksize = ceil(Int, N_sources / (2Threads.nthreads())) + chunks = chunk(N_sources, chunksize); + + Threads.@threads for i in 1:Threads.nthreads() + chunk_i = 2i + i1, i2 = chunks[chunk_i] + paint_map!(m, p, masses, redshifts, αs, δs, i1:i2, mult=mult) + end + + Threads.@threads for i in 1:Threads.nthreads() + chunk_i = 2i - 1 + i1, i2 = chunks[chunk_i] + paint_map!(m, p, masses, redshifts, αs, δs, i1:i2, mult=mult) + end +end +## +# cut = eachindex(halo_mass)[begin:end-5] + +print(maximum(ra), minimum(ra), maximum(dec), minimum(dec)) +print("Painting map.\n") +if healpix + @time chunked_paint!(m, p, halo_mass, redshift, ra, dec, mult=cutoff) # for a Healpix +else + @time chunked_paint!(m, p, psa, sitp, halo_mass, redshift, ra, dec, mult=cutoff) # for an Enmap +end + +# +write_map( + "/mnt/raid-cita/mlokken/buzzard/ymaps/ymap_buzzard_$(modeltype)_bbps_car$(beamstr)_cutoff$(cutoff).fits", + m) + +using PyPlot +plt.clf() +plt.figure() +plt.imshow(log10.(m.data')) +plt.axis("off") +plt.savefig("test_Healpix.png", bbox_inches="tight",pad_inches = 0) +plt.gcf() + + + +# ## +# m .= 0 +# profile_paint!(m, 0., 0., psa, sitp, 0.1, 1e14, π/90) + +# ## +# using Plots +# Plots.plot(log10.(m)) diff --git a/src/profiles.jl b/src/profiles.jl index 0508a33..acfe834 100644 --- a/src/profiles.jl +++ b/src/profiles.jl @@ -54,9 +54,6 @@ function Battaglia16ThermalSZProfile(; Omega_c::T=0.2589, Omega_b::T=0.0486, h:: return Battaglia16ThermalSZProfile(f_b, cosmo) end -abstract type AbstractPaintingProblem{T} end - - function BreakModel(; Omega_c::T=0.2589, Omega_b::T=0.0486, h::T=0.6774, alpha_break::T=1.5, M_break::T=2.0*10^14) where {T <: Real} #alpha_break = 1.486 from Shivam P paper by Nate's sleuthing OmegaM=Omega_b+Omega_c From 9423fd456072dba525c0d09e4a32ac2af9a0067a Mon Sep 17 00:00:00 2001 From: xzackli Date: Wed, 15 Nov 2023 09:05:27 -0800 Subject: [PATCH 4/4] working enmap painter --- src/XGPaint.jl | 1 + src/profiles_szp.jl | 95 +++++++++++++++++++++------------------------ 2 files changed, 46 insertions(+), 50 deletions(-) diff --git a/src/XGPaint.jl b/src/XGPaint.jl index cad1ef8..78fd12f 100644 --- a/src/XGPaint.jl +++ b/src/XGPaint.jl @@ -39,5 +39,6 @@ export get_cosmology, read_halo_catalog_hdf5, sort_halo_catalog export Radio_Sehgal2009, CIB_Planck2013, CIB_Scarfy, CO_CROWNED, LRG_Yuan23 export paint!, generate_sources, process_sources, profile_grid, profile_paint!, profileworkspace, paint_szp!, profile_grid_szp, profile_paint_szp!, paint_rsz!, profile_grid_rsz, profile_paint_rsz! export build_interpolator, Battaglia16ThermalSZProfile, Battaglia16RelativisticSZProfile, build_interpolator_szp, build_interpolator_rsz +export Battaglia16SZPackProfile end # module diff --git a/src/profiles_szp.jl b/src/profiles_szp.jl index 9e4004a..5f186ed 100644 --- a/src/profiles_szp.jl +++ b/src/profiles_szp.jl @@ -7,6 +7,26 @@ function read_szpack_table(filename) return szpack_interp end +struct Battaglia16SZPackProfile{T,C,TSZ, ITP1, ITP2} <: AbstractGNFW{T} + f_b::T # Omega_b / Omega_c = 0.0486 / 0.2589 + cosmo::C + X::T # X = 2.6408 corresponding to frequency 150 GHz + 𝕡_tsz::TSZ + tsz_interp::ITP1 + szpack_interp::ITP2 + τ::T +end + +function Battaglia16SZPackProfile(𝕡_tsz, tsz_interp, filename::String, τ; + Omega_c::T=0.2589, Omega_b::T=0.0486, h::T=0.6774, x::T=2.6408) where {T <: Real} + OmegaM=Omega_b+Omega_c + f_b = Omega_b / OmegaM + cosmo = get_cosmology(T, h=h, OmegaM=OmegaM) + X = x + szpack_interp = read_szpack_table(filename) + return Battaglia16SZPackProfile(f_b, cosmo, X, 𝕡_tsz, tsz_interp, szpack_interp, τ) +end + function SZpack(𝕡, M_200, z, r, τ=0.01) """ Outputs the integrated compton-y signal calculated using SZpack along the line of sight. @@ -15,18 +35,11 @@ function SZpack(𝕡, M_200, z, r, τ=0.01) T_e = T_vir_calc(𝕡, M_200, z) θ_e = (constants.k_B*T_e)/(constants.m_e*constants.c_0^2) ω = (X*constants.k_B*T_cmb)/constants.ħ - + t = ustrip(uconvert(u"keV",T_e * constants.k_B)) nu = log(ustrip(uconvert(u"Hz",ω))) - #if t > 30 - ## t = 30 - # println(t) - # println(M_200) - # println(z) - # println(X) - #end - dI = szpack_interp(t, nu)*u"MJy/sr" - + dI = 𝕡.szpack_interp(t, nu)*u"MJy/sr" + y = XGPaint.compton_y_rsz(𝕡, M_200, z, r) I = y * (dI/(τ * θ_e)) * (2π)^4 T = I/abs((2 * constants.h^2 * ω^4 * ℯ^X)/(constants.k_B * constants.c_0^2 * T_cmb * (ℯ^X - 1)^2)) @@ -68,9 +81,10 @@ end -function profile_paint_szp!(m::Enmap{T, 2, Matrix{T}, CarClenshawCurtis{T}}, p, +function profile_paint_szp!(m::Enmap{T, 2, Matrix{T}, CarClenshawCurtis{T}}, + p::Battaglia16SZPackProfile, α₀, δ₀, psa::CarClenshawCurtisProfileWorkspace, - sitp, z, Ms, θmax) where T + z, Ms, θmax) where T # get indices of the region to work on i1, j1 = sky2pix(m, α₀ - θmax, δ₀ - θmax) i2, j2 = sky2pix(m, α₀ + θmax, δ₀ + θmax) @@ -78,14 +92,24 @@ function profile_paint_szp!(m::Enmap{T, 2, Matrix{T}, CarClenshawCurtis{T}}, p, i_stop = ceil(Int, min(max(i1, i2), size(m, 1))) j_start = floor(Int, max(min(j1, j2), 1)) j_stop = ceil(Int, min(max(j1, j2), size(m, 2))) - - X_0 = calc_null(p, Ms*M_sun, z) + + # needs mass in M_200 X = p.X - if X > X_0 - sign = 1 - else - sign = -1 - end + T_e = T_vir_calc(p, Ms * M_sun, z) + θ_e = (constants.k_B*T_e)/(constants.m_e*constants.c_0^2) + ω = (X*constants.k_B*T_cmb)/constants.ħ + t = ustrip(uconvert(u"keV",T_e * constants.k_B)) + nu = log(ustrip(uconvert(u"Hz",ω))) + + logMs = log10(Ms) + dI = p.szpack_interp(t, nu)*u"MJy/sr" + rsz_factor_I = (dI/(p.τ * θ_e)) * (2π)^4 + rsz_factor_T = abs(rsz_factor_I / ( (2 * constants.h^2 * ω^4 * ℯ^X) / + (constants.k_B * constants.c_0^2 * T_cmb * (ℯ^X - 1)^2))) + X_0 = calc_null(p, Ms*M_sun, z) + if X < X_0 + rsz_factor_T *= -1 + end x₀ = cos(δ₀) * cos(α₀) y₀ = cos(δ₀) * sin(α₀) @@ -98,42 +122,13 @@ function profile_paint_szp!(m::Enmap{T, 2, Matrix{T}, CarClenshawCurtis{T}}, p, z₁ = psa.sin_δ[i,j] d² = (x₁ - x₀)^2 + (y₁ - y₀)^2 + (z₁ - z₀)^2 θ = acos(1 - d² / 2) - m[i,j] += ifelse(θ < θmax, - sign * exp(sitp(log(θ), z, log10(Ms))), - zero(T)) + y = exp(p.tsz_interp(log(θ), z, logMs)) + m[i,j] += (θ < θmax) * rsz_factor_T * y end end end -"""Helper function to build a tSZ interpolator""" -function build_interpolator_szp(model::AbstractGNFW; cache_file::String="", - N_logθ=512, pad=256, overwrite=true, verbose=true) - - if overwrite || (isfile(cache_file) == false) - verbose && print("Building new interpolator from model.\n") - rft = RadialFourierTransform(n=N_logθ, pad=pad) - logθ_min, logθ_max = log(minimum(rft.r)), log(maximum(rft.r)) - prof_logθs, prof_redshift, prof_logMs, prof_y = profile_grid_szp(model; - N_logθ=N_logθ, logθ_min=logθ_min, logθ_max=logθ_max) - if length(cache_file) > 0 - verbose && print("Saving new interpolator to $(cache_file).\n") - save(cache_file, Dict("prof_logθs"=>prof_logθs, - "prof_redshift"=>prof_redshift, "prof_logMs"=>prof_logMs, "prof_y"=>prof_y)) - end - else - print("Found cached Battaglia profile model. Loading from disk.\n") - model_grid = load(cache_file) - prof_logθs, prof_redshift, prof_logMs, prof_y = model_grid["prof_logθs"], - model_grid["prof_redshift"], model_grid["prof_logMs"], model_grid["prof_y"] - end - - itp = Interpolations.interpolate(log.(prof_y), BSpline(Cubic(Line(OnGrid())))) - sitp = scale(itp, prof_logθs, prof_redshift, prof_logMs) - return sitp -end - - function profile_paint_szp!(m::HealpixMap{T, RingOrder}, p, α₀, δ₀, w::HealpixProfileWorkspace, z, Mh, θmax) where T ϕ₀ = α₀