text:
using VectorSphericalHarmonics
using Test
using Aqua
using LinearAlgebra
using LegendrePolynomials
using OffsetArrays
using SphericalHarmonics
using SphericalHarmonics: NorthPole, SouthPole, getY, eltypeY
using SphericalHarmonicModes
using WignerD
using HCubature
using VectorSphericalHarmonics: basisconversionmatrix, cache, cache!
using Rotations
using StaticArrays

@testset "project quality" begin
    if VERSION >= v"1.6"
        Aqua.test_all(VectorSphericalHarmonics, ambiguities = (recursive = false,))
    else
        Aqua.test_all(VectorSphericalHarmonics, ambiguities = false)
    end
end

isapproxdefault(a, b) = isapprox(a, b, atol = 1e-14, rtol = sqrt(eps(Float64)))
isapproxdefault(x::Tuple{Any, Any}) = isapproxdefault(x...)

function genspharm2(l, m, n, θ, ϕ)
    @assert abs(m) <= l "m must satisfy -l <= m <= l, received l = $l, m = $m"
    @assert abs(n) <= l "n must satisfy -l <= n <= l, received l = $l, n = $n"
    √((2l+1)/4π) * WignerD.wignerdjmn(l, m, n, θ) * cis(m * ϕ)
end

const S = cache(10);
const Ylm = getY(S);

@testset "Basis" begin
    for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
        for B2 in [HelicityCovariant(), Polar(), SphericalCovariant(), Cartesian()]
            for B1 in [HelicityCovariant(), SphericalCovariant()]
                M = basisconversionmatrix(B1, B2, θ, ϕ)
                Minv = basisconversionmatrix(B2, B1, θ, ϕ)
                # transformations must be unitary
                @test M' * M ≈ I
                @test Minv ≈ M'
                @test abs(det(M)) ≈ 1
            end
        end
    end
    @testset "multiple transformations" begin
        for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
            for B in [HelicityCovariant(), Polar(), SphericalCovariant(), Cartesian()]
                M1 = basisconversionmatrix(HelicityCovariant(), B, θ, ϕ)
                M2 = basisconversionmatrix(B, SphericalCovariant(), θ, ϕ)
                M3 = basisconversionmatrix(HelicityCovariant(), SphericalCovariant(), θ, ϕ)
                @test begin
                    res = M2 * M1 ≈ M3
                    if !res
                        println(HelicityCovariant(),"->",B,"->",SphericalCovariant(),", θ = $θ, ϕ = $ϕ")
                    end
                    res
                end

                M1 = basisconversionmatrix(SphericalCovariant(), B, θ, ϕ)
                M2 = basisconversionmatrix(B, HelicityCovariant(), θ, ϕ)
                M3 = basisconversionmatrix(SphericalCovariant(), HelicityCovariant(), θ, ϕ)
                @test begin
                    res = M2 * M1 ≈ M3
                    if !res
                        println(SphericalCovariant(),"->",B,"->",HelicityCovariant(),", θ = $θ, ϕ = $ϕ")
                    end
                    res
                end
            end
        end
    end
    @testset "basis at poles" begin
        @testset "Cartesian and polar" begin
            M = basisconversionmatrix(Cartesian(), Polar(), 0, 0)
            @test M * [1,0,0] ≈ [0,1,0] # x̂ == θ̂ at the north pole
            @test M * [0,1,0] ≈ [0,0,1] # ŷ == ϕ̂ at the north pole
            @test M * [0,0,1] ≈ [1,0,0] # ẑ == r̂ at the north pole

            M = basisconversionmatrix(Polar(), Cartesian(), pi, 0)
            @test M * [1,0,0] ≈ [0,0,-1] # r̂ == -ẑ at the south pole
            @test M * [0,1,0] ≈ [-1,0,0] # θ̂ == -x̂ at the south pole
            @test M * [0,0,1] ≈ [0,1,0]  # ϕ̂ == ŷ at the south pole
        end
    end
end

@testset "cache" begin
    θ, ϕ = pi/3, pi/4
    S2 = @inferred cache(Float64, 1)
    cache!(S2, θ, ϕ)
    @test getY(S2) == getY(cache(Float64, θ, ϕ, 1))

    S2 = @inferred cache(1)
    cache!(S2, θ, ϕ)
    @test getY(S2) == getY(cache(Float64, θ, ϕ, 1))

    @test vshbasis(PB(), Polar(), 1, 1, θ, ϕ, S2) == vshbasis(PB(), Polar(), 1, 1, θ, ϕ)
    @test genspharm(1, 1, θ, ϕ, S2) == genspharm(1, 1, θ, ϕ)
end

