jmert / cmb.jl Goto Github PK

Because of the dependency to SphericalHarmonicTransforms.jl which isn't registered

(@v1.7) pkg> add https://github.com/jmert/CMB.jl
     Cloning git-repo `https://github.com/jmert/CMB.jl`
    Updating git-repo `https://github.com/jmert/CMB.jl`
    Updating registry at `~/.julia/registries/General`
    Updating git-repo `https://github.com/JuliaRegistries/General.git`
   Resolving package versions...
ERROR: Unsatisfiable requirements detected for package SphericalHarmonicTransforms [b59b3b60]:
 SphericalHarmonicTransforms [b59b3b60] log:
 ├─SphericalHarmonicTransforms [b59b3b60] has no known versions!
 └─restricted to versions 0.1.1-0.1 by CMB [592991c5] — no versions left
   └─CMB [592991c5] log:
     ├─possible versions are: 0.3.2 or uninstalled
     └─CMB [592991c5] is fixed to version 0.3.2

Currently I have some non-general code for computing outer products of sparse vectors, but the functionality should really exist in Julia itself. I've opened JuliaLang/julia#24980 to work on this.

If[/when] accepted, the functionality should be copied/backported as a version-limited code block here to support the necessary computations even on julia v0.6.

The JuliaLang/julia#36107 PR changed how type aliases are printed, and that has the consequence that every occurrence of a doctest that shows a static vector has changed from a summary line like

julia> CMB.Sphere.ẑ
3-element SArray{Tuple{3},Int64,1,3} with indices SOneTo(3):
...

to

julia> CMB.Sphere.ẑ
3-element SVector{3,Int64} with indices SOneTo(3):
...

The nightly CI doctests are now failing due to this change, so before 1.6 gets branched this issue will need to be sorted out.

Maybe a long shot, have you ever thought about using different grids than equi-distant cylindrical (i.e. lon-lat) grids? In weather forecasting, some models also use

• reduced Gaussian grids: Fewer zonal points towards the poles to correct for the convergence of the longitudes which cause grid points to cover smaller areas towards the poles. However, I'm not sure how "fewer" is determined actually.
• octahedral grids. They are actually mathematically quite beautiful, see picture (from ECMWF)

They are constructed as follows. Get the Gaussian latitudes

julia> using FastGaussQuadrature
julia> nlat = 32
julia> latd = -asind.(FastGaussQuadrature.gausslegendre(nlat)[1])[1:nlat÷2]
16-element Vector{Float64}:
 85.7605871204438
 80.26877907225
 74.7445403686358
 69.21297616937085
 63.67863556109686
 58.14295404920328
 52.606526034345265
 47.06964205968768
 41.53246124665608
 35.99507841127159
 30.457553961152094
 24.91992862994861
 19.38223134643438
 13.84448373438486
  8.306702856518806
  2.76890300773601

Then start with 20 points for the latitude closest to the pole and increase by 4 points for every latitude towards the Equator

julia> lons = [0:360/(20+4i):360-180/(20+4i) for i in 0:length(latd)-1]

Some of the Fourier transforms get an ugly length

julia> using Primes
julia> nlons = [length(lon) for lon in lons]
julia> [factor(nlon) for nlon in nlons]
16-element Vector{Primes.Factorization{Int64}}:
 2^2 * 5
 2^3 * 3
 2^2 * 7
 2^5
 2^2 * 3^2
 2^3 * 5
 2^2 * 11
 2^4 * 3
 2^2 * 13
 2^3 * 7
 2^2 * 3 * 5
 2^6
 2^2 * 17
 2^3 * 3^2
 2^2 * 19
 2^4 * 5

but ECMWF uses FFTW successfully for that.

Hey Justin, many thanks for all your blog articles and AssociatedLegendrePolynomials.jl, great work! We are currently working on SpeedyWeather.jl, a spectral model of the atmosphere. You may have guessed that it also uses a spectral transform, much like you've defined the synthesize_* functions here. We have inherited the numerics from a decade old Fortran model, but want to overhaul the implementation completely. As I have learned loads from your blog articles, especially the one on "Notes on Calculating the Spherical Harmonics" I was wondering whether you have already draft/intent to write a similar one for the analysis, i.e. the inverse transform to the synthesis you've excellently described?

Through call overloading, a nicer interface to the cache-enabled Legendre calculations can be provided.

The trick should also work for the polarization covariance cache/calculations as well.

