text:
GLAbstraction.gl_convert{T}(::Type{T}, img::Images.Image) = gl_convert(T, Images.data(img))

_default{T <: Colorant, X}(main::Images.Image{T, 2, X}, s::Style, d::Dict) = _default(Images.data(main), s, d)
_default{T <: Colorant, X}(main::Signal{Images.Image{T, 2, X}}, s::Style, d::Dict) = _default(const_lift(Images.data, main), s, d)

"""
A matrix of colors is interpreted as an image
"""
_default{T <: Colorant}(main::MatTypes{T}, ::Style, data::Dict) = @gen_defaults! data begin
    image                 = main => (Texture, "image, can be a Texture or Array of colors")
    primitive::GLUVMesh2D = SimpleRectangle{Float32}(0f0, 0f0, size(value(main))...) => "the 2D mesh the image is mapped to. Can be a 2D Geometry or mesh"
    boundingbox           = GLBoundingBox(primitive)
    preferred_camera      = :orthographic_pixel
    shader                = GLVisualizeShader("fragment_output.frag", "uv_vert.vert", "texture.frag")
end
function _default{T <: Colorant}(main::VecTypes{T}, ::Style, data::Dict)
    @gen_defaults! data begin
        image                 = main => (Texture, "image, can be a Texture or Array of colors")
        primitive::GLUVMesh2D = SimpleRectangle{Float32}(0f0, 0f0, length(value(main)), 50f0) => "the 2D mesh the image is mapped to. Can be a 2D Geometry or mesh"
        boundingbox           = GLBoundingBox(primitive)
        preferred_camera      = :orthographic_pixel
        shader                = GLVisualizeShader("fragment_output.frag", "uv_vert.vert", "texture.frag")
    end
end

"""
A matrix of Intensities will result in a contourf kind of plot
"""
function _default{T <: Intensity}(main::MatTypes{T}, s::Style, data::Dict)
    main_v = value(main)
    @gen_defaults! data begin
        ranges = (0:size(main_v, 1), 0:size(main_v, 2))
    end
    x, y, xw, yh = first(ranges[1]), first(ranges[2]), last(ranges[1]), last(ranges[2])
    @gen_defaults! data begin
        intensity             = main => Texture
        color_map             = default(Vector{RGBA{U8}},s) => Texture
        primitive::GLUVMesh2D = SimpleRectangle{Float32}(x, y, xw-x, yh-y)
        color_norm            = const_lift(extrema2f0, main)
        stroke_width::Float32 = 0.05f0
        levels::Float32       = 5f0
        stroke_color          = RGBA{Float32}(1,1,1,1)
        boundingbox           = GLBoundingBox(primitive)
        shader                = GLVisualizeShader("fragment_output.frag", "uv_vert.vert", "intensity.frag")
        preferred_camera      = :orthographic_pixel
    end
end

"""
Float matrix with the style distancefield will be interpreted as a distancefield.
A distancefield is describing a shape, with positive values denoting the inside
of the shape, negative values the outside and 0 the border
"""
function _default{T <: AbstractFloat}(main::MatTypes{T}, s::style"distancefield", data::Dict)
    @gen_defaults! data begin
        distancefield = main => Texture
        shape         = DISTANCEFIELD
    end
    rect = SimpleRectangle{Float32}(0f0,0f0, size(value(main))...)
    _default((rect, Point2f0[0]), s, data)
end

export play
"""
With play, you can slice a 3D array along `timedim` at time `t`.
This can be used to treat a 3D array like a video and create an image stream from it.
"""
function play{T}(array::Array{T, 3}, timedim::Integer, t::Integer)
    index = ntuple(dim->dim==timedim ? t : Colon(), Val{3})
    array[index...]
end
"""
Turns an Image into a video stream
"""
function play{T<:Colorant, X}(img::Images.Image{T, 3, X})
    props = img.properties
    if haskey(props, "timedim")
        timedim = props["timedim"]
        return const_lift(play, img.data, timedim, loop(1:size(img, timedim)))
    end
    error("Image has no time channel")
end
"""
Plays a video stream from VideoIO.jl. You need to supply the image `buffer`,
which will be reused for better performance.
"""
function play{T}(buffer::Array{T, 2}, video_stream, t)
    eof(video_stream) && seekstart(video_stream)
    w,h     = size(buffer)
    buffer     = reinterpret(UInt8, buffer, (3, w,h))
    read!(video_stream, buffer) # looses type and shape
    return reinterpret(T, buffer, (w,h))
end

"""
Takes a 3D image and decides if it is a volume or an animated Image.
"""
function _default{T<:Colorant, X}(img::Images.Image{T, 3, X}, s::Style, data::Dict)
    props = img.properties
    if haskey(props, "timedim")
        timedim = props["timedim"]
        video_signal = const_lift(play, img.data, timedim, loop(1:size(img, timedim)))
        return _default(video_signal, s, data)
    elseif haskey(props, "pixelspacing")
        spacing = Vec3f0(map(float, img.properties["pixelspacing"]))
        pdims   = Vec3f0(size(img))
        dims    = pdims .* spacing
        dims    = dims/maximum(dims)
        data[:dimensions] = dims
    end
    _default(img.data, s, data)
end

"""
Takes a shader as a parametric function. The shader should contain a function stubb
like this:
```GLSL
uniform float arg1; // you can add arbitrary uniforms and supply them via the keyword args
float function(float x) {
 return arg1*sin(1/tan(x));
}
```
"""
_default(func::Shader, s::Style, data::Dict) = @gen_defaults! data begin
    color                 = default(RGBA, s)  => Texture
    dimensions            = (120f0,120f0)
    primitive::GLUVMesh2D = SimpleRectangle{Float32}(0f0,0f0, dimensions...)
    preferred_camera      = :orthographic_pixel
    boundingbox           = GLBoundingBox(primitive)
    shader                = GLVisualizeShader("fragment_output.frag", "parametric.vert", "parametric.frag", view=Dict(
         "function" => Compat.String(func.source)
     ))
end

#Volumes
typealias VolumeElTypes Union{Colorant, AbstractFloat}

const default_style = Style{:default}()

function _default{T<:VolumeElTypes}(a::VolumeTypes{T}, s::Style{:iso}, data::Dict)
    data = @gen_defaults! data begin
        isovalue  = 0.5f0
        algorithm = IsoValue
    end
     _default(a, default_style, data)
end

function _default{T<:VolumeElTypes}(a::VolumeTypes{T}, s::Style{:absorption}, data::Dict)
    data = @gen_defaults! data begin
        absorption = 1f0
        algorithm  = Absorption
    end
    _default(a, default_style, data)
end

immutable VolumePrerender
end
@compat function (::VolumePrerender)()
    GLAbstraction.StandardPrerender()()
    glEnable(GL_CULL_FACE)
    glCullFace(GL_FRONT)
end

_default{T<:VolumeElTypes}(main::VolumeTypes{T}, s::Style, data::Dict) = @gen_defaults! data begin
    intensities      = main => Texture
    dimensions       = Vec3f0(1)
    hull::GLUVWMesh  = AABB{Float32}(Vec3f0(0), dimensions)
    light_position   = Vec3f0(0.25, 1.0, 3.0)
    light_intensity  = Vec3f0(15.0)

    color_map        = default(Vector{RGBA}, s) => Texture
    color_norm       = color_map == nothing ? nothing : const_lift(extrema2f0, main)
    color            = color_map == nothing ? default(RGBA, s) : nothing

    algorithm        = MaximumIntensityProjection
    boundingbox      = hull
    absorption       = 1f0
    isovalue         = 0.5f0
    isorange         = 0.01f0
    shader           = GLVisualizeShader("fragment_output.frag", "util.vert", "volume.vert", "volume.frag")
    # prerender        = VolumePrerender()
    # postrender       = () -> begin
    #     glDisable(GL_CULL_FACE)
    # end
end


text:
using Term.box

@testset "\e34mBOX" begin
    testbox = Term.box.Box(
        "ASCII",
        """
        +--+
        | ||
        |-+|
        | ||
        |-+|
        |-+|
        | ||
        +--+
        """,
    )

    @test string(testbox) == "Box\e[2m(ASCII)\e[0m"

    @suppress_out begin
        @test_nothrow println(stdout, testbox)
    end

    _bstring = Term.box.fit(Term.box.ASCII, [1, 1, 1, 1, 1, 1, 1, 1])
    @test _bstring == "+---------------+\n| | | | | | | | |\n|-+-+-+-+-+-+-+-|\n| | | | | | | | |\n|-+-+-+-+-+-+-+-|\n|-+-+-+-+-+-+-+-|\n| | | | | | | | |\n+---------------+"
end

ext: jl


text:
module KLB

using FileIO
using Libdl
using UUIDs
import Base.convert

export klbheader, loadklb, loadklb!, klbtype

const KLB_DATA_DIMS = 5
const DEFAULT_NUMTHREADS = Sys.CPU_THREADS ÷ 2

# FileIO
function load(
    file::File{format"KLB"}
    ;
    bounds::Array{UInt32} = UInt32[0],
    numthreads::Int = DEFAULT_NUMTHREADS,
    inplace::AbstractArray = UInt8[0],
    nochecks::Bool = false
    )
    if length(bounds) != 1
        if length(inplace) != 1
            loadklb!(inplace, filename(file), bounds[:,1], bounds[:,2], numthreads=numthreads, nochecks=nochecks)
        else
            loadklb(filename(file), bounds[:,1], bounds[:,2], numthreads=numthreads)
        end
    else
        if length(inplace) != 1
            loadklb!(inplace, filename(file), numthreads=numthreads, nochecks=nochecks)
        else
            loadklb(filename(file), numthreads=numthreads)
        end
    end
end

# FileIO
function save(
    file::File{format"KLB"},
    A::AbstractArray
    ;
    numthreads::Int = DEFAULT_NUMTHREADS
    )
    writearray(filename(file), A, numthreads)
end

function klbheader( file::AbstractString )
    imagesize = zeros(UInt32, KLB_DATA_DIMS)
    blocksize = zeros(UInt32, KLB_DATA_DIMS)
    pixelspacing = ones(Float32, KLB_DATA_DIMS)
    ktype = Ref{Cint}(-1)
    compressiontype = Ref{Cint}(-1)
    metadata = repeat(" ", 256)

    errid = ccall( (:readKLBheader, "klb"), Cint,
        (Cstring, Ptr{UInt32}, Ref{Cint}, Ptr{Float32}, Ptr{UInt32}, Ref{Cint}, Ptr{Cchar}),
        file, imagesize, ktype, pixelspacing, blocksize, compressiontype, metadata)

    if errid != 0
        error("Could not read KLB header of file '$file'. Error code $errid")
    end

    Dict{AbstractString, Any}(
        "imagesize" => round.(Int, imagesize),
        "blocksize" => round.(Int, blocksize),
        "pixelspacing" => pixelspacing,
        "metadata" => strip(metadata),
        "datatype" => juliatype(ktype[]),
        "compressiontype" => compressiontype[],
        "spatialorder" => ["x", "y", "z"],
        "colordim" => 4,
        "timedim" => 5)
end

function loadklb(
    file::AbstractString
    ;
    numthreads::Int = DEFAULT_NUMTHREADS
    )
    header = klbheader(file)
    A = Array{header["datatype"]}(undef, header["imagesize"]...)
    loadklb!(A, file, numthreads=numthreads, nochecks=true)
    A
end

function loadklb!(
    A::AbstractArray,
    file::AbstractString
    ;
    numthreads::Int = DEFAULT_NUMTHREADS,
    nochecks::Bool = false
    )
    if !nochecks
        header = klbheader(file)
        @assert header["datatype"] == eltype(A)
        imagesize = header["imagesize"]
        for d in 1:KLB_DATA_DIMS
            @assert imagesize[d] == size(A, d)
        end
    end

    errid = ccall( (:readKLBstackInPlace, "klb"), Cint,
        (Cstring, Ptr{Cvoid}, Ref{Cint}, Cint),
        file, A, Ref{Cint}(0), numthreads)

    if errid != 0
        error("Could not read KLB file '$file'. Error code $errid")
    end
end

function loadklb(
    file::AbstractString,
    lower_bounds::Vector{UInt32},
    upper_bounds::Vector{UInt32}
    ;
    numthreads::Int = DEFAULT_NUMTHREADS
    )
    header = klbheader(file)
    roisize = 1 .+ upper_bounds .- lower_bounds
    A = Array{header["datatype"]}(undef, roisize...)
    loadklb!(A, file, lower_bounds, upper_bounds, numthreads=numthreads, nochecks=true)
    A
end

