text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TestSymmetric

using Test, LinearAlgebra, SparseArrays, Random

Random.seed!(101)

@testset "Pauli σ-matrices: $σ" for σ in map(Hermitian,
        Any[ [1 0; 0 1], [0 1; 1 0], [0 -im; im 0], [1 0; 0 -1] ])
    @test ishermitian(σ)
end

@testset "Hermitian matrix exponential/log" begin
    A1 = randn(4,4) + im*randn(4,4)
    A2 = A1 + A1'
    @test exp(A2) ≈ exp(Hermitian(A2))
    @test log(A2) ≈ log(Hermitian(A2))
    A3 = A1 * A1' # posdef
    @test exp(A3) ≈ exp(Hermitian(A3))
    @test log(A3) ≈ log(Hermitian(A3))

    A1 = randn(4,4)
    A3 = A1 * A1'
    A4 = A1 + transpose(A1)
    @test exp(A4) ≈ exp(Symmetric(A4))
    @test log(A3) ≈ log(Symmetric(A3))
    @test log(A3) ≈ log(Hermitian(A3))
end

@testset "Core functionality" begin
    n = 10
    areal = randn(n,n)/2
    aimg  = randn(n,n)/2
    @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int)
        a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
        asym = transpose(a) + a                 # symmetric indefinite
        aherm = a' + a                 # Hermitian indefinite
        apos  = a' * a                 # Hermitian positive definite
        aposs = apos + transpose(apos)        # Symmetric positive definite
        ε = εa = eps(abs(float(one(eltya))))

        x = randn(n)
        y = randn(n)
        b = randn(n,n)/2
        x = eltya == Int ? rand(1:7, n) : convert(Vector{eltya}, eltya <: Complex ? complex.(x, zeros(n)) : x)
        y = eltya == Int ? rand(1:7, n) : convert(Vector{eltya}, eltya <: Complex ? complex.(y, zeros(n)) : y)
        b = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(b, zeros(n,n)) : b)
        @testset "basic ops" begin
            @testset "constructor" begin
                @test Symmetric(Symmetric(asym, :U))     === Symmetric(asym, :U)
                @test Hermitian(Hermitian(aherm, :U))    === Hermitian(aherm, :U)
                @test Symmetric(Symmetric(asym, :U), :U) === Symmetric(asym, :U)
                @test Hermitian(Hermitian(aherm, :U), :U) === Hermitian(aherm, :U)
                @test_throws ArgumentError Symmetric(Symmetric(asym, :U), :L)
                @test_throws ArgumentError Hermitian(Hermitian(aherm, :U), :L)
                # mixed cases with Hermitian/Symmetric
                if eltya <: Real
                    @test Symmetric(Hermitian(aherm, :U))     === Symmetric(aherm, :U)
                    @test Hermitian(Symmetric(asym, :U))     === Hermitian(asym, :U)
                    @test Symmetric(Hermitian(aherm, :U), :U) === Symmetric(aherm, :U)
                    @test Hermitian(Symmetric(asym, :U), :U) === Hermitian(asym, :U)
                    @test_throws ArgumentError Symmetric(Hermitian(aherm, :U), :L)
                    @test_throws ArgumentError Hermitian(Symmetric(aherm, :U), :L)
                end
            end
            @testset "similar" begin
                @test isa(similar(Symmetric(asym)), Symmetric{eltya})
                @test isa(similar(Hermitian(aherm)), Hermitian{eltya})
                @test isa(similar(Symmetric(asym), Int), Symmetric{Int})
                @test isa(similar(Hermitian(aherm), Int), Hermitian{Int})
                @test isa(similar(Symmetric(asym), (3,2)), Matrix{eltya})
                @test isa(similar(Hermitian(aherm), (3,2)), Matrix{eltya})
                @test isa(similar(Symmetric(asym), Int, (3,2)), Matrix{Int})
                @test isa(similar(Hermitian(aherm), Int, (3,2)), Matrix{Int})
            end

            @testset "Array/Matrix constructor from Symmetric/Hermitian" begin
                @test asym  == Matrix(Symmetric(asym))  == Array(Symmetric(asym))
                @test aherm == Matrix(Hermitian(aherm)) == Array(Hermitian(aherm))
            end

            @testset "parent" begin
                @test asym === parent(Symmetric(asym))
                @test aherm === parent(Hermitian(aherm))
            end
            # Unary minus for Symmetric/Hermitian matrices
            @testset "Unary minus for Symmetric/Hermitian matrices" begin
                @test (-Symmetric(asym))::typeof(Symmetric(asym)) == -asym
                @test (-Hermitian(aherm))::typeof(Hermitian(aherm)) == -aherm
            end

            @testset "Addition and subtraction for Symmetric/Hermitian matrices" begin
                for f in (+, -)
                    @test (f(Symmetric(asym), Symmetric(aposs)))::typeof(Symmetric(asym)) == f(asym, aposs)
                    @test (f(Hermitian(aherm), Hermitian(apos)))::typeof(Hermitian(aherm)) == f(aherm, apos)
                    @test (f(Symmetric(real(asym)), Hermitian(aherm)))::typeof(Hermitian(aherm)) == f(real(asym), aherm)
                    @test (f(Hermitian(aherm), Symmetric(real(asym))))::typeof(Hermitian(aherm)) == f(aherm, real(asym))
                    @test (f(Symmetric(asym), Hermitian(aherm))) == f(asym, aherm)
                    @test (f(Hermitian(aherm), Symmetric(asym))) == f(aherm, asym)
                end
            end

            @testset "getindex and unsafe_getindex" begin
                @test aherm[1,1] == Hermitian(aherm)[1,1]
                @test asym[1,1] == Symmetric(asym)[1,1]
                @test Symmetric(asym)[1:2,1:2] == asym[1:2,1:2]
                @test Hermitian(aherm)[1:2,1:2] == aherm[1:2,1:2]
            end

            @testset "conversion" begin
                @test Symmetric(asym) == convert(Symmetric,Symmetric(asym))
                if eltya <: Real
                    typs = [Float16,Float32,Float64]
                    for typ in typs
                        @test Symmetric(convert(Matrix{typ},asym)) == convert(Symmetric{typ,Matrix{typ}},Symmetric(asym))
                    end
                end
                if eltya <: Complex
                    typs = [ComplexF32,ComplexF64]
                    for typ in typs
                        @test Symmetric(convert(Matrix{typ},asym)) == convert(Symmetric{typ,Matrix{typ}},Symmetric(asym))
                        @test Hermitian(convert(Matrix{typ},aherm)) == convert(Hermitian{typ,Matrix{typ}},Hermitian(aherm))
                    end
                end
                @test Symmetric{eltya, Matrix{eltya}}(Symmetric(asym, :U)) === Symmetric(asym, :U)
                @test Hermitian{eltya, Matrix{eltya}}(Hermitian(aherm, :U)) === Hermitian(aherm, :U)
            end

            @testset "issymmetric, ishermitian" begin
                @test issymmetric(Symmetric(asym))
                @test ishermitian(Hermitian(aherm))
                if eltya <: Real
                    @test ishermitian(Symmetric(asym))
                    @test issymmetric(Hermitian(asym))
                elseif eltya <: Complex
                    # test that zero imaginary component is
                    # handled properly
                    @test ishermitian(Symmetric(b + b'))
                end
            end

            @testset "tril/triu" begin
                for (op, validks) in (
                        (triu, (-n + 1):(n + 1)),
                        (tril, (-n - 1):(n - 1)) )
                    for di in validks
                        @test op(Symmetric(asym), di) == op(asym, di)
                        @test op(Hermitian(aherm), di) == op(aherm, di)
                        @test op(Symmetric(asym, :L), di) == op(asym, di)
                        @test op(Hermitian(aherm, :L), di) == op(aherm, di)
                    end
                end
            end

            @testset "transpose, adjoint" begin
                S = Symmetric(asym)
                H = Hermitian(aherm)
                @test transpose(S) === S == asym
                @test adjoint(H) === H == aherm
                if eltya <: Real
                    @test adjoint(S) === S == asym
                    @test  transpose(H) === H == aherm
                else
                    @test adjoint(S) ==  Symmetric(conj(asym))
                    @test transpose(H) ==  Hermitian(copy(transpose(aherm)))
                end
                @test copy(adjoint(H)) == copy(aherm)
                @test copy(transpose(S)) == copy(asym)
            end

            @testset "real, imag" begin
                S = Symmetric(asym)
                H = Hermitian(aherm)
                @test issymmetric(real(S))
                @test ishermitian(real(H))
                if eltya <: Real
                    @test real(S) === S == asym
                    @test real(H) === H == aherm
                elseif eltya <: Complex
                    @test issymmetric(imag(S))
                    @test !ishermitian(imag(H))
                end
            end

        end

        @testset "linalg unary ops" begin
            @testset "tr" begin
                @test tr(asym) == tr(Symmetric(asym))
                @test tr(aherm) == tr(Hermitian(aherm))
            end

            @testset "isposdef[!]" begin
                @test isposdef(Symmetric(asym))  == isposdef(asym)
                @test isposdef(Symmetric(aposs)) == isposdef(aposs) == true
                @test isposdef(Hermitian(aherm)) == isposdef(aherm)
                @test isposdef(Hermitian(apos))  == isposdef(apos) == true
                if eltya != Int #chol! won't work with Int
                    @test isposdef!(Symmetric(copy(asym)))  == isposdef(asym)
                    @test isposdef!(Symmetric(copy(aposs))) == isposdef(aposs) == true
                    @test isposdef!(Hermitian(copy(aherm))) == isposdef(aherm)
                    @test isposdef!(Hermitian(copy(apos)))  == isposdef(apos) == true
                end
            end

            @testset "$f" for f in (det, logdet, logabsdet)
                for uplo in (:U, :L)
                    @test all(f(apos)  .≈ f(Hermitian(apos, uplo)))
                    @test all(f(aposs) .≈ f(Symmetric(aposs, uplo)))
                    if f != logdet
                        @test all(f(aherm) .≈ f(Hermitian(aherm, uplo)))
                        @test all(f(asym)  .≈ f(Symmetric(asym, uplo)))
                    end
                end
            end

            @testset "inversion" begin
                for uplo in (:U, :L)
                    @test inv(Symmetric(asym, uplo))::Symmetric ≈ inv(asym)
                    @test inv(Hermitian(aherm, uplo))::Hermitian ≈ inv(aherm)
                    @test inv(Symmetric(a, uplo))::Symmetric ≈ inv(Matrix(Symmetric(a, uplo)))
                    if eltya <: Real
                        @test inv(Hermitian(a, uplo))::Hermitian ≈ inv(Matrix(Hermitian(a, uplo)))
                    end
                end
                if eltya <: LinearAlgebra.BlasComplex
                    @testset "inverse edge case with complex Hermitian" begin
                        # Hermitian matrix, where inv(lu(A)) generates non-real diagonal elements
                        for T in (ComplexF32, ComplexF64)
                            A = T[0.650488+0.0im 0.826686+0.667447im; 0.826686-0.667447im 1.81707+0.0im]
                            H = Hermitian(A)
                            @test inv(H) ≈ inv(A)
                            @test ishermitian(Matrix(inv(H)))
                        end
                    end
                end
            end

            # Revisit when implemented in julia
            if eltya != BigFloat
                @testset "cond" begin
                    if eltya <: Real #svdvals! has no method for Symmetric{Complex}
                        @test cond(Symmetric(asym)) ≈ cond(asym)
                    end
                    @test cond(Hermitian(aherm)) ≈ cond(aherm)
                end

                @testset "symmetric eigendecomposition" begin
                    if eltya <: Real # the eigenvalues are only real and ordered for Hermitian matrices
                        d, v = eigen(asym)
                        @test asym*v[:,1] ≈ d[1]*v[:,1]
                        @test v*Diagonal(d)*transpose(v) ≈ asym
                        @test isequal(eigvals(asym[1]), eigvals(asym[1:1,1:1])[1])
                        @test abs.(eigen(Symmetric(asym), 1:2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2)
                        @test abs.(eigen(Symmetric(asym), d[1] - 1, (d[2] + d[3])/2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2)
                        @test eigvals(Symmetric(asym), 1:2) ≈ d[1:2]
                        @test eigvals(Symmetric(asym), d[1] - 1, (d[2] + d[3])/2) ≈ d[1:2]
                        # eigen doesn't support Symmetric{Complex}
                        @test Matrix(eigen(asym)) ≈ asym
                        @test eigvecs(Symmetric(asym)) ≈ eigvecs(asym)
                    end

                    d, v = eigen(aherm)
                    @test aherm*v[:,1] ≈ d[1]*v[:,1]
                    @test v*Diagonal(d)*v' ≈ aherm
                    @test isequal(eigvals(aherm[1]), eigvals(aherm[1:1,1:1])[1])
                    @test abs.(eigen(Hermitian(aherm), 1:2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2)
                    @test abs.(eigen(Hermitian(aherm), d[1] - 1, (d[2] + d[3])/2).vectors'v[:,1:2]) ≈ Matrix(I, 2, 2)
                    @test eigvals(Hermitian(aherm), 1:2) ≈ d[1:2]
                    @test eigvals(Hermitian(aherm), d[1] - 1, (d[2] + d[3])/2) ≈ d[1:2]
                    @test Matrix(eigen(aherm)) ≈ aherm
                    @test eigvecs(Hermitian(aherm)) ≈ eigvecs(aherm)

                    # relation to svdvals
                    if eltya <: Real #svdvals! has no method for Symmetric{Complex}
                        @test sum(sort(abs.(eigvals(Symmetric(asym))))) == sum(sort(svdvals(Symmetric(asym))))
                    end
                    @test sum(sort(abs.(eigvals(Hermitian(aherm))))) == sum(sort(svdvals(Hermitian(aherm))))
                end

                @testset "rank" begin
                    let A = a[:,1:5]*a[:,1:5]'
                        # Make sure A is Hermitian even in the presence of rounding error
                        # xianyi/OpenBLAS#729
                        A = (A + A') / 2
                        @test rank(A) == rank(Hermitian(A))
                    end
                end

                @testset "pow" begin
                    # Integer power
                    @test (asym)^2   ≈ (Symmetric(asym)^2)::Symmetric
                    @test (asym)^-2  ≈ (Symmetric(asym)^-2)::Symmetric
                    @test (aposs)^2  ≈ (Symmetric(aposs)^2)::Symmetric
                    @test (aherm)^2  ≈ (Hermitian(aherm)^2)::Hermitian
                    @test (aherm)^-2 ≈ (Hermitian(aherm)^-2)::Hermitian
                    @test (apos)^2   ≈ (Hermitian(apos)^2)::Hermitian
                    # integer floating point power
                    @test (asym)^2.0   ≈ (Symmetric(asym)^2.0)::Symmetric
                    @test (asym)^-2.0  ≈ (Symmetric(asym)^-2.0)::Symmetric
                    @test (aposs)^2.0  ≈ (Symmetric(aposs)^2.0)::Symmetric
                    @test (aherm)^2.0  ≈ (Hermitian(aherm)^2.0)::Hermitian
                    @test (aherm)^-2.0 ≈ (Hermitian(aherm)^-2.0)::Hermitian
                    @test (apos)^2.0   ≈ (Hermitian(apos)^2.0)::Hermitian
                    # non-integer floating point power
                    @test (asym)^2.5   ≈ (Symmetric(asym)^2.5)::Symmetric
                    @test (asym)^-2.5  ≈ (Symmetric(asym)^-2.5)::Symmetric
                    @test (aposs)^2.5  ≈ (Symmetric(aposs)^2.5)::Symmetric
                    @test (aherm)^2.5  ≈ (Hermitian(aherm)^2.5)#::Hermitian
                    @test (aherm)^-2.5 ≈ (Hermitian(aherm)^-2.5)#::Hermitian
                    @test (apos)^2.5   ≈ (Hermitian(apos)^2.5)::Hermitian
                end
            end
        end

        @testset "linalg binary ops" begin
            @testset "mat * vec" begin
                @test Symmetric(asym)*x+y ≈ asym*x+y
                # testing fallbacks for transpose-vector * transpose(SymHerm)
                xadj = transpose(x)
                @test xadj * transpose(Symmetric(asym)) ≈ xadj * asym
                @test x' * Symmetric(asym) ≈ x' * asym

                @test Hermitian(aherm)*x+y ≈ aherm*x+y
                # testing fallbacks for adjoint-vector * SymHerm'
                xadj = x'
                @test x' * Hermitian(aherm) ≈ x' * aherm
                @test xadj * Hermitian(aherm)' ≈ xadj * aherm
            end

            @testset "mat * mat" begin
                C = zeros(eltya,n,n)
                @test Hermitian(aherm) * a ≈ aherm * a
                @test a * Hermitian(aherm) ≈ a * aherm
                @test Hermitian(aherm) * Hermitian(aherm) ≈ aherm*aherm
                @test_throws DimensionMismatch Hermitian(aherm) * Vector{eltya}(undef, n+1)
                LinearAlgebra.mul!(C,a,Hermitian(aherm))
                @test C ≈ a*aherm

                @test Symmetric(asym) * Symmetric(asym) ≈ asym*asym
                @test Symmetric(asym) * a ≈ asym * a
                @test a * Symmetric(asym) ≈ a * asym
                @test_throws DimensionMismatch Symmetric(asym) * Vector{eltya}(undef, n+1)
                LinearAlgebra.mul!(C,a,Symmetric(asym))
                @test C ≈ a*asym

                tri_b = UpperTriangular(triu(b))
                @test Array(transpose(Hermitian(aherm)) * tri_b) ≈ transpose(aherm) * Array(tri_b)
                @test Array(tri_b * transpose(Hermitian(aherm))) ≈ Array(tri_b) * transpose(aherm)
                @test Array(Hermitian(aherm)' * tri_b) ≈ aherm' * Array(tri_b)
                @test Array(tri_b * Hermitian(aherm)') ≈ Array(tri_b) * aherm'

                @test Array(transpose(Symmetric(asym)) * tri_b) ≈ transpose(asym) * Array(tri_b)
                @test Array(tri_b * transpose(Symmetric(asym))) ≈ Array(tri_b) * transpose(asym)
                @test Array(Symmetric(asym)' * tri_b) ≈ asym' * Array(tri_b)
                @test Array(tri_b * Symmetric(asym)') ≈ Array(tri_b) * asym'
            end
            @testset "solver" begin
                @test Hermitian(aherm)\x ≈ aherm\x
                @test Hermitian(aherm)\b ≈ aherm\b
                @test Symmetric(asym)\x  ≈ asym\x
                @test Symmetric(asym)\b  ≈ asym\b
            end
        end
    end
end

#Issue #7647: test xsyevr, xheevr, xstevr drivers.
@testset "Eigenvalues in interval for $(typeof(Mi7647))" for Mi7647 in
        (Symmetric(diagm(0 => 1.0:3.0)),
         Hermitian(diagm(0 => 1.0:3.0)),
         Hermitian(diagm(0 => complex(1.0:3.0))),
         SymTridiagonal([1.0:3.0;], zeros(2)))
    @test eigmin(Mi7647)  == eigvals(Mi7647, 0.5, 1.5)[1] == 1.0
    @test eigmax(Mi7647)  == eigvals(Mi7647, 2.5, 3.5)[1] == 3.0
    @test eigvals(Mi7647) == eigvals(Mi7647, 0.5, 3.5) == [1.0:3.0;]
end

@testset "Hermitian wrapper ignores imaginary parts on diagonal" begin
    A = [1.0+im 2.0; 2.0 0.0]
    @test !ishermitian(A)
    @test Hermitian(A)[1,1] == 1
end

@testset "Issue #7933" begin
    A7933 = [1 2; 3 4]
    B7933 = copy(A7933)
    C7933 = Matrix(Symmetric(A7933))
    @test A7933 == B7933
end

@testset "Issues #8057 and #8058. f=$f, A=$A" for f in
        (eigen, eigvals),
            A in (Symmetric([0 1; 1 0]), Hermitian([0 im; -im 0]))
    @test_throws ArgumentError f(A, 3, 2)
    @test_throws ArgumentError f(A, 1:4)
end

@testset "Ignore imaginary part of Hermitian diagonal" begin
    A = [1.0+im 2.0; 2.0 0.0]
    @test !ishermitian(A)
    @test diag(Hermitian(A)) == real(diag(A))
end

@testset "Issue #17780" begin
    a = randn(2,2)
    a = a'a
    b = complex.(a,a)
    c = Symmetric(b)
    @test conj(c) == conj(Array(c))
    cc = copy(c)
    @test conj!(c) == conj(Array(cc))
    c = Hermitian(b + b')
    @test conj(c) == conj(Array(c))
    cc = copy(c)
    @test conj!(c) == conj(Array(cc))
end

@testset "Issue # 19225" begin
    X = [1 -1; -1 1]
    for T in (Symmetric, Hermitian)
        Y = T(copy(X))
        _Y = similar(Y)
        copyto!(_Y, Y)
        @test _Y == Y

        W = T(copy(X), :L)
        copyto!(W, Y)
        @test W.data == Y.data
        @test W.uplo != Y.uplo

        W[1,1] = 4
        @test W == T([4 -1; -1 1])
        @test_throws ArgumentError (W[1,2] = 2)

        @test Y + I == T([2 -1; -1 2])
        @test Y - I == T([0 -1; -1 0])
        @test Y * I == Y

        @test Y .+ 1 == T([2 0; 0 2])
        @test Y .- 1 == T([0 -2; -2 0])
        @test Y * 2 == T([2 -2; -2 2])
        @test Y / 1 == Y

        @test T([true false; false true]) .+ true == T([2 1; 1 2])
    end
end

@testset "Issue #21981" begin
    B = complex(rand(4,4))
    B[4,1] += 1im;
    @test ishermitian(Symmetric(B, :U))
    @test issymmetric(Hermitian(B, :U))
    B[4,1]  = real(B[4,1])
    B[1,4] += 1im
    @test ishermitian(Symmetric(B, :L))
    @test issymmetric(Hermitian(B, :L))
end

@testset "$HS solver with $RHS RHS - $T" for HS in (Hermitian, Symmetric),
        RHS in (Hermitian, Symmetric, Diagonal, UpperTriangular, LowerTriangular),
        T   in (Float64, ComplexF64)
    D = rand(T, 10, 10); D = D'D
    A = HS(D)
    B = RHS(D)
    @test A\B ≈ Matrix(A)\Matrix(B)
end

@testset "inversion of Hilbert matrix" begin
    for T in (Float64, ComplexF64)
        H = T[1/(i + j - 1) for i in 1:8, j in 1:8]
        @test norm(inv(Symmetric(H))*(H*fill(1., 8)) .- 1) ≈ 0 atol = 1e-5
        @test norm(inv(Hermitian(H))*(H*fill(1., 8)) .- 1) ≈ 0 atol = 1e-5
    end
end

@testset "similar should preserve underlying storage type and uplo flag" begin
    m, n = 4, 3
    sparsemat = sprand(m, m, 0.5)
    for SymType in (Symmetric, Hermitian)
        symsparsemat = SymType(sparsemat)
        @test isa(similar(symsparsemat), typeof(symsparsemat))
        @test similar(symsparsemat).uplo == symsparsemat.uplo
        @test isa(similar(symsparsemat, Float32), SymType{Float32,<:SparseMatrixCSC{Float32}})
        @test similar(symsparsemat, Float32).uplo == symsparsemat.uplo
        @test isa(similar(symsparsemat, (n, n)), typeof(sparsemat))
        @test isa(similar(symsparsemat, Float32, (n, n)), SparseMatrixCSC{Float32})
    end
end

@testset "#24572: eltype(A::HermOrSym) === eltype(parent(A))" begin
    A = rand(Float32, 3, 3)
    @test_throws TypeError Symmetric{Float64,Matrix{Float32}}(A, 'U')
    @test_throws TypeError Hermitian{Float64,Matrix{Float32}}(A, 'U')
end

@testset "fill[stored]!" begin
    for uplo in (:U, :L)
        # Hermitian
        A = Hermitian(fill(1.0+0im, 2, 2), uplo)
        @test fill!(A, 2) == fill(2, 2, 2)
        @test A.data == (uplo == :U ? [2 2; 1.0+0im 2] : [2 1.0+0im; 2 2])
        @test_throws ArgumentError fill!(A, 2+im)

        # Symmetric
        A = Symmetric(fill(1.0+im, 2, 2), uplo)
        @test fill!(A, 2) == fill(2, 2, 2)
        @test A.data == (uplo == :U ? [2 2; 1.0+im 2] : [2 1.0+im; 2 2])
    end
end

@testset "#25625 recursive transposition" begin
    A = Matrix{Matrix{Int}}(undef, 2, 2)
    A[1,1] = [1 2; 2 3]
    A[1,2] = [4 5 6; 7 8 9]
    A[2,1] = [4 7; 5 8; 6 9]
    A[2,2] = [1 2; 3 4]
    for uplo in (:U, :L)
        S = Symmetric(A, uplo)
        @test S[1,1] == A[1,1]
        @test S[1,2] == transpose(S[2,1]) == A[1,2]
        @test S[2,2] == Symmetric(A[2,2], uplo)
        @test S == transpose(S) == Matrix(S) == Matrix(transpose(S)) == transpose(Matrix(S))
    end

    B = Matrix{Matrix{Complex{Int}}}(undef, 2, 2)
    B[1,1] = [1 2+im; 2-im 3]
    B[1,2] = [4 5+1im 6-2im; 7+3im 8-4im 9+5im]
    B[2,1] = [4 7-3im; 5-1im 8+4im; 6+2im 9-5im]
    B[2,2] = [1+1im 2+2im; 3-3im 4-2im]
    for uplo in (:U, :L)
        H = Hermitian(B, uplo)
        @test H[1,1] == Hermitian(B[1,1], uplo)
        @test H[1,2] == adjoint(H[2,1]) == B[1,2]
        @test H[2,1] == adjoint(H[1,2]) == B[2,1]
        @test H[2,2] == Hermitian(B[2,2], uplo)
        @test H == adjoint(H) == Matrix(H) == Matrix(adjoint(H)) == adjoint(Matrix(H))
    end
end

@testset "getindex of diagonal element (#25972)" begin
    A = rand(ComplexF64, 2, 2)
    @test Hermitian(A, :U)[1,1] == Hermitian(A, :L)[1,1] == real(A[1,1])
end

@testset "issue #29392: SymOrHerm scaled with Number" begin
    R = rand(Float64, 2, 2); C = rand(ComplexF64, 2, 2)
    # Symmetric * Real, Real * Symmetric
    A = Symmetric(R); x = 2.0
    @test (A * x)::Symmetric == (x * A)::Symmetric
    A = Symmetric(C); x = 2.0
    @test (A * x)::Symmetric == (x * A)::Symmetric
    # Symmetric * Complex, Complex * Symmetrics
    A = Symmetric(R); x = 2.0im
    @test (A * x)::Symmetric == (x * A)::Symmetric
    A = Symmetric(C); x = 2.0im
    @test (A * x)::Symmetric == (x * A)::Symmetric
    # Hermitian * Real, Real * Hermitian
    A = Hermitian(R); x = 2.0
    @test (A * x)::Hermitian == (x * A)::Hermitian
    A = Hermitian(C); x = 2.0
    @test (A * x)::Hermitian == (x * A)::Hermitian
    # Hermitian * Complex, Complex * Hermitian
    A = Hermitian(R); x = 2.0im
    @test (A * x)::Matrix == (x * A)::Matrix
    A = Hermitian(C); x = 2.0im
    @test (A * x)::Matrix == (x * A)::Matrix
    # Symmetric / Real
    A = Symmetric(R); x = 2.0
    @test (A / x)::Symmetric == Matrix(A) / x
    A = Symmetric(C); x = 2.0
    @test (A / x)::Symmetric == Matrix(A) / x
    # Symmetric / Complex
    A = Symmetric(R); x = 2.0im
    @test (A / x)::Symmetric == Matrix(A) / x
    A = Symmetric(C); x = 2.0im
    @test (A / x)::Symmetric == Matrix(A) / x
    # Hermitian / Real
    A = Hermitian(R); x = 2.0
    @test (A / x)::Hermitian == Matrix(A) / x
    A = Hermitian(C); x = 2.0
    @test (A / x)::Hermitian == Matrix(A) / x
    # Hermitian / Complex
    A = Hermitian(R); x = 2.0im
    @test (A / x)::Matrix == Matrix(A) / x
    A = Hermitian(C); x = 2.0im
    @test (A / x)::Matrix == Matrix(A) / x
end

@testset "issue #30814: Symmetric of Hermitian if diag is not real" begin
    A = [1 2; 3 4] * (1 + im)
    B = Hermitian(A)
    @test_throws ArgumentError Symmetric(B) == Symmetric(Matrix(B))
    A[1,1] = 1; A[2,2] = 4
    @test Symmetric(B) == Symmetric(Matrix(B))
end

@testset "symmetric()/hermitian() for Numbers" begin
    @test LinearAlgebra.symmetric(1, :U) == 1
    @test LinearAlgebra.symmetric_type(Int) == Int
    @test LinearAlgebra.hermitian(1, :U) == 1
    @test LinearAlgebra.hermitian_type(Int) == Int
end

end # module TestSymmetric

max_stars_repo_name: musm/julia

max_issues_repo_name: musm/julia

max_forks_repo_name: musm/julia


text:
sign(a::Float120) = sign(a.hi)
@inline signbit(a::Float120) = signbit(a.hi)

#=
    negation, abs, flipsign, copysign
=#

@inline (-)(a::Float120) = Float120(-a.hi, -a.lo)
@inline abs(a::Float120) = ifelse(signbit(a.hi), -a, a)

flipsign(a::Float120, b::Float120)   = ifelse(signbit(b.hi), -a, a)
flipsign{T<:Real}(a::Float120, b::T) = ifelse(signbit(b), -a, a)
copysign(a::Float120, b::Float120)   = ifelse(signbit(b.hi), -abs(a), abs(a))
copysign{T<:Real}(a::Float120, b::T) = ifelse(signbit(b), -abs(a), abs(a))

#=
    frexp, ldexp, zero, one
=#

function frexp(a::Float120)
    frhi, xphi = frexp(a.hi)
    frlo, xplo = frexp(a.lo)
    Float120(frhi, ldexp(frlo,xplo-xphi)), xphi
end

function ldexp(a::Float120,xp::Int)
    Float120(ldexp(a.hi,xp),ldexp(a.lo,xp))
end
ldexp{I<:Integer}(fx::Tuple{Float120,I}) = ldexp(fx...)

zero(::Type{Float120}) = Float120(0.0, 0.0)
one(::Type{Float120})  = Float120(1.0, 0.0)
zero(x::Float120) = zero(Float120)
one(x::Float120)  = one(Float120)

# Inf, NaN

Inf30(::Type{Float120}) = Float120(Inf(T), Inf(T))
NaN30(::Type{Float120}) = Float120(NaN(T), NaN(T))
Inf30(x::Float120) = Inf30(Float120)
NaN30(x::Float120) = NaN30(Float120)

ext: jl


text:
""" Hanoi water network model. """
function hanoi_network(oct::Bool = false, friction::String = "HW")
    m = JuMP.Model()
    N = 32
    m[:N] = 32 # Number of nodes
    K = 34
    m[:K] = 34

    sink = [0, 890, 850, 130, 725, 1005, 1350, 550, 525,
            525, 500, 560, 940, 615, 280, 310,
            865, 1345, 60, 1275, 930, 485, 1045, 820, 170,
            900, 370, 290, 360, 360, 105, 805]/3600.0
    source = zeros(m[:N])
    source[1] = sum(sink)
    topology_list = [[0, 1], [1, 2], [2, 3], [2, 19], [2, 18], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10],
                    [10, 11], [11, 12], [9, 13], [13, 14], [14, 15], [16, 15], [17, 16], [18, 17], [15, 26], [26, 25],
                    [25, 24], [23, 24], [22, 23], [22, 27], [27, 28], [28, 29], [29, 30], [24, 31], [31, 30], [19, 20],
                    [20, 21], [19, 22]]
    topology_dict = Dict(i => topology_list[i] .+ 1 for i = 1:m[:K])
    coordinates = Dict(1 => (5021.20, 1582.17), 2 => (5025.20, 2585.42), 3 => (5874.22, 2588.30), 4 => (6873.11, 2588.30),
                5 => (8103.51, 2585.42), 6 => (8103.51, 3234.67), 7 => (8106.66, 4179.28), 8 => (8106.66, 5133.78),
                9 => (7318.64, 5133.78), 10 => (7319.94, 5831.65), 11 => (7319.94, 6671.19), 12 => (5636.76, 6676.24),
                13 => (6530.63, 5133.78), 14 => (5676.02, 5133.78), 15 => (5021.20, 5133.78), 16 => (5021.20, 4412.36),
                17 => (5021.20, 3868.52), 18 => (5021.20, 3191.49), 19 => (3587.87, 2588.30), 20 => (3587.87, 1300.84),
                21 => (3587.87, 901.29), 22 => (1978.55, 2588.30), 23 => (1975.58, 4084.35), 24 => (1980.46, 5137.63),
                25 => (3077.46, 5137.63), 26 => (3933.52, 5133.78), 27 => (846.04, 2588.20), 28 => (-552.41, 2588.20),
                29 => (-552.38, 4369.06), 30 => (-549.36, 5137.63), 31 => (536.45, 5137.63), 0 => (5360.71, 1354.05))
    L = [100, 1350, 900, 1150, 1450, 450, 850, 850, 800, 950, 1200, 3500, 800, 500, 550, 2730, 1750, 800, 400, 2200,
        1500, 500, 2650, 1230, 1300, 850, 300, 750, 1500, 2000, 1600, 150, 860, 950]

    # Parameters
    rho = 1000 # Density (kg/m^3)
    mu =  8.9e-4 # Viscosity, kg/m/s
    g = 9.81 # Acceleration of gravity, m/s^2
    H_min = 30
    D_min = 0.3048 # Minimum pipe diameter, m
    D_max = 1.016 # Maximum pipe diameter, m
    V_max = 100

    if friction == "HW"
        rough = 130 .* ones(K)
    elseif friction == "DW"
        rough = 0.26e-6 .* ones(K)
    end

    @variable(m, H[1:N] >= H_min)             # Head
    @constraint(m, H[1] == 100)
    @variable(m, pipeCost[1:K] >= 0)      # Pipe cost
    @variable(m, D_max >= D[1:K] >= D_min)             # Pipe diameter
    @variable(m, 6 >= flow[1:K] >= 0)
    # @variable(m, flow_dir[1:K], Bin)      # Flow direction
    @variable(m, H_loss[1:K] >= 0)        # Head loss
    @variable(m, V_max >= V[1:K] >= 0)             # Flow velocity (m/s)

    # Flow directionality
    flow_dir = ones(K)

    if friction == "DW"
        @variable(m, relRough[1:K] >= 0)    # Relative pipe roughness
        @variable(m, Re[1:K] >= 0)          # Reynolds number
        @variable(m, 1 >= f[1:K] >= 0)           # Friction factor
    end

    @objective(m, Min, sum(pipeCost))

    for i = 1:N
        @constraint(m, source[i] - sink[i] +
            sum(flow[pipe_index] for (pipe_index, node) in topology_dict if node[2] == i) -
            sum(flow[pipe_index] for (pipe_index, node) in topology_dict if node[1] == i) >= 0)
    end

    for (pipe_index, node) in topology_dict
        @constraint(m, H_loss[pipe_index] <= H[node[1]] - H[node[2]]) # TODO: add directionality
    end

    # Now for the nonlinear constraints
    if oct
        gm = GlobalModel(model = m, name = "Hanoi")
        @variable(gm.model, Ls[1:K])
        set_lower_bound.(Ls, L)
        set_upper_bound.(Ls, L)
        for i = 1:K
            add_nonlinear_constraint(gm, :((D,Ls) -> 1.1 * D ^ 1.5 * Ls), vars = [D[i], Ls[i]], dependent_var = pipeCost[i])
            add_nonlinear_constraint(gm, :((flow, D) -> 4 * flow / (pi*D^2)), vars = [flow[i], D[i]], dependent_var = V[i], equality = true)
        end

        for pipe_index = 1:K
            if friction == "HW"
                add_nonlinear_constraint(gm, :((flow, D, Ls) -> $(flow_dir[pipe_index]) * 10.67 * Ls * flow ^ 1.852 /
                ($(rough[pipe_index])^1.852*D^4.8704)),
                vars = [flow[pipe_index], D[pipe_index], Ls[pipe_index]], dependent_var = H_loss[pipe_index])
            elseif friction == "DW"
                throw(OCTHaGOnException("DW water network formulation incomplete. "))
                # @NLconstraint(m, H_loss[pipe_index] == f[pipe_index] * L[pipe_index] * V[pipe_index] ^ 2 / (
                #                         2 * D[pipe_index] * g))
                # @NLconstraint(m, relRough[pipe_index] == rough[pipe_index] / D[pipe_index])
                # @NLconstraint(m, Re[pipe_index] == rho * V[pipe_index] * D[pipe_index] / mu)
                # # From frictionFactorFitting.py
                # @NLconstraint(m, f[pipe_index] ^ 2.39794 >= 3.26853e-06 * Re[pipe_index] ^ 0.0574443 *
                #             relRough[pipe_index] ^ 0.364794 + 0.0001773 * Re[pipe_index] ^ -0.529499 *
                #             relRough[pipe_index] ^ -0.0810121 + 0.00301918 * Re[pipe_index] ^ -0.0220498 *
                #             relRough[pipe_index] ^ 1.73526 + 0.0734922 * Re[pipe_index] ^ -1.13629 *
                #             relRough[pipe_index] ^ 0.0574655 + 0.000214297 * Re[pipe_index] ^ 0.00035242 *
                #             relRough[pipe_index] ^ 0.823896)
            end
        end
        return gm
    else
        for i = 1:K
            @NLconstraint(m, pipeCost[i] >= 1.1*D[i]^1.5*L[i])
            @NLconstraint(m, V[i] == 4 * flow[i] / (pi*D[i]^2))
        end

        for pipe_index = 1:K
            if friction == "HW"
                @NLconstraint(m, H_loss[pipe_index] == flow_dir[pipe_index] * 10.67 * L[pipe_index] * flow[pipe_index] ^ 1.852 /
                                                    (rough[pipe_index]^1.852*D[pipe_index]^4.8704))
                            # S (hydraulic slope H/L) = 10.67*Q^1.852 (volumetric flow rate) /
                            # C^1.852 (pipe roughness) /d^4.8704 (pipe diameter)
            elseif friction == "DW"
                @NLconstraint(m, H_loss[pipe_index] == f[pipe_index] * L[pipe_index] * V[pipe_index] ^ 2 / (
                                        2 * D[pipe_index] * g))
                @NLconstraint(m, relRough[pipe_index] == rough[pipe_index] / D[pipe_index])
                @NLconstraint(m, Re[pipe_index] == rho * V[pipe_index] * D[pipe_index] / mu)
                # From frictionFactorFitting.py
                @NLconstraint(m, f[pipe_index] ^ 2.39794 >= 3.26853e-06 * Re[pipe_index] ^ 0.0574443 *
                            relRough[pipe_index] ^ 0.364794 + 0.0001773 * Re[pipe_index] ^ -0.529499 *
                            relRough[pipe_index] ^ -0.0810121 + 0.00301918 * Re[pipe_index] ^ -0.0220498 *
                            relRough[pipe_index] ^ 1.73526 + 0.0734922 * Re[pipe_index] ^ -1.13629 *
                            relRough[pipe_index] ^ 0.0574655 + 0.000214297 * Re[pipe_index] ^ 0.00035242 *
                            relRough[pipe_index] ^ 0.823896)
            end
        end
        return m
    end
end

# # Plotting results
# using LightGraphs
# using GraphPlot

# G = Graph(N)
# for i = 1:K
#     add_edge!(G, topology_dict[1][1], topology_dict[1][2])
# end

# gplot(G, nodelabel = 1:N)


text:
module TestConvert
using Test
using CategoricalArrays
using CategoricalArrays: DefaultRefType, level, reftype, leveltype

@testset "convert() for CategoricalPool{Int, DefaultRefType} and values" begin
    pool = CategoricalPool([1, 2, 3])
    @test convert(CategoricalPool{Int, DefaultRefType}, pool) === pool
    @test convert(CategoricalPool{Int}, pool) === pool
    @test convert(CategoricalPool, pool) === pool
    convert(CategoricalPool{Float64, UInt8}, pool)
    convert(CategoricalPool{Float64}, pool)
    convert(CategoricalPool, pool)

    v1 = CategoricalValue(1, pool)
    v2 = CategoricalValue(2, pool)
    v3 = CategoricalValue(3, pool)
    @test eltype(v1) === Any
    @test eltype(typeof(v1)) === Any
    @test leveltype(v1) === Int
    @test leveltype(typeof(v1)) === Int
    @test reftype(v1) === DefaultRefType
    @test reftype(typeof(v1)) === DefaultRefType
    @test v1 isa CategoricalArrays.CategoricalValue{Int, DefaultRefType}

    @test convert(Int, v1) === 1
    @test convert(Int, v2) === 2
    @test convert(Int, v3) === 3

    @test convert(Int32, v1) === Int32(1)
    @test convert(Int32, v2) === Int32(2)
    @test convert(Int32, v3) === Int32(3)

    @test convert(UInt8, v1) === 0x01
    @test convert(UInt8, v2) === 0x02
    @test convert(UInt8, v3) === 0x03

    @test convert(CategoricalValue, v1) === v1
    @test convert(CategoricalValue{Int}, v1) === v1
    @test convert(CategoricalValue{Int, DefaultRefType}, v1) === v1

    @test convert(Any, v1) === v1
    @test convert(Any, v2) === v2
    @test convert(Any, v3) === v3

    for T in (typeof(v1), CategoricalValue{Int}, CategoricalValue), U in (Missing, Nothing)
        @test convert(Union{T, U}, v1) === v1
        @test convert(Union{T, U}, v2) === v2
        @test convert(Union{T, U}, v3) === v3
    end

    @test get(v1) === 1
    @test get(v2) === 2
    @test get(v3) === 3

    @test promote(1, v1) === (1, 1)
    @test promote(1.0, v1) === (1.0, 1.0)
    @test promote(0x1, v1) === (1, 1)
end

@testset "promote_type" begin
    @test promote_type(CategoricalValue{Int}, CategoricalValue{Float64}) ===
        CategoricalValue{Float64}
    @test promote_type(CategoricalValue{Int, UInt8}, CategoricalValue{Float64, UInt32}) ===
        CategoricalValue{Float64, UInt32}
    @test promote_type(CategoricalValue{Int, UInt8}, CategoricalValue{Float64}) ===
        CategoricalValue{Float64}
    # Tests that return Any before Julia 1.3 are due to JuliaLang/julia#29348
    if VERSION >= v"1.3.0-DEV"
        @test promote_type(CategoricalValue{Int},
                           Union{CategoricalValue{Float64}, Missing}) ===
            Union{Missing, Float64}
        @test promote_type(CategoricalValue{Int, UInt8},
                           Union{CategoricalValue{Float64, UInt32}, Missing}) ===
            Union{CategoricalValue{Float64, UInt32}, Missing}
        @test promote_type(Union{CategoricalValue{Int}, Missing},
                           CategoricalValue{Float64}) ===
            Union{Missing, Float64}
        @test promote_type(Union{CategoricalValue{Int, UInt8}, Missing},
                           CategoricalValue{Float64, UInt32}) ===
            Union{CategoricalValue{Float64, UInt32}, Missing}
        @test promote_type(Union{CategoricalValue{Int}, Missing},
                           Union{CategoricalValue{Float64}, Missing}) ===
            Union{Missing, Float64}
    else
        @test promote_type(CategoricalValue{Int},
                           Union{CategoricalValue{Float64}, Missing}) ===
            Any
        @test promote_type(CategoricalValue{Int, UInt8},
                           Union{CategoricalValue{Float64, UInt32}, Missing}) ===
            Union{CategoricalValue{Float64, UInt32}, Missing}
        @test promote_type(Union{CategoricalValue{Int}, Missing},
                           CategoricalValue{Float64}) ===
            Any
        @test promote_type(Union{CategoricalValue{Int, UInt8}, Missing},
                        CategoricalValue{Float64, UInt32}) ===
            Union{CategoricalValue{Float64, UInt32}, Missing}
        @test promote_type(Union{CategoricalValue{Int}, Missing},
                        Union{CategoricalValue{Float64}, Missing}) ===
            Any
    end
    @test promote_type(Union{CategoricalValue{Int, UInt8}, Missing},
                       Union{CategoricalValue{Float64, UInt32}, Missing}) ===
        Union{CategoricalValue{Float64, UInt32}, Missing}

    @test promote_type(CategoricalValue, Missing) === Union{CategoricalValue, Missing}
    @test promote_type(CategoricalValue{Int}, Missing) === Union{CategoricalValue{Int}, Missing}
    @test promote_type(CategoricalValue{Int, UInt32}, Missing) ===
        Union{CategoricalValue{Int, UInt32}, Missing}
    @test promote_type(CategoricalValue{Int, UInt32}, Any) === Any
end

@testset "convert() preserves `ordered`" begin
    pool = CategoricalPool([1, 2, 3], true)
    @test convert(CategoricalPool{Float64, UInt8}, pool).ordered === true
end

@testset "convert() with Union{T, Nothing}" begin
    pool = CategoricalPool([nothing, 2, 3])
    v1 = CategoricalValue(1, pool)
    v2 = CategoricalValue(2, pool)
    @test convert(Union{Int, Nothing}, v1) === nothing
    @test convert(Union{Int, Nothing}, v2) === 2
    @test convert(Union{Float64, Nothing}, v2) === 2.0
end

@testset "levelcode" begin
    pool = CategoricalPool{Int,UInt8}([3, 1, 2])
    levels!(pool, [2, 1, 3])
    for i in 1:3
        v = CategoricalValue(i, pool)
        @test levelcode(v) isa Int16
        @test levels(pool)[levelcode(v)] == v
    end

    @test levelcode(missing) === missing
end

end


text:
# Autogenerated wrapper script for CFITSIO_jll for i686-w64-mingw32
export libcfitsio

using LibCURL_jll
JLLWrappers.@generate_wrapper_header("CFITSIO")
JLLWrappers.@declare_library_product(libcfitsio, "libcfitsio.dll")
function __init__()
    JLLWrappers.@generate_init_header(LibCURL_jll)
    JLLWrappers.@init_library_product(
        libcfitsio,
        "bin\\libcfitsio.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "GMM" begin
    Random.seed!(28374544);

    # generate dataset
    d = 2
    theta = [0,3,5,7] * (pi/4)
    k = length(theta)
    X = hcat([MultivariateMixtures.clusters(t, 250, d) for t in theta]...)

    G1 = fit_mle(MvNormal, X[:,1:250])
    GMM = fit_mm(FullNormal, X, k, tol=1e-10)
    @test Distributions.component_type(GMM) <: FullNormal
    @test ncomponents(GMM) == k
    @test probs(GMM) ≈ fill(1/4, k) atol=1e-2
    G2 = component(GMM,2)
    @test mean(G1) ≈ mean(G2) atol=1e-2
    @test cov(G1) ≈ cov(G2) atol=1e-2
    @test length(GMM) == d
    @test size(GMM) == (d,)
    @test mean(GMM) ≈ mean(X, dims=2) atol=1e-2
    @test cov(GMM) ≈ cov(X, dims=2) atol=1e-1

    GMM = fit_mm(FullNormal, X, k, tol=1e-5, homoscedastic=true)
    G3 = component(GMM,1)
    @test mean(G1) ≈ mean(G3) atol=1e-2
    @test cov(G1) ≈ cov(G3) atol=0.5
    @test cov(G3) == cov(component(GMM,1))

    GMM = fit_mm(DiagNormal, X, k, tol=1e-4)
    @test Distributions.component_type(GMM) <: DiagNormal
    @test ncomponents(GMM) == k
    @test probs(GMM) ≈ fill(1/4, k) atol=1e-2

    GMM = fit_mm(IsoNormal, X, k, tol=1e-4)
    @test Distributions.component_type(GMM) <: IsoNormal

    X = hcat([MultivariateMixtures.clusters(t, 250, 5) for t in theta]...)
    GMM = fit_mm(FullNormal, X, k, tol=1e-3)
    @test length(GMM) == 5

end


text:
module ExpmV

  # package code goes here
  include("expmv_fun.jl")
  include("degree_selector.jl")
  include("normAm.jl")
  include("expmv_tspan.jl")
  include("select_taylor_degree.jl")
  include("norm1est.jl")

end # module


text:
using BlackBoxOptim: general_stop_reason

@testset "general_stop_reason" begin
    problem = BlackBoxOptim.example_problems["Sphere"]

    @testset "Within fitness tolerance" begin
        # Set a low MinDeltaFitnessTolerance so it does not interfere.
        res = bboptimize(problem; NumDimensions = 1, TraceMode = :silent,
            FitnessTolerance = 1e-1, MinDeltaFitnessTolerance = 1e-100)
        @test general_stop_reason(res) == "Within fitness tolerance of optimum"
    end

    @testset "Delta fitness below tolerance" begin
        # Set a low FitnessTolerance so it does not interfere.
        res = bboptimize(problem; NumDimensions = 1, TraceMode = :silent,
            MinDeltaFitnessTolerance = 1e-1, FitnessTolerance = 1e-100)
        @test general_stop_reason(res) == "Delta fitness below tolerance"
    end
end


text:
module FerriteViz

using Makie
using Tensors
import Ferrite
import LinearAlgebra

abstract type AbstractPlotter end

include("utils.jl")
include("makieplotting.jl")

export MakiePlotter
export ferriteviewer

end


text:
using Test
using SeisSolXDMF
using Pkg.Artifacts

@testset "XDMF" begin
    test_resources_path = artifact"test_resources"

    filename = joinpath(test_resources_path, "output-surface-h5.xdmf")
    xdmf = XDMFFile(filename)

    @testset "Header" begin
        timesteps = timesteps_of(xdmf)
        @test length(timesteps) == 3
        @test timesteps == [float(i) / 10 for i ∈ 0:2]
        @test eltype(timesteps) == Float64
    end

    @testset "HDF5 - Surface" begin
        simplices, points = grid_of(xdmf)
        @test size(simplices) == (3, 172616)
        @test size(points) == (3, 86468)

        u = data_of(xdmf, 1, "u")
        @test ndims(u) == 1
        @test length(u) == size(simplices)[2]

        u = data_of(xdmf, 3, "u")
        @test minimum(u) < 0.
        @test maximum(u) > 0.
    end

    filename = joinpath(test_resources_path, "output-volume-bin.xdmf")
    xdmf = XDMFFile(filename)

    @testset "POSIX - Volume" begin
        simplices, points = grid_of(xdmf)
        @test size(simplices) == (4, 10000)
        @test size(points) == (3, 8000)

        u = data_of(xdmf, 1, "u")
        @test ndims(u) == 1
        @test length(u) == size(simplices)[2]
        @test all(tup -> begin (i, val) = tup; val == i * .001 end, enumerate(u))

        u = data_of(xdmf, 3, "u")
        @test all(tup -> begin (i, val) = tup; val ≈ 200 + i * .001 end, enumerate(u))
    end
end

max_stars_repo_name: mtronicsDev/SeisSolXDMF

max_issues_repo_name: mtronicsDev/SeisSolXDMF

max_forks_repo_name: mtronicsDev/SeisSolXDMF


text:
function build_identity_matrix2(::Type{HansenUtilities})
    return jcall(HansenUtilities, "buildIdentityMatrix2", PolynomialFunctionMatrix, ())
end

function build_identity_matrix4(::Type{HansenUtilities})
    return jcall(HansenUtilities, "buildIdentityMatrix4", PolynomialFunctionMatrix, ())
end

function build_zero_matrix2(::Type{HansenUtilities})
    return jcall(HansenUtilities, "buildZeroMatrix2", PolynomialFunctionMatrix, ())
end

function build_zero_matrix4(::Type{HansenUtilities})
    return jcall(HansenUtilities, "buildZeroMatrix4", PolynomialFunctionMatrix, ())
end


text:
@discord_object struct ChannelMention
    id::Snowflake
    guild_id::Snowflake
    type::ChannelType.ChannelTypeEnum
    name::String
end


text:
module val_

export
	val

include("../src/libiup_h.jl")

using IUP

#------------------------example html/examples/C/val.c ----------------------------------
function val()
	global lbl_v, lbl_h         # In a real case we would probably like to avoid the globals
	IupOpen()                   # Initializes IUP

	val_v = IupVal("VERTICAL")
	val_h = IupVal("HORIZONTAL")
	lbl_v = IupLabel("VALUE=")
	lbl_h = IupLabel("VALUE=")
	IupSetAttribute(lbl_v, "SIZE", "50x")
	IupSetAttribute(lbl_h, "SIZE", "50x")
	IupSetAttribute(val_v, "SHOWTICKS", "5")

	dlg_val = IupDialog(
		IupHbox(
			IupSetAttributes(IupHbox(val_v, lbl_v), "ALIGNMENT=ACENTER"),
			IupSetAttributes(IupVbox(val_h, lbl_h), "ALIGNMENT=ACENTER")
		)
	)

	IupSetCallback(val_v, "BUTTON_PRESS_CB",   cfunction(button_press, Int, (Ptr{Ihandle}, Float64)))
	IupSetCallback(val_v, "BUTTON_RELEASE_CB", cfunction(button_release, Int, (Ptr{Ihandle}, Float64)))
	IupSetCallback(val_v, "MOUSEMOVE_CB",      cfunction(mousemove, Int, (Ptr{Ihandle}, Float64)))

	IupSetCallback(val_h, "BUTTON_PRESS_CB",   cfunction(button_press, Int, (Ptr{Ihandle}, Float64)))
	IupSetCallback(val_h, "BUTTON_RELEASE_CB", cfunction(button_release, Int, (Ptr{Ihandle}, Float64)))
	IupSetCallback(val_h, "MOUSEMOVE_CB",      cfunction(mousemove, Int, (Ptr{Ihandle}, Float64)))

	IupSetAttribute(dlg_val, "TITLE", "IupVal")
	IupSetAttribute(dlg_val, "MARGIN", "10x10")
	IupShowXY(dlg_val, IUP_CENTER, IUP_CENTER)  # Displays the dialog (at screen center)

	IupMainLoop()                               # Initializes IUP main loop
	IupClose()                                  # And close it when ready
end

function button_press(c::Ptr{Ihandle}, a::Float64)
	global lbl_v, lbl_h
	type_ = bytestring(IupGetAttribute(c, "ORIENTATION"))

	if (type_[1] == 'V')
		IupSetAttribute(lbl_v, "FGCOLOR", "255 0 0")
	elseif (type_[1] == 'H')
		IupSetAttribute(lbl_h, "FGCOLOR", "255 0 0")
	end

	mousemove(c, a)

	return IUP_DEFAULT
end

function button_release(c::Ptr{Ihandle}, a::Float64)
	global lbl_v, lbl_h
	type_ = bytestring(IupGetAttribute(c, "ORIENTATION"))

	if (type_[1] == 'V')
		IupSetAttribute(lbl_v, "FGCOLOR", "0 0 0")
	elseif (type_[1] == 'H')
		IupSetAttribute(lbl_h, "FGCOLOR", "0 0 0")
	end

	mousemove(c, a)

	return IUP_DEFAULT
end

function mousemove(c::Ptr{Ihandle}, a::Float64)
	global lbl_v, lbl_h

	buffer = @sprintf("VALUE=%.2f", a)
	type_ = bytestring(IupGetAttribute(c, "ORIENTATION"))

	if (type_[1] == 'V')
		IupStoreAttribute(lbl_v, "TITLE", buffer)
	elseif (type_[1] == 'H')
		IupStoreAttribute(lbl_h, "TITLE", buffer)
	end

	return IUP_DEFAULT
end

end 	# module


text:
using Documenter
makedocs(;
    sitename="Ji"
)


text:
using RecipesBase, VisualRegressionTests, Plots, Colors
include("utils.jl")
# Note that the `false`s avoid popups and the tests directly fail.
# If you have changed something and need to recreate the reference in the test avoids to start Gtk

@testset "Recipes Test" begin
    references_folder = joinpath(@__DIR__, "assets")
    @testset "2D Recipes in GR" begin
        ENV["GKSwstype"] = "100"
        gr()
        function Hyp2PB_plot()
            M = Hyperbolic(2)
            p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 1.0]])
            p2 = convert.(Ref(PoincareBallPoint), p)
            return plot(M, p2)
        end
        @plottest Hyp2PB_plot joinpath(references_folder, "Hyp2PBPlot.png") false

        function Hyp2PB_plot_geo()
            M = Hyperbolic(2)
            p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 1.0]])
            p2 = convert.(Ref(PoincareBallPoint), p)
            return plot(M, p2; geodesic_interpolation=80)
        end
        @plottest Hyp2PB_plot_geo joinpath(references_folder, "Hyp2PBPlotGeo.png") false

        function Hyp2PB_quiver()
            M = Hyperbolic(2)
            p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 0.0], [0.0, 1.0]])
            p2 = convert.(Ref(PoincareBallPoint), p)
            X = [log(M, p2[2], p2[1]), log(M, p2[1], p2[3])]
            return plot(M, [p2[2], p2[1]], X)
        end
        @plottest Hyp2PB_quiver joinpath(references_folder, "Hyp2PBQuiver.png") false

        function Hyp2PH_plot()
            M = Hyperbolic(2)
            p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 1.0]])
            p2 = convert.(Ref(PoincareHalfSpacePoint), p)
            return plot(M, p2)
        end
        @plottest Hyp2PH_plot joinpath(references_folder, "Hyp2PHPlot.png") false

        function Hyp2PH_plot_geo()
            M = Hyperbolic(2)
            p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 1.0]])
            p2 = convert.(Ref(PoincareHalfSpacePoint), p)
            return plot(M, p2; geodesic_interpolation=80)
        end
        @plottest Hyp2PH_plot_geo joinpath(references_folder, "Hyp2PHPlotGeo.png") false

        function Hyp2PH_quiver()
            M = Hyperbolic(2)
            p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 0.0], [0.0, 1.0]])
            p2 = convert.(Ref(PoincareHalfSpacePoint), p)
            X = [log(M, p2[2], p2[1]), log(M, p2[1], p2[3])]
            return plot(M, [p2[2], p2[1]], X)
        end
        @plottest Hyp2PH_quiver joinpath(references_folder, "Hyp2PHQuiver.png") false
    end
    @testset "3D Recipes in GR" begin
        ENV["GKSwstype"] = "100"
        gr()
        function Hyp2_plot()
            M = Hyperbolic(2)
            p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 1.0]])
            return plot(M, p)
        end
        @plottest Hyp2_plot joinpath(references_folder, "Hyp2Plot.png") false

        function Hyp2_surfplot()
            M = Hyperbolic(2)
            p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 1.0]])
            return plot(M, p; surface=true)
        end
        @plottest Hyp2_surfplot joinpath(references_folder, "Hyp2SurfPlot.png") false

        function Hyp2_plot_geo()
            M = Hyperbolic(2)
            p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 1.0]])
            return plot(M, p; geodesic_interpolation=80)
        end
        @plottest Hyp2_plot_geo joinpath(references_folder, "Hyp2PlotGeo.png") false

        function Hyp2_quiver()
            M = Hyperbolic(2)
            p = Manifolds._hyperbolize.(Ref(M), [[1.0, 0.0], [0.0, 0.0], [0.0, 1.0]])
            X = [log(M, p[2], p[1]), log(M, p[1], p[3])]
            return plot(M, [p[2], p[1]], X)
        end
        @plottest Hyp2_quiver joinpath(references_folder, "Hyp2Quiver.png") false
    end
    @testset "3D Recipes in pyplot" begin
        pyplot()
        function Sphere2_plot()
            M = Sphere(2)
            pts = [[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]]
            return plot(M, pts; wireframe_color=colorant"#CCCCCC", markersize=10)
        end
        @plottest Sphere2_plot joinpath(references_folder, "Sphere2Plot.png") false

        function Sphere2_surfplot()
            M = Sphere(2)
            pts = [[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]]
            return plot(
                M,
                pts;
                surface=true,
                wireframe_color=colorant"#CCCCCC",
                markersize=10,
            )
        end
        @plottest Sphere2_surfplot joinpath(references_folder, "Sphere2SurfPlot.png") false

        function Sphere2_plot_geo()
            M = Sphere(2)
            pts = [[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.0]]
            return Plots.plot(
                M,
                pts;
                wireframe_color=colorant"#CCCCCC",
                geodesic_interpolation=80,
            )
        end
        @plottest Sphere2_plot_geo joinpath(references_folder, "Sphere2PlotGeo.png") false

        function Sphere2_quiver()
            pyplot()
            M = Sphere(2)
            pts2 = [[1.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0]]
            p3 = 1 / sqrt(3) .* [1.0, -1.0, 1.0]
            vecs = log.(Ref(M), pts2, Ref(p3))
            return plot(M, pts2, vecs; wireframe_color=colorant"#CCCCCC", linewidth=1.5)
        end
        @plottest Sphere2_quiver joinpath(references_folder, "Sphere2Quiver.png") false
    end
