text:
# -------------- Main results  --------------

include("ndid_core.jl")
colC = "black"
colD = "#c7254e"
colN = "#ffdd55"
colRe = "#3fb3b2"
colRle = "#8edb7c"

# Bifurcation analysis (increasing df from 0 to 1)

## A- All nodes are fast

### i- NDCRe
param = EcoPar()

### results for all components
bife = bif_analysis2(0.0:0.01:1.0, eco_NDCRe!, fill(0.03, 20), (0.0, 30000.0),
    25000:30000, param)
### bifurcation analysis for the all components
fig_bif(bife, [1 13 17], "figRe_N_30.svg", "df", "N (max/min)", 0, 12.5)
fig_bif(bife, [2 14 18], "figRe_R_30.svg", "df", "Re (max/min)", 0, 2, colRe, 2)
fig_bif(bife, [3 15 19], "figRe_C_30.svg", "df", "C (max/min)", 0, 4)
fig_bif(bife, [4 16 20], "figRe_D_30.svg", "df", "D (max/min)", 0, 2)
### combining 2 results (not used in the last version of the article)
fig_bif2(bife, [1 13 17], [3 15 19], [colN colC], "figRe_p1.svg", "", "",  0, 12.5)
fig_bif2(bife, [1 13 17], [4 16 20], [colN colD], "figRe_p2.svg", "", "",  0, 10.5)
### actual results
fig_bif(bife, [2 14 18], "figRe_Re.svg", "", "", 0, 2, colRe, 3)
barplotmean(bife, [1 13 17], "hist_bife_N.svg", 0, 11, colN)
barplotmean(bife, [4 16 20], "hist_bife_D.svg", 0.170, 0.180, colD)

### ii- NDCReRle
param = EcoPar()
param.a42 = 0.15
param.a43 = 0.05
# NB
# a43=.15; a42=.05 -- loss of edible R
# a43=.15; a42=0.01 -- no loss but edible weakly declines while Ri rises -- this is universal
# a42= .15; a42=0; Redible flatlines and Ri rises

### bifurcation analysis for the all components
bifle = bif_analysis2(0.0:0.01:1.0, eco_NDCReRle!, fill(0.03, 25), (0.0, 30000.0), 25000:30000, param)
### results for all components
fig_bif(bifle, [1 16 21], "figRle_N.svg", "df", "N max/min", 0, 1.5)
fig_bif(bifle, [2 17 22], "figRle_Re.svg", "df", "Re max/min", 0, .6, colRe, 2)
fig_bif(bifle, [3 18 23], "figRle_Rle.svg", "df", "Rle max/min", 0, 1.5)
fig_bif(bifle, [4 19 24], "figRle_C.svg", "df", "C max/min", 0, 2)
fig_bif(bifle, [5 20 25], "figRle_D.svg", "df", "D max/min", 0, 1.5)
### combining 2 or 3 results (not used in the last version of the article)
fig_bif2(bifle, [1 16 21], [5 20 25], [colN colD], "figRle_p2.svg", "", "",  0, 1)
fig_bif3(bifle, ([2 14 18], [4 16 20], [3 18 23]), [colRe colD colRle], "figRle_p2.svg", "", "Max/min",  0, 1.5)
### actual results
fig_bif(bifle, [2 17 22], "figRle_Re.svg", "", "", 0, .6, colRe, 3)
barplotmean(bifle, [1 16 21], "hist_bifle_N.svg", 0, .15, colN)
barplotmean(bifle, [5 20 25], "hist_bifle_D.svg", 0, 1, colD)

## B- All nodes are fast but the hub node that is slow

### i- NDCRe

param = EcoPar()

### bifurcation analysis for the all components
bife_sh = bif_analysis2(0.0:0.01:1.0, eco_NDCRe_sh!, fill(0.03, 20),
    (0.0, 30000.0), 25000:30000, param)

### results for all components
fig_bif(bife_sh, [1 13 17], "figRe_N_sh.svg", "df", "N (max/min)", 0, 20)
fig_bif(bife_sh, [2 14 18], "figRe_R_sh.svg", "df", "Re (max/min)", 0, .8)
fig_bif(bife_sh, [3 15 19], "figRe_C_sh.svg", "df", "C (max/min)", 0, 10)
fig_bif(bife_sh, [4 16 20], "figRe_D_sh.svg", "df", "D (max/min)", 0, 1)

### combining more than 2 or results
fig_bif2(bife_sh, [1 13 17], [3 15 19], [colN colD], "figRe_sh_p1.svg", "", "Max/min",  0, 9.2)
fig_bif2(bife_sh, [2 14 18], [4 16 20], [colRe colD], "figRe_sh_p2.svg", "", "Max/min",  0, .8)
fig_bif2(bife_sh, [1 13 17], [4 16 20], [colN colD], "figRe_sh_p2.svg", "", "",  0, 8.5)
### actual results
fig_bif(bife_sh, [2 14 18], "figRe_R_sh.svg", "", "", 0, .5, colRe, 3)
barplotmean(bife_sh, [1 13 17], "hist_bife_sh_N.svg", 0, 6, colN)
barplotmean(bife_sh, [4 16 20], "hist_bife_sh_D.svg", 0.17, 0.18, colD)

### ii- NDCReRle

param = EcoPar()
param.a42 = 0.15
param.a43 = 0.05

bifle_sh = bif_analysis2(0.0:0.01:1.0, eco_NDCReRle_sh!, fill(0.03, 25), (0.0, 30000.0), 25000:30000, param)

fig_bif(bifle_sh, [1 16 21], "figRle_N_sh.svg", "df", "N (max/min)", 0, 1)
fig_bif(bifle_sh, [2 17 22], "figRle_Re_sh.svg", "df", "Re (max/min)", 0, .8)
fig_bif(bifle_sh, [3 18 23], "figRle_Rle_sh.svg", "df", "Rle (max/min)", 0, 1)
fig_bif(bifle_sh, [4 19 24], "figRle_C_sh.svg", "df", "C (max/min)", 0, 2)
fig_bif(bifle_sh, [5 20 25], "figRle_D_sh.svg", "df", "D (max/min)", 0, 1)

### combining more than 2 or results
fig_bif2(bifle_sh, [1 16 21], [3 15 19], [colN colC], "figRle_sh_p1.svg", "df", "Max/min",  0, 1.45)
fig_bif3(bifle_sh, ([2 14 18], [4 16 20], [3 18 23]), [colRe colD colRle], "figRle_sh_p2.svg", "df", "Max/min",  0, 1.45)

### actual results
fig_bif(bifle_sh, [2 17 22], "figRle_sh_Re.svg", "", "", 0, .6, colRe, 3)
barplotmean(bifle_sh, [1 16 21], "hist_bifle_sh_N.svg", 0, .17, colN)
barplotmean(bifle_sh, [5 20 25], "hist_bifle_sh_D.svg", 0, 1, colD)

max_stars_repo_name: McCannLab/NDID

max_issues_repo_name: McCannLab/NDID

max_forks_repo_name: McCannLab/NDID


text:
using CSV, DataFrames, LaTeXStrings, MVDA, Plots

# Global plot options
upscale = 1.0
font_main = Plots.font("Computer Modern", 18)
font_large = Plots.font("Computer Modern", 16)
font_small = Plots.font("Computer Modern", 12)
default()
default(
    annotationfontfamily="Computer Modern",
    markerstrokewidth=0.5,
    titlefont=font_main,
    guidefont=font_large,
    tickfont=font_small,
    legendfont=font_small,
    markersize=5,
    linewidth=3,
    dpi=200,
    thickness_scaling=1.0,
    top_margin=5Plots.mm,
)
scalefontsizes(upscale)

# Helper function to add a separate x-axis that emphasizes the number of features/support vectors in a model.
function add_model_size_guide(fig, N, is_linear)
    # Compute ticks based on maximum model size N.
    if N > 16
        xticks = collect(round.(Int, N .* range(0, 1, length=11)))
    else
        xticks = collect(0:N)
    end
    sort!(xticks, rev=true)

    xlabel = ifelse(is_linear, "no. features", "no. support vectors")

    # Register figure inside main subplot and append extra x-axis.
    model_size_guide = plot(yticks=nothing, xticks=xticks, xlim=(0,N), xlabel=xlabel, xflip=true)
    full_figure = plot(fig, model_size_guide, layout=@layout [a{1.0h}; b{1e-8h}])

    return full_figure
end

function make_plots(filename)
    # Common plot options.
    w, h = default(:size)
    options = (; left_margin=5Plots.mm, size=(1.5*w, 1.5*h),)

    # Read file containing CV results
    PROGRESS_ID = "make_plot"
    println()
    @info "Reading $(filename)" _id=PROGRESS_ID
    partial_filename = first(splitext(filename))
    cv_results = CSV.read(partial_filename*".dat", DataFrame, header=true)

    example = first(cv_results.example)
    data_subsets = first(cv_results.ntrain), first(cv_results.nvalidate), first(cv_results.ntest)
    p, c = first(cv_results.nfeatures), first(cv_results.nclasses)
    is_linear = first(cv_results.kernel) == "none"
    nvars = is_linear ? p : data_subsets[1]

    titles = ["$(example)\n$(_n) samples / $(p) features / $(c) classes" for _n in data_subsets]
    metrics = (:train, :validation, :test)

    # Summarize CV results over folds + make plot.
    @info "Summarizing over folds" _id=PROGRESS_ID
    cv_paths = MVDA.cv_error(cv_results)
    for (title, metric) in zip(titles, metrics)
        fig = MVDA.plot_cv_paths(cv_paths, metric)
        plot!(fig; title=title, options...)
        fig = add_model_size_guide(fig, nvars, is_linear)
        figname = partial_filename*"-replicates=$(metric).png"
        @info "Saving $(figname)" _id=PROGRESS_ID
        savefig(fig, figname)
    end

    # Construct credible intervals for detailed summary plot.
    @info "Constructing credible intervals" _id=PROGRESS_ID
    cv_intervals = MVDA.credible_intervals(cv_paths)
    for (title, metric) in zip(titles, metrics)
        fig = MVDA.plot_credible_intervals(cv_intervals, metric)
        plot!(fig; title=title, options...)
        fig = add_model_size_guide(fig, nvars, is_linear)
        figname = partial_filename*"-summary=$(metric).png"
        @info "Saving $(figname)" _id=PROGRESS_ID
        savefig(fig, figname)
    end
end

dir = ARGS[1]
@info "Reading from directory: $(dir)"
files = readdir(dir, join=true)
filter!(contains(".dat"), files)
for filename in files
    make_plots(filename)
end

ext: jl


text:
# function warn_type(diag::Diagrams, loopVar, siteVar, evalPropagator, evalNodeFactor = nothing, root = diag.root)
#     @code_warntype evalNaive!(diag, loopVar, siteVar, evalPropagator, evalNodeFactor, root)
# end

function evalNaive!(diag::ExpressionTree, loopVar, siteVar, eval = DiagTree.eval)
    loopPool = diag.loopBasis
    tree = diag.node
    tweight = tree.current

    # calculate new loop
    if isnothing(loopPool) == false
        update(loopPool, loopVar)
    end

    #calculate diagram tree
    for (ni, node) in enumerate(tree.object)
        if isempty(node.children)
            if isnothing(loopPool) == false
                tweight[ni] = eval(node.para, current(loopPool, node.loopidx), node.siteidx, siteVar) * node.factor
            else
                tweight[ni] = eval(node.para, nothing, node.siteidx, siteVar) * node.factor
            end
        else
            if node.operation == MUL
                tweight[ni] = node.factor
                for nidx in node.children
                    tweight[ni] *= tweight[nidx]
                end

            elseif node.operation == ADD
                tweight[ni] = 0.0
                for nidx in node.children
                    tweight[ni] += tweight[nidx]
                end
                tweight[ni] *= node.factor
            else
                error("not implemented!")
            end
        end
    end
end


text:
type Arc
    initNode::Int
    termNode::Int
    capacity::Float64
    freeflowtime::Float64
    trueflow::Float64
    obsflow::Float64
end

Arc(initNode::Int, termNode::Int, capacity::Float64, freeflowtime::Float64) =
    Arc(initNode, termNode, capacity, freeflowtime, 0., 0.)

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
export
    AVRESAMPLE_MAX_CHANNELS,
    AVAudioResampleContext,
    AVMixCoeffType,
    AV_MIX_COEFF_TYPE_Q8,
    AV_MIX_COEFF_TYPE_Q15,
    AV_MIX_COEFF_TYPE_FLT,
    AV_MIX_COEFF_TYPE_NB,
    AVResampleFilterType,
    AV_RESAMPLE_FILTER_TYPE_CUBIC,
    AV_RESAMPLE_FILTER_TYPE_BLACKMAN_NUTTALL,
    AV_RESAMPLE_FILTER_TYPE_KAISER,
    AVResampleDitherMethod,
    AV_RESAMPLE_DITHER_NONE,
    AV_RESAMPLE_DITHER_RECTANGULAR,
    AV_RESAMPLE_DITHER_TRIANGULAR,
    AV_RESAMPLE_DITHER_TRIANGULAR_HP,
    AV_RESAMPLE_DITHER_TRIANGULAR_NS,
    AV_RESAMPLE_DITHER_NB

