text:
export neighborcover, neighborcover_mapreduce, isclique
export generate_set, mirrors, mirrorcover, neighbors2

function neighborcover_mapreduce(func, red, eg::EliminateGraph, vi::Int)
    pre = func(vi)
    for vj in vertices(eg)
        if isconnected(eg,vj,vi)
            # find a neighbor cover
            pre = red(pre, func(vj))
        end
    end
    return pre
end

function isclique(eg::EliminateGraph, vs)
    for vj in vs
        for vi in vs
            vi!=vj && !isconnected(eg, vi, vj) && return false
        end
    end
    return true
end

isclique(eg::EliminateGraph) = isclique(eg, vertices(eg))

"""
    generate_set(eg::EliminateGraph, set::AbstractVertexSet) -> Vector

Generate (allocate) the set of vertices.
"""
function generate_set end

# nearest neighbors
@inline function generate_set(eg::EliminateGraph, nn::NearestNeighbors{OPEN})
    return filter(j->isconnected(eg,j,nn.i), vertices(eg))
end

# 2nd nearest neighbors, nbs can be either neighbors or neighbor cover
@inline function generate_set(eg::EliminateGraph, nnn::Neighbors{OPEN,2}, nbs)
    return filter(j->j!=nnn.i && !isconnected(eg,j,nnn.i) && any(nb->isconnected(eg,j,nb), nbs), vertices(eg))
end

@inline function generate_set(eg::EliminateGraph, nnn::Neighbors{CLOSED,2}, nbs)
    return filter(j->j in nbs || any(nb->isconnected(eg,j,nb), nbs), vertices(eg))
end
generate_set(eg::EliminateGraph, nnn::Neighbors{SP,2} where SP) = generate_set(eg,nnn,generate_set(eg,Neighbors{OPEN,1}(nnn.i)))

# nearest neighbors + itself
@inline function generate_set(eg::EliminateGraph, nn::NearestNeighbors{CLOSED})
    return filter(vj->isconnected(eg,vj,nn.i) || nn.i==vj, vertices(eg))
end

# mirrors, including itself
function generate_set(eg::EliminateGraph, set::Mirrors{CLOSED}, nc1, nbs2)
    push!(generate_set(eg, Mirrors{OPEN}(set.i), nc1, nbs2), set.i)
end

function generate_set(eg::EliminateGraph, set::Mirrors{CLOSED})
    nc1 = generate_set(eg, NeighborCover(set.i))
    nbs2 = generate_set(eg, Neighbors{OPEN, 2}(set.i), nc1)
    generate_set(eg, set, nc1, nbs2)
end

# mirrors
function generate_set(eg::EliminateGraph, set::Mirrors{OPEN}, nc1, nbs2)
    filter(v->isclique(eg, setdiff(nc1,neighbors(eg,v))), nbs2)
end

function generate_set(eg::EliminateGraph, set::Mirrors{OPEN})
    nc1 = generate_set(eg, NeighborCover(set.i))
    generate_set(eg, set, nc1, generate_set(eg, Neighbors{OPEN,2}(set.i), nc1))
end

function generate_set(eg::EliminateGraph, set::Vertex)
    [set.i]
end

function generate_set(eg::EliminateGraph, set::UnionOf)
    generate_set(eg, set.A) ∪ generate_set(eg, set.B)
end

########## Short cuts #############
"""
    neighbors2(eg::EliminateGraph, i::Int)

Get 2nd nearest neighbors, i.e. distance 2 vertices.
"""
neighbors2(eg::EliminateGraph, i::Int) = generate_set(eg,Neighbors{OPEN,2}(i))

"""
    mirrors(eg::EliminateGraph, vi::Int)

Get mirrors of vertex `vi`.
"""
mirrors(eg::EliminateGraph, vi::Int) = generate_set(eg, Mirrors{OPEN}(vi))

"""
    mirrorcover(eg::EliminateGraph, vi::Int)

Get mirrors of vertex `vi`, as well as itself.
"""
mirrorcover(eg::EliminateGraph, vi::Int) = generate_set(eg, Mirrors{CLOSED}(vi))


text:
using FunManifolds
using Test
using StaticArrays

include("../utils.jl")

@testset "CurveWarpingSRSFSpace" begin
    for N in [21, 50]
        nodes = range(0.0, 1.0, length=N)
        M = CurveWarpingSpace(nodes)
        MS = CurveWarpingSRSFSpace(nodes)

        test_manifold(
            MS,
            [srsf(M, t -> t^2), srsf(M, t -> t^3), srsf(M, t -> t^(0.7))];
            exp_log_atol_multiplier=N % 2 == 1 ? 1e8 : 10,
            test_default_vector_transport=true,
            test_reverse_diff=false,
            test_forward_diff=false,
            test_injectivity_radius=false,
        )
    end

    @testset "reverse SRSF" begin
        N = 11
        nodes = range(0.0, 1.0, length=N)
        M = CurveWarpingSpace(nodes)
        MS = CurveWarpingSRSFSpace(nodes)

        γ = FunManifolds.make_warping(M, [t^2 for t in nodes])
        γSRSF = srsf(M, γ)
        @test isapprox(M, γ, reverse_srsf(M, γSRSF), atol=1e-15)
    end
end

ext: jl


text:
# Tests only a small subset of the functions in Chromosome.jl

# So far, tests only arity 2 gates and with maxarity=2
using Test
@testset "chromosome_to_int() and enumerate_circuits() " begin
# set p.numlevelsback to be large since enumerate_circuits() and chromosome_to_int() assume maximum
p = Parameters(2,1,3,6)
funcs = default_funcs(p.numinputs)
nfuncs = length(funcs)
cl=enumerate_circuits(p,funcs);
cil = map(ch->chromosome_to_int(ch,funcs),cl);
@test length(cil) == length(unique(cil))
@test length(cil) == count_circuits(p,nfuncs=nfuncs)
ri = rand(cil)
@test ri == chromosome_to_int(int_to_chromosome(ri,p))  # test that chromosome_to_int() inverts int_to_chromosome()
end # @testst


text:
@testset "Test automatic file deletion" begin
    storage = IS.Hdf5TimeSeriesStorage()
    filename = storage.file_path
    @test isfile(filename)
    storage = nothing
    GC.gc()
    @test !isfile(filename)
end


text:
@testset "1232.check-if-it-is-a-straight-line.jl" begin
    @test check_straight_line([[1,2],[2,3],[3,4],[4,5],[5,6],[6,7]]) == true
    @test check_straight_line([[1,1],[2,2],[3,4],[4,5],[5,6],[7,7]]) == false
end


text:
using Pkg
try using ArraysOfArrays
catch
    Pkg.add("ArraysOfArrays")
    using ArraysOfArrays
end
try using CSV
catch
    Pkg.add("CSV")
    using CSV
end
try using Dates
catch
    Pkg.add("Dates")
    using Dates
end
try using Distributions
catch
    Pkg.add("Distributions")
    using Distributions
end
try using DSP
catch
    Pkg.add("DSP")
    using DSP
end
# using GERDAMetadata
try using Glob
catch
    Pkg.add("Glob")
    using Glob
end
try using HDF5
catch
    Pkg.add("HDF5")
    using HDF5
end
try using Interpolations
catch
    Pkg.add("Interpolations")
    using Interpolations
end
try using JSON
catch
    Pkg.add("JSON")
    using JSON
end
try using LaTeXStrings
catch
    Pkg.add("LaTeXStrings")
    using JSON
end
try using LegendDataTypes
catch
    Pkg.add(url="https://github.com/legend-exp/LegendDataTypes.jl")
    using LegendDataTypes
end
try using LegendHDF5IO
catch
    Pkg.add(url="https://github.com/legend-exp/LegendHDF5IO.jl")
    using LegendHDF5IO
end
using LegendHDF5IO: readdata, writedata
# try using UpROOT
# catch
#     Pkg.add("UpROOT")
#     using UpROOT
# end
# try using LegendUpROOTIO
# catch
#     Pkg.add(url="https://github.com/legend-exp/LegendUpROOTIO.jl")
#     using LegendUpROOTIO
# end
# using LegendUpROOTIO: read_mgtevent, raw2mgtevent, mgdo2legend
try using LsqFit
catch
    Pkg.add("LsqFit")
    using LsqFit
end
try using Measurements
catch
    Pkg.add("Measurements")
    using Measurements
end
try using Plots
catch
    Pkg.add("Plots")
    using Plots
end
try using ProgressMeter
catch
    Pkg.add("ProgressMeter")
    using ProgressMeter
end
try using Query
catch
    Pkg.add("Query")
    using Query
end
try using RadiationDetectorSignals
catch
    Pkg.add("RadiationDetectorSignals")
    using RadiationDetectorSignals
end
try using RadiationSpectra
catch
    Pkg.add("RadiationSpectra")
    using RadiationSpectra
end
try using SpecialFunctions
catch
    Pkg.add("SpecialFunctions")
    using SpecialFunctions
end
try using SolidStateDetectors
catch
    Pkg.add("SolidStateDetectors")
    using SolidStateDetectors
end
try using StaticArrays
catch
    Pkg.add("StaticArrays")
    using StaticArrays
end
try using StatsBase
catch
    Pkg.add("StatsBase")
    using StatsBase
end
try using StatsPlots
catch
    Pkg.add("StatsPlots")
    using StatsPlots
end
try using StructArrays
catch
    Pkg.add("StructArrays")
    using StructArrays
end
try using Suppressor
catch
    Pkg.add("Suppressor")
    using Suppressor
end
try using TypedTables
catch
    Pkg.add("TypedTables")
    using TypedTables
end
try using Unitful
catch
    Pkg.add("Unitful")
    using Unitful
end
try using Random
catch
    Pkg.add("Random")
    using Random
end
try using LinearAlgebra
catch
    Pkg.add("LinearAlgebra")
    using LinearAlgebra
end
try using Statistics
catch
    Pkg.add("Statistics")
    using Statistics
end
try using BAT
catch
    Pkg.add(url="https://github.com/bat/BAT.jl")
    using BAT
end
try using IntervalSets
catch
    Pkg.add("IntervalSets")
    using IntervalSets
end
try using ValueShapes
catch
    Pkg.add("ValueShapes")
    using ValueShapes
end
try using IJuliaBell
catch
    Pkg.add("IJuliaBell")
    using IJuliaBell
end

SSD = SolidStateDetectors
T = Float64
eval_ch = [0, 1, 2, 3, 4, 11, 12, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 30, 31, 32, 33, 34, 35]


text:
using Test

include("readSims.jl")
include("runTestExamples.jl")

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


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

#### Specialized matrix types ####

## (complex) symmetric tridiagonal matrices
struct SymTridiagonal{T, DV<:AbstractVector{T}, EV<:AbstractVector{T}} <: AbstractBandedMatrix{T}
    dv::DV                        # diagonal
    ev::EV                        # superdiagonal
    function SymTridiagonal{T, DV, EV}(dv, ev) where {T, DV<:AbstractVector{T}, EV<:AbstractVector{T}}
        require_one_based_indexing(dv, ev)
        if !(length(dv) - 1 <= length(ev) <= length(dv))
            throw(DimensionMismatch("subdiagonaldataonal has wrong length. Has length $(length(ev)), but should be either $(length(dv) - 1) or $(length(dv))."))
        end
        new{T, DV, EV}(dv, ev)
    end
end

"""
    SymTridiagonal(dv::DV, ev::EV) where {DV <: AbstractVector, EV <: AbstractVector}

Construct a symmetric tridiagonal matrix from the diagonal (`dv`) and first
sub/super-diagonal (`ev`), respectively. The result is of type `SymTridiagonal`
and provides efficient specialized eigensolvers, but may be converted into a
regular matrix with [`convert(Array, _)`](@ref) (or `Array(_)` for short).

For `SymTridiagonal` block matrices, the elements of `dv` are symmetrized.
The argument `ev` is interpreted as the superdiagonal. Blocks from the
subdiagonaldataonal are (materialized) transpose of the corresponding superdiagonal blocks.

# Examples
```jldoctest
julia> dv = [1, 2, 3, 4]
4-element Vector{Int64}:
 1
 2
 3
 4

julia> ev = [7, 8, 9]
3-element Vector{Int64}:
 7
 8
 9

julia> SymTridiagonal(dv, ev)
4×4 SymTridiagonal{Int64, Vector{Int64}}:
 1  7  ⋅  ⋅
 7  2  8  ⋅
 ⋅  8  3  9
 ⋅  ⋅  9  4

julia> A = SymTridiagonal(fill([1 2; 3 4], 3), fill([1 2; 3 4], 2));

julia> A[1,1]
2×2 Symmetric{Int64, Matrix{Int64}}:
 1  2
 2  4

julia> A[1,2]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> A[2,1]
2×2 Matrix{Int64}:
 1  3
 2  4
```
"""
SymTridiagonal(dv::DV, ev::EV) where {T,DV<:AbstractVector{T},EV<:AbstractVector{T}} = SymTridiagonal{T}(dv, ev)
SymTridiagonal{T}(dv::DV, ev::EV) where {T,DV<:AbstractVector{T},EV<:AbstractVector{T}} = SymTridiagonal{T,DV,EV}(dv, ev)
function SymTridiagonal{T}(dv::AbstractVector, ev::AbstractVector) where {T}
    SymTridiagonal(convert(AbstractVector{T}, dv)::AbstractVector{T},
                   convert(AbstractVector{T}, ev)::AbstractVector{T})
end

"""
    SymTridiagonal(A::AbstractMatrix)

Construct a symmetric tridiagonal matrix from the diagonal and first superdiagonal
of the symmetric matrix `A`.

# Examples
```jldoctest
julia> A = [1 2 3; 2 4 5; 3 5 6]
3×3 Matrix{Int64}:
 1  2  3
 2  4  5
 3  5  6

julia> SymTridiagonal(A)
3×3 SymTridiagonal{Int64, Vector{Int64}}:
 1  2  ⋅
 2  4  5
 ⋅  5  6

julia> B = reshape([[1 2; 2 3], [1 2; 3 4], [1 3; 2 4], [1 2; 2 3]], 2, 2);

julia> SymTridiagonal(B)
2×2 SymTridiagonal{Matrix{Int64}, Vector{Matrix{Int64}}}:
 [1 2; 2 3]  [1 3; 2 4]
 [1 2; 3 4]  [1 2; 2 3]
```
"""
function SymTridiagonal(A::AbstractMatrix)
    if (diag(A, 1) == transpose.(diag(A, -1))) && all(issymmetric.(diag(A, 0)))
        SymTridiagonal(diag(A, 0), diag(A, 1))
    else
        throw(ArgumentError("matrix is not symmetric; cannot convert to SymTridiagonal"))
    end