function loadklb!(
    A::AbstractArray,
    file::AbstractString,
    lower_bounds::Union{ Vector{UInt32}, NTuple{KLB_DATA_DIMS} },
    upper_bounds::Union{ Vector{UInt32}, NTuple{KLB_DATA_DIMS} }
    ;
    numthreads::Int = DEFAULT_NUMTHREADS,
    nochecks::Bool = false
    )
    lb = lower_bounds .- UInt32(1)
    ub = upper_bounds .- UInt32(1)

    if !nochecks
        header = klbheader(file)
        @assert header["datatype"] == eltype(A)
        roisize = 1 .+ ub .- lb
        for d in 1:5
            @assert roisize[d] == size(A, d)
        end
    end

    errid = _readKLBroiInPlace(A, file, lb, ub, numthreads)

    if errid != 0
        error("Could not read KLB file '$file'. Error code $errid")
    end
    A
end

# Vector UInt32
function _readKLBroiInPlace(
    A::AbstractArray,
    file::AbstractString,
    lb::Vector{UInt32},
    ub::Vector{UInt32},
    numthreads::Int
)
    ccall( (:readKLBroiInPlace, "klb"), Cint,
        (Cstring, Ptr{Cvoid}, Ptr{UInt32}, Ptr{UInt32}, Cint),
        file, A, lb, ub, numthreads)
end

# NTuple{KLB_DATA_DIMS,T}
function _readKLBroiInPlace(
    A::AbstractArray,
    file::AbstractString,
    lb::NTuple{KLB_DATA_DIMS,T},
    ub::NTuple{KLB_DATA_DIMS,T},
    numthreads::Int
) where T
    ccall( (:readKLBroiInPlace, "klb"), Cint,
        (Cstring, Ptr{Cvoid}, NTuple{5,UInt32}, NTuple{5,UInt32}, Cint),
        file, A, (lb...,), (ub...,), numthreads)
end

function writearray(
    filepath::AbstractString,
    A::AbstractArray,
    numthreads::Int = DEFAULT_NUMTHREADS
    ;
    pixelspacing = C_NULL,
    blocksize = C_NULL,
    compressiontype::Int = 1,
    metadata = C_NULL
    )
    ktype = klbtype(eltype(A))
    imagesize = UInt32[i for i in size(A)]
    while length(imagesize) < 5
        push!(imagesize, 1)
    end

    errid = ccall( (:writeKLBstack, "klb"), Cint,
        (Ptr{Cvoid}, Cstring, Ptr{UInt32}, Cint, Cint, Ptr{Float32}, Ptr{UInt32}, Cint, Ptr{Cchar}),
        A, filepath, imagesize, ktype, numthreads, pixelspacing, blocksize, compressiontype, metadata)

    if errid != 0
        error("Could not write KLB file '$filepath'. Error code $errid")
    end
end

@enum klbtype::Cint begin
    KLB_UInt8   = 0
    KLB_UInt16  = 1
    KLB_UInt32  = 2
    KLB_UInt64  = 3
    KLB_Int8    = 4
    KLB_Int16   = 5
    KLB_Int32   = 6
    KLB_Int64   = 7
    KLB_Float32 = 8
    KLB_Float64 = 9
end

function juliatype( _klbtype )
    _klbtype = klbtype(_klbtype)
    if     _klbtype == KLB_UInt8  # 0
        return UInt8
    elseif _klbtype == KLB_UInt16 # 1
        return UInt16
    elseif _klbtype == KLB_UInt32 # 2
        return UInt32
    elseif _klbtype == KLB_UInt64 # 3
        return UInt64
    elseif _klbtype == KLB_Int8   # 4
        return Int8
    elseif _klbtype == KLB_Int16  # 5
        return Int16
    elseif _klbtype == KLB_Int32  # 6
        return Int32
    elseif _klbtype == KLB_Int64  # 7
        return Int64
    elseif _klbtype == KLB_Float32 # 8
        return Float32
    elseif _klbtype == KLB_Float64 # 9
        return Float64
    end
    error( "Unknown or unsupported data type of KLB array: $klbtype" )
end

klbtype(::Type{UInt8}  )  = KLB_UInt8   # 0
klbtype(::Type{UInt16} )  = KLB_UInt16  # 1
klbtype(::Type{UInt32} )  = KLB_UInt32  # 2
klbtype(::Type{UInt64} )  = KLB_UInt64  # 3
klbtype(::Type{Int8}   )  = KLB_Int8    # 4
klbtype(::Type{Int16}  )  = KLB_Int16   # 5
klbtype(::Type{Int32}  )  = KLB_Int32   # 6
klbtype(::Type{Int64}  )  = KLB_Int64   # 7
klbtype(::Type{Float32})  = KLB_Float32 # 8
klbtype(::Type{Float64})  = KLB_Float64 # 9
klbtype(k::klbtype) = k

function klbtype( juliatype::Type )
    #=
    if juliatype == UInt8
        return 0
    elseif juliatype == UInt16
        return 1
    elseif juliatype == UInt32
        return 2
    elseif juliatype == UInt64
        return 3
    elseif juliatype == Int8
        return 4
    elseif juliatype == Int16
        return 5
    elseif juliatype == Int32
        return 6
    elseif juliatype == Int64
        return 7
    elseif juliatype == Float32
        return 8
    elseif juliatype == Float64
        return 9
    end
    =#
    error( "Unknown or unsupported data type of KLB array: $juliatype" )
end

convert(::DataType, klb::klbtype ) = juliatype(klb)

function __init__()
    klb_deps_folder = joinpath(dirname(dirname(@__FILE__)), "deps")
    global library_location
    library_location = Libdl.find_library(["klb","libklb"],[klb_deps_folder])
    if isempty(library_location)
        error("[lib]klb.$(Libdl.dlext) library could not be located in $klb_deps_folder or elsewhere.")
    else
        @info "$library_location.$(Libdl.dlext)"
        Libdl.dlopen(library_location)
    end
    add_format(format"KLB", (), ".klb", [:KLB => UUID("8bb66c0d-974d-412d-8f46-f9b8d1ef37d0")])
end

end # module


text:
@testset "DemoSection" begin
    root = joinpath("assets", "section")

    # default section behavior
    sec = DemoSection(joinpath(root, "default"))
    @test sec.title == "Default"
    @test sec.cards == []
    subsec1, subsec2 = sec.subsections
    @test subsec1.title == "Subsection 1"
    @test subsec2.title == "Subsection 2"
    @test subsec1.subsections == []
    @test subsec2.subsections == []
    card1, card2 = subsec1.cards
    card3, card4 = subsec2.cards
    @test card1.path == joinpath(subsec1.root, "card_1.md")
    @test card2.path == joinpath(subsec1.root, "card_2.md")
    @test card3.path == joinpath(subsec2.root, "card_3.md")
    @test card4.path == joinpath(subsec2.root, "card_4.md")

    # only config subsection_1
    sec = DemoSection(joinpath(root, "partial_config"))
    @test sec.title == "Partial Config"
    @test sec.cards == []
    subsec1, subsec2 = sec.subsections
    @test subsec1.title == "Subsection 1"
    @test subsec2.title == "Subsection 2"
    @test subsec1.subsections == []
    @test subsec2.subsections == []
    card1, card2 = subsec1.cards
    card3, card4 = subsec2.cards
    @test card1.path == joinpath(subsec1.root, "card_2.md")
    @test card2.path == joinpath(subsec1.root, "card_1.md")
    @test card3.path == joinpath(subsec2.root, "card_3.md")
    @test card4.path == joinpath(subsec2.root, "card_4.md")

    # invalid cases
    @test_throws ArgumentError @suppress_err DemoSection(joinpath(root, "partial_order"))
    @test_throws ArgumentError @suppress_err DemoSection(joinpath(root, "cards_and_subsections"))

    @testset "unmatched dirty files" begin
        msg = @capture_err DemoSection(joinpath(root, "dirty"))
        @test occursin("skip unmatched file: \"data.csv\", \"script.py\"", msg)
        sec = @suppress_err DemoSection(joinpath(root, "dirty"))

        # unmatched files are excluded during the structure building stage.
        @test length(sec.cards) == 1
    end
end


text:
AccFree_MTiXsens_sync=Dict()
Xsensor_reduce_sync=Dict()
IBS_sync = Dict()

for experiment=1:size(filenameMTiXsens,1)
    AccFree_MTiXsens_sync[filenameMTiXsens[experiment,1]]=Dict()
    # @info experiment
    for n in sensorXsensNames
        arr = AccFree_MTiXsens[filenameMTiXsens[experiment,1]][n][:,"PacketCounter"]
        ligne = findall(x -> x == sync.MTiXsens[experiment], arr)
        AccFree_MTiXsens_sync[filenameMTiXsens[experiment,1]][n]=AccFree_MTiXsens[filenameMTiXsens[experiment,1]][n][ligne[1]:end,:]
    end
end
for experiment=1:size(filenameXsensor,1)
    Xsensor_reduce_sync[filenameXsensor[experiment,1]]=Xsensor_reduce[filenameXsensor[experiment,1]][sync.Xsensor[experiment]:end,:]
end
for experiment=1:size(filenameIBS,1)
    IBS_sync[filenameIBS[experiment,1]] = IBS[filenameIBS[experiment,1]][sync.IBS[experiment]:end,:]
end

println("Synchro data finished")


text:
module MyModule
export mycoolfunction

function mycoolfunction()
   println("this is my cool function")
end
function mysecretfunction()
   println("this is my secret function")
end

end

max_stars_repo_name: nnhjy/julia-introduction

max_issues_repo_name: nnhjy/julia-introduction

max_forks_repo_name: nnhjy/julia-introduction


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

"""
Internal abstract type.
"""
abstract type AbstractSimulation end

"""
Internal unimplemented exception type.
"""
Base.@kwdef struct UnimplementedError <: Exception
    msg::String = "Required function has not been implemented."
end
Base.showerror(io::IO, e::UnimplementedError) = print(io, e.msg)
unimplemented() = throw(UnimplementedError())

"""
Set of environment variables and their distributions.
"""
const Environment = Dict{Symbol, Sampleable}

"""
Samples from environment.
"""
const EnvironmentValue = Dict{Symbol, Any}

Base.zero(val::EnvironmentValue) = EnvironmentValue(k => zero(v) for (k, v) in val)

"""
Environment that will be flattened upon sampling.
"""
mutable struct Flat
    env::Environment
end

"""
Inferred properties of single environment variable.
"""
struct VariableInfo
	n::Int64	# environment variable dimensionality
	t::Type		# environment variable type
end

"""
Properties of environment variables.
"""
const EnvironmentInfo = Dict{Symbol, VariableInfo}

"""
Abstract type for AST state.
"""
abstract type State end

"""
State type for observable simulation.
"""
struct ObservableState <: State end

"""
State type for unobservable simulation.
"""
struct UnobservableState <: State end

"""
Abstract type for AST action.
"""
abstract type Action end

"""
Action type corresponding to instantiation of stochastic environment.
"""
struct SampleAction <: Action
    sample::EnvironmentValue
end

"""
Action type corresponding to setting random seed of simulation.
"""
struct SeedAction <: Action
    seed::UInt32
end

"""
Abstract AST MDP type.
"""
abstract type AbstractASTMDP{S<:State, A<:Action} <: CommonRLInterface.AbstractEnv end

"""
Return state type.
"""
state_type(::AbstractASTMDP{S, A}) where {S, A} = S

"""
Return action type.
"""
action_type(::AbstractASTMDP{S, A}) where {S, A} = A

"""
Abstract AST reward structure.
"""
abstract type AbstractReward end

"""
Standard MDP object for AST. Wraps simulation and contains auxiliary information and parameters.
"""
Base.@kwdef mutable struct ASTMDP{S<:State, A<:Action} <: AbstractASTMDP{S, A}
	sim::AbstractSimulation						      # simulation wrapping system under test
    reward::AbstractReward    = Reward()              # reward structure
    episodic::Bool            = false                 # episodic evaluation
    num_steps::Int64          = 0                     # number of steps (zero if unknown)
	env_info::EnvironmentInfo = EnvironmentInfo()	  # inferred environment properties
    flatten::Bool             = false                 # whether environment is flattened
    rng::AbstractRNG          = Random.default_rng()  # RNG used for simulation
end


text:
###
### Hamiltonian Monte Carlo samplers.
###

"""
    HMC(n_iters::Int, ϵ::Float64, n_leapfrog::Int)

Hamiltonian Monte Carlo sampler.

Arguments:

- `n_iters::Int` : The number of samples to pull.
- `ϵ::Float64` : The leapfrog step size to use.
- `n_leapfrog::Int` : The number of leapfrop steps to use.

Usage:

```julia
HMC(1000, 0.05, 10)
```

Tips:

- If you are receiving gradient errors when using `HMC`, try reducing the
`step_size` parameter.

```julia
# Original step_size
sample(gdemo([1.5, 2]), HMC(1000, 0.1, 10))

# Reduced step_size.
sample(gdemo([1.5, 2]), HMC(1000, 0.01, 10))
```
"""
mutable struct HMC{AD, space, metricT <: AHMC.AbstractMetric} <: StaticHamiltonian{AD}
    n_iters     ::  Int       # number of samples
    ϵ           ::  Float64   # leapfrog step size
    n_leapfrog  ::  Int       # leapfrog step number