end


text:
struct SortExpr
  expr::Any
end

# Base.isless(e1::Expr, e2::Expr) = _isless_sortexpr(e1, e2)
Base.isless(e1::SortExpr, e2::SortExpr) = _isless_sortexpr(e1.expr, e2.expr)
function _isless_sortexpr(e1::Expr, e2::Expr)
  @iftrue e1.head < e2.head
  @iftrue length(e1.args) < length(e2.args)
  all(zip(e1.args, e2.args)) do (a1, a2)
    _isless_sortexpr(a1, a2)
  end
end
_isless_sortexpr(a1::T, a2::T) where T = a1 < a2
function _isless_sortexpr(l1::LineNumberNode, l2::LineNumberNode)
  l1.line < l2.line && l1.file < l2.file
end
_isless_sortexpr(a1, a2) = true

function sortexpr(a::Vector)
  [s.expr for s in sort(SortExpr.(a))]
end


text:
"""
    function expintegrator(A, t::Number, u₀, u₁, …; kwargs...)
    function expintegrator(A, t::Number, (u₀, u₁, …); kwargs...)
    function expintegrator(A, t::Number, (u₀, u₁, …), algorithm)

Compute ``y = ϕ₀(t*A)*u₀ + t*ϕ₁(t*A)*u₁ + t^2*ϕ₂(t*A)*u₂ + …``, where `A` is a general linear map, i.e. a `AbstractMatrix` or just a general function or callable object and `u₀`, `u₁` are of any Julia type with vector like behavior. Here, ``ϕ₀(z) = exp(z)`` and ``ϕⱼ₊₁ = (ϕⱼ(z) - 1/j!)/z``. In particular, ``y = x(t)`` represents the solution of the ODE
``x' = A*x + ∑ⱼ t^j/j! uⱼ₊₁`` with ``x(0) = u₀``.

!!! note
    When there are only input vectors `u₀` and `u₁`, `t` can equal `Inf`, in which the
    algorithm tries to evolve all the way to the fixed point `y = - A \\ u₁ + P₀ u₀` with
    `P₀` the projector onto the eigenspace of eigenvalue zero (if any) of `A`. If `A` has
    any eigenvalues with real part larger than zero, however, the solution to the ODE will
    diverge, i.e. the fixed point is not stable.

!!! warning
    The returned solution might be the solution of the ODE integrated up to a smaller time ``t̃ = sign(t) * |t̃|`` with ``|t̃| < |t|``, when the required precision could not be attained. Always check `info.converged > 0` or `info.residual == 0` (see below).

### Arguments:
The linear map `A` can be an `AbstractMatrix` (dense or sparse) or a general function or
callable object that implements the action of the linear map on a vector. If `A` is an
`AbstractMatrix`, `x` is expected to be an `AbstractVector`, otherwise `x` can be of any
type that behaves as a vector and supports the required methods (see KrylovKit docs).

The time parameter `t` can be real or complex, and it is better to choose `t` e.g. imaginary
and `A` hermitian, then to absorb the imaginary unit in an antihermitian `A`. For the
former, the Lanczos scheme is used to built a Krylov subspace, in which an approximation to
the exponential action of the linear map is obtained. The arguments `u₀`, `u₁`, … can be
of any type and should be in the domain of `A`.

### Return values:
The return value is always of the form `y, info = expintegrator(...)` with
*   `y`: the result of the computation, i.e.
    ``y = ϕ₀(t̃*A)*u₀ + t̃*ϕ₁(t̃*A)*u₁ + t̃^2*ϕ₂(t̃*A)*u₂ + …``
    with ``t̃ = sign(t) * |t̃|`` with ``|t̃| <= |t|``, such that the accumulated error in
    `y` per unit time is at most equal to the keyword argument `tol`
*   `info`: an object of type [`ConvergenceInfo`], which has the following fields
    -   `info.converged::Int`: 0 or 1 if the solution `y` was evolved all the way up to the
        requested time `t`.
    -   `info.residual`: there is no residual in the conventional sense, however, this
        value equals the residual time `t - t̃`, i.e. it is zero if `info.converged == 1`
    -   `info.normres::Real`: a (rough) estimate of the total error accumulated in the
        solution, should be smaller than `tol * |t̃|`
    -   `info.numops::Int`: number of times the linear map was applied, i.e. number of times
        `f` was called, or a vector was multiplied with `A`
    -   `info.numiter::Int`: number of times the Krylov subspace was restarted (see below)

### Keyword arguments:
Keyword arguments and their default values are given by:
*   `verbosity::Int = 0`: verbosity level, i.e. 0 (no messages), 1 (single message
    at the end), 2 (information after every iteration), 3 (information per Krylov step)
*   `krylovdim = 30`: the maximum dimension of the Krylov subspace that will be constructed.
    Note that the dimension of the vector space is not known or checked, e.g. `x₀` should
    not necessarily support the `Base.length` function. If you know the actual problem
    dimension is smaller than the default value, it is useful to reduce the value of
    `krylovdim`, though in principle this should be detected.
*   `tol = 1e-12`: the requested accuracy per unit time, i.e. if you want a certain
    precision `ϵ` on the final result, set `tol = ϵ/abs(t)`. If you work in e.g. single
    precision (`Float32`), you should definitely change the default value.
*   `maxiter::Int = 100`: the number of times the Krylov subspace can be rebuilt; see below
    for further details on the algorithms.
*   `issymmetric`: if the linear map is symmetric, only meaningful if `T<:Real`
*   `ishermitian`: if the linear map is hermitian
The default value for the last two depends on the method. If an `AbstractMatrix` is used,
`issymmetric` and `ishermitian` are checked for that matrix, ortherwise the default values
are `issymmetric = false` and `ishermitian = T <: Real && issymmetric`.

### Algorithm
The last method, without keyword arguments and the different vectors `u₀`, `u₁`, … in a
tuple, is the one that is finally called, and can also be used directly. Here, one
specifies the algorithm explicitly as either [`Lanczos`](@ref), for real symmetric or
complex hermitian linear maps, or [`Arnoldi`](@ref), for general linear maps. Note that
these names refer to the process for building the Krylov subspace, and that one can still
use complex time steps in combination with e.g. a real symmetric map.
"""
function expintegrator end