end

SymTridiagonal{T,DV,EV}(S::SymTridiagonal{T,DV,EV}) where {T,DV<:AbstractVector{T},EV<:AbstractVector{T}} = S
SymTridiagonal{T,DV,EV}(S::SymTridiagonal) where {T,DV<:AbstractVector{T},EV<:AbstractVector{T}} =
    SymTridiagonal(convert(DV, S.dv)::DV, convert(EV, S.ev)::EV)
SymTridiagonal{T}(S::SymTridiagonal{T}) where {T} = S
SymTridiagonal{T}(S::SymTridiagonal) where {T} =
    SymTridiagonal(convert(AbstractVector{T}, S.dv)::AbstractVector{T},
                   convert(AbstractVector{T}, S.ev)::AbstractVector{T})
SymTridiagonal(S::SymTridiagonal) = S
SymTridiagonal(dv::AbstractVector{U}, ev::AbstractVector{V}) where {U,V} = SymTridiagonal{promote_type(U,V)}(dv, ev)

LinearAlgebra.SymTridiagonal{T,V}(S::SymTridiagonal) where {T,V} = LinearAlgebra.SymTridiagonal{T,V}(convert(V,S.dv),convert(V,S.ev))

AbstractMatrix{T}(S::SymTridiagonal) where {T} =
    SymTridiagonal(convert(AbstractVector{T}, S.dv)::AbstractVector{T},
                   convert(AbstractVector{T}, S.ev)::AbstractVector{T})
function Matrix{T}(M::SymTridiagonal) where T
    n = size(M, 1)
    Mf = zeros(T, n, n)
    if n == 0
        return Mf
    end
    @inbounds begin
        @simd for i = 1:n-1
            Mf[i,i] = M.dv[i]
            Mf[i+1,i] = M.ev[i]
            Mf[i,i+1] = M.ev[i]
        end
        Mf[n,n] = M.dv[n]
    end
    return Mf
end
Matrix(M::SymTridiagonal{T}) where {T} = Matrix{T}(M)
Array(M::SymTridiagonal) = Matrix(M)

size(A::SymTridiagonal) = (length(A.dv), length(A.dv))
function size(A::SymTridiagonal, d::Integer)
    if d < 1
        throw(ArgumentError("dimension must be ≥ 1, got $d"))
    elseif d<=2
        return length(A.dv)
    else
        return 1
    end
end

# For S<:SymTridiagonal, similar(S[, neweltype]) should yield a SymTridiagonal matrix.
# On the other hand, similar(S, [neweltype,] shape...) should yield a sparse matrix.
# The first method below effects the former, and the second the latter.
similar(S::SymTridiagonal, ::Type{T}) where {T} = SymTridiagonal(similar(S.dv, T), similar(S.ev, T))
# The method below is moved to SparseArrays for now
# similar(S::SymTridiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...)

#Elementary operations
for func in (:conj, :copy, :real, :imag)
    @eval ($func)(M::SymTridiagonal) = SymTridiagonal(($func)(M.dv), ($func)(M.ev))
end

transpose(S::SymTridiagonal) = S
adjoint(S::SymTridiagonal{<:Real}) = S
adjoint(S::SymTridiagonal) = Adjoint(S)
Base.copy(S::Adjoint{<:Any,<:SymTridiagonal}) = SymTridiagonal(map(x -> copy.(adjoint.(x)), (S.parent.dv, S.parent.ev))...)
Base.copy(S::Transpose{<:Any,<:SymTridiagonal}) = SymTridiagonal(map(x -> copy.(transpose.(x)), (S.parent.dv, S.parent.ev))...)

function diag(M::SymTridiagonal{<:Number}, n::Integer=0)
    # every branch call similar(..., ::Int) to make sure the
    # same vector type is returned independent of n
    absn = abs(n)
    if absn == 0
        return copyto!(similar(M.dv, length(M.dv)), M.dv)
    elseif absn == 1
        return copyto!(similar(M.ev, length(M.ev)), M.ev)
    elseif absn <= size(M,1)
        return fill!(similar(M.dv, size(M,1)-absn), 0)
    else
        throw(ArgumentError(string("requested diagonal, $n, must be at least $(-size(M, 1)) ",
            "and at most $(size(M, 2)) for an $(size(M, 1))-by-$(size(M, 2)) matrix")))
    end
end
function diag(M::SymTridiagonal, n::Integer=0)
    # every branch call similar(..., ::Int) to make sure the
    # same vector type is returned independent of n
    if n == 0
        return copyto!(similar(M.dv, length(M.dv)), symmetric.(M.dv, :U))
    elseif n == 1
        return copyto!(similar(M.ev, length(M.ev)), M.ev)
    elseif n == -1
        return copyto!(similar(M.ev, length(M.ev)), transpose.(M.ev))
    elseif n <= size(M,1)
        throw(ArgumentError("requested diagonal contains undefined zeros of an array type"))
    else
        throw(ArgumentError(string("requested diagonal, $n, must be at least $(-size(M, 1)) ",
            "and at most $(size(M, 2)) for an $(size(M, 1))-by-$(size(M, 2)) matrix")))
    end
end

+(A::SymTridiagonal, B::SymTridiagonal) = SymTridiagonal(A.dv+B.dv, A.ev+B.ev)
-(A::SymTridiagonal, B::SymTridiagonal) = SymTridiagonal(A.dv-B.dv, A.ev-B.ev)
-(A::SymTridiagonal) = SymTridiagonal(-A.dv, -A.ev)
*(A::SymTridiagonal, B::Number) = SymTridiagonal(A.dv*B, A.ev*B)
*(B::Number, A::SymTridiagonal) = A*B
/(A::SymTridiagonal, B::Number) = SymTridiagonal(A.dv/B, A.ev/B)
\(B::Number, A::SymTridiagonal) = A/B
==(A::SymTridiagonal, B::SymTridiagonal) = (A.dv==B.dv) && (A.ev==B.ev)

#tril and triu

function istriu(M::SymTridiagonal, k::Integer=0)
    if k <= -1
        return true
    elseif k == 0
        return iszero(M.ev)
    else # k >= 1
        return iszero(M.ev) && iszero(M.dv)
    end
end
istril(M::SymTridiagonal, k::Integer) = istriu(M, -k)
iszero(M::SymTridiagonal) = iszero(M.ev) && iszero(M.dv)
isone(M::SymTridiagonal) = iszero(M.ev) && all(isone, M.dv)
isdiag(M::SymTridiagonal) = iszero(M.ev)

function tril!(M::SymTridiagonal, k::Integer=0)
    n = length(M.dv)
    if !(-n - 1 <= k <= n - 1)
        throw(ArgumentError(string("the requested diagonal, $k, must be at least ",
            "$(-n - 1) and at most $(n - 1) in an $n-by-$n matrix")))
    elseif k < -1
        fill!(M.ev,0)
        fill!(M.dv,0)
        return Tridiagonal(M.ev,M.dv,copy(M.ev))
    elseif k == -1
        fill!(M.dv,0)
        return Tridiagonal(M.ev,M.dv,zero(M.ev))
    elseif k == 0
        return Tridiagonal(M.ev,M.dv,zero(M.ev))
    elseif k >= 1
        return Tridiagonal(M.ev,M.dv,copy(M.ev))
    end
end

function triu!(M::SymTridiagonal, k::Integer=0)
    n = length(M.dv)
    if !(-n + 1 <= k <= n + 1)
        throw(ArgumentError(string("the requested diagonal, $k, must be at least ",
            "$(-n + 1) and at most $(n + 1) in an $n-by-$n matrix")))
    elseif k > 1
        fill!(M.ev,0)
        fill!(M.dv,0)
        return Tridiagonal(M.ev,M.dv,copy(M.ev))
    elseif k == 1
        fill!(M.dv,0)
        return Tridiagonal(zero(M.ev),M.dv,M.ev)
    elseif k == 0
        return Tridiagonal(zero(M.ev),M.dv,M.ev)
    elseif k <= -1
        return Tridiagonal(M.ev,M.dv,copy(M.ev))
    end
end

###################
# Generic methods #
###################

# Implements the determinant using principal minors
# a, b, c are assumed to be the subdiagonaldataonal, diagonal, and superdiagonal of
# a tridiagonal matrix.
#Reference:
#    R. Usmani, "Inversion of a tridiagonal Jacobi matrix",
#    Linear Algebra and its Applications 212-213 (1994), pp.413-414
#    doi:10.1016/0024-3795(94)90414-6
function det_usmani(a::V, b::V, c::V, shift::Number=0) where {T,V<:AbstractVector{T}}
    require_one_based_indexing(a, b, c)
    n = length(b)
    θa = oneunit(T)+zero(shift)
    if n == 0
        return θa
    end
    θb = b[1]+shift
    for i in 2:n
        θb, θa = (b[i]+shift)*θb - a[i-1]*c[i-1]*θa, θb
    end
    return θb
end

# det with optional diagonal shift for use with shifted Hessenberg factorizations
det(A::SymTridiagonal; shift::Number=false) = det_usmani(A.ev, A.dv, A.ev, shift)
logabsdet(A::SymTridiagonal; shift::Number=false) = logabsdet(ldlt(A; shift=shift))

function getindex(A::SymTridiagonal{T}, i::Integer, j::Integer) where T
    if !(1 <= i <= size(A,2) && 1 <= j <= size(A,2))
        throw(BoundsError(A, (i,j)))
    end
    if i == j
        return symmetric(A.dv[i], :U)::symmetric_type(eltype(A.dv))
    elseif i == j + 1
        return copy(transpose(A.ev[j])) # materialized for type stability
    elseif i + 1 == j
        return A.ev[i]
    else
        return zero(T)
    end
end

