text:
push!(LOAD_PATH, "Modules")
using Colors
using ColorVectorSpace
using IntervalSets
using Revise
using JuliRay
using BasicBSpline
Base.@irrational ° 0.0174532925199432957692369076848861271344 (big(pi)/big(180))

function Smooth(x::Real)
    f(x) = ifelse(x>0,exp(-1/x),0)
    return f(2x/√3)/(f(2x/√3)+f(2/√3-2x/√3))
end
function Smooth(a::Real, b::Real, x::Real)
    return Smooth((x-a)/(b-a))
end
function Smooth(dict::Dict{A,B} where A <: Real where B <: Any, t::Real)
    T=sort(collect(keys(dict)))
    L=length(T)
    V=[dict[T[i]] for i in 1:L]

    if t ≤ T[1]
        return V[1]
    end
    for i in 2:L
        if t ≤ T[i]
            return V[i-1]+(V[i]-V[i-1])*Smooth(T[i-1],T[i],t)
        end
    end
    return V[end]
end

_p = 3
_k = Knots(1:13)
_P = BSplineSpace(_p,_k)
_a = [[[rand()-1/2,rand()-1/2,rand()-1/2] for _ in 1:dim(_P)-_p] for i in 1:4, j in 1:4]
_a = [vcat(_a[i,j],_a[i,j][1:_p]) for i in 1:4, j in 1:4]
_M = [BSplineManifold([_P], _a[i,j]) for i in 1:4, j in 1:4]

N=20
for i in 0:N-1
    t = _k[1+_p] + (_k[end-_p] - _k[1+_p]) * i/N
    k = 4*Knots([0,1])
    p = 3
    P = BSplineSpace(p,k)
    a0 = [[i-2.5,j-2.5,0]+mapping(_M[i,j],[t]) for i in 1:4, j in 1:4]
    a2 = [[0,0,Float64((i∈2:3) & (j∈2:3))] for i in 1:4, j in 1:4]
    a = a0+0*a2
    M = BSplineManifold([P,P],a)

    𝒑(u) = mapping(M,u)
    D1 = 0..1
    D2 = 0..1
    mesh = (40,40)
    nn = 5 # checker mesh
    rgb(u) = Float64(xor(sin(2π*nn*u[1])<0,sin(2π*nn*u[2])>0))*RGB(1,0.6,0.2) + Float64(xor(sin(2π*nn*u[1])>0,sin(2π*nn*u[2])>0))*RGB(1,0.2,0.2)
    ft(u) = FT(0.1,0.1)
    surface = ColoredParametricSurface(𝒑, D1, D2, mesh=mesh, rgb=rgb, ft=ft, smooth = false)

    ctrlpts = rgbColor(csgUnion([Sphere(p,0.05) for p in [a...]]),RGB(0.02,0.02,0.02))
    ctrlsgts1 = rgbColor(csgUnion([
        Cylinder(a[1,1],a[2,1],0.02),Cylinder(a[2,1],a[3,1],0.02),Cylinder(a[3,1],a[4,1],0.02),
        Cylinder(a[1,2],a[2,2],0.02),Cylinder(a[2,2],a[3,2],0.02),Cylinder(a[3,2],a[4,2],0.02),
        Cylinder(a[1,3],a[2,3],0.02),Cylinder(a[2,3],a[3,3],0.02),Cylinder(a[3,3],a[4,3],0.02),
        Cylinder(a[1,4],a[2,4],0.02),Cylinder(a[2,4],a[3,4],0.02),Cylinder(a[3,4],a[4,4],0.02),
        ]), RGB(0.1,0.1,0.1))
    ctrlsgts2 = rgbColor(csgUnion([
        Cylinder(a[1,1],a[1,2],0.02),Cylinder(a[1,2],a[1,3],0.02),Cylinder(a[1,3],a[1,4],0.02),
        Cylinder(a[2,1],a[2,2],0.02),Cylinder(a[2,2],a[2,3],0.02),Cylinder(a[2,3],a[2,4],0.02),
        Cylinder(a[3,1],a[3,2],0.02),Cylinder(a[3,2],a[3,3],0.02),Cylinder(a[3,3],a[3,4],0.02),
        Cylinder(a[4,1],a[4,2],0.02),Cylinder(a[4,2],a[4,3],0.02),Cylinder(a[4,3],a[4,4],0.02),
        ]), RGB(0.1,0.1,0.1))
    ctrlsgts = csgUnion(ctrlsgts1, ctrlsgts2)

    lng = 20°
    lat = 30°
    lookat = [0,0,0.25]*0
    zoom = 0.25
    render(csgUnion(surface,ctrlpts,ctrlsgts),camera=LngLatCamera(lng=lng,lat=lat,pers=0.2,zoom=zoom,lookat=lookat, width=500,height=300),name="DancingBezierSurface",index=i+1)
end

T = 0.2

N=40
for i in 0:N
    t = i/N
    tt = _k[1+_p] + (_k[end-_p] - _k[1+_p]) * T
    k = 4*Knots([0,1])
    p = 3
    P = BSplineSpace(p,k)
    a0 = [[i-2.5,j-2.5,0]+0.6*mapping(_M[i,j],[tt]) for i in 1:4, j in 1:4]
    a2 = [[0,0,Float64((i∈2:3) & (j∈2:3))] for i in 1:4, j in 1:4]
    a = a0+1.5*a2
    M = BSplineManifold([P,P],a)

    𝒑(u) = mapping(M,u)
    D1 = 0..1
    D2 = 0..1
    mesh = (40,40)
    nn = 5 # checker mesh
    rgb(u) = Float64(xor(sin(2π*nn*u[1])<0,sin(2π*nn*u[2])>0))*RGB(0.3,0.3,0.3) + Float64(xor(sin(2π*nn*u[1])>0,sin(2π*nn*u[2])>0))*RGB(0.5,0.5,0.5)
    ft(u) = FT(0.1,0.1)
    surface = ColoredParametricSurface(𝒑, D1, D2, mesh=mesh, rgb=rgb, ft=ft, smooth = false)

    ctrlpts = rgbColor(csgUnion([Sphere(p,0.05) for p in [a...]]),RGB(0,0.01,0.02))
    ctrlsgts1 = rgbColor(csgUnion([
        Cylinder(a[1,1],a[2,1],0.02),Cylinder(a[2,1],a[3,1],0.02),Cylinder(a[3,1],a[4,1],0.02),
        Cylinder(a[1,2],a[2,2],0.02),Cylinder(a[2,2],a[3,2],0.02),Cylinder(a[3,2],a[4,2],0.02),
        Cylinder(a[1,3],a[2,3],0.02),Cylinder(a[2,3],a[3,3],0.02),Cylinder(a[3,3],a[4,3],0.02),
        Cylinder(a[1,4],a[2,4],0.02),Cylinder(a[2,4],a[3,4],0.02),Cylinder(a[3,4],a[4,4],0.02),
        ]), RGB(0,0.1,0.2))
    ctrlsgts2 = rgbColor(csgUnion([
        Cylinder(a[1,1],a[1,2],0.02),Cylinder(a[1,2],a[1,3],0.02),Cylinder(a[1,3],a[1,4],0.02),
        Cylinder(a[2,1],a[2,2],0.02),Cylinder(a[2,2],a[2,3],0.02),Cylinder(a[2,3],a[2,4],0.02),
        Cylinder(a[3,1],a[3,2],0.02),Cylinder(a[3,2],a[3,3],0.02),Cylinder(a[3,3],a[3,4],0.02),
        Cylinder(a[4,1],a[4,2],0.02),Cylinder(a[4,2],a[4,3],0.02),Cylinder(a[4,3],a[4,4],0.02),
        ]), RGB(0,0.1,0.2))
    ctrlsgts = csgUnion(ctrlsgts1)

    _k_ = 4*Knots(0,1)
    _p_ = 3
    _P_ = BSplineSpace(_p_, _k_)
    _M1_ = [BSplineManifold([_P_], [a[1,i],a[2,i],a[3,i],a[4,i]]) for i in 1:4]
    c1 = rgbColor(csgUnion(ParametricCurve(t->mapping(_M1_[i],[t]),0..1,0.02) for i in 1:4), RGB(0,0.8,1))

    la1 = [mapping(_M1_[i],[t]) for i in 1:4]
    lM1 = BSplineManifold([_P_], la1)
    ca1 = rgbColor(csgUnion([Sphere(p,0.05) for p in [la1...]]),RGB(0.02,0.01,0))
    ctr1 = rgbColor(csgUnion([
        Cylinder(la1[1],la1[2],0.02),Cylinder(la1[2],la1[3],0.02),Cylinder(la1[3],la1[4],0.02),
        ]), RGB(0.2,0.1,0))
    lc1 = rgbColor(ParametricCurve(t->mapping(lM1,[t]),0..1,0.02), RGB(1,0.5,0.0))

    lng = 20°
    lat = 30°
    lookat = [0,0,0.1]
    zoom = 0.25
    render(csgUnion(surface,ctrlpts,ctrlsgts, c1, lc1, ca1, ctr1),camera=LngLatCamera(lng=lng,lat=lat,pers=0.2,zoom=zoom,lookat=lookat, width=600,height=400),name="BezierSurfaceLocus1",index=i+1)
end

for i in 0:N
    t = i/N
    tt = _k[1+_p] + (_k[end-_p] - _k[1+_p]) * T
    k = 4*Knots([0,1])
    p = 3
    P = BSplineSpace(p,k)
    a0 = [[i-2.5,j-2.5,0]+0.6*mapping(_M[i,j],[tt]) for i in 1:4, j in 1:4]
    a2 = [[0,0,Float64((i∈2:3) & (j∈2:3))] for i in 1:4, j in 1:4]
    a = a0+1.5*a2
    M = BSplineManifold([P,P],a)

    𝒑(u) = mapping(M,u)
    D1 = 0..1
    D2 = 0..1
    mesh = (40,40)
    nn = 5 # checker mesh
    rgb(u) = Float64(xor(sin(2π*nn*u[1])<0,sin(2π*nn*u[2])>0))*RGB(0.3,0.3,0.3) + Float64(xor(sin(2π*nn*u[1])>0,sin(2π*nn*u[2])>0))*RGB(0.5,0.5,0.5)
    ft(u) = FT(0.1,0.1)
    surface = ColoredParametricSurface(𝒑, D1, D2, mesh=mesh, rgb=rgb, ft=ft, smooth = false)

    ctrlpts = rgbColor(csgUnion([Sphere(p,0.05) for p in [a...]]),RGB(0,0.01,0.02))
    ctrlsgts1 = rgbColor(csgUnion([
        Cylinder(a[1,1],a[2,1],0.02),Cylinder(a[2,1],a[3,1],0.02),Cylinder(a[3,1],a[4,1],0.02),
        Cylinder(a[1,2],a[2,2],0.02),Cylinder(a[2,2],a[3,2],0.02),Cylinder(a[3,2],a[4,2],0.02),
        Cylinder(a[1,3],a[2,3],0.02),Cylinder(a[2,3],a[3,3],0.02),Cylinder(a[3,3],a[4,3],0.02),
        Cylinder(a[1,4],a[2,4],0.02),Cylinder(a[2,4],a[3,4],0.02),Cylinder(a[3,4],a[4,4],0.02),
        ]), RGB(0,0.1,0.2))
    ctrlsgts2 = rgbColor(csgUnion([
        Cylinder(a[1,1],a[1,2],0.02),Cylinder(a[1,2],a[1,3],0.02),Cylinder(a[1,3],a[1,4],0.02),
        Cylinder(a[2,1],a[2,2],0.02),Cylinder(a[2,2],a[2,3],0.02),Cylinder(a[2,3],a[2,4],0.02),
        Cylinder(a[3,1],a[3,2],0.02),Cylinder(a[3,2],a[3,3],0.02),Cylinder(a[3,3],a[3,4],0.02),
        Cylinder(a[4,1],a[4,2],0.02),Cylinder(a[4,2],a[4,3],0.02),Cylinder(a[4,3],a[4,4],0.02),
        ]), RGB(0,0.1,0.2))
    ctrlsgts = csgUnion(ctrlsgts2)

    _k_ = 4*Knots(0,1)
    _p_ = 3
    _P_ = BSplineSpace(_p_, _k_)
    _M1_ = [BSplineManifold([_P_], [a[i,1],a[i,2],a[i,3],a[i,4]]) for i in 1:4]
    c1 = rgbColor(csgUnion(ParametricCurve(t->mapping(_M1_[i],[t]),0..1,0.02) for i in 1:4), RGB(0,0.8,1))

    la1 = [mapping(_M1_[i],[t]) for i in 1:4]
    lM1 = BSplineManifold([_P_], la1)
    ca1 = rgbColor(csgUnion([Sphere(p,0.05) for p in [la1...]]),RGB(0.02,0.01,0))
    ctr1 = rgbColor(csgUnion([
        Cylinder(la1[1],la1[2],0.02),Cylinder(la1[2],la1[3],0.02),Cylinder(la1[3],la1[4],0.02),
        ]), RGB(0.2,0.1,0))
    lc1 = rgbColor(ParametricCurve(t->mapping(lM1,[t]),0..1,0.02), RGB(1,0.5,0.0))

    lng = 20°
    lat = 30°
    lookat = [0,0,0.1]
    zoom = 0.25
    render(csgUnion(surface,ctrlpts,ctrlsgts, c1, lc1, ca1, ctr1),camera=LngLatCamera(lng=lng,lat=lat,pers=0.2,zoom=zoom,lookat=lookat, width=600,height=400),name="BezierSurfaceLocus2",index=i+1)
end

