eschnett/SHTOOLS.jl - Change X5ISFPUVKYET33KOY3IN3FSI54BWY2JDPLGBAK5ENTAGBH54NK4QC

Implement all spherical harmonic transform functions

Created by Erik Schnetter on February 17, 2021

X5ISFPUVKYET33KOY3IN3FSI54BWY2JDPLGBAK5ENTAGBH54NK4QC

Dependencies

In channels

main

Change contents

Insertion in Project.toml at line 7 [3.3694]
[3.3835]
[3.3835]
```
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
```

Insertion in README.md at line 57 [3.3301]

[2.39]

[4.281]


SHExpandLSQ
MakeGrid2d
MakeGridPoint
MakeGridPointC
SHMultiply

Insertion in src/SHTOOLS.jl at line 741 [3.421]

[2.6439]

end
# Other routines
export SHExpandLSQ!
function SHExpandLSQ!(cilm::AbstractArray{Cdouble,3},
                      d::AbstractVector{Cdouble}, lat::AbstractVector{Cdouble},
                      lon::AbstractVector{Cdouble}, nmax::Integer,
                      lmax::Integer; norm::Integer=1, csphase::Integer=1,
                      weights::Vector{Cdouble},
                      exitstatus::Optional{Ref{<:Integer}}=nothing)
    @assert lmax ≥ 0
    @assert size(cilm, 1) == 2
    @assert size(cilm, 2) ≥ lmax + 1
    @assert size(cilm, 3) == size(cilm, 2)
    @assert nmax ≥ 0
    @assert nmax ≥ (lmax + 1)^2
    @assert length(d) ≥ nmax
    @assert length(lat) ≥ nmax
    @assert length(lon) ≥ nmax
    @assert norm ∈ (1, 2, 3, 4)
    @assert length(weights) == nmax
    chi2′ = Ref{Cdouble}()
    exitstatus′ = Ref{Cint}()
    ccall((:SHExpandLSQ, libSHTOOLS), Cvoid,
          (Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cint,
           Cint, Ref{Cint}, Ref{Cdouble}, Ref{Cint}, Ptr{Cdouble}, Ref{Cint}),
          cilm, size(cilm, 2), d, lat, lon, nmax, lmax, norm, chi2′, csphase,
          weights, exitstatus′)
    if exitstatus === nothing
        exitstatus′[] ≠ 0 && error("SHExpandLSQ!: Error code $(exitstatus′[])")
    else
        exitstatus[] = exitstatus′[]
    end
    return cilm, Cdouble(chi2′[])
end
export SHExpandLSQ
function SHExpandLSQ(d::AbstractVector{Cdouble}, lat::AbstractVector{Cdouble},
                     lon::AbstractVector{Cdouble}, nmax::Integer, lmax::Integer;
                     norm::Integer=1, csphase::Integer=1,
                     weights::Vector{Cdouble},
                     exitstatus::Optional{Ref{<:Integer}}=nothing)
    @assert lmax ≥ 0
    cilm = Array{Cdouble}(undef, 2, lmax + 1, lmax + 1)
    _, chi2 = SHExpandLSQ!(cilm, d, lat, lon, nmax, lmax; norm=norm,
                           csphase=csphase, weights=weights, exitstatus)
    return cilm, chi2
