Edit on GitHub

Corrections for Point Sources (corrections.jl)

We need to estimate a correction to the covariance matrix due to insufficiently apodized point source holes in the likelihood mask, and we had generated spectra in signalsim. Now we will estimate a smooth correction from these spectra, by calculating the analytic covariance and then using a Gaussian process to smoothly interpolate the deviation from 1 between the sample covariance and the analytic covariance.

$ julia corrections.jl example.toml 143 143 TT

Configuration

As usual, we set up the various information we'll need later.

configfile, freq1, freq2, spec = ARGS

freqs = [freq1, freq2]
splits = ["1", "2"]
print("$(freqs[1]), $(freqs[2]), $(spec)\n")

# We choose to never do ET, we'll just simulate TE but swap the maps.
transposed = (spec == "ET")
if transposed
    spec = "TE"
    splits = reverse(splits)
    freqs = reverse(freqs)
end

using Plots
using Healpix
using PowerSpectra
using TOML
using UUIDs, JLD2, FileIO
using Statistics
using GaussianProcesses

include("util.jl")

config = TOML.parsefile(configfile)
nside = config["general"]["nside"]
run_name = config["general"]["name"]
simpath = joinpath(config["scratch"], "signalsims",
    "$(freq1)_$(freq2)_$(spec)")
lmax = nside2lmax(nside)
lmax_planck = min(2508, lmax)
767

We'll load the existing simulated spectra from disk.

function readsims(path)
    files = [f for f in readdir(simpath; join=false) if f != "summary.jld2"]
    test_spec_ = load(joinpath(path,files[1]), "cl")
    cl_array = zeros((length(test_spec_), length(files)))
    for (i,f) in enumerate(files)
        cl_array[:,i] .= parent(load(joinpath(path,f), "cl"))
    end
    return cl_array
end


@time cl_array = readsims(simpath)
cl_mean = mean(cl_array, dims=2)
cl_var = var(cl_array, dims=2)
@save joinpath(simpath, "summary.jld2") cl_mean cl_var

signal11, th = signal_and_theory(freq1, freq1, config)
signal12, th = signal_and_theory(freq1, freq2, config)
signal21, th = signal_and_theory(freq2, freq1, config)
signal22, th = signal_and_theory(freq2, freq2, config)
(Dict("TE" => [0.0, 0.0, 5.94969686942062, 2.7509157103624298, 1.474579711071198, 0.8484170459518878, 0.5133839618849767, 0.3259405744091539, 0.21810157931696741, 0.15439615404796334  …  -3.371934880663754e-6, -3.372986717016007e-6, -3.3745077231190236e-6, -3.376587070494609e-6, -3.3790630196926563e-6, -3.38218504386031e-6, -3.385651323499353e-6, -3.389691055273055e-6, -3.3940830169080422e-6, -3.3989262130620895e-6], "ET" => [0.0, 0.0, 5.94969686942062, 2.7509157103624298, 1.474579711071198, 0.8484170459518878, 0.5133839618849767, 0.3259405744091539, 0.21810157931696741, 0.15439615404796334  …  -3.371934880663754e-6, -3.372986717016007e-6, -3.3745077231190236e-6, -3.376587070494609e-6, -3.3790630196926563e-6, -3.38218504386031e-6, -3.385651323499353e-6, -3.389691055273055e-6, -3.3940830169080422e-6, -3.3989262130620895e-6], "EE" => [0.0, 0.0, 1.259572042012656, 0.48455705177483216, 0.24335423862793176, 0.14094070573132123, 0.08980646519341885, 0.06148186807264802, 0.044446533787336574, 0.033457540787315714  …  8.457772862047944e-6, 8.454204955605177e-6, 8.450619967918835e-6, 8.447078141478032e-6, 8.44344905119932e-6, 8.439803028508405e-6, 8.43606004534395e-6, 8.43239035465071e-6, 8.428573937589084e-6, 8.424740954615954e-6], "TT" => [0.0, 0.0, 1303.6832491814655, 588.6901144664341, 326.77976561475975, 205.70604152979342, 141.0400847789298, 102.84165006948633, 78.4819516937242, 62.05775766063433  …  0.00010789386244007327, 0.00010757534769610053, 0.00010725981794470818, 0.00010694525973150039, 0.00010663297428151154, 0.00010632285628185678, 0.00010601590174247024, 0.00010570960313242752, 0.00010540485854981548, 0.00010510376143781316]), Dict("TE" => [0.0, 0.0, 2.7410710061836303, 1.538364618634338, 0.8666585989752018, 0.49257031215634367, 0.28364841670661556, 0.167248294900734, 0.10292206599010562, 0.06757831087917948  …  -2.958512853138577e-6, -2.9598952948706985e-6, -2.961746509944255e-6, -2.9641556705145465e-6, -2.966961037763696e-6, -2.970412085469816e-6, -2.9742069947643985e-6, -2.978574962939165e-6, -2.9832947683479433e-6, -2.988465416274458e-6], "ET" => [0.0, 0.0, 2.7410710061836303, 1.538364618634338, 0.8666585989752018, 0.49257031215634367, 0.28364841670661556, 0.167248294900734, 0.10292206599010562, 0.06757831087917948  …  -2.958512853138577e-6, -2.9598952948706985e-6, -2.961746509944255e-6, -2.9641556705145465e-6, -2.966961037763696e-6, -2.970412085469816e-6, -2.9742069947643985e-6, -2.978574962939165e-6, -2.9832947683479433e-6, -2.988465416274458e-6], "EE" => [0.0, 0.0, 0.03234028781433917, 0.020781792483129162, 0.01083730084967642, 0.004836816993417869, 0.0019374949892939108, 0.0007854228473397264, 0.0003928177640878598, 0.00025154872003633637  …  2.9108320503513206e-6, 2.9116999222465746e-6, 2.912545394219194e-6, 2.9134287172579747e-6, 2.914219474762186e-6, 2.914988006623316e-6, 2.9156542932289212e-6, 2.916388595955796e-6, 2.9169709043795734e-6, 2.9175313873555934e-6], "TT" => [0.0, 0.0, 1064.7171662281169, 504.6062772110218, 286.70425884072733, 183.1500345955494, 126.93664956689628, 93.35938208515364, 71.75877586343385, 57.093419930748745  …  7.920122279308883e-5, 7.890565298286018e-5, 7.861304065330096e-5, 7.832137239398163e-5, 7.803194947381507e-5, 7.774466662370739e-5, 7.746051989729621e-5, 7.717700180727358e-5, 7.689501049525576e-5, 7.661663944786449e-5]))