function expintegrator(A, t::Number, u₀, us...; kwargs...)
    alg = eigselector(A, promote_type(typeof(t), eltype(u₀), eltype.(us)...); kwargs...)
    expintegrator(A, t, (u₀, us...), alg)
end

function expintegrator(A, t::Number, u::Tuple, alg::Union{Lanczos,Arnoldi})
    length(u) == 1 && return expintegrator(A, t, (u[1], rmul!(similar(u[1]), false)), alg)

    p = length(u) - 1

    # process initial vector and determine result type
    u₀ = first(u)
    β₀ = norm(u₀)
    Au₀ = apply(A, u₀) # used to determine return type
    numops = 1
    T = promote_type(promote_type(eltype(Au₀), typeof(β₀), typeof(t)),
                        promote_type(eltype.(u)...))
    S = real(T)
    w₀ = copyto!(similar(u₀, T), u₀)

    # krylovdim and related allocations
    krylovdim = alg.krylovdim
    K = krylovdim
    HH = zeros(T, (krylovdim+p+1, krylovdim+p+1))

    # time step parameters
    η::S = alg.tol # tol is per unit time
    totalerr = zero(η)
    sgn = sign(t)
    τ::S = abs(t)
    Δτ::S = one(τ) # don't try any clever initial guesses, rely on correction mechanism
    τ₀ = zero(τ)

    # safety factors
    δ::S = 1.2
    γ::S = 0.8

    # initial vectors
    w = Vector{typeof(w₀)}(undef, p+1)
    w[1] = w₀
    for j = 1:p
        w[j+1] = apply(A, w[j])
        numops += 1
        lfac = 1
        for l = 0:p-j
            w[j+1] = axpy!((sgn*τ₀)^l/lfac, u[j+l+1], w[j+1])
            lfac *= l+1
        end
    end
    v = similar(w₀)
    β = norm(w[p+1])
    if β < alg.tol && p == 1
        if alg.verbosity > 0
            @info """expintegrate finished after 0 iterations, converged to fixed point up to error = $β"""
        end
        return w₀, ConvergenceInfo(1, zero(τ), β, 0, numops)
    end
    mul!(v, w[p+1], 1/β)

    # initialize iterator
    if alg isa Lanczos
        iter = LanczosIterator(A, w[p+1], alg.orth)
    else
        iter = ArnoldiIterator(A, w[p+1], alg.orth)
    end
    fact = initialize(iter; verbosity = alg.verbosity - 2)
    numops += 1
    sizehint!(fact, krylovdim)

    # start outer iteration loop
    maxiter = alg.maxiter
    numiter = 0
    while true
        if β < alg.tol && p == 1 # w₀ is fixed point of ODE
            if alg.verbosity > 0
                @info """expintegrate finished after $numiter iterations, converged to fixed point up to error = $β"""
            end
            return w₀, ConvergenceInfo(1, zero(τ), β, numiter, numops)
        end

        numiter += 1
        Δτ = min(Δτ, τ-τ₀)

        # Lanczos or Arnoldi factorization
        while normres(fact) > eps() && length(fact) < krylovdim
            fact = expand!(iter, fact; verbosity = alg.verbosity-2)
            numops += 1
        end
        K = fact.k # current Krylov dimension
        V = basis(fact)

        # Small matrix exponential and error estimation
        H = fill!(view(HH, 1:K+p+1, 1:K+p+1), zero(T))
        mul!(view(H, 1:K, 1:K), rayleighquotient(fact), sgn*Δτ)
        H[1, K+1] = 1
        for i = 1:p
            H[K+i, K+i+1] = 1
        end
        expH = LinearAlgebra.exp!(H)
        ϵ = abs(Δτ^p * β * normres(fact) * expH[K,K+p+1])
        ω = ϵ / (Δτ * η)

        q = K/2
        while ω > one(ω)
            ϵ_prev = ϵ
            Δτ_prev = Δτ
            Δτ *= (γ/ω)^(1/(q+1))
            H = fill!(view(HH, 1:K+p+1, 1:K+p+1), zero(T))
            mul!(view(H, 1:K, 1:K), rayleighquotient(fact), sgn*Δτ)
            H[1, K+1] = 1
            for i = 1:p
                H[K+i, K+i+1] = 1
            end
            expH = LinearAlgebra.exp!(H)
            ϵ = abs(Δτ^p * β * normres(fact) * expH[K,K+p+1])
            ω = ϵ / (Δτ * η)
            q = max(zero(q),  log(ϵ / ϵ_prev)/log(Δτ / Δτ_prev)-1)
        end

        # take time step
        totalerr += ϵ
        τ₀ += Δτ
        jfac = 1
        for j = 1:p-1
            w₀ = axpy!((sgn*Δτ)^j/jfac, w[j+1], w₀)
            jfac *= (j+1)
        end
        w[p+1] = mul!(w[p+1], basis(fact), view(expH, 1:K, K+p))
        # add first correction
        w[p+1] = axpy!(expH[K,K+p+1], residual(fact), w[p+1])
        w₀ = axpy!(β*(sgn*Δτ)^p, w[p+1], w₀)

        # increase time step for next iteration:
        if ω < γ
            Δτ *= (γ/ω)^(1/(q+1))
        end

        if alg.verbosity > 1
            msg = "expintegrate in iteration $numiter: "
            msg *= "reached time " * @sprintf("%.2e", τ₀)
            msg *= ", total error = " * @sprintf("%.4e", totalerr)
            @info msg
        end

        if τ₀ >= τ
            if alg.verbosity > 0
                @info """expintegrate finished after $numiter iterations: total error = $totalerr"""
            end
            return w₀, ConvergenceInfo(1, zero(τ), totalerr, numiter, numops)
        elseif numiter == maxiter
            if alg.verbosity > 0
                @warn """expintegrate finished without convergence after $numiter iterations:
                total error = $totalerr, residual time = $(τ - τ₀)"""
            end
            return w₀, ConvergenceInfo(0, τ-τ₀, totalerr, numiter, numops)
        else
            for j = 1:p
                w[j+1] = apply(A, w[j])
                numops += 1
                lfac = 1
                for l = 0:p-j
                    w[j+1] = axpy!((sgn*τ₀)^l/lfac, u[j+l+1], w[j+1])
                    lfac *= l+1
                end
            end
            β = norm(w[p+1])
            mul!(v, w[p+1], 1/β)
            if alg isa Lanczos
                iter = LanczosIterator(A, w[p+1], alg.orth)
            else
                iter = ArnoldiIterator(A, w[p+1], alg.orth)
            end
            fact = initialize!(iter, fact; verbosity = alg.verbosity-2)
            numops += 1
        end
    end