end
HMC(args...) = HMC{ADBackend()}(args...)
function HMC{AD}(n_iters::Int, ϵ::Float64, n_leapfrog::Int, ::Type{metricT}, space::Tuple) where {AD, metricT <: AHMC.AbstractMetric}
    return HMC{AD, space, metricT}(n_iters, ϵ, n_leapfrog)
end
function HMC{AD}(
    n_iters::Int,
    ϵ::Float64,
    n_leapfrog::Int,
    ::Tuple{};
    kwargs...
) where AD
    return HMC{AD}(n_iters, ϵ, n_leapfrog; kwargs...)
end
function HMC{AD}(
    n_iters::Int,
    ϵ::Float64,
    n_leapfrog::Int,
    space::Symbol...;
    metricT=AHMC.UnitEuclideanMetric
) where AD
    return HMC{AD}(n_iters, ϵ, n_leapfrog, metricT, space)
end

"""
    HMCDA(n_iters::Int, n_adapts::Int, δ::Float64, λ::Float64; init_ϵ::Float64=0.1)

Hamiltonian Monte Carlo sampler with Dual Averaging algorithm.

Usage:

```julia
HMCDA(1000, 200, 0.65, 0.3)
```

Arguments:

- `n_iters::Int` : Number of samples to pull.
- `n_adapts::Int` : Numbers of samples to use for adaptation.
- `δ::Float64` : Target acceptance rate. 65% is often recommended.
- `λ::Float64` : Target leapfrop length.
- `init_ϵ::Float64=0.1` : Inital step size; 0 means automatically search by Turing.

For more information, please view the following paper ([arXiv link](https://arxiv.org/abs/1111.4246)):

- Hoffman, Matthew D., and Andrew Gelman. "The No-U-turn sampler: adaptively
  setting path lengths in Hamiltonian Monte Carlo." Journal of Machine Learning
  Research 15, no. 1 (2014): 1593-1623.
"""
mutable struct HMCDA{AD, space, metricT <: AHMC.AbstractMetric} <: AdaptiveHamiltonian{AD}
    n_iters     ::  Int       # number of samples
    n_adapts    ::  Int       # number of samples with adaption for ϵ
    δ           ::  Float64   # target accept rate
    λ           ::  Float64   # target leapfrog length
    init_ϵ      ::  Float64
end
HMCDA(args...; kwargs...) = HMCDA{ADBackend()}(args...; kwargs...)
function HMCDA{AD}(n_iters::Int, n_adapts::Int, δ::Float64, λ::Float64, init_ϵ::Float64, ::Type{metricT}, space::Tuple) where {AD, metricT <: AHMC.AbstractMetric}
    return HMCDA{AD, space, metricT}(n_iters, n_adapts, δ, λ, init_ϵ)
end

function HMCDA{AD}(
    n_iters::Int,
    δ::Float64,
    λ::Float64;
    init_ϵ::Float64=0.0,
    metricT=AHMC.UnitEuclideanMetric
) where AD
    n_adapts_default = Int(round(n_iters / 2))
    n_adapts = n_adapts_default > 1000 ? 1000 : n_adapts_default
    return HMCDA{AD}(n_iters, n_adapts, δ, λ, init_ϵ,metricT, ())
end

function HMCDA{AD}(
    n_iters::Int,
    n_adapts::Int,
    δ::Float64,
    λ::Float64,
    ::Tuple{};
    kwargs...
) where AD
    return HMCDA{AD}(n_iters, n_adapts, δ, λ; kwargs...)
end

function HMCDA{AD}(
    n_iters::Int,
    n_adapts::Int,
    δ::Float64,
    λ::Float64,
    space::Symbol...;
    init_ϵ::Float64=0.0,
    metricT=AHMC.UnitEuclideanMetric
) where AD
    return HMCDA{AD}(n_iters, n_adapts, δ, λ, init_ϵ, metricT, space)
end

"""
    NUTS(n_iters::Int, n_adapts::Int, δ::Float64; max_depth::Int=5, Δ_max::Float64=1000.0, init_ϵ::Float64=0.1)

No-U-Turn Sampler (NUTS) sampler.

Usage:

```julia
NUTS(1000, 200, 0.6j_max)
```

Arguments:

- `n_iters::Int` : The number of samples to pull.
- `n_adapts::Int` : The number of samples to use with adapatation.
- `δ::Float64` : Target acceptance rate.
- `max_depth::Float64` : Maximum doubling tree depth.
- `Δ_max::Float64` : Maximum divergence during doubling tree.
- `init_ϵ::Float64` : Inital step size; 0 means automatically search by Turing.

"""
mutable struct NUTS{AD, space, metricT <: AHMC.AbstractMetric} <: AdaptiveHamiltonian{AD}
    n_iters     ::  Int       # number of samples
    n_adapts    ::  Int       # number of samples with adaption for ϵ
    δ           ::  Float64   # target accept rate
    max_depth   ::  Int
    Δ_max       ::  Float64
    init_ϵ      ::  Float64
end

NUTS(args...; kwargs...) = NUTS{ADBackend()}(args...; kwargs...)
function NUTS{AD}(
    n_iters::Int,
    n_adapts::Int,
    δ::Float64,
    max_depth::Int,
    Δ_max::Float64,
    init_ϵ::Float64,
    ::Type{metricT},
    space::Tuple
) where {AD, metricT}
    return NUTS{AD, space, metricT}(n_iters, n_adapts, δ, max_depth, Δ_max, init_ϵ)
end

function NUTS{AD}(
    n_iters::Int,
    n_adapts::Int,
    δ::Float64,
    ::Tuple{};
    kwargs...
) where AD
    NUTS{AD}(n_iters, n_adapts, δ; kwargs...)
end

function NUTS{AD}(
    n_iters::Int,
    n_adapts::Int,
    δ::Float64,
    space::Symbol...;
    max_depth::Int=10,
    Δ_max::Float64=1000.0,
    init_ϵ::Float64=0.0,
    metricT=AHMC.DiagEuclideanMetric
) where AD
    NUTS{AD}(n_iters, n_adapts, δ, max_depth, Δ_max, init_ϵ, metricT, space)
end

function NUTS{AD}(
    n_iters::Int,
    δ::Float64;
    max_depth::Int=10,
    Δ_max::Float64=1000.0,
    init_ϵ::Float64=0.0,
    metricT=AHMC.DiagEuclideanMetric
) where AD
    n_adapts_default = Int(round(n_iters / 2))
    NUTS{AD}(n_iters, n_adapts_default > 1000 ?
        1000 : n_adapts_default, δ, max_depth, Δ_max, init_ϵ, metricT, ())
end

for alg in (:HMC, :HMCDA, :NUTS)
    @eval getmetricT(::$alg{<:Any, <:Any, metricT}) where {metricT} = metricT
end

####
#### Sampler construction
####

# Sampler(alg::Hamiltonian) =  Sampler(alg, AHMCAdaptor())
function Sampler(alg::Hamiltonian, s::Selector=Selector())
    info = Dict{Symbol, Any}()

    info[:eval_num] = 0
    info[:i] = 0

    Sampler(alg, info, s)
end

function sample(
    model::Model,
    alg::Hamiltonian;
    save_state=false,                                   # flag for state saving
    resume_from=nothing,                                # chain to continue
    reuse_spl_n=0,                                      # flag for spl re-using
    adaptor=AHMCAdaptor(alg),
    init_theta::Union{Nothing,Array{<:Any,1}}=nothing,
    rng::AbstractRNG=GLOBAL_RNG,
    discard_adapt::Bool=true,
    verbose::Bool=true,
    progress::Bool=false,
    kwargs...
)
    # Create sampler
    spl = reuse_spl_n > 0 ? resume_from.info[:spl] : Sampler(alg)
    @assert isa(spl.alg, Hamiltonian) "[Turing] alg type mismatch; please use resume() to re-use spl"

    # Resume selector
    resume_from != nothing && (spl.selector = resume_from.info[:spl].selector)

    # TODO: figure out what does this line do
    n = reuse_spl_n > 0 ? reuse_spl_n : alg.n_iters

    # Init samples
    samples = Vector{Sample}(undef, n)
    weight = 1 / n
    for i = 1:n
        samples[i] = Sample(weight, Dict{Symbol, Any}())
    end

    # Create VarInfo
    vi = if resume_from == nothing
        VarInfo(model)
    else
        deepcopy(resume_from.info[:vi])
    end

    # Get `init_theta`
    if init_theta != nothing
        verbose && @info "Using passed-in initial variable values" init_theta
        # Convert individual numbers to length 1 vector; `ismissing(v)` is needed as `size(missing)` is undefined`
        init_theta = [ismissing(v) || size(v) == () ? [v] : v for v in init_theta]
        # Flatten `init_theta`
        init_theta_flat = foldl(vcat, map(vec, init_theta))
        # Create a mask to inidicate which values are not missing
        theta_mask = map(x -> !ismissing(x), init_theta_flat)
        # Get all values
        theta = vi[spl]
        @assert length(theta) == length(init_theta_flat) "Provided initial value doesn't match the dimension of the model"
        # Update those which are provided (i.e. not missing)
        theta[theta_mask] .= init_theta_flat[theta_mask]
        # Update in `vi`
        vi[spl] = theta
    end

    # Convert to transformed space
    if spl.selector.tag == :default
        link!(vi, spl)
        runmodel!(model, vi, spl)
    end

    # Init h, prop and adaptor
    step(model, spl, vi, Val(true); rng=rng, adaptor=adaptor)

    # Sampling using AHMC and store samples in `samples`
    steps!(model, spl, vi, samples; rng=rng, verbose=verbose, progress=progress)

    # Concatenate samples
    if resume_from != nothing
        pushfirst!(samples, resume_from.info[:samples]...)
    end

    # Wrap the result by Chain
    c = typeof(alg) <: AdaptiveHamiltonian && discard_adapt ?
        Chain(0.0, samples[(alg.n_adapts+1):end]) :
        Chain(0.0, samples)

    # Save state
    if save_state
        # Convert vi back to X if vi is required to be saved.
        spl.selector.tag == :default && invlink!(vi, spl)
        c = save(c, spl, model, vi, samples)
    end

    return c
end

####
#### Transition / step functions for HMC samplers.
####

# Init for StaticHamiltonian
function step(
    model,
    spl::Sampler{<:StaticHamiltonian},
    vi::VarInfo,
    is_first::Val{true};
    kwargs...
)
    ∂logπ∂θ = gen_∂logπ∂θ(vi, spl, model)
    logπ = gen_logπ(vi, spl, model)
    spl.info[:h] = AHMC.Hamiltonian(getmetricT(spl.alg)(length(vi[spl])), logπ, ∂logπ∂θ)
    spl.info[:traj] = gen_traj(spl.alg, spl.alg.ϵ)
    return vi, true
end

# Init for AdaptiveHamiltonian
function step(
    model,
    spl::Sampler{<:AdaptiveHamiltonian},
    vi::VarInfo,
    is_first::Val{true};
    adaptor=AHMCAdaptor(spl.alg),
    rng::AbstractRNG=GLOBAL_RNG,
    kwargs...
)
    spl.selector.tag != :default && link!(vi, spl)

    ∂logπ∂θ = gen_∂logπ∂θ(vi, spl, model)
    logπ = gen_logπ(vi, spl, model)

    θ_init = Vector{Float64}(vi[spl])
    metric = getmetricT(spl.alg)(length(θ_init))
    h = AHMC.Hamiltonian(metric, logπ, ∂logπ∂θ)
    init_ϵ = spl.alg.init_ϵ

    # Find good eps if not provided one
    if init_ϵ == 0.0
        init_ϵ = AHMC.find_good_eps(rng, h, θ_init)
        @info "Found initial step size" init_ϵ
    end
    if AHMC.getϵ(adaptor) == 0.0
        adaptor = AHMCAdaptor(spl.alg; init_ϵ=init_ϵ)
    end

    spl.info[:h] = h
    spl.info[:traj] = gen_traj(spl.alg, init_ϵ)
    spl.info[:adaptor] = adaptor

    spl.selector.tag != :default && invlink!(vi, spl)
    return vi, true
end