@testset "VSH orientations" begin
    θ1, ϕ1 = pi/3, pi/3;
    θNP, ϕNP = 0, 0;
    @testset "PB" begin
        @testset "HelicityCovariant" begin
            cache!(S, θ1, ϕ1, 0)
            Y = @inferred vshbasis(PB(), HelicityCovariant(), 0, 0, θ1, ϕ1, S)
            # for l=0 there is only one non-zero vector Ylm0, which is radial
            for n in axes(Y,2), basisind in axes(Y,1)
                if !(n == 0 && basisind == 0 #= e₀′ (radial) component =#)
                    @test isapproxdefault(0, Y[basisind, n])
                end
            end
        end
        @testset "Polar" begin
            cache!(S, θ1, ϕ1, 0)
            Y = @inferred vshbasis(PB(), Polar(), 0, 0, θ1, ϕ1, S)
            # for l=0 there is only one non-zero vector Ylm0, which is radial
            for n in axes(Y,2), basisind in axes(Y,1)
                if !(n == 0 && basisind == 1 #= radial component =#)
                    @test isapproxdefault(0, Y[basisind, n])
                end
            end
        end
        @testset "SphericalCovariant" begin
            # north pole
            cache!(S, θNP, ϕNP, 0)
            Y = @inferred vshbasis(PB(), SphericalCovariant(), 0, 0, θNP, ϕNP, S)
            # for l=0 there is only one non-zero vector Ylm0, which is radial
            # at north pole r̂ == ẑ. so there is only one non-zero component of Ylm0
            for n in axes(Y,2), basisind in axes(Y,1)
                if !(n == 0 && basisind == 0 #= e₀ (z) component =#)
                    @test isapproxdefault(0, Y[basisind, n])
                end
            end
        end
        @testset "Cartesian" begin
            # north pole
            cache!(S, θNP, ϕNP, 0)
            Y = @inferred vshbasis(PB(), Cartesian(), 0, 0, θNP, ϕNP, S)
            # for l=0 there is only one non-zero vector Ylm0, which is radial
            # at north pole r̂ == ẑ. so there is only one non-zero component of Ylm0
            for n in axes(Y,2), basisind in axes(Y,1)
                if !(n == 0 && basisind == 3 #= z component =#)
                    @test isapproxdefault(0, Y[basisind, n])
                end
            end
        end
    end
    @testset "Hansen" begin
        @testset "HelicityCovariant" begin
            cache!(S, θ1, ϕ1, 0)
            Y = @inferred vshbasis(Hansen(), HelicityCovariant(), 0, 0, θ1, ϕ1, S)
            # for l=0 there is only one non-zero vector Ylm⁽⁻¹⁾, which is radial
            for n in axes(Y,2), basisind in axes(Y,1)
                if !(n == -1 && basisind == 0 #= e₀′ (radial) component =#)
                    @test isapproxdefault(0, Y[basisind, n])
                end
            end
        end
        @testset "Polar" begin
            cache!(S, θ1, ϕ1, 0)
            Y = @inferred vshbasis(Hansen(), Polar(), 0, 0, θ1, ϕ1, S)
            # for l=0 there is only one non-zero vector Ylm⁽⁻¹⁾, which is radial
            for n in axes(Y,2), basisind in axes(Y,1)
                if !(n == -1 && basisind == 1 #= radial component =#)
                    @test isapproxdefault(0, Y[basisind, n])
                end
            end
        end
        @testset "SphericalCovariant" begin
            # north pole
            cache!(S, θNP, ϕNP, 0)
            Y = @inferred vshbasis(Hansen(), SphericalCovariant(), 0, 0, θNP, ϕNP, S)
            # for l=0 there is only one non-zero vector Ylm⁽⁻¹⁾, which is radial
            # at north pole r̂ == ẑ. so there is only one non-zero component of Ylm0
            for n in axes(Y,2), basisind in axes(Y,1)
                if !(n == -1 && basisind == 0 #= e₀ (z) component =#)
                    @test isapproxdefault(0, Y[basisind, n])
                end
            end
        end
        @testset "Cartesian" begin
            # north pole
            cache!(S, θNP, ϕNP, 0)
            Y = @inferred vshbasis(Hansen(), Cartesian(), 0, 0, θNP, ϕNP, S)
            # for l=0 there is only one non-zero vector Ylm⁽⁻¹⁾⁾, which is radial
            # at north pole r̂ == ẑ. so there is only one non-zero component of Ylm0
            for n in axes(Y,2), basisind in axes(Y,1)
                if !(n == -1 && basisind == 3 #= z component =#)
                    @test isapproxdefault(0, Y[basisind, n])
                end
            end
        end
    end
    @testset "Irreducible" begin
        cache!(S, θ1, ϕ1, 0)
        Y = @inferred vshbasis(Irreducible(), Polar(), 0, 0, θ1, ϕ1, S)
        @test all(iszero, Y[:, -1:0])
        @test all(x -> isapproxdefault(x, 0), Y[2:3, 1])
        @test Y[1, 1] ≈ -1/√(4pi)

        for m = -1:1
            Y = @inferred vshbasis(Irreducible(), SphericalCovariant(), 1, m, θ1, ϕ1, S)
            @test Y[m, 0] ≈ 1/√(4pi)
            @test all(iszero, Y[-1:m-1, 0])
            @test all(iszero, Y[m+1:1, 0])
        end

        cache!(S, θNP, ϕNP, 0)
        Y = @inferred vshbasis(Irreducible(), Cartesian(), 0, 0, θNP, ϕNP, S)
        @test all(iszero, Y[:, -1:0])
        @test all(iszero, Y[1:2, 1])
        @test Y[3, 1] ≈ -1/√(4pi)
    end
end

@testset "norm in different bases" begin
    lmax = 4
    for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
        cache!(S, θ, ϕ)
        for l = 0:lmax, m = -l:l
            Y = vshbasis(Irreducible(), SphericalCovariant(), l, m, θ, ϕ, S)
            Yp = parent(Y)
            for B in [HelicityCovariant(), Polar(), Cartesian()]
                Y2 = vshbasis(Irreducible(), B, l, m, θ, ϕ, S)
                Y2p = parent(Y2)
                for l in 1:3
                    @test norm(Y2p[:, l]) ≈ norm(Yp[:, l])
                end
            end
        end
    end
end

@testset "VSH components" begin
    lmax = 5
    @testset "HelicityCovariant" begin
        @testset "PB" begin
            for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
                cache!(S, θ, ϕ)
                for l in 0:lmax, m in -l:l
                    Y = vshbasis(PB(), HelicityCovariant(), l, m, θ, ϕ, S)
                    @test isapproxdefault(Y[0,0], Ylm[(l,m)])
                    for n in -1:1
                        if n in -min(1, l):min(1, l)
                            @test isapproxdefault(Y[n,n], genspharm2(l, m, n, θ, ϕ))
                        else
                            @test isapproxdefault(Y[n,n], 0)
                        end
                        for basisind in -1:1
                            basisind == n && continue
                            @test isapproxdefault(Y[basisind,n], 0)
                        end
                    end
                end
            end
        end
        @testset "Hansen" begin
            # l = 0 has only one component
            for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
                cache!(S, θ, ϕ, 0)
                Y = vshbasis(Hansen(), HelicityCovariant(), 0, 0, θ, ϕ, S)
                @test isapproxdefault(Y[ 0, -1],  √(1/4pi) * WignerD.wignerDjmn(0, 0, 0, 0, θ, ϕ))
                for n in -1:1, basisind in -1:1
                    basisind == 0 && n == -1 && continue
                    @test isapproxdefault(Y[basisind, n], 0)
                end
            end
            for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
                cache!(S, θ, ϕ)
                for l in 1:lmax, m in -l:l
                    norm = √((2l+1)/8pi)
                    sqrt2norm = √((2l+1)/4pi)
                    Y = vshbasis(Hansen(), HelicityCovariant(), l, m, θ, ϕ, S)
                    @test isapproxdefault(Y[-1, +1],  norm * WignerD.wignerDjmn(l,  1, -m, 0, θ, ϕ))
                    @test isapproxdefault(Y[ 0, +1],  0)
                    @test isapproxdefault(Y[+1, +1],  norm * WignerD.wignerDjmn(l, -1, -m, 0, θ, ϕ))
                    @test isapproxdefault(Y[-1,  0],  norm * WignerD.wignerDjmn(l,  1, -m, 0, θ, ϕ))
                    @test isapproxdefault(Y[ 0,  0],  0)
                    @test isapproxdefault(Y[+1,  0], -norm * WignerD.wignerDjmn(l, -1, -m, 0, θ, ϕ))
                    @test isapproxdefault(Y[-1, -1],  0)
                    @test isapproxdefault(Y[ 0, -1],  sqrt2norm * WignerD.wignerDjmn(l, 0, -m, 0, θ, ϕ))
                    @test isapproxdefault(Y[+1, -1],  0)
                end
            end
        end
        @testset "Irreducible" begin
            # j = 0 has only one component
            for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
                cache!(S, θ, ϕ)
                Y = vshbasis(Irreducible(), HelicityCovariant(), 0, 0, θ, ϕ, S)
                @test isapproxdefault(Y[0, +1], -√(1/4pi) * WignerD.wignerDjmn(0, 0, 0, 0, θ, ϕ))
                for n in -1:1, basisind in -1:1
                    basisind == 0 && n == 1 && continue
                    @test isapproxdefault(Y[basisind, n], 0)
                end
            end
            for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
                cache!(S, θ, ϕ)
                for j in 1:lmax, m in -j:j
                    Y = vshbasis(Irreducible(), HelicityCovariant(), j, m, θ, ϕ, S)
                    @test isapproxdefault(Y[-1, j+1], √(j/8pi) * WignerD.wignerDjmn(j,  1, -m, 0, θ, ϕ))
                    @test isapproxdefault(Y[ 0, j+1], -√((j+1)/4pi) * WignerD.wignerDjmn(j, 0, -m, 0, θ, ϕ))
                    @test isapproxdefault(Y[+1, j+1], √(j/8pi) * WignerD.wignerDjmn(j, -1, -m, 0, θ, ϕ))
                    @test isapproxdefault(Y[-1,  j], √((2j+1)/8pi) * WignerD.wignerDjmn(j,  1, -m, 0, θ, ϕ))
                    @test isapproxdefault(Y[ 0,  j], 0)
                    @test isapproxdefault(Y[+1,  j], -√((2j+1)/8pi) * WignerD.wignerDjmn(j, -1, -m, 0, θ, ϕ))
                    @test isapproxdefault(Y[-1, j-1],  √((j+1)/8pi) * WignerD.wignerDjmn(j,  1, -m, 0, θ, ϕ))
                    @test isapproxdefault(Y[ 0, j-1],  √(j/4pi) * WignerD.wignerDjmn(j, 0, -m, 0, θ, ϕ))
                    @test isapproxdefault(Y[+1, j-1],  √((j+1)/8pi) * WignerD.wignerDjmn(j, -1, -m, 0, θ, ϕ))
                end
            end
        end
    end
    @testset "Polar" begin
        @testset "Hansen" begin
            # l = 0 has only one component
            for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
                cache!(S, θ, ϕ, 0)
                Y = vshbasis(Hansen(), Polar(), 0, 0, θ, ϕ, S)
                @test begin
                    res = isapproxdefault(Y[ 1, -1], Ylm[(0,0)])
                    if !res
                        @show θ, ϕ
                    end
                    res
                end
                for n in axes(Y,2), basisind in axes(Y,1)
                    basisind == 1 && n == -1 && continue
                    @test isapproxdefault(Y[basisind, n], 0)
                end
            end
            for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
                cache!(S, θ, ϕ)
                for l in 1:lmax, m in -l:l
                    norml⁻ = 1/2*√(((l-m)*(l+m+1))/(l*(l+1)))
                    norml⁺ = 1/2*√(((l+m)*(l-m+1))/(l*(l+1)))
                    Y = vshbasis(Hansen(), Polar(), l, m, θ, ϕ, S)
                    @test isapproxdefault(Y[1, +1],  0)
                    @test isapproxdefault(Y[2, +1],  (m < l ? norml⁻*cis(-ϕ)*Ylm[(l, m+1)] : 0) -
                                          (m > -l ? norml⁺*cis(ϕ)*Ylm[(l, m-1)] : 0)
                        )
                    if sin(θ) != 0
                        @test isapproxdefault(Y[3, +1],  im*m/√(l*(l+1))*1/sin(θ)*Ylm[(l,m)])
                    end
                    @test isapproxdefault(Y[1,  0],  0)
                    if sin(θ) != 0
                        @test isapproxdefault(Y[2,  0],  -m/√(l*(l+1))*1/sin(θ)*Ylm[(l,m)])
                    end
                    @test isapproxdefault(Y[3, 0],  (m < l ? -im*norml⁻*cis(-ϕ)*Ylm[(l, m+1)] : 0) +
                                          (m > -l ? im*norml⁺*cis(ϕ)*Ylm[(l, m-1)] : 0)
                        )
                    @test isapproxdefault(Y[1, -1],  Ylm[(l,m)])
                    @test isapproxdefault(Y[2, -1],  0)
                    @test isapproxdefault(Y[3, -1],  0)

                    if m != 0 && θ == 0
                        # n.Y⁽¹⁾ and n.Y⁽⁰⁾ are zero, and n.Y⁽⁻¹⁾ == YLM which is zero at poles for M != 0
                        @test all(x -> isapproxdefault(x, 0), Y[1, :])
                    end
                end
            end
        end
        @testset "Irreducible" begin
            @testset "Pole" begin
                θ = 0; ϕ = 0
                cache!(S, θ, ϕ)
                for l in 1:lmax, m in -l:l
                    Y = vshbasis(Irreducible(), Polar(), l, m, θ, ϕ, S)
                    # n.Yᴶ⁺¹ and n.Yᴶ⁻¹ ∝ YLM which is zero at poles for M != 0. n.Yᴶ is zero
                    m == 0 && continue
                    @test all(x -> isapproxdefault(x, 0), Y[1, :])
                end
            end
            @testset "J = 0, M = 0" begin
                for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
                    cache!(S, θ, ϕ)
                    Y = vshbasis(Irreducible(), Polar(), 0, 0, θ, ϕ, S)
                    @test all(iszero, Y[:,-1:0])
                    @test Y[1, 1] ≈ -1/√(4pi)
                    @test begin
                        res = all(x -> isapproxdefault(0, x), Y[2:3, 1])
                        if !res
                            @show θ, ϕ
                        end
                        res
                    end
                end
            end
        end
    end
    @testset "SphericalCovariant" begin
        @testset "Irreducible" begin
            # For l = 0 only Y₀₀¹ is non-zero
            for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
                cache!(S, θ, ϕ, 0)
                Y = vshbasis(Irreducible(), SphericalCovariant(), 0, 0, θ, ϕ, S)
                @test all(iszero, Y[:, -1:0])
                @test isapproxdefault(Y[1, 1], 1/√3 * Ylm[(1,-1)])
                @test isapproxdefault(Y[0, 1], -1/√3 * Ylm[(1, 0)])
                @test isapproxdefault(Y[-1, 1], 1/√3 * Ylm[(1, 1)])
            end
            for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
                cache!(S, θ, ϕ)
                for j in 1:lmax, m in -j:j
                    Y = vshbasis(Irreducible(), SphericalCovariant(), j, m, θ, ϕ, S)
                    @test isapproxdefault(Y[ 1, j+1],  √((j-m+1)*(j-m+2)/(2*(j+1)*(2j+3))) * Ylm[(j+1,m-1)])
                    @test isapproxdefault(Y[ 0, j+1], -√((j-m+1)*(j+m+1)/((j+1)*(2j+3)))   * Ylm[(j+1,m)])
                    @test isapproxdefault(Y[-1, j+1],  √((j+m+1)*(j+m+2)/(2*(j+1)*(2j+3))) * Ylm[(j+1,m+1)])
                    @test isapproxdefault(Y[ 1, j], m > -j ? -√((j+m)*(j-m+1)/(2j*(j+1))) * Ylm[(j,m-1)] : 0)
                    @test isapproxdefault(Y[ 0, j],  m/√(j*(j+1)) * Ylm[(j,m)])
                    @test isapproxdefault(Y[-1, j], m < j ? √((j-m)*(j+m+1)/(2j*(j+1))) * Ylm[(j,m+1)] : 0)
                    @test isapproxdefault(Y[ 1, j-1], abs(m - 1) <= j-1 ? √((j+m)*(j+m-1)/(2j*(2j-1))) * Ylm[(j-1,m-1)] : 0)
                    @test isapproxdefault(Y[ 0, j-1], abs(m) <= j-1 ? √((j-m)*(j+m)/(j*(2j-1))) * Ylm[(j-1,m)] : 0)
                    @test isapproxdefault(Y[-1, j-1], abs(m+1) <= j-1 ? √((j-m)*(j-m-1)/(2j*(2j-1))) * Ylm[(j-1,m+1)] : 0)
                end
            end
            @testset "Pole" begin
                @testset "NorthPole" begin
                    cache!(S, NorthPole(), 0, lmax)
                    for j in 1:lmax, m in -j:j
                        # test that θ = 0 and θ = NorthPole() lead to identical results
                        Y0 = vshbasis(Irreducible(), SphericalCovariant(), j, m, 0, 0)
                        Y = vshbasis(Irreducible(), SphericalCovariant(), j, m, NorthPole(), 0, S)
                        @test isapproxdefault(Y0, Y)
                        if abs(m) > 1
                            @test all(iszero, Y)
                        elseif m == 1 || m == -1
                            @test Y[m,j+1] ≈ √(j/8pi)
                            @test Y[m,j] ≈ -m*√((2j+1)/8pi)
                            @test Y[m,j-1] ≈ √((j+1)/8pi)
                        elseif m == 0
                            @test Y[0,j+1] ≈ -√((j+1)/4pi)
                            @test isapproxdefault(Y[0,j], 0)
                            @test Y[0,j-1] ≈ √(j/4pi)
                        end

                        Y = vshbasis(Hansen(), SphericalCovariant(), j, m, NorthPole(), 0, S)
                        if abs(m) > 1
                            @test all(iszero, Y)
                        elseif m == 1 || m == -1
                            @test Y[m, 1] ≈ √((2j+1)/8pi)
                            @test Y[m, 0] ≈ -m*√((2j+1)/8pi)
                            @test isapproxdefault(Y[m,-1], 0)
                        elseif m == 0
                            @test isapproxdefault(Y[0, 1], 0)
                            @test isapproxdefault(Y[0, 0], 0)
                            @test Y[0,-1] ≈ √((2j+1)/4pi)
                        end
                    end
                end
                @testset "SouthPole" begin
                    cache!(S, SouthPole(), 0, lmax)
                    for j in 1:lmax, m in -j:j
                        # test that θ = pi and θ = SouthPole() lead to identical results
                        Y0 = vshbasis(Irreducible(), SphericalCovariant(), j, m, pi, 0)
                        Y = vshbasis(Irreducible(), SphericalCovariant(), j, m, SouthPole(), 0, S)
                        @test isapproxdefault(Y0, Y)
                        if abs(m) > 1
                            @test all(iszero, Y)
                        elseif m == 1 || m == -1
                            @test Y[m,j+1] ≈ (-1)^(j-1)*√(j/8pi)
                            @test Y[m,j] ≈ -(-1)^(j)*m*√((2j+1)/8pi)
                            @test Y[m,j-1] ≈ (-1)^(j-1)*√((j+1)/8pi)
                        elseif m == 0
                            @test Y[0,j+1] ≈ -(-1)^(j-1)*√((j+1)/4pi)
                            @test isapproxdefault(Y[0,j], 0)
                            @test Y[0,j-1] ≈ (-1)^(j-1)*√(j/4pi)
                        end

                        Y = vshbasis(Hansen(), SphericalCovariant(), j, m, SouthPole(), 0, S)
                        if abs(m) > 1
                            @test all(iszero, Y)
                        elseif m == 1 || m == -1
                            @test Y[m, 1] ≈ (-1)^(j-1)*√((2j+1)/8pi)
                            @test Y[m, 0] ≈ -(-1)^(j)*m*√((2j+1)/8pi)
                            @test isapproxdefault(Y[m,-1], 0)
                        elseif m == 0
                            @test isapproxdefault(Y[0, 1], 0)
                            @test isapproxdefault(Y[0, 0], 0)
                            @test Y[0,-1] ≈ (-1)^(j-1)*√((2j+1)/4pi)
                        end
                    end
                end
            end
        end
        @testset "Hansen" begin
            for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
                cache!(S, θ, ϕ, 0)
                Y = vshbasis(Hansen(), SphericalCovariant(), 0, 0, θ, ϕ, S)
                @test all(iszero, Y[:, 0:1])
                @test isapproxdefault(Y[ 1, -1], -1/√3 * Ylm[(1,-1)])
                @test isapproxdefault(Y[ 0, -1], 1/√3 * Ylm[(1, 0)])
                @test isapproxdefault(Y[-1, -1], -1/√3 * Ylm[(1, 1)])
            end
            for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
                cache!(S, θ, ϕ)
                for l in 1:lmax, m in -l:l
                    Y = vshbasis(Hansen(), SphericalCovariant(), l, m, θ, ϕ, S)
                    @test isapproxdefault(Y[ 1, 1],  √((l-m+1)*(l-m+2)*l/(2*(l+1)*(2l+1)*(2l+3))) * Ylm[(l+1,m-1)] +
                                                    (abs(m - 1) <= l-1 ? √((l+m)*(l+m-1)*(l+1)/(2l*(2l-1)*(2l+1))) * Ylm[(l-1,m-1)] : 0))
                    @test isapproxdefault(Y[ 0, 1], -√((l-m+1)*(l+m+1)*l/((l+1)*(2l+1)*(2l+3)))   * Ylm[(l+1,m)] +
                                                    (abs(m) <= l-1 ? √((l-m)*(l+m)*(l+1)/(l*(2l-1)*(2l+1))) * Ylm[(l-1,m)] : 0))
                    @test isapproxdefault(Y[-1, 1],  √((l+m+1)*(l+m+2)*l/(2*(l+1)*(2l+1)*(2l+3))) * Ylm[(l+1,m+1)] +
                                                    (abs(m+1) <= l-1 ? √((l-m)*(l-m-1)*(l+1)/(2l*(2l-1)*(2l+1))) * Ylm[(l-1,m+1)] : 0))
                    @test isapproxdefault(Y[ 1, 0], m > -l ? -√((l+m)*(l-m+1)/(2l*(l+1))) * Ylm[(l,m-1)] : 0)
                    @test isapproxdefault(Y[ 0, 0],  m/√(l*(l+1)) * Ylm[(l,m)])
                    @test isapproxdefault(Y[-1, 0], m < l ? √((l-m)*(l+m+1)/(2l*(l+1))) * Ylm[(l,m+1)] : 0)
                    @test isapproxdefault(Y[ 1,-1], (abs(m - 1) <= l-1 ? √((l+m)*(l+m-1)/(2*(2l-1)*(2l+1))) * Ylm[(l-1,m-1)] : 0) -
                                                    √((l-m+1)*(l-m+2)/(2*(2l+1)*(2l+3))) * Ylm[(l+1,m-1)])
                    @test isapproxdefault(Y[ 0,-1], (abs(m) <= l-1 ? √((l-m)*(l+m)/((2l-1)*(2l+1))) * Ylm[(l-1,m)] : 0) +
                                                    √((l-m+1)*(l+m+1)/((2l+1)*(2l+3))) * Ylm[(l+1,m)])
                    @test isapproxdefault(Y[-1,-1], (abs(m+1) <= l-1 ? √((l-m)*(l-m-1)/(2*(2l-1)*(2l+1))) * Ylm[(l-1,m+1)] : 0) -
                                                    √((l+m+1)*(l+m+2)/(2*(2l+1)*(2l+3))) * Ylm[(l+1,m+1)])
                end
            end
        end
    end
    @testset "matrix and vectors" begin
        for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
            cache!(S, θ, ϕ)
            for YT in [PB(), Irreducible(), Hansen()],
                B in [Cartesian(), Polar(), SphericalCovariant(), HelicityCovariant()]

                for l in 0:lmax, m in -l:l
                    Y = vshbasis(YT, B, l, m, θ, ϕ, S)
                    Yp = parent(Y)
                    for n in axes(Y, 2)
                        Yn = vshbasis(YT, B, l, m, n, θ, ϕ, S)
                        @test begin
                            res = parent(Y[:, n]) ≈ parent(Yn)
                            if !res
                                @show YT, B, l, m, n
                            end
                            res
                        end
                    end
                end
            end
        end

        @testset "genspharm" begin
            for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
                cache!(S, θ, ϕ)
                G2 = genspharm(ML(0:lmax), θ, ϕ, S)
                for l in 0:lmax, m in -l:l
                    Y = vshbasis(PB(), HelicityCovariant(), l, m, θ, ϕ, S)
                    G = genspharm(l, m, θ, ϕ, S)
                    @test isapproxdefault(G2[(l,m)], G)
                    nmax = min(1, l)
                    for n in -nmax:nmax
                        Ylmn = genspharm2(l, m, n, θ, ϕ)
                        @test isapproxdefault(Y[n,n], G[n])
                        @test isapproxdefault(G[n], Ylmn)
                    end
                end
            end
        end
    end
end

@testset "one vs multiple" begin
    lmax = 2
    modes = LM(0:lmax)
    modes_high = LM(2:lmax, 2:2)
    S2 = VectorSphericalHarmonics.cache(Float64, lmax)
    for YT in [Irreducible(), PB(), Hansen()], B in [Cartesian(), Polar(), SphericalCovariant(), HelicityCovariant()]
        @testset "$YT $B" begin
            for m in [modes, modes_high], (θ, ϕ) in ((pi/3, pi/3), (0,0))
                cache!(S2, θ, ϕ)
                Y_multiple = vshbasis(YT, B, m, θ, ϕ, S2)
                Y_multiple2 = vshbasis(YT, B, m, θ, ϕ)
                for (ind, (j, m)) in enumerate(m)
                    Y = vshbasis(YT, B, j, m, θ, ϕ, S2)
                    @test isapproxdefault(Y_multiple[ind], Y)
                    @test isapproxdefault(Y_multiple2[ind], Y)
                end
            end
        end
    end

    @testset "genspharm" begin
        for m in [modes, modes_high], (θ, ϕ) in ((pi/3, pi/3), (0,0))
            cache!(S2, θ, ϕ)
            Y_multiple = genspharm(m, θ, ϕ, S2)
            Y_multiple2 = genspharm(m, θ, ϕ)
            for (ind, (j, m)) in enumerate(m)
                Y = genspharm(j, m, θ, ϕ, S2)
                @test isapproxdefault(Y_multiple[ind], Y)
                @test isapproxdefault(Y_multiple2[ind], Y)
            end
        end
    end
end

@testset "Orthogonal" begin
    lmax = 10
    @testset "Locally orthogonal" begin
        @testset "PB" begin
            # The PB basis vectors are directed along the Helicity basis vectors, so they are orthogonal at each point
            for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
                cache!(S, θ, ϕ);
                for l in 0:lmax, m in -l:l, B in [SphericalCovariant(), Polar(), HelicityCovariant(), Cartesian()]
                    Y = vshbasis(PB(), B, l, m, θ, ϕ, S)
                    Yp = parent(Y)
                    M = Yp'*Yp
                    for j in axes(M, 2), i in axes(M, 1)
                        i == j && continue
                        @test isapprox(M[i, j], 0, atol=1e-13, rtol=1e-8)
                    end
                end
            end
        end
        @testset "Hansen" begin
            # For Hansen basis, the local orthogonality uses transpose instead of adjoint.
            for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
                cache!(S, θ, ϕ);
                for l in 0:lmax, m in -l:l
                    for B in [Polar(), Cartesian()]
                        Y = vshbasis(Hansen(), B, l, m, θ, ϕ, S)
                        Yp = parent(Y)
                        M = transpose(Yp)*Yp
                        for j in axes(M, 2), i in axes(M, 1)
                            i == j && continue
                            @test begin
                                res = isapprox(M[i, j], 0, atol=1e-13, rtol=1e-8)
                                if !res
                                    @show B, l, m, θ, ϕ
                                end
                                res
                            end
                        end
                    end
                    # For the complex bases, we need to use eμ⋅eν = (-1)^μ δμ,-ν
                    for B in [SphericalCovariant(), HelicityCovariant()]
                        Y = vshbasis(Hansen(), B, l, m, θ, ϕ, S)
                        for j in axes(Y, 2), i in axes(Y, 2)
                            i == j && continue
                            M = sum(Y[b,j] * Y[-b,i] * (-1)^b for b = -1:1)
                            @test begin
                                res = isapprox(M, 0, atol=1e-13, rtol=1e-8)
                                if !res
                                    @show B, l, m, θ, ϕ
                                end
                                res
                            end
                        end
                    end
                end
            end
        end
    end
    @testset "Orthonormal" begin
        lmax = 3
        function f(θ, ϕ, YT, B, l, m)
            cache!(S, θ, ϕ)
            Y = vshbasis(YT, B, l, m, θ, ϕ, S)
            Yp = parent(Y)
            (Yp'*Yp) .* sin(θ)
        end
        @testset "PB" begin
            for l in 0:lmax, m in -l:l, B in [SphericalCovariant(), Polar(), HelicityCovariant(), Cartesian()]
                M = hcubature(x -> f(x..., PB(), B, l, m), [0, 0], [pi, 2pi])[1]
                if l == 0
                    for I in CartesianIndices(M)
                        i, j = Tuple(I)
                        i == 2 && j == 2 && continue
                        @test isapproxdefault(M[i,j], 0)
                    end
                    @test M[2,2] ≈ 1
                else
                    @test M ≈ I
                end
            end
        end
        @testset "Irreducible" begin
            for l in 0:lmax, m in -l:l, B in [SphericalCovariant(), Polar(), HelicityCovariant(), Cartesian()]
                M = hcubature(x -> f(x..., Irreducible(), B, l, m), [0, 0], [pi, 2pi])[1]
                if l == 0
                    for I in CartesianIndices(M)
                        i, j = Tuple(I)
                        i == 3 && j == 3 && continue
                        @test isapproxdefault(M[i,j], 0)
                    end
                    @test M[3,3] ≈ 1
                else
                    @test M ≈ I
                end
            end
        end
        @testset "Hansen" begin
            for l in 0:lmax, m in -l:l, B in [SphericalCovariant(), Polar(), HelicityCovariant(), Cartesian()]
                M = hcubature(x -> f(x..., Hansen(), B, l, m), [0, 0], [pi, 2pi])[1]
                if l == 0
                    for I in CartesianIndices(M)
                        i, j = Tuple(I)
                        i == 1 && j == 1 && continue
                        @test isapproxdefault(M[i,j], 0)
                    end
                    @test M[1,1] ≈ 1
                else
                    @test M ≈ I
                end
            end
        end
    end
end

@testset "Addition theorem" begin
    lmax = 10
    @testset "same point" begin
        @testset "Irreducible" begin
            θ, ϕ = pi/3, pi/3
            for B in [Polar(), SphericalCovariant(), Cartesian(), HelicityCovariant()]
                Y1_all = vshbasis(Irreducible(), B, ML(0:lmax), θ, ϕ);
                for j in 1:lmax, l1 in j-1:j+1, l2 in j-1:j+1
                    S = zero(Float64)
                    for m in -j:j
                        Yjm = Y1_all[(j,m)]
                        Yjml1 = parent(Yjm[:,l1])
                        Yjml2 = parent(Yjm[:,l2])
                        S += dot(Yjml1, Yjml2)
                    end
                    if l1 == l2
                        @test isapproxdefault(S, (2j+1)/4pi)
                    else
                        @test isapproxdefault(S, 0)
                    end
                end
            end
        end
        @testset "Hansen" begin
            θ, ϕ = pi/3, pi/3
            for B in [Polar(), SphericalCovariant(), Cartesian(), HelicityCovariant()]
                Y1_all = vshbasis(Hansen(), B, ML(0:lmax), θ, ϕ);
                for j in 1:lmax, λ1 in -1:+1, λ2 in -1:+1
                    S = zero(Float64)
                    for m in -j:j
                        Yjm = Y1_all[(j,m)]
                        Yjmλ1 = parent(Yjm[:,λ1])
                        Yjmλ2 = parent(Yjm[:,λ2])
                        S += dot(Yjmλ1, Yjmλ2)
                    end
                    if λ1 == λ2
                        @test isapproxdefault(S, (2j+1)/4pi)
                    else
                        @test isapproxdefault(S, 0)
                    end
                end
            end
        end
    end
    @testset "different points" begin
        @testset "Irreducible" begin
            θ1, ϕ1 = pi/3, pi/3
            θ2, ϕ2 = pi/4, pi/4
            cosω12 = cos(θ1)cos(θ2) + sin(θ1)sin(θ2)cos(ϕ1 - ϕ2)
            P = LegendrePolynomials.collectPl(cosω12, lmax = lmax + 1);
            allmodes = ML(0:lmax)
            # We only test in bases that are constant, so they are identical at both the points
            for B in [SphericalCovariant(), Cartesian()]
                Y1_all = vshbasis(Irreducible(), B, allmodes, θ1, ϕ1);
                Y2_all = vshbasis(Irreducible(), B, allmodes, θ2, ϕ2);
                for j in 1:lmax, l1 in j-1:j+1, l2 in j-1:j+1
                    x = zero(Float64)
                    for m in -j:j
                        Y1jm = Y1_all[(j,m)]
                        Y2jm = Y2_all[(j,m)]
                        Y1jml1 = parent(Y1jm[:,l1])
                        Y2jml2 = parent(Y2jm[:,l2])
                        x += dot(Y1jml1, Y2jml2)
                    end
                    @test isapproxdefault(imag(x), 0)
                    if l1 == l2
                        @test begin
                            res = isapproxdefault(x, (2j+1)/4pi * P[l1])
                            if !res
                                @show j, l1, B, x, (2j+1)/4pi * P[l1]
                            end
                            res
                        end
                    else
                        @test begin
                            res = isapproxdefault(x, 0)
                            if !res
                                @show j, l1, l2, B, x
                            end
                            res
                        end
                    end
                end
            end
        end
    end
end

@testset "Complex conjugation" begin
    @testset "Irreducible" begin
        # For l = 0 Y₀₀¹ is purely real
        for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
            cache!(S, θ, ϕ, 0)
            Y = vshbasis(Irreducible(), Polar(), 0, 0, θ, ϕ, S)
            @test all(x -> isapproxdefault(imag(x), 0), Y)
        end
        for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
            cache!(S, θ, ϕ)
            for j in 1:5, m in -j:j
                Ym = vshbasis(Irreducible(), Polar(), j, m, θ, ϕ, S)
                Ymconj = conj.(Ym)
                Y₋m = vshbasis(Irreducible(), Polar(), j, -m, θ, ϕ, S)
                for l in axes(Ym, 2)
                    Ymnconj = Ymconj[:, l]
                    Y₋mn = Y₋m[:, l]
                    phase = (-1)^(j + m + l + 1)
                    @test all(isapproxdefault, zip(Ymnconj, phase .* Y₋mn))
                end
            end
        end
    end
    @testset "Hansen" begin
        # For l = 0 Y₀₀⁽⁻¹⁾ is purely real
        for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
            cache!(S, θ, ϕ, 0)
            Y = vshbasis(Hansen(), Polar(), 0, 0, θ, ϕ, S)
            @test all(x -> isapproxdefault(imag(x), 0), Y)
        end
        for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
            cache!(S, θ, ϕ)
            for l in 1:5, m in -l:l
                Ym = vshbasis(Hansen(), Polar(), l, m, θ, ϕ, S)
                Ymconj = conj.(Ym)
                Y₋m = vshbasis(Hansen(), Polar(), l, -m, θ, ϕ, S)
                for λ in axes(Ym, 2)
                    Ymnconj = Ymconj[:, λ]
                    Y₋mn = Y₋m[:, λ]
                    phase = (-1)^(m + λ + 1)
                    @test all(isapproxdefault.(Ymnconj, phase .* Y₋mn))
                end
            end
        end
    end
    @testset "PB" begin
        for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
            cache!(S, θ, ϕ)
            for j in 1:5, m in -j:j
                Ym = vshbasis(PB(), Polar(), j, m, θ, ϕ, S)
                Ymconj = conj.(Ym)
                Y₋m = vshbasis(PB(), Polar(), j, -m, θ, ϕ, S)
                for n in axes(Ym, 2)
                    Y₋m₋n = Y₋m[:, -n]
                    Ymn = Ym[:, n]
                    phase = (-1)^m
                    @test begin
                        res = all(isapproxdefault, zip(conj.(Y₋m₋n), phase .* Ymn))
                        if !res
                            @show j, m, n, θ, ϕ
                        end
                        res
                    end
                end
            end
        end
        for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
            cache!(S, θ, ϕ)
            for j in 1:5, m in -j:j
                Ym = vshbasis(PB(), HelicityCovariant(), j, m, θ, ϕ, S)
                Ymconj = conj.(Ym)
                Y₋m = vshbasis(PB(), HelicityCovariant(), j, -m, θ, ϕ, S)
                for n in axes(Ym, 2)
                    Y₋m₋n = Y₋m[-n, -n]
                    Ymn = Ym[n, n]
                    phase = (-1)^(m + n)
                    @test begin
                        res = all(isapproxdefault, zip(conj.(Y₋m₋n), phase .* Ymn))
                        if !res
                            @show j, m, n, θ, ϕ
                        end
                        res
                    end
                end
            end
        end
    end
end

@testset "Coordinate inversion" begin
    jmax = 4
    @testset "Irreducible" begin
        for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
            cache!(S, θ, ϕ, jmax)
            for j in 0:jmax, m in -j:j
                for B in [Cartesian(), SphericalCovariant()]
                    Yjm = vshbasis(Irreducible(), B, j, m, θ, ϕ, S)
                    Yjmrot = vshbasis(Irreducible(), B, j, m, pi - θ, pi + ϕ)
                    for l in axes(Yjm, 2)
                        Yjml = Yjm[:, l] * (-1)^(l+1)
                        Yjmlrot = -Yjmrot[:, l]
                        @test begin
                            res = all(isapproxdefault, zip(Yjml, Yjmlrot))
                            if !res
                                @show j, m, θ, ϕ
                            end
                            res
                        end
                    end
                end
            end
        end
    end
    @testset "Hansen" begin
        for θ in LinRange(0, pi, 10), ϕ in LinRange(0, 2pi, 10)
            cache!(S, θ, ϕ, jmax)
            for j in 0:4, m in -j:j
                for B in [Cartesian(), SphericalCovariant()]
                    Yjm = vshbasis(Hansen(), B, j, m, θ, ϕ, S)
                    Yjmrot = vshbasis(Hansen(), B, j, m, pi - θ, pi + ϕ)
                    for λ in axes(Yjm, 2)
                        Yjml = Yjm[:, λ] * (-1)^(j + λ +1)
                        Yjmlrot = -Yjmrot[:, λ]
                        @test all(isapproxdefault, zip(Yjml, Yjmlrot))
                    end
                end
            end
        end
    end
end

@testset "Coordinate rotation" begin
    lmax = 3
    modes = ML(0:lmax)

    S2 = VectorSphericalHarmonics.cache(0, 0, lmax);
    S′ = VectorSphericalHarmonics.cache(0, 0, lmax);
    Dvec = OffsetArray([zeros(ComplexF64, 2l+1, 2l+1) for l in 0:lmax], 0:lmax);

    @testset "arbitrary point" begin
        cartvec(θ, ϕ) = SVector{3}(sin(θ)cos(ϕ), sin(θ)sin(ϕ), cos(θ))
        function polcoords(n)
            nx, ny, nz = normalize(n)
            θ = acos(nz)
            if abs(nz) == 1
                # ϕ is not defined uniquely at the poles
                return promote(θ, zero(θ))
            else
                ϕ = atan(ny, nx)
                return promote(θ, ϕ)
            end
        end

        @testset "genspharm" begin
            for θ in LinRange(0, pi, 6), ϕ in LinRange(0, 2pi, 6)
                VectorSphericalHarmonics.cache!(S2, θ, ϕ);
                Y = genspharm(modes, θ, ϕ, S2);
                n = cartvec(θ, ϕ)
                Un = basisconversionmatrix(Cartesian(), HelicityCovariant(), θ, ϕ);

                for θ′ in LinRange(0, pi, 6), ϕ′ in LinRange(0, 2pi, 6)
                    VectorSphericalHarmonics.cache!(S′, θ′, ϕ′);
                    α, β, γ = ϕ, θ-θ′, -ϕ′
                    R = RotZYZ(α, β, γ)
                    R⁻¹ = inv(R)
                    # this specific rotation satisfies the condition U′n * R' * Un' == I,
                    # so the basis rotation matrix may be left out
                    # This happens because Un for polar == RotYZ(-θ, -ϕ) (albeit with permuted rows), so
                    # U′n * R' * Un' = RotYZ(-θ′, -ϕ′) * RotZYZ(ϕ′, θ′ - θ, -ϕ) * RotZY(ϕ, θ)  ≈ I
                    # Un for HelicityCovariant is related to that for the polar matrix through a unitary
                    # transformation, so the same relation holds
                    U′n = basisconversionmatrix(Cartesian(), HelicityCovariant(), θ′, ϕ′);
                    @test U′n * R⁻¹ * Un' ≈ I
                    for l in 0:lmax
                        Dp = wignerD!(Dvec[l], l, α, β, γ);
                        D = OffsetArray(Dp, -l:l, -l:l)
                        for m in -l:l
                            Ylθϕrot_m = sum(D[m′, m] * Y[(l, m′)] for m′ in -l:l)
                            Ylmθ′ϕ′ = genspharm(l, m, θ′, ϕ′, S′)
                            @test isapprox(Ylmθ′ϕ′, Ylθϕrot_m, atol = 1e-13, rtol = 1e-8)
                        end
                    end
                end

                # obtain the second point from an arbitrary rotation
                # In general we need to retain the basis conversion matrix to compute the components correctly
                for α in LinRange(0, 2pi, 6), β in LinRange(0, pi, 6), γ in LinRange(0, 2pi, 6)
                    # the rotation that transforms between frames S′ = RS
                    R = RotZYZ(α, β, γ)
                    R⁻¹ = inv(R)
                    n′ = R⁻¹ * n
                    θ′, ϕ′ = polcoords(n′)
                    VectorSphericalHarmonics.cache!(S′, θ′, ϕ′);
                    U′n = basisconversionmatrix(Cartesian(), HelicityCovariant(), θ′, ϕ′);
                    M⁻¹ = U′n * R⁻¹ * Un'
                    @test M⁻¹ ≈ Diagonal(M⁻¹)
                    for l in 0:lmax
                        Dp = wignerD!(Dvec[l], l, α, β, γ);
                        D = OffsetArray(Dp, -l:l, -l:l)
                        for m in -l:l
                            Ylθϕrot_m = M⁻¹ * sum(D[m′, m] * parent(Y[(l, m′)]) for m′ in -l:l)
                            Ylmθ′ϕ′ = parent(genspharm(l, m, θ′, ϕ′, S′))
                            @test isapprox(Ylmθ′ϕ′, Ylθϕrot_m, atol = 1e-13, rtol = 1e-8)
                        end
                    end
                end
            end
        end

        @testset "VSH" begin
            for θ in LinRange(0, pi, 6), ϕ in LinRange(0, 2pi, 6)
                n = cartvec(θ, ϕ)
                VectorSphericalHarmonics.cache!(S2, θ, ϕ);

                for YT in [Irreducible(), Hansen(), PB()], B in [Cartesian(), SphericalCovariant(), HelicityCovariant(), Polar()]
                    Un = basisconversionmatrix(Cartesian(), B, θ, ϕ)
                    Y = vshbasis(YT, B, modes, θ, ϕ, S2);

                    # We evaluate the rotation of the dot product between VSH and a vector field
                    # In this case we choose a constant vector field directed along x̂ at each point,
                    # where x̂ is directed along the x axis of the frame S (which differs from that of S′ in general)
                    # We evaluate the components in the basis B at both points
                    x_n = Un * SVector{3}(1,0,0)

                    if B ∈ (Polar(), HelicityCovariant())
                        for θ′ in LinRange(0, pi, 6), ϕ′ in LinRange(0, 2pi, 6)
                            VectorSphericalHarmonics.cache!(S′, θ′, ϕ′);
                            α, β, γ = ϕ, θ-θ′, -ϕ′
                            R = RotZYZ(α, β, γ)
                            # this specific rotation satisfies the condition U′n * R' * Un' == I for the
                            # Polar and HelicityCovariant bases, so the basis rotation matrix may be left out
                            # This happens because Un for polar == RotYZ(-θ, -ϕ) (albeit with permuted rows), so
                            # U′n * R' * Un' = RotYZ(-θ′, -ϕ′) * RotZYZ(ϕ′, θ′ - θ, -ϕ) * RotZY(ϕ, θ) ≈ I
                            # Un for HelicityCovariant is related to that for the polar matrix through a unitary
                            # transformation, so the same relation holds
                            U′n = basisconversionmatrix(Cartesian(), B, θ′, ϕ′);
                            if B === Polar()
                                # We go from (r,θ,ϕ) to (θ,ϕ,r)
                                @test Un[[2,3,1],:] ≈ RotYZ(-θ, -ϕ)
                                @test U′n[[2,3,1],:] ≈ RotYZ(-θ′, -ϕ′)
                            end
                            @test U′n * R' * Un' ≈ I
                            for l in 0:lmax
                                Dp = wignerD!(Dvec[l], l, α, β, γ);
                                D = OffsetArray(Dp, -l:l, -l:l)
                                for m in -l:l
                                    Ylθϕrot_m = sum(D[m′, m] * Y[(l, m′)] for m′ in -l:l)
                                    Ylmθ′ϕ′ = vshbasis(YT, B, l, m, θ′, ϕ′, S′)
                                    @test isapprox(Ylmθ′ϕ′, Ylθϕrot_m, atol = 1e-13, rtol = 1e-8)
                                end
                            end
                        end
                    end

                    # obtain the second point from an arbitrary rotation
                    # In general we need to retain the basis conversion matrix to compute the components correctly
                    for α in LinRange(0, 2pi, 6), β in LinRange(0, pi, 6), γ in LinRange(0, 2pi, 6)
                        # the rotation that transforms between frames S′ = RS
                        R = RotZYZ(α, β, γ)
                        R⁻¹ = inv(R)
                        n′ = R⁻¹ * n
                        θ′, ϕ′ = polcoords(n′)
                        VectorSphericalHarmonics.cache!(S′, θ′, ϕ′);
                        U′n = basisconversionmatrix(Cartesian(), B, θ′, ϕ′);

                        M⁻¹ = (U′n * R⁻¹ * Un')
                        if B === HelicityCovariant()
                            @test M⁻¹ ≈ Diagonal(M⁻¹)
                        end
                        M = inv(M⁻¹)

                        x_n′ = M⁻¹ * x_n

                        for l in 0:lmax
                            Dp = wignerD!(Dvec[l], l, α, β, γ);
                            D = OffsetArray(Dp, -l:l, -l:l)
                            for m in -l:l
                                Ylθϕrot_m = M⁻¹ * sum(D[m′, m] * parent(Y[(l, m′)]) for m′ in -l:l)
                                Ylmθ′ϕ′ = parent(vshbasis(YT, B, l, m, θ′, ϕ′, S′))
                                @test isapprox(Ylmθ′ϕ′, Ylθϕrot_m, atol = 1e-13, rtol = 1e-8)
                                # (M⁻¹ x)⋅Ylm′(n′) = (M⁻¹ x)⋅M⁻¹∑Dlmm′ Ylm(n) = x⋅∑Dlmm′ Ylm(n)
                                Ylθϕrot_m = x_n' * sum(D[m′, m] * parent(Y[(l, m′)]) for m′ in -l:l)
                                Ylmθ′ϕ′ = x_n′' * parent(vshbasis(YT, B, l, m, θ′, ϕ′, S′))
                                @test isapprox(Ylmθ′ϕ′, Ylθϕrot_m, atol = 1e-13, rtol = 1e-8)
                                # x.Ylm′(n′) = x.M⁻¹∑Dlmm′ Ylm(n) = (Mx)⋅∑Dlmm′ Ylm(n)
                                Ylθϕrot_m = (M*x_n)' * sum(D[m′, m] * parent(Y[(l, m′)]) for m′ in -l:l)
                                Ylmθ′ϕ′ = x_n' * parent(vshbasis(YT, B, l, m, θ′, ϕ′, S′))
                                @test isapprox(Ylmθ′ϕ′, Ylθϕrot_m, atol = 1e-13, rtol = 1e-8)
                            end
                        end
                    end
                end
            end
        end
    end
end

@testset "VSHcache" begin
    modes = LM(0:1, 0:1)
    lmax = 2
    θ1, ϕ1 = pi/3, pi/2
    θ2, ϕ2 = pi/4, pi/6
    @testset "vshbasis!" begin
        V = @inferred VectorSphericalHarmonics.VSHCache(Float64, Irreducible(), Polar(), θ1, ϕ1, modes);
        M = vshbasis(Irreducible(), Polar(), modes, θ1, ϕ1)
        @test getY(V) == M
        @test eltypeY(V) == eltype(getY(V))

        for θ in Any[θ2, NorthPole(), SouthPole()]
            vshbasis!(V, Irreducible(), Polar(), θ, ϕ2)
            M = vshbasis(Irreducible(), Polar(), modes, θ, ϕ2)
            @test getY(V) == M
        end
    end
    @testset "genspharm!" begin
        V = @inferred VectorSphericalHarmonics.VSHCache(Float64, θ1, ϕ1, modes);
        M = genspharm(modes, θ1, ϕ1)
        @test getY(V) == M

        for θ in Any[θ2, NorthPole(), SouthPole()]
            genspharm!(V, θ, ϕ2)
            M = genspharm(modes, θ, ϕ2)
            @test getY(V) == M
        end
    end
    @testset "uninitialized" begin
        V = VectorSphericalHarmonics.VSHCache(PB(), Cartesian(), LM(0:1))
        @test eltype(VectorSphericalHarmonics.eltypeY(V)) == ComplexF64
        @test all(iszero, VectorSphericalHarmonics.getY(V))

        V = VectorSphericalHarmonics.VSHCache(LM(0:1))
        @test eltype(VectorSphericalHarmonics.eltypeY(V)) == ComplexF64
        @test all(iszero, VectorSphericalHarmonics.getY(V))
    end
end

@testset "broadcast tags" begin
    function components_of_x̂(Basis, (θ,ϕ))
        B = VectorSphericalHarmonics.basisconversionmatrix(Cartesian(), Basis, θ, ϕ)
        B * SVector{3}(1,0,0)
    end
    for B in [Cartesian(), Polar(), SphericalCovariant(), HelicityCovariant()]
        x1, x2 = components_of_x̂.(B, ((pi/2, 0), (pi/2, pi/4)))
        @test x1 == components_of_x̂(B, (pi/2, 0))
        @test x2 == components_of_x̂(B, (pi/2, pi/4))
    end
    V = vshbasis.(Hansen(), Polar(), 1:2, -1, pi/2, 0)
    @test V[1] == vshbasis(Hansen(), Polar(), 1, -1, pi/2, 0)
    @test V[2] == vshbasis(Hansen(), Polar(), 2, -1, pi/2, 0)
    V = vshbasis.(Hansen(), Polar(), 1, -1, (pi/2, pi/4), 0)
    @test V[1] == vshbasis(Hansen(), Polar(), 1, -1, pi/2, 0)
    @test V[2] == vshbasis(Hansen(), Polar(), 1, -1, pi/4, 0)
end


text:
using EffectiveWaves, Test
using LinearAlgebra

@testset "high frequency effective" begin
        medium = Medium(1.0,1.0+0.0im)
        # Large weak scatterers with low volume fraciton
        species = [
            Specie(ρ=10.,r=1.9, c=12., volfrac=0.04),
            Specie(ρ=3., r=0.7, c=2.0, volfrac=0.02)
        ]
        ωs2 = [120.]

        tol = 1e-6
        k_eff_φs = wavenumber_low_volfrac(ωs2, medium, species; tol=tol)
        k_effs = [wavenumbers(ω, medium, species; tol=tol, num_wavenumbers=1) for ω in ωs2]
        inds = [argmin(abs.(k_effs[i] .- k_eff_φs[i])) for i in eachindex(ωs2)]
        k_effs2 = [k_effs[i][inds[i]] for i in eachindex(inds)]

        @test norm(k_effs2 - k_eff_φs)/norm(k_effs2) < tol

        Rs = map(eachindex(ωs2)) do i
            wave = EffectiveWave(ωs2[i], k_effs2[i], medium, species)
            reflection_coefficient(ωs2[i], wave, medium, species)
        end
        # warning is expected, as k_eff_φs are assymptotic approximations.
        Rs_φs = map(eachindex(ωs2)) do i
            wave = EffectiveWave(ωs2[i], k_eff_φs[i], medium, species)
            reflection_coefficient(ωs2[i], wave, medium, species)
        end
        Rs_φs2 = reflection_coefficient_low_volfrac(ωs2, medium, species)

        # the incident wave has amplitude 1, so this is already a relative difference
        @test maximum(abs.(Rs_φs - Rs)) < tol
        @test maximum(abs.(Rs_φs2 - Rs)) < tol
end

ext: jl


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8,6))
plot(curves.parameter_values, curves.measurements)
xlabel("Number of trees", fontsize=14)
xticks([10, 250, 500, 750, 1000])
ylim([4, 5])

savefig(joinpath(@OUTPUT, "A-ensembles-2-curves.svg")) # hide

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
module TestFunctions

using ModiaMedia
using ModiaMedia.Test  # included via ModiaMedia, to avoid requirement to add it in the standard environment

const mediaNames = ["ConstantPropertyLiquidWater", "SimpleAir", "N2", "MoistAir"]

function testMediumFunctions(mediumName; figure=1)
   Medium = getMedium(mediumName)
   infos  = Medium.infos
   state  = setState_pTX(Medium, infos.reference_p,
                                 infos.reference_T,
                                 infos.reference_X)

   p   = pressure(state)
   T   = temperature(state)
   h   = specificEnthalpy(state)
   d   = density(state)
   u   = specificInternalEnergy(state)
   #Cp  = specificHeatCapacityCp(state)
   #eta = dynamicViscosity(state)

   #println("... ", mediumName, " (figure=",figure,"): p = ", p, ", T = ", T, ", h = ", h, ", d = ", d, ", u = ", u, ", Cp = ", Cp, ", eta = ", eta)
   println("... ", mediumName, " (figure=",figure,"): p = ", p, ", T = ", T, ", h = ", h, ", d = ", d, ", u = ", u)

   state_b = isenthalpicState(state, 0.1e5)
   h_b     = specificEnthalpy(state_b)
   @test isapprox(h,h_b)

   state_b = setState_phX(Medium, p, h, infos.reference_X)
   T_b     = temperature(state_b)
   @test isapprox(T,T_b)

   if !infos.singleState
      state_b = setState_dTX(Medium, d, T, infos.reference_X)
      h_b     = specificEnthalpy(state_b)
      @test isapprox(h,h_b)
   end

   state_c = setState_pTX(Medium, 0.9*infos.reference_p,
                                  0.9*infos.reference_T,
                                  0.9*infos.reference_X)
   state_d = deepcopy(state_c)
   setState_pTX!(state_d, p,T,infos.reference_X)
   @test isapprox(p, pressure(state_d))
   @test isapprox(T, temperature(state_d))

   state_d = deepcopy(state_c)
   setState_phX!(state_d, p,h,infos.reference_X)
   @test isapprox(p, pressure(state_d))
   @test isapprox(h, specificEnthalpy(state_d))

   if !infos.singleState
      state_d = deepcopy(state_c)
      setState_dTX!(state_d, d,T,infos.reference_X)
      @test isapprox(d, density(state_d))
      @test isapprox(T, temperature(state_d))
   end

   if typeof(Medium) <: PureSubstance
      state_b = setState_pT(Medium, p, T)
      T_b     = temperature(state_b)
      @test isapprox(T,T_b)

      state_b = setState_ph(Medium, p, h)
      T_b     = temperature(state_b)
      @test isapprox(T,T_b)

      state_d = deepcopy(state_c)
      setState_pT!(state_d, p,T)
      @test isapprox(p, pressure(state_d))
      @test isapprox(T, temperature(state_d))

      state_d = deepcopy(state_c)
      setState_ph!(state_d, p,h)
      @test isapprox(p, pressure(state_d))
      @test isapprox(h, specificEnthalpy(state_d))

      d2 = density_ph(Medium,p,h)
      T2 = temperature_ph(Medium,p,h)
      @test isapprox(d,d2)
      @test isapprox(T,T2)

      if !infos.singleState
         state_b = setState_dT(Medium, d, T)
         h_b     = specificEnthalpy(state_b)
         @test isapprox(h,h_b)

         state_d = deepcopy(state_c)
         setState_dT!(state_d, d,T)
         @test isapprox(d, density(state_d))
         @test isapprox(T, temperature(state_d))

         p2 = pressure_dT(Medium,d,T)
         h2 = specificEnthalpy_dT(Medium,d,T)
         @test isapprox(p,p2)
         @test isapprox(h,h2)
      end
   end

   ModiaMedia.standardPlot(Medium; figure=figure)
end

function testAllMediumFunctions(mediaNames)
   println("\n... Test thermodynamic property functions:")
   i = 0
   for name in mediaNames
      i += 1
      testMediumFunctions(name; figure=i)
   end
end

@testset "ModiaMedia/test/TestFunctions.jl:" begin
   testAllMediumFunctions(mediaNames)
end

end


text:
module Minimap

using minimap2_jll

import Base: show
import ..Utility: Alignment, Hit

export align

# ------------------------------------------------------------------------
# globals

const NO_DIAG = 0x001
const NO_DUAL = 0x002
const CIGAR   = 0x004

# ------------------------------------------------------------------------
# functions

const idx_str  = (:mm_idx_str, minimap2)
const idx_stat = (:mm_idx_stat, minimap2)
const idx_free = (:mm_idx_destroy, minimap2)

const tbuf_init = (:mm_tbuf_init, minimap2)
const tbuf_free = (:mm_tbuf_destroy, minimap2)

const set_opt    = (:mm_set_opt, minimap2)
const update_opt = (:mm_mapopt_update, minimap2)

const minimap  = (:mm_map, minimap2)
const map_frag = (:mm_map_frag, minimap2)

const divergence = (:mm_event_identity, minimap2)

# ------------------------------------------------------------------------
# types

"""
	struct Index
		handle :: Ptr{Cvoid}
	end

Store the untyped address to a minimap2 sequence index (set of minimizers).
"""
struct Index
    handle :: Ptr{Cvoid}
end

"""
	struct Buffer
		handle :: Ptr{Cvoid}
	end

Store the untyped address to a minimap2 sequence buffer (working space).
"""
struct Buffer
    handle :: Ptr{Cvoid}
end

# copied from minimap2.h
"""
	struct Extra
		capacity :: UInt32
		dp_score :: Int32
		dp_max   :: Int32
		dp_max2  :: Int32
		packed   :: UInt32 # n_ambi (30 bits) / strand (2 bits)
		n_cigar  :: UInt32
		cigar UInt32[] (variable length array)
	end

Copied from minimap2.h. See mm_extra_t.
"""
struct Extra
    capacity :: UInt32
    dp_score :: Int32
    dp_max   :: Int32
    dp_max2  :: Int32
    packed   :: UInt32 # n_ambi (30 bits) / strand (2 bits)
    n_cigar  :: UInt32
    # cigar UInt32[] (variable length array)
end

"""
	struct Record
		id :: Int32; cnt :: Int32; rid :: Int32; score :: Int32
		qs :: Int32; qe  :: Int32; rs  :: Int32; re :: Int32

		parent :: Int32; subsc :: Int32

		as     :: Int32
		mlen   :: Int32; blen :: Int32
		nsub   :: Int32;
		score0 :: Int32;

		pack :: UInt32 # mapq(8) ; split(2); rev(1); inv(1); sam_pri(1); proper_frag(1); pe_thru(1); seg_split(1); seg_id(8); split_inv(1); is_alt(1); dummy(6)
		hash :: UInt32
		div  :: Cfloat
		p    :: Ptr{Extra}
	end

Copied from minimap2.h. See mm_reg1_t.
"""
struct Record
    id :: Int32; cnt :: Int32; rid :: Int32; score :: Int32
    qs :: Int32; qe  :: Int32; rs  :: Int32; re :: Int32

    parent :: Int32; subsc :: Int32

    as     :: Int32
    mlen   :: Int32; blen :: Int32
    nsub   :: Int32;
    score0 :: Int32;

    pack :: UInt32 # mapq(8) ; split(2); rev(1); inv(1); sam_pri(1); proper_frag(1); pe_thru(1); seg_split(1); seg_id(8); split_inv(1); is_alt(1); dummy(6)
    hash :: UInt32
    div  :: Cfloat
    p    :: Ptr{Extra}
end

"""
	struct IndexOptions
		k    :: Cshort
		w    :: Cshort
		flag :: Cshort
		bucket_bits :: Cshort
		mini_batch_size :: Int64
		batch_size :: UInt64

		IndexOptions() = new()
	end

Copied from minimap2.h. See mm_idxopt_t.
"""
struct IndexOptions
    k    :: Cshort
    w    :: Cshort
    flag :: Cshort
    bucket_bits :: Cshort
    mini_batch_size :: Int64
    batch_size :: UInt64

    IndexOptions() = new()
end

"""
	mutable struct MapOptions
		flag :: Int64
		seed :: Cint
		sdust_thres :: Cint

		max_qlen :: Cint

		bw :: Cint; bw_long :: Cint
		max_gap :: Cint; max_gap_ref :: Cint
		max_frag_len :: Cint
		max_chain_skip :: Cint; max_chain_iter :: Cint
		min_cnt :: Cint
		min_chain_score :: Cint
		chain_gap_scale :: Cfloat
		rmq_size_cap :: Cint; rmq_inner_dist :: Cint
		rmq_rescue_size :: Cint
		rmq_rescue_ratio :: Cfloat

		mask_level :: Cfloat
		mask_len :: Cint
		pri_ratio :: Cfloat
		best_n :: Cint

		alt_drop :: Cfloat

		a :: Cint; b :: Cint; q :: Cint; e :: Cint; q2 :: Cint; e2 :: Cint
		sc_ambi :: Cint
		noncan :: Cint
		junc_bonus :: Cint
		zdrop :: Cint; zdrop_inv :: Cint
		end_bonus :: Cint
		min_dp_max :: Cint
		min_ksw_len :: Cint
		anchor_ext_len :: Cint; anchor_ext_shift :: Cint
		max_clip_ratio :: Cfloat

		rank_min_len :: Cint
		rank_frac :: Cfloat

		pe_ori :: Cint; pe_bonus :: Cint

		mid_occ_frac :: Cfloat
		min_mid_occ :: Int32; max_mid_occ :: Int32
		mid_occ :: Int32
		max_occ :: Int32; max_max_occ :: Int32; occ_dist :: Int32;
		mini_batch_size :: Int64
		max_sw_mat :: Int64
		cap_kalloc :: Int64

		split_prefix :: Cstring
	end

Copied from minimap2.h. See mm_mapopt_t.
"""
mutable struct MapOptions
	flag :: Int64
	seed :: Cint
	sdust_thres :: Cint

	max_qlen :: Cint

	bw :: Cint; bw_long :: Cint
	max_gap :: Cint; max_gap_ref :: Cint
	max_frag_len :: Cint
	max_chain_skip :: Cint; max_chain_iter :: Cint
	min_cnt :: Cint
	min_chain_score :: Cint
	chain_gap_scale :: Cfloat
	rmq_size_cap :: Cint; rmq_inner_dist :: Cint
	rmq_rescue_size :: Cint
	rmq_rescue_ratio :: Cfloat

	mask_level :: Cfloat
	mask_len :: Cint
	pri_ratio :: Cfloat
	best_n :: Cint

	alt_drop :: Cfloat

	a :: Cint; b :: Cint; q :: Cint; e :: Cint; q2 :: Cint; e2 :: Cint
	sc_ambi :: Cint
	noncan :: Cint
	junc_bonus :: Cint
	zdrop :: Cint; zdrop_inv :: Cint
	end_bonus :: Cint
	min_dp_max :: Cint
	min_ksw_len :: Cint
	anchor_ext_len :: Cint; anchor_ext_shift :: Cint
	max_clip_ratio :: Cfloat

	rank_min_len :: Cint
	rank_frac :: Cfloat

	pe_ori :: Cint; pe_bonus :: Cint

	mid_occ_frac :: Cfloat
	min_mid_occ :: Int32; max_mid_occ :: Int32
	mid_occ :: Int32
	max_occ :: Int32; max_max_occ :: Int32; occ_dist :: Int32;
	mini_batch_size :: Int64
	max_sw_mat :: Int64
	cap_kalloc :: Int64

	split_prefix :: Cstring

    function MapOptions(idx::Ref{IndexOptions}, minblock::Int, preset::String)
        opt = Ref{MapOptions}(new())

        # TODO: error handling
        err = ccall(set_opt, Cint, (Cstring, Ptr{Cvoid}, Ptr{Cvoid}),
              C_NULL, idx, opt
        )
        err = ccall(set_opt, Cint, (Cstring, Ptr{Cvoid}, Ptr{Cvoid}),
              preset, idx, opt
        )

        # manual overrides
        opt[].min_dp_max = minblock-10
        opt[].flag |= (NO_DIAG|CIGAR)

        return opt
    end
end

"""
	struct PanContigs
		name     :: T
		sequence :: T
	end

A synonym for a consensus sequence of `Block`.
"""
struct PanContigs{T <: AbstractArray{S} where S <: AbstractString}
    name     :: T
    sequence :: T
end

# ------------------------------------------------------------------------
# wrapper

# reference index
"""
	makeindex(w, k, names, sequence; bucketbits::Int=14)

Given a window size `w` and kmer length `k`, return a handle to a minimizer index
for sequences `sequence`.
"""
function makeindex(w, k, names, sequence; bucketbits::Int=14)
    @assert length(names) == length(sequence)

    return Index(
        ccall(idx_str, Ptr{Cvoid}, (Cint, Cint, Cint, Cint, Cint, Ptr{Cstring}, Ptr{Cstring}),
          w, k, 0, bucketbits, length(names), sequence, names
        )
    )
end

function show(ref::Index)
    ccall(idx_stat, Cvoid, (Ptr{Cvoid},), ref.handle)
end

"""
	freeindex()

Free memory associated to an opaque handle to a sequence index.
"""
function freeindex(ref::Index)
    ccall(idx_free, Cvoid, (Ptr{Cvoid},), ref.handle)
end

# thread buffer
"""
	makebuffer()

Return an opaque handle thread buffer
"""
function makebuffer()
    return Buffer(ccall(tbuf_init, Ptr{Cvoid}, ()))
end

"""
	freebuffer()

Free memory associated to an opaque handle thread buffer
"""
function freebuffer(buf::Buffer)
    return ccall(tbuf_free, Cvoid, (Ptr{Cvoid},), buf.handle)
end

# cigar string
const cigartypes = ['M','I','D','N','S','H','P','E','X']
unpackcigar(base::Ptr{UInt32}, n) = [let
    cg = unsafe_load(base,i)
    (Int(cg>>4), cigartypes[cg&0xf+1])
end for i in 1:n]

# alignment
"""
	align(ref::PanContigs, qry::PanContigs, minblock::Int, preset::String)

Call into minimap to align the set of blocks `qry` to blocks `ref`.
Preset should be a string ∈ ["asm5","asm10","asm20"]. See minimap2 manual for details.
This is probably the function you want.
If you call into the function specifically, all memory management is taken care of for you.
"""
function align(ref::PanContigs, qry::PanContigs, minblock::Int, preset::String)
    iopt = Ref{IndexOptions}(IndexOptions())
    mopt = MapOptions(iopt, minblock, preset)
    len  = Ref{Cint}(0)

    index  = makeindex(iopt[].w, iopt[].k, ref.name, ref.sequence)

    buffer = makebuffer()

    ccall(update_opt, Cvoid, (Ptr{MapOptions}, Ptr{Cvoid}), mopt, index.handle)

    homologs = Alignment[]
    for (i,seq) in enumerate(qry.sequence)
        hits = ccall(minimap, Ptr{Record},
                (Ptr{Cvoid}, Cint, Ptr{Cchar}, Ptr{Cint}, Ptr{Cvoid}, Ptr{Cvoid}, Cstring),
                index.handle, length(seq), seq, len, buffer.handle, mopt, qry.name[i]
        )

        for n in 1:len[]
            hit = unsafe_load(hits, n)
            aln = unsafe_load(hit.p, 1)

            j = hit.rid+1 # minimap2 0 indexes
            push!(homologs, Alignment(
                Hit(qry.name[i], length(seq), hit.qs+1, hit.qe, nothing),
                Hit(ref.name[j], length(ref.sequence[j]), hit.rs+1, hit.re, nothing),
                Int(hit.mlen),
                Int(hit.blen),
                Int(hit.pack & 0xff),
                ((hit.pack >> 10) & 1) == 0,
                unpackcigar(Ptr{UInt32}(hit.p+Core.sizeof(Extra)), aln.n_cigar),
                1 - ccall(divergence, Cdouble, (Ptr{Record},), Ref(hit)),
                nothing
            ))

            Libc.free(hit.p)
        end
        Libc.free(hits)
    end

    freebuffer(buffer)
    freeindex(index)

    return homologs
end

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
@testset "Elements" begin

  element_result = Element("El 1")
  @test element_result.Name == "El 1"
  @test name(element_result) == "El 1"

  el_2 = JSON.parse("""
    {
      "Name": "El 2",
      "UniqueName": "[El 2]",
      "Attributes": {
        "Caption": "El 2"
      }
    }
""")

  @test name(Element(el_2)) == "El 2"
  @test Element(el_2) != element_result

end


text:
module Inference

using ..Core, ..Core.RandomVariables, ..Utilities
using ..Core.RandomVariables: Metadata, _tail
using Distributions, Libtask, Bijectors
using ProgressMeter, LinearAlgebra
using ..Turing: PROGRESS, CACHERESET, AbstractSampler
using ..Turing: Model, runmodel!, get_pvars, get_dvars,
    Sampler, SampleFromPrior, SampleFromUniform,
    Selector
using ..Turing: in_pvars, in_dvars, Turing
using StatsFuns: logsumexp
using Random: GLOBAL_RNG, AbstractRNG
import AdvancedHMC; const AHMC = AdvancedHMC

import ..Turing: getspace
import Distributions: sample
import ..Core: getchunksize, getADtype
import ..Utilities: Sample, save, resume

export  InferenceAlgorithm,
        Hamiltonian,
        AbstractGibbs,
        GibbsComponent,
        StaticHamiltonian,
        AdaptiveHamiltonian,
        SampleFromUniform,
        SampleFromPrior,
        MH,
        Gibbs,      # classic sampling
        HMC,
        SGLD,
        SGHMC,
        HMCDA,
        NUTS,       # Hamiltonian-like sampling
        DynamicNUTS,
        IS,
        SMC,
        CSMC,
        PG,
        PIMH,
        PMMH,
        IPMCMC,  # particle-based sampling
        assume,
        observe,
        step,
        WelfordVar,
        WelfordCovar,
        NaiveCovar,
        get_var,
        get_covar,
        add_sample!,
        reset!

#######################
# Sampler abstraction #
#######################
abstract type AbstractAdapter end
abstract type InferenceAlgorithm end
abstract type Hamiltonian{AD} <: InferenceAlgorithm end
abstract type StaticHamiltonian{AD} <: Hamiltonian{AD} end
abstract type AdaptiveHamiltonian{AD} <: Hamiltonian{AD} end

getchunksize(::T) where {T <: Hamiltonian} = getchunksize(T)
getchunksize(::Type{<:Hamiltonian{AD}}) where AD = getchunksize(AD)
getADtype(alg::Hamiltonian) = getADtype(typeof(alg))
getADtype(::Type{<:Hamiltonian{AD}}) where {AD} = AD

"""
    mh_accept(H::T, H_new::T, log_proposal_ratio::T) where {T<:Real}

Peform MH accept criteria with log acceptance ratio. Returns a `Bool` for acceptance.

Note: This function is only used in PMMH.
"""
function mh_accept(H::T, H_new::T, log_proposal_ratio::T) where {T<:Real}
    return log(rand()) + H_new < H + log_proposal_ratio, min(0, -(H_new - H))
end

# Concrete algorithm implementations.
include("sghmc.jl")
include("hmc.jl")
include("mh.jl")
include("is.jl")
include("AdvancedSMC.jl")
include("gibbs.jl")

getspace(::IS) = ()
getspace(::Type{<:IS}) = ()
for alg in (:SMC, :PG, :PMMH, :IPMCMC, :MH)
    @eval getspace(::$alg{space}) where {space} = space
    @eval getspace(::Type{<:$alg{space}}) where {space} = space
    @eval getspace(::Type{<:$alg{space}}) where {space} = space
end
for alg in (:HMC, :HMCDA, :NUTS, :SGLD, :SGHMC)
    @eval getspace(::$alg{<:Any, space}) where {space} = space
    @eval getspace(::Type{<:$alg{<:Any, space}}) where {space} = space
end

@inline floatof(::Type{T}) where {T <: Real} = typeof(one(T)/one(T))
@inline floatof(::Type) = Real

@inline Turing.Core.get_matching_type(spl::Turing.Sampler, vi::Turing.RandomVariables.VarInfo, ::Type{T}) where {T <: AbstractFloat} = floatof(eltype(vi, spl))
@inline Turing.Core.get_matching_type(spl::Turing.Sampler{<:Hamiltonian}, vi::Turing.RandomVariables.VarInfo, ::Type{TV}) where {T, N, TV <: Array{T, N}} = Array{Turing.Core.get_matching_type(spl, vi, T), N}
@inline Turing.Core.get_matching_type(spl::Turing.Sampler{<:Union{PG, SMC}}, vi::Turing.RandomVariables.VarInfo, ::Type{TV}) where {T, N, TV <: Array{T, N}} = TArray{T, N}

## Fallback functions

# utility funcs for querying sampler information
require_gradient(spl::Sampler) = false
require_particles(spl::Sampler) = false

assume(spl::Sampler, dist::Distribution) =
error("Turing.assume: unmanaged inference algorithm: $(typeof(spl))")

observe(spl::Sampler, weight::Float64) =
error("Turing.observe: unmanaged inference algorithm: $(typeof(spl))")

## Default definitions for assume, observe, when sampler = nothing.
function assume(spl::A,
    dist::Distribution,
    vn::VarName,
    vi::VarInfo) where {A<:Union{SampleFromPrior, SampleFromUniform}}

    if haskey(vi, vn)
        r = vi[vn]
    else
        r = isa(spl, SampleFromUniform) ? init(dist) : rand(dist)
        push!(vi, vn, r, dist, spl)
    end
    # NOTE: The importance weight is not correctly computed here because
    #       r is genereated from some uniform distribution which is different from the prior
    # acclogp!(vi, logpdf_with_trans(dist, r, istrans(vi, vn)))

    r, logpdf_with_trans(dist, r, istrans(vi, vn))

end

function assume(spl::A,
    dists::Vector{T},
    vn::VarName,
    var::Any,
    vi::VarInfo) where {T<:Distribution, A<:Union{SampleFromPrior, SampleFromUniform}}

    @assert length(dists) == 1 "Turing.assume only support vectorizing i.i.d distribution"
    dist = dists[1]
    n = size(var)[end]

    vns = map(i -> VarName(vn, "[$i]"), 1:n)

    if haskey(vi, vns[1])
        rs = vi[vns]
    else
        rs = isa(spl, SampleFromUniform) ? init(dist, n) : rand(dist, n)

        if isa(dist, UnivariateDistribution) || isa(dist, MatrixDistribution)
            for i = 1:n
                push!(vi, vns[i], rs[i], dist, spl)
            end
            @assert size(var) == size(rs) "Turing.assume: variable and random number dimension unmatched"
            var = rs
        elseif isa(dist, MultivariateDistribution)
            for i = 1:n
                push!(vi, vns[i], rs[:,i], dist, spl)
            end
            if isa(var, Vector)
                @assert length(var) == size(rs)[2] "Turing.assume: variable and random number dimension unmatched"
                for i = 1:n
                    var[i] = rs[:,i]
                end
            elseif isa(var, Matrix)
                @assert size(var) == size(rs) "Turing.assume: variable and random number dimension unmatched"
                var = rs
            else
                @error("Turing.assume: unsupported variable container"); error()
            end
        end
    end

    # acclogp!(vi, sum(logpdf_with_trans(dist, rs, istrans(vi, vns[1]))))

    var, sum(logpdf_with_trans(dist, rs, istrans(vi, vns[1])))

end

observe(::Nothing,
        dist::T,
        value::Any,
        vi::VarInfo) where T = observe(SampleFromPrior(), dist, value, vi)

function observe(spl::A,
    dist::Distribution,
    value::Any,
    vi::VarInfo) where {A<:Union{SampleFromPrior, SampleFromUniform}}

    vi.num_produce += one(vi.num_produce)
    Turing.DEBUG && @debug "dist = $dist"
    Turing.DEBUG && @debug "value = $value"

    # acclogp!(vi, logpdf(dist, value))
    logpdf(dist, value)

end

function observe(spl::A,
    dists::Vector{T},
    value::Any,
    vi::VarInfo) where {T<:Distribution, A<:Union{SampleFromPrior, SampleFromUniform}}

    @assert length(dists) == 1 "Turing.observe only support vectorizing i.i.d distribution"
    dist = dists[1]
    @assert isa(dist, UnivariateDistribution) || isa(dist, MultivariateDistribution) "Turing.observe: vectorizing matrix distribution is not supported"
    if isa(dist, UnivariateDistribution)  # only univariate distributions support broadcast operation (logpdf.) by Distributions.jl
        # acclogp!(vi, sum(logpdf.(Ref(dist), value)))
        sum(logpdf.(Ref(dist), value))
    else
        # acclogp!(vi, sum(logpdf(dist, value)))
        sum(logpdf(dist, value))
    end

end

##############
# Utilities  #
##############

# VarInfo to Sample
Sample(vi::VarInfo) = Sample(0.0, todict(vi))
function todict(vi::VarInfo)
    value = todict(vi.metadata, vi)
    value[:lp] = getlogp(vi)
    return value
end
function todict(md::Metadata, vi::VarInfo)
    value = Dict{Symbol, Any}() # value is named here because of Sample has a field called value
    for vn in keys(md.idcs)
        value[Symbol(vn)] = vi[vn]
    end
    return value
end
function todict(metadata::NamedTuple{names}, vi::VarInfo) where {names}
    length(names) === 0 && return Dict{Symbol, Any}()
    f = names[1]
    mdf = getfield(metadata, f)
    return merge(todict(mdf, vi), todict(_tail(metadata), vi))
end

# VarInfo, combined with spl.info, to Sample
function Sample(vi::AbstractVarInfo, spl::Sampler)
    s = Sample(vi)
    if haskey(spl.info, :adaptor)
        s.value[:lf_eps] = AHMC.getϵ(spl.info[:adaptor])
    end
    if haskey(spl.info, :eval_num)
        s.value[:eval_num] = spl.info[:eval_num]
    end
    return s
end

getspace(spl::Sampler) = getspace(typeof(spl))
getspace(::Type{<:Sampler{Talg}}) where {Talg} = getspace(Talg)

end # module


text:
"""
CSCI 480 Final Project
Authors: Josh Myers-Dean, Mina Shin, Ryan Wells
Winter 2020
"""

module WWURay

export main

using FileIO
using Images
using StaticArrays
using LinearAlgebra

push!(LOAD_PATH, pwd())
include("GfxBase.jl")
include("Lights.jl")
include("Materials.jl")
include("OBJMeshes.jl")
include("Scenes.jl")
include("Cameras.jl")
include("TestScenes.jl")

using .GfxBase
using .Lights
using .Materials

import .Scenes
import .Scenes.Scene
import .Scenes.HitRecord
import .Cameras
import .TestScenes

# Ray-Scene intersection:
""" Find the closest intersection point among all objects in the scene
along a ray, constraining the search to values of t between tmin and tmax. """
function closest_intersect(objects::Array{Any, 1}, ray::Ray, tmin, tmax)
    closestHit = nothing
    # various cases for error handling
    for i ∈ objects
        hit = Scenes.ray_intersect(ray, i)
        if isnothing(hit) && isnothing(closestHit)
            continue
        elseif !isnothing(hit) && isnothing(closestHit)
            if in_bounds(hit.t, tmin, tmax)
                closestHit = hit
            else
                continue
            end
        elseif isnothing(hit) && !isnothing(closestHit)
            if in_bounds(closestHit.t, tmin, tmax)
                closestHit = closestHit
            else
                continue
            end
        elseif hit.t < closestHit.t
            if in_bounds(hit.t, tmin, tmax)
                closestHit = hit
            else
                continue
            end
        end
    end
    return closestHit
end

# helper function to check if a t value is valid
function in_bounds(t, tmin, tmax)
    return (t >= tmin) && (t <= tmax)
end

""" Trace a ray from orig along ray through scene, using Whitted recursive raytracing
limited to rec_depth recursive calls. """
function traceray(scene::Scene, ray::Ray, tmin, tmax, rec_depth=1)
    # find the closest intersection
    closest_hitrec = closest_intersect(scene.objects, ray, tmin, tmax)
    if isnothing(closest_hitrec)
        return scene.background
    end
    # hitrec vlaues
    object = closest_hitrec.object
    point = closest_hitrec.intersection
    normal = closest_hitrec.normal
    material = object.material
    shader = material.shading_model

    # determine the color of the pixel
    local_color = determine_color(shader, object.material, ray, closest_hitrec, scene)
    mirror_co = object.material.mirror_coeff
    # ray direction as per the lecture
    ray_dir = 2* normal*dot(normal, -ray.direction)+ray.direction
    mirror_ray = Ray(point, normalize(ray_dir))

    # recurse until rec depth is hit or if there is no mirror coef
    if (mirror_co == 0) || (rec_depth <= 0)
        return local_color
    end
    if scene.id == 1
        tpRay = Vec3(0.545, 0.0, -0.545)
        #rotate mirror_ray's direction (90 degrees clockwise) without changing it's origin
        new_dir = Vec3(-ray_dir[3], -ray_dir[2], ray_dir[1])

                  #new origin      new direction
        mirror_ray = Ray(point - tpRay, -normalize(new_dir))
    elseif scene.id == 4
            tpRay = Vec3(-2.5,0,0)
            new_dir = Vec3(ray_dir[1], -ray_dir[2], -ray_dir[3])
            mirror_ray = Ray(point - tpRay, -normalize(new_dir))
    end
    reflected_color = traceray(scene, mirror_ray, 1e-8, Inf, rec_depth-1)
    return (local_color*(1-mirror_co))+(reflected_color*mirror_co)

end

""" Determine the color of interesction point described by hitrec
Flat shading - just color the pixel the material's diffuse color """
function determine_color(shader::Flat, material::Material, ray::Ray, hitrec::HitRecord, scene::Scene)
    get_diffuse(material, hitrec.uv)
end
""" Normal shading - color-code pixels according to their normals """
function determine_color(shader::Normal, material::Material, ray::Ray, hitrec::HitRecord, scene::Scene)
    normal_color = normalize(hitrec.normal) / 2 .+ 0.5
    RGB{Float32}(normal_color...)
end

""" Determine the color of a physical (Lambertian, BlinnPhong, etc.) surface """
function determine_color(shader::PhysicalShadingModel, material::Material, ray::Ray, hitrec::HitRecord, scene::Scene)
    # start with black color, increment if point is not shadowed
    lightAmount = RGB(0.0,0.0,0.0)
    for i ∈  scene.lights
        # call light source depending on type
        if !is_shadowed(scene, i, hitrec.intersection)
            lightAmount += shade_light(shader, material, ray, hitrec, i, scene)
        end
    end
    return lightAmount

end

""" shade_light(shader, material, ray, hitrec, light, scene)
Determine the color contribution of the given light along the given ray.
Color depends on the material, the shading model (shader), properties of the intersection
given in hitrec, """
# just diffuse shading
function shade_light(shader::Lambertian, material::Material, ray::Ray, hitrec, light, scene)
    l = get_diffuse(material, hitrec.uv) * light.intensity*max(0, dot(hitrec.normal, normalize(Lights.light_direction(light, hitrec.intersection))))
end

""" Blinn-Phong surface shading """
function shade_light(shader::BlinnPhong, material::Material, ray::Ray, hitrec, light, scene)
    # diffuse coefficent
    diffuse = get_diffuse(material, hitrec.uv) * light.intensity*max(0, dot(hitrec.normal, normalize(Lights.light_direction(light, hitrec.intersection))))
    # half vector
    h = normalize((-ray.direction)+Lights.light_direction(light, hitrec.intersection))
    # specular reflection
    spec = shader.specular_color*light.intensity*(max(0, dot(hitrec.normal, h))^(shader.specular_exp))
    return diffuse+spec
end

""" Determine whether point is in shadow wrt light """
function is_shadowed(scene, light::DirectionalLight, point::Vec3)
    if closest_intersect(scene.objects, Ray(point, Lights.light_direction(light,point)), 1e-8, Inf) != nothing
        return true
    else
        return false
    end
end

function is_shadowed(scene, light::PointLight, point::Vec3)
    if closest_intersect(scene.objects, Ray(point, Lights.light_direction(light, point)), 1e-8, 1) != nothing
        return true
    else
        return false
    end
end

# Main loop:
function main(scene, height, width, out)
    if scene == 1
	    scene = TestScenes.portalScene1(height, width)[1]
        camera = TestScenes.portalScene1(height, width)[2]
        getImages(scene, camera, height, width, out, 1, 0)
    elseif scene == 2
	    scene = TestScenes.portalScene2(height, width)[1]
        camera = TestScenes.portalScene2(height, width)[2]
        getImages(scene, camera, height, width, out, 1, -1)
    elseif scene == 3
	    scene = TestScenes.portalScene3(height, width)[1]
        camera = TestScenes.portalScene3(height, width)[2]
        getImages(scene, camera, height, width, out, 1, 0)
    elseif scene == 4
        scene = TestScenes.portalScene4(height, width)[1]
        camera = TestScenes.portalScene4(height, width)[2]
        getImages(scene, camera, height, width, out, 1, 0)
   end
end

# gets 60 images with varying camera position
function getImages(scene, camera, height, width, out, amountx, amountz)
    RANGLE = 45         # Rotate angle size
    NUMS = 60           # Number of images
    for idx = 1:NUMS
        # change camera
        T = [
            1.0 0.0 0.0 inv(tan(RANGLE/NUMS))/30 * amountx;
            0.0 1.0 0.0 0.0;
            0.0 0.0 1.0 inv(tan(RANGLE/NUMS))/30 * amountz;
            0.0 0.0 0.0 1.0
        ]

        R = TestScenes.rotate(RANGLE/NUMS, 0, 1, 0);
        newcam = T * Vec4(camera.eye[1], camera.eye[2], camera.eye[3], 1)
        camera.eye = Vec3(newcam[1], newcam[2], newcam[3])

        #get objects
        objs = scene.objects
        # Create a blank canvas to store the image:
        canvas = zeros(RGB{Float32}, height, width)

        # for each pixel iterate over each object
        # color the canavs accordingly
        for i = 1:height
            for j = 1:width
                ray:: Ray = Cameras.pixel_to_ray(camera, i, j)
                color = traceray(scene, ray, 1, Inf, 8)
                canvas[i,j] = color
            end
        end

        # clamp canvas to valid range:
        clamp01!(canvas)

        outfolder = string("video/", out)
        outfname = string(idx, ".png")
        save(File(format"PNG", string(outfolder, outfname)), colorview(RGB, canvas))

    end # end of idx loop
end

end # module WWURay


text:
using Colors

@testset "din99" begin
    # Test data from the DIN 6176 specification
    testdata = [
        ((50,  10,  10), (61.43,   9.70,   3.76)),
        ((50,  50,  50), (61.43,  28.64,  11.11)),
        ((50, -10,  10), (61.43,  -5.57,   7.03)),
        ((50, -50,  50), (61.43, -17.22,  21.75)),
        ((50, -10, -10), (61.43,  -9.70,  -3.76)),
        ((50, -50, -50), (61.43, -28.64, -11.11)),
        ((50,  10, -10), (61.43,   5.57,  -7.03)),
        ((50,  50, -50), (61.43,  17.22, -21.75)),
        (( 0,   0,   0), ( 0,      0,      0)),
        ((100,  0,   0), (100,     0,      0))]

    # A high error threshold has been chosen because converting from DIN99
    # to CIELAB with only two decimal places of accuracy yields fairly inaccurate
    # results due to the exponentiation.
    conveps = 0.05
    diffeps = 0.01
    let a, b, metric
        for (i, (a, b)) in enumerate(testdata)
            converted = convert(DIN99, Lab(a...))
            test = DIN99(b...)

            @test (abs(converted.l - test.l) < conveps)
            @test (abs(converted.a - test.a) < conveps)
            @test (abs(converted.b - test.b) < conveps)

            converted = convert(Lab, DIN99(b...))
            test = Lab(a...)

            @test (abs(converted.l - test.l) < conveps)
            @test (abs(converted.a - test.a) < conveps)
            @test (abs(converted.b - test.b) < conveps)

            # This is not a real test of the color difference metric, but at least
            # makes sure it isn't doing anything really crazy.
            metric = DE_DIN99()
            @test (abs(colordiff(convert(DIN99, Lab(a...)), DIN99(b...); metric=metric)) < diffeps)
        end
    end

    # From #256
    for c in (RGB(1, 0.5, 0),)
        @test colordiff(convert(RGB,convert(DIN99d, c)), c) < 1e-3
    end
    for c in (RGBA(1,0.5,0,0.4),)
        cc = convert(RGBA,convert(DIN99dA, c))
        @test colordiff(color(cc), color(c)) < 1e-3
        @test cc.alpha ≈ c.alpha atol=1e-3
    end
end  # @testset


text:
module TestExamples

import ShowGraphviz
import TaskDAGAnalyzers
import TaskDAGRecorders
using TaskDAGAnalyzers: Examples
using TaskDAGAnalyzers.Internal: SpawnNode, SequentialNode, DurationNS, foreachnode
using Test

function check_smoke(f)
    ctx = TaskDAGRecorders.sample(f)
    dag = TaskDAGAnalyzers.dag(ctx)
    dot = ShowGraphviz.DOT(dag)
    src = sprint(show, "text/plain", dot)
    @test occursin("digraph", src)
    txt = sprint(show, "text/plain", TaskDAGAnalyzers.summary(ctx))
    @test occursin("work", txt)
    @test occursin("span", txt)
    @test occursin("parallelism", txt)
end

function test_smoke()
    @testset for f in [Examples.dac, Examples.dac4]
        check_smoke(f)
    end
end

function set_uniform_unit_work!(dag, unit)
    foreachnode(dag) do node
        if node isa SequentialNode
            node.data = unit
        elseif node isa SpawnNode
            node.data = DurationNS(typemax(Int))  # this should be ignored
        end
    end
end

# TODO: don't depend on noisefloor
function check_unit_span(f)
    ctx = TaskDAGRecorders.sample(f, 4)[end]
    dag = TaskDAGAnalyzers.dag(ctx)
    unit = DurationNS(123_456_789)
    set_uniform_unit_work!(dag, unit)
    @test TaskDAGAnalyzers.span(dag) == unit
end

function test_unit_span()
    @testset for f in [Examples.dac, Examples.dac4]
        check_unit_span(f)
    end
end

function test_dac_work()
    ctx = TaskDAGRecorders.sample(Examples.dac, 4)[end]
    dag = TaskDAGAnalyzers.dag(ctx)
    unit = DurationNS(123_456_789)
    set_uniform_unit_work!(dag, unit)
    expected_work = DurationNS(unit.value * 4)
    @test TaskDAGAnalyzers.work(dag) == expected_work
end

function test_dac4_work()
    ctx = TaskDAGRecorders.sample(Examples.dac4, 4)[end]
    dag = TaskDAGAnalyzers.dag(ctx)
    unit = DurationNS(123_456_789)
    set_uniform_unit_work!(dag, unit)
    expected_work = DurationNS(unit.value * 16)
    @test TaskDAGAnalyzers.work(dag) == expected_work
end

end  # module


text:
using Documenter

makedocs(
    sitename = "llir/llvm",
    format = Documenter.HTML(),
    pages = [
        "index.md",
        "User Guide" => [
            "user-guide/basic.md",
            "user-guide/control.md",
            "More Function" => [
                "user-guide/func/linkage.md",
                "user-guide/func/vaarg.md",
                "user-guide/func/overloading.md",
                "user-guide/func/closure.md",
                "user-guide/func/exception.md",
            ],
            "user-guide/types.md",
            "user-guide/support.md"
        ]
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#

max_stars_repo_name: dannypsnl/document

max_issues_repo_name: dannypsnl/document

max_forks_repo_name: dannypsnl/document


text:
@testset "Parsing" begin
    @testset "tryparse" begin
        @test tryparse(c, Int, 123) == (123, nothing)
        @test tryparse(c, Vector{UInt16}, Int[1, 2, 3]) == (UInt16[1, 2, 3], nothing)
        val, e = @test_logs (:error, r"Parsing failed") tryparse(c, Int, Dict())
        @test val === nothing
        @test e isa MethodError
    end

    @testset "Snowflake" begin
        # Snowflakes usually come in as strings.
        s = snowflake(string(typemax(UInt64)))
        @test isa(s, Snowflake)
        # Snowflake === UInt64.
        @test s == typemax(UInt64)

        # https://discordapp.com/developers/docs/reference#snowflakes-snowflake-id-format-structure-left-to-right
        s = Snowflake(0x06ecefa78e42000c)
        @test snowflake2datetime(s) == DateTime(2018, 10, 9, 1, 55, 31, 1)
        @test worker_id(s) == 0x01
        @test process_id(s) == 0x00
        @test increment(s) == 0x0c
    end

    @testset "DateTime" begin
        # Discord sends dates in some weird, inconsistent ways.
        d = datetime("2018-10-08T05:20:22.782643+00:00")
        @test d == DateTime(2018, 10, 8, 5, 20, 22, 782)
        d = datetime("2018-10-08T05:20:22.782Z")
        @test d == DateTime(2018, 10, 8, 5, 20, 22, 782)

        # But sometimes they also send nice millisecond Unix timestamps.
        d = datetime(1541288588543)
        @test d == DateTime(2018, 11, 3, 23, 43, 08, 543)
    end
end


text:
"""
    glaser_liu_rokhlin_gauss_nodes(π::P, x0=pqr_start(P,degree(π));
                                         symmetry=pqr_symmetry(P),
                                         m=5)

Find Gauss nodes and weights of the  basis  polynomial  `π = basis(P,n)`. The nodes  are the  roots of the  polynomial and the weights are computed  from known formulas.

The method from  [Glaser,  Liu, and Rokhlin](DOI: 10.1137/06067016X) is  used. This is  an  O(n)  method (whereas the method  basedon the Jacobi  matrix is O(n^2)). This method fits easily  into the framework  provided  through  the `AbstractCCOP` types.  The [FastGaussQuadrature](https://github.com/JuliaApproximation/FastGaussQuadrature.jl) package provides even more efficient  algorithms and  pays attention to numerical  issues, examples of  which  can  be found in  [Hale and Townsend](https://core.ac.uk/download/pdf/9403469.pdf). The `FastGuassQuadrature` package is used when  available.

An example:

```
function  gauss_nodes_weights(P::Type{<:Jacobi{α,β}}, n)  where {α,β}
    # we don't  have a  good starting point  unless α=β
    if α == β
        xs, ws = glaser_liu_rokhlin_gauss_nodes(basis(MonicJacobi{α,β},n))
        λ = kn(P,n)^2
        xs, ws/λ
    else
        J = jacobi_matrix(P, n)
        eig = eigen(J, extrema(P)...)
        wts =  π̃βn(P,0) * (eig.vectors[1,:]).^2
        eig.values,  wts
    end

end

pqr_start(P::Type{MonicJacobi{α,α}}, n) where {α} =  0.0 # zero(eltype(P))
pqr_symmetry(P::Type{<:MonicJacobi{α,α}}) where {α} = true
function pqr_weight(P::Type{<:MonicJacobi{α,α}}, n, x, dπx) where {α}
    # https://www.chebfun.org/publications/HaleTownsend2013a.pdf
    # Eqn  (1.4)
    # Cnαβ should be  computed  using asymptotic formula  for larger n (§3.2.3)
    # XXX TThere i ss ome constant  that  makes this not work....
    β = α
    Cnαβ = 2^(α+β+1)
    Cnαβ *= gamma(n + α +  1) * gamma(n + β + 1) / gamma(n + α + β + 1)
    Cnαβ /= gamma(1 + n)
    val = Cnαβ / (1-x^2) / dπx^2
    val
end
```

"""
glaser_liu_rokhlin_gauss_nodes()

"""
    pqr

Compute  `p,q,r,p',q',r'` where `p(x)y'' + q(x)y' + r(x)  = 0`. Uses
the `a,b,c,d,e` values characterizing the  system.

"""
function pqr(P::Type{<:AbstractCCOP}, n, x)
    a,b,c,d,e = abcde(P)

    p =  a*x^2 + b*x  +  c
    dp = 2a*x + b
    q = d*x +  e
    dq = d
    r  = -(a*n*(n-1)  + d*n)
    dr = 0
    (p,q,r,dp,dq,dr)
end

"""
    pqr_start(p::P, n)

A starting  value for finding the gauss nodes
"""
pqr_start(::Type{P}, n) where {P <: AbstractCCOP} = 0

"""
    pqr_symmetry(p::P)

Boolean to specify if symmetry should apply to output
"""
pqr_symmetry(::Type{P}) where {P <: AbstractCCOP} = false

"""
    pqr_weight(p::P, x, dx)

Compute weight from x, dπx, x a node and dπx the derivative's value at the node
"""
pqr_weight(p::P, n, x, dx) where {P} = throw(MethodError())

# We just use clenshaw to  evaluate  `p(x), dp(x)` for Newton's method below,
# but implementing the following can be more efficient
# ## Evaluate basis(P,n)  and its  derivative using the three-point recursion representation
# function orthogonal_polyval_derivative(p::P, x::S) where {P <: AbstractOrthogonalPolynomial, S}

#     T = eltype(p)
#     oS = one(x)
#     R = eltype(one(T) * oS)

#     length(p) == 0 && return (zero(R), zero(R))
#     d = length(p) - 1

#     # P0,  P_{-1} case
#     pn,  pn_1 = P0(P,x), zero(x)
#     dpn, dpn_1  = dP0(P,x), zero(x)

#     ptot = p[0]*P0(P,x)
#     dptot  = p[0]*dP0(P,x)

#     ## we compute  forward here, finding pi, dp_i,  p_{i+1}, dp_{i+1}, ...
#     for i in 0:d-1
#         an, bn, cn=  An(p,i), Bn(p,i), Cn(p,i)
#         dpn_1, dpn = dpn,  an*pn + (an*x + bn)*dpn + cn*dpn_1
#         pn_1,  pn =  pn, (an*x + bn)*pn + cn*pn_1
#         ptot += p[i+1]*pn
#         dptot += p[i+1]*dpn
#     end
#     return (ptot, dptot)
# end

# run Newton's method to find zero of p
function newton(p::P, x0::S) where  {P <: AbstractCCOP, S}
    maxsteps = 25
    dp  = derivative(p)
    while maxsteps > 0
        px, dpx = p(x0),  dp(x0) #orthogonal_polyval_derivative(p, x0)
        Δ = px/dpx
        isnan(Δ) && return (x0, dpx)
        if isinf(Δ)
            x0 += sqrt(eps())  # give a nudge from minimum or maximum
            continue
        end
        x0 -= Δ
        min(abs(Δ),abs(px)) <= sqrt(eps(S))/100 && return  x0, dp(x0)# orthogonal_polyval_derivative(p, x0)[2]
        maxsteps -= 1
    end
    @warn "Newton's method did not converge from $x0"
    NaN
end

# Second-order Runge-Kutta with global truncation error of h^2
function RK(t0,  x0, F, h, n)
    k0 = h*F(t0,x0)
    local x1
    for i in 1:n
        t1 = t0 + h
        k1 = h * F(t0+h,x0+k0)
        x1 = x0 + 1/2 * (k0+k1)
        t0, x0, k0 = t1, x1, k1
    end
    x1
end

function prufer(Π::Type{P},n) where {P <: AbstractCCOP}
    (θ, x) -> begin
        dom = domain(Π)
        a, b = first(dom)+eps(), last(dom)-eps()
        x = clamp(x, a, b)
        p,q,r,dp,dq,dr = pqr(Π, n, x)
        -inv(sqrt(r/p) + (dr*p - dp*r + 2r*q)/(2r*p) * sin(2θ)/2)
    end
end

## Compute gauss node using  algorithm from
## A Fast Algorithm for the Calculation of the Roots of Special Functions
## by Glaser, Liu, Rokhlin
## DOI: 10.1137/06067016X
## specialized to the orthogonal polynomial case
function glaser_liu_rokhlin_gauss_nodes(π::P, x0=pqr_start(P,degree(π));
                                        symmetry=pqr_symmetry(P),
                                        m=5) where {P <: AbstractCCOP}

    n = degree(π)
    F = prufer(P, n)
    x = float(x0)

    # step 1 to get initial might be newton or might be RK
    if symmetry && iseven(n)
        # a max at π(x0)
        x = RK(0, x0, F, -(pi/2)/m, m)
    end
    x1, dπ1 = newton(π, x)
    rts, dπrts = [x1], [dπ1]
    N = symmetry ? div(n+1,2) : n
    for i in 2:N
        xx = RK(pi/2, rts[end], F, -pi/m, m)
        x, dπx = newton(π, xx)
        push!(rts, x)
        push!(dπrts, dπx)
    end
    if symmetry
        if iseven(n)
            rts = vcat(-reverse(rts), rts)
            dπrts = vcat(-reverse(dπrts), dπrts)
        else
            rts = vcat(-reverse(rts[2:end]), rts)
            dπrts = vcat(-reverse(dπrts[2:end]), dπrts)
        end
    end
    # get weights
    weights = pqr_weight.(P, n, rts, dπrts)
    rts,  weights
end