const AVRESAMPLE_MAX_CHANNELS = 32

typealias AVAudioResampleContext Void

# begin enum AVMixCoeffType
typealias AVMixCoeffType UInt32
const AV_MIX_COEFF_TYPE_Q8 = UInt32(0)
const AV_MIX_COEFF_TYPE_Q15 = UInt32(1)
const AV_MIX_COEFF_TYPE_FLT = UInt32(2)
const AV_MIX_COEFF_TYPE_NB = UInt32(3)
# end enum AVMixCoeffType

# begin enum AVResampleFilterType
typealias AVResampleFilterType UInt32
const AV_RESAMPLE_FILTER_TYPE_CUBIC = UInt32(0)
const AV_RESAMPLE_FILTER_TYPE_BLACKMAN_NUTTALL = UInt32(1)
const AV_RESAMPLE_FILTER_TYPE_KAISER = UInt32(2)
# end enum AVResampleFilterType

# begin enum AVResampleDitherMethod
typealias AVResampleDitherMethod UInt32
const AV_RESAMPLE_DITHER_NONE = UInt32(0)
const AV_RESAMPLE_DITHER_RECTANGULAR = UInt32(1)
const AV_RESAMPLE_DITHER_TRIANGULAR = UInt32(2)
const AV_RESAMPLE_DITHER_TRIANGULAR_HP = UInt32(3)
const AV_RESAMPLE_DITHER_TRIANGULAR_NS = UInt32(4)
const AV_RESAMPLE_DITHER_NB = UInt32(5)
# end enum AVResampleDitherMethod


text:
include("../src/utilities.jl")

@testset "Covariance estimator stats massaging" begin
    wanted = [:white, :white , :bogus , :nw]
    needed = ([:white], [:nw])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = [:hc0, :hc1, :nw]
    needed = ([:hc0, :hc1], [:nw])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = ["White", "HC0" , "nw" , "nothing"]
    needed = ([:white, :hc0], [:nw])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = ["Hc1", "HC0" , "nothing"]
    needed = ([:hc0 , :hc1], [])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = "Nw"
    needed = ([], [:nw])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = :hc0
    needed = ([:hc0], [])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = :hc0
    needed = ([:hc0], [])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = []
    needed = ([], [])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = :none
    needed = ([], [])
    @test needed == get_needed_robust_cov_stats(wanted)

    wanted = :all
    needed = ([:white, :hc0, :hc1, :hc2, :hc3], [:nw])
    @test needed == get_needed_robust_cov_stats(wanted)

end

@testset "model stats massaging" begin
    wanted = ["r2", "rmse"]
    needed = Set([:r2, :rmse, :coefs, :sse, :mse, :sst])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["aic"]
    needed = Set([:coefs, :mse, :sse, :aic])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["Ci"]
    needed = Set([:coefs, :mse, :sse, :ci, :sigma, :stderror, :t_statistic])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["Adjr2"]
    needed = Set([:coefs, :mse, :sse, :r2, :adjr2, :sst])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["stderror"]
    needed = Set([:coefs, :mse, :sse, :stderror, :sigma])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["t_values"]
    needed = Set([:coefs, :mse, :sse, :t_values, :stderror, :sigma])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["p_values"]
    needed = Set([:coefs, :mse, :sse, :p_values, :t_values, :stderror, :sigma])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["ci", "P_values"]
    needed = Set([:coefs, :mse, :sse, :p_values, :t_values, :stderror, :sigma, :ci, :t_statistic])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["none"]
    needed = Set([:coefs, :mse, :sse])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["all"]
    needed = Set([:coefs, :sse, :mse, :sst, :rmse, :aic, :sigma, :t_statistic, :vif, :r2, :adjr2, :stderror, :t_values, :p_values, :ci])
    @test needed == get_needed_model_stats(wanted)
    wanted = [ ]
    needed = Set([:coefs, :mse, :sse])
    @test needed == get_needed_model_stats(wanted)
    wanted = [ "bogus"]
    needed = Set([:coefs, :mse, :sse])
    @test needed == get_needed_model_stats(wanted)
    wanted = ["stderror", "Bogus"]
    needed = Set([:coefs, :mse, :sse, :stderror, :sigma])
    @test needed == get_needed_model_stats(wanted)
    wanted = Set([:stderror, :Bogus])
    needed = Set([:coefs, :mse, :sse, :stderror, :sigma])
    @test needed == get_needed_model_stats(wanted)
    wanted = Set()
    needed = Set([:coefs, :mse, :sse])
    @test needed == get_needed_model_stats(wanted)

end

@testset "prediction stats massaging" begin
    wanted = ["none"]
    target_need = Set([:predicted])
    target_present = Set([:predicted])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["all"]
    target_need = Set([:predicted, :residuals, :leverage, :stdp, :stdi, :stdr, :student, :rstudent, :lcli, :ucli, :lclp, :uclp, :press, :cooksd])
    target_present = Set([:predicted, :residuals, :leverage, :stdp, :stdi, :stdr, :student, :rstudent, :lcli, :ucli, :lclp, :uclp, :press, :cooksd])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["CooksD"]
    target_need = Set([:predicted, :residuals, :leverage, :stdp, :stdr, :student, :cooksd])
    target_present = Set([:predicted, :cooksd])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["CooksD", "lcli"]
    target_need = Set([:predicted, :residuals, :leverage, :stdp, :stdr, :student, :cooksd, :stdi, :lcli])
    target_present = Set([:predicted, :cooksd, :lcli])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["CooksD", "Bogus", "lcli"]
    target_need = Set([:predicted, :residuals, :leverage, :stdp, :stdr, :student, :cooksd, :stdi, :lcli])
    target_present = Set([:predicted, :cooksd, :lcli])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["PRESS"]
    target_need = Set([:predicted, :residuals, :leverage, :press])
    target_present = Set([:predicted, :press])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["lclp"]
    target_need = Set([:predicted, :stdp, :leverage, :lclp])
    target_present = Set([:predicted, :lclp])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["uclp"]
    target_need = Set([:predicted, :stdp, :leverage, :uclp])
    target_present = Set([:predicted, :uclp])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["ucli"]
    target_need = Set([:predicted, :stdi, :leverage, :ucli])
    target_present = Set([:predicted, :ucli])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["ucli", "LCLP"]
    target_need = Set([:predicted, :stdi, :leverage, :ucli, :lclp, :stdp])
    target_present = Set([:predicted, :ucli, :lclp])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = ["student"]
    target_need = Set([:predicted, :residuals, :leverage, :stdr, :student])
    target_present = Set([:predicted, :student])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = [:rstudent]
    target_need = Set([:predicted, :residuals, :leverage, :stdr, :student, :rstudent])
    target_present = Set([:predicted, :rstudent])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = [:stdp]
    target_need = Set([:predicted, :leverage, :stdp])
    target_present = Set([:predicted, :stdp])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = [:stdp, :stdi]
    target_need = Set([:predicted, :leverage, :stdp, :stdi])
    target_present = Set([:predicted, :stdp, :stdi])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = [:stdpp, :stdr]
    target_need = Set([:predicted, :leverage, :stdr])
    target_present = Set([:predicted, :stdr])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = []
    target_need = Set([:predicted])
    target_present = Set([:predicted])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

    wanted = Set()
    target_need = Set([:predicted])
    target_present = Set([:predicted])
    need, present = get_prediction_stats(wanted)
    @test target_need == need
    @test target_present == present

end


text:
import Base: *, eltype, size
import LinearAlgebra: adjoint, mul!, rmul!, lmul!

# The reflector! code in stdlib has no underflow or accuracy protection

# These are translations of xLARFG from LAPACK
# LAPACK Copyright:
# Univ. of Tennessee
# Univ. of California Berkeley
# Univ. of Colorado Denver
# NAG Ltd.
function _reflector!(x::AbstractVector{T}) where {T<:Real}
    require_one_based_indexing(x)
    n = length(x)
    n <= 1 && return zero(T)
    sfmin = 2floatmin(T) / eps(T)
    @inbounds begin
        α = x[1]
        xnorm = _norm2(view(x,2:n))
        if iszero(xnorm)
            return zero(T)
        end
        β = -copysign(hypot(α, xnorm), α)
        kount = 0
        smallβ = abs(β) < sfmin
        if smallβ
            # recompute xnorm and β if needed for accuracy
            rsfmin = one(T) / sfmin
            while smallβ
                kount += 1
                for j in 2:n
                    x[j] *= rsfmin
                end
                β *= rsfmin
                α *= rsfmin
                # CHECKME: is 20 adequate for BigFloat?
                smallβ = (abs(β) < sfmin) && (kount < 20)
            end
            # now β ∈ [sfmin,1]
            xnorm = _norm2(view(x,2:n))
            β = -copysign(hypot(α, xnorm), α)
        end
        τ = (β - α) / β
        t = one(T) / (α - β)
        for j in 2:n
            x[j] *= t
        end
        for j in 1:kount
            β *= sfmin
        end
        x[1] = β
    end
    return τ
end

function _reflector!(x::AbstractVector{T}) where {T<:Complex}
    require_one_based_indexing(x)
    n = length(x)
    # we need to make subdiagonals real so the n=1 case is nontrivial for complex eltype
    n < 1 && return zero(T)
    RT = real(T)
    sfmin = floatmin(RT) / eps(RT)
    @inbounds begin
        α = x[1]
        αr, αi = reim(α)
        xnorm = _norm2(view(x,2:n))
        if iszero(xnorm) && iszero(αi)
            return zero(T)
        end
        β = -copysign(_hypot3(αr, αi, xnorm), αr)
        kount = 0
        smallβ = abs(β) < sfmin
        if smallβ
            # recompute xnorm and β if needed for accuracy
            rsfmin = one(real(T)) / sfmin
            while smallβ
                kount += 1
                for j in 2:n
                    x[j] *= rsfmin
                end
                β *= rsfmin
                αr *= rsfmin
                αi *= rsfmin
                smallβ = (abs(β) < sfmin) && (kount < 20)
            end
            # now β ∈ [sfmin,1]
            xnorm = _norm2(view(x,2:n))
            α = complex(αr, αi)
            β = -copysign(_hypot3(αr, αi, xnorm), αr)
        end
        τ = complex((β - αr) / β, -αi / β)
        t = one(T) / (α - β)
        for j in 2:n
            x[j] *= t
        end
        for j in 1:kount
            β *= sfmin
        end
        x[1] = β
    end
    return τ
end

# copied from Andreas Noack's GenericLinearAlgebra.jl, with trivial mods

"""
a Householder reflection represented as the essential part of the
vector and the normalizing factor
"""
struct Householder{T,S<:StridedVector}
    v::S
    τ::T
end

# warning; size -> length(v) in GLA, but this makes more sense to us:
size(H::Householder) = (length(H.v)+1, length(H.v)+1)
size(H::Householder, i::Integer) = i <= 2 ? length(H.v)+1 : 1

eltype(H::Householder{T})      where T = T

adjoint(H::Householder{T})      where {T} = Adjoint{T,typeof(H)}(H)

function lmul!(H::Householder, A::StridedMatrix)
    m, n = size(A)
    size(H,1) == m || throw(DimensionMismatch(""))
    v = view(H.v, 1:m - 1)
    τ = H.τ
    for j = 1:n
        va = A[1,j]
        Aj = view(A, 2:m, j)
        va += dot(v, Aj)
        va = τ*va
        A[1,j] -= va
        axpy!(-va, v, Aj)
    end
    A
end

function rmul!(A::StridedMatrix, H::Householder)
    m, n = size(A)
    size(H,1) == n || throw(DimensionMismatch(""))
    v = view(H.v, :)
    τ = H.τ
    a1 = view(A, :, 1)
    A1 = view(A, :, 2:n)
    x = A1*v
    axpy!(one(τ), a1, x)
    axpy!(-τ, x, a1)
    rankUpdate!(-τ, x, v, A1)
    A
end

function lmul!(adjH::Adjoint{<:Any,<:Householder}, A::StridedMatrix)
    H = parent(adjH)
    m, n = size(A)
    size(H,1) == m || throw(DimensionMismatch("A: $m,$n H: $(size(H))"))
    v = view(H.v, 1:m - 1)
    τ = H.τ
    for j = 1:n
        va = A[1,j]
        Aj = view(A, 2:m, j)
        va += dot(v, Aj)
        va = τ'va
        A[1,j] -= va
        axpy!(-va, v, Aj)
    end
    A
end

Base.convert(::Type{Matrix}, H::Householder{T}) where {T} = lmul!(H, Matrix{T}(I, size(H, 1), size(H, 1)))
Base.convert(::Type{Matrix{T}}, H::Householder{T}) where {T} = lmul!(H, Matrix{T}(I, size(H, 1), size(H, 1)))


text:
# Translation of sqltypes.h; C typealiases for SQL functions
# http://msdn.microsoft.com/en-us/library/windows/desktop/ms716298(v=vs.85).aspx
# http://msdn.microsoft.com/en-us/library/windows/desktop/aa383751(v=vs.85).aspx
const SQLCHAR =       UInt8
const SQLSMALLINT =   Cshort
const SQLUSMALLINT =  Cushort