# Single step for Gibbs compatible HMC sampling.
function step(
    model,
    spl::Sampler{<:Hamiltonian},
    vi::VarInfo,
    is_first::Val{false}
)
    # Get step size
    ϵ = :adaptor in keys(spl.info) ? AHMC.getϵ(spl.info[:adaptor]) : spl.alg.ϵ

    spl.info[:i] += 1
    spl.info[:eval_num] = 0

    Turing.DEBUG && @debug "current ϵ: $ϵ"

    Turing.DEBUG && @debug "X-> R..."
    if spl.selector.tag != :default
        link!(vi, spl)
        runmodel!(model, vi, spl)
    end

    grad_func = gen_∂logπ∂θ(vi, spl, model)
    lj_func = gen_logπ(vi, spl, model)
    metric = gen_metric(length(vi[spl]), spl)

    θ, lj = vi[spl], vi.logp

    θ_new, lj_new, is_accept, α = hmc_step(θ, lj_func, grad_func, ϵ, spl.alg, metric)

    Turing.DEBUG && @debug "decide whether to accept..."
    if is_accept
        vi[spl] = θ_new
        setlogp!(vi, lj_new)
    else
        vi[spl] = θ
        setlogp!(vi, lj)
    end

    if PROGRESS[] && spl.selector.tag == :default
        haskey(spl.info, :progress) && ProgressMeter.update!(
            spl.info[:progress],
            spl.info[:progress].counter;
            showvalues = [(:ϵ, ϵ), (:α, α), (:metric, metric)],
        )
    end

    if spl.alg isa AdaptiveHamiltonian
        if spl.info[:i] <= spl.alg.n_adapts
            AHMC.adapt!(spl.info[:adaptor], Vector{Float64}(vi[spl]), α)
        end
    end

    Turing.DEBUG && @debug "R -> X..."
    spl.selector.tag != :default && invlink!(vi, spl)

    return vi, is_accept
end

# Efficient multiple step sampling for adaptive HMC.
function steps!(
    model,
    spl::Sampler{<:AdaptiveHamiltonian},
    vi,
    samples;
    rng::AbstractRNG=GLOBAL_RNG,
    verbose::Bool=true,
    progress::Bool=false
)
    ahmc_samples, stats = AHMC.sample(
        rng,
        spl.info[:h],
        spl.info[:traj],
        Vector{Float64}(vi[spl]),
        spl.alg.n_iters,
        spl.info[:adaptor],
        spl.alg.n_adapts;
        verbose=verbose,
        progress=progress
    )
    for i = 1:length(samples)
        vi[spl] = ahmc_samples[i]
        samples[i].value = Sample(vi, spl).value
        foreach(name -> samples[i].value[name] = stats[i][name], typeof(stats[i]).names)
    end
end

# Efficient multiple step sampling for static HMC.
function steps!(
    model,
    spl::Sampler{<:HMC},
    vi,
    samples;
    rng::AbstractRNG=GLOBAL_RNG,
    verbose::Bool=true,
    progress::Bool=false
)
    ahmc_samples, stats = AHMC.sample(
        rng,
        spl.info[:h],
        spl.info[:traj],
        Vector{Float64}(vi[spl]),
        spl.alg.n_iters;
        verbose=verbose,
        progress=progress
    )

    for i = 1:length(samples)
        vi[spl] = ahmc_samples[i]
        samples[i].value = Sample(vi, spl).value
        foreach(name -> samples[i].value[name] = stats[i][name], typeof(stats[i]).names)
    end
end

# Default multiple step sampling for all HMC samplers.
function steps!(
    model,
    spl::Sampler{<:Hamiltonian},
    vi,
    samples;
    rng::AbstractRNG=GLOBAL_RNG,
    verbose::Bool=true,
    progress::Bool=false
)
    # Init step
    time_elapsed = @elapsed vi, is_accept = step(model, spl, vi, Val(true))
    samples[1].value = Sample(vi, spl).value    # we know we always accept the init step
    samples[1].value[:elapsed] = time_elapsed
    # Rest steps
    for i = 2:length(samples)
        time_elapsed = @elapsed vi, is_accept = step(model, spl, vi, Val(false))
        if is_accept # accepted => store the new predcits
            samples[i].value = Sample(vi, spl).value
        else         # rejected => store the previous predcits
            samples[i] = samples[i - 1]
        end
        samples[i].value[:elapsed] = time_elapsed
    end
end

#####
##### HMC core functions
#####

"""
    gen_∂logπ∂θ(vi::VarInfo, spl::Sampler, model)

Generate a function that takes a vector of reals `θ` and compute the logpdf and
gradient at `θ` for the model specified by `(vi, spl, model)`.
"""
function gen_∂logπ∂θ(vi::VarInfo, spl::Sampler, model)
    function ∂logπ∂θ(x)
        return gradient_logp(x, vi, model, spl)
    end
    return ∂logπ∂θ
end

"""
    gen_logπ(vi::VarInfo, spl::Sampler, model)

Generate a function that takes `θ` and returns logpdf at `θ` for the model specified by
`(vi, spl, model)`.
"""
function gen_logπ(vi::VarInfo, spl::Sampler, model)
    function logπ(x)::Float64
        x_old, lj_old = vi[spl], vi.logp
        vi[spl] = x
        runmodel!(model, vi, spl)
        lj = vi.logp
        vi[spl] = x_old
        setlogp!(vi, lj_old)
        return lj
    end
    return logπ
end

gen_metric(dim::Int, spl::Sampler{<:Hamiltonian}) = AHMC.UnitEuclideanMetric(dim)
gen_metric(dim::Int, ::AHMC.UnitPreconditioner)   = AHMC.UnitEuclideanMetric(dim)
gen_metric(::Int, pc::AHMC.DiagPreconditioner)    = AHMC.DiagEuclideanMetric(AHMC.getM⁻¹(pc))
gen_metric(::Int, pc::AHMC.DensePreconditioner)   = AHMC.DenseEuclideanMetric(AHMC.getM⁻¹(pc))
gen_metric(dim::Int, spl::Sampler{<:AdaptiveHamiltonian}) = gen_metric(dim, spl.info[:adaptor].pc)

gen_traj(alg::HMC, ϵ) = AHMC.StaticTrajectory(AHMC.Leapfrog(ϵ), alg.n_leapfrog)
gen_traj(alg::HMCDA, ϵ) = AHMC.HMCDA(AHMC.Leapfrog(ϵ), alg.λ)
gen_traj(alg::NUTS, ϵ) = AHMC.NUTS(AHMC.Leapfrog(ϵ), alg.max_depth, alg.Δ_max)

function hmc_step(
    θ,
    logπ,
    ∂logπ∂θ,
    ϵ,
    alg::T,
    metric
) where {T<:Union{HMC,HMCDA,NUTS}}
    # Make sure the code in AHMC is type stable
    θ = Vector{Float64}(θ)

    # Build Hamiltonian type and trajectory
    h = AHMC.Hamiltonian(metric, logπ, ∂logπ∂θ)
    traj = gen_traj(alg, ϵ)

    h = AHMC.update(h, θ) # Ensure h.metric has the same dim as θ.

    # Sample momentum
    r = AHMC.rand(h.metric)

    # Build phase point
    z = AHMC.phasepoint(h, θ, r)

    # Call AHMC to make one MCMC transition
    z_new, stat = AHMC.transition(traj, h, z)

    return z_new.θ, stat.log_density, stat.is_accept, stat.acceptance_rate
end

####
#### Compiler interface, i.e. tilde operators.
####
using Tracker
function assume(spl::Sampler{<:Hamiltonian},
    dist::Distribution,
    vn::VarName,
    vi::VarInfo
)
    Turing.DEBUG && @debug "assuming..."
    updategid!(vi, vn, spl)
    r = vi[vn]
    # acclogp!(vi, logpdf_with_trans(dist, r, istrans(vi, vn)))
    # r
    Turing.DEBUG && @debug "dist = $dist"
    Turing.DEBUG && @debug "vn = $vn"
    Turing.DEBUG && @debug "r = $r" "typeof(r)=$(typeof(r))"
    r, logpdf_with_trans(dist, r, istrans(vi, vn))
end

function assume(spl::Sampler{<:Hamiltonian},
    dists::Vector{<:Distribution},
    vn::VarName,
    var::Any,
    vi::VarInfo
)
    @assert length(dists) == 1 "[observe] Turing only support vectorizing i.i.d distribution"
    dist = dists[1]
    n = size(var)[end]

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

    rs = vi[vns]  # NOTE: inside Turing the Julia conversion should be sticked to

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

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

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

observe(spl::Sampler{<:Hamiltonian},
    d::Distribution,
    value::Any,
    vi::VarInfo) = observe(nothing, d, value, vi)

observe(spl::Sampler{<:Hamiltonian},
    ds::Vector{<:Distribution},
    value::Any,
    vi::VarInfo) = observe(nothing, ds, value, vi)

####
#### Default HMC stepsize and mass matrix adaptor
####

function AHMCAdaptor(alg::AdaptiveHamiltonian; init_ϵ=alg.init_ϵ)
    p = AHMC.Preconditioner(getmetricT(alg))
    nda = AHMC.NesterovDualAveraging(alg.δ, init_ϵ)
    if getmetricT(alg) == AHMC.UnitEuclideanMetric
        adaptor = AHMC.NaiveHMCAdaptor(p, nda)
    else
        adaptor = AHMC.StanHMCAdaptor(alg.n_adapts, p, nda)
    end
    return adaptor
end

AHMCAdaptor(alg::Hamiltonian) = nothing


text:
# Testing the LeafModels struct
A = Fvcb()
g0 = 0.03; g1 = 12.0
Gs = Medlyn(g0,g1) # Instance of a Medlyn type with g0 = 0.03 and g1 = 0.1

@testset "LeafModels()" begin
    leaf = LeafModels(photosynthesis = A, stomatal_conductance = Gs)
    @test typeof(leaf) == LeafModels{Missing,Missing,Fvcb{Float64},Medlyn{Float64},
        MutableNamedTuples.MutableNamedTuple{(:PPFD, :Tₗ, :Cₛ, :A, :Gₛ, :Cᵢ, :Dₗ),
        NTuple{7,Base.RefValue{Float64}}}}

    @test typeof(leaf.photosynthesis) == Fvcb{Float64}
    @test typeof(leaf.stomatal_conductance) == Medlyn{Float64}
    @test leaf.photosynthesis.Tᵣ == 25.0
    @test leaf.stomatal_conductance.g0 ≈ g0
    @test leaf.stomatal_conductance.g1 ≈ g1
end;

@testset "init_status!" begin
    leaf = LeafModels(photosynthesis = A, stomatal_conductance = Gs)
    @test leaf.status.Tₗ == -999.99

    init_status!(leaf, Tₗ = 25.0)
    @test leaf.status.Tₗ == 25.0

    model = read_model("inputs/models/plant_coffee.yml")
    init_status!(model, Tₗ = 25.0)
    @test model["Leaf"].status.Tₗ == 25.0
end;

@testset "Vars to initialise" begin
    leaf = LeafModels(photosynthesis = A, stomatal_conductance = Gs)
    @test to_initialise(leaf) == [:PPFD, :Tₗ, :Cₛ, :Dₗ]
    @test to_initialise(leaf) == to_initialise(A,Gs)
    @test to_initialise(A) == [:PPFD, :Tₗ, :Cₛ]

    @test leaf.status.Tₗ == -999.99
    @test is_initialised(leaf) == false

    leaf = LeafModels(photosynthesis = A, stomatal_conductance = Gs,
        Tₗ = 25.0, PPFD = 1000.0, Cₛ = 400.0, Dₗ = 1.2)

    @test is_initialised(leaf) == true
end;


text:
include("utilities.jl")

using SEAL
using Printf

