text:
"""
Placeholder for a short summary about NumericalDebugging.
"""
module NumericalDebugging

using Requires: @require

export @check_if, nonfinite

"""
Whether to throw an error if a check fails (default: `false`).

Not exported, but part of the API.

# Example

```julia
julia> using NumericalDebugging

julia> NumericalDebugging.ERROR[] = true
true

julia> f() = @check_if nonfinite a = Inf
f (generic function with 1 method)

julia> f()
┌ Warning: nonfinite: Inf
│   a = Inf
└ @ Main REPL[18]:1
ERROR: nonfinite variables
```
"""
const ERROR = Ref(false)

####
#### macros
####

function _check_if(_module, _file, _line, f, label, kwargs)
    flagged_vars = collect(k for (k, v) in pairs(kwargs) if f(v))
    if !isempty(flagged_vars)
        @warn("$(label): $(flagged_vars...)",
              _module = _module, _file = _file, _line = _line,
              kwargs...)
        if ERROR[]
            error("$(label) variables")
        end
    end
    nothing
end

"""
    @check_if f arguments...

Check each argument, which is either a `variable` or a `symbol = value` pair, with `f`.

When `f` returns true for *any* argument, emit a warning with *all* values, also throw an
error when [`ERROR`](@ref) is set.

```jldoctest; filter = r"REPL.*"
julia> let a = Inf
           @check_if nonfinite a b = 9
       end
┌ Warning: nonfinite: a
│   a = Inf
│   b = 9
└ @ Main REPL[14]:2
```
"""
macro check_if(f, args...)
    file = String(__source__.file)
    line = __source__.line
    quote
        _check_if($__module__, $file, $line,
                  $f, $(string(f)), (; $(map(esc, args)...)))
    end
end

####
#### diagnostic functions
####

"""
    nonfinite(x)

Return true iff `x` is not finite, or has non-finite elements.
"""
nonfinite(x::Real) = !isfinite(x)

function nonfinite(x)
    if first(x) ≡ x             # catch Char, etc
        if x isa Real
            throw(DomainError())
        else
            false
        end
    else
        any(nonfinite, x)
    end
end

function __init__()
    @require ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" begin
        import .ForwardDiff
        function NumericalDebugging.nonfinite(x::ForwardDiff.Dual)
            nonfinite(ForwardDiff.value(x)) || nonfinite(ForwardDiff.partials(x))
        end
    end
end

end # module


text:
# BASED ON TIKZPICTURES.JL TESTS

using TikzCDs
using Test

for file in ["testCD.pdf", "testCDDoc.pdf", "testCD.tex"]
  if isfile(file)
    rm(file)
  end
end

# check that the TEX file contains the desired environments
function has_environment(content::String, environment::String)
    has_begin = occursin("\\begin{$environment}", content)
    has_end = occursin("\\end{$environment}", content)
    if has_begin && has_end
        return true # has both
    elseif !has_begin && !has_end
        return false # has neither
    else
        error("\\begin{$environment} and \\end{$environment} do not match")
    end
end

# Test tikz-cd

data = L"""
A \arrow[rd] \arrow[r, "\phi"] & B \\
                               & C
"""
tp = TikzCD(data, options="scale=0.25")

td = TikzDocument()
push!(td, tp, caption="hello")

save(TEX("testCD"), tp)
@test isfile("testCD.tex")

filecontent = join(readlines("testCD.tex", keep=true)) # read with line breaks
@test occursin(strip(data[2:end-1]), filecontent) # also check that the data is contained
@test has_environment(filecontent, "tikzcd")
@test has_environment(filecontent, "document")

if success(`lualatex -v`)
    save(PDF("testCD"), tp)
    @test isfile("testCD.pdf")

    save(PDF("testCDDoc"), td)
    @test isfile("testCDDoc.pdf")
else
    @warn "lualatex is missing; can not test compilation"
end

ext: jl


text:
function get_test_trajdata()
    roadway = get_test_roadway()
    trajdata = Trajdata(roadway)

    trajdata.vehdefs[1] = VehicleDef(1, AgentClass.CAR, 5.0, 3.0)
    trajdata.vehdefs[2] = VehicleDef(2, AgentClass.CAR, 5.0, 3.0)

    push!(trajdata.states, TrajdataState(1, VehicleState(VecSE2(0.0,0.0,0.0), roadway, 10.0))) # car 1, frame 1
    push!(trajdata.states, TrajdataState(2, VehicleState(VecSE2(3.0,0.0,0.0), roadway, 20.0))) # car 2, frame 1
    push!(trajdata.states, TrajdataState(1, VehicleState(VecSE2(1.0,0.0,0.0), roadway, 10.0))) # car 1, frame 2
    push!(trajdata.states, TrajdataState(2, VehicleState(VecSE2(5.0,0.0,0.0), roadway, 20.0))) # car 2, frame 2

    push!(trajdata.frames, TrajdataFrame(1,2,0.0))
    push!(trajdata.frames, TrajdataFrame(3,4,0.1))

    trajdata
end

let
    trajdata = get_test_trajdata()
    roadway = trajdata.roadway

    @test nframes(trajdata) == 2
    @test !frame_inbounds(trajdata, 0)
    @test frame_inbounds(trajdata, 1)
    @test frame_inbounds(trajdata, 2)
    @test !frame_inbounds(trajdata, 3)

    @test carsinframe(trajdata, 1) == 2
    @test carsinframe(trajdata, 2) == 2

    @test nth_carid(trajdata, 1) == 1
    @test nth_carid(trajdata, 1, 2) == 2
    @test nth_carid(trajdata, 2, 1) == 1
    @test nth_carid(trajdata, 2, 2) == 2

    @test get_first_frame_with_id(trajdata, 1) == 1
    @test get_first_frame_with_id(trajdata, 2) == 1
    @test get_first_frame_with_id(trajdata, -1) == -1
    @test get_last_frame_with_id(trajdata, 1) == 2
    @test get_last_frame_with_id(trajdata, 2) == 2
    @test get_last_frame_with_id(trajdata, -1) == -1

    @test sort!(get_ids(trajdata)) == [1,2]

    @test iscarinframe(trajdata, 1, 1)
    @test iscarinframe(trajdata, 1, 2)
    @test iscarinframe(trajdata, 2, 1)
    @test iscarinframe(trajdata, 2, 2)
    @test !iscarinframe(trajdata, 3, 1)

    @test isapprox(get_time(trajdata, 1), 0.0)
    @test isapprox(get_time(trajdata, 2), 0.1)

    @test isapprox(get_elapsed_time(trajdata, 1, 2),  0.1)
    @test isapprox(get_elapsed_time(trajdata, 2, 1), -0.1)

    @test isapprox(get_mean_timestep(trajdata), 0.1)

    veh = get_vehicle(trajdata, 1, 1)
    @test veh.state == VehicleState(VecSE2(0.0,0.0,0.0), roadway, 10.0)
    @test_throws ErrorException get_vehicle(trajdata, 10, 1)
    @test_throws BoundsError get_vehicle(trajdata, 1, 10)

    get_vehicle!(veh, trajdata, 1, 2)
    @test veh.state == VehicleState(VecSE2(1.0,0.0,0.0), roadway, 10.0)

    let
        iter = TrajdataVehicleIterator(trajdata, 1)
        vehs = collect(iter)
        @test length(vehs) == 2
        @test vehs[1][1] == 1
        @test vehs[1][2].state == get_vehiclestate(trajdata, 1, 1)
        @test vehs[2][1] == 2
        @test vehs[2][2].state == get_vehiclestate(trajdata, 1, 2)

        iter = TrajdataVehicleIterator(trajdata, 2)
        vehs = collect(iter)
        @test length(vehs) == 2
        @test vehs[1][1] == 1
        @test vehs[1][2].state == get_vehiclestate(trajdata, 2, 1)
        @test vehs[2][1] == 2
        @test vehs[2][2].state == get_vehiclestate(trajdata, 2, 2)
    end

    path, io = mktemp()
    write(io, trajdata)
    close(io)

    lines = open(readlines, path)

    for (line_orig, line_test) in zip(lines,
            ["TRAJDATA",
             "2",
             "2 2 5.000 3.000",
             "1 2 5.000 3.000",
             "4",
             "1 (0.0000 0.0000 0.0000e+00) (1 0.0000 2 1) (0.0000 0.0000 0.0000e+00) 10.0000",
             "2 (3.0000 0.0000 0.0000e+00) (4 0.0078 2 1) (3.0078 0.0000 0.0000e+00) 20.0000",
             "1 (1.0000 0.0000 0.0000e+00) (2 0.0000 2 1) (1.0000 0.0000 0.0000e+00) 10.0000",
             "2 (5.0000 0.0000 0.0000e+00) (1 1.0000 3 1) (1.0000 0.0000 0.0000e+00) 20.0000",
             "2",
             "1 2 0.0000",
             "3 4 0.1000"]
        )

        @test strip(line_orig) == line_test
    end

    io = open(path)
    trajdata2 = read(io, Trajdata)
    close(io)
    rm(path)

    @test nframes(trajdata2) == nframes(trajdata)
    for i in 1 : nframes(trajdata2)
        @test carsinframe(trajdata2, i) == carsinframe(trajdata, i)
        for j in 1 : carsinframe(trajdata, i)
            veh1 = get_vehicle(trajdata, j, i)
            veh2 = get_vehicle(trajdata2, j, i)
            @test veh1.def.id == veh2.def.id
            @test veh1.def.class == veh2.def.class
            @test isapprox(veh1.def.length, veh2.def.length)
            @test isapprox(veh1.def.width, veh2.def.width)

            @test isapprox(veh1.state.v, veh2.state.v)
            @test isapprox(veh1.state.posG, veh2.state.posG, atol=1e-3)
            @test isapprox(veh1.state.posF.s, veh2.state.posF.s, atol=1e-3)
            @test isapprox(veh1.state.posF.t, veh2.state.posF.t, atol=1e-3)
            @test isapprox(veh1.state.posF.ϕ, veh2.state.posF.ϕ, atol=1e-6)
            @test veh1.state.posF.roadind.tag == veh2.state.posF.roadind.tag
            @test veh1.state.posF.roadind.ind.i == veh2.state.posF.roadind.ind.i
            @test isapprox(veh1.state.posF.roadind.ind.t, veh2.state.posF.roadind.ind.t, atol=1e-3)
        end
    end

    trajdata3 = Trajdata(trajdata2, 1, nframes(trajdata2))
    @test nframes(trajdata3) == nframes(trajdata2)
    for i in 1 : nframes(trajdata3)
        @test carsinframe(trajdata3, i) == carsinframe(trajdata2, i)
        for j in 1 : carsinframe(trajdata2, i)
            veh1 = get_vehicle(trajdata2, j, i)
            veh2 = get_vehicle(trajdata3, j, i)
            @test veh1.def.id == veh2.def.id
            @test veh1.def.class == veh2.def.class
            @test isapprox(veh1.def.length, veh2.def.length)
            @test isapprox(veh1.def.width, veh2.def.width)

            @test isapprox(veh1.state.v, veh2.state.v)
            @test isapprox(veh1.state.posG, veh2.state.posG, atol=1e-3)
            @test isapprox(veh1.state.posF.s, veh2.state.posF.s, atol=1e-3)
            @test isapprox(veh1.state.posF.t, veh2.state.posF.t, atol=1e-3)
            @test isapprox(veh1.state.posF.ϕ, veh2.state.posF.ϕ, atol=1e-6)
            @test veh1.state.posF.roadind.tag == veh2.state.posF.roadind.tag
            @test veh1.state.posF.roadind.ind.i == veh2.state.posF.roadind.ind.i
            @test isapprox(veh1.state.posF.roadind.ind.t, veh2.state.posF.roadind.ind.t, atol=1e-3)
        end
    end

    trajdata3 = Trajdata(trajdata2, 1, 1)
    @test nframes(trajdata3) == 1
    let
        i = 1
        @test carsinframe(trajdata3, i) == carsinframe(trajdata2, i)
        for j in 1 : carsinframe(trajdata2, i)
            veh1 = get_vehicle(trajdata2, j, i)
            veh2 = get_vehicle(trajdata3, j, i)
            @test veh1.def.id == veh2.def.id
            @test veh1.def.class == veh2.def.class
            @test isapprox(veh1.def.length, veh2.def.length)
            @test isapprox(veh1.def.width, veh2.def.width)

            @test isapprox(veh1.state.v, veh2.state.v)
            @test isapprox(veh1.state.posG, veh2.state.posG, atol=1e-3)
            @test isapprox(veh1.state.posF.s, veh2.state.posF.s, atol=1e-3)
            @test isapprox(veh1.state.posF.t, veh2.state.posF.t, atol=1e-3)
            @test isapprox(veh1.state.posF.ϕ, veh2.state.posF.ϕ, atol=1e-6)
            @test veh1.state.posF.roadind.tag == veh2.state.posF.roadind.tag
            @test veh1.state.posF.roadind.ind.i == veh2.state.posF.roadind.ind.i
            @test isapprox(veh1.state.posF.roadind.ind.t, veh2.state.posF.roadind.ind.t, atol=1e-3)
        end
    end