for i in 0:0
    tt = _k[1+_p] + (_k[end-_p] - _k[1+_p]) * T
    k = 4*Knots([0,1])
    p = 3
    P = BSplineSpace(p,k)
    a0 = [[i-2.5,j-2.5,0]+0.6*mapping(_M[i,j],[tt]) for i in 1:4, j in 1:4]
    a2 = [[0,0,Float64((i∈2:3) & (j∈2:3))] for i in 1:4, j in 1:4]
    a = a0+1.5*a2
    M = BSplineManifold([P,P],a)

    𝒑(u) = mapping(M,u)
    D1 = 0..1
    D2 = 0..1
    mesh = (40,40)
    nn = 5 # checker mesh
    rgb(u) = Float64(xor(sin(2π*nn*u[1])<0,sin(2π*nn*u[2])>0))*RGB(0.3,0.3,0.3) + Float64(xor(sin(2π*nn*u[1])>0,sin(2π*nn*u[2])>0))*RGB(0.5,0.5,0.5)
    ft(u) = FT(0.1,0.1)
    surface = ColoredParametricSurface(𝒑, D1, D2, mesh=mesh, rgb=rgb, ft=ft, smooth = false)

    ctrlpts = rgbColor(csgUnion([Sphere(p,0.05) for p in [a...]]),RGB(0,0.01,0.02))
    ctrlsgts1 = rgbColor(csgUnion([
        Cylinder(a[1,1],a[2,1],0.02),Cylinder(a[2,1],a[3,1],0.02),Cylinder(a[3,1],a[4,1],0.02),
        Cylinder(a[1,2],a[2,2],0.02),Cylinder(a[2,2],a[3,2],0.02),Cylinder(a[3,2],a[4,2],0.02),
        Cylinder(a[1,3],a[2,3],0.02),Cylinder(a[2,3],a[3,3],0.02),Cylinder(a[3,3],a[4,3],0.02),
        Cylinder(a[1,4],a[2,4],0.02),Cylinder(a[2,4],a[3,4],0.02),Cylinder(a[3,4],a[4,4],0.02),
        ]), RGB(0,0.1,0.2))
    ctrlsgts2 = rgbColor(csgUnion([
        Cylinder(a[1,1],a[1,2],0.02),Cylinder(a[1,2],a[1,3],0.02),Cylinder(a[1,3],a[1,4],0.02),
        Cylinder(a[2,1],a[2,2],0.02),Cylinder(a[2,2],a[2,3],0.02),Cylinder(a[2,3],a[2,4],0.02),
        Cylinder(a[3,1],a[3,2],0.02),Cylinder(a[3,2],a[3,3],0.02),Cylinder(a[3,3],a[3,4],0.02),
        Cylinder(a[4,1],a[4,2],0.02),Cylinder(a[4,2],a[4,3],0.02),Cylinder(a[4,3],a[4,4],0.02),
        ]), RGB(0,0.1,0.2))
    ctrlsgts = csgUnion(ctrlsgts1, ctrlsgts2)

    lng = 20°
    lat = 30°
    lookat = [0,0,0.1]
    zoom = 0.25
    render(csgUnion(surface,ctrlpts,ctrlsgts),camera=LngLatCamera(lng=lng,lat=lat,pers=0.2,zoom=zoom,lookat=lookat, width=600,height=400),name="BezierSurfaceLocus")
end

max_stars_repo_name: hyrodium/JuliRay

max_issues_repo_name: hyrodium/JuliRay

max_forks_repo_name: hyrodium/JuliRay


text:
using Test, AutomaticDocstrings, DeepDiffs

function testdoc(input, output)
    # Test function autodoc
    path, io = mktemp()
    write(io, input)
    close(io)
    autodoc(path)
    result = reduce(*, readlines(path, keep=true))
    if !occursin(output, result)
        println("Autodoc failed, input:   =====================================")
        println(input)
        # println("Result:     ==================================================")
        # println(result)
        # println("Expected:   ==================================================")
        # println(output)
        println("Diff:       ==================================================")
        println(deepdiff(output, result))
        return false
    end

    # Test macro @autodoc
    write(io, input)
    close(io)
    include(path)
    result = reduce(*, readlines(path, keep=true))
    occursin(output, result)
end

restore_defaults()

@testset "AutomaticDocstrings" begin
    @info "Testing AutomaticDOcstrings"

@test testdoc(
"""
@autodoc
function f(a, b; c)
end
""",
"""
\"\"\"
    f(a, b; c)

DOCSTRING

# Arguments:
- `a`: DESCRIPTION
- `b`: DESCRIPTION
- `c`: DESCRIPTION
\"\"\"
function f(a, b; c)
end
""")

@test testdoc(
"""
@autodoc
function f(a,
    b;
    c)
end
""",
"""
\"\"\"
    f(a, b; c)

DOCSTRING

# Arguments:
- `a`: DESCRIPTION
- `b`: DESCRIPTION
- `c`: DESCRIPTION
\"\"\"
function f(a,
    b;
    c)
end
""")

@test testdoc(
"""
@autodoc
f(a::A, b; c) where A = 5
""",
"""
\"\"\"
    f(a::A, b; c)

DOCSTRING

# Arguments:
- `a`: DESCRIPTION
- `b`: DESCRIPTION
- `c`: DESCRIPTION
\"\"\"
f(a::A, b; c) where A = 5
""")

@test testdoc(
"""
@autodoc
function f(a::A, b; c) where A
    5
end
""",
"""
\"\"\"
    f(a::A, b; c)

DOCSTRING

# Arguments:
- `a`: DESCRIPTION
- `b`: DESCRIPTION
- `c`: DESCRIPTION
\"\"\"
function f(a::A, b; c) where A
    5
end
""")

AutomaticDocstrings.options[:full_def] = false

@test testdoc(
"""
@autodoc
function f(a=4, b=LinRange(1,2,3); c="hello")
end
""",
"""
\"\"\"
    f(a, b, c)

DOCSTRING

# Arguments:
- `a`: DESCRIPTION
- `b`: DESCRIPTION
- `c`: DESCRIPTION
\"\"\"
function f(a=4, b=LinRange(1,2,3); c="hello")
end
""")

# Fewer that min_args
@test testdoc(
"""
@autodoc
function f(a=4, b=LinRange(1,2,3))
end
""",
"""
\"\"\"
    f(a, b)

DOCSTRING
\"\"\"
function f(a=4, b=LinRange(1,2,3))
end
""")

restore_defaults()
@test AutomaticDocstrings.options[:full_def]

# Struct

@test testdoc(
"""
@autodoc
struct Workspace{T1,T2,T3,T4,T5,T6}
    simple_input::T1
    simple_result::T2
    result::T3
    buffersetter::T4
    resultsetter::T5
    f::T6
    N::Int
end
""",
"""
\"\"\"
    Workspace{T1, T2, T3, T4, T5, T6}

DOCSTRING

# Arguments:
- `simple_input::T1`: DESCRIPTION
- `simple_result::T2`: DESCRIPTION
- `result::T3`: DESCRIPTION
- `buffersetter::T4`: DESCRIPTION
- `resultsetter::T5`: DESCRIPTION
- `f::T6`: DESCRIPTION
- `N::Int`: DESCRIPTION
\"\"\"
struct Workspace{T1,T2,T3,T4,T5,T6}
    simple_input::T1
    simple_result::T2
    result::T3
    buffersetter::T4
    resultsetter::T5
    f::T6
    N::Int
end
""")

@static if VERSION >= v"1.3"
import Base.@kwdef

@test testdoc(
"""
@autodoc
@kwdef struct Workspace{T1,T2,T3,T4,T5,T6}
    simple_input::T1
    simple_result::T2
    result::T3
    buffersetter::T4
    resultsetter::T5
    f::T6
    N::Int = 2
end
""",
"""
\"\"\"
    Workspace{T1, T2, T3, T4, T5, T6}

DOCSTRING

# Arguments:
- `simple_input::T1`: DESCRIPTION
- `simple_result::T2`: DESCRIPTION
- `result::T3`: DESCRIPTION
- `buffersetter::T4`: DESCRIPTION
- `resultsetter::T5`: DESCRIPTION
- `f::T6`: DESCRIPTION
- `N::Int = 2`: DESCRIPTION
\"\"\"
@kwdef struct Workspace{T1,T2,T3,T4,T5,T6}
    simple_input::T1
    simple_result::T2
    result::T3
    buffersetter::T4
    resultsetter::T5
    f::T6
    N::Int = 2
end
""")
end

end

ext: jl


text:
abstract type AbstractAllocator <: AbstractController end
abstract type StaticAllocator <: AbstractAllocator end

include("PseudoInverseAllocator.jl")


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Deselect <: Edit
        col::Array{Symbol,1}
        operation::String
    end