"""
    example_levels()

Illustrate the concept of levels in the BFV and CKKS schemes.
This function is based on the file `native/examples/3_levels.cpp` of the original
SEAL library and should yield the exact same output.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/3_levels.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/3_levels.cpp)

See also: [`example_bfv_basics`](@ref), [`example_ckks_basics`](@ref)
"""
function example_levels()
  print_example_banner("Example: Levels")

  enc_parms = EncryptionParameters(SchemeType.bfv)

  poly_modulus_degree = 8192
  set_poly_modulus_degree!(enc_parms, poly_modulus_degree)

  set_coeff_modulus!(enc_parms, coeff_modulus_create(poly_modulus_degree, [50, 30, 30, 50, 50]))

  set_plain_modulus!(enc_parms, plain_modulus_batching(poly_modulus_degree, 20))

  context = SEALContext(enc_parms)
  print_parameters(context)
  println()

  print_line(@__LINE__)
  println("Print the modulus switching chain.")

  context_data = key_context_data(context)
  println("----> Level (chain index): ", chain_index(context_data), " ...... key_context_data()")
  print("      parms_id:")
  for parm_id in parms_id(context_data)
    @printf(" %016x", parm_id)
  end
  println()
  print("      coeff_modulus primes: ")
  for prime in coeff_modulus(parms(context_data))
    @printf("%016x ", value(prime))
  end
  println()
  println("\\")
  print(" \\-->")

  context_data = first_context_data(context)
  while !isnothing(context_data)
    print(" Level (chain index): ", chain_index(context_data))
    if parms_id(context_data) == first_parms_id(context)
      println(" ...... first_context_data()")
    elseif parms_id(context_data) == last_parms_id(context)
      println(" ...... last_context_data()")
    else
      println()
    end
    print("      parms_id:")
    for parm_id in parms_id(context_data)
      @printf(" %016x", parm_id)
    end
    println()
    print("      coeff_modulus primes: ")
    for prime in coeff_modulus(parms(context_data))
      @printf("%016x ", value(prime))
    end
    println()
    println("\\")
    print(" \\-->")

    context_data = next_context_data(context_data)
  end
  println(" End of chain reached")

  keygen = KeyGenerator(context)
  public_key_ = PublicKey()
  create_public_key!(public_key_, keygen)
  secret_key_ = secret_key(keygen)
  relin_keys_ = RelinKeys()
  create_relin_keys!(relin_keys_, keygen)

  print_line(@__LINE__)
  println("Print the parameter IDs of generated elements.")
  print("    + public_key:  ")
  for parm_id in parms_id(public_key_)
    @printf(" %016x", parm_id)
  end
  println()
  print("    + secret_key:  ")
  for parm_id in parms_id(secret_key_)
    @printf(" %016x", parm_id)
  end
  println()
  print("    + relin_keys:  ")
  for parm_id in parms_id(relin_keys_)
    @printf(" %016x", parm_id)
  end
  println()

  encryptor = Encryptor(context, public_key_)
  evaluator = Evaluator(context)
  decryptor = Decryptor(context, secret_key_)

  plain = Plaintext("1x^3 + 2x^2 + 3x^1 + 4")
  encrypted = Ciphertext()
  encrypt!(encrypted, plain, encryptor)
  print("    + plain:       ")
  for parm_id in parms_id(plain)
    @printf(" %016x", parm_id)
  end
  println(" (not set in BFV)")
  print("    + encrypted:   ")
  for parm_id in parms_id(encrypted)
    @printf(" %016x", parm_id)
  end
  println()
  println()

  print_line(@__LINE__)
  println("Perform modulus switching on encrypted and print.")
  context_data = first_context_data(context)
  print("---->")
  while !isnothing(next_context_data(context_data))
    println(" Level (chain index): ", chain_index(context_data))
    print("      parms_id of encrypted: ")
    for parm_id in parms_id(encrypted)
      @printf(" %016x", parm_id)
    end
    println()
    println("      Noise budget at this level: ", invariant_noise_budget(encrypted, decryptor),
            " bits")
    println("\\")
    print(" \\-->")
    mod_switch_to_next_inplace!(encrypted, evaluator)
    context_data = next_context_data(context_data)
  end
  println(" Level (chain index): ", chain_index(context_data))
  print("      parms_id of encrypted: ")
  for parm_id in parms_id(encrypted)
    @printf(" %016x", parm_id)
  end
  println()
  println("      Noise budget at this level: ", invariant_noise_budget(encrypted, decryptor),
          " bits")
  println("\\")
  print(" \\-->")
  println(" End of chain reached")
  println()

  print_line(@__LINE__)
  println("Decrypt still works after modulus switching.")
  decrypt!(plain, encrypted, decryptor)
  print("    + Decryption of encrypted: ", to_string(plain))
  println(" ...... Correct.")
  println()

  println("Computation is more efficient with modulus switching.")
  print_line(@__LINE__)
  println("Compute the 8th power.")
  encrypt!(encrypted, plain, encryptor)
  println("    + Noise budget fresh:                   ",
          invariant_noise_budget(encrypted, decryptor),
          " bits")
  square_inplace!(encrypted, evaluator)
  relinearize_inplace!(encrypted, relin_keys_, evaluator)
  println("    + Noise budget of the 2nd power:         ",
          invariant_noise_budget(encrypted, decryptor),
          " bits")
  square_inplace!(encrypted, evaluator)
  relinearize_inplace!(encrypted, relin_keys_, evaluator)
  println("    + Noise budget of the 4th power:         ",
          invariant_noise_budget(encrypted, decryptor),
          " bits")

  mod_switch_to_next_inplace!(encrypted, evaluator)
  println("    + Noise budget after modulus switching:  ",
          invariant_noise_budget(encrypted, decryptor),
          " bits")

  square_inplace!(encrypted, evaluator)
  relinearize_inplace!(encrypted, relin_keys_, evaluator)
  println("    + Noise budget of the 8th power:         ",
          invariant_noise_budget(encrypted, decryptor),
          " bits")
  mod_switch_to_next_inplace!(encrypted, evaluator)
  println("    + Noise budget after modulus switching:  ",
          invariant_noise_budget(encrypted, decryptor),
          " bits")

  decrypt!(plain, encrypted, decryptor)
  println("    + Decryption of the 8th power (hexadecimal) ...... Correct.")
  println("    ", to_string(plain))
  println()

  context = SEALContext(enc_parms, expand_mod_chain=false)

  println("Optionally disable modulus switching chain expansion.")
  print_line(@__LINE__)
  println("Print the modulus switching chain.")
  print("---->")
  context_data = key_context_data(context)
  while !isnothing(context_data)
    println(" Level (chain index): ", chain_index(context_data))
    print("      parms_id:")
    for parm_id in parms_id(context_data)
      @printf(" %016x", parm_id)
    end
    println()
    print("      coeff_modulus primes: ")
    for prime in coeff_modulus(parms(context_data))
      @printf("%016x ", value(prime))
    end
    println()
    println("\\")
    print(" \\-->")

    context_data = next_context_data(context_data)
  end
  println(" End of chain reached")
  println()

end


text:
module ReactiveMPModelsHGFTest

using Test, InteractiveUtils
using Rocket, ReactiveMP, GraphPPL, Distributions
using BenchmarkTools, Random, Plots, Dates, LinearAlgebra, StableRNGs

## Model definition
## -------------------------------------------- ##
@model [ default_factorisation = MeanField() ] function hgf_online_model(real_k, real_w, z_variance, y_variance)

    # Priors from previous time step for `z`
    zt_min_mean = datavar(Float64)
    zt_min_var  = datavar(Float64)

    # Priors from previous time step for `x`
    xt_min_mean = datavar(Float64)
    xt_min_var  = datavar(Float64)

    zt_min ~ NormalMeanVariance(zt_min_mean, zt_min_var)
    xt_min ~ NormalMeanVariance(xt_min_mean, xt_min_var)

    meta = GCVMetadata(GaussHermiteCubature(9))

    # Higher layer is modelled as a random walk
    zt ~ NormalMeanVariance(zt_min, z_variance) where { q = q(zt, zt_min)q(z_variance), meta = meta }

    # Lower layer is modelled with `GCV` node
    gcv_node, xt ~ GCV(xt_min, zt, real_k, real_w) where { q = q(xt, xt_min)q(zt)q(κ)q(ω) }

    # Noisy observations
    y = datavar(Float64)
    y ~ NormalMeanVariance(xt, y_variance)

    return zt, xt, y, gcv_node, xt_min_mean, xt_min_var, zt_min_mean, zt_min_var
end
## -------------------------------------------- ##
## Inference definition
## -------------------------------------------- ##
function hgf_online_inference(data, vmp_iters, real_k, real_w, z_variance, y_variance)
    n = length(data)

    # We don't want to save all marginals from all VMP iterations
    # but only last one after all VMP iterations per time step
    # Rocket.jl exports PendingScheduler() object that postpones
    # any update unless manual `resolve!()` has been called
    ms_scheduler = PendingScheduler()

    mz = keep(Marginal)
    mx = keep(Marginal)
    fe = ScoreActor(Float64)

    model, (zt, xt, y, gcv_node, xt_min_mean, xt_min_var, zt_min_mean, zt_min_var) = hgf_online_model(real_k, real_w, z_variance, y_variance)

    # Initial priors
    current_zt_mean, current_zt_var = 0.0, 10.0
    current_xt_mean, current_xt_var = 0.0, 10.0

    s_mz = subscribe!(getmarginal(zt) |> schedule_on(ms_scheduler), mz)
    s_mx = subscribe!(getmarginal(xt) |> schedule_on(ms_scheduler), mx)
    s_fe = subscribe!(score(Float64, BetheFreeEnergy(), model), fe)

    # Initial marginals to start VMP procedire
    setmarginal!(gcv_node, :y_x, MvNormalMeanCovariance([ 0.0, 0.0 ], [ 5.0, 5.0 ]))
    setmarginal!(gcv_node, :z, NormalMeanVariance(0.0, 5.0))

    # For each observations we perofrm `vmp_iters` VMP iterations
    for i in 1:n

        for _ in 1:vmp_iters
            update!(y, data[i])
            update!(zt_min_mean, current_zt_mean)
            update!(zt_min_var, current_zt_var)
            update!(xt_min_mean, current_xt_mean)
            update!(xt_min_var, current_xt_var)
        end

        # After all VMP iterations we release! `PendingScheduler`
        # as well as release! `ScoreActor` to indicate new time step
        release!(ms_scheduler)
        release!(fe)

        current_zt_mean, current_zt_var = mean_var(last(mz))::Tuple{Float64, Float64}
        current_xt_mean, current_xt_var = mean_var(last(mx))::Tuple{Float64, Float64}
    end

    # It is important to unsubscribe at the end of the inference procedure
    unsubscribe!((s_mz, s_mx, s_fe))

    return map(getvalues, (mz, mx, fe))
end
## -------------------------------------------- ##

@testset "Hierarchical Gaussian Filter" begin

    @testset "Online inference" begin
        ## -------------------------------------------- ##
        ## Data creation
        ## -------------------------------------------- ##
        function generate_data(rng, k, w, zv, yv)
            z_prev = 0.0
            x_prev = 0.0

            z = Vector{Float64}(undef, n)
            v = Vector{Float64}(undef, n)
            x = Vector{Float64}(undef, n)
            y = Vector{Float64}(undef, n)

            for i in 1:n
                z[i] = rand(rng, Normal(z_prev, sqrt(zv)))
                v[i] = exp(k * z[i] + w)
                x[i] = rand(rng, Normal(x_prev, sqrt(v[i])))
                y[i] = rand(rng, Normal(x[i], sqrt(yv)))

                z_prev = z[i]
                x_prev = x[i]
            end

            return z, x, y
        end
        ## -------------------------------------------- ##
        rng = StableRNG(42)
        # Parameters of HGF process
        real_k = 1.0
        real_w = 0.0
        z_variance = abs2(0.2)
        y_variance = abs2(0.1)
        # Number of observations
        n = 2000
        z, x, y = generate_data(rng, real_k, real_w, z_variance, y_variance);
        ## -------------------------------------------- ##
        ## Inference execution
        vmp_iters = 10
        mz, mx, fe = hgf_online_inference(y, vmp_iters, real_k, real_w, z_variance, y_variance);
        ## -------------------------------------------- ##
        ## Test inference results
        @test length(mz) === n
        @test length(mx) === n
        @test length(fe) === vmp_iters
        @test abs(last(fe) - 2027.8628798126442) < 0.01
        @test all(filter(e -> abs(e) > 0.1, diff(fe)) .< 0)
        # Check if all estimates are within 6std interval
        @test all((mean.(mz) .- 6 .* std.(mz)) .< z .< (mean.(mz) .+ 6 .* std.(mz)))
        @test all((mean.(mx) .- 6 .* std.(mx)) .< x .< (mean.(mx) .+ 6 .* std.(mx)))
        # Check if more than 95% of estimates are within 3std interval
        @test (sum((mean.(mz) .- 3 .* std.(mz)) .< z .< (mean.(mz) .+ 3 .* std.(mz))) / n) > 0.95
        @test (sum((mean.(mx) .- 3 .* std.(mx)) .< x .< (mean.(mx) .+ 3 .* std.(mx))) / n) > 0.95
        @test all(var.(mx) .> 0.0)
        @test all(var.(mz) .> 0.0)
        ## -------------------------------------------- ##
        ## Form debug output
        base_output = joinpath(pwd(), "_output", "models")
        mkpath(base_output)
        timestamp        = Dates.format(now(), "dd-mm-yyyy-HH-MM")
        plot_output      = joinpath(base_output, "hgf_online_plot_$(timestamp)_v$(VERSION).png")
        benchmark_output = joinpath(base_output, "hgf_online_benchmark_$(timestamp)_v$(VERSION).txt")
        ## -------------------------------------------- ##
        ## Create output plots
        pz = plot(title = "Hidden States Z")
        px = plot(title = "Hidden States X")

        plot!(pz, 1:n, z, label = "z_i", color = :orange)
        plot!(pz, 1:n, mean.(mz), ribbon = std.(mz), label = "estimated z_i", color = :teal)

        plot!(px, 1:n, x, label = "x_i", color = :green)
        plot!(px, 1:n, mean.(mx), ribbon = std.(mx), label = "estimated x_i", color = :violet)

        pf = plot(fe, label = "Bethe Free Energy")

        p = plot(pz, px, pf, layout = @layout([ a; b; c ]))
        savefig(p, plot_output)
        ## -------------------------------------------- ##
        ## Create output benchmarks
        benchmark = @benchmark hgf_online_inference($y, $vmp_iters, $real_k, $real_w, $z_variance, $y_variance) seconds=15
        open(benchmark_output, "w") do io
            show(io, MIME("text/plain"), benchmark)
            versioninfo(io)
        end
        ## -------------------------------------------- ##
    end