function setindex!(A::SymTridiagonal, x, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    if i == j
        @inbounds A.dv[i] = x
    else
        throw(ArgumentError("cannot set off-diagonal entry ($i, $j)"))
    end
    return x
end

## Tridiagonal matrices ##
struct Tridiagonal{T,DL<:AbstractVector{T},D<:AbstractVector{T},DU<:AbstractVector{T}} <: AbstractBandedMatrix{T}
    dl::DL    # sub-diagonal
    d::D     # diagonal
    du::DU   # sup-diagonal
    function Tridiagonal{T,DL,D,DU}(dl, d, du) where {T,DL<:AbstractVector{T},D<:AbstractVector{T},DU<:AbstractVector{T}}
        require_one_based_indexing(dl, d, du)
        n = length(d)
        if (length(dl) != n-1 || length(du) != n-1) && !(length(d) == 0 && length(dl) == 0 && length(du) == 0)
            throw(ArgumentError(string("cannot construct Tridiagonal from incompatible ",
                "lengths of subdiagonaldataonal, diagonal and superdiagonal: ",
                "($(length(dl)), $(length(d)), $(length(du)))")))
        end
        new{T,DL,D,DU}(dl, d, du)
    end
end

"""
    Tridiagonal(dl::DL, d::D, du::DU) where V <: AbstractVector

Construct a tridiagonal matrix from the first subdiagonaldataonal, diagonal, and first superdiagonal,
respectively. The result is of type `Tridiagonal` and provides efficient specialized linear
solvers, but may be converted into a regular matrix with
[`convert(Array, _)`](@ref) (or `Array(_)` for short).
The lengths of `dl` and `du` must be one less than the length of `d`.

# Examples
```jldoctest
julia> dl = [1, 2, 3];

julia> du = [4, 5, 6];

julia> d = [7, 8, 9, 0];

julia> Tridiagonal(dl, d, du)
4×4 Tridiagonal{Int64, Vector{Int64}}:
 7  4  ⋅  ⋅
 1  8  5  ⋅
 ⋅  2  9  6
 ⋅  ⋅  3  0
```
"""
Tridiagonal(dl::DL, d::D, du::DU) where {T,DL<:AbstractVector{T},D<:AbstractVector{T},DU<:AbstractVector{T}} = Tridiagonal{T,DL,D,DU}(dl, d, du)
Tridiagonal{T}(dl::AbstractVector, d::AbstractVector, du::AbstractVector) where {T} = Tridiagonal(map(x->convert(AbstractVector{T}, x), (dl, d, du))...)
Tridiagonal(dl::AbstractVector{T}, d::AbstractVector{U}, du::AbstractVector{V}) where {T,U,V} = Tridiagonal{promote_type(T,U,V)}(dl, d, du)

"""
    Tridiagonal(A)

Construct a tridiagonal matrix from the first sub-diagonal,
diagonal and first super-diagonal of the matrix `A`.

# Examples
```jldoctest
julia> A = [1 2 3 4; 1 2 3 4; 1 2 3 4; 1 2 3 4]
4×4 Matrix{Int64}:
 1  2  3  4
 1  2  3  4
 1  2  3  4
 1  2  3  4

julia> Tridiagonal(A)
4×4 Tridiagonal{Int64, Vector{Int64}}:
 1  2  ⋅  ⋅
 1  2  3  ⋅
 ⋅  2  3  4
 ⋅  ⋅  3  4
```
"""
Tridiagonal(A::AbstractMatrix) = Tridiagonal(diag(A,-1), diag(A,0), diag(A,1))

Tridiagonal(A::Tridiagonal) = A
Tridiagonal{T}(A::Tridiagonal{T}) where {T} = A
function Tridiagonal{T}(A::Tridiagonal) where {T}
    dl, d, du = map(x->convert(AbstractVector{T}, x)::AbstractVector{T},
                    (A.dl, A.d, A.du))
    Tridiagonal(dl, d, du)
end

size(M::Tridiagonal) = (length(M.d), length(M.d))
function size(M::Tridiagonal, d::Integer)
    if d < 1
        throw(ArgumentError("dimension d must be ≥ 1, got $d"))
    elseif d <= 2
        return length(M.d)
    else
        return 1
    end
end

function Matrix{T}(M::Tridiagonal{T}) where T
    A = zeros(T, size(M))
    for i = 1:length(M.d)
        A[i,i] = M.d[i]
    end
    for i = 1:length(M.d)-1
        A[i+1,i] = M.dl[i]
        A[i,i+1] = M.du[i]
    end
    A
end
Matrix(M::Tridiagonal{T}) where {T} = Matrix{T}(M)
Array(M::Tridiagonal) = Matrix(M)

# For M<:Tridiagonal, similar(M[, neweltype]) should yield a Tridiagonal matrix.
# On the other hand, similar(M, [neweltype,] shape...) should yield a sparse matrix.
# The first method below effects the former, and the second the latter.
similar(M::Tridiagonal, ::Type{T}) where {T} = Tridiagonal(similar(M.dl, T), similar(M.d, T), similar(M.du, T))
# The method below is moved to SparseArrays for now
# similar(M::Tridiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...)

# Operations on Tridiagonal matrices
copyto!(dest::Tridiagonal, src::Tridiagonal) = (copyto!(dest.dl, src.dl); copyto!(dest.d, src.d); copyto!(dest.du, src.du); dest)

#Elementary operations
for func in (:conj, :copy, :real, :imag)
    @eval function ($func)(M::Tridiagonal)
        Tridiagonal(($func)(M.dl), ($func)(M.d), ($func)(M.du))
    end
end

adjoint(S::Tridiagonal) = Adjoint(S)
transpose(S::Tridiagonal) = Transpose(S)
adjoint(S::Tridiagonal{<:Real}) = Tridiagonal(S.du, S.d, S.dl)
transpose(S::Tridiagonal{<:Number}) = Tridiagonal(S.du, S.d, S.dl)
Base.copy(aS::Adjoint{<:Any,<:Tridiagonal}) = (S = aS.parent; Tridiagonal(map(x -> copy.(adjoint.(x)), (S.du, S.d, S.dl))...))
Base.copy(tS::Transpose{<:Any,<:Tridiagonal}) = (S = tS.parent; Tridiagonal(map(x -> copy.(transpose.(x)), (S.du, S.d, S.dl))...))

function diag(M::Tridiagonal, n::Integer=0)
    # every branch call similar(..., ::Int) to make sure the
    # same vector type is returned independent of n
    if n == 0
        return copyto!(similar(M.d, length(M.d)), M.d)
    elseif n == -1
        return copyto!(similar(M.dl, length(M.dl)), M.dl)
    elseif n == 1
        return copyto!(similar(M.du, length(M.du)), M.du)
    elseif abs(n) <= size(M,1)
        return fill!(similar(M.d, size(M,1)-abs(n)), 0)
    else
        throw(ArgumentError(string("requested diagonal, $n, must be at least $(-size(M, 1)) ",
            "and at most $(size(M, 2)) for an $(size(M, 1))-by-$(size(M, 2)) matrix")))
    end
end

function getindex(A::Tridiagonal{T}, i::Integer, j::Integer) where T
    if !(1 <= i <= size(A,2) && 1 <= j <= size(A,2))
        throw(BoundsError(A, (i,j)))
    end
    if i == j
        return A.d[i]
    elseif i == j + 1
        return A.dl[j]
    elseif i + 1 == j
        return A.du[i]
    else
        return zero(T)
    end
end

function setindex!(A::Tridiagonal, x, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    if i == j
        @inbounds A.d[i] = x
    elseif i - j == 1
        @inbounds A.dl[j] = x
    elseif j - i == 1
        @inbounds A.du[i] = x
    elseif !iszero(x)
        throw(ArgumentError(string("cannot set entry ($i, $j) off ",
            "the tridiagonal band to a nonzero value ($x)")))
    end
    return x
end

#tril and triu

iszero(M::Tridiagonal) = iszero(M.dl) && iszero(M.d) && iszero(M.du)
isone(M::Tridiagonal) = iszero(M.dl) && all(isone, M.d) && iszero(M.du)
function istriu(M::Tridiagonal, k::Integer=0)
    if k <= -1
        return true
    elseif k == 0
        return iszero(M.dl)
    elseif k == 1
        return iszero(M.dl) && iszero(M.d)
    else # k >= 2
        return iszero(M.dl) && iszero(M.d) && iszero(M.du)
    end
end
function istril(M::Tridiagonal, k::Integer=0)
    if k >= 1
        return true
    elseif k == 0
        return iszero(M.du)
    elseif k == -1
        return iszero(M.du) && iszero(M.d)
    else # k <= -2
        return iszero(M.du) && iszero(M.d) && iszero(M.dl)
    end
end
isdiag(M::Tridiagonal) = iszero(M.dl) && iszero(M.du)

function tril!(M::Tridiagonal, k::Integer=0)
    n = length(M.d)
    if !(-n - 1 <= k <= n - 1)
        throw(ArgumentError(string("the requested diagonal, $k, must be at least ",
            "$(-n - 1) and at most $(n - 1) in an $n-by-$n matrix")))
    elseif k < -1
        fill!(M.dl,0)
        fill!(M.d,0)
        fill!(M.du,0)
    elseif k == -1
        fill!(M.d,0)
        fill!(M.du,0)
    elseif k == 0
        fill!(M.du,0)
    end
    return M
end

function triu!(M::Tridiagonal, k::Integer=0)
    n = length(M.d)
    if !(-n + 1 <= k <= n + 1)
        throw(ArgumentError(string("the requested diagonal, $k, must be at least ",
            "$(-n + 1) and at most $(n + 1) in an $n-by-$n matrix")))
    elseif k > 1
        fill!(M.dl,0)
        fill!(M.d,0)
        fill!(M.du,0)
    elseif k == 1
        fill!(M.dl,0)
        fill!(M.d,0)
    elseif k == 0
        fill!(M.dl,0)
    end
    return M
end

###################
# Generic methods #
###################

-(A::Tridiagonal) = Tridiagonal(-A.dl, -A.d, -A.du)
+(A::Tridiagonal, B::Tridiagonal) = Tridiagonal(A.dl+B.dl, A.d+B.d, A.du+B.du)
-(A::Tridiagonal, B::Tridiagonal) = Tridiagonal(A.dl-B.dl, A.d-B.d, A.du-B.du)
*(A::Tridiagonal, B::Number) = Tridiagonal(A.dl*B, A.d*B, A.du*B)
*(B::Number, A::Tridiagonal) = A*B
/(A::Tridiagonal, B::Number) = Tridiagonal(A.dl/B, A.d/B, A.du/B)
\(B::Number, A::Tridiagonal) = A/B

==(A::Tridiagonal, B::Tridiagonal) = (A.dl==B.dl) && (A.d==B.d) && (A.du==B.du)
==(A::Tridiagonal, B::SymTridiagonal) = (A.dl==A.du==B.ev) && (A.d==B.dv)
==(A::SymTridiagonal, B::Tridiagonal) = (B.dl==B.du==A.ev) && (B.d==A.dv)

det(A::Tridiagonal) = det_usmani(A.dl, A.d, A.du)

AbstractMatrix{T}(M::Tridiagonal) where {T} = Tridiagonal{T}(M)
Tridiagonal{T}(M::SymTridiagonal{T}) where {T} = Tridiagonal(M)
function SymTridiagonal{T}(M::Tridiagonal) where T
    if M.dl == M.du
        return SymTridiagonal{T}(convert(AbstractVector{T},M.d), convert(AbstractVector{T},M.dl))
    else
        throw(ArgumentError("Tridiagonal is not symmetric, cannot convert to SymTridiagonal"))
    end
end

Base._sum(A::Tridiagonal, ::Colon) = sum(A.d) + sum(A.dl) + sum(A.du)
Base._sum(A::SymTridiagonal, ::Colon) = sum(A.dv) + 2sum(A.ev)

function Base._sum(A::Tridiagonal, dims::Integer)
    res = Base.reducedim_initarray(A, dims, zero(eltype(A)))
    n = length(A.d)
    if n == 0
        return res
    elseif n == 1
        res[1] = A.d[1]
        return res
    end
    @inbounds begin
        if dims == 1
            res[1] = A.dl[1] + A.d[1]
            for i = 2:n-1
                res[i] = A.dl[i] + A.d[i] + A.du[i-1]
            end
            res[n] = A.d[n] + A.du[n-1]
        elseif dims == 2
            res[1] = A.d[1] + A.du[1]
            for i = 2:n-1
                res[i] = A.dl[i-1] + A.d[i] + A.du[i]
            end
            res[n] = A.dl[n-1] + A.d[n]
        elseif dims >= 3
            for i = 1:n-1
                res[i,i+1] = A.du[i]
                res[i,i]   = A.d[i]
                res[i+1,i] = A.dl[i]
            end
            res[n,n] = A.d[n]
        end
    end
    res
end

function Base._sum(A::SymTridiagonal, dims::Integer)
    res = Base.reducedim_initarray(A, dims, zero(eltype(A)))
    n = length(A.dv)
    if n == 0
        return res
    elseif n == 1
        res[1] = A.dv[1]
        return res
    end
    @inbounds begin
        if dims == 1
            res[1] = A.ev[1] + A.dv[1]
            for i = 2:n-1
                res[i] = A.ev[i] + A.dv[i] + A.ev[i-1]
            end
            res[n] = A.dv[n] + A.ev[n-1]
        elseif dims == 2
            res[1] = A.dv[1] + A.ev[1]
            for i = 2:n-1
                res[i] = A.ev[i-1] + A.dv[i] + A.ev[i]
            end
            res[n] = A.ev[n-1] + A.dv[n]
        elseif dims >= 3
            for i = 1:n-1
                res[i,i+1] = A.ev[i]
                res[i,i]   = A.dv[i]
                res[i+1,i] = A.ev[i]
            end
            res[n,n] = A.dv[n]
        end
    end
    res
end

function dot(x::AbstractVector, A::Tridiagonal, y::AbstractVector)
    require_one_based_indexing(x, y)
    nx, ny = length(x), length(y)
    (nx == size(A, 1) == ny) || throw(DimensionMismatch())
    if iszero(nx)
        return dot(zero(eltype(x)), zero(eltype(A)), zero(eltype(y)))
    end
    x₀ = x[1]
    x₊ = x[2]
    dl, d, du = A.dl, A.d, A.du
    r = dot(adjoint(d[1])*x₀ + adjoint(dl[1])*x₊, y[1])
    @inbounds for j in 2:nx-1
        x₋, x₀, x₊ = x₀, x₊, x[j+1]
        r += dot(adjoint(du[j-1])*x₋ + adjoint(d[j])*x₀ + adjoint(dl[j])*x₊, y[j])
    end
    r += dot(adjoint(du[nx-1])*x₀ + adjoint(d[nx])*x₊, y[nx])
    return r
end

MemoryLayout(::Type{<:SymTridiagonal{<:Any,DV,EV}}) where {DV,EV} = symtridiagonallayout(MemoryLayout(DV), MemoryLayout(EV))
MemoryLayout(::Type{<:Tridiagonal{<:Any,DL,D,DU}}) where {DL,D,DU} = tridiagonallayout(MemoryLayout(DL), MemoryLayout(D), MemoryLayout(DU))
bandwidths(::SymTridiagonal) = (1,1)
bandwidths(::Tridiagonal) = (1,1)

Base.BroadcastStyle(::Type{SymTridiagonal{T,DV,EV}}) where {T,DV,EV} = StructuredMatrixStyle{LinearAlgebra.SymTridiagonal{T,DV}}()
Base.BroadcastStyle(::Type{Tridiagonal{T,DL,D,DU}}) where {T,DL,D,DU} = StructuredMatrixStyle{LinearAlgebra.Tridiagonal{T,D}}()

convert(::Type{LinearAlgebra.Tridiagonal}, B::Tridiagonal) = LinearAlgebra.Tridiagonal(B.dl, B.d, B.du)
convert(::Type{LinearAlgebra.SymTridiagonal}, B::SymTridiagonal) = LinearAlgebra.SymTridiagonal(B.dv, B.ev)

permutedims(S::SymTridiagonal) = S
permutedims(T::Tridiagonal) = Tridiagonal(T.du, T.d, T.dl)

diagonaldata(D::SymTridiagonal) = D.dv
diagonaldata(D::Tridiagonal) = D.d

supdiagonaldata(D::SymTridiagonal) = D.ev
subdiagonaldata(D::SymTridiagonal) = D.ev

subdiagonaldata(D::Tridiagonal) = D.dl
supdiagonaldata(D::Tridiagonal) = D.du

function BlockArrays.sizes_from_blocks(A::Tridiagonal, _)
    # for k = 1:length(A.du)
    #     size(A.du[k],1) == sz[1][k] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
    #     size(A.du[k],2) == sz[2][k+1] || throw(ArgumentError("block sizes of upper diagonal inconsisent with diagonal"))
    #     size(A.dl[k],1) == sz[1][k+1] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
    #     size(A.dl[k],2) == sz[2][k] || throw(ArgumentError("block sizes of lower diagonal inconsisent with diagonal"))
    # end
    (size.(A.d, 1), size.(A.d,2))
end

eigvals(A::SymTridiagonal{T}) where T = eigvals(LinearAlgebra.SymTridiagonal{T,Vector{T}}(A))
eigen(A::SymTridiagonal{T}) where T = eigen(LinearAlgebra.SymTridiagonal{T,Vector{T}}(A))


text:
# on-device memory management

export @LocalMemory

## thread-local array

"""
Creates a local static memory shared inside one block.
Equivalent to `__local` of OpenCL or `__shared__ (<variable>)` of CUDA.
"""
macro LocalMemory(ctx, T, N)
    id = gensym("local_memory")
    quote
        LocalMemory($(esc(ctx)), $(esc(T)), Val($(esc(N))), Val($(QuoteNode(id))))
    end
end

"""
Creates a block local array pointer with `T` being the element type
and `N` the length. Both T and N need to be static! C is a counter for
approriately get the correct Local mem id in CUDAnative.
This is an internal method which needs to be overloaded by the GPU Array backends
"""
function LocalMemory(ctx, ::Type{T}, ::Val{dims}, ::Val{id}) where {T, dims, id}
    error("Not implemented") # COV_EXCL_LINE
end


text:
export i_mean_sum, i_var_mean_sum, i_normal_logpdf, i_cor_cov

"""
    i_mean_sum(out!, sum!, x)

get the `mean` and `sum` of `x`.
"""
@i function i_mean_sum(out!, sum!, x)
    for i=1:length(x)
        sum! += x[i]
    end
    out! += sum!/length(x)
end

"""
    i_var_mean_sum(var!, varsum!, mean!, sum!, sqv)

Compute the variance, the accumulated variance, mean and sum.
"""
@i function i_var_mean_sum(var!, varsum!, mean!, sum!, v::AbstractVector{T}) where T
    i_mean_sum(mean!, sum!, v)
    for i=1:length(v)
        v[i] -= mean!
        varsum! += v[i] ^ 2
        v[i] += mean!
    end
    var! += varsum! / (length(v)-1)
end

"""
    i_normal_logpdf(out, x, μ, σ)

get the pdf of `Normal(μ, σ)` at point `x`.
"""
@i function i_normal_logpdf(out, x::T, μ, σ) where T
    @zeros T anc1 anc2 anc3

    @routine begin
        anc1 += x
        anc1 -= μ
        anc2 += anc1 / σ  # (x- μ)/σ
        anc3 += anc2^2 # (x-μ)^2/σ^2
    end

    out -= anc3 * 0.5 # -(x-μ)^2/2σ^2
    out -= log(σ) # -(x-μ)^2/2σ^2 - log(σ)
    out -= log(2π)/2 # -(x-μ)^2/2σ^2 - log(σ) - log(2π)/2

    ~@routine
end

"""
     i_cor_cov(rho!,cov!,a,b)

get Pearson correlation and covariance of two vectors `a` and `b`

"""

@i function i_cor_cov(rho!::T,cov!::T,a::AbstractVector{T},b::AbstractVector{T}) where T
    @safe @assert length(a) == length(b)
    @routine  @invcheckoff begin
        @zeros T var1 varsum1 mean1 sum1 std1
        i_var_mean_sum(var1, varsum1, mean1, sum1, a)
        std1 += sqrt(var1)
        @zeros T var2 varsum2 mean2 sum2 std2
        i_var_mean_sum(var2, varsum2, mean2, sum2, b)
        std2 += sqrt(var2)
        @zeros T anc3 anc4 anc5 anc6 anc7
        @inbounds for i=1:length(b)
            a[i] -= mean1
            anc3 += a[i]
            b[i] -= mean2
            anc4 += b[i]
            anc5 += anc3*anc4
            anc3 -= a[i]
            anc4 -= b[i]
            a[i] += mean1
            b[i] += mean2
        end
        anc6+=std1*std2
        anc7+=anc6*(length(b)-1)
    end
        cov! += anc5/(length(b)-1)
        rho! += anc5/anc7
    ~@routine
end


text:
export DiaSymSemiseparable

struct DiaSymSemiseparable <: SymSemiseparableMatrix
    n::Int64
    p::Int64
    U::AbstractArray
    V::AbstractArray
    d::AbstractArray
end

# Constuctors
function DiaSymSemiseparable(U::AbstractArray, V::AbstractArray, d::AbstractArray)
    if size(U,1) == size(V,1) && size(U,2) == size(V,2) && length(d) == size(U,1)
        return DiaSymSemiseparable(size(U,1),size(U,2),U,V,d);
    else
        error("Dimension mismatch between the generators U, V and d")
    end
end
# Mappings
mul!(y::AbstractArray, L::DiaSymSemiseparable, 		            b::AbstractArray) = dss_mul_mat!(y, L.U, L.V, L.d, b);
mul!(y::AbstractArray, L::AdjointOperator{DiaSymSemiseparable}, b::AbstractArray) = dss_mul_mat!(y, L.A.U, L.A.V, L.A.d, b);
function inv!(y, K::DiaSymSemiseparable, b::AbstractArray)
	L = DiaSymSemiseparableChol(K)
	y[:,:] = L'\(L\b)
end
function inv!(y, K::AdjointOperator{DiaSymSemiseparable}, b::AbstractArray)
	L = DiaSymSemiseparableChol(K.A)
	y[:,:] = L'\(L\b)
end

#################################################
#### Cholesky factoriaztion of:              ####
####    Higher-order quasiseparable matrices ####
#################################################

#### Matrix-matrix product ####
function dss_mul_mat!(Y::Array, U::Array, V::Array, d::Array, X::Array)
    n, m = size(U);
    mx = size(X,2);
    Vbar = zeros(m,mx);
    Ubar = U'*X;
    for i = 1:n
        tmpV = V[i,:];
        tmpU = U[i,:];
        tmpX = X[i:i,:];
        Ubar -= tmpU .* tmpX;
        Vbar += tmpV .* tmpX;
        Y[i,:] = tmpU'*Vbar + tmpV'*Ubar + d[i]*tmpX;
    end
end
#### Going from the Cholesky factorization back to matrix Σ ####
function dss_create_vd(U::Array, W::Array, dbar::Array)
    n,m = size(U)
    d = zeros(n)
    V = zeros(n,m)
    P = zeros(m,m)
    for i = 1:n
        tmpU = U[i,:]
        tmpW = W[i,:]
        d[i] = dbar[i]^2 - tmpU'*tmpW
        V[i,:] = P*tmpU
        P += tmpW*tmpW';
    end
    return V, d
end


text:
struct Hull{T<:AbstractVector}
    vertices::Vector{T}
    k::Int
    converged::Bool
end

Hull(::Type{T},k::Int) where {T<:AbstractVector} = Hull(T[],k, false)

function is_left(p0,p1,p2)
    return ((p1[1] - p0[1]) * (p2[2] - p0[2]) - (p2[1] -  p0[1]) * (p1[2] - p0[2]))
end

function in_hull(p,hull::Hull)

    v = copy(hull.vertices)
    push!(v,hull.vertices[1])
    nv = length(v)
    wn = 0
    @inbounds for i = 1:nv-1
        if v[i][2] <= p[2]
            if v[i+1][2] > p[2]
                if is_left(v[i],v[i+1],p) > 0.0
                    wn = wn +1
                end
            end
        else
            if v[i+1][2] <= p[2]
                if is_left(v[i],v[i+1],p) < 0.0
                    wn = wn - 1
                end
            end
        end
    end
    if wn == 0
        return false
    else
        return true
    end
end

function signed_area(h::Hull)
    v = h.vertices
    if v[end] == v[1]
        n = length(v)-1
    else
        n = length(v)
    end

    n < 3 && return 0.0

    A1 = 0.0
    @inbounds for i=1:n-1
        A1 = A1 + v[i][1]*v[i+1][2]
    end
    A1 = A1 + v[n][1]*v[1][2]

    A2 = 0.0
    @inbounds for i=1:n-1
        A2 = A2 + v[i+1][1]*v[i][2]
    end
    A2 = A2 + v[1][1]*v[n][2]

    return (A1 - A2)/2
end

area(h::Hull) = abs(signed_area(h))

orientation(h::Hull) = sign(signed_area(h))

function Base.show(io::IO, hull::Hull)
    println(io, typeof(hull))
    println(io, "  Number of vertices: ", length(hull.vertices))
    println(io, "  Number of neighbors used: ", hull.k)
end


text:
# Goal: Make sure that our solvers are "correct" for small steps or tolerances
# Verify this for many (ideally all) combinations
# Compare with an algorithm from OrdinaryDiffEq.jl with high precision
using ProbNumDiffEq
using Test
using OrdinaryDiffEq
using LinearAlgebra
using Statistics: mean
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems;
importodeproblems();
import DiffEqProblemLibrary.ODEProblemLibrary:
    prob_ode_lotkavoltera, prob_ode_fitzhughnagumo

import ProbNumDiffEq: remake_prob_with_jac

EK1FDB1(; kwargs...) = EK1FDB(; jac_quality=1, kwargs...)
EK1FDB2(; kwargs...) = EK1FDB(; jac_quality=2, kwargs...)
EK1FDB3(; kwargs...) = EK1FDB(; jac_quality=3, kwargs...)
ALGS = [EK0, EK1, EK1FDB1, EK1FDB2, EK1FDB3]
DIFFUSIONS =
    [FixedDiffusion(), DynamicDiffusion(), FixedMVDiffusion(), DynamicMVDiffusion()]
INITS = [TaylorModeInit(), ClassicSolverInit()]

for (prob, probname) in [
    (remake_prob_with_jac(prob_ode_lotkavoltera), "lotkavolterra"),
    (remake_prob_with_jac(prob_ode_fitzhughnagumo), "fitzhughnagumo"),
]
    @testset "Constant steps: $probname" begin
        true_sol =
            solve(remake(prob, u0=big.(prob.u0)), Tsit5(), abstol=1e-20, reltol=1e-20)

        for Alg in ALGS, diffusion in DIFFUSIONS, init in INITS, q in [2, 3, 5]
            if (diffusion isa FixedMVDiffusion || diffusion isa DynamicMVDiffusion) &&
               Alg != EK0
                continue
            end

            @testset "Constant steps: $probname; alg=$Alg, diffusion=$diffusion, init=$init, q=$q" begin
                @debug "Testing for correctness: Constant steps" probname alg diffusion q dt

                sol = solve(
                    prob,
                    Alg(order=q, diffusionmodel=diffusion),
                    adaptive=false,
                    dt=5e-3,
                )
                @test sol.u ≈ true_sol.(sol.t) rtol = 1e-5
            end
        end
    end
end

for (prob, probname) in [
    (remake_prob_with_jac(prob_ode_lotkavoltera), "lotkavolterra"),
    (remake_prob_with_jac(prob_ode_fitzhughnagumo), "fitzhughnagumo"),
]
    @testset "Adaptive steps: $probname" begin
        t_eval = prob.tspan[1]:0.01:prob.tspan[end]
        true_sol =
            solve(remake(prob, u0=big.(prob.u0)), Tsit5(), abstol=1e-20, reltol=1e-20)
        true_dense_vals = true_sol.(t_eval)

        for Alg in ALGS, diffusion in DIFFUSIONS, init in INITS, q in [2, 3, 5]
            if (diffusion isa FixedMVDiffusion || diffusion isa DynamicMVDiffusion) &&
               Alg != EK0
                continue
            end

            @testset "Adaptive steps: $probname; alg=$Alg, diffusion=$diffusion, init=$init, q=$q" begin
                @debug "Testing for correctness: Adaptive steps" probname Alg diffusion q

                sol = solve(prob, Alg(order=q, diffusionmodel=diffusion), adaptive=true)

                @test sol.u ≈ true_sol.(sol.t) rtol = 1e-3
                @test mean.(sol.(t_eval)) ≈ true_dense_vals rtol = 1e-3
            end
        end
    end
end


text:
using Yao
using QuAlgorithmZoo
using YaoBlocks.ConstGate
using Test
using LinearAlgebra

@testset "hadamard test" begin
    n = 2
    U = chain(put(n, 2=>Rx(0.2)), put(n, 1=>Rz(1.2)), put(n, 1=>phase(0.4)))
    US = chain(2n, put(2n, (3,4)=>U),
            chain(2n, [swap(2n,i,i+n) for i=1:n]))
    reg = zero_state(2n)
    reg |> repeat(2n, H, 1:n) |> chain(2n, [cnot(2n,i,i+n) for i=1:n])
    @show real(tr(mat(U)))/(1<<n)
    @test hadamard_test(US, reg, 0.0) ≈ real(tr(mat(U)))/(1<<n)
    @test hadamard_test(US, reg, -π/2) ≈ imag(tr(mat(U)))/(1<<n)
end


text:
@testset "Quadratic and hermitian forms" begin
  include("QuadForm/Basic.jl")
  include("QuadForm/Lattices.jl")
  include("QuadForm/Genus.jl")
end


text:
using Chemicals
using Test

@testset "Chemicals.jl" begin
    # Write your tests here.
end


text:
import Base.iterate

function iterate(x::AbstractVariable, s = 0)
    return s >= length(x) ? nothing : (x[s+1], s + 1)
end


text:
################################################################################
####################### Approximation error plot ###############################
################################################################################
function approx_error2(DVEC::Array{Array{Float64,1},1},
                       T::Vector{String}, L::Vector{Tuple{Symbol,Symbol}},
                       frac::Float64 = 0.30)
    # This version plots the relative L2 errors against
    # the FRACTION of coefficients retained.
    plot(xaxis = "Fraction of Coefficients Retained",
         yaxis = "Relative Approximation Error")
    for i = 1:length(DVEC)
        dvec = DVEC[i]
        N = length(dvec)
        dvec_norm = norm(dvec,2)
        dvec_sort = sort(dvec.^2) # the smallest first
        er = sqrt.(reverse(cumsum(dvec_sort)))/dvec_norm
        er[er .== 0.0] .= minimum(er[er .!= 0.0]) # avoid blowup by taking log in the plot below
        # er is the relative L^2 error of the whole thing: length(er)=N.
        p = Int64(floor(frac*N)) + 1 # upper limit
        plot!(frac*(0:(p-1))/(p-1), er[1:p], yaxis=:log, xlims = (0.,frac),
              label = T[i], line = L[i], linewidth = 2, grid = false)
    end
    display(current())
end

function approx_error3(DVEC::Array{Array{Float64,1},1},
                       T::Vector{String}, L::Vector{Tuple{Symbol,Symbol}},
                       N::Int64)
    # This version plots the relative L2 errors against
    # the NUMBER of coefficients retained.
    plot(xaxis = "Number of Coefficients Retained",
         yaxis = "Relative Approximation Error")
    for i = 1:length(DVEC)
        dvec = DVEC[i]
        dvec_norm = norm(dvec,2)
        dvec_sort = sort(dvec.^2) # the smallest first
        er = sqrt.(reverse(cumsum(dvec_sort)))/dvec_norm
        er[er .== 0.0] .= minimum(er[er .!= 0.0]) # avoid blowup by taking log in the plot below
        # er is the relative L^2 error of the whole thing: length(er)=N.
        plot!(0:N-1, er[1:N], yaxis=:log, label = T[i], line = L[i],
              linewidth = 2, grid = false)
    end
    display(current())
end

################################################################################
############ Computing a weight matrix using the Gaussian affinity #############
################################################################################
function image_Gaussian_affinity(img::Matrix{Float64}, r::Int64, σ::Float64)
# Get the neighbors for affinity matrix computation
    l = 2*r + 1
    temp_x, temp_y = fill(1,l^2), fill(1,l^2)
    temp_xy = CartesianIndices((1:l,1:l))
    for i = 1:l^2
        temp_x[i], temp_y[i] = temp_xy[i][1], temp_xy[i][2]
    end
    # (r+1, r+1) is the index of the center location.
    temp_ind = ((temp_x .- (r + 1)).^2 + (temp_y .- (r + 1)).^2) .<= r^2
    # Now, temp_ind indicates those points withinin the circle of radius r
    # from the center while neighbor_x, neighbor_y represent relative positions
    # of those points from the center.
    neighbor_x = temp_x[temp_ind] .- (r + 1)
    neighbor_y = temp_y[temp_ind] .- (r + 1)
   # So, for any index (x, y), points within (x ± neighbor_x, y ± neighbor_y) are
   # its neighbors for the purpose of calculating affinity

# Create affinity matrix
    m, n = size(img)
    sig = img[:]
    W = fill(0., (m*n,m*n))
    for i = 1:m*n
        cur = CartesianIndices((m,n))[i]
        for j = 1:length(neighbor_x) # assuming dim(neighbor_x) == dim(neighbor_y)
            if 1 <= cur[1] + neighbor_x[j] <= m && 1 <= cur[2] + neighbor_y[j] <= n
                tempd = LinearIndices((m,n))[cur[1] + neighbor_x[j], cur[2] + neighbor_y[j]]
                W[i,tempd] = exp(-(sig[i] - sig[tempd])^2/σ)
            end
        end
    end
    return sparse(W)
end # end of the image_Gaussian_affinity function

################################################################################
######## Display top 9 basis vectors of various bases for an image data ########
################################################################################
function top_vectors_plot2(dvec::Array{Float64, 2}, m::Int64, n::Int64, BS::BasisSpec, GP::GraphPart;
                           clims::Tuple{Float64,Float64} = (-0.01, 0.01))

    # Get the indices from max to min in terms of absolute value of the coefs
    # Note that m*n == length(dvec); m and n are the image size.
    sorted_ind = sortperm(abs.(dvec[:]), rev = true);

    # Set up the layout as 3 x 3 subplots
    plot(9, layout = (3,3), framestyle = :none, legend = false)

    # Do the display!
    for i=1:9
        dvecT = fill(0., size(dvec))
        dvecT[sorted_ind[i]] = 1
        f = ghwt_synthesis(dvecT, GP, BS)
        heatmap!(reshape(f, (m,n)), subplot=i, ratio=1, yaxis=:flip,
                 showaxis=false, ticks = false, color = :grays, clims = clims)
    end
    display(current())
end
################################################################################

function top_vectors_plot3(dvec::Array{Float64, 2}, m::Int64, n::Int64, BS::BasisSpec,
    GP::GraphPart; clims::Tuple{Float64,Float64} = (-0.01, 0.01), K::Int64 = 9)

    # Get the indices from max to min in terms of absolute value of the coefs
    # Note that m*n == length(dvec); m and n are the image size.
    sorted_ind = sortperm(abs.(dvec[:]), rev = true);

    # Set up the layout as 3 x 3 subplots
    plot(K, layout = (Int(sqrt(K)), Int(sqrt(K))), framestyle = :none, legend = false)

    # Do the display!
    for i=1:K
        dvecT = fill(0., size(dvec))
        dvecT[sorted_ind[i]] = 1
        f = ghwt_synthesis(dvecT, GP, BS)
        heatmap!(reshape(f, (m,n)), subplot=i, ratio=1, yaxis=:flip,
                 showaxis=false, ticks = false, color = :grays, clims = clims)
    end
    display(current())
end
################################################################################


text:
"""
	options = NewtonPar(tol = 1e-4,...)

Returns a variable containing parameters to affect the `newton` algorithm when solving `F(x) = 0`.

# Arguments (with default values):
- `tol = 1e-10`: absolute tolerance for `F(x)`
- `maxIter = 50`: number of Newton iterations
- `verbose = false`: display Newton iterations?
- `linsolver = DefaultLS()`: linear solver, must be `<: AbstractLinearSolver`
- `eigsolver = DefaultEig()`: eigen solver, must be `<: AbstractEigenSolver`

# Arguments only used in `newtonPALC`
- `linesearch = false`: use line search algorithm
- `alpha = 1.0`: alpha (damping) parameter for line search algorithm
- `almin  = 0.001 `: minimal vslue of the damping `alpha`

!!! tip "Mutating"
    For performance reasons, we decided to use an immutable structure to hold the parameters. One can use the package `Setfield.jl` to drastically simplify the mutation of different fields. See the tutorials for examples.
"""
@with_kw struct NewtonPar{T, L <: AbstractLinearSolver, E <: AbstractEigenSolver}
	tol::T			 = 1e-10
	maxIter::Int64 	 = 50
	alpha::T         = convert(typeof(tol), 1.0)        # damping
	almin::T         = convert(typeof(tol), 0.001)      # minimal damping
	verbose::Bool    = false
	linesearch::Bool = false
	linsolver::L 	 = DefaultLS()
	eigsolver::E 	 = DefaultEig()
end

####################################################################################################
"""
		newton(F, J, x0, p0, options::NewtonPar; normN = norm, callback = (x, f, J, res, iteration, itlinear, optionsN; kwargs...) -> true, kwargs...)

This is the Newton-Krylov Solver for `F(x, p0) = 0` with Jacobian w.r.t. `x` written `J(x, p0)` and initial guess `x0`. The function `normN` allows to specify a norm for the convergence criteria. It is important to set the linear solver `options.linsolver` properly depending on your problem. This linear solver is used to solve ``J(x, p0)u = -F(x, p0)`` in the Newton step. You can for example use `linsolver = DefaultLS()` which is the operator backslash: it works well for Sparse / Dense matrices. See [Linear solvers](@ref) for more informations.

# Arguments:
- `(x, p) -> F(x, p)` functional whose zeros are looked for. In particular, it is not **inplace**,
- `dF(x, p) = (x, p) -> J(x, p)` compute the jacobian of `F` at `x`. It is then passed to `options.linsolver`. The Jacobian `J(x, p)` can be a matrix or an out-of-place function.
- `x0` initial guess
- `p0` set of parameters to be passed to `F` and `J`
- `options` variable holding the internal parameters used by the `newton` method
- `callback` function passed by the user which is called at the end of each iteration. Can be used to update a preconditionner for example. The arguments passed to the callback are as follows
    - `x` current solution
    - `f` current residual
    - `J` current jacobian
    - `res` current norm of the residual
    - `iteration` current newton iteration
    - `itlinear` number of iterations to solve the linear system
    - `optionsN` a copy of the argument `options` passed to `newton`
    - `kwargs` kwargs arguments, contain your initial guess `x0`
- `kwargs` arguments passed to the callback. Useful when `newton` is called from `continuation`

# Output:
- solution:
- history of residuals
- flag of convergence
- number of iterations

# Simplified calls
When `J` is not passed. It is then computed with finite differences. The call is as follows:

	newton(Fhandle, x0, p0, options::NewtonPar; kwargs...)

You can also pass functions which do not have parameters `x ->F(x)`, `x->J(x)` as follows

	newton(F, J, x0, options::NewtonPar;  kwargs...)

or

	newton(F, x0, options::NewtonPar;  kwargs...)

# Example

```
julia> F(x, p) = x.^3 .- 1
julia> Jac(x, p) = spdiagm(0 => 3 .* x.^2) # sparse jacobian
julia> x0 = rand(1_000)
julia> opts = NewtonPar()
julia> sol, hist, flag, _ = newton(F, Jac, x0, nothing, opts, normN = x->norm(x, Inf))
```

!!! tip "Other formulation"
    If you don't have parameters, you can still use `newton` as follows `newton((x,p) -> F(x), (x,p)-> J(x), x0, nothing, options)`
"""
function newton(Fhandle, Jhandle, x0, p0, options::NewtonPar; normN = norm, callback = (x, f, J, res, iteration, itlinear, optionsN; kwargs...) -> true, kwargs...)
	# Extract parameters
	@unpack tol, maxIter, verbose, linesearch = options

	# Initialize iterations
	x = similar(x0); copyto!(x, x0) # x = copy(x0)
	f = Fhandle(x, p0)
	d = similar(f); copyto!(d, f)	# d = copy(f)

	neval = 1
	res = normN(f)
	resHist = [res]
	it = 0

	# Displaying results
	verbose && displayIteration(it, neval, res)

	# invoke callback before algo really starts
	compute = callback(x, f, nothing, res, it, 0, options; x0 = x0, kwargs...)

	# Main loop
	while (res > tol) & (it < maxIter) & compute
		J = Jhandle(x, p0)
		d, _, itlinear = options.linsolver(J, f)

		# Update solution: x .= x .- d
		minus!(x, d)

		copyto!(f, Fhandle(x, p0))
		res = normN(f)

		neval += 1
		push!(resHist, res)
		it += 1

		verbose && displayIteration(it, neval, res, itlinear)

		if callback(x, f, J, res, it, itlinear, options; x0 = x0, kwargs...) == false
			break
		end
	end
	((resHist[end] > tol) && verbose) && @error("\n--> Newton algorithm failed to converge, residual = $(res[end])")
	flag = (resHist[end] < tol) & callback(x, f, nothing, res, it, nothing, options; x0 = x0, kwargs...)
	return x, resHist, flag, it
end

# simplified call to newton when no Jacobian is passed in which case we estimate it using finiteDifferences
function newton(Fhandle, x0, p0, options::NewtonPar; kwargs...)
	Jhandle = (u, p) -> finiteDifferences(z -> Fhandle(z, p), u)
	return newton(Fhandle, Jhandle, x0, p0, options; kwargs...)
end


text:
abstract type AbstractBuffer end
abstract type AbstractNetworkBuffer <: AbstractBuffer end

"Store indexing on target device"
struct IndexingCache{IVT} <: AbstractBuffer
    index_pv::IVT
    index_pq::IVT
    index_ref::IVT
    index_generators::IVT
    index_pv_to_gen::IVT
    index_ref_to_gen::IVT
end

"""
    PolarNetworkState{VI, VT} <: AbstractNetworkBuffer

Buffer to store current values of all the variables describing
the network, in polar formulation. Attributes are:

- `vmag` (length: nbus): voltage magnitude at each bus
- `vang` (length: nbus): voltage angle at each bus
- `pinj` (length: nbus): power injection RHS. Equal to `Cg * Pg - Pd`
- `qinj` (length: nbus): power injection RHS. Equal to `Cg * Qg - Qd`
- `pg`   (length: ngen): active power of generators
- `qg`   (length: ngen): reactive power of generators
- `dx`   (length: nstates): cache the difference between two consecutive states (used in power flow resolution)
- `balance` (length: nstates): cache for current power imbalance (used in power flow resolution)
- `bus_gen` (length: ngen): generator-bus incidence matrix `Cg`

"""
struct PolarNetworkState{VI,VT} <: AbstractNetworkBuffer
    vmag::VT
    vang::VT
    pinj::VT
    qinj::VT
    pg::VT
    qg::VT
    balance::VT
    dx::VT
    bus_gen::VI   # Generator-Bus incidence matrix
end

function PolarNetworkState{VT}(nbus::Int, ngen::Int, nstates::Int, bus_gen::VI) where {VI, VT}
    # Bus variables
    pbus = xzeros(VT, nbus)
    qbus = xzeros(VT, nbus)
    vmag = xzeros(VT, nbus)
    vang = xzeros(VT, nbus)
    # Generators variables
    pg = xzeros(VT, ngen)
    qg = xzeros(VT, ngen)
    # Buffers
    balance = xzeros(VT, nstates)
    dx = xzeros(VT, nstates)
    return PolarNetworkState{VI,VT}(vmag, vang, pbus, qbus, pg, qg, balance, dx, bus_gen)
end

setvalues!(buf::PolarNetworkState, ::PS.VoltageMagnitude, values) = copyto!(buf.vmag, values)
setvalues!(buf::PolarNetworkState, ::PS.VoltageAngle, values) = copyto!(buf.vang, values)
function setvalues!(buf::PolarNetworkState, ::PS.ActivePower, values)
    pgenbus = view(buf.pinj, buf.bus_gen)
    # Remove previous values
    pgenbus .-= buf.pg
    # Add new values
    copyto!(buf.pg, values)
    pgenbus .+= buf.pg
end
function setvalues!(buf::PolarNetworkState, ::PS.ReactivePower, values)
    qgenbus = view(buf.qinj, buf.bus_gen)
    # Remove previous values
    qgenbus .-= buf.qg
    # Add new values
    copyto!(buf.qg, values)
    qgenbus .+= buf.qg
end
function setvalues!(buf::PolarNetworkState, ::PS.ActiveLoad, values)
    fill!(buf.pinj, 0)
    # Pbus = Cg * Pg - Pd
    buf.pinj .-= values
    pgenbus = view(buf.pinj, buf.bus_gen)
    pgenbus .+= buf.pg
end
function setvalues!(buf::PolarNetworkState, ::PS.ReactiveLoad, values)
    fill!(buf.qinj, 0)
    # Qbus = Cg * Qg - Qd
    buf.qinj .-= values
    qgenbus = view(buf.qinj, buf.bus_gen)
    qgenbus .+= buf.qg
end

function Base.iszero(buf::PolarNetworkState)
    return iszero(buf.pinj) &&
        iszero(buf.qinj) &&
        iszero(buf.vmag) &&
        iszero(buf.vang) &&
        iszero(buf.pg) &&
        iszero(buf.qg) &&
        iszero(buf.balance) &&
        iszero(buf.dx)
end

"Store topology of the network on target device."
struct NetworkTopology{VTI, VTD}
    # Bus admittance matrix
    ybus_re::Spmat{VTI, VTD} # nb x nb
    ybus_im::Spmat{VTI, VTD} # nb x nb
    # Branches admittance matrix
    ## Real part
    yff_re::VTD # nl
    yft_re::VTD # nl
    ytf_re::VTD # nl
    ytt_re::VTD # nl
    ## Imag part
    yff_im::VTD # nl
    yft_im::VTD # nl
    ytf_im::VTD # nl
    ytt_im::VTD # nl
    # Correspondence
    f_buses::VTI # nl
    t_buses::VTI # nl
    sortperm::VTI # nnz
end

function NetworkTopology{VTI, VTD}(pf::PS.PowerNetwork) where {VTI, VTD}
    Y = pf.Ybus
    ybus_re, ybus_im = Spmat{VTI, VTD}(Y)
    lines = pf.lines
    yff_re = real.(lines.Yff) |> VTD
    yft_re = real.(lines.Yft) |> VTD
    ytf_re = real.(lines.Ytf) |> VTD
    ytt_re = real.(lines.Ytt) |> VTD

    yff_im = imag.(lines.Yff) |> VTD
    yft_im = imag.(lines.Yft) |> VTD
    ytf_im = imag.(lines.Ytf) |> VTD
    ytt_im = imag.(lines.Ytt) |> VTD

    f = lines.from_buses |> VTI
    t = lines.to_buses   |> VTI
    i, j, _ = findnz(Y)
    sp = sortperm(i) |> VTI

    return NetworkTopology(
        ybus_re, ybus_im,
        yff_re, yft_re, ytf_re, ytt_re,
        yff_im, yft_im, ytf_im, ytt_im,
        f, t, sp,
    )
end

get(net::NetworkTopology, ::PS.BusAdmittanceMatrix) = (net.ybus_re, net.ybus_im)

struct HessianStructure{IT} <: AbstractStructure where {IT}
    map::IT
end


text:
using Test

gpx_sample_file2=joinpath(@__DIR__, gpx_sample_file)
tcx_sample_file2=joinpath(@__DIR__, tcx_sample_file)
tcx_centrypark_run_file2=joinpath(@__DIR__, tcx_centrypark_run_file)
tcx_treadmill_run_file2=joinpath(@__DIR__, tcx_treadmill_run_file)

@testset "TESTSETS: Test read TCX file" begin
    @testset "CASE: Test file not XML" begin
        err, data =  TCX.parse_tcx_file(joinpath(@__DIR__, "runtests.jl"))
        @test err == 400
    end

    @testset "CASE: Test file not exists" begin
        err, data =  TCX.parse_tcx_file("non-existing-file.tcx")
        @test err == 404
    end

    @testset "CASE: Test file exist but not a TCX file" begin
        err, _ =  TCX.parse_tcx_file(gpx_sample_file2)
        @test err == 400
    end

    @testset "CASE: Test a TCX records run." begin
        err, data =  TCX.parse_tcx_file(tcx_sample_file2)
        @test (err == 200) & (getActivityType(data) == "Running")
    end

    @testset "CASE: Test a TCX contains track points." begin
        err, data =  TCX.parse_tcx_file(tcx_sample_file2)
        @test (err == 200) & (size(getDataFrame(data), 1) == 3427)
    end

    @testset "CASE: Test a TCX treadmill run, no logitude/latitude." begin
        err, data =  TCX.parse_tcx_file(tcx_treadmill_run_file2)
        @test (err == 200) & (getDistance2(data) == 0)
    end

    @testset "CASE: Test a TCX static running distance." begin
        err, data =  TCX.parse_tcx_file(tcx_centrypark_run_file2)
        @test (err == 200) & (getDistance(data) == 27020)
    end

    @testset "CASE: Test a TCX running distance from Geodesy." begin
        err, data =  TCX.parse_tcx_file(tcx_centrypark_run_file2)
        @test (err == 200) & (getDistance2(data) > 0)
    end

    @testset "CASE: Test a TCX static average pace." begin
        err, data =  TCX.parse_tcx_file(tcx_centrypark_run_file2)
        @test (err == 200) & (getAveragePace(data) > 0 )
    end

    @testset "CASE: Test to get TCX duration." begin
        err, data =  TCX.parse_tcx_file(tcx_centrypark_run_file2)
        @test (err == 200) & (getDuration(data) > 0 )
    end
end


text:
#Ax_mul_Bxt

export Ax_mul_Bxt

"""
`Ax_mul_Bxt(A::AbstractOperator,B::AbstractOperator)`

Create an operator `P` such that:

`P == (Ax)*(Bx)'`

# Example: Matrix multiplication

```julia
julia> A,B = randn(4,4),randn(4,4);

julia> P = Ax_mul_Bxt(MatrixOp(A),MatrixOp(B))
▒*▒  ℝ^4 -> ℝ^(4, 4)

julia> x = randn(4);

julia> P*x == (A*x)*(B*x)'
true

```
"""
struct Ax_mul_Bxt{
                   L1 <: AbstractOperator,
                   L2 <: AbstractOperator,
                   C <: AbstractArray,
                   D <: AbstractArray,
                  } <: NonLinearOperator
  A::L1
  B::L2
  bufA::C
  bufB::C
  bufC::C
  bufD::D
  function Ax_mul_Bxt(A::L1, B::L2, bufA::C, bufB::C, bufC::C, bufD::D) where {L1,L2,C,D}
    if ndims(A,1) == 1
      if size(A) != size(B)
        throw(DimensionMismatch("Cannot compose operators"))
      end
    elseif ndims(A,1) == 2 && ndims(B,1) == 2 && size(A,2) == size(B,2)
      if size(A,1)[2] != size(B,1)[2]
        throw(DimensionMismatch("Cannot compose operators"))
      end
    else
      throw(DimensionMismatch("Cannot compose operators"))
    end
    new{L1,L2,C,D}(A,B,bufA,bufB,bufC,bufD)
  end
end

struct Ax_mul_BxtJac{
                      L1 <: AbstractOperator,
                      L2 <: AbstractOperator,
                      C <: AbstractArray,
                      D <: AbstractArray,
                     } <: LinearOperator
  A::L1
  B::L2
  bufA::C
  bufB::C
  bufC::C
  bufD::D
end

# Constructors
function Ax_mul_Bxt(A::AbstractOperator,B::AbstractOperator)
  s,t = size(A,1), codomainType(A)
  bufA = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)
  bufC = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)
  s,t = size(B,1), codomainType(B)
  bufB = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)
  s,t = size(A,2), domainType(A)
  bufD = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)
  Ax_mul_Bxt(A,B,bufA,bufB,bufC,bufD)
