text:
@with_kw struct DGP
    suppZ::Array{<:Real, 2}
    densZ::Array{<:Real, 1}
    pscore::Array{<:Real, 1}
    mtrs::Tuple{MTR, MTR}

    function DGP(suppZ, densZ, pscore, mtrs)
        @assert size(suppZ, 1) == length(densZ)
        @assert length(densZ) == length(pscore)

        new(suppZ, densZ, pscore, mtrs)
    end
end
export DGP

function find_pscore(z, dgp::DGP)
    idx = findall([z == dgp.suppZ[i, :] for i in 1:size(dgp.suppZ, 1)])
    @assert length(idx) == 1 # find one and only one match
    return dgp.pscore[idx][1]
end
export find_pscore

function find_density(z, dgp::DGP)
    idx = findall([z == dgp.suppZ[i, :] for i in 1:size(dgp.suppZ, 1)])
    @assert length(idx) == 1 # find one and only one match
    return dgp.densZ[idx][1]
end
export find_density

function default_mtrs()
    # Note: Coefficients here are in the Bernstein basis
    # In the paper I wrote it using the standard basis
    mtrbasis = bernstein_basis(1)
    mtr0 = MTR(mtrbasis, hcat(.5, .4))
    mtr1 = MTR(mtrbasis, hcat(.8, .4))
    return (mtr0, mtr1)
end
export default_mtrs

function default_suppZ()
    return [0 0 ;
            0 1 ;
            1 0 ;
            1 1  ]
end

function default_pscore()
    return [.3, .5, .6, .7]
end

function illustration_dgp()
    dgp = DGP(
        suppZ = default_suppZ(),
        densZ = fill(.25, 4),
        pscore = default_pscore(),
        mtrs = default_mtrs()
    )
end
export illustration_dgp

function simulation_dgp()
    dgp = DGP(
        suppZ = default_suppZ(),
        densZ = [.4, .3, .1, .2],
        pscore = default_pscore(),
        mtrs = default_mtrs()
    )
end
export simulation_dgp

function prte_dgp(; ℓ_gen = 1)
                    # E[Y | G = g]    Pr[G = g]
    n = NamedArray([    .6                .1 ;   # at
                        .6                .1 ;   # ec
                        .1                .4 ;   # 1c
                        .6                .2 ;   # 2c
                        .2                .1 ;   # rc
                        .1                .1  ], # nt
                   (OrderedDict("at" => 1, "ec" => 2, "1c" => 3,
                                "2c" => 4, "rc" => 5, "nt" => 6),
                    OrderedDict("ExpY" => 1, "Pr" => 2)),
                   ("Rows", "Cols"))
    @assert sum(n[:, "Pr"]) == 1
    # default, but writing it out here so construction below is clear
    suppZ = [0 0 ;
             0 1 ;
             1 0 ;
             1 1  ]
    pscore = [n["at", "Pr"],
              n["at", "Pr"] + n["ec", "Pr"] + n["2c", "Pr"],
              n["at", "Pr"] + n["ec", "Pr"] + n["1c", "Pr"],
              1 - n["nt", "Pr"]]
    knots = vcat(0, pscore, 1)
    @assert length(knots) == 6 # otherwise construction below will fail
    @assert pscore == sort(pscore) # also needed for construction below

    θ = [fill(NaN, (2, length(knots) - 1)) for ℓ in 1:2]
    ############################################################################
    # Mapping from choice groups and their average effects to MTEs averaged
    # over one of the five knot regions.
    #
    # e.g. ℓ = 1, given Z2 = 0
    # if you are in (0, p(0,0)]
    # then you must be an always taker
    # if you are in (p(0,0), p(1,0)]
    # then you could be an eager complier or Z1-complier
    # if you are in (p(1,0), p(1,1)]
    # you could be a reluctant complier, Z2-complier, or never-taker
    #
    # e.g. ℓ = 1, given Z2 = 1
    # if you are in (0, p(0,1)]
    # then you could be an always-taker, eager complier, or Z2-complier
    # if you are in (p(0,1), p(1,1)]
    # then you could be a reluctant complier or a Z1-complier
    # if you are in (p(1,1), 1]
    # then you must be a never-taker
    ############################################################################
    # ℓ = 1
    #
    # Z2 = 0
    θ[1][1, 1]    = n["at", "ExpY"]
    θ[1][1, 2:3] .= group_average(n, ["ec", "1c"])
    θ[1][1, 4:5] .= group_average(n, ["rc", "2c", "nt"])
    # Z2 = 1
    θ[1][2, 1:2] .= group_average(n, ["at", "2c", "ec"])
    θ[1][2, 3:4] .= group_average(n, ["1c", "rc"])
    θ[1][2, 5]    = n["nt", "ExpY"]
    # ℓ = 2
    #
    # Z1 = 0
    θ[2][1, 1]    = n["at", "ExpY"]
    θ[2][1, 2]    = group_average(n, ["ec", "2c"])
    θ[2][1, 3:5] .= group_average(n, ["rc", "1c", "nt"])
    # Z1 = 1
    θ[2][2, 1:3] .= group_average(n, ["at", "1c", "ec"])
    θ[2][2, 4]    = group_average(n, ["2c", "rc"])
    θ[2][2, 5]    = n["nt", "ExpY"]

    # Assume Y(0) = 0 to make life easy
    notℓ = 2 - (ℓ_gen - 1) # = 2 if ℓ_gen is 1, and 1 if ℓ_gen is 2
    mtr0 = MTR(switching_constantspline_basis(knots, notℓ = notℓ),
               zeros(size(θ[1])))
    mtr1 = MTR(switching_constantspline_basis(knots, notℓ = notℓ), θ[ℓ_gen])

    dgp = DGP(suppZ = default_suppZ(),
              pscore = pscore,
              densZ = [.25, 0, .5, .25], # note misnomer in "support" -- oh well
              mtrs = (mtr0, mtr1))
    return dgp
end
export prte_dgp

function group_average(n, glist)
    dot(n[glist, "ExpY"], n[glist, "Pr"])/sum(n[glist, "Pr"])
end


text:
@registerN  Charlier AbstractCDOP1 μ
export  Charlier

"""
    Charlier{μ}

References: [Koekoek and Swarttouw §1.12](https://arxiv.org/pdf/math/9602214.pdf)
"""
Charlier

abcde(::Type{<:Charlier{μ}})  where {μ} = NamedTuple{(:a,:b,:c,:d,:e)}((0,1,0,-1,μ))

basis_symbol(::Type{<:Charlier{μ}}) where {μ} = "cᵘ"
Polynomials.domain(::Type{<:Charlier{μ}}) where {μ} = Polynomials.Interval(0, Inf)
weight_function(::Type{<:Charlier{μ}}) where {μ} = x  ->  μ^x/gamma(1+x)

function k1k0(P::Type{<:Charlier{μ}}, n::Int)  where {μ}
    -one(eltype(P))/μ
end

function classical_hypergeometric(P::Type{<:Charlier{μ}}, n::Int, x) where {μ}
    pFq((-n,-x),(),-1/μ)
end

ext: jl


text:
using Test
using POMDPPolicies
using POMDPs
using FiniteHorizonPOMDPs
using BeliefUpdaters
using POMDPSimulators
using POMDPModelTools
using POMDPModels
using Random

@testset "alpha" begin
    include("test_alpha_policy.jl")
end
@testset "function" begin
    include("test_function_policy.jl")
end
@testset "stochastic" begin
    include("test_stochastic_policy.jl")
end
@testset "utility" begin
    include("test_utility_wrapper.jl")
end
@testset "vector" begin
    include("test_vector_policy.jl")
end
@testset "random" begin
    include("test_random_solver.jl")
end
@testset "pretty_printing" begin
    include("test_pretty_printing.jl")
end
@testset "exploration policies" begin
    include("test_exploration_policies.jl")
end
@testset "playback policies" begin
    include("test_playback_policy.jl")
end
@testset "test_staged_policy" begin
    include("test_staged_policy.jl")
end


text:
using ARMA
using Test

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

max_stars_repo_name: wmkouw/ARMA

max_issues_repo_name: wmkouw/ARMA

max_forks_repo_name: wmkouw/ARMA


text:
"""
    TeXElement

Abstract super type for elements of a MathTeX string.

Currently represent either
    - a char
    - a space
    - a horizontal or vertical line
    - a group of TeXElements positioned relative to each other
"""
abstract type TeXElement end

"""
    leftinkbound(elem::TeXElement)

Return the position of the leftmost position the ink of the element uses.
"""
function leftinkbound end

"""
    rightinkbound(elem::TeXElement)

Return the position of the rightmost position the ink of the element uses.
"""
function rightinkbound end

"""
    bottominkbound(elem::TeXElement)

Return the position of the lowest position the ink of the element uses.
"""
function bottominkbound end

"""
    topinkbound(elem::TeXElement)

Return the position of the highest position the ink of the element uses.
"""
function topinkbound end

"""
    advance(elem::TeXElement)

Return the horizontal distance between the origin of the element and the origin of the
next element to be drawn, without kerning.
"""
advance(x::TeXElement) = inkwidth(x)

"""
    ascender(elem::TeXElement)

Return the vertical spaced reserved above the element, measured from the
baseline.

Only makes sense for TeXChar and derived elements, is set to zero for all other
for convenience.
"""
ascender(x::TeXElement) = 0

"""
    descender(elem::TeXElement)

Return the vertical spaced reserved below the element, measured from the
baseline.

Only makes sense for TeXChar and derived elements, is set to zero for all other
for convenience.
"""
descender(x::TeXElement) = 0

"""
    xheight(elem::TeXElement)

Return the height of the 'x' character associated to the font of the element.

Only makes sense for TeXChar and derived elements, is set to zero for all other
for convenience.
"""
xheight(x::TeXElement) = 0

"""
    hmid(elem::TeXElement)

Return the horizontal middle of the element ink.
"""
hmid(x::TeXElement) = 0.5*(leftinkbound(x) + rightinkbound(x))

"""
    vmid(elem::TeXElement)

Return the vertical middle of the element ink.
"""
vmid(x::TeXElement) = 0.5*(bottominkbound(x) + topinkbound(x))

"""
    inkwidth(elem::TeXElement)

Return the total width of the element ink.
"""
inkwidth(x::TeXElement) = rightinkbound(x) - leftinkbound(x)

"""
    inkheight(elem::TeXElement)

Return the total height of the element ink.
"""
inkheight(x::TeXElement) = topinkbound(x) - bottominkbound(x)

"""
    TeXChar(char, font)

A MathTeX character with an associated font.
"""
struct TeXChar <: TeXElement
    char::Char
    font::FTFont
end

function TeXChar(char, fontset::FontSet, char_type, command=nothing)
    return TeXChar(char, get_font(fontset, char_type))
end

for inkfunc in (:leftinkbound, :rightinkbound, :bottominkbound, :topinkbound)
    @eval $inkfunc(char::TeXChar) = $inkfunc(get_extent(char.font, char.char))
end

advance(char::TeXChar) = begin
    hadvance(get_extent(char.font, char.char))
end
ascender(char::TeXChar) = ascender(char.font)
descender(char::TeXChar) = descender(char.font)
xheight(char::TeXChar) = xheight(char.font)

Base.show(io::IO, tc::TeXChar) =
    print(io, "TeXChar '$(tc.char)' [U+$(uppercase(string(codepoint(tc.char), base=16, pad=4))) in $(tc.font.family_name) - $(tc.font.style_name)]")

"""
    Space

A MathTeX space of a given width.
"""
struct Space{T} <: TeXElement
    width::T
end

leftinkbound(::Space) = 0
rightinkbound(s::Space) = s.width
bottominkbound(::Space) = 0
topinkbound(::Space) = 0

"""
    ScaledChar

A scaled TeXChar.
"""
struct ScaledChar{T} <: TeXElement
    char::TeXChar
    scale::T
end

for func in (:leftinkbound, :rightinkbound, :bottominkbound, :topinkbound,
             :advance, :ascender, :descender, :xheight)
    @eval $func(scaled::ScaledChar) = $func(scaled.char) * scaled.scale
end

"""
    Vline

A vertical line.
"""
struct VLine{T} <: TeXElement
    height::T
    thickness::T
end

VLine(height, thickness) = VLine(promote(height, thickness)...)