const SQLSCHAR =      Cchar
const SQLDATE =       Cuchar
const SQLDECIMAL =    Cuchar
const SQLDOUBLE =     Cdouble
const SQLFLOAT =      Cdouble

const SQLVARCHAR =    Cuchar
const SQLNUMERIC =    Cuchar
const SQLREAL =       Cfloat
const SQLTIME =       Cuchar
const SQLTIMESTAMP =  Cuchar

# ODBC API	64-bit platform	32-bit platform
# SQLINTEGER	32 bits	32 bits
# SQLUINTEGER	32 bits	32 bits
# SQLLEN	64 bits	32 bits
# SQLULEN	64 bits	32 bits
# SQLSETPOSIROW	64 bits	16 bits
# SQL_C_BOOKMARK	64 bits	32 bits
# BOOKMARK	64 bits	32 bits

const SQLINTEGER =  Cint
const SQLUINTEGER = Cuint
const SQLLEN = Int
const SQLULEN = UInt

# if WORD_SIZE == 64
#     const SQLINTEGER =    Cint
#     const SQLUINTEGER =   Cuint
# else
#     const SQLINTEGER =    Clong
#     const SQLUINTEGER =   Culong
# end
#
# const SQLLEN =        SQLINTEGER
# const SQLULEN =       SQLUINTEGER
const SQLSETPOSIROW = SQLUSMALLINT

const SQLROWCOUNT =   SQLULEN
const SQLROWSETSIZE = SQLULEN
const SQLTRANSID =    SQLULEN
const SQLROWOFFSET =  SQLLEN
const SQLPOINTER =    Ptr{Cvoid}
const SQLRETURN =     SQLSMALLINT
const SQLHANDLE =     Ptr{Cvoid}
const SQLHENV =       SQLHANDLE
const SQLHDBC =       SQLHANDLE
const SQLHSTMT =      SQLHANDLE
const SQLHDESC =      SQLHANDLE
const ULONG =         Cuint
const PULONG =        Ptr{ULONG}
const USHORT =        Cushort
const PUSHORT =       Ptr{USHORT}
const UCHAR =         Cuchar
const PUCHAR =        Ptr{Cuchar}
const PSZ =           Ptr{Cchar}
const SCHAR =         Cchar
const SDWORD =        Cint
const SWORD =         Cshort
const UDWORD =        Cuint
const UWORD =         Cushort
const SLONG =         Cint
const SSHORT =        Cshort
const SDOUBLE =       Cdouble
const LDOUBLE =       Cdouble
const SFLOAT =        Cfloat
const PTR =           Ptr{Cvoid}
const HENV =          Ptr{Cvoid}
const HDBC =          Ptr{Cvoid}
const HSTMT =         Ptr{Cvoid}
const RETCODE =       Cshort
const SQLHWND =       Ptr{Cvoid}

#################

# Data Type Mappings
# SQL data types are returned in resultset metadata calls (ODBCMetadata)
# C data types are used in SQLBindCol (ODBCFetch) to allocate column memory; the driver manager converts from the SQL type to this C type in memory
# Julia types indicate how julia should read the returned C data type memory from the previous step

# SQL Data Type                     C Data Type                         Julia Type
# ---------------------------------------------------------------------------------
# SQL_CHAR                          SQL_C_CHAR                          UInt8
# SQL_VARCHAR                       SQL_C_CHAR                          UInt8
# SQL_LONGVARCHAR                   SQL_C_CHAR                          UInt8
# SQL_WCHAR                         SQL_C_WCHAR                         Cwchar_t
# SQL_WVARCHAR                      SQL_C_WCHAR                         Cwchar_t
# SQL_WLONGVARCHAR                  SQL_C_WCHAR                         Cwchar_t
# SQL_DECIMAL                       SQL_C_DOUBLE                        SQLNumeric
# SQL_NUMERIC                       SQL_C_DOUBLE                        SQLNumeric
# SQL_SMALLINT                      SQL_C_SHORT                         Int16
# SQL_INTEGER                       SQL_C_LONG                          Int32
# SQL_REAL                          SQL_C_FLOAT                         Float64
# SQL_FLOAT                         SQL_C_DOUBLE                        Float64
# SQL_DOUBLE                        SQL_C_DOUBLE                        Float64
# SQL_BIT                           SQL_C_BIT                           Int8
# SQL_TINYINT                       SQL_C_TINYINT                       Int8
# SQL_BIGINT                        SQL_C_BIGINT                        Int64
# SQL_BINARY                        SQL_C_BINARY                        UInt8
# SQL_VARBINARY                     SQL_C_BINARY                        UInt8
# SQL_LONGVARBINARY                 SQL_C_BINARY                        UInt8
# SQL_TYPE_DATE                     SQL_C_TYPE_DATE                     SQLDate
# SQL_TYPE_TIME                     SQL_C_TYPE_TIME                     SQLTime
# SQL_TYPE_TIMESTAMP                SQL_C_TYPE_TIMESTAMP                SQLTimestamp
# SQL_INTERVAL_MONTH                SQL_C_INTERVAL_MONTH                UInt8
# SQL_INTERVAL_YEAR                 SQL_C_INTERVAL_YEAR                 UInt8
# SQL_INTERVAL_YEAR_TO_MONTH        SQL_C_INTERVAL_YEAR_TO_MONTH        UInt8
# SQL_INTERVAL_DAY                  SQL_C_INTERVAL_DAY                  UInt8
# SQL_INTERVAL_HOUR                 SQL_C_INTERVAL_HOUR                 UInt8
# SQL_INTERVAL_MINUTE               SQL_C_INTERVAL_MINUTE               UInt8
# SQL_INTERVAL_SECOND               SQL_C_INTERVAL_SECOND               UInt8
# SQL_INTERVAL_DAY_TO_HOUR          SQL_C_INTERVAL_DAY_TO_HOUR          UInt8
# SQL_INTERVAL_DAY_TO_MINUTE        SQL_C_INTERVAL_DAY_TO_MINUTE        UInt8
# SQL_INTERVAL_DAY_TO_SECOND        SQL_C_INTERVAL_DAY_TO_SECOND        UInt8
# SQL_INTERVAL_HOUR_TO_MINUTE       SQL_C_INTERVAL_HOUR_TO_MINUTE       UInt8
# SQL_INTERVAL_HOUR_TO_SECOND       SQL_C_INTERVAL_HOUR_TO_SECOND       UInt8
# SQL_INTERVAL_MINUTE_TO_SECOND     SQL_C_INTERVAL_MINUTE_TO_SECOND     UInt8
# SQL_GUID                          SQL_C_GUID                          UUID

# SQL Data Type Definitions
const SQL_NULL_DATA     = -1
const SQL_NO_TOTAL      = -4
const SQL_CHAR          = Int16(  1) # Character string of fixed string length n.
const SQL_VARCHAR       = Int16( 12) # Variable-length character string with a maximum string length n.
const SQL_LONGVARCHAR   = Int16( -1) # Variable length character data. Maximum length is data source–dependent.
const SQL_WCHAR         = Int16( -8) # Unicode character string of fixed string length n
const SQL_WVARCHAR      = Int16( -9) # Unicode variable-length character string with a maximum string length n
const SQL_WLONGVARCHAR  = Int16(-10) # Unicode variable-length character data. Maximum length is data source–dependent
const SQL_DECIMAL       = Int16(  3)
const SQL_NUMERIC       = Int16(  2)
const SQL_SMALLINT      = Int16(  5) # Exact numeric value with precision 5 and scale 0 (signed: –32,768 <= n <= 32,767, unsigned: 0 <= n <= 65,535)
const SQL_INTEGER       = Int16(  4) # Exact numeric value with precision 10 and scale 0 (signed: –2[31] <= n <= 2[31] – 1, unsigned: 0 <= n <= 2[32] – 1)
const SQL_REAL          = Int16(  7) # Signed, approximate, numeric value with a binary precision 24 (zero or absolute value 10[–38] to 10[38]).
const SQL_FLOAT         = Int16(  6) # Signed, approximate, numeric value with a binary precision of at least p. (The maximum precision is driver-defined.)
const SQL_DOUBLE        = Int16(  8) # Signed, approximate, numeric value with a binary precision 53 (zero or absolute value 10[–308] to 10[308]).
const SQL_BIT           = Int16( -7) # Single bit binary data.
const SQL_TINYINT       = Int16( -6) # Exact numeric value with precision 3 and scale 0 (signed: –128 <= n <= 127, unsigned: 0 <= n <= 255)
const SQL_BIGINT        = Int16( -5) # Exact numeric value with precision 19 (if signed) or 20 (if unsigned) and scale 0 (signed: –2[63] <= n <= 2[63] – 1, unsigned: 0 <= n <= 2[64] – 1)
const SQL_BINARY        = Int16( -2) # Binary data of fixed length n.
const SQL_VARBINARY     = Int16( -3) # Variable length binary data of maximum length n. The maximum is set by the user.
const SQL_LONGVARBINARY = Int16( -4) # Variable length binary data. Maximum length is data source–dependent.
const SQL_TYPE_DATE     = Int16( 91) # Year, month, and day fields, conforming to the rules of the Gregorian calendar.
const SQL_TYPE_TIME     = Int16( 92) # Hour, minute, and second fields, with valid values for hours of 00 to 23,
                                     # valid values for minutes of 00 to 59, and valid values for seconds of 00 to 61. Precision p indicates the seconds precision.
const SQL_TYPE_TIMESTAMP = Int16( 93) # Year, month, day, hour, minute, and second fields, with valid values as defined for the DATE and TIME data types.
# SQL Server specific
const SQL_SS_TIME2       = Int16(-154)
const SQL_SS_TIMESTAMPOFFSET = Int16(-155)

#const SQL_INTERVAL_MONTH            = Int16(102)
#const SQL_INTERVAL_YEAR             = Int16(101)
#const SQL_INTERVAL_YEAR_TO_MONTH    = Int16(107)
#const SQL_INTERVAL_DAY              = Int16(103)
#const SQL_INTERVAL_HOUR             = Int16(104)
#const SQL_INTERVAL_MINUTE           = Int16(105)
#const SQL_INTERVAL_SECOND           = Int16(106)
#const SQL_INTERVAL_DAY_TO_HOUR      = Int16(108)
#const SQL_INTERVAL_DAY_TO_MINUTE    = Int16(109)
#const SQL_INTERVAL_DAY_TO_SECOND    = Int16(110)
#const SQL_INTERVAL_HOUR_TO_MINUTE   = Int16(111)
#const SQL_INTERVAL_HOUR_TO_SECOND   = Int16(112)
#const SQL_INTERVAL_MINUTE_TO_SECOND = Int16(113)
const SQL_GUID                      = Int16(-11) # Fixed length GUID.

# C Data Types
const SQL_C_CHAR      = Int16(  1)
const SQL_C_WCHAR     = Int16( -8)
const SQL_C_DOUBLE    = Int16(  8)
const SQL_C_SHORT     = Int16(  5)
const SQL_C_LONG      = Int16(  4)
const SQL_C_FLOAT     = Int16(  7)
const SQL_C_NUMERIC   = Int16(  2)
const SQL_C_BIT       = Int16( -7)
const SQL_C_TINYINT   = Int16( -6)
const SQL_C_BIGINT    = Int16(-27)
const SQL_C_BINARY    = Int16( -2)
const SQL_C_TYPE_DATE = Int16( 91)
const SQL_C_TYPE_TIME = Int16( 92)
const SQL_C_TYPE_TIMESTAMP = Int16( 93)
const SQL_SIGNED_OFFSET = Int16(-20)
const SQL_UNSIGNED_OFFSET = Int16(-22)
const SQL_C_SLONG = (SQL_C_LONG+SQL_SIGNED_OFFSET)
const SQL_C_SSHORT = (SQL_C_SHORT+SQL_SIGNED_OFFSET)
const SQL_C_STINYINT = (SQL_TINYINT+SQL_SIGNED_OFFSET)
const SQL_C_SBIGINT = (SQL_BIGINT+SQL_SIGNED_OFFSET)
const SQL_C_ULONG = (SQL_C_LONG+SQL_UNSIGNED_OFFSET)
const SQL_C_USHORT = (SQL_C_SHORT+SQL_UNSIGNED_OFFSET)
const SQL_C_UTINYINT = (SQL_TINYINT+SQL_UNSIGNED_OFFSET)
const SQL_C_UBIGINT = (SQL_BIGINT+SQL_UNSIGNED_OFFSET)