end

# Jacobian
function Jacobian(P::Ax_mul_Bxt{L1,L2,C,D}, x::AbstractArray) where {L1,L2,C,D}
  JA, JB = Jacobian(P.A, x), Jacobian(P.B, x)
  Ax_mul_BxtJac{typeof(JA),typeof(JB),C,D}(JA,JB,P.bufA,P.bufB,P.bufC,P.bufD)
end

# Mappings
function mul!(y, P::Ax_mul_Bxt{L1,L2,C,D}, b) where {L1,L2,C,D}
  mul!(P.bufA,P.A,b)
  mul!(P.bufB,P.B,b)
  mul!(y,P.bufA, P.bufB')
end

function mul!(y, J::AdjointOperator{Ax_mul_BxtJac{L1,L2,C,D}}, b) where {L1,L2,C,D}
  #y .= J.A.A'*(b*(J.A.bufB)) + J.A.B'*(b'*(J.A.bufA))
  mul!(J.A.bufC, b, J.A.bufB)
  mul!(y, J.A.A', J.A.bufC)
  mul!(J.A.bufB, b', J.A.bufA)
  mul!(J.A.bufD, J.A.B', J.A.bufB)
  y .+= J.A.bufD
  return y
end

size(P::Union{Ax_mul_Bxt,Ax_mul_BxtJac}) = ((size(P.A,1)[1],size(P.B,1)[1]),size(P.A,2))

fun_name(L::Union{Ax_mul_Bxt,Ax_mul_BxtJac}) = fun_name(L.A)*"*"*fun_name(L.B)

domainType(L::Union{Ax_mul_Bxt,Ax_mul_BxtJac})   = domainType(L.A)
codomainType(L::Union{Ax_mul_Bxt,Ax_mul_BxtJac}) = codomainType(L.A)

# utils
function permute(P::Ax_mul_Bxt{L1,L2,C,D},
                 p::AbstractVector{Int}) where {L1,L2,C,D <:ArrayPartition}
  Ax_mul_Bxt(permute(P.A,p),permute(P.B,p),P.bufA,P.bufB,P.bufC,ArrayPartition(P.bufD.x[p]) )
end

remove_displacement(P::Ax_mul_Bxt) =
Ax_mul_Bxt(remove_displacement(P.A), remove_displacement(P.B), P.bufA, P.bufB, P.bufC, P.bufD)


text:
# Autogenerated wrapper script for Osi_jll for armv7l-linux-gnueabihf-cxx03
export libOsi, libOsiCommonTests

using CoinUtils_jll
using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Osi")
JLLWrappers.@declare_library_product(libOsi, "libOsi.so.1")
JLLWrappers.@declare_library_product(libOsiCommonTests, "libOsiCommonTests.so.1")
function __init__()
    JLLWrappers.@generate_init_header(CoinUtils_jll, OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libOsi,
        "lib/libOsi.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libOsiCommonTests,
        "lib/libOsiCommonTests.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
header("LABELS")
@testset "labels" begin
    # Exersize the :labels of a PnmlDict

    d = PNML.PnmlDict(:labels => PNML.PnmlLabel[])
    reg = PNML.IDRegistry()
    for i in 1:4
        x = i<3 ? 1 : 2
        node = root(EzXML.parsexml("<test$x> $i </test$x>"))
        n = PNML.add_label!(d, node, PnmlCore(); reg)
        @test_call PNML.add_label!(d, node, PnmlCore(); reg)
        @test length(PNML.labels(d)) == i
        @test d[:labels] == n # Returned value is the vector
        @test PNML.labels(d) == n # Returned value is the vector
    end
    printnode(d)
    @test_call PNML.labels(d)
    # labels(d) infers as Any
    @test PNML.labels(d) isa Vector{PNML.PnmlLabel}
    @test length(d[:labels]) == 4
    @test length(labels(d)) == 4
    foreach(PNML.labels(d)) do l
        @test_call tag(l)
        @test tag(l) === :test1 || tag(l) === :test2
        @test xmlnode(l) isa Maybe{EzXML.Node}
    end
    @show typeof(d), collect(keys(d))
    @show typeof(labels(d))

    @test_call PNML.has_label(d, :test1)
    @test_call PNML.get_label(d, :test1)
    @test_call PNML.get_labels(d, :test1)

    @test PNML.has_label(d, :test1)
    @test !PNML.has_label(d, :bumble)

    @show v = PNML.get_label(d, :test2)
    @test v.dict[:content] == "3"

    @testset "label $labeltag" for labeltag in [:test1, :test2]
        @show v = PNML.get_labels(d, labeltag)
        @test length(v) == 2
        for l in v
            @test tag(l) === labeltag
        end
    end
end


text:
# Efficient implementation a' * Q * b that avoid unnecessary type promotion as
# well as unnecessary allocation
function quad_form(a::AbstractVector{<:Real},
                   Q::SymMatrix{<:Real},
                   b::AbstractVector{<:Real})
    n = length(a)
    @assert n == LinearAlgebra.checksquare(Q)
    @assert n == length(b)
    out = zero(typeof(zero(eltype(a)) * zero(eltype(Q)) * zero(eltype(b))))
    k = 0
    for j in 1:n
        α = a[j]
        β = b[j]
        for i in 1:(j-1)
            k += 1
            out += a[i] * Q.Q[k] * β
            out += α * Q.Q[k] * b[i]
        end
        k += 1
        out += a[j] * Q.Q[k] * b[j]
    end
    return out
end
function quad_form(a::AbstractVector,
                   Q::AbstractMatrix,
                   b::AbstractVector)
    (n, m) = size(Q)
    n == length(a) || throw(DimensionMismatch())
    m == length(b) || throw(DimensionMismatch())
    U = MA.promote_operation(*, eltype(a), eltype(Q), eltype(b))
    out = zero(MA.promote_operation(+, U, U))
    for j in 1:m
        for i in 1:n
            out = MA.add_mul!!(out, a[i], Q[i, j], b[j])
        end
    end
    return out
end
# Same as quad_form(a, Q, a)
function quad_form(Q::Symmetric,
                   a::AbstractVector)
    return quad_form(a, Q, a)
end
function quad_form(Q::Symmetric{JuMP.VariableRef}, a::AbstractVector{<:Real})
    return quad_form(a, Q, a)
end

# We have x' Q x and we want y' Q y where y is obtained by substituting
# vars for A * new_vars in x. We want to compute the matrix M such that
# y = M z where z is the vector of monomials of degree d of new_vars
# Then we will have y' Q y = z' M' Q M z
struct GramTransformation{T, MT <: MP.AbstractMonomial,
                          MVT <: AbstractVector{MT}}
    M::Matrix{T}
    monos::MVT
end

function apply_transformation(p::SumOfSquares.GramMatrix,
                              t::GramTransformation) where T
    new_n = length(t.monos)
    new_Q = [quad_form(t.M[:, i], p.Q, t.M[:, j]) for j in 1:new_n for i in 1:j]
    return GramMatrix(SymMatrix(new_Q, new_n), t.monos)
end

function transformation(old_monos, A::AbstractMatrix, new_vars, d)
    new_monos = monomials(new_vars, d)
    new_n = length(new_monos)
    M = zeros(eltype(A), length(old_monos), new_n)
    mapped_vars = A * new_vars
    # Cache the result of mapped_vars[i]^n
    powers = [Union{Nothing, eltype(mapped_vars)}[mapped_vars[i]]
              for i in eachindex(mapped_vars)]
    # Compute mapped_vars[i]^n by "Power by Squaring" and cache it in `powers`
    function _power(i, n)
        @assert n > 0
        while n > length(powers[i])
            push!(powers[i], nothing)
        end
        if powers[i][n] === nothing
            if isodd(n)
                powers[i][n] = mapped_vars[i] * _power(i, n - 1)
            else
                p_2 = _power(i, div(n, 2))
                powers[i][n] = p_2 * p_2
            end
        end
        return powers[i][n]::eltype(mapped_vars)
    end
    function _map(mono::MP.AbstractMonomial)
        exps = MP.exponents(mono)
        length(exps) == length(mapped_vars) || throw(ArgumentError("A monomial have less variables than `new_vars`"))
        cur = one(eltype(mapped_vars))
        for i in eachindex(exps)
            if exps[i] > 0
                cur *= _power(i, exps[i])
            end
        end
        return cur
    end
    for i in eachindex(old_monos)
        y = _map(old_monos[i])
        j = 1
        for term in MP.terms(y)
            mono = MP.monomial(term)
            while j <= length(new_monos) && mono < new_monos[j]
                j += 1
            end
            M[i, j] = MP.coefficient(term)
        end
    end
    return GramTransformation(M, new_monos)
end

# computes p ∘ A or more precisely p(variables(p) => A * new_vars)
function apply_matrix(p::SumOfSquares.GramMatrix{T, <:MonomialBasis},
                      A::AbstractMatrix, new_vars, d) where T
    return apply_transformation(p, transformation(p.basis.monomials, A, new_vars, d))
end

# computes A # μ or more precisely p(variables(p) => A * new_vars)

function psd_constraint(Q::Symmetric)
    n = LinearAlgebra.checksquare(Q)
    q = [Q[i, j] for j in 1:n for i in 1:j]
    # For n == 0, it will create no constraint, for n == 1, it will simply
    # be a Nonnegatives constraint and for n == 2 it will be a rotated SOC.
    set = SumOfSquares.matrix_cone(MOI.PositiveSemidefiniteConeTriangle, n)
    return PolyJuMP.bridgeable(JuMP.build_constraint(error, q, set),
                               JuMP.moi_function_type(typeof(q)), typeof(set))
end

function psd_constraint(model, Q::Symmetric)
    return JuMP.add_constraint(model, psd_constraint(Q))
end


text:
let
    G1 = MetaGraph()
    add_vertex!(G1,Dict(:x=>0,:y=>0))
    add_vertex!(G1,Dict(:x=>1,:y=>0))
    add_vertex!(G1,Dict(:x=>2,:y=>0))
    add_vertex!(G1,Dict(:x=>3,:y=>0))
    add_vertex!(G1,Dict(:x=>4,:y=>0))
    # add_vertex!(G1,Dict(:x=>1,:y=>1))
    # add_vertex!(G1,Dict(:x=>3,:y=>1))
    add_edge!(G1,1,2)
    add_edge!(G1,2,3)
    add_edge!(G1,3,4)
    add_edge!(G1,4,5)
    # add_edge!(G1,2,6)
    # add_edge!(G1,4,7)
    start_vtxs = [1, 5]
    goal_vtxs = [2, 3]

    T = 2 # time horizon to extend graph
    # Gx, model, x = construct_ILP_flow_model(G1,T,start_vtxs,goal_vtxs,milp_optimizer())
    # optimize!(model)
    Gx = CRCBS.FlowProblems.extend_G_through_time(G1,T)
    @test nv(Gx) == nv(G1)*(T+1) + ne(G1)*4*T
    x = [
        1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0,
        0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
        1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0
        ]
    paths = CRCBS.FlowProblems.extract_solution_from_flow_matrix(Gx,x,start_vtxs,goal_vtxs)
    true_paths = CRCBS.FlowProblems.extract_paths_from_flow_solution(G1,paths)
    for (i,path) in enumerate(true_paths)
        @test path[1] == start_vtxs[i]
        @test path[end] in goal_vtxs
    end
end


text:
# Autogenerated wrapper script for Julia_jll for aarch64-linux-gnu-libgfortran4-cxx11
export julia, libjulia

using Zlib_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `julia`
const julia_splitpath = ["bin", "julia"]

# This will be filled out by __init__() for all products, as it must be done at runtime
julia_path = ""

# julia-specific global declaration
function julia(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(julia_path)
    end
end

# Relative path to `libjulia`
const libjulia_splitpath = ["lib", "libjulia.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libjulia_path = ""

# libjulia-specific global declaration
# This will be filled out by __init__()
libjulia_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libjulia = "libjulia.so.1"

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"Julia")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (Zlib_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (Zlib_jll.LIBPATH_list,))

    global julia_path = normpath(joinpath(artifact_dir, julia_splitpath...))

    push!(PATH_list, dirname(julia_path))
    global libjulia_path = normpath(joinpath(artifact_dir, libjulia_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libjulia_handle = dlopen(libjulia_path)
    push!(LIBPATH_list, dirname(libjulia_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()


text:
module ExampleMethods

export rk4_table, extrapolated_euler_table,
    rkck5_table, dopri5_table, rkf8_table

# This is the classical Runge-Kutta method, for which no introduction should
# be needed. (If you don't know what this is, you probably need to do some more
# background reading before using this package!)
rk4_table(::Type{T}) where {T <: Real} = T[
    inv(T(2)), T(0), inv(T(2)), T(0), T(0), T(1),
    inv(T(6)), inv(T(3)), inv(T(3)), inv(T(6))]

# This is the Runge-Kutta method obtained by applying Richardson extrapolation
# to n independent executions of Euler's method using step sizes h, h/2, h/3,
# ..., h/n, yielding a method of order n. The methods obtained in this fashion
# are not practically useful, but provide a simple proof that Runge-Kutta
# methods exist of any order using a quadratic number of stages.
function extrapolated_euler_table(order::Int)::Vector{Rational{BigInt}}
    result = Vector{Rational{BigInt}}[]
    skip = 0
    for i = 2 : order
        for j = 1 : i - 1
            push!(result, vcat(
                [Rational{BigInt}(1, i)],
                zeros(Rational{BigInt}, skip),
                [Rational{BigInt}(1, i) for _ = 1 : j - 1]))
        end
        skip += i - 1
    end
    num_stages = div(order * (order - 1), 2) + 1
    mat = zeros(Rational{BigInt}, order, num_stages)
    skip = 0
    for i = 1 : order
        mat[i, 1] = Rational{BigInt}(1, i)
        for j = skip + 2 : skip + i
            mat[i, j] = Rational{BigInt}(1, i)
        end
        skip += i - 1
    end
    lhs = [Rational{BigInt}(1, j)^i for i = 0 : order - 1, j = 1 : order]
    rhs = [Rational{BigInt}(i == 1, 1) for i = 1 : order]
    push!(result, transpose(mat) * (lhs \ rhs))
    vcat(result...)
end

function extrapolated_euler_table(::Type{T},
        order::Int)::Vector{T} where {T <: Real}
    T.(extrapolated_euler_table(order))
end

# This is the 5th-order method presented on p. 206 of the following paper.
# Interestingly, it is presented not as a standalone method or in an embedded
# pair, but as the highest-order component of an embedded quintuple of orders
# 5(4,3,2,1). Cash and Karp themselves named it "RKFNC"; I'm not sure why.
# Numerical Recipes calls this method "RKCK."
# Cash and Karp 1990, "A Variable Order Runge-Kutta Method for Initial Value
#                      Problems with Rapidly Varying Right-Hand Sides"
rkck5_table(::Type{T}) where {T <: Real} = T[
    inv(T(5)), T(3)/T(40), T(9)/T(40), T(3)/T(10), T(-9)/T(10), T(6)/T(5),
    T(-11)/T(54), T(5)/T(2), T(-70)/T(27), T(35)/T(27), T(1631)/T(55296),
    T(175)/T(512), T(575)/T(13824), T(44275)/T(110592), T(253)/T(4096),
    T(37)/T(378), T(0), T(250)/T(621), T(125)/T(594), T(0), T(512)/T(1771)]

# This is the higher-order component of the 5(4) embedded pair presented on
# p. 23 of the following paper. Dormand and Prince call this method "RK5(4)7M,"
# but it has become commonly known as "DOPRI5" following a popular Fortran
# implementation.
# Dormand and Prince 1980, "A family of embedded Runge-Kutta formulae"
dopri5_table(::Type{T}) where {T <: Real} = T[
    inv(T(5)), T(3)/T(40), T(9)/T(40), T(44)/T(45), T(-56)/T(15), T(32)/T(9),
    T(19372)/T(6561), T(-25360)/T(2187), T(64448)/T(6561), T(-212)/T(729),
    T(9017)/T(3168), T(-355)/T(33), T(46732)/T(5247), T(49)/T(176),
    T(-5103)/T(18656), T(35)/T(384), T(0), T(500)/T(1113), T(125)/T(192),
    T(-2187)/T(6784), T(11)/T(84), T(35)/T(384), T(0), T(500)/T(1113),
    T(125)/T(192), T(-2187)/T(6784), T(11)/T(84), T(0)]

# This is the 8th-order method presented on p. 65 of Fehlberg's 1968 NASA
# paper as the higher-order component of the embedded pair RK7(8).
# Fehlberg 1968, "Classical Fifth-, Sixth-, Seventh-, and Eighth-Order
#                 Runge-Kutta Formulas with Stepsize Control"
rkf8_table(::Type{T}) where {T <: Real} = T[
    T(2)/T(27), inv(T(36)), inv(T(12)), inv(T(24)), T(0), inv(T(8)),
    T(5)/T(12), T(0), T(-25)/T(16), T(25)/T(16), inv(T(20)), T(0), T(0),
    inv(T(4)), inv(T(5)), T(-25)/T(108), T(0), T(0), T(125)/T(108),
    T(-65)/T(27), T(125)/T(54), T(31)/T(300), T(0), T(0), T(0), T(61)/T(225),
    T(-2)/T(9), T(13)/T(900), T(2)/T(1), T(0), T(0), T(-53)/T(6), T(704)/T(45),
    T(-107)/T(9), T(67)/T(90), T(3)/T(1), T(-91)/T(108), T(0), T(0),
    T(23)/T(108), T(-976)/T(135), T(311)/T(54), T(-19)/T(60), T(17)/T(6),
    T(-1)/T(12), T(2383)/T(4100), T(0), T(0), T(-341)/T(164), T(4496)/T(1025),
    T(-301)/T(82), T(2133)/T(4100), T(45)/T(82), T(45)/T(164), T(18)/T(41),
    T(3)/T(205), T(0), T(0), T(0), T(0), T(-6)/T(41), T(-3)/T(205),
    T(-3)/T(41), T(3)/T(41), T(6)/T(41), T(0), T(-1777)/T(4100), T(0), T(0),
    T(-341)/T(164), T(4496)/T(1025), T(-289)/T(82), T(2193)/T(4100),
    T(51)/T(82), T(33)/T(164), T(12)/T(41), T(0), T(1), T(0), T(0), T(0), T(0),
    T(0), T(34)/T(105), T(9)/T(35), T(9)/T(35), T(9)/T(280), T(9)/T(280), T(0),
    T(41)/T(840), T(41)/T(840)]

end # module ExampleMethods


text:
# Julia wrapper for header: p6est.h
# Automatically generated using Clang.jl

function p6est_connectivity_new(conn4, top_vertices, height)
    ccall((:p6est_connectivity_new, p4est_lib), Ptr{p6est_connectivity_t}, (Ptr{p4est_connectivity_t}, Ptr{Cdouble}, Ptr{Cdouble}), conn4, top_vertices, height)
end

function p6est_connectivity_destroy(conn)
    ccall((:p6est_connectivity_destroy, p4est_lib), Cvoid, (Ptr{p6est_connectivity_t},), conn)
end

function p6est_tree_get_vertices(conn, which_tree, vertices)
    ccall((:p6est_tree_get_vertices, p4est_lib), Cvoid, (Ptr{p6est_connectivity_t}, p4est_topidx_t, Ptr{Cdouble}), conn, which_tree, vertices)
end

function p6est_qcoord_to_vertex(connectivity, treeid, x, y, z, vxyz)
    ccall((:p6est_qcoord_to_vertex, p4est_lib), Cvoid, (Ptr{p6est_connectivity_t}, p4est_topidx_t, p4est_qcoord_t, p4est_qcoord_t, p4est_qcoord_t, Ptr{Cdouble}), connectivity, treeid, x, y, z, vxyz)
end

function p6est_new(mpicomm, connectivity, data_size, init_fn, user_pointer)
    ccall((:p6est_new, p4est_lib), Ptr{p6est_t}, (MPI_Comm, Ptr{p6est_connectivity_t}, Csize_t, p6est_init_t, Ptr{Cvoid}), mpicomm, connectivity, data_size, init_fn, user_pointer)
end

function p6est_new_from_p4est(p4est, top_vertices, height, min_zlevel, data_size, init_fn, user_pointer)
    ccall((:p6est_new_from_p4est, p4est_lib), Ptr{p6est_t}, (Ptr{p4est_t}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Csize_t, p6est_init_t, Ptr{Cvoid}), p4est, top_vertices, height, min_zlevel, data_size, init_fn, user_pointer)
end

function p6est_destroy(p6est)
    ccall((:p6est_destroy, p4est_lib), Cvoid, (Ptr{p6est_t},), p6est)
end

function p6est_copy(input, copy_data)
    ccall((:p6est_copy, p4est_lib), Ptr{p6est_t}, (Ptr{p6est_t}, Cint), input, copy_data)
end

function p6est_reset_data(p6est, data_size, init_fn, user_pointer)
    ccall((:p6est_reset_data, p4est_lib), Cvoid, (Ptr{p6est_t}, Csize_t, p6est_init_t, Ptr{Cvoid}), p6est, data_size, init_fn, user_pointer)
end

function p6est_refine_columns(p6est, refine_recursive, refine_fn, init_fn)
    ccall((:p6est_refine_columns, p4est_lib), Cvoid, (Ptr{p6est_t}, Cint, p6est_refine_column_t, p6est_init_t), p6est, refine_recursive, refine_fn, init_fn)
end

function p6est_refine_layers(p6est, refine_recursive, refine_fn, init_fn)
    ccall((:p6est_refine_layers, p4est_lib), Cvoid, (Ptr{p6est_t}, Cint, p6est_refine_layer_t, p6est_init_t), p6est, refine_recursive, refine_fn, init_fn)
end

function p6est_coarsen_columns(p6est, coarsen_recursive, coarsen_fn, init_fn)
    ccall((:p6est_coarsen_columns, p4est_lib), Cvoid, (Ptr{p6est_t}, Cint, p6est_coarsen_column_t, p6est_init_t), p6est, coarsen_recursive, coarsen_fn, init_fn)
end

function p6est_coarsen_layers(p6est, coarsen_recursive, coarsen_fn, init_fn)
    ccall((:p6est_coarsen_layers, p4est_lib), Cvoid, (Ptr{p6est_t}, Cint, p6est_coarsen_layer_t, p6est_init_t), p6est, coarsen_recursive, coarsen_fn, init_fn)
end

function p6est_balance(p6est, btype, init_fn)
    ccall((:p6est_balance, p4est_lib), Cvoid, (Ptr{p6est_t}, p8est_connect_type_t, p6est_init_t), p6est, btype, init_fn)
end

function p6est_partition(p6est, weight_fn)
    ccall((:p6est_partition, p4est_lib), p4est_gloidx_t, (Ptr{p6est_t}, p6est_weight_t), p6est, weight_fn)
end

function p6est_partition_correct(p6est, num_layers_in_proc)
    ccall((:p6est_partition_correct, p4est_lib), Cvoid, (Ptr{p6est_t}, Ptr{p4est_locidx_t}), p6est, num_layers_in_proc)
end

function p6est_partition_to_p4est_partition(p6est, num_layers_in_proc, num_columns_in_proc)
    ccall((:p6est_partition_to_p4est_partition, p4est_lib), Cvoid, (Ptr{p6est_t}, Ptr{p4est_locidx_t}, Ptr{p4est_locidx_t}), p6est, num_layers_in_proc, num_columns_in_proc)
end

function p4est_partition_to_p6est_partition(p6est, num_columns_in_proc, num_layers_in_proc)
    ccall((:p4est_partition_to_p6est_partition, p4est_lib), Cvoid, (Ptr{p6est_t}, Ptr{p4est_locidx_t}, Ptr{p4est_locidx_t}), p6est, num_columns_in_proc, num_layers_in_proc)
end

function p6est_partition_for_coarsening(p6est, num_layers_in_proc)
    ccall((:p6est_partition_for_coarsening, p4est_lib), p4est_gloidx_t, (Ptr{p6est_t}, Ptr{p4est_locidx_t}), p6est, num_layers_in_proc)
end

function p6est_partition_given(p6est, num_layers_in_proc)
    ccall((:p6est_partition_given, p4est_lib), p4est_gloidx_t, (Ptr{p6est_t}, Ptr{p4est_locidx_t}), p6est, num_layers_in_proc)
end

function p6est_checksum(p6est)
    ccall((:p6est_checksum, p4est_lib), UInt32, (Ptr{p6est_t},), p6est)
end

function p6est_save(filename, p6est, save_data)
    ccall((:p6est_save, p4est_lib), Cvoid, (Cstring, Ptr{p6est_t}, Cint), filename, p6est, save_data)
end

function p6est_load(filename, mpicomm, data_size, load_data, user_pointer, connectivity)
    ccall((:p6est_load, p4est_lib), Ptr{p6est_t}, (Cstring, MPI_Comm, Csize_t, Cint, Ptr{Cvoid}, Ptr{Ptr{p6est_connectivity_t}}), filename, mpicomm, data_size, load_data, user_pointer, connectivity)
end

function p6est_compress_columns(p6est)
    ccall((:p6est_compress_columns, p4est_lib), Cvoid, (Ptr{p6est_t},), p6est)
end

function p6est_update_offsets(p6est)
    ccall((:p6est_update_offsets, p4est_lib), Cvoid, (Ptr{p6est_t},), p6est)
end
# Julia wrapper for header: p6est_extended.h
# Automatically generated using Clang.jl

function p6est_new_ext(mpicomm, connectivity, min_quadrants, min_level, min_zlevel, num_zroot, fill_uniform, data_size, init_fn, user_pointer)
    ccall((:p6est_new_ext, p4est_lib), Ptr{p6est_t}, (MPI_Comm, Ptr{p6est_connectivity_t}, p4est_locidx_t, Cint, Cint, Cint, Cint, Csize_t, p6est_init_t, Ptr{Cvoid}), mpicomm, connectivity, min_quadrants, min_level, min_zlevel, num_zroot, fill_uniform, data_size, init_fn, user_pointer)
end

function p6est_copy_ext(input, copy_data, duplicate_mpicomm)
    ccall((:p6est_copy_ext, p4est_lib), Ptr{p6est_t}, (Ptr{p6est_t}, Cint, Cint), input, copy_data, duplicate_mpicomm)
end

function p6est_save_ext(filename, p6est, save_data, save_partition)
    ccall((:p6est_save_ext, p4est_lib), Cvoid, (Cstring, Ptr{p6est_t}, Cint, Cint), filename, p6est, save_data, save_partition)
end

function p6est_load_ext(filename, mpicomm, data_size, load_data, autopartition, broadcasthead, user_pointer, connectivity)
    ccall((:p6est_load_ext, p4est_lib), Ptr{p6est_t}, (Cstring, MPI_Comm, Csize_t, Cint, Cint, Cint, Ptr{Cvoid}, Ptr{Ptr{p6est_connectivity_t}}), filename, mpicomm, data_size, load_data, autopartition, broadcasthead, user_pointer, connectivity)
end

function p6est_refine_columns_ext(p6est, refine_recursive, maxlevel, refine_fn, init_fn, replace_fn)
    ccall((:p6est_refine_columns_ext, p4est_lib), Cvoid, (Ptr{p6est_t}, Cint, Cint, p6est_refine_column_t, p6est_init_t, p6est_replace_t), p6est, refine_recursive, maxlevel, refine_fn, init_fn, replace_fn)
end

function p6est_refine_layers_ext(p6est, refine_recursive, maxlevel, refine_fn, init_fn, replace_fn)
    ccall((:p6est_refine_layers_ext, p4est_lib), Cvoid, (Ptr{p6est_t}, Cint, Cint, p6est_refine_layer_t, p6est_init_t, p6est_replace_t), p6est, refine_recursive, maxlevel, refine_fn, init_fn, replace_fn)
end

function p6est_coarsen_columns_ext(p6est, coarsen_recursive, callback_orphans, coarsen_fn, init_fn, replace_fn)
    ccall((:p6est_coarsen_columns_ext, p4est_lib), Cvoid, (Ptr{p6est_t}, Cint, Cint, p6est_coarsen_column_t, p6est_init_t, p6est_replace_t), p6est, coarsen_recursive, callback_orphans, coarsen_fn, init_fn, replace_fn)
end

function p6est_coarsen_layers_ext(p6est, coarsen_recursive, callback_orphans, coarsen_fn, init_fn, replace_fn)
    ccall((:p6est_coarsen_layers_ext, p4est_lib), Cvoid, (Ptr{p6est_t}, Cint, Cint, p6est_coarsen_layer_t, p6est_init_t, p6est_replace_t), p6est, coarsen_recursive, callback_orphans, coarsen_fn, init_fn, replace_fn)
end

function p6est_partition_ext(p6est, partition_for_coarsening, weight_fn)
    ccall((:p6est_partition_ext, p4est_lib), p4est_gloidx_t, (Ptr{p6est_t}, Cint, p6est_weight_t), p6est, partition_for_coarsening, weight_fn)
end

function p6est_balance_ext(p6est, btype, max_diff, min_diff, init_fn, replace_fn)
    ccall((:p6est_balance_ext, p4est_lib), Cvoid, (Ptr{p6est_t}, p8est_connect_type_t, Cint, Cint, p6est_init_t, p6est_replace_t), p6est, btype, max_diff, min_diff, init_fn, replace_fn)
end
# Julia wrapper for header: p6est_ghost.h
# Automatically generated using Clang.jl

function p6est_ghost_memory_used()
    ccall((:p6est_ghost_memory_used, p4est_lib), Cint, ())
end

function p6est_ghost_new(p4est, btype)
    ccall((:p6est_ghost_new, p4est_lib), Ptr{p6est_ghost_t}, (Ptr{p6est_t}, p4est_connect_type_t), p4est, btype)
end

function p6est_ghost_expand(p6est, ghost)
    ccall((:p6est_ghost_expand, p4est_lib), Cvoid, (Ptr{p6est_t}, Ptr{p6est_ghost_t}), p6est, ghost)
end

function p6est_ghost_destroy(ghost)
    ccall((:p6est_ghost_destroy, p4est_lib), Cvoid, (Ptr{p6est_ghost_t},), ghost)
end

function p6est_ghost_bsearch(ghost, which_proc, which_tree, column, layer)
    ccall((:p6est_ghost_bsearch, p4est_lib), ssize_t, (Ptr{p6est_ghost_t}, Cint, p4est_topidx_t, Ptr{p4est_quadrant_t}, Ptr{p2est_quadrant_t}), ghost, which_proc, which_tree, column, layer)
end

function p6est_ghost_contains(ghost, which_proc, which_tree, column, layer)
    ccall((:p6est_ghost_contains, p4est_lib), ssize_t, (Ptr{p6est_ghost_t}, Cint, p4est_topidx_t, Ptr{p4est_quadrant_t}, Ptr{p2est_quadrant_t}), ghost, which_proc, which_tree, column, layer)
end

function p6est_layer_exists(p6est, ghost, treeid, column, layer, exists_arr, rproc_arr, rquad_arr)
    ccall((:p6est_layer_exists, p4est_lib), Cint, (Ptr{p6est_t}, Ptr{p6est_ghost_t}, p4est_topidx_t, Ptr{p4est_quadrant_t}, Ptr{p2est_quadrant_t}, Ptr{sc_array_t}, Ptr{sc_array_t}, Ptr{sc_array_t}), p6est, ghost, treeid, column, layer, exists_arr, rproc_arr, rquad_arr)
end

function p6est_is_balanced(p6est, btype)
    ccall((:p6est_is_balanced, p4est_lib), Cint, (Ptr{p6est_t}, p8est_connect_type_t), p6est, btype)
end

function p6est_ghost_checksum(p6est, ghost)
    ccall((:p6est_ghost_checksum, p4est_lib), UInt32, (Ptr{p6est_t}, Ptr{p6est_ghost_t}), p6est, ghost)
end


text:
using Game2048, Flux
using BenchmarkTools
using Flux: relu, softmax
using StatsBase, Statistics

function mirrors_rotations(board)
    tboard = transpose(board) |> collect

    return board,  rotl90(board),  rot180(board),  rotr90(board),
           tboard, rotl90(tboard), rot180(tboard), rotr90(tboard)
end

# simulate the game current that

# do a policy gradient approach
# for each step update the policy
model = Chain(
    normalised_board -> reshape(normalised_board, 4, 4, 1, 1),
    Conv((3,3), 1=>256, relu),
    Conv((3,3), 256=>128, relu),
    x->reshape(x, 128),
    Dense(128, 4),
    softmax
)

loss(board, means_scores) = crossentropy(model(board), softmax(means_scores))

params = Flux.params(model);

opt = ADAM()

function simulate_game_for_mean_score(board, dir)::Tuple{Float32, Bool}
    pts, possible = move!(deepcopy(board), dir)

    if possible
        return mean(simulate_games(1000, board, pts)), possible
    else
        return 0.0, possible
    end
end

function return_best_dir(board)
    ms_pos = simulate_game_for_mean_score.(Ref(board), DIRS)
    mean_scores = [ms for (ms, _) in ms_pos]
    possible = [p for (_, p) in ms_pos]
    (mean_scores .- mean(mean_scores))./std(mean_scores), possible
end

function play_game_with_policy_update!(loss, params, model, opt)
    board = initboard()
    most_sure_move_pct = 0.25

    while true
        mean_scores, possible = return_best_dir(board)

        if all(!, possible)
            return board
        end

        board_float32 = board ./ maximum(board)
        # training_data = [(b, mean_scores) for b in mirrors_rotations(board_float32)]

        training_data = [(board_float32, mean_scores)]

        Flux.train!(loss, params, training_data, opt)
        prob = model(board_float32)
        if maximum(prob) > most_sure_move_pct
            max_prob, i = findmax(prob)
            most_sure_move_pct = max_prob
            println("it is very sure - '$(DIRS[i])' $most_sure_move_pct")
            display(board)
        end
        # sample a direction from policy
        dir = wsample(DIRS[possible], prob[possible])
        # println(dir)
        move!(board, dir)
        gen_new_tile!(board)
        # display(board)
    end
end

function play_game_with_policy(model)
    board = initboard()

    while true
        possible = [p for (_, p) in move.(Ref(board), DIRS)]
        if all(!, possible)
            return board
        end

        board_float32 = board ./ maximum(board)
        prob = model(board_float32)
        # println(maximum(prob))

        # sample a direction from policy
        dir = wsample(DIRS[possible], prob[possible])
        # println(dir)
        move!(board, dir)
        gen_new_tile!(board)
        # display(board)
    end
end

function train_many_times!(res, times, loss, params, model, opt)
    for i in 1:times
        finalised_board = play_game_with_policy_update!(loss, params, model, opt)
        println("******************************self play report trends:***********************")
        display(finalised_board)
        play_max = [maximum(play_game_with_policy(model)) for i in 1:100]
        cm = countmap(play_max)
        display(sort(collect(cm), by = x->x[1]))
        push!(res, mean(play_max))
        println("$(mean(res)) $(res[end-7:end])")
        println("*****************************************************")
    end
end

# res = [6.0 for i in 1:8]
@time train_many_times!(res, 1200, loss, params, model, opt)