# Arguments
- `col::Array{Symbol,1}`: name of column containing data to remove
- `operation::String`: how to determine what to drop
"""
mutable struct Deselect <: Edit
    "name of column containing data to remove"
    col::Array{Symbol,1}
    "how to determine what to drop"
    operation::String
end

function Deselect(; col, operation, )
    Deselect(col, operation, )
end

"""Get [`Deselect`](@ref) `col`."""
get_col(value::Deselect) = value.col
"""Get [`Deselect`](@ref) `operation`."""
get_operation(value::Deselect) = value.operation

"""Set [`Deselect`](@ref) `col`."""
set_col!(value::Deselect, val) = value.col = val
"""Set [`Deselect`](@ref) `operation`."""
set_operation!(value::Deselect, val) = value.operation = val

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
function promote_utype(t1::TypeOfUnknown, t2::TypeOfUnknown)
    return OpaqueType()
end

function compatible(s1::Nullable{Int64}, s2::Nullable{Int64})
    if isnull(s1) || isnull(s2)
        return true
    end

    s1 = get(s1)
    s2 = get(s2)

    s1 == s2 || s1 == 1 || s2 == 1
end

function broadcast(s1::Nullable{Int64}, s2::Nullable{Int64})
    if isnull(s1) || isnull(s2)
        return nothing
    end

    s1v = get(s1)
    s2v = get(s2)

    if s1v == s2v
        return s1
    end

    if s1v == 1
        return s2
    else
        return s1
    end
end

function promote_utype(atype1::ArrayType, atype2::ArrayType)
    shape1 = atype1.shape
    shape2 = atype2.shape

    if isnull(shape1) || isnull(shape2)
        res_shape = Nullable{Array{Nullable{Int64}, 1}}()
    else
        shape1 = get(shape1)
        shape2 = get(shape2)

        if length(shape1) != length(shape2)
            extend_dims = abs(length(shape1) - length(shape2))
            extend_shape = [Nullable{Int64}(1) for i in 1:extend_dims]

            if length(shape1) < length(shape2)
                shape1 = [shape1..., extend_shape...]
            else
                shape2 = [shape2..., extend_shape...]
            end
        end

        @assert all(compatible(s1, s2) for (s1, s2) in zip(shape1, shape2))
        res_shape = map((x) -> broadcast(x[1], x[2]), zip(shape1, shape2))
        res_shape = Nullable{Array{Nullable{Int64}, 1}}(res_shape)
    end

    eltype1 = atype1.eltype
    eltype2 = atype2.eltype

    if isnull(eltype1) || isnull(eltype2)
        res_eltype = Nullable{DataType}()
    else
        res_eltype = Nullable{DataType}(promote_type(get(eltype1), get(eltype2)))
    end

    ArrayType(res_shape, res_eltype)
end

# Some machinery for intercepting broadcast with `ArrayExpr` objects.
# Usage of e.g. `.*` transforms into a broadcast call with an anonymous function with `*` inside.
# Since we would like to distinguish between `.*` and `*` in array expressions,
# we wrap `ArrayExpr` objects into a special object, add arithmetic methods for it,
# call the anonymous function, and unwrap it back.

struct Broadcasted{T <: ArrayExpr}
    val::T
end

to_bc(x::Broadcasted) = x
to_bc(x) = Broadcasted(x)

from_bc(x::Broadcasted) = x.val

to_aexpr(x::ArrayExpr) = x
to_aexpr(x) = convert(ArrayExpr, x)

call_node(op, a, b) = CallNode(op, [a, b], promote_utype(expr_type(a), expr_type(b)))

for op in (:+, :-, :*, :/, :^)
    @eval begin
        Base.$op(a::Broadcasted, b::Broadcasted) =
            to_bc(call_node(Symbol(".", $op), from_bc(a), from_bc(b)))
        Base.$op(a::Broadcasted, b) = $op(a, to_bc(to_aexpr(b)))
        Base.$op(a, b::Broadcasted) = $op(to_bc(to_aexpr(a)), b)
    end
end

Base.Broadcast.containertype(::ArrayExpr) = ArrayExpr
Base.Broadcast.promote_containertype(::Type{ArrayExpr}, ::Type{ArrayExpr}) = ArrayExpr
# We would want to only handle the numbers and arrays of numbers here.
# But currently:
# - `containertype()` of a scalar type results in `Any`, so we cannot add a
#   `promote_containertype()` method for `ArrayExpr` and `Number` - it just will not be called.
# - there is already a `promote_containertype()` for `Any` and `Array` which would
#   intercept anything more specialized than the methods below.
# That's why constructions like e.g. `[x, y] .== x` where `x` is of some `ArrayExpr` subtype
# will fall into our custom broadcast and fail miserably.
Base.Broadcast.promote_containertype(::Type{ArrayExpr}, ::Type{Any}) = ArrayExpr
Base.Broadcast.promote_containertype(::Type{Any}, ::Type{ArrayExpr}) = ArrayExpr
Base.Broadcast.promote_containertype(::Type{ArrayExpr}, ::Type{Array}) = ArrayExpr
Base.Broadcast.promote_containertype(::Type{Array}, ::Type{ArrayExpr}) = ArrayExpr

# Custom broadcast on our container type
function Base.Broadcast.broadcast_c(f, ::Type{ArrayExpr}, as...)
    from_bc(f(to_bc.(to_aexpr.(as))...))
end

for op in (:+, :-, :*, :/, :^)
    @eval begin
        Base.$op(x::ArrayExpr, y::ArrayExpr) = call_node(Symbol($op), x, y)
    end

    for ext_type in (Number, AbstractArray)
        @eval begin
            Base.$op(x::ArrayExpr, y::$ext_type) = ($op)(x, to_aexpr(y))
            Base.$op(x::$ext_type, y::ArrayExpr) = ($op)(to_aexpr(x), y)
        end
    end
end

import Base: size, reshape, permutedims, repeat

function size_expr(arr::ArrayExpr)
    CallNode(:size, [arr], ArrayType([1], Int64))
end

function size_expr(arr::ArrayExpr, dim::ArrayExpr)
    CallNode(:size, [arr, dim], ArrayType([], Int64))
end

function size(arr::ArrayExpr)
    atype = expr_type(arr)
    if !isnull(atype.shape)
        shape = get(atype.shape)

        if !any(map(isnull, shape))
            shape = map(get, shape)
            return tuple(shape...)
        end
    end

    size_expr(arr)
end

function size(arr::ArrayExpr, dim)
    atype = expr_type(arr)
    if !isnull(atype.shape)
        shape = get(atype.shape)
        if !isnull(shape[dim])
            return get(shape[dim])
        end
    end
    size_expr(arr, to_aexpr(dim))
end

function size(arr, dim::ArrayExpr)
    size_expr(to_aexpr(arr), dim)
end

function reshape(arr::ArrayExpr, dims::ArrayExpr)
    CallNode(:reshape, [arr, dims], ArrayType(nothing, expr_type(arr).eltype))
end

normalize_dims(dim::Number) = [dim]
normalize_dims(dim::AbstractArray) = dim
normalize_dims{T <: Number,N}(dim::Tuple{Vararg{T,N}}) = collect(dim)
normalize_dims(dims...) = collect(dims)

function reshape(arr::ArrayExpr, dims...)
    dims = normalize_dims(dims)
    atype = expr_type(arr)

    if !isnull(atype.shape)
        shape = get(atype.shape)
        if !any(map(isnull, shape))
            shape = map(get, shape)
            @assert prod(shape) == prod(dims)
        end
    end

    CallNode(:reshape, [arr, to_aexpr(dims)], ArrayType(dims, atype.eltype))
end

function reshape(arr, dims::ArrayExpr)
    CallNode(:reshape, [to_aexpr(arr), dims], ArrayType(nothing, atype.eltype))
end

function permutedims(arr::ArrayExpr, perm)
    perm = collect(perm)
    @assert ndims(perm) == 1
    atype = expr_type(arr)

    if !isnull(atype.shape)
        shape = get(atype.shape)
        @assert length(shape) == length(perm)
        new_shape = shape[perm]
    else
        new_shape = nothing
    end

    CallNode(:permutedims, [arr, to_aexpr(perm)], ArrayType(new_shape, atype.eltype))
end

function permutedims(arr, perm::ArrayExpr)
    CallNode(:permutedims, [to_aexpr(arr), perm], ArrayType(nothing, atype.eltype))
end

normalize_outer(rep::Number) = [rep]
normalize_outer(rep::AbstractArray) = rep
normalize_outer{T <: Number,N}(rep::Tuple{Vararg{T,N}}) = collect(rep)

function repeat(arr::ArrayExpr; outer=nothing)
    if outer === nothing
        return arr
    end

    outer = normalize_outer(outer)
    atype = expr_type(arr)

    if !isnull(atype.shape)
        shape = get(atype.shape)
        @assert length(shape) == length(outer)

        new_shape = [isnull(x) ? nothing : get(x) * rep for (x, rep) in zip(shape, outer)]
    else
        new_shape = nothing
    end

    CallNode(:repeat_outer, [arr, to_aexpr(outer)], ArrayType(new_shape, atype.eltype))
end

for builtin_function in (:cosh, :fft, :ifft, :abs, :sqrt, :sum)
    @eval begin
        import Base: $builtin_function

        function $builtin_function(arg::ArrayExpr, args...)
            # This check is not necessary, [x, args...] should work anyway,
            # but it currently doesn't because of bug #17003
            if length(args) > 0
                args = [convert(ArrayExpr, arg) for arg in args]
                args = [arg, args...]
            else
                args = [arg]
            end
            CallNode($(Meta.quot(builtin_function)), args, expr_type(arg))
        end
    end
end

function unpack(aexpr::ArrayExpr)
    etype = expr_type(aexpr)
    if typeof(etype) != TupleType
        return aexpr
    else
        return tuple((
            CallNode(:unpack, [aexpr, convert(ArrayExpr, i)], typ)
            for (i, typ) in enumerate(etype.elements))...)
    end
end

as_aexpr(val) = convert(ArrayExpr, val)

function random_uniform(state::ArrayExpr, min, max, dims...)
    atype = ArrayType(collect(dims), Float64)
    dims = map(as_aexpr, dims)
    arr, new_state = unpack(
        CallNode(:random_uniform, [state, dims...], TupleType([atype, RandomStateType()])))
    arr .* (max .- min) .+ min, new_state
end

function random_normal(state::ArrayExpr, mean, std, dims...)
    atype = ArrayType(collect(dims), Float64)
    dims = map(as_aexpr, dims)
    arr, new_state = unpack(
        CallNode(:random_normal, [state, dims...], TupleType([atype, RandomStateType()])))
    arr .* std .+ mean, new_state
end


text:
using Documenter
import Pkg

# Fix for https://github.com/trixi-framework/Trixi.jl/issues/668
if (get(ENV, "CI", nothing) != "true") && (get(ENV, "TRIXI_DOC_DEFAULT_ENVIRONMENT", nothing) != "true")
    push!(LOAD_PATH, dirname(@__DIR__))
end

using Trixi
using Trixi2Vtk

# Get Trixi root directory
trixi_root_dir = dirname(@__DIR__)

include(joinpath(trixi_root_dir, "docs", "literate", "make.jl"))

# Copy list of authors to not need to synchronize it manually
authors_text = read(joinpath(trixi_root_dir, "AUTHORS.md"), String)
authors_text = replace(authors_text, "in the [LICENSE.md](LICENSE.md) file" => "under [License](@ref)")
write(joinpath(@__DIR__, "src", "authors.md"), authors_text)

# Define module-wide setups such that the respective modules are available in doctests
DocMeta.setdocmeta!(Trixi,     :DocTestSetup, :(using Trixi);     recursive=true)
DocMeta.setdocmeta!(Trixi2Vtk, :DocTestSetup, :(using Trixi2Vtk); recursive=true)

# Create tutorials for the following files:
# Normal structure: "title" => "filename.jl"
# If there are several files for one topic and one folder, the structure is:
#   "title" => ["subtitle 1" => ("folder 1", "filename 1.jl"),
#               "subtitle 2" => ("folder 2", "filename 2.jl")]
files = [
    "Adding a new equation" => ["Scalar conservation law" => ("adding_new_equations", "cubic_conservation_law.jl"),
                                "Nonconservative equation" => ("adding_new_equations", "nonconservative_advection.jl")],
    "Differentiable programming" => "differentiable_programming.jl",
    "Unstructured meshes with HOHQMesh.jl" => "hohqmesh_tutorial.jl",
    ]
tutorials = create_tutorials(files)

# Make documentation
makedocs(
    # Specify modules for which docstrings should be shown
    modules = [Trixi, Trixi2Vtk],
    # Set sitename to Trixi
    sitename="Trixi.jl",
    # Provide additional formatting options
    format = Documenter.HTML(
        # Disable pretty URLs during manual testing
        prettyurls = get(ENV, "CI", nothing) == "true",
        # Explicitly add favicon as asset
        assets = ["assets/favicon.ico"],
        # Set canonical URL to GitHub pages URL
        canonical = "https://trixi-framework.github.io/Trixi.jl/stable"
    ),
    # Explicitly specify documentation structure
    pages = [
        "Home" => "index.md",
        "Getting started" => [
            "Overview" => "overview.md",
            "Visualization" => "visualization.md",
        ],
        "Tutorials" => tutorials,
        "Basic building blocks" => [
            "Meshes" => [
                "Tree mesh" => joinpath("meshes", "tree_mesh.md"),
                "Structured mesh" => joinpath("meshes", "structured_mesh.md"),
                "Unstructured mesh" => joinpath("meshes", "unstructured_quad_mesh.md"),
                "P4est-based mesh" => joinpath("meshes", "p4est_mesh.md"),
                "Simplicial mesh" => joinpath("meshes", "mesh_data_meshes.md"),
            ],
            "Time integration" => "time_integration.md",
            "Callbacks" => "callbacks.md",
        ],
        "Advanced topics & developers" => [
            "Conventions" =>"conventions.md",
            "Development" => "development.md",
            "GitHub & Git" => "github-git.md",
            "Style guide" => "styleguide.md",
            "Testing" => "testing.md",
            "Performance" => "performance.md",
            "Parallelization" => "parallelization.md",
        ],
        "Troubleshooting and FAQ" => "troubleshooting.md",
        "Reference" => [
                        "Trixi.jl" => "reference-trixi.md",
                        "Trixi2Vtk.jl" => "reference-trixi2vtk.md"
                       ],
        "Authors" => "authors.md",
        "Contributing" => "contributing.md",
        "License" => "license.md"
    ],
    strict = true # to make the GitHub action fail when doctests fail, see https://github.com/neuropsychology/Psycho.jl/issues/34
)

deploydocs(
    repo = "github.com/trixi-framework/Trixi.jl",
    devbranch = "main",
    push_preview = true
)


text:
module TestSymmetry

using LinearAlgebra, SparseArrays, Test
using SumOfSquares

function test_linsolve()
    x = [1, 2]
    for A in [
        [1 0 2 3
         0 1 3 -2],
        [1 2 0 3
         0 3 1 -2],
    ]
        b = A' * x
        @test Certificate.Symmetry._linsolve(A, b, Symmetry._RowEchelonMatrix()) ≈ x
        B = float.(A)
        for i in 1:size(B, 1)
            B[i, :] = normalize(B[i, :])
        end
        b = B' * x
        @test Certificate.Symmetry.__linsolve(B, b, Symmetry._OrthogonalMatrix()) ≈ x
        @test Certificate.Symmetry.__linsolve(sparse(B), b, Symmetry._OrthogonalMatrix()) ≈ x
    end
end

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$name", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
end

end

TestSymmetry.runtests()

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
using QCDNUM

export QCDNUMGrid, QCDNUMParameters
export SPLINTParameters

"""
    QCDNUMGrid

Struct for holding the QCDNUM grid parameters.
"""
@with_kw struct QCDNUMGrid
    "x boundaries of grid"
    x_min::Vector{Float64}
    "x grid weights"
    x_weights::Vector{Int32}
    "number of x grid boundaries"
    x_num_bounds::Integer = size(x_min)[1]
    "number of x grid points"
    nx::Integer
    "qq boundaries of grid"
    qq_bounds::Vector{Float64}
    "qq grid weights"
    qq_weights::Vector{Float64}
    "number of qq grid boundaries"
    qq_num_bounds::Integer = size(qq_bounds)[1]
    "number of qq grid points"
    nq::Integer
    "degree of spline interpolation used"
    spline_interp::Integer = 3
end

"""
   QCDNUMParameters

Struct for holding all QCDNUM Parameters.
"""
@with_kw struct QCDNUMParameters
    "order of evolution in pQCD"
    order::Integer
    "coupling constant at starting scale"
    α_S::Float64
    "starting scale"
    q0::Float64
    "QCDNUMGrid"
    grid::QCDNUMGrid
    "number of fixed flavours in FFNS"
    n_fixed_flav::Integer = 5
    "charm threshold as index of qq grid (VFNS only)"
    iqc::Integer = 1
    "bottom threshold as index of qq grid (VFNS only)"
    iqb::Integer = 1
    "top threshold as index of qq grid (VFNS only)"
    iqt::Integer = 1
    "weight table type (1<=>unpolarised)"
    weight_type::Integer = 1
end

"""
    SplineAddresses

Lookup table for addresses of different
structure function splines.
"""
@with_kw struct SplineAddresses
    F2up::Integer = 1
    F2dn::Integer = 2
    F3up::Integer = 3
    F3dn::Integer = 4
    FLup::Integer = 5
    FLdn::Integer = 6
    F_eP::Integer = 7
    F_eM::Integer = 8
end

"""
    SPLINTParameters