end
export MakeGrid2d!
function MakeGrid2d!(grid::AbstractArray{Cdouble,2},
                     cilm::AbstractArray{Cdouble,3}, lmax::Integer,
                     interval::Union{AbstractFloat,Integer}; norm::Integer=1,
                     csphase::Integer=1,
                     f::Optional{Union{AbstractFloat,Integer}}=nothing,
                     a::Optional{Union{AbstractFloat,Integer}}=nothing,
                     north::Union{AbstractFloat,Integer}=90.0,
                     south::Union{AbstractFloat,Integer}=-90.0,
                     east::Union{AbstractFloat,Integer}=360.0,
                     west::Union{AbstractFloat,Integer}=0.0,
                     dealloc::Bool=false,
                     exitstatus::Optional{Ref{<:Integer}}=nothing)
    @assert interval > 0
    @assert size(grid, 1) ≥ floor(Int, 180 / interval + 1)
    @assert size(grid, 2) ≥ floor(Int, 360 / interval + 1)
    @assert size(cilm, 1) == 2
    @assert lmax ≥ 0
    @assert size(cilm, 2) ≥ lmax + 1
    @assert size(cilm, 3) == size(cilm, 2)
    @assert (f !== nothing) == (a !== nothing)
    nlat = Ref{Cint}()
    nlong = Ref{Cint}()
    exitstatus′ = Ref{Cint}()
    ccall((:MakeGrid2d, libSHTOOLS), Cvoid,
          (Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Cint, Cint, Cdouble,
           Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cdouble},
           Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
           Ref{Cint}, Ref{Cint}), grid, size(grid, 1), size(grid, 2), cilm,
          size(cilm, 2), lmax, interval, nlat, nlong, norm, csphase,
          optional(f, Ptr{Cdouble}()), optional(a, Ptr{Cdouble}()), north,
          south, east, west, dealloc, exitstatus′)
    if exitstatus === nothing
        exitstatus′[] ≠ 0 && error("MakeGrid2d!: Error code $(exitstatus′[])")
    else
        exitstatus[] = exitstatus′[]
    end
    return grid, Int(nlat[]), Int(nlong[])
end
export MakeGrid2d
function MakeGrid2d(cilm::AbstractArray{Cdouble,3}, lmax::Integer,
                    interval::Union{AbstractFloat,Integer}; norm::Integer=1,
                    csphase::Integer=1,
                    f::Optional{Union{AbstractFloat,Integer}}=nothing,
                    a::Optional{Union{AbstractFloat,Integer}}=nothing,
                    north::Union{AbstractFloat,Integer}=90.0,
                    south::Union{AbstractFloat,Integer}=-90.0,
                    east::Union{AbstractFloat,Integer}=360.0,
                    west::Union{AbstractFloat,Integer}=0.0, dealloc::Bool=false,
                    exitstatus::Optional{Ref{<:Integer}}=nothing)
    @assert interval > 0
    grid = Array{Cdouble}(undef, floor(Int, 180 / interval + 1),
                          floor(Int, 360 / interval + 1))
    MakeGrid2d!(grid, cilm, lmax, interval; norm=norm, csphase=csphase, f=f,
                a=a, north=north, south=south, east=east, west=west,
                dealloc=dealloc, exitstatus=exitstatus)
    return grid
end
export MakeGridPoint
function MakeGridPoint(cilm::AbstractArray{Cdouble,3}, lmax::Integer,
                       lat::Union{AbstractFloat,Integer},
                       lon::Union{AbstractFloat,Integer}; norm::Integer=1,
                       csphase::Integer=1, dealloc::Bool=false)
    @assert lmax ≥ 0
    @assert size(cilm, 1) == 2
    @assert size(cilm, 2) ≥ lmax + 1
    @assert size(cilm, 3) ≥ size(cilm, 2)
    @assert norm ∈ (1, 2, 3, 4)
    @assert csphase ∈ (1, -1)
    value = ccall((:MakeGridPoint, libSHTOOLS), Cdouble,
                  (Ptr{Cdouble}, Cint, Cint, Cdouble, Cdouble, Ref{Cint},
                   Ref{Cint}, Ref{Cint}), cilm, size(cilm, 2), lmax, lat, lon,
                  norm, csphase, dealloc)
    return Float64(value)

Insertion in src/SHTOOLS.jl at line 869 [3.421]

[2.6444]

export MakeGridPointC
function MakeGridPointC(cilm::AbstractArray{Complex{Cdouble},3}, lmax::Integer,
                        lat::Union{AbstractFloat,Integer},
                        lon::Union{AbstractFloat,Integer}; norm::Integer=1,
                        csphase::Integer=1, dealloc::Bool=false)
    @assert lmax ≥ 0
    @assert size(cilm, 1) == 2
    @assert size(cilm, 2) ≥ lmax + 1
    @assert size(cilm, 3) ≥ size(cilm, 2)
    @assert norm ∈ (1, 2, 3, 4)
    @assert csphase ∈ (1, -1)
    value = ccall((:MakeGridPointC, libSHTOOLS), Complex{Cdouble},
                  (Ptr{Complex{Cdouble}}, Cint, Cint, Cdouble, Cdouble,
                   Ref{Cint}, Ref{Cint}, Ref{Cint}), cilm, size(cilm, 2), lmax,
                  lat, lon, norm, csphase, dealloc)
    return Complex{Float64}(value)