end


text:
using .Graphs

export ContMarkovChain, add_state!, add_transition!, state_count

struct ContMarkovChain
    state_graph::Digraph
    trans_rates::Vector{Float64}
    function ContMarkovChain()
        new(Digraph(), Vector{Float64}())
    end
    function ContMarkovChain(nstates::Integer)
        new(Digraph(nstates), Vector{Float64}())
    end
end

struct Transition
    src::Int
    dst::Int
    rate::Float64
end

add_state!(chain::ContMarkovChain) = add_node!(chain.state_graph)

function add_transition!(chain::ContMarkovChain, src::Integer, dst::Integer, rate::Real)
    idx = add_edge!(chain.state_graph, src, dst)
    push!(chain.trans_rates, rate)
    return idx
end

state_count(chain) = node_count(chain.state_graph)
states(chain) = nodes(chain.state_graph)
function transitions(chain)
    ts = Vector{Transition}()
    for eid in edges(chain.state_graph)
        ep = chain.state_graph.edge_props[eid]
        push!(ts, Transition(ep.src, ep.dst, chain.trans_rates[eid]))
    end
    return ts
end


text:
using MPI
using CLIMA
using Logging
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using CLIMA.MPIStateArrays
using LinearAlgebra
using Printf
using Dates
using CLIMA.GenericCallbacks: EveryXWallTimeSeconds, EveryXSimulationSteps
using CLIMA.ODESolvers
using CLIMA.VTK: writevtk, writepvtu
using CLIMA.Mesh.Grids: min_node_distance