end


text:
using Documenter, COBS

makedocs(
    modules = [COBS],
    format = Documenter.HTML(),
    checkdocs = :exports,
    sitename = "COBS.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/yakir12/COBS.jl.git",
)


text:
"""
    metropolis_sampler(;target_pdf, known_pdf, known_sampler, starting_state, n_iter, n_burn = 100, n_thinning = 100)

Implement a metropolis independent sampler for standard boson sampling following. The burn
in period `n_burn`and the thinning interval `n_thinning` both have default value of 100.

!!! note "Reference"
    [https://arxiv.org/abs/1705.00686](https://arxiv.org/abs/1705.00686): As the paper is limited to collinionless events, we keep track of this thanks
    to [`iterate_until_collisionless`](@ref).

!!! warning
    Burn in perdiod and thinning interval must have the same value.
"""
function metropolis_sampler(;target_pdf, known_pdf, known_sampler, starting_state, n_iter, n_burn = 100, n_thinning = 100)

    function transition_probability(target_pdf, known_pdf, new_state, previous_state)

        """metropolis_independant_sampler transition probability, see eq 1 in https://arxiv.org/abs/1705.00686"""
        min(1, target_pdf(new_state) * known_pdf(previous_state)/(target_pdf(previous_state) * known_pdf(new_state)))
    end

    function do_with_probability(p)

        """returns true with probability p, false with (1-p)"""

        rand() < p ? true : false

    end

    function update_state_parameters(target_pdf, known_pdf, new_state, previous_state)

        if do_with_probability(transition_probability(target_pdf, known_pdf, new_state, previous_state))

            return new_state
        else
            return previous_state
        end

    end

    function metropolis_iteration!(target_pdf, known_pdf, known_sampler, new_state, previous_state)
        new_state[:] = known_sampler()[:]
        new_state[:] = update_state_parameters(target_pdf, known_pdf, new_state, previous_state)[:]
        previous_state[:] = new_state[:]
    end

    samples = Array{eltype(starting_state)}(undef, n_iter, size(starting_state,1))

    previous_state = similar(starting_state) #x
    new_state = similar(previous_state)
    previous_state = starting_state

    if n_burn == n_thinning
        @showprogress for iter in 1:n_iter
            for i in 1:n_thinning
                metropolis_iteration!(target_pdf, known_pdf, known_sampler, new_state, previous_state)
            end
            samples[iter,:] = new_state[:]
        end
    else
        throw(Exception("n_burn != n_thinning not implemented"))
    end

    samples

end


text:
#
# Casting to normal matrix types.
#
# This implements the casts described in README.md and the cast from
# BoundaryPaddedArray to Array.
#

function Base.copyto!(L::AbstractMatrix{T}, A::DerivativeOperator{T}, N::Int) where T
    bl = A.boundary_point_count
    stencil_length = A.stencil_length
    stencil_pivot = use_winding(A) ? (1 + stencil_length%2) : div(stencil_length,2)
    bstl = A.boundary_stencil_length

    coeff   = A.coefficients
    get_coeff = if coeff isa AbstractVector
        i -> coeff[i]
    elseif coeff isa Number
        i -> coeff
    else
        i -> true
    end

    for i in 1:bl
        cur_coeff   = get_coeff(i)
        cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.low_boundary_coefs[i]) : A.low_boundary_coefs[i]
        L[i,1:bstl] = cur_coeff * cur_stencil
    end

    for i in bl+1:N-bl
        cur_coeff   = get_coeff(i)
        stencil     = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i-bl] : A.stencil_coefs
        cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(stencil) : stencil
        L[i,i+1-stencil_pivot:i-stencil_pivot+stencil_length] = cur_coeff * cur_stencil
    end

    for i in N-bl+1:N
        cur_coeff   = get_coeff(i)
        cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
        L[i,N-bstl+3:N+2] = cur_coeff * cur_stencil
    end

    L
end

LinearAlgebra.Array(A::DerivativeOperator{T}, N::Int=A.len) where T =
    copyto!(zeros(T, N, N+2), A, N)

function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.len) where T
    bl = A.boundary_point_count
    stencil_length = A.stencil_length
    stencil_pivot = use_winding(A) ? (1 + stencil_length%2) : div(stencil_length,2)
    bstl = A.boundary_stencil_length

    coeff   = A.coefficients
    get_coeff = if coeff isa AbstractVector
        i -> coeff[i]
    elseif coeff isa Number
        i -> coeff
    else
        i -> true
    end

    Is = Int[]
    Js = Int[]
    Vs = T[]

    nvalues = 2*bl * bstl + (N - 2*bl) * stencil_length
    sizehint!(Is, nvalues)
    sizehint!(Js, nvalues)
    sizehint!(Vs, nvalues)

    for i in 1:bl
        cur_coeff   = get_coeff(i)
        cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.low_boundary_coefs[i]) : A.low_boundary_coefs[i]
        append!(Is, ((i for j in 1:bstl)...))
        append!(Js, 1:bstl)
        append!(Vs, cur_coeff * cur_stencil)
    end

    for i in bl+1:N-bl
        cur_coeff   = get_coeff(i)
        stencil     = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i-bl] : A.stencil_coefs
        cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(stencil) : stencil
        append!(Is, ((i for j in 1:stencil_length)...))
        append!(Js, i+1-stencil_pivot:i-stencil_pivot+stencil_length)
        append!(Vs, cur_coeff * cur_stencil)
    end

    for i in N-bl+1:N
        cur_coeff   = get_coeff(i)
        cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
        append!(Is, ((i for j in N-bstl+3:N+2)...))
        append!(Js, N-bstl+3:N+2)
        append!(Vs, cur_coeff * cur_stencil)
    end

   # ensure efficient allocation
    @assert length(Is) == nvalues
    @assert length(Js) == nvalues
    @assert length(Vs) == nvalues

    return sparse(Is, Js, Vs, N, N+2)
end

SparseArrays.sparse(A::DerivativeOperator{T}, N::Int=A.len) where T = SparseMatrixCSC(A,N)

Base.copyto!(L::AbstractSparseArray{T}, A::DerivativeOperator{T}, N::Int) where T =
    copyto!(L, sparse(A))

function BandedMatrices.BandedMatrix(A::DerivativeOperator{T}, N::Int=A.len) where T
    stencil_length = A.stencil_length
    bstl = A.boundary_stencil_length
    L = BandedMatrix{T}(Zeros(N, N+2), (bstl-3,bstl-1))
    copyto!(L, A, N)
end

Base.convert(::Type{Mat}, A::DerivativeOperator) where {Mat<:Union{Array,SparseMatrixCSC,BandedMatrix}} =
    Mat(A)

Base.convert(::Type{AbstractMatrix},A::DerivativeOperator) =
    BandedMatrix(A)

################################################################################
# Boundary Padded Array concretizations
################################################################################

function LinearAlgebra.Array(Q::BoundaryPaddedArray{T,D,N,M,V,B}) where {T,D,N,M,V,B}
    S = size(Q)
    out = zeros(T, S...)
    dim = D
    ulowview = selectdim(out, dim, 1)
    uhighview = selectdim(out, dim, S[dim])
    uview = selectdim(out, dim, 2:(S[dim]-1))
    ulowview .= Q.lower
    uhighview .= Q.upper
    uview .= Q.u
    return out
end

function LinearAlgebra.Array(Q::ComposedBoundaryPaddedArray{T,N,M,V,B}) where {T,N,M,V,B}
    S = size(Q)
    out = zeros(T, S...)
    dimset = 1:N
    uview = out
    for dim in dimset
        ulowview = selectdim(out, dim, 1)
        uhighview = selectdim(out, dim, S[dim])
        uview = selectdim(uview, dim, 2:(S[dim]-1))
        for (index, otherdim) in enumerate(setdiff(dimset, dim))
            ulowview = selectdim(ulowview, index, 2:(S[otherdim]-1))
            uhighview = selectdim(uhighview, index, 2:(S[otherdim]-1))
        end
        ulowview .= Q.lower[dim]
        uhighview .= Q.upper[dim]
    end
    uview .= Q.u
    return out
end

function Base.convert(::Type{AbstractArray}, A::AbstractBoundaryPaddedArray)
    Array(A)
end

################################################################################
# Boundary Condition Operator concretizations
################################################################################

#Atomic BCs
function LinearAlgebra.Array(Q::AffineBC{T}, N::Int) where {T}
    Q_L = [transpose(Q.a_l) transpose(zeros(T, N-length(Q.a_l))); Diagonal(ones(T,N)); transpose(zeros(T, N-length(Q.a_r))) transpose(Q.a_r)]
    Q_b = [Q.b_l; zeros(T,N); Q.b_r]
    return (Array(Q_L), Q_b)
end

function SparseArrays.SparseMatrixCSC(Q::AffineBC{T}, N::Int) where {T}
    Q_l = spzeros(T, N+2, N)
    for i in 1:N
        Q_l[i+1,i] = one(T)
    end

    a_r = Q.a_r[findfirst(!iszero, Q.a_r) isa Nothing ? something(end+1:end) : something(findfirst(!iszero, Q.a_r):end)]
    a_l = Q.a_l[findfirst(!iszero, Q.a_l) isa Nothing ? something(end+1:end) : something(findfirst(!iszero, Q.a_l):end)]
    for (j,e) ∈ enumerate(a_l)
        BandedMatrices.inbands_setindex!(Q_l, e, 1, j)
    end
    for (j,e) ∈ enumerate(a_r)
        BandedMatrices.inbands_setindex!(Q_l, e, N+2, N-length(a_r)+j)
    end

    Q_b = spzeros(T,N+2)
    Q_b[1] = Q.b_l
    Q_b[N+2] = Q.b_r
    return (Q_l,Q_b)
end

function BandedMatrices.BandedMatrix(Q::AffineBC{T}, N::Int) where {T}
    a_r = Q.a_r[findfirst(!iszero, Q.a_r) isa Nothing ? something(end+1:end) : something(findfirst(!iszero, Q.a_r):end)]
    a_l = Q.a_l[findfirst(!iszero, Q.a_l) isa Nothing ? something(end+1:end) : something(findfirst(!iszero, Q.a_l):end)]

    l = max(count(!iszero, a_r)+1, 1)
    u = max(count(!iszero, a_l)-1, -1)

    Q_l = BandedMatrix((-1 => ones(T,N),), (N+2,N), (l, u))
    for (j,e) ∈ enumerate(a_l)
        BandedMatrices.inbands_setindex!(Q_l, e, 1, j)
    end
    for (j,e) ∈ enumerate(a_r)
        BandedMatrices.inbands_setindex!(Q_l, e, N+2, N-length(a_r)+j)
    end
    Q_b = [Q.b_l; zeros(T,N); Q.b_r]
    return (Q_l, Q_b)
end

function SparseArrays.sparse(Q::AffineBC{T}, N::Int) where {T}
    SparseMatrixCSC(Q,N)
end

LinearAlgebra.Array(Q::PeriodicBC{T}, N::Int) where T = (Array([transpose(zeros(T, N-1)) one(T); Diagonal(ones(T,N)); one(T) transpose(zeros(T, N-1))]), zeros(T, N+2))
function SparseArrays.SparseMatrixCSC(Q::PeriodicBC{T}, N::Int) where T
    Q_l = spzeros(T,N+2,N)
    for i in 1:N
       Q_l[i+1,i] = one(T)
    end
    Q_l[1, end] = one(T)
    Q_l[end, 1] = one(T)
    Q_b = spzeros(T,N+2)
    return (Q_l,Q_b)
end
SparseArrays.sparse(Q::PeriodicBC{T}, N::Int) where T = SparseMatrixCSC(Q,N)
function BandedMatrices.BandedMatrix(Q::PeriodicBC{T}, N::Int) where T #Not recommended!
    Q_array = BandedMatrix{T}((-1 => ones(T,N),), (N+2, N), (N+1,N+1))
    Q_array[1, end] = one(T)
    Q_array[end, 1] = one(T)
    return (Q_array, zeros(T, N+2))
end

function LinearAlgebra.Array(Q::BoundaryPaddedVector)
    return [Q.l; Q.u; Q.r]
end

function Base.convert(::Type{Array},A::AbstractBC{T}) where T
    Array(A)
end

function Base.convert(::Type{SparseMatrixCSC},A::AbstractBC{T}) where T
    SparseMatrixCSC(A)
end

function Base.convert(::Type{AbstractMatrix},A::AbstractBC{T}) where T
    SparseMatrixCSC(A)
end

# Multidimensional BC operators
_concretize(Q::MultiDimDirectionalBC, M) = _concretize(Q.BCs, M)

function _concretize(Q::AbstractArray{T,N}, M) where {T,N}
    return (stencil.(Q, fill(M,size(Q))), affine.(Q))
end