Struct for storage of parameters used
with SPLINT package of QCDNUM.
"""
@with_kw struct SPLINTParameters
    "number of words in memory for user space"
    nuser::Integer = 10
    "number of steps in x"
    nsteps_x::Integer = 5
    "number of steps in qq"
    nsteps_q::Integer = 10
    "number of nodes in x"
    nnodes_x::Integer = 100
    "number of nodes in qq"
    nnodes_q::Integer = 100
    "rs constraint"
    rs::Float64 = 318.0
    "cut on rs"
    rscut::Float64 = 370.0
    "spline addresses"
    spline_addresses::SplineAddresses = SplineAddresses()
end


text:
# implementation of the GPUCompiler interfaces for generating GCN code

## target

export GCNCompilerTarget

Base.@kwdef struct GCNCompilerTarget <: AbstractCompilerTarget
    dev_isa::String
end

llvm_triple(::GCNCompilerTarget) = "amdgcn-amd-amdhsa"

function llvm_machine(target::GCNCompilerTarget)
    triple = llvm_triple(target)
    t = Target(triple=triple)

    cpu = target.dev_isa
    feat = ""
    optlevel = LLVM.API.LLVMCodeGenLevelDefault
    reloc = LLVM.API.LLVMRelocPIC
    tm = TargetMachine(t, triple, cpu, feat, optlevel, reloc)
    asm_verbosity!(tm, true)

    return tm
end

## job

# TODO: encode debug build or not in the compiler job
#       https://github.com/JuliaGPU/CUDAnative.jl/issues/368
runtime_slug(job::CompilerJob{GCNCompilerTarget}) = "gcn-$(job.target.dev_isa)"

const gcn_intrinsics = () # TODO: ("vprintf", "__assertfail", "malloc", "free")
isintrinsic(::CompilerJob{GCNCompilerTarget}, fn::String) = in(fn, gcn_intrinsics)

function process_kernel!(job::CompilerJob{GCNCompilerTarget}, mod::LLVM.Module, kernel::LLVM.Function)
    kernel = wrap_entry!(job, mod, kernel)
    # AMDGPU kernel calling convention
    callconv!(kernel, LLVM.API.LLVMCallConv(91))
    kernel
end

function add_lowering_passes!(job::CompilerJob{GCNCompilerTarget}, pm::LLVM.PassManager)
    add!(pm, ModulePass("LowerThrowExtra", lower_throw_extra!))
    add!(pm, FunctionPass("FixAllocaAddrspace", fix_alloca_addrspace!))
end

function lower_throw_extra!(mod::LLVM.Module)
    job = current_job::CompilerJob
    ctx = context(mod)
    changed = false
    @timeit_debug to "lower throw (extra)" begin

    throw_functions = [
        r"julia_bounds_error.*",
        r"julia_throw_boundserror.*",
        r"julia_error_if_canonical_getindex.*",
        r"julia_error_if_canonical_setindex.*",
        r"julia___subarray_throw_boundserror.*",
    ]

    for f in functions(mod)
        f_name = LLVM.name(f)
        for fn in throw_functions
            if occursin(fn, f_name)
                for use in uses(f)
                    call = user(use)::LLVM.CallInst

                    # replace the throw with a trap
                    let builder = Builder(ctx)
                        position!(builder, call)
                        emit_exception!(builder, f_name, call)
                        dispose(builder)
                    end

                    # remove the call
                    call_args = collect(operands(call))[1:end-1] # last arg is function itself
                    unsafe_delete!(LLVM.parent(call), call)

                    # HACK: kill the exceptions' unused arguments
                    for arg in call_args
                        # peek through casts
                        if isa(arg, LLVM.AddrSpaceCastInst)
                            cast = arg
                            arg = first(operands(cast))
                            isempty(uses(cast)) && unsafe_delete!(LLVM.parent(cast), cast)
                        end

                        if isa(arg, LLVM.Instruction) && isempty(uses(arg))
                            unsafe_delete!(LLVM.parent(arg), arg)
                        end
                    end

                    changed = true
                end

                @compiler_assert isempty(uses(f)) job
            end
        end
    end

    end
    return changed
end
function fix_alloca_addrspace!(fn::LLVM.Function)
    changed = false
    alloca_as = 5
    ctx = context(fn)

    for bb in blocks(fn)
        for insn in instructions(bb)
            if isa(insn, LLVM.AllocaInst)
                ty = llvmtype(insn)
                ety = eltype(ty)
                addrspace(ty) == alloca_as && continue

                new_insn = nothing
                Builder(ctx) do builder
                    position!(builder, insn)
                    _alloca = alloca!(builder, ety, name(insn))
                    new_insn = addrspacecast!(builder, _alloca, ty)
                end
                replace_uses!(insn, new_insn)
                unsafe_delete!(LLVM.parent(insn), insn)
            end
        end
    end

    return changed
end

function emit_trap!(job::CompilerJob{GCNCompilerTarget}, builder, mod, inst)
    ctx = context(mod)
    trap = if haskey(functions(mod), "llvm.trap")
        functions(mod)["llvm.trap"]
    else
        LLVM.Function(mod, "llvm.trap", LLVM.FunctionType(LLVM.VoidType(ctx)))
    end
    if Base.libllvm_version < v"9"
        rl_ft = LLVM.FunctionType(LLVM.Int32Type(ctx),
                                  [LLVM.Int32Type(ctx)])
        rl = if haskey(functions(mod), "llvm.amdgcn.readfirstlane")
            functions(mod)["llvm.amdgcn.readfirstlane"]
        else
            LLVM.Function(mod, "llvm.amdgcn.readfirstlane", rl_ft)
        end
        # FIXME: Early versions of the AMDGPU target fail to skip machine
        # blocks with certain side effects when EXEC==0, except when certain
        # criteria are met within said block. We emit a v_readfirstlane_b32
        # instruction here, as that is sufficient to trigger a skip. Without
        # this, the target will only attempt to do a "masked branch", which
        # only works on vector instructions (trap is a scalar instruction, and
        # therefore it is executed even when EXEC==0).
        rl_val = call!(builder, rl, [ConstantInt(Int32(32), ctx)])
        rl_bc = inttoptr!(builder, rl_val, LLVM.PointerType(LLVM.Int32Type(ctx)))
        store!(builder, rl_val, rl_bc)
    end
    call!(builder, trap)
end

# manual implementation of byval, as the backend doesn't support it for kernel args
# https://reviews.llvm.org/D79744
function wrapper_type(julia_t::Type, codegen_t::LLVMType)::LLVMType
    if !isbitstype(julia_t)
        # don't pass jl_value_t by value; it's an opaque structure
        return codegen_t
    elseif isa(codegen_t, LLVM.PointerType) && !(julia_t <: Ptr)
        # we didn't specify a pointer, but codegen passes one anyway.
        # make the wrapper accept the underlying value instead.
        return eltype(codegen_t)
    else
        return codegen_t
    end
end
# generate a kernel wrapper to fix & improve argument passing
function wrap_entry!(job::CompilerJob, mod::LLVM.Module, entry_f::LLVM.Function)
    ctx = context(mod)
    entry_ft = eltype(llvmtype(entry_f)::LLVM.PointerType)::LLVM.FunctionType
    @compiler_assert return_type(entry_ft) == LLVM.VoidType(ctx) job

    # filter out types which don't occur in the LLVM function signatures
    sig = Base.signature_type(job.source.f, job.source.tt)::Type
    julia_types = Type[]
    for dt::Type in sig.parameters
        if !isghosttype(dt) && (VERSION < v"1.5.0-DEV.581" || !Core.Compiler.isconstType(dt))
            push!(julia_types, dt)
        end
    end

    # generate the wrapper function type & definition
    wrapper_types = LLVM.LLVMType[wrapper_type(julia_t, codegen_t)
                                  for (julia_t, codegen_t)
                                  in zip(julia_types, parameters(entry_ft))]
    wrapper_fn = LLVM.name(entry_f)
    LLVM.name!(entry_f, wrapper_fn * ".inner")
    wrapper_ft = LLVM.FunctionType(LLVM.VoidType(ctx), wrapper_types)
    wrapper_f = LLVM.Function(mod, wrapper_fn, wrapper_ft)

    # emit IR performing the "conversions"
    let builder = Builder(ctx)
        entry = BasicBlock(wrapper_f, "entry", ctx)
        position!(builder, entry)

        wrapper_args = Vector{LLVM.Value}()

        # perform argument conversions
        codegen_types = parameters(entry_ft)
        wrapper_params = parameters(wrapper_f)
        param_index = 0
        for (julia_t, codegen_t, wrapper_t, wrapper_param) in
            zip(julia_types, codegen_types, wrapper_types, wrapper_params)
            param_index += 1
            if codegen_t != wrapper_t
                # the wrapper argument doesn't match the kernel parameter type.
                # this only happens when codegen wants to pass a pointer.
                @compiler_assert isa(codegen_t, LLVM.PointerType) job
                @compiler_assert eltype(codegen_t) == wrapper_t job

                # copy the argument value to a stack slot, and reference it.
                ptr = alloca!(builder, wrapper_t)
                if LLVM.addrspace(codegen_t) != 0
                    ptr = addrspacecast!(builder, ptr, codegen_t)
                end
                store!(builder, wrapper_param, ptr)
                push!(wrapper_args, ptr)
            else
                push!(wrapper_args, wrapper_param)
                for attr in collect(parameter_attributes(entry_f, param_index))
                    push!(parameter_attributes(wrapper_f, param_index), attr)
                end
            end
        end

        call!(builder, entry_f, wrapper_args)

        ret!(builder)

        dispose(builder)
    end

    # early-inline the original entry function into the wrapper
    push!(function_attributes(entry_f), EnumAttribute("alwaysinline", 0, ctx))
    linkage!(entry_f, LLVM.API.LLVMInternalLinkage)

    fixup_metadata!(entry_f)
    ModulePassManager() do pm
        always_inliner!(pm)
        run!(pm, mod)
    end

    return wrapper_f
end
# HACK: get rid of invariant.load and const TBAA metadata on loads from pointer args,
#       since storing to a stack slot violates the semantics of those attributes.
# TODO: can we emit a wrapper that doesn't violate Julia's metadata?
function fixup_metadata!(f::LLVM.Function)
    for param in parameters(f)
        if isa(llvmtype(param), LLVM.PointerType)
            # collect all uses of the pointer
            worklist = Vector{LLVM.Instruction}(user.(collect(uses(param))))
            while !isempty(worklist)
                value = popfirst!(worklist)

                # remove the invariant.load attribute
                md = metadata(value)
                if haskey(md, LLVM.MD_invariant_load)
                    delete!(md, LLVM.MD_invariant_load)
                end
                if haskey(md, LLVM.MD_tbaa)
                    delete!(md, LLVM.MD_tbaa)
                end

                # recurse on the output of some instructions
                if isa(value, LLVM.BitCastInst) ||
                   isa(value, LLVM.GetElementPtrInst) ||
                   isa(value, LLVM.AddrSpaceCastInst)
                    append!(worklist, user.(collect(uses(value))))
                end

                # IMPORTANT NOTE: if we ever want to inline functions at the LLVM level,
                # we need to recurse into call instructions here, and strip metadata from
                # called functions (see CUDAnative.jl#238).
            end
        end
    end
end


text:
#unique!(push!(LOAD_PATH, "~/Documents/UGent/PhD/Code_cMPS/julia1.0/CMPSKit.jl/src/"))
unique!(push!(LOAD_PATH, joinpath(pwd(), "src")))

using Revise
using CMPSKit
using KrylovKit
using OptimKit
using LinearAlgebra
using JLD2
using TensorOperations
using Plots

D = 4
k = 1.
μ = 2.3
c = 1.5
g = 2.

Id = 1*Matrix(I,D,D)
KL = Constant(randn(D^2,D^2))
KL = 0.5*(KL-KL')
R1 = Constant(kron(Id,randn(D,D)))
R2 = Constant(kron(randn(D,D),Id))
RLs = (R1,R2)
QL = KL
for R in RLs
    mul!(QL, R', R, -1/2, 1)
end
#Put them in cMPS form
Ψ = InfiniteCMPS(QL, (R1,R2); gauge = :left)

#h = k * (∂ψ[1]'*∂ψ[1] + ∂ψ[2]'*∂ψ[2]) - μ * (ψ[1]'*ψ[1] + ψ[2]'*ψ[2]) + c * ((ψ[1]')^2*ψ[1]^2 + (ψ[2]')^2*ψ[2]^2) + g * (ψ[1]'*ψ[2]'*ψ[1]*ψ[2])
#H = ∫(h, (-Inf,+Inf))

h = k * (∂ψ[1]'*∂ψ[1] + ∂ψ[2]'*∂ψ[2]) - μ * (ψ[1]'*ψ[1] + ψ[2]'*ψ[2]) + c * ((ψ[1]')^2*ψ[1]^2 + (ψ[2]')^2*ψ[2]^2)
H = ∫(h, (-Inf,+Inf))

alg1 = LBFGS(; verbosity = 2, maxiter = 1000000, gradtol = 1e-4);

#Ψ, ρR, E, e, normgrad, numfg, history = groundstate4_constrained(H, Ψ, (0.63,0.63); optalg = alg1, linalg = linalg)
αs,fs, dfs1, dfs2 = groundstate4_constrained(H, Ψ, (0.63,0.63); optalg = alg1, linalg = GMRES(; tol = 1e-5))
αs = (αs[1:end-1] + αs[2:end])/2
push!(αs,0.1)
display(plot(αs,[dfs1,dfs2]))
gui()


text:
### Semi-infinite models ###

function _get_fit_model!(Rt::Array{Float64}, β::Array{Float64,1}, data::fitDTOF, mod::Val{refl_DA_semiinf_TD})
    Rt = refl_DA_semiinf_TD(data.t, β, data.ndet, data.nmed)
end

# fit DTOF to fluence with EBPC conditions
function _get_fit_model!(Rt::Array{Float64}, β::Array{Float64,1}, data::fitDTOF, mod::Val{fluence_DA_semiinf_TD})
    Rt = fluence_DA_semiinf_TD(data.t, β, data.ndet, data.nmed, 0.0)
    return Rt./(2*get_afac(data.nmed/data.ndet))
end

### Slab models ###

#=
function _get_fit_model!(Rt::Array{Float64}, β::Array{Float64,1}, data::fitDTOF, mod::Val{refl_DA_slab_TD})
    Rt = refl_DA_slab_TD(data.t, β, data.ndet, data.nmed, data.s)
end

### Parallelepiped models ###

function _get_fit_model!(Rt::Array{Float64}, β::Array{Float64,1}, data::fitDTOF, mod::Val{refl_DA_paralpip_TD})
    Rt = refl_DA_paralpip_TD(data.t, β, data.ndet, data.nmed, data.rd, data.rs, data.L)
end
=#

max_issues_repo_name: heltonmc/LightPropagation


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 9daa72fd-1145-4733-a21e-bd209367380a
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ 947fbdc9-3c48-44f3-b4d3-183f70e2194e
using LinearAlgebra

# ╔═╡ a465636e-bc12-4e8f-9db2-1a20ea771504
begin
	𝐱 = rand(100, 3)
	𝐱 = hcat(ones(size(𝐱, 1)), 𝐱)
	y = 𝐱 * [1, 2, 3, 4] + 1e-3randn(size(𝐱, 1))
end;

# ╔═╡ 833c56c2-d2ab-4f94-968b-b5771c9fa945
function irls(𝐱, y; p=2, n_iter=100)
	w = ones(size(𝐱, 1))
	𝐰 = diagm(w)
	β = inv(𝐱' * 𝐰 * 𝐱) * 𝐱' * 𝐰 * y

	if p == 2
		return β
	end

	for _ in 1:n_iter
		w = abs.(y - 𝐱 * β).^(p-2)
		𝐰 = diagm(w)
		β = inv(𝐱' * 𝐰 * 𝐱) * 𝐱' * 𝐰 * y
	end

	return β
end

# ╔═╡ 2c663777-b6e3-42b2-a279-d7b42dd27d70
irls(𝐱, y, p=4, n_iter=10)

# ╔═╡ Cell order:
# ╟─9daa72fd-1145-4733-a21e-bd209367380a
# ╠═947fbdc9-3c48-44f3-b4d3-183f70e2194e
# ╠═a465636e-bc12-4e8f-9db2-1a20ea771504
# ╠═833c56c2-d2ab-4f94-968b-b5771c9fa945
# ╠═2c663777-b6e3-42b2-a279-d7b42dd27d70


text:
using ITensors

function ctmrg(T::ITensor, Cₗᵤ::ITensor, Aₗ::ITensor; χmax::Int, cutoff=0.0, nsteps::Int)
  sₕ = commonind(T, Aₗ)
  sᵥ = uniqueind(T, Aₗ, Aₗ'; plev=0)
  lᵥ = commonind(Cₗᵤ, Aₗ)
  lₕ = uniqueind(Cₗᵤ, Aₗ)
  Aᵤ = replaceinds(Aₗ, lᵥ => lₕ, lᵥ' => lₕ', sₕ => sᵥ)
  Cₗᵤ = dense(Cₗᵤ)
  for i in 1:nsteps
    ## Get the grown corner transfer matrix (CTM)
    Cₗᵤ⁽¹⁾ = Aₗ * Cₗᵤ * Aᵤ * T

    ## Diagonalize the grown CTM
    # TODO: replace with
    # eigen(Cₗᵤ⁽¹⁾, "horiz" => "vert"; tags = "horiz" => "vert", kwargs...)
    Cₗᵤ, Uᵥ = eigen(
      Cₗᵤ⁽¹⁾,
      (lₕ', sₕ'),
      (lᵥ', sᵥ');
      ishermitian=true,
      cutoff,
      maxdim=χmax,
      lefttags=tags(lₕ),
      righttags=tags(lᵥ),
    )
    Cₗᵤ = dense(Cₗᵤ)
    lᵥ = commonind(Cₗᵤ, Uᵥ)
    lₕ = uniqueind(Cₗᵤ, Uᵥ)

    # The renormalized CTM is the diagonal matrix of eigenvalues
    # Normalize the CTM
    Cₗ = Cₗᵤ * prime(dag(Cₗᵤ), lₕ)
    normC = (Cₗ * dag(Cₗ))[]^(1 / 4)
    Cₗᵤ = Cₗᵤ / normC

    # Calculate the renormalized half row transfer matrix (HRTM)
    Uᵥ = noprime(Uᵥ)
    Aₗ = Aₗ * Uᵥ * T * dag(Uᵥ')
    Aₗ = replaceinds(Aₗ, sₕ' => sₕ)

    # Normalize the HRTM
    ACₗ = Aₗ * Cₗᵤ * prime(dag(Cₗᵤ))
    normA = √((ACₗ * dag(ACₗ))[])
    Aₗ = Aₗ / normA
    Aᵤ = replaceinds(Aₗ, lᵥ => lₕ, lᵥ' => lₕ', sₕ => sᵥ)
  end
  return Cₗᵤ, Aₗ
end


text:
@testset "get_growingseason" begin
    rng = StableRNG(815)
    GPPd = @pipe sin.((1:365).*π./365) .+ 0.5 .* rand(rng,365) |> allowmissing(_)
    GPPd[100:120] .= missing
    #plot(1:365, GPPd)
    growseas = @test_logs (:warn, r"gap in 'GPPd'") get_growingseason(GPPd, 0.5)
    #plot(growseas)
    @test rle(growseas.is_growingseason) == (Bool[0, 1, 0], [51, 258, 56])
    growseas = get_growingseason(GPPd, 0.5; min_int = 56, warngap = false)
    @test rle(growseas.is_growingseason) == (Bool[1, 0], [309, 56])
end

@testset "get_growingseason with few data" begin
    rng = StableRNG(815)
    GPPd = @pipe sin.((1:365).*π./365) .+ 0.5 .* rand(rng,365) |> allowmissing(_)
    GPPd[1:200] .= missing
    @test_throws ErrorException get_growingseason(GPPd, 0.5)
end

@testset "setinvalid_nongrowingseason!" begin
    rng = StableRNG(815)
    nday = 365
    GPPd = @pipe sin.((1:nday).*π./365) .+ 0.5 .* rand(rng,365) |> allowmissing(_)
    GPPd[100:120] .= missing
    df = DataFrame(datetime = DateTime(2021) .+ Hour(1).+ Day.(0:nday-1), GPP = GPPd)
    df2 = copy(df)
    df2a = setinvalid_nongrowingseason!(df2, 0.5; warngap = false)
    @test df2a === df2
    @test propertynames(df2a) == union(propertynames(df), SA[:valid])
    @test rle(df2a.valid) == (Bool[0, 1, 0], [51, 258, 56])
    #
    df2.valid[80:90] .= false
    df2a = setinvalid_nongrowingseason!(df2, 0.5; warngap = false, update_GPPd_smoothed = true)
    @test df2a === df2
    @test propertynames(df2a) == union(propertynames(df), SA[:valid, :GPPd_smoothed])
    @test rle(df2a.valid) == (Bool[0, 1, 0, 1, 0], [51, 28, 11, 219, 56])
end

@testset "setinvalid_qualityflag!" begin
    df = DataFrame(
        NEE = 1:3,
        GPP = 10:10:30,
        NEE_qc = [1,1,2],
        GPP_qc = [1,missing,1],
    )
    df2 = copy(df)
    df2a = setinvalid_qualityflag!(df2; vars = SA["NEE", "GPP"], setvalmissing = false)
    @test df2a === df2
    @test df2a.valid == [true, false, false]
    @test df2a.NEE == 1:3 # not set to missing
    # test with existing :valid column, where already false entries are kept
    df2.valid[1] = false
    df2a = setinvalid_qualityflag!(df2; vars = SA["NEE", "GPP"], setvalmissing = false)
    @test df2a === df2
    @test df2a.valid == [false, false, false]
    @test df2a.NEE == 1:3 # not set to missing
    #
    # test with setting values to missing
    df2 = copy(df)
    df2a = setinvalid_qualityflag!(df2; vars = SA["NEE", "GPP"], setvalmissing = true)
    @test df2a === df2
    @test df2a.valid == [true, false, false]
    @test isequal(df2a.NEE, [1,2,missing])
    @test isequal(df2a.GPP, [10,missing,30])
    df2 = copy(df)
    df2a = setinvalid_qualityflag!(df2; vars = SA["NEE", "GPP"], setvalmissing = true, missing_qc_as_bad = false)
    @test df2a === df2
    @test df2a.valid == [true, true, false]
    @test isequal(df2a.NEE, [1,2,missing])
    @test isequal(df2a.GPP, [10,20,30])
end

@testset "setinvalid_range!" begin
    df = DataFrame(
        NEE = 1:3,
        GPP = 10:10:30,
    )
    allowmissing!(df, :NEE)
    var_ranges = [:NEE => (-2.0,4.0), :GPP => (8.0,28.0)]
    df2 = copy(df)
    df3 = copy(df2)
    df2a = setinvalid_range!(df2, var_ranges...; setvalmissing = false)
    # @btime setinvalid_range!($df2, $(var_ranges)...; setvalmissing = false)
    @test df2a === df2
    @test df2a.valid == [true, true, false]
    df3a = setinvalid_range!(df3, var_ranges...)
    @test df3a === df3
    @test isequal(df3a.valid, df2a.valid)
    @test ismissing(df3a.GPP[3])
    @test all(.!ismissing.(df3a.GPP[1:2]))
    #
    df2.NEE[1] = missing
    df3 = copy(df2)
    df2b = setinvalid_range!(df2, var_ranges...; setvalmissing = false)
    @test df2b === df2
    @test df2b.valid == [false, true, false]
    df3a = setinvalid_range!(df3, var_ranges...)
    @test isequal(df3a.valid, df2a.valid)
    @test ismissing(df3a.GPP[3])
    @test all(.!ismissing.(df3a.GPP[1:2]))
    @test all(.!ismissing.(df3a.NEE[2:3]))
    #
    df2.valid[2] = false
    df3 = copy(df2)
    df2b = setinvalid_range!(df2, var_ranges...; setvalmissing = false)
    @test df2b === df2
    @test df2b.valid == [false, false, false]
    df3a = setinvalid_range!(df3, var_ranges...)
    @test isequal(df3a.valid, df2a.valid)
end

@testset "setinvalid_afterprecip!" begin
    min_precip = 0.01
    hours_after = 1.0
    dfo = DataFrame(datetime = DateTime(2021) .+ Hour.(1:13), precip = 0.0)
    dfo.precip[3] = 0.001
    dfo.precip[10:12] .= 2.0
    dfo.precip[6:8] .= 2.0
    allowmissing!(dfo, :precip); dfo.precip[2] = missing
    #
    df = copy(dfo)
    df2 = setinvalid_afterprecip!(df; min_precip, hours_after)
    @test df2 === df
    @test all(.!df.valid[6:9])
    @test all(.!df.valid[10:13])
    @test all(df.valid[Not(vcat(6:9,10:13))])
end


text:
using LinearAlgebra
include("findlocalmaxima.jl")

function MUSIC(Y,K,m::Int=0,f_steps=1000)
    # Estimation of signal parameters via rotational invariance techniques
    # (ESPRIT) algorithm for line spectral estimation from measurements.
    # Can handle multiple snapshots by constructing a low rank Hankel matrix
    # from the data.
    #
    # Inputs:
    # Y: Data matrix
    # K: number of sinusoids
    # m: Hankel parameter of the Hankel matrix of the data
    # f_steps: number of data points of the MUSIC pseudo spectrum
    #
    # Parameters:
    # N = size(Y,1): number of measurements per snapshot
    # L = size(Y,2): number of snapshots
    #
    # Output:
    # f_hat: estimate of the frequencies of the K sinusoids in [rad]

(N,L) = size(Y)
f_range = range(0, 1, length=f_steps)

if m == 0
m = ceil(Int,N/2)
end

C = zeros(m,m)
for j = 1:L, i = m:N
    C = C + Y[i-m+1:i,j]*Y[i-m+1:i,j]'/(L*N)
end

    evdC = eigen(C)
    p = sortperm(real(evdC.values))[1:m-K]
    Un = evdC.vectors[:,p]

    a(f) = exp.(2*pi*1im*f*(0:m-1))
    spec(f) = 1/norm(Un'*a(f))^2
    Mspec = spec.(f_range)/maximum(spec.(f_range))
    localmaxima = findlocalmaxima(Mspec)
    q = localmaxima[sortperm(Mspec[localmaxima],rev=true)[1:min(K,length(localmaxima))]]
    f_hat = f_range[q]
    return f_hat #, Mspec

end


text:
using GeoLearning
using Meshes
using GeoStatsBase
using CategoricalArrays
using Test, Random, Plots
using ReferenceTests, ImageIO

# learning models from MLJ
using MLJ: @load
dtree = @load DecisionTreeClassifier pkg=DecisionTree verbosity=0

# workaround GR warnings
ENV["GKSwstype"] = "100"

# environment settings
isCI = "CI" ∈ keys(ENV)
islinux = Sys.islinux()
visualtests = !isCI || (isCI && islinux)
datadir = joinpath(@__DIR__,"data")

# list of tests
testfiles = [
  "pointwise.jl"
]

@testset "GeoLearning.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
module TestUnit

using Test
using SimpleMock
using Documenter
using Trixi

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

# Run various unit (= non-elixir-triggered) tests
@testset "Unit tests" begin
  @testset "SerialTree" begin
    @testset "constructors" begin
      @test_nowarn Trixi.SerialTree(Val(1), 10, 0.0, 1.0)
    end

    @testset "helper functions" begin
      t = Trixi.SerialTree(Val(1), 10, 0.0, 1.0)
      @test_nowarn display(t)
      @test Trixi.ndims(t) == 1
      @test Trixi.has_any_neighbor(t, 1, 1) == true
      @test Trixi.isperiodic(t, 1) == true
      @test Trixi.n_children_per_cell(t) == 2
      @test Trixi.n_children_per_cell(2) == 4
      @test Trixi.n_directions(t) == 2
    end

    @testset "refine!/coarsen!" begin
      t = Trixi.SerialTree(Val(1), 10, 0.0, 1.0)
      @test Trixi.refine!(t) == [1]
      @test Trixi.coarsen!(t) == [1]
      @test Trixi.refine!(t) == [1]
      @test Trixi.coarsen!(t, 1) == [1]
      @test Trixi.coarsen!(t) == Int[] # Coarsen twice to check degenerate case of single-cell tree
      @test Trixi.refine!(t) == [1]
      @test Trixi.refine!(t) == [2,3]
      @test Trixi.coarsen_box!(t, [-0.5], [0.0]) == [2]
      @test Trixi.coarsen_box!(t, 0.0, 0.5) == [3]
      @test isnothing(Trixi.reset_data_structures!(t))
    end
  end

  @testset "ParallelTree" begin
    @testset "constructors" begin
      @test_nowarn Trixi.ParallelTree(Val(1), 10, 0.0, 1.0)
    end

    @testset "helper functions" begin
      t = Trixi.ParallelTree(Val(1), 10, 0.0, 1.0)
      @test isnothing(display(t))
      @test isnothing(Trixi.reset_data_structures!(t))
    end
  end

  @testset "TreeMesh" begin
    @testset "constructors" begin
      @test TreeMesh{1, Trixi.SerialTree{1}}(1, 5.0, 2.0) isa TreeMesh
    end
  end

  @testset "ParallelTreeMesh" begin
    @testset "partition!" begin
      @testset "mpi_nranks() = 2" begin
        mock((Trixi.mpi_nranks) => () -> 2) do _
          mesh = TreeMesh{2, Trixi.ParallelTree{2}}(30, (0.0, 0.0), 1)
          # Refine twice
          Trixi.refine!(mesh.tree)
          Trixi.refine!(mesh.tree)

          # allow_coarsening = true
          Trixi.partition!(mesh)
          # Use parent for OffsetArray
          @test parent(mesh.n_cells_by_rank) == [11, 10]
          @test mesh.tree.mpi_ranks[1:21] ==
              [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
          @test parent(mesh.first_cell_by_rank) == [1, 12]

          # allow_coarsening = false
          Trixi.partition!(mesh; allow_coarsening=false)
          @test parent(mesh.n_cells_by_rank) == [11, 10]
          @test mesh.tree.mpi_ranks[1:21] ==
              [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
          @test parent(mesh.first_cell_by_rank) == [1, 12]
        end
      end

      @testset "mpi_nranks() = 3" begin
        mock((Trixi.mpi_nranks) => () -> 3) do _
          mesh = TreeMesh{2, Trixi.ParallelTree{2}}(100, (0.0, 0.0), 1)
          # Refine twice
          Trixi.refine!(mesh.tree)
          Trixi.refine!(mesh.tree)

          # allow_coarsening = true
          Trixi.partition!(mesh)
          # Use parent for OffsetArray
          @test parent(mesh.n_cells_by_rank) == [11, 5, 5]
          @test mesh.tree.mpi_ranks[1:21] ==
              [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2]
          @test parent(mesh.first_cell_by_rank) == [1, 12, 17]

          # allow_coarsening = false
          Trixi.partition!(mesh; allow_coarsening=false)
          @test parent(mesh.n_cells_by_rank) == [9, 6, 6]
          @test mesh.tree.mpi_ranks[1:21] ==
              [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2]
          @test parent(mesh.first_cell_by_rank) == [1, 10, 16]
        end
      end

      @testset "mpi_nranks() = 9" begin
        mock((Trixi.mpi_nranks) => () -> 9) do _
          mesh = TreeMesh{2, Trixi.ParallelTree{2}}(1000, (0.0, 0.0), 1)
          # Refine twice
          Trixi.refine!(mesh.tree)
          Trixi.refine!(mesh.tree)
          Trixi.refine!(mesh.tree)
          Trixi.refine!(mesh.tree)

          # allow_coarsening = true
          Trixi.partition!(mesh)
          # Use parent for OffsetArray
          @test parent(mesh.n_cells_by_rank) == [44, 37, 38, 37, 37, 37, 38, 37, 36]
          @test parent(mesh.first_cell_by_rank) == [1, 45, 82, 120, 157, 194, 231, 269, 306]
        end
      end

      @testset "mpi_nranks() = 3 non-uniform" begin
        mock((Trixi.mpi_nranks) => () -> 3) do _
          mesh = TreeMesh{2, Trixi.ParallelTree{2}}(100, (0.0, 0.0), 1)
          # Refine whole tree
          Trixi.refine!(mesh.tree)
          # Refine left leaf
          Trixi.refine!(mesh.tree, [2])

          # allow_coarsening = true
          Trixi.partition!(mesh)
          # Use parent for OffsetArray
          @test parent(mesh.n_cells_by_rank) == [6, 1, 2]
          @test mesh.tree.mpi_ranks[1:9] == [0, 0, 0, 0, 0, 0, 1, 2, 2]
          @test parent(mesh.first_cell_by_rank) == [1, 7, 8]

          # allow_coarsening = false
          Trixi.partition!(mesh; allow_coarsening=false)
          @test parent(mesh.n_cells_by_rank) == [5, 2, 2]
          @test mesh.tree.mpi_ranks[1:9] == [0, 0, 0, 0, 0, 1, 1, 2, 2]
          @test parent(mesh.first_cell_by_rank) == [1, 6, 8]
        end
      end

      @testset "not enough ranks" begin
        mock((Trixi.mpi_nranks) => () -> 3) do _
          mesh = TreeMesh{2, Trixi.ParallelTree{2}}(100, (0.0, 0.0), 1)

          # Only one leaf
          @test_throws AssertionError(
            "Too many ranks to properly partition the mesh!") Trixi.partition!(mesh)

          # Refine to 4 leaves
          Trixi.refine!(mesh.tree)

          # All four leaves will need to be on one rank to allow coarsening
          @test_throws AssertionError(
            "Too many ranks to properly partition the mesh!") Trixi.partition!(mesh)
          @test_nowarn Trixi.partition!(mesh; allow_coarsening=false)
        end
      end
    end
  end

  @testset "interpolation" begin
    @testset "nodes and weights" begin
      @test Trixi.gauss_nodes_weights(1) == ([0.0], [2.0])
    end

    @testset "multiply_dimensionwise" begin
      nodes_in  = [0.0, 0.5, 1.0]
      nodes_out = [0.0, 1/3, 2/3, 1.0]
      matrix = Trixi.polynomial_interpolation_matrix(nodes_in, nodes_out)
      data_in = [3.0 4.5 6.0]
      @test isapprox(Trixi.multiply_dimensionwise(matrix, data_in), [3.0 4.0 5.0 6.0])

      n_vars   = 3
      size_in  = 4
      size_out = 5
      matrix   = randn(size_out, size_in)
      # 1D
      data_in  = randn(n_vars, size_in)
      data_out = Trixi.multiply_dimensionwise_naive(matrix, data_in)
      @test isapprox(data_out, Trixi.multiply_dimensionwise(matrix, data_in))
      # 2D
      data_in  = randn(n_vars, size_in, size_in)
      data_out = Trixi.multiply_dimensionwise_naive(matrix, data_in)
      @test isapprox(data_out, Trixi.multiply_dimensionwise(matrix, data_in))
      # 3D
      data_in  = randn(n_vars, size_in, size_in, size_in)
      data_out = Trixi.multiply_dimensionwise_naive(matrix, data_in)
      @test isapprox(data_out, Trixi.multiply_dimensionwise(matrix, data_in))
    end
  end

  @testset "L2 projection" begin
    @testset "calc_reverse_upper for LGL" begin
      @test isapprox(Trixi.calc_reverse_upper(2, Val(:gauss_lobatto)), [[0.25, 0.25] [0.0, 0.5]])
    end
    @testset "calc_reverse_lower for LGL" begin
      @test isapprox(Trixi.calc_reverse_lower(2, Val(:gauss_lobatto)), [[0.5, 0.0] [0.25, 0.25]])
    end
  end

  @testset "containers" begin
    # Set up mock container
    mutable struct MyContainer <: Trixi.AbstractContainer
      data::Vector{Int}
      capacity::Int
      length::Int
      dummy::Int
    end
    function MyContainer(data, capacity)
      c = MyContainer(Vector{Int}(undef, capacity+1), capacity, length(data), capacity+1)
      c.data[1:length(data)] .= data
      return c
    end
    MyContainer(data::AbstractArray) = MyContainer(data, length(data))
    Trixi.invalidate!(c::MyContainer, first, last) = (c.data[first:last] .= 0; c)
    function Trixi.raw_copy!(target::MyContainer, source::MyContainer, first, last, destination)
      Trixi.copy_data!(target.data, source.data, first, last, destination)
      return target
    end
    Trixi.move_connectivity!(c::MyContainer, first, last, destination) = c
    Trixi.delete_connectivity!(c::MyContainer, first, last) = c
    Trixi.reset_data_structures!(c::MyContainer) = (c.data = Vector{Int}(undef, c.capacity+1); c)
    function Base.:(==)(c1::MyContainer, c2::MyContainer)
      return (c1.capacity == c2.capacity &&
              c1.length == c2.length &&
              c1.dummy == c2.dummy &&
              c1.data[1:c1.length] == c2.data[1:c2.length])
    end

    @testset "size" begin
      c = MyContainer([1, 2, 3])
      @test size(c) == (3,)
    end

    @testset "resize!" begin
      c = MyContainer([1, 2, 3])
      @test length(resize!(c, 2)) == 2
    end

    @testset "copy!" begin
      c1 = MyContainer([1, 2, 3])
      c2 = MyContainer([4, 5])
      @test Trixi.copy!(c1, c2, 2, 1, 2) == MyContainer([1, 2, 3]) # no-op

      c1 = MyContainer([1, 2, 3])
      c2 = MyContainer([4, 5])
      @test Trixi.copy!(c1, c2, 1, 2, 2) == MyContainer([1, 4, 5])

      c1 = MyContainer([1, 2, 3])
      @test Trixi.copy!(c1, c2, 1, 2) == MyContainer([1, 4, 3])

      c1 = MyContainer([1, 2, 3])
      @test Trixi.copy!(c1, 2, 3, 1) == MyContainer([2, 3, 3])

      c1 = MyContainer([1, 2, 3])
      @test Trixi.copy!(c1, 1, 3) == MyContainer([1, 2, 1])
    end

    @testset "move!" begin
      c = MyContainer([1, 2, 3])
      @test Trixi.move!(c, 1, 1) == MyContainer([1, 2, 3]) # no-op

      c = MyContainer([1, 2, 3])
      @test Trixi.move!(c, 1, 2) == MyContainer([0, 1, 3])
    end

    @testset "swap!" begin
      c = MyContainer([1,2])
      @test Trixi.swap!(c, 1, 1) == MyContainer([1, 2]) # no-op

      c = MyContainer([1,2])
      @test Trixi.swap!(c, 1, 2) == MyContainer([2,1])
    end

    @testset "erase!" begin
      c = MyContainer([1, 2])
      @test Trixi.erase!(c, 2, 1) == MyContainer([1, 2]) # no-op

      c = MyContainer([1, 2])
      @test Trixi.erase!(c, 1) == MyContainer([0, 2])
    end

    @testset "remove_shift!" begin
      c = MyContainer([1, 2, 3, 4])
      @test Trixi.remove_shift!(c, 2, 1) == MyContainer([1, 2, 3, 4]) # no-op

      c = MyContainer([1, 2, 3, 4])
      @test Trixi.remove_shift!(c, 2, 2) == MyContainer([1, 3, 4], 4)

      c = MyContainer([1, 2, 3, 4])
      @test Trixi.remove_shift!(c, 2) == MyContainer([1, 3, 4], 4)
    end

    @testset "remove_fill!" begin
      c = MyContainer([1, 2, 3, 4])
      @test Trixi.remove_fill!(c, 2, 1) == MyContainer([1, 2, 3, 4]) # no-op

      c = MyContainer([1, 2, 3, 4])
      @test Trixi.remove_fill!(c, 2, 2) == MyContainer([1, 4, 3], 4)
    end

    @testset "reset!" begin
      c = MyContainer([1, 2, 3])
      @test Trixi.reset!(c, 2) == MyContainer(Int[], 2)
    end
  end

  @testset "example elixirs" begin
    @test basename(examples_dir()) == "examples"
    @test !isempty(get_examples())
    @test endswith(default_example(), "elixir_advection_basic.jl")
  end

  @testset "DG L2 mortar container debug output" begin
    c2d = Trixi.L2MortarContainer2D{Float64, 1, 1}(1)
    @test isnothing(display(c2d))
    c3d = Trixi.L2MortarContainer3D{Float64, 1, 1}(1)
    @test isnothing(display(c3d))
  end

  @testset "Printing indicators/controllers" begin
    # OBS! Constructing indicators/controllers using the parameters below doesn't make sense. It's
    # just useful to run basic tests of `show` methods.

    c = ControllerThreeLevelCombined(1, 2, 3, 10.0, 11.0, 12.0, "primary", "secondary", "cache")
    @test_nowarn show(stdout, c)

    indicator_hg = IndicatorHennemannGassner(1.0, 0.0, true, "variable", "cache")
    @test_nowarn show(stdout, indicator_hg)

    indicator_loehner = IndicatorLöhner(1.0, "variable", (; cache=nothing))
    @test_nowarn show(stdout, indicator_loehner)

    indicator_max = IndicatorMax("variable", (; cache=nothing))
    @test_nowarn show(stdout, indicator_max)
  end

  @testset "LBM 2D constructor" begin
    # Neither Mach number nor velocity set
    @test_throws ErrorException LatticeBoltzmannEquations2D(Ma=nothing, Re=1000)
    # Both Mach number and velocity set
    @test_throws ErrorException LatticeBoltzmannEquations2D(Ma=0.1, Re=1000, u0=1)
    # Neither Reynolds number nor viscosity set
    @test_throws ErrorException LatticeBoltzmannEquations2D(Ma=0.1, Re=nothing)
    # Both Reynolds number and viscosity set
    @test_throws ErrorException LatticeBoltzmannEquations2D(Ma=0.1, Re=1000, nu=1)

    # No non-dimensional values set
    @test LatticeBoltzmannEquations2D(Ma=nothing, Re=nothing, u0=1, nu=1) isa LatticeBoltzmannEquations2D
  end

  @testset "LBM 3D constructor" begin
    # Neither Mach number nor velocity set
    @test_throws ErrorException LatticeBoltzmannEquations3D(Ma=nothing, Re=1000)
    # Both Mach number and velocity set
    @test_throws ErrorException LatticeBoltzmannEquations3D(Ma=0.1, Re=1000, u0=1)
    # Neither Reynolds number nor viscosity set
    @test_throws ErrorException LatticeBoltzmannEquations3D(Ma=0.1, Re=nothing)
    # Both Reynolds number and viscosity set
    @test_throws ErrorException LatticeBoltzmannEquations3D(Ma=0.1, Re=1000, nu=1)

    # No non-dimensional values set
    @test LatticeBoltzmannEquations3D(Ma=nothing, Re=nothing, u0=1, nu=1) isa LatticeBoltzmannEquations3D
  end

  @testset "LBM 2D functions" begin
    # Set up LBM struct and dummy distribution
    equation = LatticeBoltzmannEquations2D(Ma=0.1, Re=1000)
    u = Trixi.equilibrium_distribution(1, 2, 3, equation)

    # Component-wise velocity
    @test isapprox(Trixi.velocity(u, 1, equation), 2)
    @test isapprox(Trixi.velocity(u, 2, equation), 3)
  end

  @testset "LBM 3D functions" begin
    # Set up LBM struct and dummy distribution
    equation = LatticeBoltzmannEquations3D(Ma=0.1, Re=1000)
    u = Trixi.equilibrium_distribution(1, 2, 3, 4, equation)

    # Component-wise velocity
    @test isapprox(velocity(u, 1, equation), 2)
    @test isapprox(velocity(u, 2, equation), 3)
    @test isapprox(velocity(u, 3, equation), 4)
  end

  @testset "LBMCollisionCallback" begin
    # Printing of LBM collision callback
    callback = LBMCollisionCallback()
    @test_nowarn show(stdout, callback)
    println()
    @test_nowarn show(stdout, "text/plain", callback)
    println()
  end

  # Test docstrings
  DocMeta.setdocmeta!(Trixi, :DocTestSetup, :(using Trixi); recursive=true)
  doctest(Trixi, manual=false)
end

end #module


text:
#################################################################
# This file implements the grid world problem as an MDP.
# In the problem, the agent is tasked with navigating in a
# stochatic environemnt. For example, when the agent chooses
# to go right, it may not always go right, but may go up, down
# or left with some probability. The agent's goal is to reach the
# reward states. The states with a positive reward are terminal,
# while the states with a negative reward are not.
#################################################################

#################################################################
# States and Actions
#################################################################
# state of the agent in grid world
struct GridWorldState # this is not immutable because of how it is used in transition(), but maybe it should be
    x::Int64 # x position
    y::Int64 # y position
    done::Bool # entered the terminal reward state in previous step - there is only one terminal state
    GridWorldState(x,y,done) = new(x,y,done)
    GridWorldState() = new()
end

# simpler constructors
GridWorldState(x::Int64, y::Int64) = GridWorldState(x,y,false)
# for state comparison
function ==(s1::GridWorldState,s2::GridWorldState)
    if s1.done && s2.done
        return true
    elseif s1.done || s2.done
        return false
    else
        return posequal(s1, s2)
    end
end
# for hashing states in dictionaries in Monte Carlo Tree Search
posequal(s1::GridWorldState, s2::GridWorldState) = s1.x == s2.x && s1.y == s2.y
function hash(s::GridWorldState, h::UInt64 = zero(UInt64))
    if s.done
        return hash(s.done, h)
    else
        return hash(s.x, hash(s.y, h))
    end
end
Base.copy!(dest::GridWorldState, src::GridWorldState) = (dest.x=src.x; dest.y=src.y; dest.done=src.done; return dest)

# action taken by the agent indicates desired travel direction
const GridWorldAction = Symbol # deprecated - this is here so that other people's code won't break

#################################################################
# Grid World MDP
#################################################################
# the grid world mdp type
mutable struct LegacyGridWorld <: MDP{GridWorldState, Symbol}
    size_x::Int64 # x size of the grid
    size_y::Int64 # y size of the grid
    reward_states::Vector{GridWorldState} # the states in which agent recieves reward
    reward_values::Vector{Float64} # reward values for those states
    bounds_penalty::Float64 # penalty for bumping the wall (will be added to reward)
    tprob::Float64 # probability of transitioning to the desired state
    terminals::Set{GridWorldState}
    discount_factor::Float64 # disocunt factor
end
# we use key worded arguments so we can change any of the values we pass in
function LegacyGridWorld(sx::Int64, # size_x
                   sy::Int64; # size_y
                   rs::Vector{GridWorldState}=[GridWorldState(4,3), GridWorldState(4,6), GridWorldState(9,3), GridWorldState(8,8)],
                   rv::Vector{Float64}=[-10.,-5,10,3],
                   penalty::Float64=0.0, # penalty for trying to go out of bounds  (will be added to reward)
                   tp::Float64=0.7, # tprob
                   discount_factor::Float64=0.95,
                   terminals=Set{GridWorldState}([rs[i] for i in filter(i->rv[i]>0.0, 1:length(rs))]))
    return LegacyGridWorld(sx, sy, rs, rv, penalty, tp, Set{GridWorldState}(terminals), discount_factor)
end

LegacyGridWorld(;sx::Int64=10, sy::Int64=10, kwargs...) = LegacyGridWorld(sx, sy; kwargs...)

#################################################################
# State and Action Spaces
#################################################################
# This could probably be implemented more efficiently without vectors

function states(mdp::LegacyGridWorld)
    s = vec(collect(GridWorldState(x, y, false) for x in 1:mdp.size_x, y in 1:mdp.size_y))
    push!(s, GridWorldState(0, 0, true))
    return s
end

actions(mdp::LegacyGridWorld) = [:up, :down, :left, :right]

function reward(mdp::LegacyGridWorld, state::GridWorldState, action::Symbol)
    if state.done
        return 0.0
    end
    r = static_reward(mdp, state)
    if !inbounds(mdp, state, action)
        r += mdp.bounds_penalty
    end
  return r
end

"""
    static_reward(mdp::LegacyGridWorld, state::GridWorldState)