leftinkbound(line::VLine) = -line.thickness/2
rightinkbound(line::VLine) = line.thickness/2
bottominkbound(line::VLine{T}) where T = min(line.height, zero(T))
topinkbound(line::VLine{T}) where T = max(line.height, zero(T))

"""
    Hline

A horizontal line.
"""
struct HLine{T} <: TeXElement
    width::T
    thickness::T
end

HLine(height, thickness) = HLine(promote(height, thickness)...)

leftinkbound(line::HLine{T}) where T = min(line.width, zero(T))
rightinkbound(line::HLine{T}) where T = max(line.width, zero(T))
bottominkbound(line::HLine) = -line.thickness/2
topinkbound(line::HLine) = line.thickness/2

"""
    Group

A group of TeXElements.

Positions and scales are relative to the group.
"""
struct Group{T} <: TeXElement
    elements::Vector{<:TeXElement}
    positions::Vector{Point2f0}
    scales::Vector{T}
end

Group(elements, positions) = Group(elements, positions, ones(length(elements)))

xpositions(g::Group) = [p[1] for p in g.positions]
ypositions(g::Group) = [p[2] for p in g.positions]

function leftinkbound(g::Group)
    lefts = leftinkbound.(g.elements) .* g.scales .+ xpositions(g)
    return minimum(lefts)
end

function rightinkbound(g::Group)
    rights = rightinkbound.(g.elements) .* g.scales .+ xpositions(g)
    return maximum(rights)
end

function bottominkbound(g::Group)
    bottoms = bottominkbound.(g.elements) .* g.scales .+ ypositions(g)
    return minimum(bottoms)
end

function topinkbound(g::Group)
    tops = topinkbound.(g.elements) .* g.scales .+ ypositions(g)
    return maximum(tops)
end

function advance(g::Group)
    adv = xpositions(g) .+ advance.(g.elements) .* g.scales
    return maximum(adv)
end

function ascender(g::Group)
    asc = ypositions(g) .+ ascender.(g.elements) .* g.scales
    return maximum(asc)
end

function descender(g::Group)
    des = ypositions(g) .+ descender.(g.elements) .* g.scales
    return minimum(des)
end

xheight(g::Group) = maximum(xheight.(g.elements) .* g.scales)


text:
printstyled("Making plots ... "; color = :green);println(" ")
##################################################################
using PyPlot

figure(1)
plot(J[1,1,:],Energy)
legend(["a","ab","abc"]);
xlabel("Exchange interaction J (meV)");
ylabel("Energy (meV)")

figure(2)
plot(J[1,1,:],Energy.-Energy[:,1])
xlabel("Exchange interaction J (meV)");
ylabel("Excitation Energy (meV)")
##################################################################
println("Success making plots!");


text:
###

println("Low level methods")

decode_mode(d, 0b00000001)
println("decode_mode for digit 1 only\n")
sleep(5)

decode_mode(d)
println("decode_mode for all digits\n")
sleep(5)

decode_mode(d, 0b0)
println("decode_mode for none\n")
sleep(5)

println("change intensity")
for i in 1:16
    set_intensity(d, i)
    println("$i/16")
    sleep(0.5)
end
for i in collect(16:(-1):1)
    set_intensity(d, i)
    println("$i/16")
    sleep(0.5)
end
sleep(1)

println("change limit")
for i in collect(8:(-1):1)
    set_limit(d, i)
    println("$i/8")
    sleep(0.5)
end

for i in 1:8
    set_limit(d, i)
    println("$i/8")
    sleep(0.5)
end
sleep(1)

println("shutdown_mode_on()")
shutdown_mode_on(d)
sleep(5)

println("shutdown_mode_off()")
shutdown_mode_off(d)
sleep(5)

println("test_mode_on()")
test_mode_on(d)
sleep(5)

println("test_mode_off()")
test_mode_off(d)
sleep(5)

print("switch to decode and set 3,3,3,3 at positions 8,7,6,5")
decode_mode(d)
write_digit(d, UInt8(3), 8)
sleep(1)
write_digit(d, UInt8(3), 7)
sleep(1)
write_digit(d, UInt8(3), 6)
sleep(1)
write_digit(d, UInt8(3), 5)
sleep(1)


text:
module NativeHTML

import Base.div

struct HTML
    data::Vector{UInt8}
end

Base.showable(::MIME"text/html", _::HTML) = true

Base.showable(::MIME{Symbol("juliavscode/html")}, _::HTML) = true

function Base.show(io::IO, html::HTML)
    write(io, String(copy(html.data)))
end

function Base.show(io::IO, ::MIME"text/plain", html::HTML)
    (isdefined(Main, :IJulia) && Main.IJulia.inited ||
     isdefined(Main, :Juno) && Main.Juno.isactive()) && return
    filename = "NativeHTML.html"
    open(filename, "w") do io
        write(io, html.data)
    end
    if Sys.isapple()
        run(`open $(filename)`)
    elseif Sys.iswindows()
        cmd = get(ENV, "COMSPEC", "cmd")
        run(`$(ENV["COMSPEC"]) /c start $(filename)`)
    elseif Sys.isunix()
        run(`xdg-open $(filename)`)
    end
    return filename
end

function Base.show(io::IO, ::MIME"text/html", html::HTML)
    write(io, html.data)
end

function Base.show(io::IO, ::MIME{Symbol("juliavscode/html")}, html::HTML)
    write(io, html.data)
end

function Base.write(io::IO, html::HTML)
    write(io, html.data)
end

function html(f::Function)
    io = IOBuffer()
    println(io, "<!DOCTYPE html>")
    println(io, "<html>")
    task_local_storage(:htmlio, io)
    task_local_storage(:level, 1)
    f()
    println(io, "</html>")
    HTML(take!(io))
end

export html

function fragment(f::Function; level::Int = 0)
    io = IOBuffer()
    task_local_storage(:htmlio, io)
    task_local_storage(:level, level)
    f()
    HTML(take!(io))
end

export fragment

const PRIMITIVES = Dict(
    # metadata
    :head => true,
    :title => true,
    :base => false,
    :link => false,
    :meta => false,
    :style => true,
    # sections
    :body => true,
    :article => true,
    :section => true,
    :nav => true,
    :aside => true,
    :h1 => true,
    :h2 => true,
    :h3 => true,
    :h4 => true,
    :h5 => true,
    :h6 => true,
    :hgroup => true,
    :header => true,
    :footer => true,
    :address => true,
    # grouping
    :p => true,
    :hr => false,
    :pre => true,
    :blockquote => true,
    :ol => true,
    :ul => true,
    :menu => true,
    :li => true,
    :dl => true,
    :dt => true,
    :dd => true,
    :figure => true,
    :figcaption => true,
    :main => true,
    :div => true,
    # text-level
    :a => true,
    :em => true,
    :strong => true,
    :small => true,
    :s => true,
    :cite => true,
    :q => true,
    :dfn => true,
    :abbr => true,
    :ruby => true,
    :rt => true,
    :rp => true,
    :data => true,
    :time => true,
    :code => true,
    :var => true,
    :samp => true,
    :kbd => true,
    :sub => true,
    :sup => true,
    :i => true,
    :b => true,
    :u => true,
    :mark => true,
    :bdi => true,
    :bdo => true,
    :span => true,
    :br => false,
    :wbr => false,
    # edits
    :ins => true,
    :del => true,
    # embedded
    :picture => true,
    :source => false,
    :img => false,
    :iframe => true,
    :embed => false,
    :object => true,
    :param => false,
    :video => true,
    :audio => true,
    :track => false,
    :map => true,
    :area => false,
    # tables
    :table => true,
    :caption => true,
    :colgroup => true,
    :col => false,
    :tbody => true,
    :thead => true,
    :tfoot => true,
    :tr => true,
    :td => true,
    :th => true,
    # forms
    :form => true,
    :label => true,
    :input => false,
    :button => true,
    :select => true,
    :datalist => true,
    :optgroup => true,
    :option => true,
    :textarea => true,
    :output => true,
    :progress => true,
    :meter => true,
    :fieldset => true,
    :legend => true,
    # interactive
    :details => true,
    :summary => true,
    # scripting
    :script => true,
    :noscript => true,
    :template => true,
    :slot => true,
    :canvas => true,
)

for (primitive, notvoid) in PRIMITIVES
    eval(quote
        function $primitive(txt::String = ""; kwargs...)
            let io = task_local_storage(:htmlio), level = task_local_storage(:level)
                print(io, " " ^ level, "<", $primitive)
                for (arg, val) in kwargs
                    print(io, " ", replacenotallowed(arg))
                    val === nothing || print(io, "=\"", val, "\"")
                end
                if txt === ""
                    $notvoid ? println(io, "/>") : println(io, ">")
                else
                    println(io, ">", txt, "</", $primitive, ">")
                end
            end
        end
    end)
end

for primitive in keys(PRIMITIVES)
    eval(quote
        export $primitive
    end)
end

for primitive in keys(filter(d->last(d), PRIMITIVES))
    eval(quote
        function $primitive(f::Function; kwargs...)
            let io = task_local_storage(:htmlio), level = task_local_storage(:level)
                print(io, " " ^ level, "<", $primitive)
                for (arg, val) in kwargs
                    print(io, " ", replacenotallowed(arg))
                    val === nothing || print(io, "=\"", val, "\"")
                end
                println(io, ">")
                task_local_storage(:level, level + 1)
                f()
                task_local_storage(:level, level)
                println(io, " " ^ level, "</", $primitive, ">")
            end
        end
    end)
end

function text(txt::String)
    let io = task_local_storage(:htmlio), level = task_local_storage(:level)
        println(io, " " ^ level, txt)
    end
end

export text

function cdata(txt::String)
    let io = task_local_storage(:htmlio), level = task_local_storage(:level)
        println(io, " " ^ level, "<![CDATA[", txt, "]]>")
    end
end

export cdata

function replacenotallowed(key::Symbol)
    str = replace(String(key), '_' => '-')
    str[begin] === '-' && return str[begin+1:end]
    str
end

end


text:
"""
    using Plots
    plot(pi::ParamInterval)

Plots profile `L(theta)` for parameter `theta_num`,
`identifiability level`, `identifiability interval`.
Use `update_profile_points!(pi::ProfileInterval)` function to refine
profile points and make your plot more smooth
"""
@recipe function f(pi::ParamInterval)
    # get points for plot
    xs, ys = get_grid(pi)
    xguide --> "theta"
    yguide --> "L(theta)"

    if length(xs) > 1
        @series begin
            label --> "CI endpoint"
            seriestype --> :line
            linestyle --> :dot
            markershape --> :circle
            markercolor --> :green
            line := (2.2,:green)
            (xs, ys)
        end
    end

    # critical level subplot
    @series begin
        label --> "Identifiability level"
        seriestype --> :hline
        line := (2.5, :purple)
        [pi.input.options[:loss_crit]]
    end

    if !isa(pi.result[1].value, Nothing) || !isa(pi.result[2].value, Nothing)
        @series begin
            label --> "Identifiability inteval"
            seriestype --> :vline
            line := (2.1, :pink)
            [
                !isa(pi.result[1].value, Nothing) ? pi.result[1].value : NaN,
                !isa(pi.result[2].value, Nothing) ? pi.result[2].value : NaN
            ]

            #[pi.result[1].value,pi.result[2].value]
        end
    end

    @series begin
        label --> "Initial point"
        seriestype --> :scatter
        markershape --> :diamond
        markercolor --> :orange
        [(pi.input.theta_init[pi.input.theta_num], pi.loss_init)]
    end

end

function get_grid(pi::ParamInterval)
    xs = Vector{Float64}()
    ys = Vector{Float64}()
    #ns = Vector{String}()

    ep_start = pi.input.theta_init[pi.input.theta_num]
    push!(xs, ep_start)
    push!(ys, pi.loss_init)
    #push!(ns, "1")

    for ep in pi.result
        x_pps, y_pps = get_pps(ep)
        append!(xs, x_pps)
        append!(ys, y_pps)
        #append!(ns, n_pps)
    end
    return (xs, ys)
end

function get_pps(ep::EndPoint)
    l = length(ep.profilePoints)

    xs = Vector{Float64}(undef, l)
    ys = Vector{Float64}(undef, l)
    #ns = Vector{String}(undef, l)

    for i in 1:l
        xs[i] = ep.profilePoints[i].value
        ys[i] = ep.profilePoints[i].loss
        #ns[i] = ep.profilePoints[i].counter != nothing ?
        #    string(ep.profilePoints[i].counter) : ""
    end
    (xs,ys)