Insertion in src/SHTOOLS.jl at line 886 [3.421]

[2.6448]


export SHMultiply!
function SHMultiply!(cilmout::AbstractArray{Cdouble,3},
                     cilm1::AbstractArray{Cdouble,3}, lmax1::Integer,
                     cilm2::AbstractArray{Cdouble,3}, lmax2::Integer;
                     precomp::Bool=false, norm::Integer=1, csphase::Integer=1,
                     exitstatus::Optional{Ref{<:Integer}}=nothing)
    @assert lmax1 ≥ 0
    @assert lmax2 ≥ 0
    @assert size(cilmout, 1) == 2
    @assert size(cilmout, 2) ≥ lmax1 + lmax2 + 1
    @assert size(cilmout, 3) == size(cilmout, 2)
    @assert size(cilm1, 1) == 2
    @assert size(cilm1, 2) ≥ lmax1 + 1
    @assert size(cilm1, 3) == size(cilm1, 2)
    @assert size(cilm2, 1) == 2
    @assert size(cilm2, 2) ≥ lmax2 + 1
    @assert size(cilm2, 3) == size(cilm2, 2)
    @assert norm ∈ (1, 2, 3, 4)
    @assert csphase ∈ (1, -1)
    exitstatus′ = Ref{Cint}()
    ccall((:SHMultiply, libSHTOOLS), Cvoid,
          (Ptr{Cdouble}, Cint, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Cint,
           Cint, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}), cilmout,
          size(cilmout, 2), cilm1, size(cilm1, 2), lmax1, cilm2, size(cilm2, 2),
          lmax2, precomp, norm, csphase, exitstatus′)
    if exitstatus === nothing
        exitstatus′[] ≠ 0 && error("SHMultiply!: Error code $(exitstatus′[])")
    else
        exitstatus[] = exitstatus′[]
    end
    return cilmout
end
export SHMultiply
function SHMultiply(cilm1::AbstractArray{Cdouble,3}, lmax1::Integer,
                    cilm2::AbstractArray{Cdouble,3}, lmax2::Integer;
                    precomp::Bool=false, norm::Integer=1, csphase::Integer=1,
                    exitstatus::Optional{Ref{<:Integer}}=nothing)
    @assert lmax1 ≥ 0
    @assert lmax2 ≥ 0
    cilmout = Array{Cdouble}(undef, 2, lmax1 + lmax2 + 1, lmax1 + lmax2 + 1)
    SHMultiply!(cilmout, cilm1, lmax1, cilm2, lmax2; precomp=precomp, norm=norm,
                csphase=csphase, exitstatus=exitstatus)
    return cilmout
end
end

Insertion in test/runtests.jl at line 1 [3.22]
[3.22]
[3.23]
```
using Random
```
Insertion in test/runtests.jl at line 116 [3.22]
[5.27]
[4.10347]
```
    Random.seed!(0)
```

Replacement in test/runtests.jl at line 124 [3.22]

B:BD[4.10443] → [5.59:114]

    # Go back to the grid (it'll be low-pass filtered)

[4.10443]

    # Go back to the grid (it will be low-pass filtered)

Insertion in test/runtests.jl at line 141 [3.22]
[6.4633]
[6.4633]
```
    Random.seed!(0)
```

Replacement in test/runtests.jl at line 149 [3.22]

B:BD[6.4778] → [6.4778:4833]

    # Go back to the grid (it'll be low-pass filtered)

[6.4778]

[6.4833]

    # Go back to the grid (it will be low-pass filtered)