Return the reward for being in the state (the reward not including bumping)
"""
function static_reward(mdp::LegacyGridWorld, state::GridWorldState)
    r = 0.0
    n = length(mdp.reward_states)
    for i = 1:n
        if posequal(state, mdp.reward_states[i])
            r += mdp.reward_values[i]
        end
    end
    return r
end

#checking boundries- x,y --> points of current state
inbounds(mdp::LegacyGridWorld,x::Int64,y::Int64) = 1 <= x <= mdp.size_x && 1 <= y <= mdp.size_y
inbounds(mdp::LegacyGridWorld,state::GridWorldState) = inbounds(mdp, state.x, state.y)

"""
    inbounds(mdp::LegacyGridWorld, s::GridWorldState, a::Symbol)

Return false if `a` is trying to go out of bounds, true otherwise.
"""
function inbounds(mdp::LegacyGridWorld, s::GridWorldState, a::Symbol)
    xdir = s.x
    ydir = s.y
    if a == :right
        xdir += 1
    elseif a == :left
        xdir -= 1
    elseif a == :up
        ydir += 1
    else
        # @assert a == :down
        ydir -= 1
    end
    return inbounds(mdp, GridWorldState(xdir, ydir, s.done))
end

function fill_probability!(p::AbstractVector{Float64}, val::Float64, index::Int64)
    for i = 1:length(p)
        if i == index
            p[i] = val
        else
            p[i] = 0.0
        end
    end
end

function transition(mdp::LegacyGridWorld, state::GridWorldState, action::Symbol)

    a = action
    x = state.x
    y = state.y

    neighbors = MVector(
        GridWorldState(x+1, y, false), # right
        GridWorldState(x-1, y, false), # left
        GridWorldState(x, y-1, false), # down
        GridWorldState(x, y+1, false), # up
        GridWorldState(x, y, false)    # stay
    )

    probability = MVector{5, Float64}(undef)
    fill!(probability, 0.0)

    if state.done
        fill_probability!(probability, 1.0, 5)
        neighbors[5] = GridWorldState(x, y, true)
        return SparseCat(neighbors, probability)
    end

    reward_states = mdp.reward_states
    reward_values = mdp.reward_values
    n = length(reward_states)
    if state in mdp.terminals
        fill_probability!(probability, 1.0, 5)
        neighbors[5] = GridWorldState(x, y, true)
        return SparseCat(neighbors, probability)
    end

    # The following match the definition of neighbors
    # given above
    target_neighbor = 0
    if a == :right
        target_neighbor = 1
    elseif a == :left
        target_neighbor = 2
    elseif a == :down
        target_neighbor = 3
    elseif a == :up
        target_neighbor = 4
    end
    # @assert target_neighbor > 0

    if !inbounds(mdp, neighbors[target_neighbor])
        # If would transition out of bounds, stay in
        # same cell with probability 1
        fill_probability!(probability, 1.0, 5)
    else
        probability[target_neighbor] = mdp.tprob

        oob_count = 0 # number of out of bounds neighbors

        for i = 1:length(neighbors)
            if !inbounds(mdp, neighbors[i])
                oob_count += 1
                @assert probability[i] == 0.0
            end
        end

        new_probability = (1.0 - mdp.tprob)/(3-oob_count)

        for i = 1:4 # do not include neighbor 5
            if inbounds(mdp, neighbors[i]) && i != target_neighbor
                probability[i] = new_probability
            end
        end
    end

    return SparseCat(neighbors, probability)
end

function actionindex(mdp::LegacyGridWorld, a::Symbol)
    # lazy, replace with switches when they arrive
    if a == :up
        return 1
    elseif a == :down
        return 2
    elseif a == :left
        return 3
    elseif a == :right
        return 4
    else
        error("Invalid action symbol $a")
    end
end

function stateindex(mdp::LegacyGridWorld, s::GridWorldState)
    return s2i(mdp, s)
end

function s2i(mdp::LegacyGridWorld, state::GridWorldState)
    if state.done
        return mdp.size_x*mdp.size_y + 1
    else
        return LinearIndices((mdp.size_x, mdp.size_y))[state.x, state.y]
    end
end

#=
function i2s(mdp::GridWorld, i::Int)
end
=#

isterminal(mdp::LegacyGridWorld, s::GridWorldState) = s.done

discount(mdp::LegacyGridWorld) = mdp.discount_factor

convert_s(::Type{A}, s::GridWorldState, mdp::LegacyGridWorld) where A<:AbstractArray = Float64[s.x, s.y, s.done]
convert_s(::Type{GridWorldState}, s::AbstractArray, mdp::LegacyGridWorld) = GridWorldState(s[1], s[2], s[3])

function a2int(a::Symbol, mdp::LegacyGridWorld)
    if a == :up
        return 0
    elseif a == :down
        return 1
    elseif a == :left
        return 2
    elseif a == :right
        return 3
    else
        throw("Action $a is invalid")
    end
end

function int2a(a::Int, mdp::LegacyGridWorld)
    if a == 0
        return :up
    elseif a == 1
        return :down
    elseif a == 2
        return :left
    elseif a == 3
        return :right
    else
        throw("Action $a is invalid")
    end
end

convert_a(::Type{A}, a::Symbol, mdp::LegacyGridWorld) where A<:AbstractArray = [Float64(a2int(a, mdp))]
convert_a(::Type{Symbol}, a::A, mdp::LegacyGridWorld) where A<:AbstractArray = int2a(Int(a[1]), mdp)

initialstate(mdp::LegacyGridWorld, rng::AbstractRNG) = GridWorldState(rand(rng, 1:mdp.size_x), rand(rng, 1:mdp.size_y))

initialstate(mdp::LegacyGridWorld) = ImplicitDistribution(mdp) do m, rng
    GridWorldState(rand(rng, 1:m.size_x), rand(rng, 1:m.size_y))
end

# Visualization
#=
function colorval(val, brightness::Real = 1.0)
    val = convert(Vector{Float64}, val)
    x = 255 .- min.(255, 255 * (abs.(val) ./ 10.0) .^ brightness)
    r = 255 * ones(size(val))
    g = 255 * ones(size(val))
    b = 255 * ones(size(val))
    r[val .>= 0] = x[val .>= 0]
    b[val .>= 0] = x[val .>= 0]
    g[val .< 0] = x[val .< 0]
    b[val .< 0] = x[val .< 0]
    return (r, g, b)
end

function plot(g::GridWorld, f::Function)
    V = map(f, states(g))
    plot(g, V)
end

function plot(mdp::GridWorld, V::Vector, state=GridWorldState(0,0,true))
    o = IOBuffer()
    sqsize = 1.0
    twid = 0.05
    (r, g, b) = colorval(V)
    for s in states(mdp)
        if !s.done
            (xval, yval) = (s.x, mdp.size_y-s.y+1)
            i = state_index(mdp, s)
            yval = 10 - yval
            println(o, "\\definecolor{currentcolor}{RGB}{$(r[i]),$(g[i]),$(b[i])}")
            println(o, "\\fill[currentcolor] ($((xval-1) * sqsize),$((yval) * sqsize)) rectangle +($sqsize,$sqsize);")
            if s == state
                println(o, "\\fill[orange] ($((xval-1) * sqsize),$((yval) * sqsize)) rectangle +($sqsize,$sqsize);")
            end
            vs = @sprintf("%0.2f", V[i])
            println(o, "\\node[above right] at ($((xval-1) * sqsize), $((yval) * sqsize)) {\$$(vs)\$};")
        end
    end
    println(o, "\\draw[black] grid(10,10);")
    tikzDeleteIntermediate(false)
    TikzPicture(String(take!(o)), options="scale=1.25")
end

function plot(mdp::GridWorld, state=GridWorldState(0,0,true))
    plot(mdp, zeros(n_states(mdp)), state)
end

function plot(g::GridWorld, f::Function, policy::Policy, state=GridWorldState(0,0,true))
    V = map(f, states(g))
    plot(g, V, policy, state)
end

function plot(mdp::GridWorld, V::Vector, policy::Policy, state=GridWorldState(0,0,true))
    o = IOBuffer()
    sqsize = 1.0
    twid = 0.05
    (r, g, b) = colorval(V)
    for s in states(mdp)
        if !s.done
            (xval, yval) = (s.x, mdp.size_y-s.y+1)
            i = state_index(mdp, s)
            yval = 10 - yval
            println(o, "\\definecolor{currentcolor}{RGB}{$(r[i]),$(g[i]),$(b[i])}")
            println(o, "\\fill[currentcolor] ($((xval-1) * sqsize),$((yval) * sqsize)) rectangle +($sqsize,$sqsize);")
            if s == state
                println(o, "\\fill[orange] ($((xval-1) * sqsize),$((yval) * sqsize)) rectangle +($sqsize,$sqsize);")
            end
        end
    end
    println(o, "\\begin{scope}[fill=gray]")
    for s in states(mdp)
        if !s.done
            (xval, yval) = (s.x, mdp.size_y-s.y+1)
            i = state_index(mdp, s)
            yval = 10 - yval + 1
            c = [xval, yval] * sqsize .- sqsize / 2
            C = [c'; c'; c']'
            RightArrow = [0 0 sqsize/2; twid -twid 0]
            dir = action(policy, s)
            if dir == :left
                A = [-1 0; 0 -1] * RightArrow + C
                println(o, "\\fill ($(A[1]), $(A[2])) -- ($(A[3]), $(A[4])) -- ($(A[5]), $(A[6])) -- cycle;")
            end
            if dir == :right
                A = RightArrow + C
                println(o, "\\fill ($(A[1]), $(A[2])) -- ($(A[3]), $(A[4])) -- ($(A[5]), $(A[6])) -- cycle;")
            end
            if dir == :up
                A = [0 -1; 1 0] * RightArrow + C
                println(o, "\\fill ($(A[1]), $(A[2])) -- ($(A[3]), $(A[4])) -- ($(A[5]), $(A[6])) -- cycle;")
            end
            if dir == :down
                A = [0 1; -1 0] * RightArrow + C
                println(o, "\\fill ($(A[1]), $(A[2])) -- ($(A[3]), $(A[4])) -- ($(A[5]), $(A[6])) -- cycle;")
            end

            vs = @sprintf("%0.2f", V[i])
            println(o, "\\node[above right] at ($((xval-1) * sqsize), $((yval-1) * sqsize)) {\$$(vs)\$};")
        end
    end
    println(o, "\\end{scope}");
    println(o, "\\draw[black] grid(10,10);");
    TikzPicture(String(take!(o)), options="scale=1.25")
end
=#


text:
#Download aima-data
using LibGit2
using Pkg

println("Downloading AIMACore from GitHub")
Pkg.clone("https://github.com/sambitdash/AIMACore.jl")

!isdir("downloads") && mkdir("downloads")
cd("downloads")

println("Downloading aima-data from GitHub")
if !isdir("aima-data")
    repo_path = joinpath("aima-data")
    repo_url  = "https://github.com/aimacode/aima-data.git"
    LibGit2.clone(repo_url, repo_path)
else
    repo = LibGit2.GitRepo("aima-data")
    LibGit2.fetch(repo)
end
println("aima-data - downloaded")


text:
@enum Status begin
    unsolved_status
    optimal_status
    unbounded_status
    infeasible_status
    max_iters_status
    numerical_error_status
end

@with_kw mutable struct Solution{T <: AbstractFloat}
    status::Status = unsolved_status
    optval::T = NaN
    x::Vector{T}
    y::Vector{T}
    z::Vector{T}
    s::Vector{T}
end

Solution{T}(n::Int64, m::Int64, p::Int64) where T = Solution{T}(unsolved_status, NaN, fill(NaN, n),
                                                                fill(NaN, m), fill(NaN, p), fill(NaN, p))
Solution(prob::Problem{T}) where T = Solution{T}(get_prob_sizes(prob)...)

mutable struct Stats{T <: AbstractFloat}
    iters::Int64
    residual::Vector{T}
    best::Vector{Bool}
    σ::Vector{T}
    μ::Vector{T}
    κ::Vector{T}
    τ::Vector{T}
    α_aff::Vector{T}
    α_comb::Vector{T}
    ir_iters_cache::Vector{Int64}
    ir_iters_aff::Vector{Int64}
    ir_iters_comb::Vector{Int64}
    regularization::Vector{T}
end

function Stats(settings::Settings)
    max_iters = settings.max_iters
    max_iters_ir = settings.max_iters_ir

    iters = 0
    residual = fill(NaN, max_iters)
    best = fill(false, max_iters)
    σ = fill(NaN, max_iters)
    μ = fill(NaN, max_iters)
    κ = fill(NaN, max_iters)
    τ = fill(NaN, max_iters)
    α_aff = fill(NaN, max_iters)
    α_comb = fill(NaN, max_iters)
    ir_iters_cache = fill(0, max_iters)
    ir_iters_aff = fill(0, max_iters)
    ir_iters_comb = fill(0, max_iters)
    regularization = fill(NaN, max_iters)
    return Stats(iters, residual, best, σ, μ, κ, τ, α_aff, α_comb,
                 ir_iters_cache, ir_iters_aff, ir_iters_comb, regularization)
end

# TODO can this be shortened?
function truncate_stats!(stats::Stats)
    i = stats.iters
    stats.residual = stats.residual[1:i]
    stats.best = stats.best[1:i]
    stats.σ = stats.σ[1:i]
    stats.μ = stats.μ[1:i]
    stats.κ = stats.κ[1:i]
    stats.τ = stats.τ[1:i]
    stats.α_aff = stats.α_aff[1:i]
    stats.α_comb = stats.α_comb[1:i]
    stats.ir_iters_cache = stats.ir_iters_cache[1:i]
    stats.ir_iters_aff = stats.ir_iters_aff[1:i]
    stats.ir_iters_comb = stats.ir_iters_comb[1:i]
    stats.regularization = stats.regularization[1:i]
end

function register_aff_dir_stats!(stats::Stats, ir_iters::Int64)
    stats.ir_iters_aff[stats.iters] = ir_iters
end

function register_comb_dir_stats!(stats::Stats{T}, ir_iters::Int64, μ_hat::T) where T
    stats.ir_iters_comb[stats.iters] = ir_iters
    stats.μ[stats.iters] = μ_hat
end

function register_cache_stats!(stats::Stats{T}, ir_iters::Int64, reg::T) where T
    stats.ir_iters_cache[stats.iters] = ir_iters
    stats.regularization[stats.iters] = reg
end


text:
using Revise

using RepresentativeDaysFinders
using JuMP
#using GLPK
using Gurobi
import DataFrames
##################################################################################
# Specify location of config-file
##################################################################################
config_file = normpath(joinpath(@__DIR__, "scenarios", "investment_model_BE.yaml"))

##
using Cbc

dft = findRepresentativeDays(config_file, with_optimizer(Cbc.Optimizer; seconds = 20))
# findRepresentativeDays(config_file, with_optimizer(GLPK.Optimizer; presolve=true, msg_lev=GLPK.MSG_ALL, tm_lim=180*1000))
# Juno.@enter findRepresentativeDays(config_file, with_optimizer(GLPK.Optimizer; presolve=true, msg_lev=GLPK.MSG_ALL, tm_lim=180*1000))

#RepresentativeDaysFinders.create_plots(dft)

df_dv = DataFrames.DataFrame(
            periods     = sort([k for k in dft.periods]),
            weights     = [dft.w[k] for k in sort([k for k in keys(dft.w)])],
            used_days   = [dft.u[k] for k in sort([k for k in keys(dft.u)])])
#CSV.write(joinpath(result_dir, "decision_variables.csv"), df_dv, delim=';')

df_dv_s = deepcopy(df_dv[df_dv[:weights] .> 0, :])

function ParseRepresentativeDaysResults(df_dv_s::DataFrames.DataFrame)
    selected_representativeDays = zeros(size(df_dv_s,1),2)
    for row in 1:size(df_dv_s,1)
        day = parse(Int,df_dv_s[row,:][:periods][2:4])
        weight = df_dv_s[row,:][:weights]
        println([day weight])
        selected_representativeDays[row,:] = [day weight]
        println([day,weight])
    end
    return selected_representativeDays
end
selected_representativeDays = ParseRepresentativeDaysResults(df_dv_s)


text:
"""
    Info <: AbstractDict{Symbol,String}