#const SQL_C_INTERVAL_MONTH            = Int16(102)
#const SQL_C_INTERVAL_YEAR             = Int16(101)
#const SQL_C_INTERVAL_YEAR_TO_MONTH    = Int16(107)
#const SQL_C_INTERVAL_DAY              = Int16(103)
#const SQL_C_INTERVAL_HOUR             = Int16(104)
#const SQL_C_INTERVAL_MINUTE           = Int16(105)
#const SQL_C_INTERVAL_SECOND           = Int16(106)
#const SQL_C_INTERVAL_DAY_TO_HOUR      = Int16(108)
#const SQL_C_INTERVAL_DAY_TO_MINUTE    = Int16(109)
#const SQL_C_INTERVAL_DAY_TO_SECOND    = Int16(110)
#const SQL_C_INTERVAL_HOUR_TO_MINUTE   = Int16(111)
#const SQL_C_INTERVAL_HOUR_TO_SECOND   = Int16(112)
#const SQL_C_INTERVAL_MINUTE_TO_SECOND = Int16(113)
const SQL_C_GUID                      = Int16(-11)

"Convenience mapping of SQL types to their string representation"
const SQL_TYPES = Dict(
      1 => "SQL_CHAR",
     12 => "SQL_VARCHAR",
     -1 => "SQL_LONGVARCHAR",
     -8 => "SQL_WCHAR",
     -9 => "SQL_WVARCHAR",
    -10 => "SQL_WLONGVARCHAR",
      3 => "SQL_DECIMAL",
      2 => "SQL_NUMERIC",
      5 => "SQL_SMALLINT",
      4 => "SQL_INTEGER",
      7 => "SQL_REAL",
      6 => "SQL_FLOAT",
      8 => "SQL_DOUBLE",
     -7 => "SQL_BIT",
     -6 => "SQL_TINYINT",
     -5 => "SQL_BIGINT",
     -2 => "SQL_BINARY",
     -3 => "SQL_VARBINARY",
     -4 => "SQL_LONGVARBINARY",
     91 => "SQL_TYPE_DATE",
     92 => "SQL_TYPE_TIME",
     93 => "SQL_TYPE_TIMESTAMP",
   -154 => "SQL_SS_TIME2",
   -155 => "SQL_SS_TIMESTAMPOFFSET",
    102 => "SQL_INTERVAL_MONTH",
    101 => "SQL_INTERVAL_YEAR",
    107 => "SQL_INTERVAL_YEAR_TO_MONTH",
    103 => "SQL_INTERVAL_DAY",
    104 => "SQL_INTERVAL_HOUR",
    105 => "SQL_INTERVAL_MINUTE",
    106 => "SQL_INTERVAL_SECOND",
    108 => "SQL_INTERVAL_DAY_TO_HOUR",
    109 => "SQL_INTERVAL_DAY_TO_MINUTE",
    110 => "SQL_INTERVAL_DAY_TO_SECOND",
    111 => "SQL_INTERVAL_HOUR_TO_MINUTE",
    112 => "SQL_INTERVAL_HOUR_TO_SECOND",
    113 => "SQL_INTERVAL_MINUTE_TO_SECOND",
    -11 => "SQL_GUID")

"Convenience mapping of SQL types to their C-type equivalent as a string"
const C_TYPES = Dict(
  SQL_C_CHAR => "SQL_C_CHAR",
  SQL_C_WCHAR => "SQL_C_WCHAR",
  SQL_C_DOUBLE => "SQL_C_DOUBLE",
  SQL_C_SHORT => "SQL_C_SHORT",
  SQL_C_LONG => "SQL_C_LONG",
  SQL_C_FLOAT => "SQL_C_FLOAT",
  SQL_C_NUMERIC => "SQL_C_NUMERIC",
  SQL_C_BIT => "SQL_C_BIT",
  SQL_C_TINYINT => "SQL_C_TINYINT",
  SQL_C_BIGINT => "SQL_C_BIGINT",
  SQL_C_BINARY => "SQL_C_BINARY",
  SQL_C_TYPE_DATE => "SQL_C_TYPE_DATE",
  SQL_C_TYPE_TIME => "SQL_C_TYPE_TIME",
  SQL_C_TYPE_TIMESTAMP => "SQL_C_TYPE_TIMESTAMP",
  SQL_SIGNED_OFFSET => "SQL_SIGNED_OFFSET",
  SQL_UNSIGNED_OFFSET => "SQL_UNSIGNED_OFFSET",
  SQL_C_SLONG => "SQL_C_SLONG",
  SQL_C_SSHORT => "SQL_C_SSHORT",
  SQL_C_STINYINT => "SQL_C_STINYINT",
  SQL_C_ULONG => "SQL_C_ULONG",
  SQL_C_USHORT => "SQL_C_USHORT",
  SQL_C_UTINYINT => "SQL_C_UTINYINT",
  SQL_C_GUID => "SQL_C_GUID",
  SQL_C_SBIGINT => "SQL_C_SBIGINT",
)

# success codes
const SQL_SUCCESS           = Int16(0)
const SQL_SUCCESS_WITH_INFO = Int16(1)
const SQL_STILL_EXECUTING   = Int16(2)

# error codes
const SQL_ERROR             = Int16(-1)
const SQL_INVALID_HANDLE    = Int16(-2)
const SQL_NTS               = -3

# status codes
const SQL_STILL_EXECUTING   = Int16(2)
const SQL_NO_DATA           = Int16(100)

const RETURN_VALUES = Dict(SQL_ERROR   => "SQL_ERROR",
                           SQL_NO_DATA => "SQL_NO_DATA",
                           SQL_INVALID_HANDLE  => "SQL_INVALID_HANDLE",
                           SQL_STILL_EXECUTING => "SQL_STILL_EXECUTING")

const SQL_NO_NULLS = Int16(0)


text:
# Autogenerated wrapper script for Gloo_jll for x86_64-w64-mingw32-cxx03
export libgloo