function LinearAlgebra.Array(Q::MultiDimDirectionalBC{T, B, D, N, L}, s::NTuple{N,G}) where {T, B, D,N,L, G<:Int}
    @assert size(Q.BCs) == perpindex(s, D) "The size of the BC array in Q, $(size(Q.BCs)) is incompatible with s, $s"
    blip = zeros(Int64, N)
    blip[D] = 2
    s_pad = s.+ blip # extend s in the right direction
    Q = _concretize(Q.BCs, s[D])
    ē = unit_indices(N)
    QL = zeros(T, prod(s_pad), prod(s))
    Qb = zeros(T, prod(s_pad))
    ranges = Union{typeof(1:10), Int64}[1:s[i] for i in 1:N]
    ranges[D] = ranges[D] .+ 1

    interior = CartesianIndices(Tuple(ranges))
    I1 = CartesianIndex(Tuple(ones(Int64, N)))
    for I in interior
        i = cartesian_to_linear(I, s_pad)
        j = cartesian_to_linear(I-ē[D], s)
        QL[i,j] = one(T)
    end
    ranges[D] = 1
    lower = CartesianIndices((Tuple(ranges)))
    ranges[D] = s_pad[D]
    upper = CartesianIndices((Tuple(ranges)))
    for K in CartesianIndices(upper)
        I = CartesianIndex(Tuple(K)[setdiff(1:N, D)])
        il = cartesian_to_linear(lower[K], s_pad)
        iu = cartesian_to_linear(upper[K], s_pad)
        Qb[il] = Q[2][I][1]
        Qb[iu] = Q[2][I][2]
        for k in 0:s[D]-1
            j = cartesian_to_linear(lower[K] + k*ē[D], s)
            QL[il, j] = Q[1][I][1][k+1]
            QL[iu, j] = Q[1][I][2][k+1]
        end
    end

    return (QL, Qb)
end

"""
This is confusing, but it does work.
"""
function LinearAlgebra.Array(Q::ComposedMultiDimBC{T, B, N,M} , s::NTuple{N,G}) where {T, B, N, M, G<:Int}
    for d in 1:N
        @assert size(Q.BCs[d]) == perpindex(s, d) "The size of the BC array in Q along dimension $d, $(size(Q.BCs[d])) is incompatible with s, $s"
    end
    s_pad = s.+2
    Q = Tuple(_concretize.(Q.BCs, s)) #essentially, finding the first and last rows of the matrix part and affine part for every atomic BC

    QL = zeros(T, prod(s_pad), prod(s))
    Qb = zeros(T, prod(s_pad))

    ranges = Union{typeof(1:10), Int64}[2:s_pad[i]-1 for i in 1:N] #Set up indices corresponding to the interior
    interior = CartesianIndices(Tuple(ranges))

    ē = unit_indices(N) #setup unit indices in each direction
    I1 = CartesianIndex(Tuple(ones(Int64, N))) #set up the ones index
    for I in interior #loop over interior
        i = cartesian_to_linear(I, s_pad) #find the index on the padded side
        j = cartesian_to_linear(I-I1, s)  #find the index on the unpadded side
        QL[i,j] = one(T)  #create a padded identity matrix
    end
    for dim in 1:N #Loop over boundaries
        r_ = deepcopy(ranges)
        r_[dim] = 1
        lower = CartesianIndices((Tuple(r_))) #set up upper and lower indices
        r_[dim] = s_pad[dim]
        upper = CartesianIndices((Tuple(r_)))
        for K in CartesianIndices(upper) #for every element of the boundaries
            I = CartesianIndex(Tuple(K)[setdiff(1:N, dim)]) #convert K to 2D index for indexing the BC arrays
            il = cartesian_to_linear(lower[K], s_pad) #Translate to linear indices
            iu = cartesian_to_linear(upper[K], s_pad) # ditto
            Qb[il] = Q[dim][2][I][1] #store the affine parts in indices corresponding with the lower index boundary
            Qb[iu] = Q[dim][2][I][2] #ditto with upper index
            for k in 1:s[dim] #loop over the direction orthogonal to the boundary
                j = cartesian_to_linear(lower[K] + k*ē[dim]-I1, s) #Find the linear index this element of the boundary stencil should be at on the unpadded side
                QL[il, j] = Q[dim][1][I][1][k]
                QL[iu, j] = Q[dim][1][I][2][k]
            end
        end
    end

    return (QL, Qb)
end

"""
See comments on the `Array` method for this type for an idea of what is going on.
"""
function SparseArrays.SparseMatrixCSC(Q::MultiDimDirectionalBC{T, B, D, N, L}, s::NTuple{N,G}) where {T, B, D,N,L, G<:Int}
    @assert size(Q.BCs) == perpindex(s, D) "The size of the BC array in Q, $(size(Q.BCs)) is incompatible with s, $s"
    blip = zeros(Int64, N)
    blip[D] = 2
    s_pad = s.+ blip # extend s in the right direction
    Q = _concretize(Q.BCs, s[D])
    ē = unit_indices(N)
    QL = spzeros(T, prod(s_pad), prod(s))
    Qb = spzeros(T, prod(s_pad))
    ranges = Union{typeof(1:10), Int64}[1:s[i] for i in 1:N]
    ranges[D] = ranges[D] .+ 1

    interior = CartesianIndices(Tuple(ranges))
    I1 = CartesianIndex(Tuple(ones(Int64, N)))
    for I in interior
        i = cartesian_to_linear(I, s_pad)
        j = cartesian_to_linear(I-ē[D], s)
        QL[i,j] = one(T)
    end
    ranges[D] = 1
    lower = CartesianIndices((Tuple(ranges)))
    ranges[D] = s_pad[D]
    upper = CartesianIndices((Tuple(ranges)))
    for K in CartesianIndices(upper)
        I = CartesianIndex(Tuple(K)[setdiff(1:N, D)])
        il = cartesian_to_linear(lower[K], s_pad)
        iu = cartesian_to_linear(upper[K], s_pad)
        Qb[il] = Q[2][I][1]
        Qb[iu] = Q[2][I][2]
        for k in 0:s[D]-1
            j = cartesian_to_linear(lower[K] + k*ē[D], s)
            QL[il, j] = Q[1][I][1][k+1]
            QL[iu, j] = Q[1][I][2][k+1]
        end
    end

    return (QL, Qb)
end

function SparseArrays.SparseMatrixCSC(Q::ComposedMultiDimBC{T, B, N,M}, s::NTuple{N,G}) where {T, B, N, M, G<:Int}
    for d in 1:N
        @assert size(Q.BCs[d]) == perpindex(s, d) "The size of the BC array in Q along dimension $d, $(size(Q.BCs[d])) is incompatible with s, $s"
    end
    s_pad = s.+2
    Q = Tuple(_concretize.(Q.BCs, s)) #essentially, finding the first and last rows of the matrix part and affine part for every atomic BC

    QL = spzeros(T, prod(s_pad), prod(s))
    Qb = spzeros(T, prod(s_pad))

    ranges = Union{typeof(1:10), Int64}[2:s_pad[i]-1 for i in 1:N] #Set up indices corresponding to the interior
    interior = CartesianIndices(Tuple(ranges))

    ē = unit_indices(N) #setup unit indices in each direction
    I1 = CartesianIndex(Tuple(ones(Int64, N))) #set up the ones index
    for I in interior #loop over interior
        i = cartesian_to_linear(I, s_pad) #find the index on the padded side
        j = cartesian_to_linear(I-I1, s)  #find the index on the unpadded side
        QL[i,j] = one(T)  #create a padded identity matrix
    end
    for dim in 1:N #Loop over boundaries
        r_ = deepcopy(ranges)
        r_[dim] = 1
        lower = CartesianIndices((Tuple(r_))) #set up upper and lower indices
        r_[dim] = s_pad[dim]
        upper = CartesianIndices((Tuple(r_)))
        for K in CartesianIndices(upper) #for every element of the boundaries
            I = CartesianIndex(Tuple(K)[setdiff(1:N, dim)]) #convert K to 2D index for indexing the BC arrays
            il = cartesian_to_linear(lower[K], s_pad) #Translate to linear indices
            iu = cartesian_to_linear(upper[K], s_pad) # ditto
            Qb[il] = Q[dim][2][I][1] #store the affine parts in indices corresponding with the lower index boundary
            Qb[iu] = Q[dim][2][I][2] #ditto with upper index
            for k in 1:s[dim] #loop over the direction orthogonal to the boundary
                j = cartesian_to_linear(lower[K] + k*ē[dim]-I1, s) #Find the linear index this element of the boundary stencil should be at on the unpadded side
                QL[il, j] = Q[dim][1][I][1][k]
                QL[iu, j] = Q[dim][1][I][2][k]
            end
        end
    end

    return (QL, Qb)
end

SparseArrays.sparse(Q::MultiDimDirectionalBC, s) = SparseMatrixCSC(Q, s)
SparseArrays.sparse(Q::ComposedMultiDimBC, s) = SparseMatrixCSC(Q, s)

function BandedMatrices.BandedMatrix(Q:: MultiDimensionalBC, M) where {T, B, D,N,K}
    throw("Banded Matrix cocnretization not yet supported for MultiDimensionalBCs")
end

################################################################################
# Higher-Dimensional DerivativeOperator Concretizations
################################################################################
# Higher-Dimensional Concretizations. The following concretizations return two-dimensional arrays
# which operate on flattened vectors. Mshape is the size of the unflattened array on which A is operating.

function LinearAlgebra.Array(A::DerivativeOperator{T,N}, Mshape) where {T,N}
    # Case where A is not differentiating along the first dimension
    if N != 1
        n = prod(Mshape[1:N-1])
        B = Kron(Array(A), Eye(n))
        if N != length(Mshape)
            n = prod(Mshape[N+1:end])
            B = Kron(Eye(n), B)
        end

    # Case where A is differentiating along the first dimension
    else
        n = prod(Mshape[2:end])
        B = Kron(Eye(n), Array(A))
    end
    return Array(B)
end

function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N}, Mshape) where {T,N}
    # Case where A is not differentiating along the first dimension
    if N != 1
        n = prod(Mshape[1:N-1])
        B = Kron(sparse(A), sparse(I,n,n))
        if N != length(Mshape)
            n = prod(Mshape[N+1:end])
            B = Kron(sparse(I,n,n), B)
        end

    # Case where A is differentiating along the first dimension
    else
        n = prod(Mshape[2:end])
        B = Kron(sparse(I,n,n), sparse(A))
    end
    return sparse(B)
end

function SparseArrays.sparse(A::DerivativeOperator{T,N}, Mshape) where {T,N}
    return SparseMatrixCSC(A,Mshape)
end

function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N}, Mshape) where {T,N}
    # Case where A is not differentiating along the first dimension
    if N != 1
        n = prod(Mshape[1:N-1])
        B = Kron(BandedMatrix(A), Eye(n))
        if N != length(Mshape)
            n = prod(Mshape[N+1:end])
            B = Kron(Eye(n), B)
        end

    # Case where A is differentiating along the first dimension
    else
        n = prod(Mshape[2:end])
        B = Kron(BandedMatrix(Eye(n)), BandedMatrix(A))
    end
    return BandedMatrix(B)
end

function BlockBandedMatrices.BandedBlockBandedMatrix(A::DerivativeOperator{T,N}, Mshape) where {T,N}
    # Case where A is not differentiating along the first dimension
    if N != 1
        n = prod(Mshape[1:N-1])
        B = Kron(BandedMatrix(A), Eye(n))
        if N != length(Mshape)
            n = prod(Mshape[N+1:end])
            B = Kron(Eye(n), B)
        end

    # Case where A is differentiating along the first dimension
    else
        n = prod(Mshape[2:end])
        B = Kron(BandedMatrix(Eye(n)), BandedMatrix(A))
    end
    return BandedBlockBandedMatrix(B)
end

################################################################################
# Upwind Operator Concretization
################################################################################

# Array Concretizations
# Uniform grid case
function LinearAlgebra.Array(A::DerivativeOperator{T,N,true}, len::Int=A.len) where {T,N}
    L = zeros(T, len, len+2)
    bpc = A.boundary_point_count
    stl = A.stencil_length
    bstl = A.boundary_stencil_length
    coeff   = A.coefficients
    offside = A.offside

    # downwind stencils at low boundary
    downwind_stencils = A.low_boundary_coefs
    # upwind stencils at upper boundary
    upwind_stencils = A.high_boundary_coefs
    # interior stencils
    stencils = A.stencil_coefs

    for i in 1:bpc
        cur_coeff   = coeff[i]
        if cur_coeff >= 0 && (i+stl <= len+2) && i >= offside
            cur_stencil = stencils
            L[i,i+1-offside:i+stl-offside] = cur_coeff*cur_stencil
        else
            cur_stencil = downwind_stencils[i]
            L[i,1:bstl] = cur_coeff * cur_stencil
        end
    end

    for i in bpc+1 : len-bpc-offside
        cur_coeff   = coeff[i]
        cur_stencil = stencils
        cur_stencil = cur_coeff >= 0 ? cur_stencil : ((-1)^A.derivative_order)*reverse(cur_stencil)
        if cur_coeff >= 0
            L[i,i+1-offside:i+stl-offside] = cur_coeff * cur_stencil
        else
            L[i,i-stl+2+offside:i+1+offside] = cur_coeff * cur_stencil
        end
    end

    for i in 1 : bpc+offside
        cur_coeff   = coeff[len-bpc+i-offside]
        if cur_coeff < 0 && (len+2-stl-bpc+i >= 1) && (i <= bpc+1)
            cur_stencil = stencils
            cur_stencil = ((-1)^A.derivative_order)*reverse(cur_stencil)
            L[len-bpc+i-offside,i+len-stl-bpc+2:i+len-bpc+1] = cur_coeff * cur_stencil
        elseif cur_coeff < 0 && (len+2-stl-bpc+i >= 1) && (i > bpc + 1)
            cur_stencil = upwind_stencils[bpc + offside + 1 - i]
            cur_stencil = ((-1)^A.derivative_order)*reverse(cur_stencil)
            L[len-bpc+i-offside ,len-bstl+3:len+2] = cur_coeff * cur_stencil
        elseif cur_coeff >=0 && i < offside + 1
            cur_stencil = stencils
            L[len-bpc+i-offside,len-stl+i-offside+3:len+i-offside+2] = cur_coeff * cur_stencil
        else
            cur_stencil = upwind_stencils[i]
            L[len-bpc+i-offside,len-bstl+3:len+2] = cur_coeff * cur_stencil
        end
    end
    return L