`MPI.Info` objects store key-value pairs, and are typically used for passing optional
arguments to MPI functions.

# Usage

These will typically be hidden from user-facing APIs by splatting keywords, e.g.
```julia
function f(args...; kwargs...)
    info = Info(kwargs...)
    # pass `info` object to `ccall`
end
```

For manual usage, `Info` objects act like Julia `Dict` objects:
```julia
info = Info(init=true) # keyword argument is required
info[key] = value
x = info[key]
delete!(info, key)
```

If `init=false` is used in the costructor (the default), a "null" `Info` object will be
returned: no keys can be added to such an object.
"""
@mpi_handle Info <: AbstractDict{Symbol,String} MPI_Info

const INFO_NULL = _Info(MPI_INFO_NULL)

function Info(;init=false)
    info = Info(INFO_NULL.val)
    if init
        @mpichk ccall((:MPI_Info_create, libmpi), Cint, (Ptr{MPI_Info},), info)
        refcount_inc()
        finalizer(free, info)
    end
    return info
end

function free(info::Info)
    if info.val != INFO_NULL.val
        @mpichk ccall((:MPI_Info_free, libmpi), Cint, (Ptr{MPI_Info},), info)
        refcount_dec()
    end
    return nothing
end

function Base.setindex!(info::Info, value::AbstractString, key::Symbol)
    skey = String(key)
    @assert isascii(skey) && isascii(value) &&
        length(skey) <= MPI_MAX_INFO_KEY && length(value) <= MPI_MAX_INFO_VAL
    @mpichk ccall((:MPI_Info_set, libmpi), Cint,
          (MPI_Info, Cstring, Cstring), info, skey, value)
end

Base.setindex!(info::Info, value::Any, key::Symbol) = info[key] = infoval(value)

"""
    infoval(x)