const output = parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_OUTPUT","false")))

if !@isdefined integration_testing
  const integration_testing =
    parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
end

include("hyperdiffusion_model.jl")

struct ConstantHyperDiffusion{dim, dir, FT} <: HyperDiffusionProblem
  D::SMatrix{3, 3, FT, 9}
end

function init_hyperdiffusion_tensor!(problem::ConstantHyperDiffusion, aux::Vars,
                                     geom::LocalGeometry)
  aux.D = problem.D
end

function initial_condition!(problem::ConstantHyperDiffusion{dim, dir}, state, aux, x, t) where {dim, dir}
  @inbounds begin
    k = SVector(1, 2, 3)
    kD = k * k' .* problem.D
    if dir === EveryDirection()
      c = sum(abs2, k[SOneTo(dim)]) * sum(kD[SOneTo(dim), SOneTo(dim)])
    elseif dir === HorizontalDirection()
      c = sum(abs2, k[SOneTo(dim - 1)]) * sum(kD[SOneTo(dim - 1), SOneTo(dim - 1)])
    elseif dir === VerticalDirection()
      c = k[dim] ^ 2 * kD[dim, dim]
    end
    state.ρ = sin(dot(k[SOneTo(dim)], x[SOneTo(dim)])) * exp(-c * t)
  end