Let's take a look at the results of the simulations compared to the signal.

if "--plot" ∈ ARGS
    plot(cl_mean[1:lmax_planck+1] .* collect(0:lmax_planck).^2,
        yerr=(sqrt.(cl_var[1:lmax_planck+1])  .* collect(0:lmax_planck).^2),
        label="sims", alpha=0.5)
    plot!(signal12[spec] .* collect(0:length(signal11[spec])-1).^2, xlim=(0,2nside),
        label="injected signal")
end

Next, we'll estimate the analytic signal-only covariance matrix. This requires the masks, the signal spectrum, and some zero arrays for the pixel variances.

X, Y = Symbol(spec[1]), Symbol(spec[2])
run_name = config["general"]["name"]
masktype1 = (X == :T) ? "T" : "P"
masktype2 = (Y == :T) ? "T" : "P"
mapid1 = "P$(freq1)hm$(splits[1])"
mapid2 = "P$(freq2)hm$(splits[2])"

maskfileT₁ = joinpath(config["scratch"], "masks", "$(run_name)_$(mapid1)_maskT.fits")
maskfileP₁ = joinpath(config["scratch"], "masks", "$(run_name)_$(mapid1)_maskP.fits")
maskfileT₂ = joinpath(config["scratch"], "masks", "$(run_name)_$(mapid2)_maskT.fits")
maskfileP₂ = joinpath(config["scratch"], "masks", "$(run_name)_$(mapid2)_maskP.fits")

maskT₁ = readMapFromFITS(maskfileT₁, 1, Float64)
maskP₁ = readMapFromFITS(maskfileP₁, 1, Float64)
maskT₂ = readMapFromFITS(maskfileT₂, 1, Float64)
maskP₂ = readMapFromFITS(maskfileP₂, 1, Float64);

zero_var_component = HealpixMap{Float64, RingOrder}(zeros(nside2npix(nside)))
zero_var = PolarizedHealpixMap(zero_var_component, zero_var_component, zero_var_component);

m1_signal = CovField(mapid1, maskT₁, maskP₁, zero_var)
m2_signal = CovField(mapid2, maskT₂, maskP₂, zero_var)


# convert a Vector (starting from 0) into a SpectralVector (i.e. a 0-indexed vector)
function format_signal(v::Vector{T}, nside) where T
    result = SpectralVector(zeros(T, nside2lmax(nside)+1))
    last_ind = min(length(v), length(result)) - 1
    for ℓ in 0:last_ind
        result[ℓ] = v[ℓ + 1]
    end
    return result
end


spectra = Dict{SpectrumName, SpectralVector{Float64, Vector{Float64}}}(
    (:TT, mapid1, mapid1) => format_signal(signal11["TT"], nside),
    (:TT, mapid1, mapid2) => format_signal(signal12["TT"], nside),
    (:TT, mapid2, mapid1) => format_signal(signal21["TT"], nside),
    (:TT, mapid2, mapid2) => format_signal(signal22["TT"], nside),

    (:EE, mapid1, mapid1) => format_signal(signal11["EE"], nside),
    (:EE, mapid1, mapid2) => format_signal(signal12["EE"], nside),
    (:EE, mapid2, mapid1) => format_signal(signal21["EE"], nside),
    (:EE, mapid2, mapid2) => format_signal(signal22["EE"], nside),

    (:TE, mapid1, mapid1) => format_signal(signal11["TE"], nside),
    (:TE, mapid1, mapid2) => format_signal(signal12["TE"], nside),
    (:TE, mapid2, mapid1) => format_signal(signal21["TE"], nside),
    (:TE, mapid2, mapid2) => format_signal(signal22["TE"], nside),
);