Convert Julia object `x` to a string representation for storing in an [`Info`](@ref) object.

The MPI specification allows passing strings, Boolean values, integers, and lists.
"""
infoval(x::String) = x
infoval(x::Bool) = string(x) # "true"/"false"
infoval(x::Integer) = string(x) # decimal values
infoval(xs::AbstractVector) = join(map(infoval, xs), ',') # separated by commas
infoval(xs::Tuple) = join(map(infoval, xs), ',') # separated by commas

function Info(kvs::Pair...)
    info = Info(;init=true)
    for (k,v) in kvs
        info[k] = v
    end
    return info
end

function Base.getindex(info::Info, key::Symbol)
    skey = String(key)
    @assert isascii(skey) && length(skey) <= MPI_MAX_INFO_KEY
    valuelen = Ref{Cint}()
    flag = Ref{Cint}()
    # int MPI_Info_get_valuelen(MPI_Info info, const char *key, int *valuelen, int *flag)
    @mpichk ccall((:MPI_Info_get_valuelen, libmpi), Cint,
          (MPI_Info, Cstring, Ptr{Cint}, Ptr{Cint}),
          info, skey, valuelen, flag)

    if flag[] == 0
        throw(KeyError(key))
    end

    # According to the MPI standard:
    #   "`valuelen` should be one less than the amount of allocated
    #   space to allow for the null terminator."
    # But MS-MPI will insists on setting the `n`th character as NUL,
    # so we simply pad both to avoid problems.
    n = valuelen[]
    buffer = Vector{UInt8}(undef, n+2)
    # int MPI_Info_get(MPI_Info info, const char *key, int valuelen, char *value, int *flag)
    @mpichk ccall((:MPI_Info_get, libmpi), Cint,
          (MPI_Info, Cstring, Cint, Ptr{UInt8}, Ptr{Cint}),
          info, skey, n+1, buffer, flag)
    return String(resize!(buffer,n))
end

function Base.delete!(info::Info,key::Symbol)
    skey = String(key)
    @assert isascii(skey) && length(skey) <= MPI_MAX_INFO_KEY
    @mpichk ccall((:MPI_Info_delete, libmpi), Cint,
          (MPI_Info, Cstring), info, skey)
end

function Base.length(info::Info)
    if info.val == INFO_NULL.val
        return 0
    end
    nkeys = Ref{Cint}()
    @mpichk ccall((:MPI_Info_get_nkeys, libmpi), Cint,
                   (MPI_Info, Ptr{Cint}), info, nkeys)
    return Int(nkeys[])
end

function nthkey(info::Info, n::Integer)
    buffer = Vector{UInt8}(undef, MPI_MAX_INFO_KEY+1)
    @mpichk ccall((:MPI_Info_get_nthkey, libmpi), Cint,
                  (MPI_Info, Cint, Ptr{UInt8}), info, n, buffer)
    i = findfirst(isequal(UInt8(0)), buffer)
    if i !== nothing
        resize!(buffer, i-1)
    end
    Symbol(buffer)
end

function Base.iterate(infokeys::Base.KeySet{Symbol,Info}, i=0)
    info = infokeys.dict
    if i < length(info)
        key = nthkey(info, i)
        return key, i+1
    else
        return nothing
    end
end

function Base.iterate(info::Info, i=0)
    if i < length(info)
        key = nthkey(info, i)
        val = info[key]
        return key=>val, i+1
    else
        return nothing
    end
end


text:
get_parameter(sdf::SDFFile, p::Symbol) = getindex(sdf.param[], p)
get_parameter(sdf::SDFFile, p::Symbol, c::Symbol) = getindex(get_parameter(sdf, p), c)

get_time(sdf::SDFFile) = sdf.header.time * u"s"
get_npart(sdf::SDFFile, species) = sdf.blocks[Symbol("grid/"*species)].np

function Base.ndims(sdf::SDFFile)
    ks = string.(keys(sdf))
    i = findfirst(k->startswith(k, "grid"), ks)
    ndims(sdf.blocks[Symbol(ks[i])])
end

function domain_length(sdf::SDFFile, dir::Symbol)
    d_min = Symbol(string(dir)*"_min")
    d_max = Symbol(string(dir)*"_max")

    get_parameter(sdf, d_max) - get_parameter(sdf, d_min)
end

function domain_volume(sdf::SDFFile)
    nd = ndims(sdf)

    Δx = domain_length(sdf, :x)

    if ndims == 1
        Δx
    elseif ndims == 2
        Δy = domain_length(sdf, :y)
        Δx * Δy
    else
        Δy = domain_length(sdf, :y)
        Δz = domain_length(sdf, :z)
        Δx * Δy * Δz
    end
end

function cell_volume(sdf::SDFFile)
    nd = ndims(sdf)
    vol = domain_volume(sdf)
    nx = get_parameter(sdf, :nx)

    if ndims == 1
        vol / nx
    elseif ndims == 2
        ny = get_parameter(sdf, :ny)
        vol / (nx * ny)
    else
        ny = get_parameter(sdf, :ny)
        nz = get_parameter(sdf, :nz)
        vol / (nx * ny * nz)
    end
end


text:
module XCSP

end


text:
"""
    set_cmd!(cmds::Dict, cmd)