end

# Non-uniform grid case
function LinearAlgebra.Array(A::DerivativeOperator{T,N,true,M}, len::Int=A.len) where {T,N,M<:AbstractArray{T}}
    L = zeros(T, len, len+2)
    bpc = A.boundary_point_count
    stl = A.stencil_length
    bstl = A.boundary_stencil_length
    coeff   = A.coefficients
    offside = A.offside

    for i in 1:bpc
        cur_coeff   = coeff[i]
        if cur_coeff >= 0 && offside == 0
            L[i,i+1:i+stl] = cur_coeff * A.low_boundary_coefs[1,i]
        elseif cur_coeff >= 0 && i < offside
            L[i,1:stl] = cur_coeff * A.low_boundary_coefs[1,i]
        elseif cur_coeff >= 0 && i >= offside
            L[i,i+1-offside:i+stl-offside] = cur_coeff * A.low_boundary_coefs[1,i]
        else
            L[i,1:bstl] = cur_coeff * A.low_boundary_coefs[2,i]
        end
    end

    for i in bpc+1:len-bpc-offside
        cur_coeff   = coeff[i]
        if cur_coeff >= 0
            L[i,i+1-offside:i+stl-offside] = cur_coeff * A.stencil_coefs[1,i-bpc]
        else
            L[i,i-stl+2+offside:i+1+offside] = cur_coeff * A.stencil_coefs[2,i-bpc]
        end
    end

    for i in len-bpc+1-offside:len
        cur_coeff   = coeff[i]
        if cur_coeff < 0 && i <= len - offside
            L[i,i-stl+2+offside:i+1+offside] = cur_coeff * A.high_boundary_coefs[2,i-len+bpc+offside]
        elseif cur_coeff < 0
            L[i,len-stl+3:len+2] = cur_coeff * A.high_boundary_coefs[2,i-len+bpc+offside]
        elseif cur_coeff >=0 && i <= len - bpc
            L[i,i-bstl+2+offside:i+1+offside] = cur_coeff * A.stencil_coefs[1,i-bpc]
        else
            L[i,len-bstl+3:len+2] = cur_coeff * A.high_boundary_coefs[1,i-len+bpc+offside]
        end
    end
    return L
end

# Sparse Concretizations
# Uniform grid case
function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true}, len::Int=A.len) where {T,N}
    L = spzeros(T, len, len+2)
    bpc = A.boundary_point_count
    stl = A.stencil_length
    bstl = A.boundary_stencil_length
    coeff   = A.coefficients
    offside = A.offside

    # downwind stencils at low boundary
    downwind_stencils = A.low_boundary_coefs
    # upwind stencils at upper boundary
    upwind_stencils = A.high_boundary_coefs
    # interior stencils
    stencils = A.stencil_coefs

    for i in 1:bpc
        cur_coeff   = coeff[i]
        if cur_coeff >= 0 && (i+stl <= len+2) && i >= offside
            cur_stencil = stencils
            L[i,i+1-offside:i+stl-offside] = cur_coeff*cur_stencil
        else
            cur_stencil = downwind_stencils[i]
            L[i,1:bstl] = cur_coeff * cur_stencil
        end
    end

    for i in bpc+1 : len-bpc-offside
        cur_coeff   = coeff[i]
        cur_stencil = stencils
        cur_stencil = cur_coeff >= 0 ? cur_stencil : ((-1)^A.derivative_order)*reverse(cur_stencil)
        if cur_coeff >= 0
            L[i,i+1-offside:i+stl-offside] = cur_coeff * cur_stencil
        else
            L[i,i-stl+2+offside:i+1+offside] = cur_coeff * cur_stencil
        end
    end

    for i in 1 : bpc+offside
        cur_coeff   = coeff[len-bpc+i-offside]
        if cur_coeff < 0 && (len+2-stl-bpc+i >= 1) && (i <= bpc+1)
            cur_stencil = stencils
            cur_stencil = ((-1)^A.derivative_order)*reverse(cur_stencil)
            L[len-bpc+i-offside,i+len-stl-bpc+2:i+len-bpc+1 ] = cur_coeff * cur_stencil
        elseif cur_coeff < 0 && (len+2-stl-bpc+i >= 1) && (i > bpc + 1)
            cur_stencil = upwind_stencils[bpc + offside + 1 - i]
            cur_stencil = ((-1)^A.derivative_order)*reverse(cur_stencil)
            L[len-bpc+i-offside,len-bstl+3:len+2] = cur_coeff * cur_stencil
        elseif cur_coeff >=0 && i < offside + 1
            cur_stencil = stencils
            L[len-bpc+i-offside,len-stl+i-offside+3:len+i-offside+2] = cur_coeff * cur_stencil
        else
            cur_stencil = upwind_stencils[i]
            L[len-bpc+i-offside,len-bstl+3:len+2] = cur_coeff * cur_stencil
        end
    end
    return L
end

# Non-uniform grid case
function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T,N,true,M}, len::Int=A.len) where {T,N,M<:AbstractArray{T}}
    L = spzeros(T, len, len+2)
    bpc = A.boundary_point_count
    stl = A.stencil_length
    bstl = A.boundary_stencil_length
    coeff   = A.coefficients
    offside = A.offside

    for i in 1:bpc
        cur_coeff   = coeff[i]
        if cur_coeff >= 0 && offside == 0
            L[i,i+1:i+stl] = cur_coeff * A.low_boundary_coefs[1,i]
        elseif cur_coeff >= 0 && i < offside
            L[i,1:stl] = cur_coeff * A.low_boundary_coefs[1,i]
        elseif cur_coeff >= 0 && i >= offside
            L[i,i+1-offside:i+stl-offside] = cur_coeff * A.low_boundary_coefs[1,i]
        else
            L[i,1:bstl] = cur_coeff * A.low_boundary_coefs[2,i]
        end
    end

    for i in bpc+1:len-bpc-offside
        cur_coeff   = coeff[i]
        if cur_coeff >= 0
            L[i,i+1-offside:i+stl-offside] = cur_coeff * A.stencil_coefs[1,i-bpc]
        else
            L[i,i-stl+2+offside:i+1+offside] = cur_coeff * A.stencil_coefs[2,i-bpc]
        end
    end

    for i in len-bpc+1-offside:len
        cur_coeff   = coeff[i]
        if cur_coeff < 0 && i <= len - offside
            L[i,i-stl+2+offside:i+1+offside] = cur_coeff * A.high_boundary_coefs[2,i-len+bpc+offside]
        elseif cur_coeff < 0
            L[i,len-stl+3:len+2] = cur_coeff * A.high_boundary_coefs[2,i-len+bpc+offside]
        elseif cur_coeff >=0 && i <= len - bpc
            L[i,i-bstl+2+offside:i+1+offside] = cur_coeff * A.stencil_coefs[1,i-bpc]
        else
            L[i,len-bstl+3:len+2] = cur_coeff * A.high_boundary_coefs[1,i-len+bpc+offside]
        end
    end
    return L
end

# Banded Concretizations
# Uniform grid case
function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true}, len::Int=A.len) where {T,N}
    bpc = A.boundary_point_count
    stl = A.stencil_length
    bstl = A.boundary_stencil_length
    coeff   = A.coefficients
    offside = A.offside

    L = BandedMatrix{T}(Zeros(len, len+2), (stl-2, stl))

    # downwind stencils at low boundary
    downwind_stencils = A.low_boundary_coefs
    # upwind stencils at upper boundary
    upwind_stencils = A.high_boundary_coefs
    # interior stencils
    stencils = A.stencil_coefs

    for i in 1:bpc
        cur_coeff   = coeff[i]
        if cur_coeff >= 0 && (i+stl <= len+2) && i >= offside
            cur_stencil = stencils
            L[i,i+1-offside:i+stl-offside] = cur_coeff*cur_stencil
        else
            cur_stencil = downwind_stencils[i]
            L[i,1:bstl] = cur_coeff * cur_stencil
        end
    end

    for i in bpc+1 : len-bpc-offside
        cur_coeff   = coeff[i]
        cur_stencil = stencils
        cur_stencil = cur_coeff >= 0 ? cur_stencil : ((-1)^A.derivative_order)*reverse(cur_stencil)
        if cur_coeff >= 0
            L[i,i+1-offside:i+stl-offside] = cur_coeff * cur_stencil
        else
            L[i,i-stl+2+offside:i+1+offside] = cur_coeff * cur_stencil
        end
    end

    for i in 1 : bpc+offside
        cur_coeff   = coeff[len-bpc+i-offside]
        if cur_coeff < 0 && (len+2-stl-bpc+i >= 1) && (i <= bpc+1)
            cur_stencil = stencils
            cur_stencil = ((-1)^A.derivative_order)*reverse(cur_stencil)
            L[len-bpc+i-offside,i+len-stl-bpc+2:i+len-bpc+1] = cur_coeff * cur_stencil
        elseif cur_coeff < 0 && (len+2-stl-bpc+i >= 1) && (i > bpc + 1)
            cur_stencil = upwind_stencils[bpc + offside + 1 - i]
            cur_stencil = ((-1)^A.derivative_order)*reverse(cur_stencil)
            L[len-bpc+i-offside,len-bstl+3:len+2] = cur_coeff * cur_stencil
        elseif cur_coeff >=0 && i < offside + 1
            cur_stencil = stencils
            L[len-bpc+i-offside,len-stl+i-offside+3:len+i-offside+2] = cur_coeff * cur_stencil
        else
            cur_stencil = upwind_stencils[i]
            L[len-bpc+i-offside,len-bstl+3:len+2] = cur_coeff * cur_stencil
        end
    end
    return L
end

# Non-uniform grid case
function BandedMatrices.BandedMatrix(A::DerivativeOperator{T,N,true,M}, len::Int=A.len) where {T,N,M<:AbstractArray{T}}
    bpc = A.boundary_point_count
    stl = A.stencil_length
    bstl = A.boundary_stencil_length
    coeff   = A.coefficients
    offside = A.offside
    L = BandedMatrix{T}(Zeros(len, len+2), (stl-2, stl))

    for i in 1:bpc
        cur_coeff   = coeff[i]
        if cur_coeff >= 0 && offside == 0
            L[i,i+1:i+stl] = cur_coeff * A.low_boundary_coefs[1,i]
        elseif cur_coeff >= 0 && i < offside
            L[i,1:stl] = cur_coeff * A.low_boundary_coefs[1,i]
        elseif cur_coeff >= 0 && i >= offside
            L[i,i+1-offside:i+stl-offside] = cur_coeff * A.low_boundary_coefs[1,i]
        else
            L[i,1:bstl] = cur_coeff * A.low_boundary_coefs[2,i]
        end
    end

    for i in bpc+1:len-bpc-offside
        cur_coeff   = coeff[i]
        if cur_coeff >= 0
            L[i,i+1-offside:i+stl-offside] = cur_coeff * A.stencil_coefs[1,i-bpc]
        else
            L[i,i-stl+2+offside:i+1+offside] = cur_coeff * A.stencil_coefs[2,i-bpc]
        end
    end

    for i in len-bpc+1-offside:len
        cur_coeff   = coeff[i]
        if cur_coeff < 0 && i <= len - offside
            L[i,i-stl+2+offside:i+1+offside] = cur_coeff * A.high_boundary_coefs[2,i-len+bpc+offside]
        elseif cur_coeff < 0
            L[i,len-stl+3:len+2] = cur_coeff * A.high_boundary_coefs[2,i-len+bpc+offside]
        elseif cur_coeff >=0 && i <= len - bpc
            L[i,i-bstl+2+offside:i+1+offside] = cur_coeff * A.stencil_coefs[1,i-bpc]
        else
            L[i,len-bstl+3:len+2] = cur_coeff * A.high_boundary_coefs[1,i-len+bpc+offside]
        end
    end
    return L