end

"""
    update_profile_points!(
        pi::ParamInterval;
        max_recursions::Int = 2
        )

Refines profile points to make your plot more smooth. Internally uses
`adapted_grid` to compute additional profile points.
See `PlotUtils.adapted_grid`.

## Arguments
- `pi` : `ParamInterval` structure to update.
- `max_recursions` : how many times each interval is allowed to be refined (default: 2).
"""
function update_profile_points!(
    pi::ParamInterval;
    max_recursions::Int = 2
)
    for ep in pi.result
        if ep.status != :SCAN_BOUND_REACHED
            ep_start = pi.input.theta_init[pi.input.theta_num]
            update_profile_endpoint!(
                ep.profilePoints,
                ep.direction == :left ?
                    (ep.value, ep_start) :
                    (ep_start, ep.value),
                pi.input.theta_init,
                pi.input.theta_num,
                pi.input.loss_func,
                pi.input.options[:local_alg],
                pi.input.options[:theta_bounds],
                pi.input.options[:loss_tol];
                max_recursions
            )
        else
            @info "$(string(ep.direction)) - a half-open interval"
        end
    end
    return nothing
end

function update_profile_endpoint!(
    pps_arr::Vector{ProfilePoint},
    interval::Tuple{Float64,Float64},
    init_params::Vector{Float64},
    id::Int,
    loss_func::Function,
    fit_alg::Symbol,
    bounds::Vector{Tuple{Float64,Float64}},
    tol::Float64;
    max_recursions::Int = 2
)
    params = copy(init_params)

    # bounds
    lb = minimum.(bounds)
    ub = maximum.(bounds)

    # function to be used in adapted_grid
    function profile_func(x::Float64)
        if length(params) == 1
            return loss_func([x])
        else
            # ! params_intervals_one_side can be used here
            fit_params_func = (p,g) -> loss_func(p)
            opt = Opt(fit_alg, length(params))

            params[id] = x
            min_objective!(opt, fit_params_func)

            # excluding params[id] from optimization
            lb[id] = x
            ub[id] = x
            lower_bounds!(opt,lb)
            upper_bounds!(opt,ub)
            ftol_abs!(opt,tol)

            (loss,minx,ret) = optimize(opt,params)

            loss = loss_func(minx)
            push!(pps_arr, ProfilePoint(x, loss, minx, ret, nothing))
            return loss
        end
    end

    # adapted_grid
    adapted_grid(profile_func,interval; max_recursions = max_recursions)
    return nothing
end

# deprecated
#=
function get_adapted_grid(
    interval::Tuple{Float64,Float64},
    init_params::Vector{Float64},
    id::Int64,
    loss_func::Function,
    fit_alg::Symbol,
    bounds::Vector{Tuple{Float64,Float64}},
    tol::Float64;
    max_recursions::Int = 1
)
    params = copy(init_params)
    # bounds
    lb = minimum.(bounds)
    ub = maximum.(bounds)

    # function to be used in adapted_grid
    function profile_func(x::Float64)
        if length(params) == 1
            return loss_func([x])
        else
            # ! params_intervals_one_side can be used here
            fit_params_func = (p,g) -> loss_func(p)
            opt = Opt(fit_alg, length(params))

            params[id] = x
            min_objective!(opt, fit_params_func)

            # excluding params[id] from optimization
            lb[id] = x
            ub[id] = x
            lower_bounds!(opt,lb)
            upper_bounds!(opt,ub)
            ftol_abs!(opt,tol)

            (loss,minx,ret) = optimize(opt,params)

            return loss
        end
    end

    # adapted_grid
    adapted_grid2(profile_func,interval; max_recursions = max_recursions)
end
=#


text:
using Documenter, TextParse

makedocs(
    modules = [TextParse],
    sitename = "TextParse.jl",
    authors = "Shashi Gowda",
    pages = Any["Home" => "index.md"],
)

deploydocs(
    repo = "github.com/JuliaComputing/TextParse.jl.git"
)


text:
module HealthMLBase

end # module


text:
# constructs a spiking layer out of a parametric type, like HHNeuron
function SpikingLayer{T <: AbstractSpikingNeuron}(::Type{T}, name::ASCIIString, N::Int64)
	units = [T() for i=1:N]
	SpikingLayer(name,N,units)
end

function update!(layer::SpikingLayer{SimpleSpikingNeuron}, pres::Vector{SpikingConnection}, doSpike::BitArray{1}, t::Float64)
	for i=1:layer.N
		update!(layer.units[i], doSpike[i], t)
	end
end

function update!(layer::SpikingLayer{HHNeuron}, pres::Vector{SpikingConnection}, Iinj::Vector{Float64}, t::Float64)
	Ge_syn = zeros(layer.N)
	Gi_syn = zeros(layer.N)

#	@bp layer.name == "LGN_L"

	# integrate all synaptic conductances
	for sc in pres
		for syn in sc.synapses
			if syn.stype == EXC
				Ge_syn[syn.t] += syn.w * syn.g
			else
				Gi_syn[syn.t] += syn.w * syn.g
			end
		end
	end

	for i=1:layer.N
		update!(layer.units[i], Ge_syn[i], Gi_syn[i], Iinj[i], t)
	end
end

function update!(layer::SpikingLayer{HHNeuron}, pres::Vector{SpikingConnection}, Iinj::Nothing, t::Float64)
	update!(layer,pres,zeros(layer.N),t)
end

function didSpike(layer::SpikingLayer)
	[unit.didSpike for unit in layer.units]
end

function reset!(layer::SpikingLayer)
	for unit in layer.units
		reset!(unit)
	end
end


text:
##  Copyright 2022 Myrl Marmarelis
##
##  Licensed under the Apache License, Version 2.0 (the "License");
##  you may not use this file except in compliance with the License.
##  You may obtain a copy of the License at
##
##    http://www.apache.org/licenses/LICENSE-2.0
##
##  Unless required by applicable law or agreed to in writing, software
##  distributed under the License is distributed on an "AS IS" BASIS,
##  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##  See the License for the specific language governing permissions and
##  limitations under the License.

no_dynamics(T) = (_...) -> (T(-Inf), T(Inf))

lipschitz_condition(slope) =
  ((prev_point, prev_weight), point, (succ_point, succ_weight)) -> begin
    prev_range = slope * (point - prev_point)
    succ_range = slope * (succ_point - point)
    if isnan(prev_point)
      (succ_weight - succ_range, succ_weight + succ_range)
    elseif isnan(succ_point)
      (prev_weight - prev_range, prev_weight + prev_range)
    else
      upper = min(prev_weight + prev_range, succ_weight + succ_range)
      lower = max(prev_weight - prev_range, succ_weight - succ_range)
      (lower, upper)
    end
  end

# `values` need not be monotonic on `points`
# on arbitrary closures, so that it's recompiled for each specialized call
function bound_expectation(static_bounds::Function,
    points::AbstractVector{T}, values::AbstractVector{T};
    dynamic_bounds::Function=no_dynamics(T))::NTuple{2,T} where T
  ordering = sortperm(values)
  points = points[ordering] # method to ensure this is done in place?
  values = values[ordering]
  maximized = maximize_expectation(static_bounds, points, values;
    dynamic_bounds).expectation
  reverse!(points)
  reverse!(values)
  values .*= -1
  minimized = -maximize_expectation(static_bounds, points, values;
    dynamic_bounds).expectation
  (minimized, maximized)
end

# assumes `points` are sorted along `values` here
function maximize_expectation(static_bounds::Function,
    points::AbstractVector{T}, values::AbstractVector{T};
    dynamic_bounds::Function=no_dynamics(T)) where T
  n_draws = length(points)
  weights = [ static_bounds(point)[2] for point in points ]
  point_ordering = sortperm(points)
  point_neighbors = Dict( ( point_ordering[p] =>
      (get(point_ordering, p-1, 0), get(point_ordering, p+1, 0)) )
    for p in 1:n_draws )
  n_iterations = 0
  mutation = true
  while mutation
    n_iterations += 1
    mutation = false
    for index in 1:n_draws
      (this_point, this_value) = (points[index], values[index])
      # all this noise in order to enforce `points` continuity for non-monotonic `values`
      (prev_index, succ_index) = point_neighbors[index]
      previous = (prev_index > 0 ?
        (points[prev_index], weights[prev_index]) : T.((NaN, NaN)) )
      successive = (succ_index > 0 ?
        (points[succ_index], weights[succ_index]) : T.((NaN, NaN)) )
      # we call the below only on legal neighbor assignments.
      (dynamic_lo, dynamic_up) = dynamic_bounds(previous, this_point, successive)
      (static_lo, static_up) = static_bounds(this_point)
      unscaled_derivative = sum(1:n_draws) do other
        weights[other] * (this_value - values[other])
      end
      if unscaled_derivative < 0
        new_weight = max(static_lo, dynamic_lo)
        if !isapprox(new_weight, weights[index])
          mutation = true
        end
        weights[index] = new_weight
      else
        # done in this one iteration; move along
        break
      end
    end
  end
  expectation = sum(weights) \ sum(1:n_draws) do index
    weights[index] * values[index]
  end
  (; expectation, weights, n_iterations )
end

max_stars_repo_name: marmarelis/TreatmentCurves

max_issues_repo_name: marmarelis/TreatmentCurves

max_forks_repo_name: marmarelis/TreatmentCurves


text:
"""
	add_data(d::ConjugatePostDistribution, X)

Add data to Posterior Distribution.
"""
function add_data(d::ConjugatePostDistribution, X)
	dist = deepcopy(d)
	add_data!(dist, X)
	return dist
end

"""
	add_data!(d::ConjugatePostDistribution, X)

Add data to Posterior Distribution (inplace).
"""
function add_data!(d::MultivariateConjugatePostDistribution, X::AbstractMatrix)
    @simd for i in 1:size(X, 1)
			@inbounds add_data!(d, @view X[i,:])
		end
end

function add_data!(d::UnivariateConjugatePostDistribution, X::AbstractVector)
	@simd for i in 1:size(X, 1)
		@inbounds add_data!(d, X[i])
	end
end

function add_data!(d::WishartGaussian, X::AbstractVector)
	@assert length(X) == d.D

	d.n += 1
  d.sums += X
  d.ssums += X * X'
end

function add_data!(d::GaussianDiagonal, X::AbstractVector)
	for (dim, dist) in enumerate(d.dists)
		add_data!(dist, X[dim])
	end
end

function add_data!(d::ContinuousUnivariateConjugatePostDistribution, X::AbstractFloat)
	d.n += 1
	d.sums += X
	d.ssums += X^2
end

function add_data!(d::DirichletMultinomial, X::AbstractVector)

    # process samples
    Dmin = minimum(X)
    Dmax = maximum(X)

    if Dmax > d.D
        throw(ArgumentError("Value of X and is larger than Multinomial Distribution!"))
    end

    if Dmin < 1
        throw(ArgumentError("Value of X and is smaller than 1!"))
    end

    d.n += length(X)

    for x in X
        d.counts[x] += 1
    end

    d.dirty = true
end

function add_data!(d::DirichletMultinomial, X::SparseVector{Int, Int})
		@assert length(X) == d.D

    d.n += 1
    d.counts += X
    d.dirty = true
end

function add_data!(d::BetaBernoulli, X::Integer)
	d.successes += X
	d.n += 1
end

"""
	remove_data(d::ConjugatePostDistribution, X::AbstractArray)

Remove data from Posterior Distribution.
"""
function remove_data(d::ConjugatePostDistribution, X)
		dist = deepcopy(d)
		remove_data!(dist, X)
		return dist
end

"""
	remove_data!(d::ConjugatePostDistribution, X)

Remove data from Posterior Distribution (inplace).
"""
function remove_data!(d::MultivariateConjugatePostDistribution, X::AbstractMatrix)
    @simd for i in 1:size(X, 1)
			@inbounds remove_data!(d, @view X[i,:])
		end
end

function remove_data!(d::UnivariateConjugatePostDistribution, X::AbstractVector)
	@simd for i in 1:size(X, 1)
		@inbounds remove_data!(d, X[i])
	end
end