register `cmd` in the command registry `cmds`, which usually is
a constant `CASTED_COMMANDS` under given module.
"""
function set_cmd!(cmds::Dict, cmd, name = cmd_name(cmd))
    if haskey(cmds, name)
        @warn "replacing command $name in the registry"
    end

    return cmds[name] = cmd
end

"""
    command(module; name="", doc=docstring(m))

Convert a module to a CLI command [`NodeCommand`](@ref).

    command(f::Function, args, kwargs; name="")

Convert a function to a CLI command [`LeafCommand`](@ref).
It requires a `Vector` of function arguments and
a `Vector` of function kwargs.

The element of arguments vector is a tuple
(name, type, require):

- `name::String` is the name of argument
- `type::DataType` is the type of this argument, can be `Any` if you don't want to specify
- `require::Bool` indicates the whether this argument is required.
- `vararg::Bool` indicates whether this argument is a vararg.

The element of kwargs vector is also a tuple
(name, type, short):

- `name::String` is the name of kwarg
- `type::DataType` is the type of this kwarg, can be `Any` if you don't want to specify
- `short::Bool` is a flag that indicates whether this kwarg has short option.
"""
function command end

function command(m::Module; name = "", doc = docstring(m), line = LineNumberNode(0))
    if isempty(name) # force to have a valid name
        name = default_name(m)
    end
    NodeCommand(name, line, doc, collect(values(m.CASTED_COMMANDS)))
end

function command(f::Function, args, kwargs; name = "", line = LineNumberNode(0))
    if isempty(name) # force to have a valid name
        name = default_name(f)
    end
    md = Base.Docs.doc(f)
    intro, doc_args, flag_docs, option_docs = read_doc(md)
    args = create_args(args, doc_args, line)
    options, flags = create_options_and_flags(kwargs, option_docs, flag_docs, line)
    LeafCommand(
        f;
        name = name,
        args = args,
        options = options,
        flags = flags,
        doc = intro,
        line = line,
    )
end

function create_args(args, docs, line)
    ret = Arg[]
    for (name, type, require, vararg, default) in args
        push!(
            ret,
            Arg(
                name;
                type = type,
                require = require,
                doc = get(docs, name, ""),
                line = line,
                vararg = vararg,
                default = default,
            ),
        )
    end
    return ret
end

function create_options_and_flags(kwargs, option_docs, flag_docs, line)
    flags, options = Flag[], Option[]
    for (name, type, isflag, default) in kwargs
        if isflag
            push!(flags, create_flag(name, type, flag_docs, line))
        else
            push!(options, create_option(name, type, option_docs, line, default))
        end
    end
    return options, flags
end

function create_flag(name::String, type, flag_docs, line)
    if haskey(flag_docs, name)
        doc, short = flag_docs[name]
        return Flag(name, line, doc, short)
    else
        return Flag(name; line = line)
    end
end

function create_option(name::String, type, option_docs, line, default)
    if haskey(option_docs, name)
        arg, doc, short = option_docs[name]
        return Option(
            name,
            Arg(arg; type = type, default = default);
            short = short,
            doc = doc,
            line = line,
        )
    else
        return Option(name, Arg(; type = type, line = line, default = default); line = line)
    end
end

function main(
    m::Module = Main;
    name = default_name(m),
    doc = "",
    version = get_version(m),
    line = LineNumberNode(0),
)
    cmd = NodeCommand(name, line, doc, collect(values(m.CASTED_COMMANDS)))
    return EntryCommand(cmd; version = version, line = line)
end


text:
"""
    LinearEvolution(Φ, b, Q) <: Evolution

Evolution of the law of `x -> Φ x + w` where ``w \\sim N(0, Q)``.

# Examples
```
    evolve(LinearEvolution(Φ, b, Q), 0 => Gaussian(x, P))
```
"""
struct LinearEvolution{TΦ,TQ} <: Evolution
    Φ::TΦ # dxd
    Q::TQ # dxd
end
Phi(M::LinearEvolution, ::Any) = M.Φ
Phi(M::LinearEvolution) = M.Φ

struct GenericLinearEvolution <: Evolution
end

evolve(M::LinearEvolution, u::Pair) = timelift_evolve(M, u)
evolve(M::LinearEvolution, u::Pair, c) = timelift_evolve(M, u, c)

function evolve(M::LinearEvolution, G)
    M.Φ*G ⊕ M.Q
end

function evolve(M::GenericLinearEvolution, (G,)::Control, control)
    Φ, Q = control
    G(Φ*mean(G) + mean(Q), Φ*cov(G)*Φ' + cov(Q))
end