end

end


text:
"""
```
newuser(dsn, username, name; email="", phone="", office="")
```

This function creates a new user in the `users` table of the database.
E-mail, phone, office are useful for contacting users about their measurements.
"""
function newuser(dsn, username, name; email="", phone="", office="")
    first, middle, last = "", "", ""
    names = split(name, " ")
    if length(names) == 1
        first = names[1]
    elseif length(names) == 2
        first, last == (names...)
    elseif length(names) == 3
        first, mid, last == (names...)
    else
        error("unexpected name convention.")
    end

    ph = replace(phone, r"[-\+\(\)\s]", "")
    ODBC.execute!(dsn, """
    INSERT INTO users VALUES (
        '$username', '$first', '$middle', '$last', '$email', '$phone', '$office'
    );
    """)
end

"""
```
updateuser(dsn, username; kwargs...)
```

Update an existing user in the `users` table, identified by `username`. Specify
the fields to update with keyword arguments specified in
[`ICDataServer.newuser`](@ref).
"""
function updateuser(dsn, username; kwargs...)
    isempty(kwargs) && return
    pstr = reduce((a,b)->a*","*b, "$k = '$v'" for (k,v) in kwargs)
    for (k,v) in kwargs
        ODBC.execute!(dsn, "UPDATE users SET "*pstr*" WHERE username = $username;")
    end
end

"""
```
deleteuser(dsn, username)
```

Delete a user from the `users` table by providing the username.
"""
function deleteuser(dsn, username)
    ODBC.execute!(dsn, "DELETE FROM users WHERE username = '$username';")
end

"""
```
listusers(dsn)
```

List all users in the `users` table.
"""
function listusers(dsn)
    ODBC.query(dsn, "SELECT * FROM users;")
end


text:
# VCF Header
# ==========

immutable VCFHeader
    metainfo::Vector{VCFMetaInfo}
    sampleID::Vector{String}

    function VCFHeader(metainfo::Vector{VCFMetaInfo}, sampleID::Vector{String})
        return new(metainfo, sampleID)
    end
end

"""
    VCFHeader()

Create an empty `VCFHeader` object.
"""
function VCFHeader()
    return VCFHeader(VCFMetaInfo[], String[])
end

"""
    VCFHeader(metainfo::Vector, sampleID::Vector)

Create a `VCFHeader` object with `metainfo` and `sampleID`.
"""
function VCFHeader(metainfo::Vector, sampleID::Vector)
    return VCFHeader(convert(Vector{VCFMetaInfo}, metainfo), convert(Vector{String}, sampleID))
end

function Base.eltype(::Type{VCFHeader})
    return VCFMetaInfo
end

function Base.length(header::VCFHeader)
    return length(header.metainfo)
end

function Base.start(header::VCFHeader)
    return 1
end

function Base.done(header::VCFHeader, i)
    return i > endof(header.metainfo)
end

function Base.next(header::VCFHeader, i)
    return header.metainfo[i], i + 1
end

function Base.find(header::VCFHeader, tag::AbstractString)
    return Base.filter(m -> isequaltag(m, tag), header.metainfo)
end

function Base.unshift!(header::VCFHeader, metainfo)
    unshift!(header.metainfo, convert(VCFMetaInfo, metainfo))
    return header
end

function Base.push!(header::VCFHeader, metainfo)
    push!(header.metainfo, convert(VCFMetaInfo, metainfo))
    return header
end

function Base.show(io::IO, header::VCFHeader)
    println(io, summary(header), ':')
    tags = metainfotag.(header.metainfo)
    println(io, "  metainfo tags: ", join(unique(tags), ' '))
      print(io, "     sample IDs: ", join(header.sampleID, ' '))
end

function Base.write(io::IO, header::VCFHeader)
    n = 0
    for metainfo in header.metainfo
        n += write(io, metainfo, '\n')
    end
    n += write(io, "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO")
    if !isempty(header.sampleID)
        n += write(io, '\t', "FORMAT")
        for id in header.sampleID
            n += write(io, '\t', id)
        end
    end
    n += write(io, '\n')
    return n
end


text:
###############################################################################
##### Training Data type and constructors
###############################################################################

struct TrainData
  n::Int
  p::Int
  xmin::Matrix{Float64}
  xmax::Matrix{Float64}
  # dt::ZScoreTransform{Float64}
  X::Matrix{Float64}
  ybar::Float64
  y::AbstractVector
  σhat::Float64
end

function TrainData(X::Matrix{Float64}, y::Vector{Float64})
  n = length(y)
  p = size(X, 2)
  Q = Matrix(qr(X).Q)
  ybar = mean(y)
  y = y .- ybar
  yhat = Q * Q' * y
  r = y - yhat
  σhat = p >= n ? std(y) : sqrt(dot(r, r) / (n - p))
  # dt = fit(ZScoreTransform, transpose(X))
  # X = StatsBase.transform(dt, transpose(X))
  # X = Matrix(transpose(X))
  xmin = minimum(X, dims = 1)
  xmax = maximum(X, dims = 1)
  TrainData(n, p, xmin, xmax, X, ybar, y, σhat)
end

function TrainData(X::Matrix{Float64}, y::Vector{Int})
  n = length(y)
  p = size(X, 2)
  # dt = fit(ZScoreTransform, transpose(X))
  # X = StatsBase.transform(dt, transpose(X))
  # X = Matrix(transpose(X))
  xmin = minimum(X, dims = 1)
  xmax = maximum(X, dims = 1)
  ybar = mean(y)
  TrainData(n, p, xmin, xmax, X, ybar, y, 1.0)
end

###############################################################################
##### MCMC options type
###############################################################################

struct Opts
  nchains::Int
  nburn::Int
  ndraw::Int
  nthin::Int
  S::Int64
  function Opts(;nchains = 4, nburn = 2500, ndraw = 2500, nthin = 1)
    new(nchains, nburn, ndraw, nthin, nburn + ndraw)
  end
end

###############################################################################
##### BART hyperparameters type
###############################################################################

struct Hypers
  m::Int64
  k::Int64
  ν::Float64
  δ::Float64
  q::Float64
  α::Float64
  β::Float64
  hard::Bool
  λmean::Float64
  λfix::Bool
  sigma_noninf::Bool
  τ::Float64
  init_trees::String
  init_depth::Vector
  sparse::Bool
  shape::Float64
  a::Float64
  b::Float64
  # group_idx::Vector{Int}
  # groups::Vector{Int}
  # group_size::Vector{Int}
  # scale::Int
  function Hypers(td::TrainData; m = 50, k = 2,
    ν = 3.0, q = 0.9, α = 0.95, β = 2.0,
    sigma_noninf = true,
    hard = true, λmean = 0.1, λfix = false,
    init_trees = "leaf", init_depth = ones(4),
    sparse = false, shape = 1.0, a = 0.5, b = 1.0
    # ,
    # group_idx = nothing
    )
    if !(init_trees in ["leaf", "rf"])
      throw(ArgumentError("init_tree options are leaf or rf"))
    end
    if sigma_noninf
      ν = 0
      δ = 0
    else
      δ = 1 / quantile(InverseGamma(ν / 2, ν / (2 * td.σhat^2)), q)
    end
    if isa(td.y, Vector{Int})
      τ = (3.0 / (k*sqrt(m)))^2
    else
      τ = ((maximum(td.y) - minimum(td.y)) / (2*k*sqrt(m)))^2
    end
    if hard
      λfix = true
      λmean = 1/10000
    end
    # group_idx = isa(group_idx, Nothing) ? collect(1:td.p) : group_idx
    # groups = unique(group_idx)
    # group_size = counts(group_idx)
    # scale = length(group_size)
    new(m, k, ν, δ, q, α, β, hard, λmean, λfix, sigma_noninf, τ,
      init_trees, init_depth, sparse, shape, a, b
      # ,
      # group_idx, groups, group_size, scale
      )
  end
end

###############################################################################
##### Bart model type and constructors
###############################################################################

struct BartModel
  hypers::Hypers
  opts::Opts
  td::TrainData
end

function BartModel(X::Matrix{Float64}, y::AbstractVector, opts::Opts; hyperargs...)
  td = TrainData(X, y)
  hypers = Hypers(td; hyperargs...)
  BartModel(hypers, opts, td)
end

###############################################################################
##### State of sampler
###############################################################################

abstract type BartState end

mutable struct SuffStats
  Lt::Int
  Ω::Matrix{Float64}
  rhat::Vector{Float64}
end

function suffstats(rt::Vector{Float64}, S::Matrix{Float64}, bs::BartState, bm::BartModel)
  Lt = size(S, 2)
  Ω = inv(transpose(S) * S / bs.σ^2 + I / bm.hypers.τ)
  rhat = transpose(S) * rt / bs.σ^2
  SuffStats(Lt, Ω, rhat)
end

mutable struct BartTree
  tree::Tree
  S::Matrix{Float64}
  ss::SuffStats
end

mutable struct BartEnsemble
  trees::Vector{BartTree}
end

mutable struct RegBartState <: BartState
  ensemble::BartEnsemble
  fhat::Vector{Float64}
  σ::Float64
  s::Vector{Float64}
  shape::Float64
end

function Base.convert(Node, x)
  isa(x, DecisionTree.Leaf) ? Leaf(x.majority) :
    Branch(x.featid, x.featval, convert(Node, x.left), convert(Node, x.right))
end

function RegBartState(bm::BartModel)
  states = RegBartState[]
  for c in 1:bm.opts.nchains
    if bm.hypers.init_trees == "leaf"
      trees = [Tree(Leaf(0.0), bm.hypers.λmean) for _ in 1:bm.hypers.m]
    else
      rf = DecisionTree.fit!(
        DecisionTree.RandomForestRegressor(
          n_trees = bm.hypers.m, max_depth = bm.hypers.init_depth[c],
          min_purity_increase = 5),
          bm.td.X, bm.td.y
      )
      trees = Tree.(convert.(Node, rf.ensemble.trees), bm.hypers.λmean)
    end
    S = [leafprob(bm.td.X, tree) for tree in trees]
    fhats = reduce(hcat, [S[t] * getμ(trees[t]) for t in eachindex(trees)])
    yhat = vec(sum(fhats, dims = 2))
    bt = Vector{BartTree}(undef, bm.hypers.m)
    for t in eachindex(trees)
      rt = bm.td.y - sum(fhats[:,eachindex(trees) .!= t], dims = 2)
      Ω = inv(transpose(S[t]) * S[t] / bm.td.σhat^2 + I / bm.hypers.τ)
      rhat = vec(transpose(S[t]) * rt / bm.td.σhat^2)
      bt[t] = BartTree(trees[t], S[t], SuffStats(size(S[t], 2), Ω, rhat))
    end
    push!(states,
      RegBartState(BartEnsemble(bt), yhat, bm.td.σhat, ones(bm.td.p) ./ bm.td.p, bm.hypers.shape))
  end
  states
end

mutable struct ProbitBartState <: BartState
  ensemble::BartEnsemble
  fhat::Vector{Float64}
  z::Vector{Float64}
  σ::Float64
  s::Vector{Float64}
  shape::Float64
end

function ProbitBartState(bm::BartModel)
  states = []
  z = map(y -> y == 1 ?
    rand(truncated(Normal(), 0, Inf)) :
    rand(truncated(Normal(), -Inf, 0)),
    bm.td.y
  )
  for c in 1:bm.opts.nchains
    if bm.hypers.init_trees == "leaf"
      trees = [Tree(Leaf(0.0), bm.hypers.λmean) for _ in 1:bm.hypers.m]
    else
      rf = DecisionTree.fit!(
        DecisionTree.RandomForestRegressor(
          n_trees = bm.hypers.m, max_depth = bm.hypers.init_depth[c],
          min_purity_increase = 5),
          bm.td.X, z
      )
      trees = Tree.(convert.(Node, rf.ensemble.trees), bm.hypers.λmean)
    end
    S = [leafprob(bm.td.X, tree) for tree in trees]
    fhats = reduce(hcat, [S[t] * getμ(trees[t]) for t in eachindex(trees)])
    yhat = vec(sum(fhats, dims = 2))
    bt = BartEnsemble(Vector{BartTree}(undef, bm.hypers.m))
    for t in eachindex(trees)
      rt = z - sum(fhats[:,eachindex(trees) .!= t], dims = 2)
      Ω = inv(transpose(S[t]) * S[t] + I / bm.hypers.τ)
      rhat = vec(transpose(S[t]) * rt)
      bt.trees[t] = BartTree(trees[t], S[t], SuffStats(size(S[t], 2), Ω, rhat))
    end
    push!(states,
      ProbitBartState(bt, yhat, z, 1, ones(bm.td.p) ./ bm.td.p, bm.hypers.shape))
  end
  states