function remove_data!(d::WishartGaussian, X::AbstractVector)
		@assert length(X) == d.D

    d.n -= 1
    d.sums -= X
    d.ssums -= X * X'
end

function remove_data!(d::ContinuousUnivariateConjugatePostDistribution, X::AbstractFloat)
    d.n -= 1
    d.sums -= X
    d.ssums -= X^2
end

function remove_data!(d::GaussianDiagonal, X::AbstractVector)
	for (dim, dist) in enumerate(d.dists)
		remove_data!(dist, X[dim])
	end
end

function remove_data!(d::DirichletMultinomial, X::AbstractVector)

    # process samples
    Dmin = minimum(X)
    Dmax = maximum(X)

    if Dmax > d.D
        throw(ArgumentError("Value of X and is larger than Multinomial Distribution!"))
    end

    if Dmin < 1
        throw(ArgumentError("Value of X and is smaller than 1!"))
    end

    d.n -= length(X)

    for x in X
        d.counts[x] -= 1
    end

    d.dirty = true
end

function remove_data!(d::DirichletMultinomial, X::SparseVector{Int, Int})

   # process samples
   @assert length(X) == d.D

   d.n -= 1
   d.counts -= X
   d.dirty = true
end

function remove_data!(d::BetaBernoulli, X::Bool)
		d.successes -= X
		d.n -= 1
end

"""
	isdistempty(d::ConjugatePostDistribution)

Check if distribution is empty (contains no samples).
"""
function isdistempty(d::ConjugatePostDistribution)
    return d.n <= 0
end

"""
	logpred(d::ConjugatePostDistribution, X)

Compute log posterior predictive.
"""
function logpred(d::WishartGaussian, x)

	if isdistempty(d)
		# posterior predictive of Normal Inverse Wishart is student-t Distribution

		C = d.Sigma0 * ((d.kappa0 + 1) / (d.kappa0 * (d.nu0 - d.D + 1)))

		dist = Distributions.MvTDist(d.nu0 - d.D + 1, d.mu0, C)
	  return Distributions.logpdf(dist, x)

	end

  # statistics
  sample_mu = d.sums / d.n

  # make sure values are not NaN
  sample_mu[sample_mu .!= sample_mu] = 0

  # compute posterior parameters
  kappa = d.kappa0 + d.n
  nu = d.nu0 + d.n
  mu = (d.kappa0 * d.mu0 + d.n * sample_mu) / kappa
  Sigma = d.Sigma0 + d.ssums - kappa * (mu * mu') + ( d.kappa0 * (d.mu0 * d.mu0') )

  # posterior predictive of Normal Inverse Wishart is student-t Distribution, see:
  # K. Murphy, Conjugate Bayesian analysis of the Gaussian distribution. Eq. 258

 C = Sigma * ((kappa + 1) / (kappa * (nu - d.D + 1)))

 if !isposdef(C)
	 throw(ErrorException("C is not positive definite! (C: $(C))"))
 end

  dist = Distributions.MvTDist(nu - d.D + 1, mu, C)
  return Distributions.logpdf(dist, x)

end

"Log PDF of Generalized student-t Distribution."
function tlogpdf(x, df, mean, sigma)

   function tdist_consts(df, sigma)
       hdf = 0.5 * df
       shdfhdim = hdf + 0.5
       v = lgamma(hdf + 1/2) - lgamma(hdf) - 0.5*log(df) - 0.5*log(pi) - log(sigma)
       return (shdfhdim, v)
   end

    shdfhdim, v = tdist_consts(df, sigma)

    xx = x .- mean
    xx = (xx ./ sqrt(sigma)).^2

    p = 1/df * xx

    return v - log((1 + p).^((df+1) / 2))
end

"Log PDF for NormalGamma."
function logpred(d::GammaNormal, x)

   if d.n == 0
      # posterior predictive of Normal Gamma is student-t Distribution
      df = 2 * d.α0
      mean = d.μ0
      sigma = ( d.β0 * (d.λ0 + 1) ) / (d.λ0 * d.α0)

      return tlogpdf(x, df, mean, sigma)
   end

    # statistics
    sample_mu = d.sums / d.n

    # make sure values are not NaN
    sample_mu = sample_mu != sample_mu ? 0 : sample_mu

    # compute posterior parameters
    μ = (d.λ0 * d.μ0 + d.sums) / (d.λ0 + d.n)
    λ = d.λ0 + d.n
    α = d.α0 + (d.n / 2)
    s = (d.ssums / d.n) - (sample_mu * sample_mu)
    β = d.β0 + 1/2 * (d.n * s + (d.λ0 * d.n * (sample_mu - d.μ0)^2 ) / (d.λ0 + d.n) )

    # posterior predictive of Normal Gamma is student-t Distribution
    df = 2 * α
    mean = μ
    sigma = ( β * (λ + 1) ) / (λ * α)

		return tlogpdf(x, df, mean, sigma)
end

function logpred(d::GaussianDiagonal, x)
	return vec(sum([logpred(di, x[dim, :]) for (dim, di) in enumerate(d.dists)]))
end

"Log PDF for NormalNormal."
function logpred(d::NormalNormal, x)
	return Float64[logpred(d, xi) for xi in x]
end

function logpred(d::NormalNormal, x::Number)

   if d.n == 0
		 return normlogpdf(d.μ0, d.σ0, x)
   end

  # statistics
  sample_mu = d.sums / d.n
	sample_var = (d.n * d.ssums - d.sums^2) / (d.n*(d.n-1))
	sample_var = (d.ssums - (d.sums^2)/d.n) / (d.n - 1)
	sample_var = 10.0

	# make sure values are not NaN
	sample_mu = sample_mu != sample_mu ? 0 : sample_mu

	σ = (sample_var * d.σ0) / (d.n * d.σ0 + sample_var)
  μ = σ * (d.n*sample_mu/sample_var + d.μ0/d.σ0)

	return normlogpdf(μ, sqrt(σ), x)
end

"Log PMF for MultinomialDirichlet."
function logpred(d::DirichletMultinomial, x)

    N = length(x)

	 # TODO: This is bad and slow code, improve!!!
    # construct sparse vector
    xx = spzeros(d.D, N)
    for i in 1:N
      xx[x[i]] += 1
    end

    m = sum(xx)
    mi = nnz(xx)

	 if d.dirty
		 d.Z2 = lgamma(d.alpha0 + d.n)
		 d.Z3 = lgamma( d.alpha0/d.D + d.counts )

		 d.dirty = false
	 end

	l1 = lgamma(m + 1) - sum(lgamma(mi + 1))
	l2 = d.Z2 - lgamma(d.alpha0 + d.n + m)
	l3 = sum( lgamma( d.alpha0 / d.D + d.counts + xx ) - d.Z3 )

    return [l1 + l2 + l3]
end

"Log PMF for MultinomialDirichlet."
function logpred(d::DirichletMultinomial, x::SparseVector{Int, Int})

   D = size(x, 1)
   N = size(x, 2)

   if N > 1
      throw(ErrorException("Multiple samples are not supported yet!"))
   end

   m = sum(x)
   mi = nnz(x)

   if d.dirty
      d.Z2 = lgamma(d.alpha0 + d.n)
      d.Z3 = lgamma( d.alpha0/d.D + d.counts )

      d.dirty = false
   end

	l1 = lgamma(m + 1) - sum(lgamma(mi + 1))
	l2 = d.Z2 - lgamma(d.alpha0 + d.n + m)
	l3 = sum( lgamma( d.alpha0 / d.D + d.counts + x ) - d.Z3 )

   return [l1 + l2 + l3]
end

"Log PMF for BernoulliBeta."
function logpred(d::BetaBernoulli, X::AbstractArray)
	return Float64[logpred(d, x) for x in X]
end

function logpred(d::BetaBernoulli, X::Bool)

  	# posterior
		α = d.α + d.successes
		β = d.β + d.n - d.successes

   return log(α) - log(α + β)
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using Statistics
using ArraysOfArrays, ElasticArrays, StatsBase
using StableRNGs

@testset "effective_sample_size" begin
    stblrng() = StableRNG(789990641)

    # Primitive MCMC of flat distribution between [0, 0, 0] and [10, 5, 8]:
    rng = stblrng()
    v = nestedview(ElasticArray{Float64, 2}(undef, 3, 0))
    push!(v, [0, 0, 0])
    for i in 1:(10^4 - 1)
        push!(v, clamp.(last(v) .+ randn(rng, 3), [0, 0, 0], [10, 5, 8]))
    end

    v1 = flatview(v)[1, :]

    @testset "BAT.bat_integrated_autocorr_len" begin
        @test @inferred(bat_integrated_autocorr_len(v1, GeyerAutocorLen())).result ≈ 52.2404651916953
        @test @inferred(bat_integrated_autocorr_len(v, GeyerAutocorLen())).result ≈ [52.240465191695314, 17.04353447359818, 38.393838710754345]

        @test @inferred(bat_integrated_autocorr_len(v1, SokalAutocorLen())).result ≈ 44.243392655975356
        @test @inferred(bat_integrated_autocorr_len(v, SokalAutocorLen())).result ≈ [44.243392655975356, 16.794891919657566, 31.94870020972804]
    end
end


text:
function LinearAlgebra.lu!(a::IndicesArray, args...; kwargs...)
    inner_lu = lu!(parent(a), args...; kwargs...)
    return LU(
        IndicesArray(getfield(inner_lu, :factors), axes(a)),
        getfield(inner_lu, :ipiv),
        getfield(inner_lu, :info)
       )
end

function Base.parent(fact::LU{T,<:IndicesArray}) where {T}
    return LU(
        parent(getfield(fact, :factors)),
        getfield(fact, :ipiv),
        getfield(fact, :info)
       )
end

function Base.getproperty(fact::LU{T,<:IndicesArray{L}}, d::Symbol) where {T, L}
    return IndicesArray(
        getproperty(parent(fact), d),
        lu_indices(a, d),
        AllUnique,
        LengthChecked
       )
end

lu_indices(a::AbstractArray, d::Symbol) = lu_indices(axes(a), d)
function lu_indices(a::Tuple, d::Symbol)
    if d === :L
        return (unsafe_axes(a, 1), nothing)
    elseif d === :U
        return (nothing, unsafe_axes(a, 2))
    elseif d === :P
        return (unsafe_axes(a, 1), unsafe_axes(a, 1))
    elseif d === :p
        return (unsafe_axes(a, 1),)
    end
end


text:
using BlockDiagonals
using BlockDiagonals: svd_blockwise, eigen_blockwise
using LinearAlgebra
using Random
using Test