end

# GhostDerivativeOperator Concretizations
################################################################################
function LinearAlgebra.Array(A::GhostDerivativeOperator{T, E, F},N::Int=A.L.len) where {T,E,F}
    Q_l, Q_b = Array(A.Q,N)
    return (Array(A.L,N)*Q_l, Array(A.L,N)*Q_b)
end

function LinearAlgebra.Array(A::GhostDerivativeOperator{T, E, F}, s::NTuple{N,I}) where {T,E,F,N,I<:Int}
    Q_l, Q_b = Array(A.Q,s)
    return (Array(A.L,s)*Q_l, Array(A.L,s)*Q_b)
end

function BandedMatrices.BandedMatrix(A::GhostDerivativeOperator{T, E, F},N::Int=A.L.len) where {T,E,F}
    Q_l, Q_b = BandedMatrix(A.Q,N)
    return (BandedMatrix(A.L,N)*Q_l, BandedMatrix(A.L,N)*Q_b)
end

function BandedMatrices.BandedMatrix(A::GhostDerivativeOperator{T, E, F}, s::NTuple{N,I}) where {T,E,F, N, I<:Int}
    Q_l,Q_b = BandedMatrix(A.Q,s)
    return (BandedMatrix(A.L,s)*Q_l, BandedMatrix(A.L,s)*Q_b)
end

function SparseArrays.SparseMatrixCSC(A::GhostDerivativeOperator{T, E, F},N::Int=A.L.len) where {T,E,F}
    Q_l, Q_b = SparseMatrixCSC(A.Q,N)
    return (SparseMatrixCSC(A.L,N)*Q_l, SparseMatrixCSC(A.L,N)*Q_b)
end

function SparseArrays.SparseMatrixCSC(A::GhostDerivativeOperator{T, E, F}, s::NTuple{N,I}) where {T,E,F,N,I<:Int}
    Q_l, Q_b = SparseMatrixCSC(A.Q,s)
    return (SparseMatrixCSC(A.L,s)*Q_l, SparseMatrixCSC(A.L,s)*Q_b)
end

function SparseArrays.sparse(A::GhostDerivativeOperator{T, E, F}, N::Int=A.L.len) where {T,E,F}
    return SparseMatrixCSC(A,N)
end

function SparseArrays.sparse(A::GhostDerivativeOperator{T, E, F}, s::NTuple{N,I}) where {T,E,F, N,I}
    return SparseMatrixCSC(A,s)
end

################################################################################
# Composite Operator Concretizations
################################################################################
Array(L::DiffEqScaledOperator, s) = L.coeff * Array(L.op, s)
Array(L::DiffEqOperatorCombination, s) = sum(Array.(L.ops, fill(s, length(L.ops))))
Array(L::DiffEqOperatorComposition, s) = prod(Array.(reverse(L.ops), fill(s, length(L.ops))))


text:
# benchmark suite execution and codespeed submission

using CUDA, BenchmarkTools

BenchmarkTools.DEFAULT_PARAMETERS.evals = 0     # to find untuned benchmarks

SUITE = BenchmarkGroup()

# NOTE: don't use spaces in benchmark names (tobami/codespeed#256)

include("kernel.jl")
include("array.jl")

@info "Warming up"
warmup(SUITE; verbose=false)

paramsfile = joinpath(first(DEPOT_PATH), "cache", "CUDA_benchmark_params.json")
mkpath(dirname(paramsfile))
if !isfile(paramsfile)
    @warn "No saved parameters found, will re-tune all benchmarks"
    tune!(SUITE)
else
    loadparams!(SUITE, BenchmarkTools.load(paramsfile)[1], :evals, :samples)

    # find untuned benchmarks for which we have the default evals==0
    function find_untuned(group::BenchmarkGroup, untuned=Dict(), prefix="")
        for (name, b) in group
            find_untuned(b, untuned, isempty(prefix) ? name : "$prefix/$name")
        end
        return untuned
    end
    function find_untuned(b::BenchmarkTools.Benchmark, untuned=Dict(), prefix="")
        if params(b).evals == 0
            untuned[prefix] = b
        end
        return untuned
    end
    untuned = find_untuned(SUITE)

    if !isempty(untuned)
        @info "Tuning parameters: $(join(keys(untuned), ", "))"
        foreach(tune!, values(untuned))
        BenchmarkTools.save(paramsfile, params(SUITE))
    end
end

@info "Running benchmarks"
results = run(SUITE, verbose=true)
println(results)

## submission

using JSON, HTTP

if get(ENV, "CODESPEED_BRANCH", nothing) == "master"
    @info "Submitting to Codespeed..."

    basedata = Dict(
        "branch"        => ENV["CODESPEED_BRANCH"],
        "commitid"      => ENV["CODESPEED_COMMIT"],
        "project"       => ENV["CODESPEED_PROJECT"],
        "environment"   => ENV["CODESPEED_ENVIRONMENT"],
        "executable"    => ENV["CODESPEED_EXECUTABLE"]
    )

    # convert nested groups of benchmark to flat dictionaries of results
    flat_results = []
    function flatten(results, prefix="")
        for (key,value) in results
            if value isa BenchmarkGroup
                flatten(value, "$prefix$key/")
            else
                @assert value isa BenchmarkTools.Trial

                # codespeed reports maxima, but those are often very noisy.
                # get rid of measurements that unnecessarily skew the distribution.
                rmskew!(value)

                push!(flat_results,
                    Dict(basedata...,
                        "benchmark" => "$prefix$key",
                        "result_value" => median(value).time / 1e9,
                        "min" => minimum(value).time / 1e9,
                        "max" => maximum(value).time / 1e9))
            end
        end
    end
    flatten(results)

    HTTP.post("$(ENV["CODESPEED_SERVER"])/result/add/json/",
                ["Content-Type" => "application/x-www-form-urlencoded"],
                HTTP.URIs.escapeuri(Dict("json" => JSON.json(flat_results))))
end


text:
import RowEchelon
import LinearAlgebra

#=
Homework 6
Section 2.5
Problems 17
=#

# Create the matrix
A = [2 -6 4 ; -4 8 0 ; 0 -4 6]
LU = LinearAlgebra.lu(A)
print("L = \t")
println(LU.L)
print("U = \t")
println(LU.U)

# Create inverses of LU
invL = inv(LU.L)
invU = inv(LU.U)
print("invL = \t")
println(invL)
print("invU = \t")
println(invU)

# Solve for invA
print("invA = \t")
println(invU*invL)


text:
using Clang

includes = ["/usr/local/include",
            "/usr/include",
            "/usr/lib/gcc/x86_64-linux-gnu/4.8/include",
            "/usr/lib/gcc/x86_64-linux-gnu/4.8/include-fixed"]
headers = ["/usr/local/cuda/include/nvrtc.h"]

# Customize how functions, constants, and structs are written
const skip_expr = []
const skip_error_check = []
function rewriter(ex::Expr)
    if in(ex, skip_expr)
        return :()
    end
    # Empty types get converted to Void
    if ex.head == :type
        a3 = ex.args[3]
        if isempty(a3.args)
            objname = ex.args[2]
            return :(typealias $objname Void)
        end
    end
    ex.head == :function || return ex
    decl, body = ex.args[1], ex.args[2]
    # omit types from function prototypes
    for i = 2:length(decl.args)
        a = decl.args[i]
        if a.head == :(::)
            decl.args[i] = a.args[1]
        end
    end
    # Error-check functions that return a cudaError_t (with some omissions)
    ccallexpr = body.args[1]
    if ccallexpr.head != :ccall
      error("Unexpected body expression: ", body)
    end
    rettype = ccallexpr.args[2]
    if rettype == :nvrtcResult
      fname = decl.args[1]
      if !in(fname, skip_error_check)
        body.args[1] = Expr(:call, :check_nvrtcresult, deepcopy(ccallexpr))
      end
    end
    ex
end

rewriter(A::Array) = [rewriter(a) for a in A]
rewriter(s::Symbol) = string(s)
rewriter(arg) = arg

context = wrap_c.init(output_file = "libnvrtc.jl",
                      common_file = "libnvrtc_types.jl",
                      header_library = _->"libnvrtc",
                      headers = headers,
                      clang_includes = includes,
                      clang_diagnostics = true,
                      header_wrapped=(x,y)->contains(y,"nvrtc"),
                      rewriter = rewriter)

context.options = wrap_c.InternalOptions(true,true)  # wrap structs, too

# Execute the wrap
run(context)


text:
function AdapterDetector(arg0::EventDetector)
    return AdapterDetector((EventDetector,), arg0)
end

function event_occurred(obj::AdapterDetector, arg0::SpacecraftState, arg1::jboolean)
    return jcall(obj, "eventOccurred", Action, (SpacecraftState, jboolean), arg0, arg1)
end