Hey Justin, was just wondering whether it makes sense to move your existing code to define the HEALPix grid and functions for it to a separate package, say HEALPix.jl. I see that you are defining a lot of functions for the grid, but I was wondering whether it also makes sense to create a struct HEALPixGrid{T,H,K} that efficiently stores data on that grid, something like

abstract type AbstractHEALPixGrid end    # not sure whether it makes sense to subtype it from any other existing array

struct HEALPixGrid{T,H,K} <: AbstractHEALPixGrid
    data::Vector{T}
    ...
end

# make H=4, K=3 the default
HEALPixGrid{T}(data::AbstractVector,...) = HEALPixGrid{eltype{data},4,3}(data,...)

then we can more easily dispatch functions and clearly define projections onto that grid and back. What do you think?

Take the following setup code:

using PyPlot
import CMB.Legendre: LegendreSphereCoeff

function integrate_simpson(f, xs)
    N = length(xs)
    I = copy(f(xs[1]))
    I .+= f(xs[end])
    for x in xs[2:2:end-1]
        I .+= oftype(x, 2)/3 .* f(x)
    end
    for x in xs[3:2:end-1]
        I .+= oftype(x, 4)/3 .* f(x)
    end
    I .*= step(xs)
    return I
end

# _sinc(x) == sinc(x/π) without the extra division
_sinc(x) = iszero(x) ? one(sin(zero(x))) : sin(x) / x

function pixel_window_function(θmin, θmax, φmin, φmax;
        lmax::Int=720, nstep::Int=100, eltype::Type=Float64)

    nstep += mod(nstep, 2) # Forces nstep to be even
    ell = 0:lmax
    zmin, zmax = cos.(convert.(eltype, (θmax, θmin)))
    zs = range(zmin, stop=zmax, length=nstep)
    Δz, Δφ = eltype.((zmax - zmin, φmax - φmin))

    Λfn(z) = let Λtmp = zeros(eltype, lmax+1, lmax+1),
                 Λtab! = LegendreSphereCoeff{eltype}(lmax)
        Λtab!(Λtmp, z)
    end
    Λ = integrate_simpson(Λfn, zs)
    Λ[:,2:end] .*= sqrt(2) # double to account for m ≤ -1

    Φ(m) = (Δφ * _sinc(m * Δφ/2))^2

    wlm2 = squeeze(sum(Λ.^2 .* Φ.(ell'), dims=2), dims=2)
    norm = 4π / (Δφ * Δz)^2
    Wl = sqrt.(norm ./ (2 .* ell .+ 1) .* wlm2)
    return (ell, Wl, Λ)
end

Computing up to ℓ_max = 2160 for the pixel settings as follow, the answer starts diverging somwhere in the range ℓ ∈ [2120, 2125].

Δδ, Δλ = (0.250,   0.466)
δ₀, λ₀ = (-57.625, 0.000)

θlim = deg2rad.(90.0 .- (δ₀ .+ (Δδ .* (1, -1))))
φlim = deg2rad.(λ₀ .+ (Δλ .* (-1, 1)))
(l, Wl, Λ) = pixel_window_function(θlim..., φlim..., lmax=3*720)

semilogy(l, Wl)

Looking at the Legendre integral term Λ directly, for ℓ = 2110 and ℓ = 2125 (x-axis being across m) we can see what is probably contributing to the divergence as a hump beyond the final "ring-down".

plot(l, Λ[2110, :], label="ℓ = 2110")
plot(l, Λ[2125, :], label="ℓ = 2125")
legend()
xlabel(raw"$m$")
ylabel(raw"$\int\lambda_\ell^m(z)\,dz$")

Repeating the same exercises with eltype = BigFloat do not show the same behavior, so the reasonable conclusion is that there's either round-off or term-by-term cancellation erros which are contributing to the breakdown for double-float precisions.

Other observations:

• Switching from Simpson's rule to trapezoidal rule integration has no impact, so I'm interpreting this as numerical effects from the Legendre polynomials, not overshoot/"ringing" induced by the integration method.

The parallelism check coded in

In commits 668184d and b4fb922, the limiting behavior of the vector-based bearing angle calculations were updated to give continuous answers as (1) r₁ approached the poles, and (2) as r₁ and r₂ become (anti)-parallel.

The angle-based calculations have not been updated, to agree with this behavior. Fix this. Solutions could be:

• Update the limiting-case logic to handle these two cases.
• Get rid of the implementation and simply use the vector-based calculations by converting angles to unit vectors internally.