end

###############################################################################
##### Posterior draws from BART model
###############################################################################

abstract type BartPosterior end

struct RegBartPosterior <: BartPosterior
  sdraws::Matrix{Float64}
  σdraws::Vector{Float64}
  treedraws::Vector{Vector{Tree}}
  function RegBartPosterior(bm::BartModel)
    new(
      Matrix{Float64}(undef, bm.td.p, bm.opts.ndraw),
      Vector{Float64}(undef, bm.opts.ndraw),
      Vector{Vector{Tree}}(undef, bm.opts.ndraw)
    )
  end
end

struct ProbitBartPosterior <: BartPosterior
  sdraws::Matrix{Float64}
  zdraws::Matrix{Float64}
  treedraws::Vector{Vector{Tree}}
  function ProbitBartPosterior(bm::BartModel)
    new(
      Matrix{Float64}(undef, bm.td.p, bm.opts.ndraw),
      Matrix{Float64}(undef, bm.td.n, bm.opts.ndraw),
      Vector{Vector{Tree}}(undef, bm.opts.ndraw)
    )
  end
end

abstract type BartChain end

struct RegBartChain <: BartChain
  bm::BartModel
  init_trees
  sdraws::Array{Float64}
  treedraws::Array{Vector{Tree}}
  σdraws::Array{Float64}
end

struct ProbitBartChain <: BartChain
  bm::BartModel
  init_trees
  sdraws::Array{Float64}
  treedraws::Array{Vector{Tree}}
end

struct SurvBartChain <: BartChain
  bm::BartModel
  time::Vector{Float64}
  event::Vector{Int}
  init_trees
  sdraws::Array{Float64}
  treedraws::Array{Vector{Tree}}
  Sdraws::Array{Float64}
end

###############################################################################
##### MCMCChains
###############################################################################

function MCMCChains.Chains(bc::Union{ProbitBartChain,SurvBartChain})
  atd = mapslices(td -> map(x -> mean(depth.(x)), td), bc.treedraws, dims = 3)
  mpost = reshape(ptd(bc), bc.bm.opts.ndraw, 1, bc.bm.opts.nchains)
  monitor = hcat(mpost, atd)
  Chains(mpost, ["PTD", "ATD"])
end

function MCMCChains.Chains(bc::RegBartChain)
  atd = mapslices(td -> map(x -> mean(depth.(x)), td), bc.treedraws, dims = 3)
  mpost = reshape(ptd(bc), bc.bm.opts.ndraw, 1, bc.bm.opts.nchains)
  monitor = hcat(bc.σdraws, mpost, atd)
  Chains(monitor, ["sigma", "PTD", "ATD"])
end


text:
using Satisfactory
using Test

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


text:
# Exploration of an apparent bug in array iteration
# https://discourse.julialang.org/t/last-row-is-broken-when-iterating-an-array/15696

using DelimitedFiles

mutable struct mytbl
    source
    function mytbl(csvdata::Base.GenericIOBuffer)
        source = readdlm(csvdata, ',')
        source = convert(Array, source[2:end,:]) # clear the headers
        source = reshape(source, size(source, 1), size(source, 2))
        new(source)
    end
end

Base.length(it::mytbl) = size(it.source, 1)
function Base.iterate(it::mytbl, (el, i)=(nothing, 1))
    if i > length(it); return nothing; end
    return (it.source[i,:], (nothing, i + 1))
end

TABLE_CAST = """id,height,age,name,occupation
1,10.0,1,string1,2012-06-15 00:00:00
2,10.1,2,string2,2013-06-15 01:00:00
3,10.2,3,string3,2014-06-15 02:00:00
4,10.3,4,string4,2015-06-15 03:00:00
5,10.4,5,string5,2016-06-15 04:00:00
"""

table =  mytbl(IOBuffer(TABLE_CAST))

[ row for row in table ]
# the last element is #undef

[ row[1] for row in table ]
# Returns:
[1, 2, 3, 4, 229445824]

[ row[2] for row in table ]
# Returns:
[10.0, 10.1, 10.2, 10.3, 4.94e-324]

max_stars_repo_name: loleg/tableschema-jl

max_issues_repo_name: loleg/tableschema-jl

max_forks_repo_name: loleg/tableschema-jl


text:
using Test
using CyclicArrays

#=
face: 1:nfaces
direction:X=1 / Y=2 / Z=3,
side: - = 1 / + = 2
destiny:
1  face: 1:nfaces
2  direction:X=1 / Y=2 / Z=3,
3  side: - = 1 / + = 2
4  flip 1 / no flip 0
=#

# 3d array example

nx=4
ny=4
nz=1
nfaces=1;

faces=zeros(nfaces,3,2,4);
faces[1,1,1,:]=[1,1,2,0];
faces[1,1,2,:]=[1,1,1,0];
faces[1,2,1,:]=[1,2,2,0];
faces[1,2,2,:]=[1,2,1,0];
faces[1,3,1,:].=-1;
faces[1,3,2,:].=-1;

grid=CyclicArray(faces);

x0=reshape(collect(0:15).%4 .+1,(4,4))