end

function do_output(mpicomm, vtkdir, vtkstep, dg, Q, Qe, model, testname)
  ## name of the file that this MPI rank will write
  filename = @sprintf("%s/%s_mpirank%04d_step%04d",
                      vtkdir, testname, MPI.Comm_rank(mpicomm), vtkstep)

  statenames = flattenednames(vars_state(model, eltype(Q)))
  exactnames = statenames .* "_exact"

  writevtk(filename, Q, dg, statenames, Qe, exactnames)

  ## Generate the pvtu file for these vtk files
  if MPI.Comm_rank(mpicomm) == 0
    ## name of the pvtu file
    pvtuprefix = @sprintf("%s/%s_step%04d", vtkdir, testname, vtkstep)

    ## name of each of the ranks vtk files
    prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
      @sprintf("%s_mpirank%04d_step%04d", testname, i - 1, vtkstep)
    end

    writepvtu(pvtuprefix, prefixes, (statenames..., exactnames...))

    @info "Done writing VTK: $pvtuprefix"
  end
end

function run(mpicomm, ArrayType, dim, topl, N, timeend, FT, direction,
             D, vtkdir, outputtime)

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                         )

  dx = min_node_distance(grid)
  dt = dx ^ 4 / 25 / sum(D)
  @info "time step" dt
  dt = outputtime / ceil(Int64, outputtime / dt)

  model = HyperDiffusion{dim}(ConstantHyperDiffusion{dim, direction(), FT}(D))
  dg = DGModel(model,
               grid,
               CentralNumericalFluxNonDiffusive(),
               CentralNumericalFluxDiffusive(),
               CentralNumericalFluxGradient(),
               direction=direction())

  Q = init_ode_state(dg, FT(0))

  lsrk = LSRK54CarpenterKennedy(dg, Q; dt = dt, t0 = 0)

  eng0 = norm(Q)
  @info @sprintf """Starting
  norm(Q₀) = %.16e""" eng0

  # Set up the information callback
  starttime = Ref(now())
  cbinfo = EveryXWallTimeSeconds(60, mpicomm) do (s=false)
    if s
      starttime[] = now()
    else
      energy = norm(Q)
      @info @sprintf("""Update
                     simtime = %.16e
                     runtime = %s
                     norm(Q) = %.16e""", gettime(lsrk),
                     Dates.format(convert(Dates.DateTime,
                                          Dates.now()-starttime[]),
                                  Dates.dateformat"HH:MM:SS"),
                     energy)
    end
  end
  callbacks = (cbinfo,)
  if ~isnothing(vtkdir)
    # create vtk dir
    mkpath(vtkdir)

    vtkstep = 0
    # output initial step
    do_output(mpicomm, vtkdir, vtkstep, dg, Q, Q, model, "hyperdiffusion")

    # setup the output callback
    cbvtk = EveryXSimulationSteps(floor(outputtime/dt)) do
      vtkstep += 1
      Qe = init_ode_state(dg, gettime(lsrk))
      do_output(mpicomm, vtkdir, vtkstep, dg, Q, Qe, model,
                "hyperdiffusion")
    end
    callbacks = (callbacks..., cbvtk)
  end

  solve!(Q, lsrk; timeend=timeend, callbacks=callbacks)

  # Print some end of the simulation information
  engf = norm(Q)
  Qe = init_ode_state(dg, FT(timeend))

  engfe = norm(Qe)
  errf = euclidean_distance(Q, Qe)
  @info @sprintf """Finished
  norm(Q)                 = %.16e
  norm(Q) / norm(Q₀)      = %.16e
  norm(Q) - norm(Q₀)      = %.16e
  norm(Q - Qe)            = %.16e
  norm(Q - Qe) / norm(Qe) = %.16e
  """ engf engf/eng0 engf-eng0 errf errf / engfe
  errf