workspace_signal = CovarianceWorkspace(m1_signal, m2_signal, m1_signal, m2_signal);

@time C = coupledcov(Symbol(spec), Symbol(spec), workspace_signal, spectra)

symspec = (spec == "EE") ? :M⁺⁺ : Symbol(spec)
mask₁ = (spec[1] == 'T') ? maskT₁ : maskP₁
mask₂ = (spec[2] == 'T') ? maskT₂ : maskP₂
@time 𝐌 = mcm(symspec, mask₁, mask₂)
C_decoupled = decouple_covmat(C, 𝐌, 𝐌)
768×768 PowerSpectra.SpectralArray{Float64, 2, Matrix{Float64}} with indices 0:767×0:767:
     1.68222e5       661.633       …   7.56077e-5    7.80722e-5
   661.633         12300.7             1.48415e-5    9.78514e-6
    -5.27774e5     -1722.57           -0.000231177  -0.000238227
  -275.604        -54746.8            -6.24081e-5   -4.04927e-5
 33485.7              63.8839          2.36566e-6    9.7496e-7
   274.778          2661.79        …  -3.00607e-6   -2.36023e-6
  7624.72             51.1977          1.48637e-6    1.3182e-6
   187.203           865.944          -3.04772e-7   -2.90548e-7
   264.239            19.0491         -3.85263e-7   -4.06246e-7
   113.23             59.8835         -3.20625e-7   -2.59605e-7
     ⋮                             ⋱                
     0.000151648       2.00745e-5     -1.27006e-9   -1.19048e-9
     0.000140257       2.35652e-5  …  -2.48684e-9   -1.27128e-9
     0.000155043       2.18833e-5     -1.50434e-9   -2.48796e-9
     0.00014011        2.31245e-5     -4.95303e-8   -1.50747e-9
     0.000156289       2.27197e-5     -2.09576e-9   -4.95174e-8
     0.000159278       2.50791e-5     -2.35794e-7   -2.15994e-9
     0.000174445       1.22594e-5  …  -2.2741e-9    -2.35708e-7
     7.56077e-5        1.48415e-5      1.98377e-6   -2.57731e-9
     7.80722e-5        9.78514e-6     -2.57731e-9    1.98288e-6

Once a covariance matrix has been estimated, it needs to be corrected for mode-coupling just like the spectra. The covariance is possibly transposed, but we only care about the diagonal.

correction = SpectralVector(cl_var[:,1] ./ diag(parent(C_decoupled)))
768-element PowerSpectra.SpectralVector{Float64, Vector{Float64}} with indices 0:767:
 0.5085376068387355
 0.6276401494475211
 0.4598645614839412
 0.5910420663000585
 0.8939125978283793
 0.7087940102599101
 1.1491497212592383
 0.8982054699336028
 1.1372377130239266
 1.0526959365593267
 ⋮
 1.27757174042366
 1.2141564109948548
 1.2793711133043584
 1.4521272200022348
 1.1625663440691654
 1.2003946266076808
 1.0368921791953223
 1.023613251148139
 1.0021097983081189

Now it's finally time to fit a smooth interpolation over the noisy correction we just estimated.

lmin_cut, lmax_cut = get_plic_ellrange(Symbol(spec), freq1, freq2)
correction[lmax_cut:end] .= 1.0

fit_ells = intersect((lmin_cut-20):(lmax_cut), 0:lmax)
u = correction[fit_ells] .- 1
t = 1.0 * collect(fit_ells)
758-element Vector{Float64}:
  10.0
  11.0
  12.0
  13.0
  14.0
  15.0
  16.0
  17.0
  18.0
  19.0
   ⋮
 759.0
 760.0
 761.0
 762.0
 763.0
 764.0
 765.0
 766.0
 767.0

We use a really generic Gaussian kernel.

se = SE(0.0, 0.0)
m = MeanZero()
gp = GP(t,u,m,se)

using Optim
optimize!(gp)   # Optimise the hyperparameters
correction_gp = predict_y(gp, float.(0:(lmax_cut+20)))[1];
correction_gp[1:(lmin_cut-21)] .= 0.0

correction_path = joinpath(config["scratch"], "point_source_corrections")
mkpath(correction_path)

if transposed
    revspec = reverse(spec)
    correction_file = joinpath(correction_path, "$(freq1)_$(freq2)_$(revspec)_corr.dat")
else
    correction_file = joinpath(correction_path, "$(freq1)_$(freq2)_$(spec)_corr.dat")
end
writedlm(correction_file, correction_gp)

This smooth curve will be applied to the diagonal of the covariance matrices later.

if "--plot" ∈ ARGS
    plot(correction_gp)
end