Insertion in test/runtests.jl at line 165 [3.22]
[2.6494]
[7.5318]
```
    Random.seed!(0)
```

Replacement in test/runtests.jl at line 176 [3.22]

B:BD[7.5525] → [7.5525:5580]

    # Go back to the grid (it'll be low-pass filtered)

[7.5525]

[7.5580]

    # Go back to the grid (it will be low-pass filtered)

Replacement in test/runtests.jl at line 190 [3.22]

B:BD[2.6536] → [2.6536:6596]

    @show (length(latglq), length(longlq)) == size(gridglq)

[2.6536]

[8.6445]

    @test (length(latglq), length(longlq)) == size(gridglq)

Insertion in test/runtests.jl at line 194 [3.22]
[2.6645]
[2.6645]
```
    Random.seed!(0)
```

Replacement in test/runtests.jl at line 205 [3.22]

B:BD[2.6871] → [2.6871:6926]

    # Go back to the grid (it'll be low-pass filtered)

[2.6871]

[2.6926]

    # Go back to the grid (it will be low-pass filtered)

Insertion in test/runtests.jl at line 217 [3.22]

[2.7400]

end
@testset "SHExpandLSQ" begin
    Random.seed!(0)
    lmax = 2
    n = (lmax + 1)^2
    lat = π * rand(n)
    lon = 2π * rand(n)
    d = randn(n)
    weights = Float64[1 for i in 1:n]
    cilm, chi2 = SHExpandLSQ(d, lat, lon, n, lmax; weights=weights)
    @test size(cilm) == (2, lmax + 1, lmax + 1)
    @test abs(chi2) < 10 * eps(Float64)
end
@testset "MakeGrid2d" begin
    Random.seed!(0)
    lmax = 4
    cilm = randn(2, lmax + 1, lmax + 1)
    interval = 15                # grid spacing in degrees
    grid = MakeGrid2d(cilm, lmax, interval)
    @test size(grid) == (180 ÷ interval + 1, 360 ÷ interval + 1)
end
@testset "MakeGridPoint" begin
    Random.seed!(0)
    # Choose random values
    lmax = 2
    nmax = (lmax + 1)^2
    lat = π * rand(nmax)
    lon = 2π * rand(nmax)
    values = randn(nmax)
    weights = Float64[1 for n in 1:nmax]
    # Expand
    cilm, chi2 = SHExpandLSQ(values, lat, lon, nmax, lmax; weights=weights)
    @test size(cilm) == (2, lmax + 1, lmax + 1)
    @test abs(chi2) < 10 * eps(Float64)
    # Convert back to values (they will be low-pass filtered)
    values′ = Float64[MakeGridPoint(cilm, lmax, lat[n], lon[n]) for n in 1:nmax]
    # Expand again (result must be the same)
    cilm′, chi2 = SHExpandLSQ(values′, lat, lon, nmax, lmax; weights=weights)
    @test abs(chi2) < 10 * eps(Float64)
    @test isapprox(cilm′, cilm)
    # Convert back to values (result must be the same this time)
    values″ = Float64[MakeGridPoint(cilm′, lmax, lat[n], lon[n])
                      for n in 1:nmax]
    @test isapprox(values″, values′)

Insertion in test/runtests.jl at line 269 [3.22]

[2.7404]


@testset "MakeGridPointC" begin
    Random.seed!(0)
    # Choose random values
    lmax = 2
    cilm = randn(Complex{Float64}, 2, lmax + 1, lmax + 1)
    lat = π * rand()
    lon = 2π * rand()
    value = MakeGridPointC(cilm, lmax, lat, lon)
    @test value isa Complex{Float64}
end
@testset "SHMultiply" begin
    Random.seed!(0)
    # Choose random values
    lmax1 = 2
    lmax2 = 3
    cilm1 = randn(2, lmax1 + 1, lmax1 + 1)
    cilm2 = randn(2, lmax2 + 1, lmax2 + 1)
    cilmout = SHMultiply(cilm1, lmax1, cilm2, lmax2)
    @test size(cilmout) == (2, lmax1 + lmax2 + 1, lmax1 + lmax2 + 1)
end