end

using Test
let
  CLIMA.init()
  ArrayType = CLIMA.array_type()
  mpicomm = MPI.COMM_WORLD

  numlevels = integration_testing || CLIMA.Settings.integration_testing ? 3 : 1

  polynomialorder = 4
  base_num_elem = 4

  expected_result = Dict()
  expected_result[2, 1, Float64, EveryDirection] = 7.6772960298563120e-03
  expected_result[2, 2, Float64, EveryDirection] = 2.3268371815073617e-03
  expected_result[2, 3, Float64, EveryDirection] = 4.2641957936779901e-05
  expected_result[2, 1, Float64, HorizontalDirection] = 8.4650675812606650e-04
  expected_result[2, 2, Float64, HorizontalDirection] = 4.4626814795055979e-05
  expected_result[2, 3, Float64, HorizontalDirection] = 1.2193396277823764e-06
  expected_result[2, 1, Float64, VerticalDirection] = 6.1690465335730834e-03
  expected_result[2, 2, Float64, VerticalDirection] = 2.3407593209031621e-03
  expected_result[2, 3, Float64, VerticalDirection] = 4.3775160749787010e-05

  expected_result[3, 1, Float64, EveryDirection] = 1.7363355506160003e-01
  expected_result[3, 2, Float64, EveryDirection] = 7.3049474767548042e-02
  expected_result[3, 3, Float64, EveryDirection] = 5.8530711333407105e-04
  expected_result[3, 1, Float64, HorizontalDirection] = 1.9244127301149615e-02
  expected_result[3, 2, Float64, HorizontalDirection] = 5.8325158696244947e-03
  expected_result[3, 3, Float64, HorizontalDirection] = 1.0688753745025491e-04
  expected_result[3, 1, Float64, VerticalDirection] = 1.4412891107361228e-01
  expected_result[3, 2, Float64, VerticalDirection] = 6.3744013545812925e-02
  expected_result[3, 3, Float64, VerticalDirection] = 9.0891011404938341e-04

  numlevels = integration_testing ? 3 : 1
  for FT in (Float64,)
    D = 1 // 100 * SMatrix{3, 3, FT}(9 // 50, 3 // 50, 5  // 50,
                                     3 // 50, 7 // 50, 4  // 50,
                                     5 // 50, 4 // 50, 10 // 50)

    result = zeros(FT, numlevels)
    for dim in (2, 3)
      for direction in (EveryDirection, HorizontalDirection, VerticalDirection)
        for l = 1:numlevels
          Ne = 2^(l-1) * base_num_elem
          xrange = range(FT(0); length=Ne+1, stop=FT(2pi))
          brickrange = ntuple(j->xrange, dim)
          periodicity = ntuple(j->true, dim)
          topl = StackedBrickTopology(mpicomm, brickrange;
                                      periodicity = periodicity)
          timeend = 1
          outputtime = 1

          @info (ArrayType, FT, dim, direction)
          vtkdir = output ? "vtk_hyperdiffusion" *
                            "_poly$(polynomialorder)" *
                            "_dim$(dim)_$(ArrayType)_$(FT)_$(direction)" *
                            "_level$(l)" : nothing
          result[l] = run(mpicomm, ArrayType, dim, topl, polynomialorder,
                          timeend, FT, direction, D, vtkdir,
                          outputtime)

          @test result[l] ≈ FT(expected_result[dim, l, FT, direction])
        end
        @info begin
          msg = ""
          for l = 1:numlevels-1
            rate = log2(result[l]) - log2(result[l+1])
            msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
          end
          msg
        end
      end
    end
  end