@testset "linalg.jl" begin
    rng = MersenneTwister(123456)
    N1, N2, N3 = 3, 4, 5
    N = N1 + N2 + N3
    b1 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N3, N3)])
    b2 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N3, N3), rand(rng, N2, N2)])
    b3 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N2, N2)])
    A = rand(rng, N, N + N1)
    B = rand(rng, N + N1, N + N2)
    A′, B′ = A', B'
    a = rand(rng, N)
    b = rand(rng, N + N1)

    @testset "mul!" begin
        c = similar(b1)
        d = similar(Matrix(b1))
        mul!(c, b1, b1)
        mul!(d, Matrix(b1), Matrix(b1))
        @test c ≈ d
        if VERSION ≥ v"1.3"
            mul!(c, b1, b1, 2.0, 3.0)
            mul!(d, Matrix(b1), Matrix(b1), 2.0, 3.0)
            @test c ≈ d
        end
    end

    @testset "Unary Linear Algebra" begin
        @testset "$f" for f in (adjoint, det, diag, eigvals, inv, pinv, svdvals, transpose, tr)
            @test f(b1) ≈ f(Matrix(b1))
        end

        @testset "permute=$p, scale=$s" for p in (true, false), s in (true, false)
            @testset "$g" for g in (eigmin, eigmax)
                # `b2` has real eigenvals, required for `eigmin`, `b1` has Complex eigenvals
                @test g(b2; permute=p, scale=s) ≈ g(Matrix(b2); permute=p, scale=s)
            end
            @testset "eigvals" begin
                result = eigvals(b2; permute=p, scale=s)
                expected = eigvals(Matrix(b2), permute=p, scale=s)
                @static if VERSION < v"1.2"
                    # the `eigvals` method we hit above did not sort real eigenvalues pre-v1.2
                    # but some `eigvals` methods did, so we always sort real eigenvalues.
                    @test sort!(result) ≈ sort!(expected)
                else
                    @test result ≈ expected
                end
            end
        end

        @testset "Symmetric/Hermitian" begin
            @test Symmetric(b1) == Symmetric(Matrix(b1))
            @test Symmetric(b1, :L) == Symmetric(Matrix(b1), :L)
            @test Hermitian(b1) == Hermitian(Matrix(b1))
            @test Hermitian(b1, :L) == Hermitian(Matrix(b1), :L)
        end

        @testset "Eigen Decomposition" begin

            @testset "eigen $name" for (name, B) in [("", b1), ("symmetric", Symmetric(b1)), ("hermitian", Hermitian(b1))]
                E = eigen(B)
                evals_bd, evecs_bd = E
                evals, evecs = eigen(Matrix(B))

                @test E isa Eigen
                @test Matrix(E) ≈ B

                # There is no test like @test eigen(B) == eigen(Matrix(B))
                # 1. this fails in the complex case. Probably a convergence thing that leads to ever so slight differences
                # 2. pre version 1.2 this can't be expected to hold at all because the order of eigenvalues was random
                # so I sort the values/vectors (if needed) and then compare them via ≈

                @static if VERSION < v"1.2"
                    # pre-v1.2 we need to sort the eigenvalues consistently
                    # Since eigenvalues may be complex here, I use this function, which works for this test.
                    # This test is already somewhat fragile w. r. t. degenerate eigenvalues
                    # and this just makes this a little worse.
                    perm_bd = sortperm(real.(evals_bd) + 100*imag.(evals_bd))
                    evals_bd = evals_bd[perm_bd]
                    evecs_bd = evecs_bd[:, perm_bd]

                    perm = sortperm(real.(evals) + 100*imag.(evals))
                    evals = evals[perm]
                    evecs = evecs[:, perm]
                end

                @test evals_bd ≈ evals
                # comparing vectors is more difficult due to a sign ambiguity
                # So we try adding and subtracting the vectors, keeping the smallest magnitude for each entry
                # and compare that with something small.
                # I performed some tests and the largest deviation I found was ~5e-14
                @test all(min.(abs.(evecs_bd - evecs), abs.(evecs_bd + evecs)) .< 1e-13)
            end

            @testset "eigen_blockwise $name" for (name, B) in [("", b1), ("symmetric", Symmetric(b1)), ("hermitian", Hermitian(b1))]
                vals, vecs = eigen_blockwise(B)

                # check types. Eltype varies so not checked here (should be ComplexF64, Float64, Float64)
                @test vecs isa BlockDiagonal
                @test vals isa Vector

                # Note: /(Matrix, BlockDiagonal) fails. I think because we can't do factorize(BlockDiagonal)
                # this is why I convert to Matrix prior to /
                @test B ≈ vecs * Diagonal(vals) / Matrix(vecs)

                # check by block
                cumulative_size = 0
                for (i, block) in enumerate(blocks(B))
                    block_vals = vals[cumulative_size+1:cumulative_size+size(block,1)]
                    cumulative_size += size(block, 1)
                    # adapt eltype to the same to block_vals.
                    # block_vals's eltype is chosen to be compatible across all eigenvalues, thus it might be different
                    block = convert.(eltype(block_vals), block)

                    # from here on the code parallel to the test code above
                    E = Eigen(block_vals, blocks(vecs)[i])
                    evals_bd, evecs_bd = E
                    evals, evecs = eigen(block)

                    @test block ≈ Matrix(E)

                    @static if VERSION < v"1.2"
                        # sorting if needed
                        perm_bd = sortperm(real.(evals_bd) + 100*imag.(evals_bd))
                        evals_bd = evals_bd[perm_bd]
                        evecs_bd = evecs_bd[:, perm_bd]

                        perm = sortperm(real.(evals) + 100*imag.(evals))
                        evals = evals[perm]
                        evecs = evecs[:, perm]
                    end

                    @test evals_bd ≈ evals
                    @test all(min.(abs.(evecs_bd - evecs), abs.(evecs_bd + evecs)) .< 1e-13)
                end
            end
        end

        @testset "eigvals on LinearAlgebra types" begin
            # `eigvals` has different methods for different types, e.g. Hermitian
            b_herm = BlockDiagonal([Hermitian(rand(rng, 3, 3) + I) for _ in 1:3])
            @test eigvals(b_herm) ≈ eigvals(Matrix(b_herm))
            @test eigvals(b_herm, 1.0, 2.0) ≈ eigvals(Hermitian(Matrix(b_herm)), 1.0, 2.0)
        end

        # Requires a postive semi-definite matrix
        @testset "logdet" begin
            b̂1 = b1 * b1'
            @test logdet(b̂1) ≈ logdet(Matrix(b̂1))
        end
    end  # Unary
    @testset "Cholesky decomposition" begin
        X = [  4  12 -16
              12  37 -43
             -16 -43  98]
        U = [ 2.0 6.0 -8.0
              0.0 1.0  5.0
              0.0 0.0  3.0]

        BD = BlockDiagonal([X, X])
        C = cholesky(BD)
        @test C isa Cholesky{Float64, <:BlockDiagonal{Float64}}
        @test C.U ≈ cholesky(Matrix(BD)).U
        @test C.U ≈ BlockDiagonal([U, U])
        @test C.L ≈ BlockDiagonal([U', U'])
        @test C.UL ≈ C.U
        @test C.uplo === 'U'
        @test C.info == 0

        M = BlockDiagonal(map(Matrix, blocks(C.L)))
        C = Cholesky(M, 'L', 0)
        @test C.U ≈ cholesky(Matrix(BD)).U
        @test C.U ≈ BlockDiagonal([U, U])
        @test C.L ≈ BlockDiagonal([U', U'])
        @test C.UL ≈ C.L
        @test C.uplo === 'L'
        @test C.info == 0
    end  # Cholesky
    @testset "Singular Value Decomposition" begin
        X = [  4  12 -16
              12  37 -43
             -16 -43  98]
        BD = BlockDiagonal([X, X])

        @testset "full=$full" for full in (true, false)

            @testset "svd_blockwise" begin
                U, S, Vt = svd_blockwise(BD; full=full)
                F = SVD(U, S, Vt)
                @test BD ≈ F.U * Diagonal(F.S) * F.Vt

                # Matrices should be BlockDiagonal
                @test F isa SVD{Float64, Float64, <:BlockDiagonal{Float64}}
                @test F.U isa BlockDiagonal
                @test F.V isa BlockDiagonal
                @test F.Vt isa BlockDiagonal

                # Should have same values, but not sorted so as to keep BlockDiagonal structure
                F_ = svd(Matrix(BD), full=full)
                for fname in fieldnames(SVD)
                    @test sort(vec(getfield(F, fname))) ≈ sort(vec(getfield(F_, fname)))
                end
                # Singular values should be block-wise
                s = svd(X).S
                @test F.S == vcat(s, s)
            end

            @testset "svd" begin
                F = svd(BD; full=full)
                F_ = svd(Matrix(BD), full=full)

                @test F isa SVD
                @test BD ≈ F.U * Diagonal(F.S) * F.Vt

                @test F == F_
                for fname in fieldnames(SVD)
                    @test getfield(F, fname) ≈ getfield(F_, fname)
                end

                # Singular values should be sorted in descending order
                @test F.S == sort(F.S, rev=true)
            end
        end
    end  # SVD
    @testset "Left division" begin
        N1 = 20
        N2 = 8
        N3 = 5
        A = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2)])
        B = BlockDiagonal([rand(rng, N1, N2), rand(rng, N3, N1)])
        x = rand(rng, N1 + N2)
        y = rand(rng, N1 + N3)
        @test A \ x ≈ inv(A) * x
        @test B \ y ≈ Matrix(B) \ y
    end
end


text:
###################################
## Example code
## for Bayesian Poisson HMM

using PyPlot, PyCall
using HDF5, JLD
@pyimport matplotlib.gridspec as gspec

push!(LOAD_PATH,".")
import PoissonHMM
import PoissonMixtureModel