function g(obj::AdapterDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_detector(obj::AdapterDetector)
    return jcall(obj, "getDetector", EventDetector, ())
end

function get_max_check_interval(obj::AdapterDetector)
    return jcall(obj, "getMaxCheckInterval", jdouble, ())
end

function get_max_iteration_count(obj::AdapterDetector)
    return jcall(obj, "getMaxIterationCount", jint, ())
end

function get_threshold(obj::AdapterDetector)
    return jcall(obj, "getThreshold", jdouble, ())
end

function init(obj::AdapterDetector, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end

function reset_state(obj::AdapterDetector, arg0::SpacecraftState)
    return jcall(obj, "resetState", SpacecraftState, (SpacecraftState,), arg0)
end


text:
# Functions for plotting molecular dynamics simulation results

PLOTTING_TIME_UNIT = u"ps"
PLOTTING_TEMP_UNIT = u"K"
PLOTTING_ENERGY_UNIT = u"hartree"

# Extract (timestep, time) tuples at stride intervals from a simulation result
function plotting_time_range(result::MolecularDynamicsResult, stride::Integer)
    [(t, uconvert(PLOTTING_TIME_UNIT, time)) for (t, time) ∈ enumerate(get_time_range(result)) if (t - 1) % stride == 0]
end

"""
    plot_temperature(result::MolecularDynamicsResult, stride::Integer)::Plot

Plot the temperature of a `MolecularDynamicsResult` against time, sampling every `stride` points.
"""
function plot_temperature(result::MolecularDynamicsResult, stride::Integer)
    N = length(get_system(result))
    p = plot(
        title = "Temperature during Simulation [n = $(N)]",
        xlab = "Time",
        ylab = "Temperature",
    )
    plot_temperature!(p, result, stride, true)
end
"""
    plot_temperature!(p::Plot, result::MolecularDynamicsResult, stride::Integer, first_plot::Bool = false)::Plot

Plot the temperature of a `MolecularDynamicsResult` against time, sampling every `stride` points.
Adds the new line to an existing plot and update the legend only if it is the first plot.
If it is not the first plot, adds a vertical line to differentiate the segments of the simulation in the plot.
"""
function plot_temperature!(p::Plot, result::MolecularDynamicsResult, stride::Integer, first_plot::Bool = false)
    time_range = plotting_time_range(result, stride)
    if (!first_plot)
        vline!(
            p,
            [time_range[1][2]],
            label = false,
            color = :black,
            linestyle = :dot,
            lw = 2
        )
    end
    plot!(
        p,
        [time for (t, time) ∈ time_range],
        [uconvert(PLOTTING_TEMP_UNIT, temperature(result, t)) for (t, time) ∈ time_range],
        label = first_plot ? "Simulation Temperature" : nothing,
        color = 1,
    )
    reference_temp = uconvert(PLOTTING_TEMP_UNIT, reference_temperature(result))
    if (!ismissing(reference_temp))
        plot!(
            p,
            [time for (t, time) ∈ time_range],
            [reference_temp for (t, time) ∈ time_range],
            label = first_plot ? "Reference Temperature" : nothing,
            color = 2,
            linestyle = :dash,
            lw = 2
        )
    end
    p
end

"""
    plot_energy(result::MolecularDynamicsResult, stride::Integer)::Plot

Plot the kinetic, potential, and total energy of a `MolecularDynamicsResult` against time, sampling every `stride` points.
"""
function plot_energy(result::MolecularDynamicsResult, stride::Integer)
    N = length(get_system(result))
    p = plot(
        title = "Energy during Simulation [n = $(N)]",
        xlab = "Time",
        ylab = "Energy",
        legend = :right
    )
    plot_energy!(p, result, stride, true)
end
"""
    plot_energy!(p::Plot, result::MolecularDynamicsResult, stride::Integer, first_plot::Bool = false)::Plot

Plot the kinetic, potential, and total energy of a `MolecularDynamicsResult` against time, sampling every `stride` points.
Add the new lines to an existing plot and update the legend only if it is the first plot.
If it is not the first plot, add a vertical line to differentiate the segments of the simulation in the plot.
"""
function plot_energy!(p::Plot, result::MolecularDynamicsResult, stride::Integer, first_plot::Bool = false)
    time_range = plotting_time_range(result, stride)
    if (!first_plot)
        vline!(p, [time_range[1][2]], label = false, color = :black, linestyle = :dot, lw = 2)
    end
    plot!(
        p,
        [time for (t, time) ∈ time_range],
        [uconvert(PLOTTING_ENERGY_UNIT, kinetic_energy(result, t)) for (t, time) ∈ time_range],
        label = first_plot ? "Kinetic Energy" : nothing,
        color = 2
    )
    plot!(
        p,
        [time for (t, time) ∈ time_range],
        [uconvert(PLOTTING_ENERGY_UNIT, potential_energy(result, t)) for (t, time) ∈ time_range],
        label = first_plot ? "Potential Energy" : nothing,
        color = 1
    )
    plot!(
        p,
        [time for (t, time) ∈ time_range],
        [uconvert(PLOTTING_ENERGY_UNIT, total_energy(result, t)) for (t, time) ∈ time_range],
        label = first_plot ? "Total Energy" : nothing,
        color = 3
    )
end

"""
    plot_rdf(result::MolecularDynamicsResult, σ::Real, sample_fraction::Float64 = 1.0)::Plot

Plot the radial distribution function of a `MolecularDynamicsResult` sampling only the trailing `sample_fraction` of the timesteps.
Use σ (from Lennard Jones) as a normalization factor for the radius (assumed to be in atomic units).
"""
function plot_rdf(result::MolecularDynamicsResult, σ::Real, sample_fraction::Float64 = 1.0)
    plot_rdf(result, σ * LENGTH_UNIT, sample_fraction)
end
"""
    plot_rdf(result::MolecularDynamicsResult, σ::Unitful.Length, sample_fraction::Float64 = 1.0)::Plot

Plot the radial distribution function of a `MolecularDynamicsResult` sampling only the trailing `sample_fraction` of the timesteps.
Use σ (from Lennard Jones) as a normalization factor for the radius.
"""
function plot_rdf(result::MolecularDynamicsResult, σ::Unitful.Length, sample_fraction::Float64 = 1.0)
    @assert 0 < sample_fraction ≤ 1
    N = length(get_system(result))
    T = length(get_time_range(result)) - 1
    rs, grf = rdf(result, sample_fraction)
    @assert length(rs) == length(grf)
    plot(
        title = "Radial Distribution Function [n = $(N)] [T = $(T)]",
        xlab = "Distance r/σ",
        ylab = "Radial Distribution g(r)",
        legend = false
    )
    plot!(rs / austrip(σ), grf)
end


text:
using CSV, CategoricalArrays, Plots; pyplot()

df = CSV.read("../data/companyData.csv")
companies = levels(df.Type)

year2012 = df[df.Year .== 2012, :MarketCap]
year2016 = df[df.Year .== 2016, :MarketCap]

p1 = pie(companies, year2012, title="2012 Market Cap \n by company")
p2 = pie(companies, year2016, title="2016 Market Cap \n by company")
plot(p1, p2, size=(800, 400))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using StatsBase

"""
    AbstractPolicy
"""
abstract type AbstractPolicy end

function Base.get(π::AbstractPolicy, state_t, action_t)
    throw(DomainError("get(PolicyType, args...) not defined!"))
end

"""
    NullPolicy <: AbstractPolicy

Always returns 1.0
"""
struct NullPolicy <: AbstractPolicy
end
Base.get(π::NullPolicy, state_t, action_t) = 1.0

"""
    PersistentPolicy <: AbstractPolicy

Returns 1 of `action_t == action`
"""
struct PersistentPolicy <: AbstractPolicy
    action::Int64
end

Base.get(π::PersistentPolicy, state_t, action_t) =
    π.action == action_t ? 1 : 0

"""
    RandomPolicy{T<:AbstractFloat} <: AbstractPolicy

"""
struct RandomPolicy{T<:AbstractFloat} <: AbstractPolicy
    probabilities::Array{T,1}
    weight_vec::Weights{T, T, Array{T, 1}}
    RandomPolicy(probabilities::Array{T,1}) where {T<:AbstractFloat} = new{T}(probabilities, Weights(probabilities))
end

Base.get(π::RandomPolicy, state_t, action_t) =
    π.probabilities[action_t]

StatsBase.sample(π::RandomPolicy) = StatsBase.sample(Random.GLOBAL_RNG, π)
StatsBase.sample(rng::Random.AbstractRNG, π::RandomPolicy) = StatsBase.sample(rng, π.weight_vec)
StatsBase.sample(rng::Random.AbstractRNG, π::RandomPolicy, state) = StatsBase.sample(rng, π.weight_vec)

"""
    FunctionalPolicy{F <: Function} <: AbstractPolicy

"""
struct FunctionalPolicy{F <: Function} <: AbstractPolicy
    func::F
end

Base.get(π::FunctionalPolicy, state_t, action_t) =
    π.func(state_t, action_t)

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
require("../../src/BlackBoxOptim.jl")
require("cec_cmsa_es.jl")
require("../experiment_framework.jl")
require("../../src/experiments/parameter_experiment.jl")

using BlackBoxOptim

for(dim in [32])

  n = dim
  p = BlackBoxOptim.as_fixed_dim_problem(BlackBoxOptim.example_problems["Rastrigin"], n)
  p = BlackBoxOptim.shifted(p)
  diameter = minimum(diameters(search_space(p)))

  pe = ParameterExperiment(
  ["lambda",
   "mu",
   "sigma",
   "decompose_covar_prob",
   "tau",
   "tau_c"
  ],
  [((ds, ps) -> int(10^ds[1])),
   ((ds, ps) -> int(max(1.0, ds[2]*ps[1]))),
   ((ds, ps) -> diameter * 10^ds[3]),
   ((ds, ps) -> ds[4]),
   ((ds, ps) -> 1 / sqrt(n * ds[5])),
   ((ds, ps) -> 1 + n * (n+1) / (ds[6]*ps[2]))
  ],
  [( log10(4.0 + log(n)), log10(10*n*n)),
   ( 1/100, 1/2),
   (-3.0, 0.0),
   (0.01, 0.99),
   (0.2, 10.0),
   (0.01, 100.0)
  ]
 )

  for num_repeats in [2]
    prefix = "exp1_6_cmsa_es_6param_log10_$(name(p))$(n)_$(num_repeats)"
    explore_parameters(pe, cmsa_es, p; experiment_prefix = prefix, num_repeats = num_repeats)
  end

end


text:
using ConstructionBaseExtras, StaticArrays, Test, ConstructionBase

@testset "constructorof" begin
    sa = @SVector [2, 4, 6, 8]
    sa2 = ConstructionBase.constructorof(typeof(sa))((3.0, 5.0, 7.0, 9.0))
    @test sa2 === @SVector [3.0, 5.0, 7.0, 9.0]

    ma = @MMatrix [2.0 4.0; 6.0 8.0]
    ma2 = ConstructionBase.constructorof(typeof(ma))((1, 2, 3, 4))
    @test ma2 isa MArray{Tuple{2,2},Int,2,4}
    @test all(ma2 .=== @MMatrix [1 3; 2 4])

    sz = SizedArray{Tuple{2,2}}([1 2;3 4])
    sz2 = ConstructionBase.constructorof(typeof(sz))([:a :b; :c :d])
    @test sz2 == SizedArray{Tuple{2,2}}([:a :b; :c :d])
    @test typeof(sz2) <: SizedArray{Tuple{2,2},Symbol,2,2}
end


text:
struct TruncatedPlotContext{T₁ <: AbstractControl, T₂ <: PlotContext,   T₃ <: NestedIntervalContext} <: AbstractPlotContext
    control::T₁
    plot_context::T₂
    interval_context::T₃
end

function (context::TruncatedPlotContext{<: PlotlinesControl, <:  PlotContext,   <: NestedIntervalContext})()
        control = context.control
        plot_context = context.plot_context
        interval_context = context.interval_context
        #CImGui.Begin("Main",C_NULL, CImGui.ImGuiWindowFlags_NoBringToFrontOnFocus)
        isenabled(control) ? control(plot_context, interval_context) : nothing
        #CImGui.End()
end

function (control::PlotlinesControl)(plot_context::PlotContext, interval_context::NestedIntervalContext)
    plot_control = plot_context.control
    plot_model = plot_context.model
    plot_display_properties = plot_context.display_properties

    interval_control = interval_context.control
    interval_model = interval_context.model
    interval_display_properties = interval_context.display_properties

    data = get_data(plot_model)
    i₀ = Int(get_start(interval_model))
    i₁ = Int(get_stop(interval_model))
    plot_model₂ = PlotlinesModel(data[i₀:i₁])

    #CImGui.Begin("EDA",C_NULL, CImGui.ImGuiWindowFlags_NoBringToFrontOnFocus)
    isenabled(plot_control) ? plot_control(plot_model₂, plot_display_properties) : nothing
    CImGui.Dummy(ImVec2(0, 20))
    isenabled(interval_control) ? interval_control(interval_model, interval_display_properties) : nothing
    #CImGui.End()
end


text:
using VectorAutoregressions, Base.Test

# comparison against http://www-personal.umich.edu/~lkilian/figure9_1_chol.zip
y = readdlm("cholvar_test_data.csv",',')

V = VAR(y, 4, true)

@test V.p == 4
@test V.Σ ≈ [
        312.524637669663 0.773643992117187 0.919341486304653
        0.773643992117187 0.0515304724687740 0.0148802276681216
        0.919341486304653 0.0148802276681216 0.557021891337960
        ]

@test V.β[:,1:4] ≈ [
        -0.680676727884458 -0.00637960178073601 0.936519611007488 1.75104737240770
        -0.0152242907833718 0.00177981748889471 0.599086054531876 0.0281955980964510
        0.436329751995733 -0.00236721376013309 -0.218737236527547 0.332361646200618
        ]

@test V.ϵ[:,1:4] ≈ [
        70.3934221013155 -0.443580106581005 4.21164200500382 1.22299717886091
        0.189207494365181 0.287694213161629 0.692414550347607 0.405418261036420
        -1.39379479232130 0.277340573620367 -0.573151724739961 -0.395501452903569
        ]

mIRFa = IRFs_a(V,4,true)

@test [mIRFa.IRF[1,1:4]'; mIRFa.IRF[4,1:4]'; mIRFa.IRF[7,1:4]'] ≈ [
        17.6783663744607 0.0192642533976305 -2.85156060358692 0.594909295963393
        0.0437621879606954 0.0591478592307352 0.113331057827868 0.109932348126692
        0.0520037579735191 -0.0341368375885218 -0.111284843056732 -0.0491881671489620
        ]

mIRFb = IRFs_b(V,4,10,true)

@test [mIRFb.IRF[1,1:4]'; mIRFb.IRF[4,1:4]'; mIRFb.IRF[7,1:4]'] ≈ [
        17.6783663744607 0.0192642533976305 -2.85156060358692 0.594909295963393
        0.0437621879606954 0.0591478592307352 0.113331057827868 0.109932348126692
        0.0520037579735191 -0.0341368375885218 -0.111284843056732 -0.0491881671489620
        ]


text:
#
# Advent of Code
# Day 12: Leonardo's Monorail
#

println("Advent of Code")
println("Day 12: Leonardo's Monorail")
println("")

DEBUG = false

# Assembunny
#
# cpy x y copies x (either an integer or the value of a register) into register y.
# inc x increases the value of register x by one.
# dec x decreases the value of register x by one.
# jnz x y jumps to an instruction y away (positive means forward; negative means backward), but only if x is not zero.
#

# Data file
# Input data has been saved to this file
inputDataFileName = "day12-input.txt"

# Read the input data, parse it down to a list of "commands"
inputDataFile = open(inputDataFileName, "r")
inputData = readstring(inputDataFile)
close(inputDataFile)

# inputData = """cpy 41 a
# inc a
# inc a
# dec a
# jnz a 2
# dec a"""

inputDataLines = split(strip(inputData), '\n')

# part 1
#registers = Dict([("a",0),("b",0),("c",0),("d",0)])

# part 2
registers = Dict([("a",0),("b",0),("c",1),("d",0)])

lineCounter = 1

while lineCounter <= length(inputDataLines)
  value = 0 # scope rules are weird
  line = inputDataLines[lineCounter]
  m = match(r"(.{3}) ([\w\d]+)(?:\s?)([-\w]*)", line)
  if m != nothing
    if m.captures[1] == "cpy"
      try
        value = registers[m.captures[2]]
        #if DEBUG println("cpy value: ", value) end
      catch
        value = parse(m.captures[2])
        #if DEBUG println("cpy from register ", m.captures[2], ": ", registers[m.captures[2]]) end
      end
      registers[m.captures[3]] = value
      if DEBUG println("cpy ", value, " to ", m.captures[3]) end
      lineCounter += 1
    elseif m.captures[1] == "inc"
      registers[m.captures[2]] += 1
      if DEBUG println("inc ", m.captures[2], ": ", registers[m.captures[2]]) end
      lineCounter += 1
    elseif m.captures[1] == "dec"
      registers[m.captures[2]] -= 1
      if DEBUG println("dec ", m.captures[2], ": ", registers[m.captures[2]]) end
      lineCounter += 1
    elseif m.captures[1] == "jnz"
      try
        value = registers[m.captures[2]]
      catch
        value = parse(m.captures[2])
      end
      if DEBUG println("jnz value: ", value) end
      if value != 0
        lineCounter += parse(m.captures[3])
        if DEBUG println("jnz to ", lineCounter) end
      else
        if DEBUG println("jnz 0") end
        lineCounter += 1
      end
    end
  else
    if DEBUG println("Unable to parse: ", line) end
  end

  if DEBUG println("line: ", lineCounter) end
end

println("Value in register a: ", registers["a"])


text:
immutable Beta <: ContinuousUnivariateDistribution
    alpha::Float64
    beta::Float64
    function Beta(a::Real, b::Real)
        (a > zero(a) && b > zero(b)) || error("alpha and beta must be positive")
        new(float64(a), float64(b))
    end
end

Beta(a::Real) = Beta(a, a) # symmetric in [0, 1]
Beta() = Beta(1.0) # uniform

@_jl_dist_2p Beta beta

function entropy(d::Beta)
    o = lbeta(d.alpha, d.beta)
    o -= (d.alpha - 1.0) * digamma(d.alpha)
    o -= (d.beta - 1.0) * digamma(d.beta)
    o += (d.alpha + d.beta - 2.0) * digamma(d.alpha + d.beta)
    o
end

function kurtosis(d::Beta)
    α, β = d.alpha, d.beta
    num = 6.0 * ((α - β)^2 * (α + β + 1.0) - α * β * (α + β + 2.0))
    den = α * β * (α + β + 2.0) * (α + β + 3.0)
    num / den
end

mean(d::Beta) = d.alpha / (d.alpha + d.beta)

median(d::Beta) = quantile(d, 0.5)

function mode(d::Beta)
    α, β = d.alpha, d.beta
    if α >= 1.0
        if β > 1.0
            (α - 1.0) / (α + β - 2.0)
        elseif α == 1.0 && β == 1.0
            # Uniform[0,1]: what should be returned?
            0.5
        else
            1.0
        end
    else
        if β >= 1.0
            0.0
        else
            # not unique: return largest
            α > β ? 1.0 : 0.0
        end
    end
end

function modes(d::Beta)
    α, β = d.alpha, d.beta
    if α >= 1.0
        if β > 1.0
            [(α - 1.0) / (α + β - 2.0)]
        elseif α == 1.0 && β == 1.0
            # Uniform[0,1]: what should be returned?
            Float64[]
        else
            [1.0]
        end
    else
        if β >= 1.0
            [0.0]
        else
            [0.0,1.0]
        end
    end
end

function skewness(d::Beta)
    num = 2.0 * (d.beta - d.alpha) * sqrt(d.alpha + d.beta + 1.0)
    den = (d.alpha + d.beta + 2.0) * sqrt(d.alpha * d.beta)
    num / den
end

function var(d::Beta)
    ab = d.alpha + d.beta
    d.alpha * d.beta / (ab * ab * (ab + 1.0))
end

function gradlogpdf(d::Beta, x::Real)
  insupport(Beta, x) ? (d.alpha - 1.0) / x - (d.beta - 1.0) / (1 - x) : 0.0
end

### handling support

@continuous_distr_support Beta 0.0 1.0

## Fit model

# TODO: add MLE method (should be similar to Dirichlet)

# This is a moment-matching method (not MLE)
#
function fit(::Type{Beta}, x::Array)
    for xi in x
        insupport(Beta, xi) || error("Beta observations must be in [0,1]")
    end
    x_bar = mean(x)
    v_bar = varm(x, x_bar)
    α = x_bar * (((x_bar * (1.0 - x_bar)) / v_bar) - 1.0)
    β = (1.0 - x_bar) * (((x_bar * (1.0 - x_bar)) / v_bar) - 1.0)
    Beta(α, β)
end


text:
using Test, DiffEqBase

f123(x,y) = 1
f123(::Val{:jac}, x) = 2

g123(x,y) = 1
g123(::Val{:jac}, x::T) where {T} = 2

h123(x,y) = 1
h123(x::T) where {T} = 2

struct T123 end
(::T123)(a) = 1
(::T123)(::Val{:jac}, a) = 1
t123 = T123()

struct G123 end
(::G123)(a) = 1
(::G123)(b, a) = 1

@test DiffEqBase.__has_jac(f123)
@test DiffEqBase.__has_jac(g123)
@test !DiffEqBase.__has_jac(h123)
@test DiffEqBase.__has_jac(T123())
@test !DiffEqBase.__has_jac(G123())


text:
using SDELab2 # SDELab2.jl should be findable via LOAD_PATH
println("run_HH.jl")
#
# Hodgkin-Huxley example from
# H. Alzubaidi and T. Shardlow. Numerical simulations of SDEs and SPDEs.
# In: Stochastic Methods in Neuroscience. Ed. by C. Laing and G. Lord. OUP,
# 2009. Chap. 12, pp. 344–366.
# doi: 10.1093/acprof:oso/9780199235070.003.0012 .
#
function drift_f(y,Mu)
  g_k=36; g_Na=120; g_L=0.3
  C_m=1 # membrane capabilities
  C1=g_k/C_m; C2=g_Na/C_m; C3=g_L/C_m
  # constants of resting potential of such iterations
  v_k=-12; v_Na=115; v_L=10.613
  a_n=(10-y[1]) / (100*(exp((10-y[1])/10)-1))
  b_n=exp(-y[1]/80)/8
  a_m=(25-y[1])/(10*(exp((25-y[1])/10)-1.))
  b_m=4*exp(-y[1]/18)
  a_h=(7*exp(-y[1]/20))/100
  b_h=1/(exp((30-y[1])/10)+1)
  #
  # compute drift function for HH model
  z1=(C1*(v_k-y[1])*(y[2]^4)) + (C2*(v_Na-y[1])*(y[3]^3)*y[4]) + (C3*(v_L-y[1])) + Mu
  z2=a_n*(1-y[2])-b_n*y[2]
  z3=a_m*(1-y[3])-b_m*y[3]
  z4=a_h*(1-y[4])-b_h*y[4]
  return [z1,z2,z3,z4]
end

#
function diff_g_mat(u,dw,sigma)
  return   [sigma*dw;0;0;0]
end
#
function diff_g_vecs(sig)
  m=1
  A=Array(Function,m)
  g1(u)=[sig,0,0,0]
  A[1]=g1
  return A
end
# dimensions
d=4
m=1
noise_type=3 # commutative
# params
sig=0.1
Mu=2.5
# initial data
a0_n=0.1/(exp(1.)-1.); b0_n=0.125
a0_m=2.5/(exp(2.5)-1); b0_m=4;
a0_h=0.07; b0_h=1/(exp(3.)+1)
y0=[0; a0_n/(a0_n+b0_n); a0_m/(a0_m+b0_m); a0_h/(a0_h+b0_h)]
#
t0=0
tf=40

#
fcn=set_fcn(d,m,x->drift_f(x,Mu),(x,dw)->diff_g_mat(x,dw,sig),diff_g_vecs(sig), noise_type)
opt=set_opt(0.01,"SIBDF")
opt["do_plotting"]=true
opt["output_plot_type"]="path_plot"
opt["output_sel"]=[1,2]
#using PyPlot
println("Start sde_strong_soln")
tic();t,y,W=sde_strong_solution(fcn, linspace(t0,tf,100), y0, opt);run_time=toc()

#


text:
module EncodeTest

using Test
using ErlPort.Exceptions
using ErlPort.ErlTerms.Encode

function run()::Nothing
    # tests for supporting functions
    @testset begin
        @test charint4pack(1) == b"\x00\x00\x00\x01"
        @test charint4pack(1000) == b"\x00\x00\x03\xe8"
        @test charint4pack(72000) == b"\x00\x01\x19@"
        @test charint4pack(1000000000) == b";\x9a\xca\x00"
        @test charint4pack(4294967295) == b"\xff\xff\xff\xff"
        @test charint2pack(1) == b"\x00\x01"
        @test charint2pack(1000) == b"\x03\xe8"
        @test charint2pack(10000) == b"'\x10"
        @test charint2pack(50092) == b"\xc3\xac"
        @test charint2pack(65535) == b"\xff\xff"
    end

    # tests for atoms
    @testset begin
        @test encode(:erlang) == b"\x83d\0\6erlang"
        @test encode(:test) == b"\x83d\0\4test"
        @test encode(Symbol("")) == b"\x83d\0\0"
        @test encodeterm(:erlang) == b"d\0\6erlang"
        @test encodeterm(:test) == b"d\0\4test"
        @test encodeterm(Symbol("")) == b"d\0\0"
    end

    # tests for bools
    @testset begin
        @test encode(true) == b"\x83d\0\4true"
        @test encode(:true) == b"\x83d\0\4true"
        @test encode(false) == b"\x83d\0\5false"
        @test encode(:false) == b"\x83d\0\5false"
        @test encodeterm(true) == b"d\0\4true"
        @test encodeterm(:true) == b"d\0\4true"
        @test encodeterm(false) == b"d\0\5false"
        @test encodeterm(:false) == b"d\0\5false"
    end

    # tests for nothing
    @testset begin
        @test encode(nothing) == b"\x83d\0\x09undefined"
        @test encode(:nothing) == b"\x83d\0\7nothing"
        @test encodeterm(nothing) == b"d\0\x09undefined"
        @test encodeterm(:nothing) == b"d\0\7nothing"
    end

    # tests for list
    @testset begin
        @test encode([]) == b"\x83l\0\0\0\0j"
        @test encode([1]) == b"\x83l\0\0\0\1a\1j"
        @test encode([1, 2]) == b"\x83l\0\0\0\2a\1a\2j"
        @test encode([:a]) == b"\x83l\0\0\0\1d\0\1aj"
        @test encode([:a, :ok]) == b"\x83l\0\0\0\2d\0\1ad\0\2okj"
        @test encodeterm([]) == b"l\0\0\0\0j"
        @test encodeterm([1]) == b"l\0\0\0\1a\1j"
        @test encodeterm([1, 2]) == b"l\0\0\0\2a\1a\2j"
        @test encodeterm([:a]) == b"l\0\0\0\1d\0\1aj"
        @test encodeterm([:a, :ok]) == b"l\0\0\0\2d\0\1ad\0\2okj"
    end

    # tests for binary data (vector of UInt8's)
    @testset begin
        @test encode(b"") == b"\x83m\0\0\0\0"
        @test encode(b"something") == b"\x83m\0\0\0\x09something"
        @test encode(b"nothing") == b"\x83m\0\0\0\7nothing"
        @test encodeterm(b"") == b"m\0\0\0\0"
        @test encodeterm(b"something") == b"m\0\0\0\x09something"
        @test encodeterm(b"nothing") == b"m\0\0\0\7nothing"
    end

    # tests for string
    @testset begin
        @test encode("") == b"\x83m\0\0\0\0"
        @test encode("something") == b"\x83m\0\0\0\x09something"
        @test encode("nothing") == b"\x83m\0\0\0\7nothing"
        @test encodeterm("") == b"m\0\0\0\0"
        @test encodeterm("something") == b"m\0\0\0\x09something"
        @test encodeterm("nothing") == b"m\0\0\0\7nothing"

        @test encode(SubString("")) == b"\x83m\0\0\0\0"
        @test encode(SubString("something")) == b"\x83m\0\0\0\x09something"
        @test encode(SubString("nothing")) == b"\x83m\0\0\0\7nothing"
        @test encodeterm(SubString("")) == b"m\0\0\0\0"
        @test encodeterm(SubString("something")) == b"m\0\0\0\x09something"
        @test encodeterm(SubString("nothing")) == b"m\0\0\0\7nothing"
    end

    # tests for tuples
    @testset begin
        @test encode(()) == b"\x83h\0"
        @test encode((1,)) == b"\x83h\1a\1"
        @test encode((1, 2)) == b"\x83h\2a\1a\2"
        @test encodeterm(()) == b"h\0"
        @test encodeterm((1,)) == b"h\1a\1"
        @test encodeterm((1, 2)) == b"h\2a\1a\2"
    end

    # tests for improper list
    @testset begin
    end

    # tests for short int
    @testset begin
        @test encode(0) == b"\x83a\0"
        @test encode(1) == b"\x83a\1"
        @test encode(255) == b"\x83a\xff"
    end

    # tests for int
    @testset begin
        @test encode(-1) == b"\x83b\xff\xff\xff\xff"
        @test encode(256) == b"\x83b\0\0\1\0"
        @test encode(-2147483648) == b"\x83b\x80\0\0\0"
        @test encode(2147483647) == b"\x83b\x7f\xff\xff\xff"
    end

    # tests for long int
    @testset begin
        @test encode(-2147483649) == b"\x83n\4\1\1\0\0\x80"
        @test encode(2147483648) == b"\x83n\4\0\0\0\0\x80"
        @test encode(BigInt(2)^2040) == vcat(b"\x83o\0\0\1\0\0", zeros(255), b"\1")
        @test encode(-BigInt(2)^2040) == vcat(b"\x83o\0\0\1\0\1", zeros(255), b"\1")
    end

    # tests for int errors
    # XXX add some!

    # tests for float
    @testset begin
        @test encode(0.0) == b"\x83\x46\0\0\0\0\0\0\0\0"
        @test encode(1.5) == b"\x83\x46\x3f\xf8\0\0\0\0\0\0"
        @test encode(NaN) == b"\x83d\0\3nan"
        @test encodeterm(NaN) == b"d\0\3nan"
    end

    # tests for opaque objects
    @testset begin
    end

    # tests for compressed terms
    @testset begin
    end

    # tests for map
    @testset "encode map" begin
        # \x61 = 100 = SMALL_INT_EXT
        # \x64 = 100 = ATOM_EXT
        # \x6a = 106 = NIL_EXT
        # \x6c = 108 = LIST_EXT
        # \x6d = 109 = BINARY_EXT
        # \x74 = 116 = MAP_EXT
        @test encode(Dict()) == b"\x83\x74\0\0\0\0"
        @test encode(Dict(:a => 2)) == b"\x83\x74\0\0\0\1\x64\0\1a\x61\2"
        @test encode(Dict(:a => Dict(:a => 2))) ==
              b"\x83\x74\0\0\0\1\x64\0\1a\x74\0\0\0\1\x64\0\1a\x61\2"

        # In the BERT version, either element of the map can come first.
        @test encode(Dict(b"b" => [], :a => 2)) in
              (b"\x83\x74\0\0\0\2\x6d\0\0\0\1b\x6c\0\0\0\0\x6a\x64\0\1a\x61\2",
               b"\x83\x74\0\0\0\2\x64\0\1a\x61\2\x6d\0\0\0\1b\x6c\0\0\0\0\x6a")

        @test encodeterm(Dict()) == b"\x74\0\0\0\0"
        @test encodeterm(Dict(:a => 2)) == b"\x74\0\0\0\1d\0\1a\x61\2"
        @test encodeterm(Dict(:a => Dict(:a => 2))) ==
              b"\x74\0\0\0\1\x64\0\1a\x74\0\0\0\1\x64\0\1a\x61\2"

        @test encodeterm(Dict(b"b" => [], :a => 2)) in
              (b"\x74\0\0\0\2\x6d\0\0\0\1b\x6c\0\0\0\0\x6a\x64\0\1a\x61\2",
               b"\x74\0\0\0\2\x64\0\1a\x61\2\x6d\0\0\0\1b\x6c\0\0\0\0\x6a")
    end

    return nothing
end

end # module


text:
using IrrationalExpressions
using Base.Test

@test BigFloat(-π) == -BigFloat(π)
@test BigFloat(-(-π)) == BigFloat(π)
@test BigFloat(2π+1//3) == BigFloat(2)*BigFloat(π)+BigFloat(1//3)

@test isa(π, Irrational)
@test isa(π+1, IrrationalExpressions.IrrationalExpr)
@test isa(π+1.0, Float64)
@test isa(2π+1.0, Float64)
@test isa(π+BigFloat(1), BigFloat)
@test +(1+π) == 1+π
@test string(1+pi) == "1 + π = 4.1415926535897..."
@test string(2π+0.0) == "6.283185307179586"
@test string(2π+BigFloat(0)) == "6.283185307179586476925286766559005768394338798750211641949889184615632812572396"


text:
using DemoPackageXL
using Test

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


text:
using DSGE, Test, ModelConstructors

m = AnSchorfheide()
if haskey(ENV, "FRED_API_KEY") || isfile(joinpath(ENV["HOME"],".freddatarc"))
    load_data(m)

    @testset "Test util functions" begin
        @test get_class(:histobs) == :obs
        @test get_class(:histpseudo) == :pseudo
        @test get_class(:histstates) == :states
        @test get_class(:histstdshock) == :stdshocks
        @test get_class(:histshocks) == :shocks
        @test_throws ErrorException get_class(:bad)

        @test get_product(:bddhistforecast4qobs) == :bddhistforecast4q
        @test get_product(:dettrendobs) == :dettrend
        @test_throws ErrorException get_product(:bad)

        @test DSGE.get_class_longname(:pseudo)==:pseudoobservable
        @test DSGE.get_class_longname(:obs)==:observable
        @test DSGE.get_class_longname(:states)==:state
        @test DSGE.get_class_longname(:shocks)==:shock
        @test DSGE.parse_transform(Symbol("DSGE.loggrowthtopct_annualized"))([1.0]) ==
            DSGE.loggrowthtopct_annualized([1.0])

        m <= Setting(:use_population_forecast, false)
        m <= Setting(:hpfilter_population, true)
        DSGE.load_population_growth(m)
        m <= Setting(:hpfilter_population, false)
        DSGE.load_population_growth(m)
        m <= Setting(:use_population_forecast, true)
        @test_throws ErrorException DSGE.load_population_growth(m)

        @test DSGE.get_y0_index(m, :forecast) == 241
        @test DSGE.get_y0_index(m, :forecast4q) == 238
        @test DSGE.get_y0_index(m, :shockdec) == 2
        @test DSGE.get_y0_index(m, :hist) == 2
        @test DSGE.get_y0_index(m, :irf) == -1
        @test_throws ErrorException DSGE.get_y0_index(m, :bad)

        @test DSGE.check_consistent_order([1; 2; 3], [4; 5; 6])
        @test_throws AssertionError DSGE.check_consistent_order([1; 2; 3], [5; 4; 6])
        #would be nice to test:
        # resize population_forecast
        # get_population_series
        # get_mb_populuation_series

    end
else
        @warn "Skipping fred_data test because FRED_API_KEY not present"
end

@testset "Test prior_table works" begin
    DSGE.prior_table(m)

    DSGE.posterior_table(m, ones(16), ones(16, 2))

    DSGE.prior_posterior_moments_table(m, ones(16), ones(16, 2))
    DSGE.prior_posterior_table(m, ones(16))

end

@testset "Test finding density bands" begin
    @test find_density_bands(ones(100, 100), 0.95) == ones(2, 100)
    @test find_density_bands(ones(1, 100), 0.95) == ones(2, 100)
    @test find_density_bands(ones(100, 100), 0.95, minimize = true) == ones(2, 100)
    dens = find_density_bands(ones(100, 100), [.8, .9])
    @test names(dens) == [Symbol("80.0% UB"), Symbol("80.0% LB"), Symbol("90.0% UB"), Symbol("90.0% LB")]
    @test Matrix(dens) == ones(100,4)
end


text:
using LinearAlgebra

function err_estimation(xfull, xhalf, order::Integer)
    return (1<<order)/(1<<order-1)*norm(xhalf .- xfull)
end

function richardson_extrapolation(xfull, xhalf, order::Integer)
    return ((1<<order)*xhalf .- xfull)/(1<<order-1)
end

function mixed_step_ode_solver(
    f,
    t0::Number,
    x0,
    stop_cond;
    method = mod_euler_step,
    min_step::AbstractFloat = 1e-3,
    max_step::AbstractFloat = 10,
    target_err::AbstractFloat = 0.1,
    step0::AbstractFloat = 1.0,
    maxit::Integer = 10000,
    debug_dict = nothing
)
    t = Array{typeof(t0 + step0), 1}()
    x = Array{typeof(x0 + step0*f(t0, x0)), 1}()
    step = step0
    push!(t, t0)
    push!(x, x0)
    i = 1
    for j = 1:maxit
        aux = last(x) .+ method(f, last(t), last(x), [step, step/2])
        xfull = aux[:,1]
        xhalf = aux[:,2] + method(f, last(t)+step/2, aux[:,2], step/2)

        err_full = err_estimation(xfull, xhalf, method(:order))
        if err_full < target_err
            push!(t, last(t)+step)
            push!(x, richardson_extrapolation(xfull, xhalf, method(:order)))
            if stop_cond(t, x)
                if debug_dict != nothing
                    debug_dict[:numit] = j
                end
                break
            end
        end

        q = (target_err*step/err_full) ^ (1/(method(:order)))
        q = min(4, max(0.1, q))
        step = min(q*step, max_step)
        if step < min_step
            throw(ErrorException("step went below the min_step($(min_step)) value"))
        end
    end
    return t, x
end

struct ExtButcherTableau
    a
    c
    b1
    b2
    order

    ExtButcherTableau(; a, c, b1, b2, order::Integer) = new(a, c, b1, b2, order)
end

DORPRI5_tableau = ExtButcherTableau(
    a = [          0  1/5        3/40    44/45     19372/6561   9017/3168     35/384;
                   0    0        9/40   -56/15    -25360/2187     -355/33          0;
                   0    0           0     32/9     64448/6561  46732/5247   500/1113;
                   0    0           0        0       -212/729      49/176    125/192;
                   0    0           0        0              0 -5103/18656 -2187/6784;
                   0    0           0        0              0           0      11/84;
                   0    0           0        0              0           0          0],
    c = [          0, 1/5,       3/10,     4/5,           8/9,          1,         1],

    b1 = [    35/384,   0,   500/1113, 125/192,    -2187/6784,      11/84,         0],
    b2 = [5179/57600,   0, 7571/16695, 393/640, -92097/339200,   187/2100,      1/40],
    order = 5
)

RKFehlberg_tableau = ExtButcherTableau(
    a = [       0  1/4        3/32    1932/2197   439/216      -8/27;
                0    0        9/32   -7200/2197        -8          2;
                0    0           0    7296/2197  3680/513 -3544/2565;
                0    0           0            0 -845/4104  1859/4104;
                0    0           0            0         0     -11/40;
                0    0           0            0         0          0],
    c = [       0, 1/4,        3/8,       12/13,        1,       1/2],

    b1 = [ 16/135,   0, 6656/12825, 28561/56430,    -9/50,      2/55],
    b2 = [ 25/216,   0,  1408/2565,   2197/4104,     -1/5,         0],
    order = 4
)

function double_step(tab::ExtButcherTableau, f, t, x, step)
    k = Array{typeof(x + step*f(t, x)), 1}(undef, length(tab.c))
    k[1] = step*f(t + step*tab.c[1], x)
    for j = 2:length(k)
        k[j] = step*f(t + step*tab.c[j], x + tab.a[1:(j-1), j]' * k[1:(j-1)])
    end
    return (x + tab.b1' * k, x + tab.b2' * k)
end

# function double_step(

function adaptative_rk_ode_solver(
    f,
    t0::Number,
    x0,
    stop_cond;
    tableau::ExtButcherTableau = RKFehlberg_tableau,
    min_step::AbstractFloat = 1e-3,
    max_step::AbstractFloat = 10,
    target_err::AbstractFloat = 0.1,
    step0::AbstractFloat = 1.0,
    maxit::Integer = 10000,
    debug_dict = nothing
)
    t = Array{typeof(t0 + step0), 1}()
    x = Array{typeof(x0 + step0*f(t0, x0)), 1}()
    step = step0
    push!(t, t0)
    push!(x, x0)
    i = 1
    for j = 1:maxit
        x1, x2 = double_step(tableau, f, last(t), last(x), step)

        err = norm(x1 .- x2)
        if err < target_err
            push!(t, last(t)+step)
            push!(x, x1)
            if stop_cond(t, x)
                if debug_dict != nothing
                    debug_dict[:numit] = j
                end
                break
            end
        end

        q = (target_err*step/err) ^ (1/(tableau.order))
        q = min(4, max(0.1, q))
        step = min(q*step, max_step)
        if step < min_step
            throw(ErrorException("step went below the min_step($(min_step)) value"))
        end
    end
    return t, x
end

max_stars_repo_name: useredsa/numeric-differential-equations

max_issues_repo_name: useredsa/numeric-differential-equations

max_forks_repo_name: useredsa/numeric-differential-equations


text:
using MAEBOL
using Base.Test

# write your own tests here
@test 1 == 1

# Testing the coefficients
m = SearchDomain(9)
@assert h_ij(0, 0, 9) == 9.0
@assert h_ij(1, 1, 9) == 4.5
@assert h_ij(2, 2, 9) == 4.5
@assert h_ij(3, 3, 9) == 4.5


text:
export shift_graph

"""
    shift_graph(S::Set{String}, len::Int = 0, trim::Bool = false)
Create a word shift graph. For example, *stranger* and *estrange* are
adjacent.

We use only words of length `len`, unless this is set to `0`
in which case all words are used.

If `trim` is set to `true` then vertices of degree zero are deleted.
"""
function shift_graph(S::Set{String}, len::Int = 0, trim::Bool = false)::SimpleGraph{String}

    G = _bare_graph(S, len)

    PM = Progress(NV(G))
    for w in G.V
        next!(PM)
        tail = w[2:end]
        for c in _A2Z
            ww = tail * c
            if has(G, ww)
                add!(G, w, ww)
            end
        end
    end

    # remove vertices of degree zero (if trim is true)
    if trim
        for w in G.V
            if deg(G, w) == 0
                delete!(G, w)
            end
        end
    end

    if len < 1
        name(G, "Shift Graph on all words")
    else
        name(G, "Shift Graph on $len-letter words")
    end

    return G

end