JLLWrappers.@generate_wrapper_header("Gloo")
JLLWrappers.@declare_library_product(libgloo, "libgloo.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgloo,
        "bin\\libgloo.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
The `Romberg` module provides the function `romberg` that
performs a Romberg integration of discrete 1-dimensional data `y` sampled at
equally spaced points over `x`.

Romberg integration combines the trapezoid method of integration with Richardson
extrapolation. Significantly lower error estimates of definite integrals can be
obtained using Romberg integration over the trapezoid method alone, at only
slightly higher computational cost.
"""
module Romberg

export romberg

import Primes, Richardson, LinearAlgebra

@views function _romberg(Δx, y, endsum, factors, numfactors; kws...)
    b, e = firstindex(y), lastindex(y)
    v = (Δx * (sum(y[b+1:e-1]) + endsum), Δx)
    vals = Vector{typeof(v)}(undef, numfactors+1)
    i = numfactors+1
    vals[i] = v
    step = 1
    for (f,K) in factors
        @inbounds for k = 1:K
            step *= f
            sΔx = step*Δx
            if i == 2 # last iteration (empty sum)
                vals[1] = (sΔx * endsum, sΔx)
            else
                vals[i -= 1] = (sΔx * (sum(y[b+step:step:e-step]) + endsum), sΔx)
            end
        end
    end
    return Richardson.extrapolate!(vals, power=2; kws...)
end

function romberg(Δx::Real, y::AbstractVector; kws...)
    n = length(y)
    endsum = n ≥ 2 ? (first(y)+last(y))/2 : (n == 1 ? zero(first(y))/1 : zero(eltype(y))/1)
    Δxf = float(Δx)
    if n <= 2
        I = endsum * Δxf
        return (I, LinearAlgebra.norm(I))
    end
    m = n - 1
    if ispow2(m)
        # fast path: no need to allocate factors array
        k = 0
        while m > 1
            k += 1
            m >>= 1
        end
        return _romberg(Δxf, y, endsum, (2=>k,), k; kws...)
    else
        factors = Primes.factor(m)
        return _romberg(Δxf, y, endsum, factors, sum(values(factors)); kws...)
    end
end

romberg(x::AbstractRange, y::AbstractVector; kws...) = romberg(step(x), y; kws...)

"""
    romberg(Δx::Real, y::AbstractVector)
    romberg(x::AbstractRange, y::AbstractVector)

Integrate `y` sampled at points with equal spacing `Δx`, or equivalently
a range `x` (`Δx=step(x)`) by Romberg integration, which corresponds
to Richardson extrapolation of a trapezoidal rule to smaller and smaller `Δx`.

The return value is a tuple `(I, E)` of the estimated integral `I` and
an estimated upper bound `E` on the error in `I`.

The algorithm is most effective if `length(x) - 1` has many small factors,
ideally being a power of two, but it can handle any length.  If `length(x) - 1`
is a prime number, it is nearly equivalent to the trapezoidal rule without extrapolation.

# Examples

```jldoctest
julia> x = range(0, stop=π, length=2^8+1);

julia> romberg(x, sin.(x))
(2.0000000000000018, 1.9984014443252818e-15)
```
"""
romberg

end # module


text:
number = 21

include("myLibraries/myLibrary.jl")
include("task_$number.jl")

r = Robot("startPositions/pos_task_$number/untitled.sit", animate = true)

numb!(r)

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs


text:
module Reiwa

import Unicode

export ReiwaString, @reiwa_str

const SQUARE_ERA_NAME_REIWA = '\u32ff'
const COMPAT_STRING_REIWA = "令和"
# const IS_UNICODE_V12_1 = Base.Unicode.isassigned(SQUARE_ERA_NAME_REIWA)

"""
    ReiwaString(str::AbstractString)

Create a new `ReiwaString` object that wraps an `str::AbstractString` and
behaves as ordinal UTF-8 `String`.
Most commonly constructed using the `@reiwa_str` macro.

The character `'\\u32ff'` (`'㋿'`: U+32FF SQUARE ERA NAME REIWA) is stored
and shown **AS IS** (not escaped to the Unicode codepoints form),
and normalized to the compatible String `"令和"` by `Unicode.normalize()`.

# Example

```jldoctest
julia> reiwa"\\u337e\\u337d\\u337c\\u337b\\u32ff"
reiwa"㍾㍽㍼㍻㋿"

julia> using Unicode; Unicode.normalize(reiwa"\\u337e\\u337d\\u337c\\u337b\\u32ff", :NFKC)
"明治大正昭和平成令和"
```
"""
struct ReiwaString <: AbstractString
    str::AbstractString
end

"""
    @reiwa_str

Create a new `ReiwaString` object.

# Example

```jldoctest
julia> reiwa"\\u337e\\u337d\\u337c\\u337b\\u32ff"
reiwa"㍾㍽㍼㍻㋿"
```
"""
macro reiwa_str(src::String)
    ReiwaString(unescape_string(src))
end

# io
function Base.show(io::IO, str::ReiwaString)
    components = split(str.str, SQUARE_ERA_NAME_REIWA)
    print(io, "reiwa\"", escape_string(components[1]))
    for cmp = components[2:end]
        print(io, SQUARE_ERA_NAME_REIWA, escape_string(cmp))
    end
    print(io, "\"")
end

(::Type{String})(str::ReiwaString) = String(str.str)

function tocompatstring(io::IO, str::ReiwaString)
    components = split(str.str, SQUARE_ERA_NAME_REIWA)
    print(io, components[1])
    for cmp = components[2:end]
        print(io, COMPAT_STRING_REIWA, cmp)
    end
end
@inline tocompatstring(str::ReiwaString) = sprint(tocompatstring, str, sizehint=lastindex(str)+1)

# interfaces
Base.ncodeunits(str::ReiwaString) = ncodeunits(str.str)
Base.isvalid(str::ReiwaString, i::Integer) = isvalid(str.str, i)

Base.iterate(str::ReiwaString) = iterate(str.str)
Base.iterate(str::ReiwaString, i::Integer) = iterate(str.str, i)

# Unicode
Unicode.isassigned(c::Char) = c == SQUARE_ERA_NAME_REIWA ? true : Base.Unicode.isassigned(c)
Unicode.isassigned(c::Integer) = c == Int32(SQUARE_ERA_NAME_REIWA) ? true : Base.Unicode.isassigned(c)

function Unicode.normalize(str::ReiwaString, nf::Symbol)
    Unicode.normalize((nf ∈ (:NFKC, :NFKD) ? tocompatstring(str) : str.str), nf)
end
function Unicode.normalize(str::ReiwaString; compat::Bool=false, kwargs...)
    Unicode.normalize((compat ? tocompatstring(str) : str.str); compat=compat, kwargs...)
end

end # module


text:
using Documenter
using Merlin

makedocs(
    modules = [Merlin],
    format = :html,
    sitename = "Merlin.jl",
    pages = [
        "Home" => "index.md",
        "Var" => "var.md",
        "Functions" => "functions.md",
        "Graph" => "graph.md",
        "Initializaters" => "initializers.md",
        "Optimizers" => "optimizers.md",
        "Datasets" => "datasets.md",
        "Save and Load" => "save_load.md",
    ]
)

deploydocs(
    # deps = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
    # deps = Deps.pip("pygments", "mkdocs", "python-markdown-math", "mkdocs-material"),
    repo = "github.com/hshindo/Merlin.jl.git",
    julia  = "0.6",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
export HarmonicAngle

"""
    HarmonicAngle(; th0, cth)

A harmonic bond angle between three atoms.
"""
struct HarmonicAngle{D, K} <: SpecificInteraction
    th0::D
    cth::K
end

HarmonicAngle(; th0, cth) = HarmonicAngle{typeof(th0), typeof(cth)}(th0, cth)

@inline @inbounds function force(a::HarmonicAngle, coords_i, coords_j, coords_k, box_size)
    # In 2D we use then eliminate the cross product
    ba = vector_pad3D(coords_j, coords_i, box_size)
    bc = vector_pad3D(coords_j, coords_k, box_size)
    cross_ba_bc = ba × bc
    pa = normalize(trim3D( ba × cross_ba_bc, box_size))
    pc = normalize(trim3D(-bc × cross_ba_bc, box_size))
    angle_term = -a.cth * (acosbound(dot(ba, bc) / (norm(ba) * norm(bc))) - a.th0)
    fa = (angle_term / norm(ba)) * pa
    fc = (angle_term / norm(bc)) * pc
    fb = -fa - fc
    return SpecificForce3Atoms(fa, fb, fc)
end

@inline @inbounds function potential_energy(a::HarmonicAngle, coords_i, coords_j,
                                            coords_k, box_size)
    θ = bond_angle(coords_i, coords_j, coords_k, box_size)
    return (a.cth / 2) * (θ - a.th0) ^ 2
end


text:
using Playground
using Test
using LinearAlgebra

using Distances
using IterTools
using LightGraphs
using SimpleWeightedGraphs
using NearestNeighbors
using StaticArrays
using Plots
using RecipesBase

n = 1000
# Datasets:
# Two perpendicular lines at 0.5 distance apart.
const two_lines = (vcat([@SVector[0.0, rand()] for _ in 1:n÷2],
                        [@SVector[rand(), 0.5] for _ in 1:n÷2]),
                   Euclidean(),
                   0.4)
# 2-Torus
const torus = let θ = rand(n), φ = rand(n), R = 3, r = 1
    ([@SVector[(R+r*cospi(2θ[i]))cospi(2φ[i]),
               (R+r*cospi(2θ[i]))sinpi(2φ[i]),
               r*sinpi(2θ[i])] for i in 1:n],
     Euclidean(),
     1.0)
end
# 3-sphere
const sphere3 = ([normalize(@SVector randn(4)) for _ in 1:n],
                 Euclidean(),
                 1.0)
# Integer points on a grid, manhattan distance.
const grid2 = ([@SVector[i, j, k] for i in 1:10 for j in 1:10 for k in 1:10],
               Cityblock(),
               5)

const datasets = (two_lines, torus, sphere3, grid2)

include("equilateraliterator.jl")
include("geodesiccomplex.jl")


text:
##
cd(@__DIR__())
using FFTViews
using Base.Iterators: cycle, takewhile, filter, accumulate
mappa = reduce(hcat, Vector{Char}.(eachline("input.txt"))) |> FFTView;
prod(((1,1), (3,1), (5, 1), (7, 1), (1, 2)) .|> CartesianIndex) do Δ
    pos = takewhile(c -> c.I[2] < size(mappa, 2), accumulate(+, cycle([Δ]); init=CartesianIndex(0,0)))
    filter(c -> mappa[c] == '#', pos) |> collect |> length
end

##
sum(enumerate(eachline("input.txt"))) do (i, line)
    map([(1,1), (3,1), (5, 1), (7, 1), (1, 2)]) do (r,d)
        (i-1) % d == 0 || return false
        line[(((i-1) ÷ d) * r) % length(line) + 1] == '#'
    end
end |> prod


text:
@doc raw"""
conductance based channel, simple format

$$i = \bar{g} m^a h^b (V - E)$$

$$i = g m^a h^b (V - E)$$
$$\dot{m} = \frac{m_\infty(V) - m}{\tau_m(V)}$$
$$\dot{h} = \frac{h_\infty(V) - h}{\tau_h(V)}$$

activation and inactivation function would be encoded
as Boltzmann function:

$$m_\infty(V) = \frac{1 - \zeta}{1 + \exp(\frac{V_{\text{half}} - V}{k})} + \zeta$$

the time constant would be encoded as guassian function:

$$\tau(V) = C_{\text{base}} + C_{\text{amp}} \exp(\frac{- (V_{\text{max}} - V) ^ 2}{\sigma ^ 2})$$
"""
mutable struct SimpleIonChannel <: AbstractIonChannel
    name::String
    ion::Symbol
    g::Real

    m::Kinetics
    h::Kinetics
end

function update!(ch::SimpleIonChannel)
    if ~isnothing(ch.m.Vhalf)
        ch.m.infty = (V) -> (1 - ch.m.zeta) / (1 + exp((ch.m.Vhalf - V) / ch.m.k)) + ch.m.zeta
    end

    if ~isnothing(ch.m.Cbase)
        ch.m.tau = (V) -> ch.m.Cbase + ch.m.Camp * exp(- ((ch.m.Vmax - V) / ch.m.sigma)^2)
    end

    if ~isnothing(ch.h.Vhalf)
        ch.h.infty = (V) -> (1 - ch.h.zeta) / (1 + exp(-(ch.h.Vhalf - V) / ch.h.k)) + ch.h.zeta
    end

    if ~isnothing(ch.h.Cbase)
        ch.h.tau = (V) -> ch.h.Cbase + ch.h.Camp * exp(- ((ch.h.Vmax - V) / ch.h.sigma)^2)
    end

    ch
end

function dof(channel::SimpleIonChannel)
    _m = channel.m._type == :evolving ? 1 : 0
    _h = channel.h._type == :evolving ? 1 : 0
    [_m, _h]
end

function step(ch::SimpleIonChannel; V::Real, var::Vector{T}, E::Real) where {T <: Real}
    _var_idx = 1
    i = ch.g * (V - E)

    if ch.m._type == :evolving
        m = var[_var_idx]
        _var_idx += 1

        i = i * m ^ ch.m.n
        dm = (ch.m.infty(V) - m) / ch.m.tau(V)
    elseif ch.m._type == :instantaneous
        i = i * ch.m.infty(V) ^ ch.m.n
        dm = nothing
    else
        dm = nothing
    end

    if ch.h._type == :evolving
        h = var[_var_idx]

        i = i * h ^ ch.h.n
        dh = (ch.h.infty(V) - h) / ch.h.tau(V)
    elseif ch.h._type == :instantaneous
        i = i * ch.h.infty(V) ^ ch.h.n
        dh = nothing
    else
        dh = nothing
    end

    return (i, [dm, dh])
end

function itr_kinetics(ch::SimpleIonChannel)
    [ch.m, ch.h]
end

function current(ch::SimpleIonChannel; V::Vector{T}, var::Array{T, 2}, E::T) where {T <: Real}
    _var_idx = 1
    i = ch.g .* (V .- E)

    if ch.m._type == :evolving
        m = var[_var_idx, :]
        _var_idx += 1

        i = i .* m .^ ch.m.n
    elseif ch.m._type == :instantaneous
        i = i .* ch.m.infty.(V) .^ ch.m.n
    else
        nothing
    end

    if ch.h._type == :evolving
        h = var[_var_idx, :]

        i = i .* h .^ ch.h.n
    elseif ch.h._type == :instantaneous
        i = i .* ch.h.infty.(V) .^ ch.h.n
    else
        nothing
    end

    return i
end


text:
using FileTools
using Base.Test

# hello.txt
@test TXTFile("hello.txt") == f"hello.txt"
@test TXTFile == typeof(f"hello.txt")
@test 2 == length(f"hello.txt")
@test 12 == filesize(f"hello.txt")
@test f"hello.txt"[1] == split("""hello\nworld""",'\n')[1]
@test f"hello.txt"[1:end] == split("""hello\nworld""",'\n')

# int-array.txt
save(f"int-array.txt", "1\n2\n3")
@test ["1","2","3"] == parse(f"int-array.txt")

open(f"int-array.txt") do a
  push!(a, "5")
  push!(a, "6")
end
@test ["1","2","3","5","6"] == parse(f"int-array.txt")

push!(f"int-array.txt", "7")
@test ["1","2","3","5","6","7"] == parse(f"int-array.txt")

# 0.5
# for a ∈ f"hello.txt"
for a in f"hello.txt"
  (word,) = a
  println(word)
  break
end


text:
# load packages
using RPOMDPs, RPOMDPModels, RPOMDPToolbox, RobustValueIteration
using SimpleProbabilitySets, DataFrames, Query, CSV, ProgressMeter
const rpbvi = RobustValueIteration

# setup results structures
# solutions: holds policies and full sim results
# data: holds solution values and copmutation time
factors = DataFrame(ID = Int[], Problem = String[], Short_Name = String[],
    Solution = String[], Reward = String[], Information_Function = String[],
    Uncertainty_Size = Float64[], Dynamics = String[])
probnames = ["Rock Diagnosis"]
shortnames = ["rock"]
soltypes = ["Standard", "Robust"]
rewardtypes = ["Standard", "Information"]
infofuncs = ["Simple", "Approximate Entropy"]
uncsizes = [0.001, 0.1, 0.2, 0.3, 0.4]
dyntypes = ["Nominal", "Ambiguous"]
respsols = ["Policy", "Simulation Values"]
headerfacts = ["ID","Problem", "Short Name", "Solution", "Reward",
    "Information Function", "Uncertainty Size", "Dynamics"]
headersols = vcat(factors, respsols)

ind = 0
for pname in probnames, sol in soltypes, r in rewardtypes, i in infofuncs, u in uncsizes, d in dyntypes
    ind += 1
    sname = shortnames[findin(probnames,[pname])[1]]
    push!(factors, [ind, pname, sname, sol, r, i, u, d])
end

for i in 1:size(factors,1)
    if factors[:Reward][i] == "Standard"
        factors[:Information_Function][i] = "NA"
    end
    if (factors[:Solution][i] == "Standard") && (factors[:Dynamics][i] == "Nominal")
        factors[:Uncertainty_Size][i] = 0.0
    end
end
new_factors = unique(factors[:, 2:8])

# simple info reward function
ralphas1 = [[1.0, -1/3, -1/3, -1/3],
         [-1/3, 1.0, -1/3, -1/3],
         [-1/3, -1/3, 1.0, -1/3],
         [-1/3, -1/3, -1/3, 1.0]]

# complicated info reward function
vhi = -10.0
vlo = 0.1
ralphas2 = [[1.0, vhi, vhi, vhi],
        [vhi, 1.0, vhi, vhi],
        [vhi, vhi, 1.0, vhi],
        [vhi, vhi, vhi, 1.0],
        [vlo, -vlo/3, -vlo/3, -vlo/3],
        [-vlo/3, vlo, -vlo/3, -vlo/3],
        [-vlo/3, -vlo/3, vlo, -vlo/3],
        [-vlo/3, -vlo/3, -vlo/3, vlo]]

# set up problems - WIP
function build(sname, robust, info, err, rewardfunc)
    prob = nothing
    if rewardfunc == "Simple"
        r = ralphas1
    elseif rewardfunc == "Approximate Entropy"
        r = ralphas2
    else
        r = nothing
    end
    if sname == "baby"
        if robust == "Robust"
            prob = (info == "Standard") ? BabyRPOMDP(err) : BabyInfoRPOMDP(err, r)
        else
            prob = (info == "Standard") ? BabyPOMDP() : BabyInfoPOMDP(r)
        end
    elseif sname == "maze1x3"
        if robust == "Robust"
            prob = info == "Standard" ? Maze13RPOMDP(err) : Maze13RIPOMDP(err, r)
        else
            prob = info == "Standard" ? Maze13POMDP() : Maze13IPOMDP(r)
        end
    elseif sname == "maze4x3"
        if robust == "Robust"
            prob = info == "Standard" ? Maze43RPOMDP(err) : Maze43RIPOMDP(err, r)
        else
            prob = info == "Standard" ? Maze43POMDP() : Maze43IPOMDP(r)
        end
    elseif sname == "aircraftid"
        if robust == "Robust"
            prob = info == "Standard" ? AircraftRPOMDP(err) : AircraftRIPOMDP(err, r)
        else
            prob = info == "Standard" ? AircraftPOMDP() : AircraftIPOMDP(r)
        end
    elseif sname == "paint"
        if robust == "Robust"
            prob = info == "Standard" ? PaintRPOMDP(err) :PaintRIPOMDP(err, r)
        else
            prob = info == "Standard" ? PaintPOMDP() : PaintIPOMDP(r)
        end
    elseif sname == "shuttle"
        if robust == "Robust"
            prob = info == "Standard" ? ShuttleRPOMDP(err) : ShuttleRIPOMDP(err, r)
        else
            prob = info == "Standard" ? ShuttlePOMDP() : ShuttleIPOMDP(r)
        end
    elseif sname == "tiger"
        if robust == "Robust"
            prob = info == "Standard" ? TigerRPOMDP(0.95, err) : TigerInfoRPOMDP(0.95, err, r)
        else
            prob = info == "Standard" ? TigerPOMDP(0.95) : TigerInfoPOMDP(0.95, r)
        end
    elseif sname == "simplebaby"
        if robust == "Robust"
            prob = info == "Standard" ? SimpleBaby2RPOMDP(-15.0, -10.0, 0.9, err) : BabyInfoRPOMDP(err, r)
        else
            prob = info == "Standard" ? Baby2POMDP(-15.0, -10.0, 0.9) : BabyInfoPOMDP(r)
        end
    elseif sname == "simpletiger"
        if robust == "Robust"
            prob = info == "Standard" ? SimpleTigerRPOMDP(0.95, err) : TigerInfoRPOMDP(0.95, err, r)
        else
            prob = info == "Standard" ? TigerPOMDP(0.95) : TigerInfoPOMDP(0.95, r)
        end
    elseif sname == "rock"
        if robust == "Robust"
            prob = info == "Standard" ? nothing : RockRIPOMDP(r, 0.95, err)
        else
            prob = info == "Standard" ? nothing : RockIPOMDP(r, 0.95)
        end
    end
    prob
end

##################################################

sname = "rock"
sversion = "9.0"
nreps = 100

# info sname runs
dfexp = @from run in new_factors begin
            @where run.Short_Name == sname && run.Reward == "Information"
            @select run
            @collect DataFrame
end

nrows = size(dfexp,1)
probs = Array{Union{POMDP,IPOMDP,RPOMDP,RIPOMDP}}(nrows)
simprobs = Array{Union{POMDP,IPOMDP,RPOMDP,RIPOMDP}}(nrows)
# problems and sim problems
for i = 1:nrows
    sname = dfexp[:Short_Name][i]
    robust = dfexp[:Solution][i]
    info = dfexp[:Reward][i]
    err = dfexp[:Uncertainty_Size][i]
    rewardfunc = dfexp[:Information_Function][i]
    probs[i] = build(sname, robust, info , err, rewardfunc)
    simrobust = (dfexp[:Dynamics][i] == "Nominal") ? "Standard" : "Robust"
    simprobs[i] = build(sname, simrobust, info, err, rewardfunc)
end

nsteps = 150
bdelt = 0.05
nb2 = 20
bs1 = [vcat([b, 1-b],zeros(2)) for b in 0.0:bdelt:1.0]
bs2 = [vcat(zeros(2), [b, 1-b]) for b in 0.0:bdelt:1.0]
bs3 = [psample(zeros(4), ones(4)) for i = 1:nb2]
bs = vcat(bs1, bs2, bs3)
solver = RPBVISolver(beliefpoints = bs, max_iterations = nsteps)
policies = Array{AlphaVectorPolicy}(nrows)
soldynamics = Array{AlphaVectorPolicy}(nrows) # worst-case dynamics
policies[1] = rpbvi.solve(solver, probs[1], verbose = true)
(dfexp[:Dynamics][1] == "Ambiguous") && (soldynamics[1] = rpbvi.solve(solver, simprobs[1]))

# calculate policies
for i = 2:nrows
    println("\rPolicy $i")
    if dfexp[:Solution][i] == "Standard"
        if dfexp[:Information_Function][i] == dfexp[:Information_Function][i-1]
            policies[i] = policies[i-1]
        else
            policies[i] = rpbvi.solve(solver, probs[i], verbose = true)
        end
    elseif dfexp[:Solution][i] == "Robust"
        if dfexp[i, 5:6] == dfexp[i-1, 5:6]
            policies[i] = policies[i-1]
        else
            policies[i] = rpbvi.solve(solver, probs[i], verbose = true)
        end
    else
        policies[i] = rpbvi.solve(solver, probs[i], verbose = true)
    end
end

function findpolicy(dfexp, row::Int, policies)
    if dfexp[:Dynamics][row] == "Nominal"
        if dfexp[:Information_Function][row] == "Simple"
            return policies[1]
        else
            return policies[7]
        end
    else
        if dfexp[:Information_Function][row] == "Simple"
            if dfexp[:Uncertainty_Size][row] == uncsizes[1]
                return policies[13]
            elseif dfexp[:Uncertainty_Size][row] == uncsizes[2]
                return policies[15]
            elseif dfexp[:Uncertainty_Size][row] == uncsizes[3]
                return policies[17]
            elseif dfexp[:Uncertainty_Size][row] == uncsizes[4]
                return policies[19]
            elseif dfexp[:Uncertainty_Size][row] == uncsizes[5]
                return policies[21]
            else
                return nothing
            end
        else
            if dfexp[:Uncertainty_Size][row] == uncsizes[1]
                return policies[23]
            elseif dfexp[:Uncertainty_Size][row] == uncsizes[2]
                return policies[25]
            elseif dfexp[:Uncertainty_Size][row] == uncsizes[3]
                return policies[27]
            elseif dfexp[:Uncertainty_Size][row] == uncsizes[4]
                return policies[29]
            elseif dfexp[:Uncertainty_Size][row] == uncsizes[5]
                return policies[31]
            else
                return nothing
            end
        end
    end
end

for i = 1:nrows
    soldynamics[i] = findpolicy(dfexp, i, policies)
end

function meanci(data::Vector{Float64})
    n = length(data)
    m = mean(data)
    s = std(data)
    tstar = 1.962
    hw = tstar * s / sqrt(n)
    (m - hw, m + hw)
end

psim = RolloutSimulator(max_steps = nsteps)
simvals = [Vector{Float64}(nreps) for _ in 1:nrows] # simulated values
simps = [Vector{Float64}(nreps) for _ in 1:nrows] # simulated percent correct
ves = Array{Float64}(nrows) # expected values
vms = Array{Float64}(nrows) # mean of sim values
vss = Array{Float64}(nrows) # std dev of sim values
vci = Array{Tuple{Float64,Float64}}(nrows) # 95% ci of mean of sim values
pms = Array{Float64}(nrows) # mean percent correct
pss = Array{Float64}(nrows) # std percent correct
pci = Array{Tuple{Float64,Float64}}(nrows) # 95% ci of mean of sim percent correct

# rollout generates state/obs/rewards based on the sim problem, and beliefs based on the solution policy
srand(923475)
for i in 1:nrows
    println("Run $i")
    println("Solution Type: ", new_factors[:Solution][i], " Dynamics, ", new_factors[:Reward][i], " Reward")
    prob = probs[i]
    simprob = simprobs[i]
    sol = policies[i]
    ves[i] = value(sol, initial_belief(prob))
    bu = updater(sol)
    binit = initial_belief_distribution(prob)
    sinit = rand(psim.rng, initial_state_distribution(simprob))
    println("Simulation Type: ", new_factors[:Dynamics][i], " ", typeof(simprob))
    @showprogress 1 "Simulating value..." for j = 1:nreps
        if typeof(simprob) <: Union{RPOMDP,RIPOMDP}
            sv, sp = simulate_worst(psim, simprob, sol, bu, binit, sinit, soldynamics[i].alphas)
        else
            sv, sp = simulate(psim, simprob, sol, bu, binit, sinit)
        end
        simvals[i][j] = sv
        simps[i][j] = sp
    end
    vms[i] = mean(simvals[i])
    vss[i] = std(simvals[i])
    vci[i] = meanci(simvals[i])
    pms[i] = mean(simps[i])
    pss[i] = std(simps[i])
    pci[i] = meanci(simps[i])
end

ids = collect(1:size(dfexp,1))
dfexp[:ID] = ids
rdata = DataFrame(ID = ids, ExpectedValue = ves, SimMean = vms, SimStd = vss,
    SimCI = vci, CorrectMean = pms, CorrectStd = pss, CorrectCI = pci)
ndec = 3
rdata_round = DataFrame(ID = ids, ExpectedValue = round.(ves, ndec),
    SimMean = round.(vms, ndec),
    SimStd = round.(vss, ndec),
    SimCI = [round.(vci[i], ndec) for i in 1:length(vci)],
    CorrectMean = round.(pms, ndec) .* 100,
    CorrectStd = round.(pss, ndec) .* 100,
    CorrectCI = [round.(pci[i], ndec) .* 100 for i in 1:length(pci)])
dfraw = join(dfexp, rdata_round, on = :ID)
simdata = hcat(ids, hcat(simvals...)') |> DataFrame
@show rdata_round
df = sort(dfraw, [order(:Information_Function), order(:Solution),
    order(:Dynamics), order(:Uncertainty_Size)])

path = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data")
fnresults = string("exp_results_", sname, "_", sversion, ".csv")
fnsim = string("exp_sim_values_", sname, "_", sversion, ".csv")
CSV.write(joinpath(path, fnresults), df)
CSV.write(joinpath(path, fnsim), simdata)

using Plots; gr()
rnum = 1
valfunc1 = [maximum(dot(policies[rnum].alphas[i], [b, 1-b, 0.0, 0.0]) for i = 1:length(policies[rnum].alphas))
        for b = 0.0:0.01:1.0]
valfunc13 = [maximum(dot(policies[13].alphas[i], [b, 1-b, 0.0, 0.0]) for i = 1:length(policies[24].alphas))
        for b = 0.0:0.01:1.0]
valfunc15 = [maximum(dot(policies[15].alphas[i], [b, 1-b, 0.0, 0.0]) for i = 1:length(policies[26].alphas))
        for b = 0.0:0.01:1.0]
valfunc17 = [maximum(dot(policies[17].alphas[i], [b, 1-b, 0.0, 0.0]) for i = 1:length(policies[28].alphas))
        for b = 0.0:0.01:1.0]
valfunc19 = [maximum(dot(policies[19].alphas[i], [b, 1-b, 0.0, 0.0]) for i = 1:length(policies[19].alphas))
        for b = 0.0:0.01:1.0]
valfunc21 = [maximum(dot(policies[21].alphas[i], [b, 1-b, 0.0, 0.0]) for i = 1:length(policies[21].alphas))
        for b = 0.0:0.01:1.0]
x = collect(0:0.01:1.0)
vplot = plot(x, valfunc1, xticks = 0:0.1:1, label = "Nominal",
        # title = "Tiger POMDP Value Function",
        xlab = "Belief, P(State = Bad | Position = 1)",
        # xlab = "Belief, P(State = Hungry)",
        ylab = "Expected Total Discounted Reward",
        legend = :top,
        # legend = :topright,
        line = :dash,
        linealpha = 1,
        linewidth = 2)
plot!(x, valfunc13, linewidth = 2, color = :red, linealpha = 1, label = "Robust: 0.001")
plot!(x, valfunc15, linewidth = 2, color = :red, linealpha = 0.7, label = "Robust: 0.1")
plot!(x, valfunc17, linewidth = 2, color = :red, linealpha = 0.6, label = "Robust: 0.2")
plot!(x, valfunc19, linewidth = 2, color = :red, linealpha = 0.5, label = "Robust: 0.3")
plot!(x, valfunc21, linewidth = 2, color = :red, linealpha = 0.3, label = "Robust: 0.4")
fn = string("value_function_", sname, "_simple_3.pdf")
path = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data\\figures\\",fn)
savefig(vplot, path)

using Plots; gr()
rnum = 7
valfunc7 = [maximum(dot(policies[rnum].alphas[i], [b, 1-b, 0.0, 0.0]) for i = 1:length(policies[rnum].alphas))
        for b = 0.0:0.01:1.0]
valfunc23 = [maximum(dot(policies[23].alphas[i], [b, 1-b, 0.0, 0.0]) for i = 1:length(policies[23].alphas))
        for b = 0.0:0.01:1.0]
valfunc215 = [maximum(dot(policies[25].alphas[i], [b, 1-b, 0.0, 0.0]) for i = 1:length(policies[25].alphas))
        for b = 0.0:0.01:1.0]
valfunc27 = [maximum(dot(policies[27].alphas[i], [b, 1-b, 0.0, 0.0]) for i = 1:length(policies[27].alphas))
        for b = 0.0:0.01:1.0]
valfunc29 = [maximum(dot(policies[29].alphas[i], [b, 1-b, 0.0, 0.0]) for i = 1:length(policies[29].alphas))
        for b = 0.0:0.01:1.0]
valfunc31 = [maximum(dot(policies[31].alphas[i], [b, 1-b, 0.0, 0.0]) for i = 1:length(policies[31].alphas))
        for b = 0.0:0.01:1.0]
x = collect(0:0.01:1.0)
vplot = plot(x, valfunc7, xticks = 0:0.1:1, label = "Nominal",
        # title = "Tiger POMDP Value Function",
        xlab = "Belief, P(State = Bad | Position = 1)",
        # xlab = "Belief, P(State = Hungry)",
        ylab = "Expected Total Discounted Reward",
        legend = :top,
        # legend = :topright,
        line = :dash,
        linealpha = 1,
        linewidth = 2)
plot!(x, valfunc23, linewidth = 2, color = :red, linealpha = 1, label = "Robust: 0.001")
plot!(x, valfunc25, linewidth = 2, color = :red, linealpha = 0.7, label = "Robust: 0.1")
plot!(x, valfunc27, linewidth = 2, color = :red, linealpha = 0.6, label = "Robust: 0.2")
plot!(x, valfunc29, linewidth = 2, color = :red, linealpha = 0.5, label = "Robust: 0.3")
plot!(x, valfunc31, linewidth = 2, color = :red, linealpha = 0.3, label = "Robust: 0.4")
fn = string("value_function_", sname, "_complex_3.pdf")
path = joinpath(homedir(),".julia\\v0.6\\RobustInfoPOMDP\\data\\figures\\",fn)
savefig(vplot, path)

max_stars_repo_name: ajkeith/RobustInfoPOMDP

max_issues_repo_name: ajkeith/RobustInfoPOMDP

max_forks_repo_name: ajkeith/RobustInfoPOMDP


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ECOS_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ECOS")
JLLWrappers.@generate_main_file("ECOS", UUID("c2c64177-6a8e-5dca-99a7-64895ad7445f"))
end  # module ECOS_jll


text:
"""
A BlockOperator has a block matrix structure, where each block is an operator,
and it acts on multiple sets.

A BlockOperator is row-like if it only has one row of blocks. In that case, the
destination of the operator is not necessarily a multidict.

A BlockOperator is column-like if it only has one column of blocks. In that case,
the source set of the operator is not necessarily a multidict.
"""
struct BlockOperator{T} <: DictionaryOperator{T}
    operators   ::  AbstractArray{DictionaryOperator{T}, 2}
    src         ::  Dictionary
    dest        ::  Dictionary

    # scratch_src and scratch_dest hold scratch memory for each subset of the
    # source and destination sets, for allocation-free implementation of the
    # action of the operator
    scratch_src
    scratch_dest

    function BlockOperator{T}(operators, src, dest) where {T}
        # Avoid 1x1 block operators
        @assert size(operators,1) + size(operators,2) > 2
        scratch_src = zeros(src)
        scratch_dest = zeros(dest)
        new(operators, src, dest, scratch_src, scratch_dest)
    end
end

function BlockOperator(operators::AbstractArray{OP,2},
            op_src = multidict(map(src, operators[1,:])...),
            op_dest = multidict(map(dest, operators[:,1])...)) where {OP <: DictionaryOperator}
    T = promote_type(operatoreltype(op_src, op_dest), map(eltype, operators)...)
    BlockOperator{T}(operators, op_src, op_dest)
end

function ArrayOperator(A::BlockArray{T}, src::Dictionary, dest::Dictionary) where {T}
    components_src = iscomposite(src) ? components(src) : (src,)
    components_dest = iscomposite(dest) ? components(dest) : (dest,)
    BlockOperator{T}([ArrayOperator(view(A, Block(i, j)), srcj, desti) for (i,desti) in enumerate(components_dest), (j,srcj) in enumerate(components_src)], src, dest)
end

# For legacy reasons
function block_row_operator(op1::DictionaryOperator, op2::DictionaryOperator, dicts::Dictionary...)
    T = promote_type(eltype(op1), eltype(op2), operatoreltype(dicts...))
    operators = Array{DictionaryOperator{T}}(1, 2)
    operators[1] = op1
    operators[2] = op2
    BlockOperator(operators, dicts...)
end

block_row_operator(ops::DictionaryOperator...) = block_row_operator(collect(ops))
function block_row_operator(ops::AbstractArray{<:DictionaryOperator,1}, dicts...)
    T = promote_type(map(eltype, ops)...)
    operators = Array{DictionaryOperator{T}}(undef, 1, length(ops))
    operators[:] = ops
    BlockOperator(operators, dicts...)
end

block_column_operator(ops::DictionaryOperator...) = block_column_operator(collect(ops))
function block_column_operator(ops::AbstractArray{<:DictionaryOperator,1}, dicts...)
    T = promote_type(map(eltype, ops)...)
    operators = Array{DictionaryOperator{T}}(undef, length(ops), 1)
    operators[:] = ops
    BlockOperator(operators, dicts...)
end

unsafe_wrap_operator(src, dest, op::BlockOperator) = similar_operator(op, src, dest)

similar_operator(op::BlockOperator{T}, src, dest) where {T} =
    BlockOperator{T}(op.operators, src, dest)

# Return a block operator the size of the given operator, but filled with zero
# operators.
function zeros(op::BlockOperator)
    operators = zeros(DictionaryOperator{eltype(op)}, composite_size(op))
    for i in 1:composite_size(op, 1)
        for j in 1:composite_size(op, 2)
            operators[i,j] = ZeroOperator(src(op, j), dest(op, i))
        end
    end
    BlockOperator(operators)
end

# Return the source of the i-th column
src(op::BlockOperator, j) = j==1 && iscolumnlike(op) ? src(op) : component(src(op), j)

# Return the destination of the i-th row
dest(op::BlockOperator, i) = i==1 && isrowlike(op) ? dest(op) : component(dest(op), i)

component(op::BlockOperator, i::Int, j::Int) = op.operators[i,j]
components(op::BlockOperator) = op.operators

composite_size(op::BlockOperator) = size(op.operators)

composite_size(op::BlockOperator, dim) = size(op.operators, dim)

isrowlike(op::BlockOperator) = size(op.operators,1) == 1

iscolumnlike(op::BlockOperator) = size(op.operators,2) == 1

function apply!(op::BlockOperator, coef_dest, coef_src)
    if isrowlike(op)
        apply_rowoperator!(op, coef_dest, coef_src, op.scratch_src, op.scratch_dest)
    elseif iscolumnlike(op)
        apply_columnoperator!(op, coef_dest, coef_src, op.scratch_src, op.scratch_dest)
    else
        apply_block_operator!(op, coef_dest, coef_src, op.scratch_src, op.scratch_dest)
    end
    coef_dest
end

function apply_block_operator!(op::BlockOperator, coef_dest::AbstractVector, coef_src::AbstractVector, scratch_src, scratch_dest)
    delinearize_coefficients!(scratch_src, coef_src)
    apply_block_operator!(op, scratch_dest, scratch_src, scratch_src, scratch_dest)
    linearize_coefficients!(coef_dest, scratch_dest)
end

function apply_block_operator!(op::BlockOperator, coef_dest::AbstractVector, coef_src::BlockVector, scratch_src, scratch_dest)
    apply_block_operator!(op, scratch_dest, coef_src, scratch_src, scratch_dest)
    linearize_coefficients!(coef_dest, scratch_dest)
end

function apply_block_operator!(op::BlockOperator, coef_dest::BlockVector, coef_src::AbstractVector, scratch_src, scratch_dest)
    delinearize_coefficients!(scratch_src, coef_src)
    apply_block_operator!(op, coef_dest, scratch_src, scratch_src, scratch_dest)
end

function apply_block_operator!(op::BlockOperator, coef_dest::BlockVector, coef_src::BlockVector, scratch_src, scratch_dest)
    for m in 1:ncomponents(coef_dest)
        fill!(component(coef_dest, m), 0)
        for n in 1:ncomponents(coef_src)
            apply_block_element!(component(op, m, n), component(coef_dest, m),
                component(coef_src,n), component(op.scratch_dest, m))
        end
    end
    coef_dest
end

function apply_block_element!(op, coef_dest, coef_src, scratch)
    apply!(op, scratch, coef_src)
    for i in eachindex(coef_dest)
        coef_dest[i] += scratch[i]
    end
end

function apply_rowoperator!(op::BlockOperator, coef_dest, coef_src::BlockVector, scratch_src, scratch_dest)
    fill!(coef_dest, 0)
    for n in 1:ncomponents(coef_src)
        apply!(op.operators[1,n], scratch_dest, component(coef_src,n))
        for i in eachindex(coef_dest)
            coef_dest[i] += scratch_dest[i]
        end
    end
end

function apply_rowoperator!(op::BlockOperator, coef_dest, coef_src::AbstractVector, scratch_src, scratch_dest)
    delinearize_coefficients!(scratch_src, coef_src)
    apply_rowoperator!(op, coef_dest, scratch_src, scratch_src, scratch_dest)
end

function apply_columnoperator!(op::BlockOperator, coef_dest::BlockVector, coef_src, scratch_src, scratch_dest)
    for m in 1:ncomponents(coef_dest)
        apply!(op.operators[m,1], component(coef_dest, m), coef_src)
    end
end

function apply_columnoperator!(op::BlockOperator, coef_dest::AbstractVector, coef_src, scratch_src, scratch_dest)
    apply_columnoperator!(op, scratch_dest, coef_src, scratch_src, scratch_dest)
    linearize_coefficients!(coef_dest, scratch_dest)
end

hcat(ops::DictionaryOperator...) = block_row_operator(ops...)
vcat(ops::DictionaryOperator...) = block_column_operator(ops...)

function hvcat(rows::Tuple{Vararg{Int}}, ops::DictionaryOperator...)
    T = promote_type(map(eltype, ops)...)
    operators = Array{DictionaryOperator{T}}(undef, floor(Int,length(ops)/rows[1]), rows[1])
    for i in 1:length(operators)
        operators[i] = ops[i]
    end
    BlockOperator(operators)
end

hcat(op1::BlockOperator, op2::BlockOperator) = BlockOperator(hcat(op1.operators, op2.operators))
vcat(op1::BlockOperator, op2::BlockOperator) = BlockOperator(vcat(op1.operators, op2.operators))

adjoint(op::BlockOperator) = BlockOperator(Array{DictionaryOperator}(adjoint(op.operators)), dest(op), src(op))
conj(op::BlockOperator) = BlockOperator(Array{DictionaryOperator}(conj(op.operators)), src(op), dest(op))

############################
# Block diagonal operators
############################

"""
A BlockDiagonalOperator has a block matrix structure like a BlockOperator, but
with only blocks on the diagonal.
"""
struct BlockDiagonalOperator{T} <: DictionaryOperator{T}
    operators   ::  AbstractArray{DictionaryOperator{T}, 1}
    src         ::  Dictionary
    dest        ::  Dictionary
end

function BlockDiagonalOperator(operators::AbstractArray{O,1}, src, dest) where {O<:DictionaryOperator}
    T = promote_type(operatoreltype(src,dest), map(eltype, operators)...)
    BlockDiagonalOperator{T}(operators, src, dest)
end

function BlockDiagonalOperator(operators::AbstractArray{O,1}) where {O<:DictionaryOperator}
    BlockDiagonalOperator(operators, multidict(map(src, operators)), multidict(map(dest, operators)))
end

operators(op::BlockDiagonalOperator) = op.operators
components(op::BlockDiagonalOperator) = op.operators

function block_operator(op::BlockDiagonalOperator)
    ops = Array(DictionaryOperator{eltype(op)}, ncomponents(op), ncomponents(op))
    n = length(op.operators)
    for i in 1:n,j in 1:n
        ops[i,j] = component(op, i, j)
    end
    BlockOperator(ops)
end

composite_size(op::BlockDiagonalOperator) = (l = length(op.operators); (l,l))

block_diagonal_operator(op1::DictionaryOperator, op2::DictionaryOperator) =
    BlockDiagonalOperator([op1,op2])
block_diagonal_operator(op1::BlockDiagonalOperator, op2::BlockDiagonalOperator) =
    BlockDiagonalOperator(vcat(components(op1),components(op2)))
block_diagonal_operator(op1::BlockDiagonalOperator, op2::DictionaryOperator) =
    BlockDiagonalOperator(vcat(components(op1), op2))
block_diagonal_operator(op1::DictionaryOperator, op2::BlockDiagonalOperator) =
    BlockDiagonalOperator(vcat(op1,components(op2)))

function block_diagonal_operator(op1::DictionaryOperator, op2::BlockOperator)
    ops = Array(DictionaryOperator{eltype(op1)}, 1+composite_size(op2,1), 1+composite_size(op2,2))
    ops[1,1] = op1
    for i in 1:composite_size(op2,1)
        ops[1+i,1] = ZeroOperator(src(op1), dest(op2,i))
    end
    for j in 1:composite_size(op2,2)
        ops[1,1+j] = ZeroOperator(src(op2,j), dest(op1))
    end
    ops[2:end,2:end] = op2.operators
    BlockOperator(ops)
end

function block_diagonal_operator(ops::AbstractArray{DictionaryOperator{T}}) where {T}
    @assert length(ops) > 1
    if length(ops) > 2
        block_diagonal_operator(ops[1], ops[2], ops[2:end])
    else
        block_diagonal_operator(ops[1], ops[2])
    end
end

⊕(op1::DictionaryOperator, op2::DictionaryOperator) = block_diagonal_operator(op1, op2)

function component(op::BlockDiagonalOperator{T}, i::Int, j::Int) where {T}
    if i == j
        op.operators[i]
    else
        ZeroOperator{T}(component(op.src, j), component(op.dest, i))
    end
end

apply!(op::BlockDiagonalOperator, coef_dest::BlockVector, coef_src::Array{T,1}) where {T} =
    apply!(op, coef_dest, delinearize_coefficients(coef_dest, coef_src))

function apply!(op::BlockDiagonalOperator, coef_dest::BlockVector, coef_src::BlockVector)
    for i in 1:ncomponents(coef_src)
        apply!(op.operators[i], component(coef_dest, i), component(coef_src, i))
    end
    coef_dest
end

adjoint(op::BlockDiagonalOperator) = BlockDiagonalOperator(map(adjoint, operators(op)), dest(op), src(op))
conj(op::BlockDiagonalOperator) = BlockDiagonalOperator(map(conj, operators(op)), src(op), dest(op))

inv(op::BlockDiagonalOperator) = BlockDiagonalOperator(map(inv, operators(op)), dest(op), src(op))
# inv(op::BlockDiagonalOperator) = BlockDiagonalOperator(DictionaryOperator{eltype(op)}[inv(o) for o in BasisFunctions.operators(op)])

function stencilarray(op::BlockOperator)
    A = Any[]
    push!(A,"[")
    for i=1:size(components(op),1)
        i!=1 && push!(A,";\t")
        push!(A,component(op,i,1))
        for j=2:size(components(op),2)
            push!(A,", \t")
            push!(A,component(op,i,j))
        end
    end
    push!(A,"]")
    A
end

function stencilarray(op::BlockDiagonalOperator)
    A = Any[]
    push!(A,"[")
    for i=1:composite_size(op)[1]
        i!=1 && push!(A,";\t 0")
        i==1 && push!(A,component(op,1,1))
        for j=2:composite_size(op)[2]
            push!(A,", \t")
            i==j && push!(A,component(op,i,j))
            i!=j && push!(A,"0")
        end
    end
    push!(A,"]")
    A
end

Matrix(op::Union{BlockOperator,BlockDiagonalOperator}) =
    Matrix(matrix(op))

function matrix(op::Union{BlockOperator,BlockDiagonalOperator})
    cs = composite_size(op)
    dest_cs = cs[1] == 1 ? [length(dest(op))] :  [map(length, components(dest(op)))...]
    src_cs = cs[2] == 1 ? [length(src(op))] :  [map(length, components(src(op)))...]
    a = BlockArray{eltype(op)}(undef_blocks, dest_cs, src_cs)
    matrix!(op, a)
    a
end

function matrix!(op::Union{BlockOperator,BlockDiagonalOperator}, a)
    for i in 1:composite_size(op)[1]
        for j in 1:composite_size(op)[2]
            a[Block(i,j)] = Matrix(component(op, i, j))
        end
    end
end


text:
export CircularTurnBuffer, circular_RTSA_buffer, capacity, isfull, circular_PRTSA_buffer

import StatsBase: sample
using .Utils: SumTree

"""
    CircularTurnBuffer{names,types,Tbs}

Contain a collection of buffers (mainly are [`CircularArrayBuffer`](@ref), but not restriced to it) to represent the interractions between agents and environments. This struct itself is very simple, some commonly used buffers are provided by:

- [`circular_PRTSA_buffer`](@ref)
- [`circular_PRTSA_buffer`](@ref)

# Fields

- `buffers`: a tuple of inner buffers

"""
struct CircularTurnBuffer{names,types,Tbs} <: AbstractTurnBuffer{names,types}
    buffers::Tbs
end

"""
    sample(b::CircularTurnBuffer, batch_size, n_step, stack_size)

!!! note
    `stack_size` can be an `Int`, and then the size of state is expanded in the last dimension.
    For example if the original state is an image of size `(84, 84)`, then the size of new state is `(84, 84, stack_size)`.
"""
function sample(b::CircularTurnBuffer, batch_size, n_step, stack_size)
    inds = sample_indices(b, batch_size, n_step, stack_size)
    inds, consecutive_view(b, inds, n_step, stack_size)
end

sample_indices(b::CircularTurnBuffer, batch_size::Int, n_step::Int, stack_size::Nothing) = sample_indices(b, batch_size, n_step, 1)

#####
# RTSA
#####

"""
    circular_RTSA_buffer(;kwargs...) -> CircularTurnBuffer

A helper function to help generate a [`CircularTurnBuffer`](@ref) with **RTSA** (**R**eward, **T**erminal, **S**tate, **A**ction) fields as buffers.

# Keywords

## Necessary

- `capacity::Int`: the maximum length of the buffer

## Optional

- `state_eltype::Type=Int`: the type of the state field in an [`Observation`](@ref), `Int` by default
- `state_size::NTuple{N, Int}=()`: the size of the state field in an [`Observation`](@ref), the `N` must match `ndims(state_eltype)`. Since the default `state_eltype` is `Int`, it is an empty tuple here by default.
- `action_eltype::Type=Int`: similar to `state_eltype`
- `action_size::NTuple{N, Int}=()`: similar to `state_size`
- `reward_eltype::Type=Float32`: similar to `state_eltype`
- `reward_size::NTuple{N, Int}=()`: similar to `state_size`
- `terminal_eltype::Type=Bool`: similar to `state_eltype`
- `terminal_size::NTuple{N, Int}=()`: similar to `state_size`

The following picture will help you understand how the data are organized.

![](../assets/img/circular_RTSA_buffer.png)

# Examples

```@julia-repl
julia> using ReinforcementLearning

julia> b = circular_RTSA_buffer(;capacity=2, state_eltype=Array{Float32, 2}, state_size=(2, 2))
0-element CircularTurnBuffer{(:reward, :terminal, :state, :action),Tuple{Float32,Bool,Array{Float32,2},Int64},NamedTuple{(:reward, :terminal, :state, :action),Tuple{CircularArrayBuffer{Float32,Float32,1},CircularArrayBuffer{Bool,Bool,1},CircularArrayBuffer{Array{Float32,2},Float32,3},CircularArrayBuffer{Int64,Int64,1}}}}

julia> push!(b; reward = 0.0, terminal = true, state = Float32[0 0; 0 0], action = 0)

julia> length(b)
0

julia> b.buffers.reward
1-element CircularArrayBuffer{Float32,Float32,1}:
 0.0

julia> b.buffers.terminal
1-element CircularArrayBuffer{Bool,Bool,1}:
 1

julia> b.buffers.state
2×2×1 CircularArrayBuffer{Array{Float32,2},Float32,3}:
[:, :, 1] =
 0.0  0.0
 0.0  0.0

julia> b.buffers.action
1-element CircularArrayBuffer{Int64,Int64,1}:
 0

julia> push!(b; reward = 1.0, terminal = false, state = Float32[1 1; 1 1], action = 1)

julia> b
1-element CircularTurnBuffer{(:reward, :terminal, :state, :action),Tuple{Float32,Bool,Array{Float32,2},Int64},NamedTuple{(:reward, :terminal, :state, :action),Tuple{CircularArrayBuffer{Float32,Float32,1},CircularArrayBuffer{Bool,Bool,1},CircularArrayBuffer{Array{Float32,2},Float32,3},CircularArrayBuffer{Int64,Int64,1}}}}:
 (state = Float32[0.0 0.0; 0.0 0.0], action = 0, reward = 1.0f0, terminal = false, next_state = Float32[1.0 1.0; 1.0 1.0], next_action = 1)

julia> push!(b; reward = 2.0, terminal = true, state = Float32[2 2; 2 2], action = 2)

julia> b
2-element CircularTurnBuffer{(:reward, :terminal, :state, :action),Tuple{Float32,Bool,Array{Float32,2},Int64},NamedTuple{(:reward, :terminal, :state, :action),Tuple{CircularArrayBuffer{Float32,Float32,1},CircularArrayBuffer{Bool,Bool,1},CircularArrayBuffer{Array{Float32,2},Float32,3},CircularArrayBuffer{Int64,Int64,1}}}}:
 (state = Float32[0.0 0.0; 0.0 0.0], action = 0, reward = 1.0f0, terminal = false, next_state = Float32[1.0 1.0; 1.0 1.0], next_action = 1)
 (state = Float32[1.0 1.0; 1.0 1.0], action = 1, reward = 2.0f0, terminal = true, next_state = Float32[2.0 2.0; 2.0 2.0], next_action = 2)

julia> push!(b; reward = 3.0, terminal = false, state = Float32[3 3; 3 3], action = 3)

julia> b
2-element CircularTurnBuffer{(:reward, :terminal, :state, :action),Tuple{Float32,Bool,Array{Float32,2},Int64},NamedTuple{(:reward, :terminal, :state, :action),Tuple{CircularArrayBuffer{Float32,Float32,1},CircularArrayBuffer{Bool,Bool,1},CircularArrayBuffer{Array{Float32,2},Float32,3},CircularArrayBuffer{Int64,Int64,1}}}}:
 (state = Float32[1.0 1.0; 1.0 1.0], action = 1, reward = 2.0f0, terminal = true, next_state = Float32[2.0 2.0; 2.0 2.0], next_action = 2)
 (state = Float32[2.0 2.0; 2.0 2.0], action = 2, reward = 3.0f0, terminal = false, next_state = Float32[3.0 3.0; 3.0 3.0], next_action = 3)

julia> b.buffers.state
2×2×3 CircularArrayBuffer{Array{Float32,2},Float32,3}:
[:, :, 1] =
 1.0  1.0
 1.0  1.0

[:, :, 2] =
 2.0  2.0
 2.0  2.0

[:, :, 3] =
 3.0  3.0
 3.0  3.0

julia> b.buffers.reward
3-element CircularArrayBuffer{Float32,Float32,1}:
 1.0
 2.0
 3.0

julia> length(b)
2
```
"""
function circular_RTSA_buffer(
    ;
    capacity,
    state_eltype = Int,
    state_size = (),
    action_eltype = Int,
    action_size = (),
    reward_eltype = Float32,
    reward_size = (),
    terminal_eltype = Bool,
    terminal_size = (),
)
    capacity += 1  # we need to store extra dummy (reward, terminal)
    buffers = (
        reward = CircularArrayBuffer{reward_eltype}(reward_size..., capacity),
        terminal = CircularArrayBuffer{terminal_eltype}(terminal_size..., capacity),
        state = CircularArrayBuffer{state_eltype}(state_size..., capacity),
        action = CircularArrayBuffer{action_eltype}(action_size..., capacity),
    )
    CircularTurnBuffer{
        RTSA,
        Tuple{reward_eltype,terminal_eltype,state_eltype,action_eltype},
        typeof(buffers),
    }(buffers)
end

sample_indices(b::CircularTurnBuffer{RTSA}, batch_size::Int, n_step::Int, stack_size::Int) =
    rand(stack_size:length(b)-n_step, batch_size)

#####
# PRTSA
#####

"""
    circular_PRTSA_buffer(;kwargs...) -> CircularTurnBuffer

The only difference compared to [`circular_RTSA_buffer`](@ref) is that a new field named `priority` is added. Notice that the struct of `priority` is not a [`CircularArrayBuffer`](@ref) but a [`SumTree`](@ref).

# Keywords

## Necessary

- `capacity::Int`: the maximum length of the buffer

## Optional

- `state_eltype::Type=Int`: the type of the state field in an [`Observation`](@ref), `Int` by default
- `state_size::NTuple{N, Int}=()`: the size of the state field in an [`Observation`](@ref), the `N` must match `ndims(state_eltype)`. Since the default `state_eltype` is `Int`, it is an empty tuple here by default.
- `action_eltype::Type=Int`: similar to `state_eltype`
- `action_size::NTuple{N, Int}=()`: similar to `state_size`
- `reward_eltype::Type=Float32`: similar to `state_eltype`
- `reward_size::NTuple{N, Int}=()`: similar to `state_size`
- `terminal_eltype::Type=Bool`: similar to `state_eltype`
- `terminal_size::NTuple{N, Int}=()`: similar to `state_size`
- `priority_eltype::Type=Float64`: the type of the `priority` field

"""
function circular_PRTSA_buffer(
    ;
    capacity,
    state_eltype = Int,
    state_size = (),
    action_eltype = Int,
    action_size = (),
    reward_eltype = Float64,
    reward_size = (),
    terminal_eltype = Bool,
    terminal_size = (),
    priority_eltype = Float64,
)
    capacity += 1  # we need to store extra dummy (reward, terminal)
    buffers = (
        priority = SumTree(priority_eltype, capacity),
        reward = CircularArrayBuffer{reward_eltype}(reward_size..., capacity),
        terminal = CircularArrayBuffer{terminal_eltype}(terminal_size..., capacity),
        state = CircularArrayBuffer{state_eltype}(state_size..., capacity),
        action = CircularArrayBuffer{action_eltype}(action_size..., capacity),
    )
    CircularTurnBuffer{
        PRTSA,
        Tuple{priority_eltype,reward_eltype,terminal_eltype,state_eltype,action_eltype},
        typeof(buffers),
    }(buffers)
end

function sample_indices(b::CircularTurnBuffer{PRTSA}, batch_size::Int, n_step::Int, stack_size::Int)
    inds = Vector{Int}(undef, batch_size)
    for i = 1:length(inds)
        ind, p = sample(priority(b))
        while ind <= stack_size || ind > length(b) - n_step + 1
            ind, p = sample(priority(b))
        end
        inds[i] = ind - 1  # !!! left shift by 1 because we are padding for priority
    end
    inds
end


text:
function to_graphviz(d::DAG, file::AbstractString)
  isfile(file) && rm(file)
  io = open(file, "w")
  gs = topological_sort(d.a)
  write(io, "digraph $(d.name) {\n")
  vars = names(gs, 1)
  for var in vars
    for (ind, entry) in enumerate(gs[:, var])
      if entry == 1
        write(io, "  $(var) -> $(vars[ind]);\n")
      end
    end
  end
  write(io, "}\n")
  close(io)
end

export
  to_graphviz


text:
module FlicsFlow

# Write your package code here.
export readopenfoammesh
include("Mesh/MeshReader.jl")

end