"""
Simple comparison between HMM and mixture model.
"""
function test_comparison()
    #########################
    ## load data
    file_name = "../data/timeseries.jld"
    X = load(file_name)["obs"]
    N = length(X)

    #########################
    ## Poison HMM

    ## set model
    K = 2 #  number of mixture components
    alpha_phi = 10.0 * ones(K)
    alpha_A = 100.0 * eye(K) + 1.0*ones(K, K)
    cmp = [PoissonHMM.Gam(1.0, 0.01), PoissonHMM.Gam(1.0, 0.01)]
    bhmm = PoissonHMM.BHMM(K, alpha_phi, alpha_A, cmp)

    ## inference
    max_iter = 100
    tic()
    Z_est_hmm, post_bhmm = PoissonHMM.learn_VI(X, bhmm, max_iter)
    toc()

    #########################
    ## Poison Mixture Model

    ## set model
    K = 2 #  number of mixture components
    alpha_phi = 10.0 * ones(K)
    cmp = [PoissonMixtureModel.Gam([1.0], 0.01), PoissonMixtureModel.Gam([1.0], 0.01)]
    bpmm = PoissonMixtureModel.BPMM(1, K, alpha_phi, cmp)

    ## inference
    max_iter = 100
    tic()
    Z_est_pmm, post_bpmm = PoissonMixtureModel.learn_VI(reshape(X, 1, N), bpmm, max_iter)
    toc()

    #########################
    ## Compare results
    figure("Hidden Markov Model vs Mixture Model")
    subplot(3,1,1);plot(X);ylabel("data")
    subplot(3,1,2);fill_between(1:N, reshape(Z_est_hmm[1,:]', N), zeros(N));ylim([0.0, 1.0]);ylabel("S (PHMM)")
    subplot(3,1,3);fill_between(1:N, reshape(Z_est_pmm[1,:]', N), zeros(N));ylim([0.0, 1.0]);ylabel("S (PMM)")
    show()
end

test_comparison()


text:
# Chapter 10. Arrays
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#chap10

# An Array is a Sequence
[10, 20, 30]
["crunchy frog", "ram bladder"]

# heterogeneous array
["spam", 2.0, 5, [10, 20]]

# [] empty array
typeof([])

cheeses = ["Cheddar", "Edam", "Gouda"];
numbers = [42, 123]
typeof(cheeses)

# Arrays Are Mutable
numbers[2] = 5
println(numbers)

# IN operator
"Edam" ∈ cheeses
"Brie" ∈ cheeses

# Traversing an Array
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_traversing_an_array
for cheese in cheeses
    println(cheese)
end

# with index
for i in eachindex(cheeses)
    println(cheeses[i])
end

# Array Slices
t = ['a', 'b', 'c', 'd', 'e', 'f']
print(t[1:3])
print(t[3:end])

t[2:3] = ['x', 'y']
println(t)

# Array Library

# push!
t = ['a', 'b', 'c']
push!(t, 'd')
println(t)

# append!
t1 = ['a', 'b', 'c']
t2 = ['d', 'e']
append!(t1, t2)
println(t1)

# sort! & sort
t = ['d', 'c', 'e', 'b', 'a']
tsorted = sort(t)
println(t)       # ['d', 'c', 'e', 'b', 'a']
println(tsorted) # ['a', 'b', 'c', 'd', 'e']

sort!(t)
println(t) # ['a', 'b', 'c', 'd', 'e']

# ----------------------
# Map, Filter and Reduce
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_map_filter_and_reduce
# ----------------------

# Reduce: An operation like this that combines a sequence of elements into
#         a single value is sometimes called a reduce operation.

# Instead of doing this
function addall(t)
    sum = 0
    for v in t
        sum += v
    end
    return sum
end
addall(1:0.5:3)

# sum(itr)
sum(1:0.5:3)

# Map: “maps” a function (in this case uppercase) onto each of the elements in a sequence.
"""
    capitalizeall(itr)

Capitalize each element of itr.
"""
function capitalizeall(itr) # ::1-D array
    res = []
    for v in itr
        push!(res, uppercase(v))
    end

    return res
end
capitalizeall(['a', 'b'])

# Filter: Selects some of the elements and filters out the others.
"""
    onlyupper(itr)

Returns only upper case elements of itr
"""
function onlyupper(itr)
    ret = []
    for v in itr
        if v == uppercase(v)
            push!(ret, v)
        end
    end

    return ret
end
line = onlyupper(['a', 'B', 'C'])
println(line) # Any['B', 'C']

# ----------
# Dot Syntax
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_dot_syntax
# ----------

# use on operator
t = [1, 2, 3]
println(t.^2) # [1, 4, 9]

# use on function
t = uppercase.(["abc", "def", "ghi"])
println(t) # ["ABC", "DEF", "GHI"]

# -----------------------------
# Deleting (Inserting) Elements
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_deleting_inserting_elements
# -----------------------------

# splice!
t = ['a', 'b', 'c']
splice!(t, 2)
println(t) # ['a', 'c']

# pop!
t = ['a', 'b', 'c']
pop!(t)
println(t) # ['a', 'b']

# popfirst!
t = ['a', 'b', 'c']
popfirst!(t)
println(t) # ['b', 'c']

# push!
t = ['a', 'b']
push!(t, 'd')
println(t) # ['a', 'b', 'd']

# pushfirst!
t = ['a', 'b']
pushfirst!(t, 'X')
println(t) # ['X', 'a', 'b']

# deleteat!
t = ['a', 'b']
deleteat!(t, 2)
println(t) # ['a']

# insert!
t = ['a', 'b', 'c']
insert!(t, 2, 'X')
println(t) # ['a', 'X', 'b', 'c']

# ---------------
# Array Arguments
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_array_arguments
# ---------------

# vcat
t = vcat(1,2)
println(t)

c = ([1 2 3], [4 5 6])
t = vcat(c...) # splat operator, see help for ...
println(t)


text:
using Plots
using Interpolations
include("../src/functions.jl")

pyplot()

function test_bilinear()
    # Same as example on wikipedia
    # https://en.wikipedia.org/wiki/Bilinear_interpolation

    # Unit square
    x_bounds = (0, 1)
    y_bounds = (0, 1)

    #   [f(x1, y1)    f(x1, y2)
    #    f(x2, y1)    f(x2, y2)]

    values = [0. 1.
              1. 0.5]

    x = collect(LinRange(0,1,100))
    y = collect(LinRange(0,1,100))

    F = Matrix{Float64}(undef, (length(x), length(y)))

    for i in eachindex(x)
        for j in eachindex(y)
            F[i, j] = bilinear(x[i], y[j], x_bounds, y_bounds, values)
        end
    end

    heatmap(x, y, F, c=:jet, dpi=500, title="Test Bilinear",
            xaxis="x", yaxis="y", aspect_ratio=:equal)

    savefig("bilinear_test")
end

function test_trilinear()
    # Unit cube
    x_bounds = [0.0, 1.0]
    y_bounds = [0.0, 1.0]
    z_bounds = [0.0, 1.0]

    #   [[f(x1, y1)    f(x1, y2)
    #    f(x2, y1)    f(x2, y2)]

    values = Array{Float64, 3}(undef, (2, 2, 2))
    values[:, 1, 1] .= 0.0
    values[:, 1, 2] .= 1.0
    values[:, 2, 1] .= 1.0
    values[:, 2, 2] .= 0.5

    x = collect(LinRange(0.1, 0.99, 100))
    y = collect(LinRange(0.1, 0.99, 100))
    z = collect(LinRange(0.1, 0.99, 100))

    F = Array{Float64, 3}(undef, (length(x), length(y), length(z)))

    for k in eachindex(z)
        for i in eachindex(x)
            for j in eachindex(y)
                F[k, i, j] = trilinear(z[k], x[i], y[j], z_bounds, x_bounds, y_bounds, values)
            end
        end
    end

    for k in [1, 50, 100]
        heatmap(x, y, F[k, :, :], c=:jet, dpi=500, title="Test Trilinear $k",
                xaxis="x", yaxis="y")

        savefig("trilinear_test_$k")
    end
end

function trilinear(z_mrk::Float64, x_mrk::Float64, y_mrk::Float64,
                   z::Vector{Float64}, x::Vector{Float64}, y::Vector{Float64},
                   vals::AbstractArray)
    # Returns the index of the first value in a greater than or equal to x
    # Subtract by 1 to get coordinate of lower corner
    idz = searchsortedfirst(z, z_mrk) - 1
    idx = searchsortedfirst(x, x_mrk) - 1
    idy = searchsortedfirst(y, y_mrk) - 1

    # bounding corner coordinates
    z0 = z[idz]; z1 = z[idz+1]
    x0 = x[idx]; x1 = x[idx+1]
    y0 = y[idy]; y1 = y[idy+1]

    # difference between coordinates and interpolation point
    x_d = (x_mrk - x0)/(x1 - x0)
    y_d = (y_mrk - y0)/(y1 - y0)
    z_d = (z_mrk - z0)/(z1 - z0)

    # values at each corner (z is first index in data array)
    c000 = vals[idz, idx, idy]
    c010 = vals[idz, idx, idy+1]
    c100 = vals[idz, idx+1, idy]
    c110 = vals[idz, idx+1, idy+1]
    c001 = vals[idz+1, idx, idy]
    c011 = vals[idz+1, idx, idy+1]
    c101 = vals[idz+1, idx+1, idy]
    c111 = vals[idz+1, idx+1, idy+1]

    # interpolate in x direction
    c00 = c000*(1 - x_d) + c100*x_d
    c01 = c001*(1 - x_d) + c101*x_d
    c10 = c010*(1 - x_d) + c110*x_d
    c11 = c011*(1 - x_d) + c111*x_d

    # interpolate in y direction
    c0 = c00*(1 - y_d) + c10*y_d
    c1 = c01*(1 - y_d) + c11*y_d

    # intepolate in z direction
    c = c0*(1 - z_d) + c1*z_d
    return c
end

test_bilinear()
test_trilinear()

max_stars_repo_name: meudnaes/VoronoiRT

max_issues_repo_name: meudnaes/VoronoiRT

max_forks_repo_name: meudnaes/VoronoiRT


text:
"""
    abundances(at::AbstractAbundanceTable)

Get the underlying sparse matrix of an `AbstractAbundanceTable`.
Note that this does not copy - any modifications to this matrix
will update the parent.
"""
const abundances = EcoBase.occurrences

const nfeatures = EcoBase.nthings

"""
    featurenames(at::AbstractAbundanceTable)

Get a vector of feature names from `at`, equivalent to `name.(features(at))`
"""
const featurenames = EcoBase.thingnames
const getfeature = EcoBase.thingoccurrences

const nsamples = EcoBase.nplaces

"""
    samplenames(at::AbstractAbundanceTable)

Get a vector of sample names from `at`, equivalent to `name.(samples(at))`
"""
const samplenames = EcoBase.placenames
const getsample = EcoBase.placeoccurrences

abstract type AbstractFeature <: EcoBase.AbstractThings end
abstract type AbstractSample <: EcoBase.AbstractPlaces{Nothing} end


text:
using Test, RandomProjectionTree

using Logging
using Base.CoreLogging

logger = ConsoleLogger(stdout, CoreLogging.Debug)
global_logger(logger)

include("testrandom.jl")
@testset "random" begin
    @test check_for_random_data()
end


text:
function noop_recommender(request; environment=nothing)
    environment == nothing && @error "No search environment provided for noop recommender."
    @debug "Noop recommender, returning empty result..."
    empty_recommendation = build_result_from_ids(environment.dbdata,
                                                 Int[],
                                                 environment.id_key,
                                                 make_id();
                                                 id_key=environment.id_key,
                                                 max_matches=request.max_matches)
    return [empty_recommendation]
end


text:
@json_read mutable struct WorkspaceEditCapabilities
    documentChanges::Union{Nothing,Bool}
end

@json_read mutable struct Capabilities
    dynamicRegistration::Union{Nothing,Bool}
end

@json_read mutable struct SymbolKindCapabilities
    valueSet::Union{Nothing,Vector{Int}}
end

@json_read mutable struct SymbolCapabilities
    dynamicRegistration::Union{Nothing,Bool}
    symbolKind::Union{Nothing,SymbolKindCapabilities}
end

@json_read mutable struct SynchronizationCapabilities
    dynamicRegistration::Union{Nothing,Bool}
    willSave::Union{Nothing,Bool}
    willSaveWaitUntil::Union{Nothing,Bool}
    didSave::Union{Nothing,Bool}
end

@json_read mutable struct CompletionItemCapabilities
    snippetSupport::Union{Nothing,Bool}
    commitCharactersSupport::Union{Nothing,Bool}
    documentationFormat::Union{Nothing,Vector{String}}
    deprecatedSupport::Union{Nothing,Bool}
end

@json_read mutable struct CompletionItemKindCapabilities
    valueSet::Union{Nothing,Vector{Int}}
end

@json_read mutable struct CompletionCapabilities
    dynamicRegistration::Union{Nothing,Bool}
    completionItem::CompletionItemCapabilities
    completionItemKind::Union{Nothing,CompletionItemKindCapabilities}
    contextSupport::Union{Nothing,Bool}
end

@json_read mutable struct HoverCapabilities
    dynamicRegistration::Union{Nothing,Bool}
    contentFormat::Union{Nothing,Vector{String}}
end

@json_read mutable struct SignatureInformationCapabilities
    documentationFormat::Union{Nothing,Vector{String}}
end

@json_read mutable struct SignatureCapabilities
    dynamicRegistration::Union{Nothing,Bool}
    signatureInformation::Union{Nothing,SignatureInformationCapabilities}
end

@json_read mutable struct DocumentSymbolCapabilities
    dynamicRegistration::Union{Nothing,Bool}
    symbolKind::Union{Nothing,SymbolKindCapabilities}
end

@json_read mutable struct CodeActionKindCapabilities
    valueSet::Vector{String}
end

@json_read mutable struct CodeActionLiteralCapabilities
    codeActionKind::CodeActionKindCapabilities
end

@json_read mutable struct CodeActionCapabilities
    dynamicRegistration::Union{Nothing,Bool}
    codeActionLiteralSupport::Union{Nothing,CodeActionLiteralCapabilities}
end

@json_read mutable struct PublishDiagnosticsCapabilities
    relatedInformation::Union{Nothing,Bool}
end

@json_read mutable struct WorkspaceClientCapabilities
    applyEdit::Union{Nothing,Bool}
    workspaceEdit::Union{Nothing,WorkspaceEditCapabilities}
    didChangeConfiguration::Union{Nothing,Capabilities}
    didChangeWatchedFiles::Union{Nothing,Capabilities}
    symbol::Union{Nothing,SymbolCapabilities}
    executeCommand::Union{Nothing,Capabilities}
    workspaceFolders::Union{Nothing,Bool}
    configuration::Union{Nothing,Bool}
end

@json_read mutable struct TextDocumentClientCapabilities
    synchronization::Union{Nothing,SynchronizationCapabilities}
    completion::Union{Nothing,CompletionCapabilities}
    hover::Union{Nothing,HoverCapabilities}
    signatureHelp::Union{Nothing,SignatureCapabilities}
    references::Union{Nothing,Capabilities}
    documentHighlight::Union{Nothing,Capabilities}
    documentSymbol::Union{Nothing,DocumentSymbolCapabilities}
    formatting::Union{Nothing,Capabilities}
    rangeFormatting::Union{Nothing,Capabilities}
    onTypeFormatting::Union{Nothing,Capabilities}
    definition::Union{Nothing,Capabilities}
    typeDefinition::Union{Nothing,Capabilities}
    implementation::Union{Nothing,Capabilities}
    codeAction::Union{Nothing,CodeActionCapabilities}
    CodeLens::Union{Nothing,Capabilities}
    documentLink::Union{Nothing,Capabilities}
    colorProvider::Union{Nothing,Capabilities}
    rename::Union{Nothing,Capabilities}
    publishDiagnostics::Union{Nothing,PublishDiagnosticsCapabilities}
end

@json_read mutable struct ClientCapabilities
    workspace::Union{Nothing,WorkspaceClientCapabilities}
    textDocument::Union{Nothing,TextDocumentClientCapabilities}
    experimental::Union{Nothing,Any}
end

@json_read mutable struct InitializeParams
    processId::Union{Nothing,Int}
    rootPath::Union{Nothing,DocumentUri}
    rootUri::Union{Nothing,DocumentUri}
    initializationOptions::Union{Nothing,Any}
    capabilities::ClientCapabilities
    trace::Union{Nothing,String}
    workspaceFolders::Union{Nothing,Vector{WorkspaceFolder}}
end

# Server Response
mutable struct CompletionOptions
    resolveProvider::Bool
    triggerCharacters::Vector{String}
end

mutable struct SignatureHelpOptions
    triggerCharacters::Vector{String}
end

mutable struct CodeLensOptions
    resolveProvider::Bool
end

mutable struct DocumentOnTypeFormattingOptions
    firstTriggerCharacter::String
    moreTriggerCharacters::Vector{String}
end

mutable struct DocumentLinkOptions
    resolveProvider::Bool
end

mutable struct ExecuteCommandOptions
    commands::Vector{String}
end

mutable struct SaveOptions
    includeText::Bool
end

mutable struct TextDocumentSyncOptions
    openClose::Bool
    change::Int
    willSave::Bool
    willSaveWaitUntil::Bool
    save::SaveOptions
end

mutable struct StaticRegistrationOptions
    id::String
end

mutable struct WorkspaceFoldersOptions
    supported::Bool
    changeNotifications::Union{Bool,String}
end

mutable struct WorkspaceOptions
    workspaceFolders::WorkspaceFoldersOptions
end

mutable struct ServerCapabilities
    textDocumentSync::Union{TextDocumentSyncOptions,Int}
    hoverProvider::Bool
    completionProvider::CompletionOptions
    signatureHelpProvider::SignatureHelpOptions
    definitionProvider::Bool
    typeDefinitionProvider::Bool
    implementationProvider::Bool
    referencesProvider::Bool
    documentHighlightProvider::Bool
    documentSymbolProvider::Bool
    workspaceSymbolProvider::Bool
    codeActionProvider::Bool
    # codeLensProvider::CodeLensOptions
    documentFormattingProvider::Bool
    documentRangeFormattingProvider::Bool
    # documentOnTypeFormattingProvider::DocumentOnTypeFormattingOptions
    renameProvider::Bool
    # documentLinkProvider::DocumentLinkOptions
    colorProvider::Bool
    executeCommandProvider::ExecuteCommandOptions
    workspace::WorkspaceOptions
    experimental::Any
end

mutable struct InitializeResult
    capabilities::ServerCapabilities
end

mutable struct Registration
    id::String
    method::String
    registerOptions::Any
end

mutable struct RegistrationParams
    registrations::Vector{Registration}
end

mutable struct Unregistration
    id::String
    method::String
end

mutable struct UnregistrationParams
    unregistrations::Vector{Unregistration}
end

mutable struct DidChangeConfiguration
    settings::Any
end

mutable struct ConfigurationItem
    scopeUri::Union{Nothing,String}
    section::Union{Nothing,String}
end

mutable struct ConfigurationParams
    items::Vector{ConfigurationItem}
end


text:
using Test
using PackedVectorsOfVectors

using Documenter
DocMeta.setdocmeta!(
    PackedVectorsOfVectors,
    :DocTestSetup,
    quote
        using PackedVectorsOfVectors
        using PackedVectorsOfVectors: cumsum1
    end
)
doctest(PackedVectorsOfVectors, testset="docs")

@testset "pack" begin
    pv = pack([[1,2],[3]])
    @test pv.p == [1,3,4]
    @test pv.v == [1,2,3]

    pv = pack([])
    @test pv.p == [1]
    @test pv.v == []

    pv = pack([[]])
    @test pv.p == [1,1]
    @test pv.v == []

    pv = pack([[],[]])
    @test pv.p == [1,1,1]
    @test pv.v == []

    pv = pack([[1],[]])
    @test pv.p == [1,2,2]
    @test pv.v == [1]

    pv = pack([[],[1]])
    @test pv.p == [1,1,2]
    @test pv.v == [1]
end

@testset "allocate_packed" begin
    pv = allocate_packed(Int, undef, [1,2,3])
    @test pv.p == [1,2,4,7]
    @test length(pv.v) == 6
end

@testset "$pf" for (f,pf) in ((zeros,packed_zeros),(ones,packed_ones))
    n = [2,3]
    pv = pf(n)
    @test eltype(eltype(pv)) == Float64
    for k = 1:length(n)
        @test pv[k] == f(n[k])
    end

    n = [2,3]
    pv = pf(Int,n)
    @test eltype(eltype(pv)) == Int
    for k = 1:length(n)
        @test pv[k] == f(Int,n[k])
    end
end

@testset "packed_fill" begin
    x = 42
    n = [2,3]
    pv = packed_fill(x,n)
    @test eltype(eltype(pv)) == typeof(x)
    for k = 1:length(n)
        @test pv[k] == fill(x,n[k])
    end
end

@testset "vector interface" begin
    vv = [[1,2],[3]]
    @test collect(pack(vv)) == vv
end

@testset "transpose" begin
    @test collect(transpose(pack([[1,2],[1]]))) == pack([[1,2],[1]])
    @test collect(transpose(pack([[1,2],[2]]))) == pack([[1],[1,2]])
    @test collect(transpose(pack([[1,2],[3]]))) == pack([[1],[1],[2]])

    @test collect(transpose(pack(Vector{Int}[]))) == pack([])
    @test collect(transpose(pack(Vector{Int}[[]]))) == pack([])
    @test collect(transpose(pack(Vector{Int}[[1],[]]))) == pack([[1]])
    @test collect(transpose(pack(Vector{Int}[[],[2]]))) == pack([[],[2]])
end


text:
#Linear Kalman filter example
#This example closely follows that given on "Greg Czerniak's Website". Namely
#the voltage example on this page: http://greg.czerniak.info/guides/kalman1/

#Let's import the modules required to execute the Kalman Filter
using StateSpace
using Distributions
using DataFrames
using Gadfly
using Colors

#To use the Kalman Filter we need three things:
#1) Some observations
#2) The parameters of the Kalman Filter model describing the system process and
#the observation model.
#3) An initial guess of the state.
#
#The following code is split into those three sections to begin with. We finish
#by performing the Kalman Filter algorithm and finally displaying the results.

################################################################################
#Section: Generate noisy Observations
#-------------------------------------------------------------------------------
#Let's assume that the true voltage value does not change throughout the
#experiment. So we can give the true voltage a value
true_voltage = 1.25

#However our measurements of the true voltage are noisy. Here we assume that the
#noise is normally distributed with zero mean but non-zero variance:
measurement_noise_variance = 0.1

#Given this information we can generate some noisy data (observations). We'll
#also choose the number of observations
number_of_observations = 60
observations = randn(number_of_observations) * sqrt(measurement_noise_variance) + true_voltage

#Note that the observation list has to be transposed because each column of the
#matrix is considered as a new observation.
observations = observations'
#End Section: Generate noisy Observations
################################################################################

################################################################################
#Section: Describe Kalman Filter parameters
#-------------------------------------------------------------------------------
#Now that we have some observations we need to describe the Kalman Filter
#parameters:
#Process matrix - Note that since this is a one-dimensional system I've given
#a floating point value rather than explicitly declaring a matrix with one
#element, although you can do it that way. The value of 1.0 here means that the
#process model suggests that the previous state is the same as the next state.
#i.e. the voltage doesn't change.
process_matrix = 1.0

#Process Covariance - Since it's a single value, the covariance is the same as
#variance. We've made it small here because we are fairly sure of our process
#model that the state doesn't change. If we were less sure that the model was
#correct then we would increase this value.
process_covariance = 0.00001

#Observation matrix - Since we measure the voltage directly (we don't measure
#some quantity related by a function - i.e. we don't measure twice the voltage)
#the observation matrix is equal to 1.0
observation_matrix = 1.0

#Observation Covariance - we set the variance of the measurement noise earlier
#so we'll give it the same value.
observation_covariance = measurement_noise_variance

#Now we can create our Linear State Space Model
linSSM = LinearGaussianSSM(process_matrix, process_covariance, observation_matrix, observation_covariance)
#End Section: Describe Kalman Filter parameters
################################################################################

################################################################################
#Section: Set initial guess of the state
#-------------------------------------------------------------------------------
#We'll begin with a guess that's quite off. We'll assume an initial mean state
#of 3 (i.e. we think we have 3 volts) and give the variance as 1
initial_guess = MvNormal([3.0], [1.0])
#End Section: Set initial guess of the state
################################################################################

################################################################################
#Section: Execute Kalman Filter
#-------------------------------------------------------------------------------
#Now that we have some noisy observations, the Kalman filter parameters and an
#intial guess for the state of the system, we can run the Kalman Filter
filtered_state = filter(linSSM, observations, initial_guess)
#End Section: Execute Kalman Filter
################################################################################

################################################################################
#Section: Plot Filtered results
#-------------------------------------------------------------------------------
#Here we are plotting the filtered results with Gadfly. See the Gadfly
#documentation for information about how plotting works if you are unfamiliar.
#Website: http://gadflyjl.org/

x_data = 1:number_of_observations
state_array = Vector{Float64}(number_of_observations+1)
confidence_array = Vector{Float64}(number_of_observations+1)
state_array[1] = initial_guess.μ[1]
confidence_array[1] = 2*sqrt(initial_guess.Σ.diag[1])
for i in x_data
    current_state = filtered_state.state[i]
    state_array[i+1] = current_state.μ[1]
    confidence_array[i+1] = 2*sqrt(current_state.Σ.mat[1])
end
df_fs = DataFrame(
    x = [0;x_data],
    y = state_array,
    ymin = state_array - confidence_array,
    ymax = state_array + confidence_array,
    f = "Filtered values"
    )

n = 3
getColors = distinguishable_colors(n, Color[LCHab(70, 60, 240)],
                                   transform=c -> deuteranopic(c, 0.5),
                                   lchoices=Float64[65, 70, 75, 80],
                                   cchoices=Float64[0, 50, 60, 70],
                                   hchoices=linspace(0, 330, 24))
filtered_state_plot = plot(
    layer(x=x_data, y=filtered_state.observations, Geom.point, Theme(default_color=getColors[2])),
    layer(x=[0;x_data], y=ones(number_of_observations+1)*true_voltage, Geom.line, Theme(default_color=getColors[3])),
    layer(df_fs, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon),
    Guide.xlabel("Measurement Number"), Guide.ylabel("Voltage (Volts)"),
    Guide.manual_color_key("Colour Key",["Filtered Estimate", "Measurements","True Value "],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Linear Kalman Filter Example")
    )
display(filtered_state_plot)
#End Section: Plot Filtered results
################################################################################

################################################################################
#Section: Execute Kalman Smoother
#-------------------------------------------------------------------------------
#Everything can be performed much better in hindsight and that's essentially
#where the Kalman smoother comes in. Given all of the observations and the
#filtered state (including the Kalman filter parameters, the smoother attempts
#to give better estimates of the system's state. Here we perform smoothing on
#the filtered data.
smoothed_state = smooth(linSSM, filtered_state)
#End Section: Execute Kalman Filter
################################################################################

################################################################################
#Section: Plot Smoothed Results
#-------------------------------------------------------------------------------
#Here we are plotting the smoothed results with Gadfly. See the Gadfly
#documentation for information about how plotting works if you are unfamiliar.
#Website: http://gadflyjl.org/
state_array = Vector{Float64}(number_of_observations)
confidence_array = Vector{Float64}(number_of_observations)
for i in x_data
    current_state = smoothed_state.state[i]
    state_array[i] = current_state.μ[1]
    confidence_array[i] = 2*sqrt(current_state.Σ.mat[1])
end
df_ss = DataFrame(
    x = x_data,
    y = state_array,
    ymin = state_array - confidence_array,
    ymax = state_array + confidence_array,
    f = "Filtered values"
    )

smoothed_state_plot = plot(
    layer(x=x_data, y=smoothed_state.observations, Geom.point, Theme(default_color=getColors[2])),
    layer(x=x_data, y=ones(number_of_observations)*true_voltage, Geom.line, Theme(default_color=getColors[3])),
    layer(df_ss, x=:x, y=:y, ymin=:ymin, ymax=:ymax, Geom.line, Geom.ribbon),
    Guide.xlabel("Measurement Number"), Guide.ylabel("Voltage (Volts)"),
    Guide.manual_color_key("Colour Key",["Smoothed Estimate", "Measurements","True Value "],[getColors[1],getColors[2],getColors[3]]),
    Guide.title("Linear Kalman Smoother Example")
    )
display(smoothed_state_plot)
#End Section: Plot Filtered results
################################################################################


text:
using ADCME
function buildops(dirname)
    PWD = pwd()
    cd(joinpath("$(@__DIR__)/CustomOp/", dirname))
    rm("build", force=true, recursive=true)
    if !isdir("build")
        mkdir("build")
        cd("build")
        ADCME.cmake()
        run(`$(ADCME.MAKE) -j`)
    end
    cd(PWD)
end

buildops("FintComp")
buildops("SymOp")
buildops("OrthotropicOp")
buildops("SPDOp")
buildops("CholOp")
buildops("CholOrthOp")

# # build ADLaw
# PWD = pwd()
# cd(joinpath("$(@__DIR__)/../", "deps/CustomOp/ADLaw"))
# if !isdir("Adept-2")
#     run(`sh download_lib.sh`)
# end
# if !isdir("build")
#     mkdir("build")
#     cd("build")
#     ADCME.cmake()
#     run(`$(ADCME.MAKE) -j`)
# end
# cd(PWD)


text:
# Import libraries.
using Turing, Turing: Variational
using Distributions
using AbstractGPs, KernelFunctions
import Random
import LinearAlgebra

# NOTE: Import other libraries ...

# NOTE: Read data ...

# Define a kernel.
sqexpkernel(alpha::Real, rho::Real) =
    alpha^2 * transform(SqExponentialKernel(), 1/(rho*sqrt(2)))

# Define model.
@model GPRegression(y, X) = begin
    # Priors.
    alpha ~ LogNormal(0.0, 0.1)
    rho ~ LogNormal(0.0, 1.0)
    sigma ~ LogNormal(0.0, 1.0)

    # Covariance function.
    kernel = sqexpkernel(alpha, rho)

    # GP (implicit zero-mean).
    gp = GP(kernel)

    # Sampling Distribution (MvNormal likelihood).
    y ~ gp(X, sigma^2 + 1e-6)  # add 1e-6 for numerical stability.
end;

# Set random number generator seed.
Random.seed!(0)

# Model creation.
m = GPRegression(y, X)

# Fit via ADVI.
q0 = Variational.meanfield(m)  # initialize variational distribution (optional)
num_elbo_samples, max_iters = (1, 2000)
# Run optimizer.
@time q = vi(m, ADVI(num_elbo_samples, max_iters), q0,
             optimizer=Flux.ADAM(1e-1));

# Fit via HMC.
burn = 1000
nsamples = 1000
@time chain = sample(m, HMC(0.01, 100), burn + nsamples)  # start sampling.

# Fit via NUTS.
@time chain = begin
    nsamples = 1000  # number of MCMC samples
    nadapt = 1000  # number of iterations to adapt tuning parameters in NUTS
    iterations = nsamples + nadapt
    target_accept_ratio = 0.8

    sample(m, NUTS(nadapt, target_accept_ratio, max_depth=10), iterations);
end

max_stars_repo_name: luiarthur/TuringBnpBenchmarks

max_issues_repo_name: luiarthur/TuringBnpBenchmarks

max_forks_repo_name: luiarthur/TuringBnpBenchmarks


text:
# Generate quadrature information for B-spline product.
# Returns weights and nodes for integration in [-1, 1].
#
# See https://en.wikipedia.org/wiki/Gaussian_quadrature.
#
# Some notes:
#
# - On each interval between two neighbouring knots, each B-spline is a
#   polynomial of degree (k - 1). Hence, the product of two B-splines has degree
#   (2k - 2).
#
# - The Gauss--Legendre quadrature rule is exact for polynomials of degree
#   <= (2n - 1), where n is the number of nodes and weights.
#
# - Conclusion: on each knot interval, `k` nodes should be enough to get an
#   exact integral. (I verified that results don't change when using more than
#   `k` nodes.)
#
# Here, p is the polynomial order (p = 2k - 2 for the product of two B-splines).
function _quadrature_prod(::Val{p}) where {p}
    n = cld(p + 1, 2)
    _gausslegendre(Val(n))
end

# Compile-time computation of Gauss--Legendre nodes and weights.
# If the @generated branch is taken, then the computation is done at compile
# time, making sure that no runtime allocations are performed.
function _gausslegendre(::Val{n}) where {n}
    if @generated
        vecs = _gausslegendre_impl(Val(n))
        :( $vecs )
    else
        _gausslegendre_impl(Val(n))
    end
end

function _gausslegendre_impl(::Val{n}) where {n}
    data = gausslegendre(n)
    map(SVector{n}, data)
end

# Metric for integration on interval [a, b].
# This is to transform from integration on the interval [-1, 1], in which
# quadratures are defined.
# See https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval
struct QuadratureMetric{T}
    α :: T
    β :: T

    function QuadratureMetric(a::T, b::T) where {T}
        α = (b - a) / 2
        β = (a + b) / 2
        new{T}(α, β)
    end
end

# Apply metric to normalised coordinate (x ∈ [-1, 1]).
Base.:*(M::QuadratureMetric, x::Real) = M.α * x + M.β
Broadcast.broadcastable(M::QuadratureMetric) = Ref(M)

# Evaluate elements of basis `B` (given by indices `is`) at points `xs`.
# The length of `xs` is assumed static.
# The length of `is` is generally equal to the B-spline order, but may me
# smaller near the boundaries (this is not significant if knots are "augmented",
# as is the default).
# TODO implement evaluation of all B-splines at once (should be much faster...)
function eval_basis_functions(B, is, xs, args...)
    k = order(B)
    @assert length(is) ≤ k
    ntuple(Val(k)) do n
        # In general, `is` will have length equal `k`.
        # If its length is less than `k` (may happen near boundaries, if knots
        # are not "augmented"), we repeat values for the first B-spline, just
        # for type stability concerns. These values are never used.
        i = n > length(is) ? is[1] : is[n]
        map(xs) do x
            B[i](x, args...)
        end
    end
end


text:
# if it's one symbol, set the guide and return the column
function handle_dfs(df::AbstractDataFrame, d::KW, letter, sym::Symbol)
    get!(d, Symbol(letter * "guide"), string(sym))
    collect(df[sym])
end

# if it's an array of symbols, set the labels and return a Vector{Any} of columns
function handle_dfs(df::AbstractDataFrame, d::KW, letter, syms::AbstractArray{Symbol})
    get!(d, :label, reshape(syms, 1, length(syms)))
    vec(Any[collect(df[s]) for s in syms])
end

# for anything else, no-op
function handle_dfs(df::AbstractDataFrame, d::KW, letter, anything)
    anything
end

# handle grouping by DataFrame column
function Plots.extractGroupArgs(group::Symbol, df::AbstractDataFrame, args...)
    Plots.extractGroupArgs(collect(df[group]))
end

# allows the passing of expressions including DataFrame columns as symbols
function processExpr!(expr, df)
    arg = map(expr.args) do x
        isa(x, Expr) && return processExpr!(x, df)

        if isa(x, QuoteNode) && isa(x.value, Symbol)
            return :(collect($(df)[$x]))
        end
        x
    end
    expr.args = arg
    return expr
end

# if a DataFrame is the first arg, lets swap symbols out for columns
@recipe function f(df::AbstractDataFrame, args...)
    # if any of these attributes are symbols, swap out for the df column
    for k in (:fillrange, :line_z, :marker_z, :markersize, :ribbon, :weights, :xerror, :yerror, :hover)
        if haskey(d, k)
            if isa(d[k], Expr)
                d[k] = eval(processExpr!(d[k], df))
            end

            if isa(d[k], Symbol)
                d[k] = collect(df[d[k]])
            end
        end
    end

    # return a list of new arguments
    tuple(Any[handle_dfs(df, d, (i==1 ? "x" : i==2 ? "y" : "z"), arg) for (i,arg) in enumerate(args)]...)
end


text:
# __BEGIN_LICENSE__
#
# ThreeDeconv.jl
#
# Copyright (c) 2018, Stanford University
#
# All rights reserved.
#
# Redistribution and use in source and binary forms for academic and other
# non-commercial purposes with or without modification, are permitted provided
# that the following conditions are met:
#
# * Redistributions of source code, including modified source code, must retain
#   the above copyright notice, this list of conditions and the following
#   disclaimer.
#
# * Redistributions in binary form or a modified form of the source code must
#   reproduce the above copyright notice, this list of conditions and the
#   following disclaimer in the documentation and/or other materials provided with
#   the distribution.
#
# * Neither the name of The Leland Stanford Junior University, any of its
#   trademarks, the names of its employees, nor contributors to the source code
#   may be used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# * Where a modified version of the source code is redistributed publicly in
#   source or binary forms, the modified source code must be published in a freely
#   accessible manner, or otherwise redistributed at no charge to anyone
#   requesting a copy of the modified source code, subject to the same terms as
#   this agreement.
#
# THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF THE LELAND STANFORD JUNIOR
# UNIVERSITY "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE LELAND STANFORD JUNIOR
# UNIVERSITY OR ITS TRUSTEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
# OF SUCH DAMAGE.
#
# __END_LICENSE__

mutable struct ADMM <: Solver
    ρ::Float64
    ϵ_rel_diff::Float64
end

function ADMM(;
        ρ::Real = NaN,
        ϵ_rel_diff::Real = 1e-3) # before: 1e-3, for accuracy: 1e-5
    return ADMM(ρ, ϵ_rel_diff)
end

struct ADMMmetric <: SolverMetric
    iteration::Int
    rel_norm_diff::Float64
    elapsed_time::Float64
end

mutable struct ADMMstate{T<:AbstractArray{Float32,3}} <: SolverState
    x::T
    z::Vector{T}
    λ::Vector{T}
    Kx::Vector{T}
    x_prev::T
end

function ADMMstate(input_shape::NTuple, output_shape::Vector{Tuple})
    if CUDA.functional()
        x = CUDA.zeros(Float32, input_shape)
        z = [CUDA.zeros(Float32, sh) for sh in output_shape]
        λ = [CUDA.zeros(Float32, sh) for sh in output_shape]
        Kx = [CUDA.zeros(Float32, sh) for sh in output_shape]
        # Allocate an array for convergence check
        x_prev = CUDA.zeros(Float32, input_shape)
    else
        x = zeros(Float32, input_shape)
        z = [zeros(Float32, sh) for sh in output_shape]
        λ = [zeros(Float32, sh) for sh in output_shape]
        Kx = [zeros(Float32, sh) for sh in output_shape]
        # Allocate an array for convergence check
        x_prev = zeros(Float32, input_shape)
    end
    return ADMMstate(x, z, λ, Kx, x_prev)
end

function print_header(method::ADMM)
    @printf "Iter     rel_diff      time  \n"
    @printf "------   ---------   --------\n"
end

function Base.show(io::IO, s::ADMMmetric)
    @printf io "%6d   %9.3g   %8.2f\n" s.iteration s.rel_norm_diff s.elapsed_time
end

function initialize_optimizer(method::ADMM,
                              K::StackedLinearOperator,
                              update_x!::Function,
                              update_z!::Function)

    initial_state = ADMMstate(K.input_shape, K.output_shape)

    function update_λ!( Kx::Vector{T}, z::Vector{T}, λ::Vector{T}) where T<:AbstractArray{Float32,3}
        for (Kx_j, z_j, λ_j) in zip(Kx, z, λ)
            λ_j .+= Kx_j .- z_j
        end
    end

    function update_state!(state::ADMMstate)
        state.x_prev .= state.x
        update_x!(state.x, state.z, state.λ)
        K.forward!(state.x, state.Kx)
        update_z!(state.Kx, state.z, state.λ)
        update_λ!(state.Kx, state.z, state.λ)
    end

    function check_state(state::ADMMstate, iter::Int, tic::Float64)
        rel_norm_diff = norm(state.x - state.x_prev) / norm(state.x_prev)
        toc = time()
        return ADMMmetric(iter, rel_norm_diff, toc - tic)
    end

    function check_convergence(metric::ADMMmetric)
        return metric.rel_norm_diff < method.ϵ_rel_diff
    end

    return Optimizer(method, initial_state, update_state!, check_state, check_convergence)
end


text:
module AccurateArithmetic

include("EFT.jl")
using .EFT
export two_sum, two_prod

include("Summation.jl")
using .Summation
export sum_naive, sum_kbn, sum_oro, sum_mixed
export dot_naive, dot_oro, dot_mixed

include("Test.jl")

end # module