x=reshape(x0,(1,4,4));
y=reshape(x0',(1,4,4));

x1=CyclicArray(x,grid)
y1=CyclicArray(y,grid)

i=5;j=1;k=1;
x1[k,j,i],y1[k,j,i]

x1[k,j,i]=2

@test x1[1,5,1]==2
@test [x1[1,2,2],y1[1,2,2]]==[2,2]
@test [x1[1,3,2],y1[1,3,2]]==[3,2]
@test [x1[1,5,2],y1[1,5,2]]==[1,2]
@test [x1[1,2,5],y1[1,2,5]]==[2,1]

# 2d array example

nx=5
ny=5
nfaces=1;

faces=zeros(nfaces,2,2,4);
faces[1,1,1,:]=[1,1,2,0];
faces[1,1,2,:]=[1,1,1,0];
faces[1,2,1,:]=[1,2,2,0];
faces[1,2,2,:]=[1,2,1,0];

grid=CyclicArray(faces);

x=reshape(collect(0:15).%4 .+1,(4,4));
y=x'

x1=CyclicArray(x,grid)
y1=CyclicArray(y,grid)

i=5;j=1;
x1[i,j],y1[i,j]

x1[i,j]=2

@test x1[5,1]==2
@test [x1[2,2],y1[2,2]]==[2,2]
@test [x1[3,2],y1[3,2]]==[3,2]
@test [x1[5,2],y1[5,2]]==[1,2]
@test [x1[2,5],y1[2,5]]==[2,1]

# 1d array example

faces=zeros(1,1,2,4);
faces[1,1,1,:]=[1,1,2,0];
faces[1,1,2,:]=[1,1,1,0];

grid=CyclicArray(faces);
x=[0:3;]
x1=CyclicArray(x,grid)

@test x1[5]==0
@test x1[0]==3


text:
mutable struct BSONArray
    _wrap_::VoidPtr
    _ref_::Any

    function BSONArray()
        _wrap_ = ccall((:bson_new, libbson), VoidPtr, ())
        bsonArray = new(_wrap_, Union{})
        finalizer(destroy, bsonArray)
        bsonArray
    end

    function BSONArray(vec::Vector)
        bsonArray = BSONArray()
        for ele in vec
            append(bsonArray, ele)
        end
        bsonArray
    end

    function BSONArray(data::BytePtr, length::Integer, _ref_::Any)
        buffer = alloc_buf(128)
        ccall((:bson_init_static, libbson), Bool, (VoidPtr, BytePtr, UInt32),
              buffer, data, length) || error("bson_init_static: failure")
        b = Base.unsafe_convert(VoidPtr, buffer)
        new(b, (_ref_, buffer))
    end

    BSONArray(_wrap_::VoidPtr, _ref_::Any) = new(_wrap_, _ref_)
end
export BSONArray

@static if VERSION < v"0.7-"
Base.iteratoreltype(::Type{BSONArray}) = Base.EltypeUnknown()
else
Base.IteratorEltype(::Type{BSONArray}) = Base.EltypeUnknown()
end

function Base.convert(::Type{AbstractString}, bsonArray::BSONArray)
    cstr = ccall((:bson_array_as_json, libbson), BytePtr, (VoidPtr, BytePtr),
                 bsonArray._wrap_, C_NULL)
    cstr == C_NULL && return ""
    result = unsafe_string(cstr)
    ccall((:bson_free, libbson), Cvoid, (VoidPtr,), cstr)
    result
end

Base.convert(::Type{Array}, b::BSONArray) = collect(b)
Base.convert(::Type{Array{T}}, b::BSONArray) where {T} = collect(T, b)
Base.convert(::Type{Vector{T}}, b::BSONArray) where {T} = collect(T, b)

Base.string(bsonArray::BSONArray) = convert(AbstractString, bsonArray)

Base.show(io::IO, bsonArray::BSONArray) =
    print(io, "BSONArray($(convert(AbstractString, bsonArray)))")

Base.length(bsonArray::BSONArray) =
    ccall((:bson_count_keys, libbson), UInt32, (VoidPtr,), bsonArray._wrap_)

function append(bsonArray::BSONArray, val::Bool)
    keyCStr = string(length(bsonArray))
    ccall((:bson_append_bool, libbson), Bool, (VoidPtr, BytePtr, Cint, Bool),
          bsonArray._wrap_, keyCStr, length(keyCStr), val) || bson_ovf()
end

function append(bsonArray::BSONArray, val::Real)
    keyCStr = string(length(bsonArray))
    ccall((:bson_append_double, libbson), Bool, (VoidPtr, BytePtr, Cint, Cdouble),
          bsonArray._wrap_, keyCStr, length(keyCStr), val) || bson_ovf()
end

function append(bsonArray::BSONArray, val::DateTime)
    keyCStr = string(length(bsonArray))
    ts = round(Int64, datetime2unix(val)*1000)
    ccall((:bson_append_date_time, libbson), Bool, (VoidPtr, BytePtr, Cint, Clonglong),
          bsonArray._wrap_, keyCStr, length(keyCStr), ts) || bson_ovf()
end

append(bsonArray::BSONArray, val::Date) = append(bsonArray, DateTime(val))

function append(bsonArray::BSONArray, val::BSONArray)
    keyCStr = string(length(bsonArray))
    childBuffer = alloc_buf(128)
    ccall((:bson_append_array_begin, libbson), Bool, (VoidPtr, BytePtr, Cint, VoidPtr),
          bsonArray._wrap_, keyCStr, length(keyCStr), childBuffer) ||
              error("bson_append_array_begin: failure")
    childBSONArray = BSONArray(Base.unsafe_convert(VoidPtr, childBuffer), childBuffer)
    for element in val
        append(childBSONArray, element)
    end
    ccall((:bson_append_array_end, libbson), Bool, (VoidPtr, VoidPtr),
          bsonArray._wrap_, childBuffer) || error("bson_append_array_end: failure")
end

function append(bsonArray::BSONArray, val::Union{Int8, UInt8, Int16, UInt16, Int32, UInt32})
    keyCStr = string(length(bsonArray))
    ccall((:bson_append_int32, libbson), Bool, (VoidPtr, BytePtr, Cint, Int32),
          bsonArray._wrap_, keyCStr, length(keyCStr), val) || bson_ovf()
end

function append(bsonArray::BSONArray, val::Union{Int64, UInt64})
    keyCStr = string(length(bsonArray))
    ccall((:bson_append_int64, libbson), Bool, (VoidPtr, BytePtr, Cint, Int64),
          bsonArray._wrap_, keyCStr, length(keyCStr), val) || bson_ovf()
end

function append(bsonArray::BSONArray, val::BSONOID)
    keyCStr = string(length(bsonArray))
    ccall((:bson_append_oid, libbson), Bool, (VoidPtr, BytePtr, Cint, BytePtr),
          bsonArray._wrap_, keyCStr, length(keyCStr), val._wrap_) || bson_ovf()
end

function append(bsonArray::BSONArray, val::AbstractString)
    keyCStr = string(length(bsonArray))
    valUTF8 = String(val)
    ccall((:bson_append_utf8, libbson), Bool, (VoidPtr, BytePtr, Cint, BytePtr, Cint),
          bsonArray._wrap_, keyCStr, length(keyCStr), valUTF8, sizeof(valUTF8)) || bson_ovf()
end

append(bsonArray::BSONArray, val::Cvoid) = append_null(bsonArray)

function append(bsonArray::BSONArray, val::Symbol)
    if val == :null
        append_null(bsonArray)
    elseif val == :minkey
        append_minkey(bsonArray)
    elseif val == :maxkey
        append_maxkey(bsonArray)
    else
        append(bsonArray, string(val))
    end
end

function append(bsonArray::BSONArray, val::AbstractDict)
    keyCStr = string(length(bsonArray))
    childBuffer = alloc_buf(128)
    ccall((:bson_append_document_begin, libbson), Bool, (VoidPtr, BytePtr, Cint, VoidPtr),
          bsonArray._wrap_, keyCStr, length(keyCStr), childBuffer) ||
              error("bson_append_document_begin: failure")
    childBSONObject = BSONObject(Base.unsafe_convert(VoidPtr, childBuffer), childBuffer)
    for (k, v) in val
        append(childBSONObject, k, v)
    end
    ccall((:bson_append_document_end, libbson), Bool, (VoidPtr, VoidPtr),
          bsonArray._wrap_, childBuffer) || error("bson_append_document_end: failure")
end

function append(bsonArray::BSONArray, val::Vector)
    keyCStr = string(length(bsonArray))
    childBuffer = alloc_buf(128)
    ccall((:bson_append_array_begin, libbson), Bool, (VoidPtr, BytePtr, Cint, VoidPtr),
          bsonArray._wrap_, keyCStr, length(keyCStr), childBuffer) ||
              error("bson_append_array_begin: failure")
    childBSONArray = BSONArray(Base.unsafe_convert(VoidPtr, childBuffer), childBuffer)
    for element in val
        append(childBSONArray, element)
    end
    ccall((:bson_append_array_end, libbson), Bool, (VoidPtr, VoidPtr),
          bsonArray._wrap_, childBuffer) || error("bson_append_array_end: failure")
end

export append

function append_null(bsonArray::BSONArray)
    keyCStr = string(length(bsonArray))
    ccall((:bson_append_null, libbson), Bool, (VoidPtr, BytePtr, Cint),
          bsonArray._wrap_, keyCStr, length(keyCStr)) || bson_ovf()
end
export append_null

function append_minkey(bsonArray::BSONArray)
    keyCStr = string(length(bsonArray))
    ccall((:bson_append_minkey, libbson), Bool, (VoidPtr, BytePtr, Cint),
          bsonArray._wrap_, keyCStr, length(keyCStr)) || bson_ovf()
end
export append_minkey

function append_maxkey(bsonArray::BSONArray)
    keyCStr = string(length(bsonArray))
    ccall((:bson_append_maxkey, libbson), Bool, (VoidPtr, BytePtr, Cint),
          bsonArray._wrap_, keyCStr, length(keyCStr)) || bson_ovf()
end
export append_maxkey

function vector(bsonArray::BSONArray)
    vec = Any[]
    for v in bsonArray
        if isa(v, BSONObject)
            push!(vec, dict(v))
        elseif isa(v, BSONArray)
            push!(vec, vector(v))
        else
            push!(vec, v)
        end
    end
    vec
end
export vector

# Private

destroy(bsonArray::BSONArray) =
    ccall((:bson_destroy, libbson), Cvoid, (VoidPtr,), bsonArray._wrap_)


text:
import .Unitful

tweight(A::AbstractArray{T}) where T <: Unitful.Quantity = tweight(Unitful.ustrip(A))


text:
using Documenter, ImageGather

makedocs(sitename="Image gather tools",
         doctest=false, clean=true,
         authors="Mathias Louboutin",
         pages = Any["Home" => "index.md"])

deploydocs(repo="github.com/slimgroup/ImageGather.jl",
           devbranch="main")


text:
Zero = [1, 0]
One = [0, 1]

I = [1 0
     0 1]

X = [0 1
     1 0]

Y = [0 -im
     im 0]

Z = [1  0
     0 -1]

H = [1  1
     1 -1]/√2

S = classical((a, b) -> (b, a), 2)

control(U) = Zero*Zero' ⊗ eye(U) + One*One' ⊗ U

CX = control(X)

T = control(CX)


text:
##############
#  Results   #
##############

struct Results
  samples::Int
  min_time::Float64
  timings::Array{Float64, 1}
end

function Results(n::Int, timings::Array{Float64, 1})
  samples = n
  min_time = minimum(timings)

  return Results(samples, min_time, timings)
end


text:
export init, ok, spin, spinOnce, now
export info, info_stream, warn, warn_stream, error, error_stream, debug, debug_stream, fatal, fatal_stream

"""
    now()

Retrieve the current time. If ROS clock time is in use, this returns the time according to the ROS clock. Otherwise returns the current wall clock time.
"""
now() = @cxx ros::Time::now()

"""
    spin()

Enter simple event loop.
"""
spin() = @cxx ros::spin()

"""
    spinOnce()

Process a single round of callbacks.
"""
spinOnce() = @cxx ros::spinOnce()

"""
    ok()

Check whether it's time to exit.
"""
ok() = @cxx ros::ok()

"""
    init(node_name::String)

ROS initialization function.
"""
init(node_name::String) = @cxx ros::init(Int32(length(ARGS)+1), pointer(pointer.(pushfirst!(ARGS,PROGRAM_FILE))), pointer(node_name))

"""
    info(n)

Print `n` with a ROS info logging level.
"""
info(n) = icxx"ROS_INFO_STREAM($n);"

"""
    info_stream(n)

Print `n` with a ROS info logging level.
"""
info_stream(n) = info(n)

"""
    warn(n)

Print `n` with a ROS warn logging level.
"""
warn(n) = icxx"ROS_WARN_STREAM($n);"

"""
    warn_stream(n)

Print `n` with a ROS warn logging level.
"""
warn_stream(n) = warn(n)

"""
    error(n)

Print `n` with a ROS error logging level.
"""
error(n) = icxx"ROS_ERROR_STREAM($n);"

"""
    error_stream(n)

Print `n` with a ROS error logging level.
"""
error_stream(n) = error(n)

"""
    debug(n)

Print `n` with a ROS debug logging level.
"""
debug(n) = icxx"ROS_DEBUG_STREAM($n);"

"""
    debug_stream(n)

Print `n` with a ROS debug logging level.
"""
debug_stream(n) = debug(n)

"""
    fatal(n)

Print `n` with a ROS fatal logging level.
"""
fatal(n) = icxx"ROS_FATAL_STREAM($n);"

"""
    fatal_stream(n)

Print `n` with a ROS fatal logging level.
"""
fatal_stream(n) = fatal(n)

function ros♥julia()
    try
        ros_root = normpath(joinpath(ENV["ROS_ROOT"], "../.."))
        ros_header_dir = ros_root * "include"
        ros_lib_dir = ros_root * "lib"
        addHeaderDir(ros_header_dir, kind = C_System)
        package_path = String.(split(ENV["ROS_PACKAGE_PATH"], ":"))
        for path in package_path
            path = normpath(path, "..")
            if path != ros_root
                addHeaderDir(joinpath(path, "devel", "include"), kind = C_System)
            end
        end
        # Base ROS
        cxxinclude(joinpath(ros_header_dir,"ros","ros.h"))
        cxxinclude(joinpath(ros_header_dir,"ros","package.h"))
        Libdl.dlopen(joinpath(ros_lib_dir, "libroscpp.so"), Libdl.RTLD_GLOBAL)
        Libdl.dlopen(joinpath(ros_lib_dir, "libroslib.so"), Libdl.RTLD_GLOBAL)
        # Actionlib
        cxxinclude(joinpath(ros_header_dir,"actionlib","server","simple_action_server.h"))
        cxxinclude(joinpath(ros_header_dir,"actionlib","client","simple_action_client.h"))
        cxxinclude(joinpath(ros_header_dir,"actionlib","client","simple_client_goal_state.h"))
        cxxinclude(joinpath(ros_header_dir,"actionlib","client","terminal_state.h"))
        Libdl.dlopen(joinpath(ros_lib_dir, "libactionlib.so"), Libdl.RTLD_GLOBAL)
        # tf
        cxxinclude(joinpath(ros_header_dir,"tf2_ros","transform_broadcaster.h"))
        cxxinclude(joinpath(ros_header_dir,"tf2_ros","static_transform_broadcaster.h"))
        cxxinclude(joinpath(ros_header_dir,"tf2_ros","transform_listener.h"))
        cxxinclude(joinpath(ros_header_dir,"tf2_geometry_msgs","tf2_geometry_msgs.h"))
        cxxinclude(joinpath(ros_header_dir, "tf2", "LinearMath", "Quaternion.h"))
        Libdl.dlopen(joinpath(ros_lib_dir, "libtf2.so"), Libdl.RTLD_GLOBAL)
        Libdl.dlopen(joinpath(ros_lib_dir, "libtf2_ros.so"), Libdl.RTLD_GLOBAL)
        #  println(pushfirst!(ARGS,PROGRAM_FILE))
    catch
        throw(ErrorException("ROS.jl cannot find your ROS installation! ¯\\_(ツ)_/¯"))
    end
end


text:
using Test

@testset "CitZAF" begin
    # Mostly the testing is performed in testagainstheinrich.jmd
    # These tests just ensure stablity.
    θ = deg2rad(40.0)
    algs = zafcorrection(CitZAF, ReedFluorescence, NullCoating, mat"NaAlSi3O8", mat"NaF", n"Na K", 15.0e3)
    @test isapprox(A(algs..., n"Na K-L3",θ,θ), 1.0404, atol=0.0001)
    @test isapprox(Z(algs...), 0.96463, atol=0.0001)
    algs = zafcorrection(CitZAF, ReedFluorescence, NullCoating, mat"NaAlSi3O8", mat"Al2O3", n"Al K", 15.0e3)
    @test isapprox(A(algs..., n"Al K-L3",θ,θ), 0.89712, atol=0.0001)
    @test isapprox(Z(algs...), 0.9896, atol=0.0001)
    algs = zafcorrection(CitZAF, ReedFluorescence, NullCoating, mat"NaAlSi3O8", mat"Al2O3", n"O K", 15.0e3)
    @test isapprox(A(algs..., n"O K-L3",θ,θ), 0.94331, atol=0.0001)
    @test isapprox(Z(algs...), 0.98962, atol=0.0001)
    algs = zafcorrection(CitZAF, ReedFluorescence, NullCoating, mat"NaAlSi3O8", mat"SiO2", n"Si K", 15.0e3)
    @test isapprox(A(algs..., n"Si K-L3",θ,θ), 0.84532, atol=0.0001)
    @test isapprox(Z(algs...), 1.0065, atol=0.0001)
end

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
module TimeIt

if VERSION >= v"0.7.0"
    using Printf
end

export @timeit

macro timen(ex, n)
    quote
        local t0 = time_ns()
        for i = 1:$(esc(n))
            local val = $(esc(ex))
        end
        local t1 = time_ns()
        (t1 - t0) / 1.e9
    end
end

macro timeit(ex)
    quote
        local val = $(esc(ex))  # Warm up
        t = zeros(3)

        # Determine number of loops so that total time > 0.1s.
        n = 1
        for i = 0:9
            n = 10^i
            t[1] = @timen $(esc(ex)) n
            if t[1] > 0.1
                break
            end
        end

        # Two more production runs.
        for i = 2:3
            t[i] = @timen $(esc(ex)) n
        end
        best = minimum(t) / n

        # Format to nano-, micro- or milliseconds.
        factor = 1.
        if best < 1e-6
            factor = 1e9
            pre = "n"
        elseif best < 1e-3
            factor = 1e6
            pre = "\u00b5"
        elseif best < 1.
            factor = 1e3
            pre = "m"
        else
            factor = 1.
            pre = ""
        end
        @printf "%d loops, best of 3: %4.2f %ss per loop\n" n best*factor pre
        best
    end
end

end  # module


text:
module GraphFlux

using CUDA
using Flux
using Functors
using ChainRulesCore
using Setfield

using DifferentialEquations
using DiffEqFlux

export
    GraphTuple,
    BatchGraphTuple,
    batchgraphs,

    GCN,
    GCNₑ,
    NodeLayer,
    graphmeanpool,

    GraphNeuralODE

include("dtypes/graphtuple.jl")

include("ops/messageprop.jl")

# GPU ops
include("cuda/ops/messageprop.jl")

# layers
include("layers/gcn.jl")
include("layers/graphode.jl")

# ogbhiv
include("features/molecules/ogbhiv.jl")

end