end

nothing

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
export BrickletAirQualityAllValues
struct BrickletAirQualityAllValues
    iaq_index::Integer
    iaq_index_accuracy::Integer
    temperature::Integer
    humidity::Integer
    air_pressure::Integer
end

export BrickletAirQualityAllValuesCallbackConfiguration
struct BrickletAirQualityAllValuesCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletAirQualityIAQIndex
struct BrickletAirQualityIAQIndex
    iaq_index::Integer
    iaq_index_accuracy::Integer
end

export BrickletAirQualityIAQIndexCallbackConfiguration
struct BrickletAirQualityIAQIndexCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletAirQualityTemperatureCallbackConfiguration
struct BrickletAirQualityTemperatureCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletAirQualityHumidityCallbackConfiguration
struct BrickletAirQualityHumidityCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletAirQualityAirPressureCallbackConfiguration
struct BrickletAirQualityAirPressureCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletAirQualitySPITFPErrorCount
struct BrickletAirQualitySPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletAirQualityIdentity
struct BrickletAirQualityIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletAirQuality
"""
Measures IAQ index, temperature, humidity and air pressure
"""
mutable struct BrickletAirQuality <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletAirQuality(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_air_quality")
        deviceInternal = package.BrickletAirQuality(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_all_values
"""
    $(SIGNATURES)

Returns all values measured by the Air Quality Bricklet. The values are
IAQ (Indoor Air Quality) Index (higher value means greater level of air pollution), IAQ Index Accuracy, Temperature, Humidity and
Air Pressure.

.. image:: /Images/Misc/bricklet_air_quality_iaq_index.png
   :scale: 100 %
   :alt: Air Quality Index description
   :align: center
   :target: ../../_images/Misc/bricklet_air_quality_iaq_index.png
"""
function get_all_values(device::BrickletAirQuality)
    return device.deviceInternal.get_all_values()
end

export set_temperature_offset
"""
    $(SIGNATURES)

Sets a temperature offset. A offset of 10 will decrease the measured temperature by 0.1 °C.

If you install this Bricklet into an enclosure and you want to measure the ambient
temperature, you may have to decrease the measured temperature by some value to
compensate for the error because of the heating inside of the enclosure.

We recommend that you leave the parts in the enclosure running for at least
24 hours such that a temperature equilibrium can be reached. After that you can measure
the temperature directly outside of enclosure and set the difference as offset.

This temperature offset is used to calculate the relative humidity and
IAQ index measurements. In case the Bricklet is installed in an enclosure, we
recommend to measure and set the temperature offset to improve the accuracy of
the measurements.
"""
function set_temperature_offset(device::BrickletAirQuality, offset)
    device.deviceInternal.set_temperature_offset(offset)
end

export get_temperature_offset
"""
    $(SIGNATURES)

Returns the temperature offset as set by
:func:`Set Temperature Offset`.
"""
function get_temperature_offset(device::BrickletAirQuality)
    return device.deviceInternal.get_temperature_offset()
end

export set_all_values_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`All Values`
callback is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after at least one of the values has changed. If the values didn't
change within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_all_values_callback_configuration(device::BrickletAirQuality, period, value_has_to_change)
    device.deviceInternal.set_all_values_callback_configuration(period, value_has_to_change)
end

export get_all_values_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by
:func:`Set All Values Callback Configuration`.
"""
function get_all_values_callback_configuration(device::BrickletAirQuality)
    return device.deviceInternal.get_all_values_callback_configuration()
end

export get_iaq_index
"""
    $(SIGNATURES)

Returns the IAQ index and accuracy. The higher the IAQ index, the greater the level of air pollution.

.. image:: /Images/Misc/bricklet_air_quality_iaq_index.png
   :scale: 100 %
   :alt: IAQ index description
   :align: center
   :target: ../../_images/Misc/bricklet_air_quality_iaq_index.png

If you want to get the value periodically, it is recommended to use the
:cb:`IAQ Index` callback. You can set the callback configuration
with :func:`Set IAQ Index Callback Configuration`.
"""
function get_iaq_index(device::BrickletAirQuality)
    return device.deviceInternal.get_iaq_index()
end

export set_iaq_index_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`IAQ Index`
callback is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after at least one of the values has changed. If the values didn't
change within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_iaq_index_callback_configuration(device::BrickletAirQuality, period, value_has_to_change)
    device.deviceInternal.set_iaq_index_callback_configuration(period, value_has_to_change)
end

export get_iaq_index_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by
:func:`Set IAQ Index Callback Configuration`.
"""
function get_iaq_index_callback_configuration(device::BrickletAirQuality)
    return device.deviceInternal.get_iaq_index_callback_configuration()
end

export get_temperature
"""
    $(SIGNATURES)

Returns temperature.

If you want to get the value periodically, it is recommended to use the
:cb:`Temperature` callback. You can set the callback configuration
with :func:`Set Temperature Callback Configuration`.
"""
function get_temperature(device::BrickletAirQuality)
    return device.deviceInternal.get_temperature()
end

export set_temperature_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Temperature` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Temperature` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_temperature_callback_configuration(device::BrickletAirQuality, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_temperature_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_temperature_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Temperature Callback Configuration`.
"""
function get_temperature_callback_configuration(device::BrickletAirQuality)
    return device.deviceInternal.get_temperature_callback_configuration()
end

export get_humidity
"""
    $(SIGNATURES)

Returns relative humidity.

If you want to get the value periodically, it is recommended to use the
:cb:`Humidity` callback. You can set the callback configuration
with :func:`Set Humidity Callback Configuration`.
"""
function get_humidity(device::BrickletAirQuality)
    return device.deviceInternal.get_humidity()
end

export set_humidity_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Humidity` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Humidity` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_humidity_callback_configuration(device::BrickletAirQuality, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_humidity_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_humidity_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Humidity Callback Configuration`.
"""
function get_humidity_callback_configuration(device::BrickletAirQuality)
    return device.deviceInternal.get_humidity_callback_configuration()
end

export get_air_pressure
"""
    $(SIGNATURES)

Returns air pressure.

If you want to get the value periodically, it is recommended to use the
:cb:`Air Pressure` callback. You can set the callback configuration
with :func:`Set Air Pressure Callback Configuration`.
"""
function get_air_pressure(device::BrickletAirQuality)
    return device.deviceInternal.get_air_pressure()
end

export set_air_pressure_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Air Pressure` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Air Pressure` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_air_pressure_callback_configuration(device::BrickletAirQuality, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_air_pressure_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_air_pressure_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Air Pressure Callback Configuration`.
"""
function get_air_pressure_callback_configuration(device::BrickletAirQuality)
    return device.deviceInternal.get_air_pressure_callback_configuration()
end

export remove_calibration
"""
    $(SIGNATURES)

Deletes the calibration from flash. After you call this function,
you need to power cycle the Air Quality Bricklet.

On the next power up the Bricklet will start a new calibration, as
if it was started for the very first time.

The calibration is based on the data of the last four days, so it takes
four days until a full calibration is re-established.

.. versionadded:: 2.0.3\$nbsp;(Plugin)
"""
function remove_calibration(device::BrickletAirQuality)
    device.deviceInternal.remove_calibration()
end

export set_background_calibration_duration
"""
    $(SIGNATURES)

The Air Quality Bricklet uses an automatic background calibration mechanism to
calculate the IAQ Index. This calibration mechanism considers a history of
measured data. The duration of this history can be configured to either be
4 days or 28 days.

If you keep the Bricklet mostly at one place and it does not get moved around
to different environments, we recommend that you use a duration of 28 days.

If you change the duration, the current calibration will be discarded and
the calibration will start from beginning again. The configuration of the
duration is saved in flash, so you should only have to call this function
once in the lifetime of the Bricklet.

The Bricklet has to be power cycled after this function is called
for a duration change to take effect.

Before firmware version 2.0.3 this was not configurable and the duration was
4 days.

The default value (since firmware version 2.0.3) is 28 days.

.. versionadded:: 2.0.3\$nbsp;(Plugin)
"""
function set_background_calibration_duration(device::BrickletAirQuality, duration)
    device.deviceInternal.set_background_calibration_duration(duration)
end

export get_background_calibration_duration
"""
    $(SIGNATURES)

Returns the background calibration duration as set by
:func:`Set Background Calibration Duration`.

.. versionadded:: 2.0.3\$nbsp;(Plugin)
"""
function get_background_calibration_duration(device::BrickletAirQuality)
    return device.deviceInternal.get_background_calibration_duration()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletAirQuality)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletAirQuality, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletAirQuality)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletAirQuality, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletAirQuality, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletAirQuality, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletAirQuality)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletAirQuality)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletAirQuality)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletAirQuality, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletAirQuality)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletAirQuality)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletAirQuality, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
import DataStructures
import XMLDict

function getXML(path)
    s = open(path) do file
        read(file,String)
    end
    return xml = XMLDict.xml_dict(s)

end

# get potentially missing value out of dict
function getVal(dict,key)
    try
        return val = parse(Float64,dict[key])
    catch y
        if isa(y,KeyError)
            return val = missing
        else
            throw(y)
        end
    end
end

struct XTbMLTable
    select::DataStructures.DefaultOrderedDict
    ultimate::DataStructures.DefaultOrderedDict
end

function loadXTbMLTable(path)

    x = getXML(path)
    name = x["XTbML"]["ContentClassification"]["TableName"]
    tbl = XTbMLTable(DataStructures.DefaultOrderedDict(missing),DataStructures.DefaultOrderedDict(missing))

    if isa(x["XTbML"]["Table"],Vector)
        # for a select and ultimate table, will have multiple tables
        # parsed into a vector of tables
        for ai in x["XTbML"]["Table"][1]["Values"]["Axis"]
            issue_age = parse(Int,ai[:t])
            tbl.select[issue_age] = DataStructures.DefaultOrderedDict(missing)
            for aj in ai["Axis"]["Y"]
                duration = parse(Int,aj[:t])
                rate = getVal(aj,"")
                if !ismissing(rate)
                    tbl.select[issue_age][duration]= rate
                end
            end
        end

        for ai in x["XTbML"]["Table"][2]["Values"]["Axis"]["Y"]
            age = parse(Int,ai[:t])
            rate = getVal(ai,"")
            if !ismissing(rate)
                tbl.ultimate[age]= rate
            end
        end
    elseif isa(x["XTbML"]["Table"],DataStructures.OrderedDict)
        # a table without select period will just have one set of values, which
        # are loaded into the OrderedDict

        for ai in x["XTbML"]["Table"]["Values"]["Axis"]["Y"]
            age = parse(Int,ai[:t])
            rate = getVal(ai,"")
            if !ismissing(rate)
                tbl.ultimate[age]= rate
            end
        end
    else
        error("don't know how to handle table: "  *  name)
    end

    return (tbl, name)
end