text:
# Simple test on matrix multiplication.
#
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\vvvvvv/!!!!!!!!!!!!!!!!!!!!
# !CAUTION: THIS TEST MUST OCCUR >BEFORE< BLIS WAS IMPORTED.!
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!/^^^^^^\!!!!!!!!!!!!!!!!!!!!
#
using Test
using Random
using LinearAlgebra
using LinearAlgebra: BLAS
using Statistics
Random.seed!(1234)
rtype(::Type{Complex{T}}) where {T} = T
zrtest(val, atol, label) = begin
    iszr = ≈(val, 0.0, atol=atol)
    if !iszr
        @info "`$label` test failed. Consider adding it to ~/.blis_jlbla_blacklist."
    end
    return iszr
end

@testset "BLAS level-3 LinearAlgebra interface" begin
αr = 1.1
βr = 1.1
αc = 1.1 + 0.2im
βc = 1.1 + 0.3im

χlarge = 500
χsmall = 20

Alarge_base = rand(ComplexF64, χlarge, χlarge)
Blarge_base = rand(ComplexF64, χlarge, χlarge)
Clarge_base = rand(ComplexF64, χlarge, χlarge)

Asmall_base = rand(ComplexF64, χsmall, χsmall)
Bsmall_base = rand(ComplexF64, χsmall, χsmall)
Csmall_base = rand(ComplexF64, χsmall, χsmall)

global Clarge_gemm = [zeros(T, χlarge, χlarge) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Clarge_hemm = [zeros(T, χlarge, χlarge) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Clarge_symm = [zeros(T, χlarge, χlarge) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Clarge_her2k= [zeros(T, χlarge, χlarge) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Clarge_syr2k= [zeros(T, χlarge, χlarge) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Clarge_herk = [zeros(T, χlarge, χlarge) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Clarge_syrk = [zeros(T, χlarge, χlarge) for T=(Float32, Float64, ComplexF32, ComplexF64)]
# global Clarge_trmm = [zeros(T, χlarge, χlarge) for T=(Float32, Float64, ComplexF32, ComplexF64)]
# global Clarge_trsm = [zeros(T, χlarge, χlarge) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Csmall_gemm = [zeros(T, χsmall, χsmall) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Csmall_hemm = [zeros(T, χsmall, χsmall) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Csmall_symm = [zeros(T, χsmall, χsmall) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Csmall_her2k= [zeros(T, χsmall, χsmall) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Csmall_syr2k= [zeros(T, χsmall, χsmall) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Csmall_herk = [zeros(T, χsmall, χsmall) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Csmall_syrk = [zeros(T, χsmall, χsmall) for T=(Float32, Float64, ComplexF32, ComplexF64)]

global Cst_lg_gemm = [zeros(T, χlarge÷2, χlarge÷2) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Cst_lg_hemm = [zeros(T, χlarge÷2, χlarge÷2) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Cst_lg_symm = [zeros(T, χlarge÷2, χlarge÷2) for T=(Float32, Float64, ComplexF32, ComplexF64)]
# global Cst_lg_her2k= [zeros(T, χlarge÷2, χlarge÷2) for T=(Float32, Float64, ComplexF32, ComplexF64)]
# global Cst_lg_syr2k= [zeros(T, χlarge÷2, χlarge÷2) for T=(Float32, Float64, ComplexF32, ComplexF64)]
# global Cst_lg_herk = [zeros(T, χlarge÷2, χlarge÷2) for T=(Float32, Float64, ComplexF32, ComplexF64)]
# global Cst_lg_syrk = [zeros(T, χlarge÷2, χlarge÷2) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Cst_sm_gemm = [zeros(T, χsmall÷2, χsmall÷2) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Cst_sm_hemm = [zeros(T, χsmall÷2, χsmall÷2) for T=(Float32, Float64, ComplexF32, ComplexF64)]
global Cst_sm_symm = [zeros(T, χsmall÷2, χsmall÷2) for T=(Float32, Float64, ComplexF32, ComplexF64)]

for (i, T)=zip(1:4, [Float32, Float64, ComplexF32, ComplexF64])
    Alarge = zeros(T, χlarge, χlarge)
    Blarge = zeros(T, χlarge, χlarge)
    Asmall = zeros(T, χsmall, χsmall)
    Bsmall = zeros(T, χsmall, χsmall)

    local elconv, αu, βu, αR, βR
    local locl_hemm!, locl_her2k!, locl_herk!
    if eltype(Alarge)<:Complex
        elconv = x -> x
        αu = T(αc)
        βu = T(βc)
        αR = rtype(T)(αr)
        βR = rtype(T)(βr)
        locl_hemm!  = BLAS.hemm!
        locl_herk!  = BLAS.herk!
        locl_her2k! = BLAS.her2k!
    else
        elconv = x -> real(x)
        αu = T(αr)
        βu = T(βr)
        αR = αu
        βR = βu
        locl_hemm!  = BLAS.symm!
        locl_herk!  = BLAS.syrk!
        locl_her2k! = BLAS.syr2k!
    end
    Alarge .= elconv.(Alarge_base)
    Blarge .= elconv.(Blarge_base)
    Asmall .= elconv.(Asmall_base)
    Bsmall .= elconv.(Bsmall_base)
    Clarge_gemm[i]  .= elconv.(Clarge_base)
    Clarge_hemm[i]  .= elconv.(Clarge_base)
    Clarge_symm[i]  .= elconv.(Clarge_base)
    Clarge_her2k[i] .= elconv.(Clarge_base)
    Clarge_syr2k[i] .= elconv.(Clarge_base)
    Clarge_herk[i]  .= elconv.(Clarge_base)
    Clarge_syrk[i]  .= elconv.(Clarge_base)
    Csmall_gemm[i]  .= elconv.(Csmall_base)
    Csmall_hemm[i]  .= elconv.(Csmall_base)
    Csmall_symm[i]  .= elconv.(Csmall_base)
    Csmall_her2k[i] .= elconv.(Csmall_base)
    Csmall_syr2k[i] .= elconv.(Csmall_base)
    Csmall_herk[i]  .= elconv.(Csmall_base)
    Csmall_syrk[i]  .= elconv.(Csmall_base)

    # Strided.
    Ast_lg = view(Alarge, 1:2:χlarge, 1:2:χlarge)
    Bst_lg = view(Blarge, 1:2:χlarge, 1:2:χlarge)
    Ast_sm = view(Asmall, 1:2:χsmall, 1:2:χsmall)
    Bst_sm = view(Bsmall, 1:2:χsmall, 1:2:χsmall)

    # Execute: column-major.
    BLAS.gemm!('N', 'N', αu, Alarge, Blarge, βu, Clarge_gemm[i])
    BLAS.gemm!('N', 'N', αu, Asmall, Bsmall, βu, Csmall_gemm[i])
    locl_hemm!('L', 'U', αu, Alarge, Blarge, βu, Clarge_hemm[i])
    locl_hemm!('R', 'U', αu, Asmall, Bsmall, βu, Csmall_hemm[i])
    BLAS.symm!('L', 'L', αu, Alarge, Blarge, βu, Clarge_symm[i])
    BLAS.symm!('L', 'L', αu, Asmall, Bsmall, βu, Csmall_symm[i])
    locl_her2k!('U', 'N', αu, Alarge, Blarge, βR, Clarge_her2k[i])
    locl_her2k!('U', 'N', αu, Asmall, Bsmall, βR, Csmall_her2k[i])
    BLAS.syr2k!('U', 'N', αu, Alarge, Blarge, βu, Clarge_syr2k[i])
    BLAS.syr2k!('U', 'N', αu, Asmall, Bsmall, βu, Csmall_syr2k[i])
    locl_herk!('U', 'N', αR, Alarge, βR, Clarge_herk[i])
    locl_herk!('U', 'N', αR, Asmall, βR, Csmall_herk[i])
    BLAS.syrk!('U', 'N', αu, Alarge, βu, Clarge_syrk[i])
    BLAS.syrk!('U', 'N', αu, Asmall, βu, Csmall_syrk[i])

    # Execute: generic-strided.
    Cst_lg_gemm[i] .= Ast_lg * Bst_lg
    Cst_sm_gemm[i] .= Ast_sm * Bst_sm
    Cst_lg_hemm[i] .= Hermitian(Array(Ast_lg)) * Array(Bst_lg)
    Cst_sm_hemm[i] .= Hermitian(Array(Ast_sm)) * Array(Bst_sm)
    Cst_lg_symm[i] .= Symmetric(Array(Ast_lg)) * Array(Bst_lg)
    Cst_sm_symm[i] .= Symmetric(Array(Ast_sm)) * Array(Bst_sm)
end

# Import BLIS for testing.
using BLIS
if length(BLIS.BLASInterface.blacklist) > 0
    @info "Blacklisted methods: $(BLIS.BLASInterface.blacklist)."
end

for (i, T)=zip(1:4, [Float32, Float64, ComplexF32, ComplexF64])
    Alarge = zeros(T, χlarge, χlarge)
    Blarge = zeros(T, χlarge, χlarge)
    Asmall = zeros(T, χsmall, χsmall)
    Bsmall = zeros(T, χsmall, χsmall)

    local elconv, αu, βu, αR, βR
    local locl_hemm!, locl_her2k!, locl_herk!
    if eltype(Alarge)<:Complex
        elconv = x -> x
        αu = T(αc)
        βu = T(βc)
        αR = rtype(T)(αr)
        βR = rtype(T)(βr)
        locl_hemm!  = BLAS.hemm!
        locl_herk!  = BLAS.herk!
        locl_her2k! = BLAS.her2k!
    else
        elconv = x -> real(x)
        αu = T(αr)
        βu = T(βr)
        αR = αu
        βR = βu
        locl_hemm!  = BLAS.symm!
        locl_herk!  = BLAS.syrk!
        locl_her2k! = BLAS.syr2k!
    end
    Alarge .= elconv.(Alarge_base)
    Blarge .= elconv.(Blarge_base)
    Asmall .= elconv.(Asmall_base)
    Bsmall .= elconv.(Bsmall_base)
    Clarge_gemm_cur  = T.(elconv.(Clarge_base))
    Clarge_hemm_cur  = T.(elconv.(Clarge_base))
    Clarge_symm_cur  = T.(elconv.(Clarge_base))
    Clarge_her2k_cur = T.(elconv.(Clarge_base))
    Clarge_syr2k_cur = T.(elconv.(Clarge_base))
    Clarge_herk_cur  = T.(elconv.(Clarge_base))
    Clarge_syrk_cur  = T.(elconv.(Clarge_base))
    Csmall_gemm_cur  = T.(elconv.(Csmall_base))
    Csmall_hemm_cur  = T.(elconv.(Csmall_base))
    Csmall_symm_cur  = T.(elconv.(Csmall_base))
    Csmall_her2k_cur = T.(elconv.(Csmall_base))
    Csmall_syr2k_cur = T.(elconv.(Csmall_base))
    Csmall_herk_cur  = T.(elconv.(Csmall_base))
    Csmall_syrk_cur  = T.(elconv.(Csmall_base))

    # Strided.
    Ast_lg = view(Alarge, 1:2:χlarge, 1:2:χlarge)
    Bst_lg = view(Blarge, 1:2:χlarge, 1:2:χlarge)
    Ast_sm = view(Asmall, 1:2:χsmall, 1:2:χsmall)
    Bst_sm = view(Bsmall, 1:2:χsmall, 1:2:χsmall)

    # Execute: column-major.
    BLAS.gemm!('N', 'N', αu, Alarge, Blarge, βu, Clarge_gemm_cur)
    BLAS.gemm!('N', 'N', αu, Asmall, Bsmall, βu, Csmall_gemm_cur)
    locl_hemm!('L', 'U', αu, Alarge, Blarge, βu, Clarge_hemm_cur)
    locl_hemm!('R', 'U', αu, Asmall, Bsmall, βu, Csmall_hemm_cur)
    BLAS.symm!('L', 'L', αu, Alarge, Blarge, βu, Clarge_symm_cur)
    BLAS.symm!('L', 'L', αu, Asmall, Bsmall, βu, Csmall_symm_cur)
    locl_her2k!('U', 'N', αu, Alarge, Blarge, βR, Clarge_her2k_cur)
    locl_her2k!('U', 'N', αu, Asmall, Bsmall, βR, Csmall_her2k_cur)
    BLAS.syr2k!('U', 'N', αu, Alarge, Blarge, βu, Clarge_syr2k_cur)
    BLAS.syr2k!('U', 'N', αu, Asmall, Bsmall, βu, Csmall_syr2k_cur)
    locl_herk!('U', 'N', αR, Alarge, βR, Clarge_herk_cur)
    locl_herk!('U', 'N', αR, Asmall, βR, Csmall_herk_cur)
    BLAS.syrk!('U', 'N', αu, Alarge, βu, Clarge_syrk_cur)
    BLAS.syrk!('U', 'N', αu, Asmall, βu, Csmall_syrk_cur)

    # Execute: generic-strided.
    Cst_lg_gemm_cur = Ast_lg * Bst_lg
    Cst_sm_gemm_cur = Ast_sm * Bst_sm
    Cst_lg_hemm_cur = Hermitian(Ast_lg) * Bst_lg
    Cst_sm_hemm_cur = Hermitian(Ast_sm) * Bst_sm
    Cst_lg_symm_cur = Symmetric(Ast_lg) * Bst_lg
    Cst_sm_symm_cur = Symmetric(Ast_sm) * Bst_sm

    # Check.
    @test zrtest(mean(abs.(Clarge_gemm_cur  - Clarge_gemm[i] )), 1e-6*χlarge^1.2, "gemm")
    @test zrtest(mean(abs.(Clarge_hemm_cur  - Clarge_hemm[i] )), 1e-6*χlarge^1.2, "hemm")
    @test zrtest(mean(abs.(Clarge_symm_cur  - Clarge_symm[i] )), 1e-6*χlarge^1.2, "symm")
    @test zrtest(mean(abs.(Clarge_her2k_cur - Clarge_her2k[i])), 1e-6*χlarge^1.2, "her2k")
    @test zrtest(mean(abs.(Clarge_syr2k_cur - Clarge_syr2k[i])), 1e-6*χlarge^1.2, "syr2k")
    @test zrtest(mean(abs.(Clarge_herk_cur  - Clarge_herk[i] )), 1e-6*χlarge^1.2, "herk")
    @test zrtest(mean(abs.(Clarge_syrk_cur  - Clarge_syrk[i] )), 1e-6*χlarge^1.2, "syrk")
    @test zrtest(mean(abs.(Csmall_gemm_cur  - Csmall_gemm[i] )), 1e-6*χsmall^1.2, "gemm")
    @test zrtest(mean(abs.(Csmall_hemm_cur  - Csmall_hemm[i] )), 1e-6*χsmall^1.2, "hemm")
    @test zrtest(mean(abs.(Csmall_symm_cur  - Csmall_symm[i] )), 1e-6*χsmall^1.2, "symm")
    @test zrtest(mean(abs.(Csmall_her2k_cur - Csmall_her2k[i])), 1e-6*χsmall^1.2, "her2k")
    @test zrtest(mean(abs.(Csmall_syr2k_cur - Csmall_syr2k[i])), 1e-6*χsmall^1.2, "syr2k")
    @test zrtest(mean(abs.(Csmall_herk_cur  - Csmall_herk[i] )), 1e-6*χsmall^1.2, "herk")
    @test zrtest(mean(abs.(Csmall_syrk_cur  - Csmall_syrk[i] )), 1e-6*χsmall^1.2, "syrk")

    # Check - strided.
    @test zrtest(mean(abs.(Cst_lg_gemm_cur - Cst_lg_gemm[i])), 1e-6*χlarge^1.2, "gemm")
    @test zrtest(mean(abs.(Cst_sm_gemm_cur - Cst_sm_gemm[i])), 1e-6*χsmall^1.2, "gemm")
    @test zrtest(mean(abs.(Cst_lg_hemm_cur - Cst_lg_hemm[i])), 1e-6*χlarge^1.2, "hemm")
    @test zrtest(mean(abs.(Cst_sm_hemm_cur - Cst_sm_hemm[i])), 1e-6*χsmall^1.2, "hemm")
    @test zrtest(mean(abs.(Cst_lg_symm_cur - Cst_lg_symm[i])), 1e-6*χlarge^1.2, "symm")
    @test zrtest(mean(abs.(Cst_sm_symm_cur - Cst_sm_symm[i])), 1e-6*χsmall^1.2, "symm")
end
end


text:
# Read an explanation of this implementation here:
#   https://nextjournal.com/erik-engheim/implementation-of-automatic-differentiation
export Dual, derive, gradient, Params

import Base: +, -, *, /, sin, cos, show,
             convert, promote_rule,
             setproperty!, getproperty

struct Dual{T<:Real} <: Real
    x::T
    ϵ::T
end

"""
    Dual(x)

Create a dual number of `x`. Dual numbers allow us to keep track
of derivatives when performing calculations. If `f(x) = x` then `f'(x) = 1`.
That is why `Dual(x)` is equal to `Dual(x, 1)`.
"""
Dual(x) = Dual(x, one(x))

show(io::IO, d::Dual) = print(io, d.x, " + ", d.ϵ, "ϵ")

a::Dual + b::Dual = Dual(a.x + b.x, a.ϵ + b.ϵ)
a::Dual - b::Dual = Dual(a.x - b.x, a.ϵ - b.ϵ)

a::Dual * b::Dual = Dual(a.x * b.x, b.x * a.ϵ + a.x * b.ϵ)
a::Dual / b::Dual = Dual(a.x / b.x, (a.ϵ*b.x - b.ϵ*a.x) / b.x^2)

sin(d::Dual) = Dual(sin(d.x), cos(d.x) * d.ϵ)
cos(d::Dual) = Dual(cos(d.x), sin(d.x) * (-d.ϵ))

function promote_rule(::Type{Dual{T}}, ::Type{Dual{R}}) where {T,R}
    Dual{promote_type(T,R)}
end

function promote_rule(::Type{Dual{T}}, ::Type{R}) where {T,R}
    Dual{promote_type(T,R)}
end

function convert(::Type{Dual{T}}, x::Dual) where T
    Dual(convert(T, x.x), convert(T, x.ϵ))
end

function convert(::Type{Dual{T}}, x::Real) where T
    Dual(convert(T, x), zero(T))
end

"""
    derive(f, x)

Given a function `y = f(x)` then `y` is the dependent variable and `x` is
the independent variable. `derive(f, x)` calculates the derivative of the function `f`
with respect to the indepedent variable `x`.
"""
derive(f, x) = f(Dual(x)).ϵ

"""
    gradient(f, xs...)

Calculate the gradient of the function `f` with respect to the variables `xs`.
So say we got `xs = [x₁, x₂, ..., xₙ]`. Then the gradient `gs = ∇f(x₁, x₂, ..., xₙ)`
is such that `gs[1] = ∂f(xs...)/∂x₁`, `gs[2] = ∂f(xs...)/∂x₂` etc. So the component `gs[i]`
of the derivative is partial derivative of `f` with respect to `xs[i]`.
"""
function gradient(f, xs...)
    gs = []
    for (i, x) in enumerate(xs)
        if i == 1
            push!(gs, f(Dual(x), xs[2:end]...))
        else
            push!(gs, f(xs[1:(i-1)]..., Dual(x), xs[(i+1):end]...))
        end
    end
    [g.ϵ for g in gs]
end

const Params = Dict{Symbol, Real}

function getproperty(params::Params, key::Symbol)
    if hasfield(Params, key)
        getfield(params, key)
    else
        getindex(params, key)
    end
end

function setproperty!(params::Params, key::Symbol, x::Real)
    if hasfield(Params, key)
        setfield!(params, key, x)
    else
        setindex!(params, x, key)
    end
end

"""
    gradient(f, params::Params)

Calculate the gradient of the function `f` with respect to the parameters `params`.
"""
function gradient(f, params::Params)
    gs = Params()       # To hold gradient

    for k in keys(params)
        backup = params[k]
        params[k] = Dual(backup) # So we can derivate with respect to paramer `k`
        gs[k] = f().ϵ
        params[k] = backup
    end

    gs
end

ext: jl


text:
using Oceananigans
using Oceananigans.Units
using Oceananigans.BoundaryConditions: fill_halo_regions!, ImpenetrableBoundaryCondition
using Printf
using GLMakie

# Buoyancy profile
bᵢ(x, y, z) = 1e-5 * z
b_particle = bᵢ(0, 0, 0) # particle buoyancy
ν = 1e-4
κ_particle = 1e-3
R_particle = 1.0
Δt = 1e-2

grid = RectilinearGrid(size=128, z=(-4, 4), topology=(Flat, Flat, Bounded))

no_penetration = ImpenetrableBoundaryCondition()
slip_bcs = FieldBoundaryConditions(grid, (Center, Center, Face),
                                   top=no_penetration, bottom=no_penetration)
w_slip = ZFaceField(grid, boundary_conditions=slip_bcs)
sinking = AdvectiveForcing(WENO5(), w=w_slip)

model = NonhydrostaticModel(; grid,
                            tracers = (:b, :P),
                            buoyancy = BuoyancyTracer(),
                            closure = ScalarDiffusivity(κ=(b=0.0, P=κ_particle)),
                            forcing = (; P = sinking))

Pᵢ(x, y, z) = exp(-(z + 1)^2)
set!(model, b=bᵢ, P=Pᵢ)

simulation = Simulation(model; Δt, stop_iteration=0)

b = model.tracers.b
w_slip_op = 2/9 * (b_particle - b) / ν * R_particle^2

function compute_slip_velocity!(sim)
    w_slip .= w_slip_op
    fill_halo_regions!(w_slip)
    return nothing
end

simulation.callbacks[:slip] = Callback(compute_slip_velocity!)

progress(sim) = @info @sprintf("Iter: %d, time: %.2e, max|w_slip|: %.2e",
                               iteration(sim), time(sim), maximum(abs, w_slip))
simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

z = znodes(Center, grid)
P = interior(model.tracers.P, 1, 1, :)

fig = Figure()
ax = Axis(fig[1, 1],
          xlabel = "Particle concentration",
          ylabel = "z",
          title = "Settling diffusion particle concentration at t=0")
xlims!(ax, -1, 4)

ℓ = lines!(ax, P, z)
display(fig)

function update_plot!(sim)
    ℓ.input_args[1][] = interior(sim.model.tracers.P, 1, 1, :)
    ax.title[] = @sprintf("Particle concentration at t=%.2e", time(sim))
end

#simulation.callbacks[:plot] = Callback(update_plot!, IterationInterval(1000))

record(fig, "settling_diffusion.mp4", 1:100, framerate=24) do nn
    simulation.stop_iteration += 200
    run!(simulation)
    update_plot!(simulation)
end


text:
if basename(pwd()) == "aoc"
    cd("2017/16")
end

function part1(filename::AbstractString, lastchar = 'p')
    chars = ['a':lastchar...]
    lines = split(readchomp(filename), ',')
    for line in lines
        if startswith(line, 's')
            chars = circshift(chars, parse(Int, line[2:end]))
        elseif startswith(line, 'x')
            i, j = parse.(Int, split(line[2:end], '/')) .+ 1
            chars[i], chars[j] = chars[j], chars[i]
        elseif startswith(line, 'p')
            replace!(chars, line[2] => line[4], line[4] => line[2])
        end
    end
    join(chars)
end

part1("example.txt", 'e')
part1("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
module NaiveResamplers
export naivefilt

# Naive rational resampler
function naivefilt( h::Vector, x::Vector, resamplerate::Rational = 1//1 )

    upfactor     = numerator( resamplerate )
    downfactor   = denominator( resamplerate )
    xLen         = length( x )
    xZeroStuffed = zeros( eltype(x), length( x ) * upfactor )

    for n in 0:length(x)-1
        xZeroStuffed[ n*upfactor+1 ] = x[ n+1 ]
    end

    y = Base.filt( h, 1.0, xZeroStuffed )
    y = [ y[n] for n = 1:downfactor:length( y ) ]
end

# Naive arbitrary resampler
function naivefilt( h::Vector, x::Vector, resamplerate::AbstractFloat, numfilters::Integer = 32 )
    xLen          = length( x )
    xInterpolated = naivefilt( h, x, numfilters//1 )
    xLen          = length( xInterpolated )
    yLen          = ceil( Int, xLen * resamplerate )
    y             = similar( x, yLen )
    yIdx          = 1
    xIdx          = 1
    α             = 0.0
    (δ, 𝜙Stride)  = modf( numfilters/resamplerate )
    𝜙Stride       = Int( 𝜙Stride )

    while xIdx < xLen
        yLower  = xInterpolated[xIdx]
        yUpper  = xInterpolated[xIdx+1]
        y[yIdx] = yLower + α*( yUpper - yLower )
        yIdx   += 1
        α      += δ
        xIdx   += floor( Int, α ) + 𝜙Stride
        α       = mod( α, 1.0 )
    end

    resize!( y, yIdx-1 )

    return y
end

end # module


text:
"""
Read meta data for NITF 2.1

    nitfparse21(fid::IOStream)
"""
function nitfparse21(fid::IOStream)
    nitf_meta = [];
    fields = ["FHDR"   "FileProfileNameAndVersion"          9;
              "CLEVEL"  "ComplexityLevel"                  2;
              "STYPE"   "StandardType"                     4;
              "OSTAID"  "OriginatingStationID"            10;
              "FDT"     "FileDateAndTime"                 14;
              "FTITLE"  "FileTitle"                       80;
              "FSCLAS"  "FileSecurityClassification"       1;
              "FSCLSY"  "FileSecurityClassificationSystem" 2;
              "FSCODE"  "FileCodewords"                   11;
              "FSCTLH"  "FileControlAndHandling"           2;
              "FSREL"   "FileReleasingInstructions"       20;
              "FSDCTP"  "FileDeclassificationType"         2;
              "FSDCDT"  "FileDeclassificationDate"         8;
              "FSDCXM"  "FileDeclassificationExemption"    4;
              "FSDG"    "FileDowngrade"                    1;
              "FSDGT"   "FileDowngradeDate"                8;
              "FSCLTX"  "FileClassificationText"          43;
              "FSCATP"  "FileClassificationAuthorityType"  1;
              "FSCAUT"  "FileClassificationAuthority"     40;
              "FSCRSN"  "FileClassificationReason"         1;
              "FSSRDT"  "FileSecuritySourceDate"           8;
              "FSCTLN"  "FileSecurityControlNumber"       15;
              "FSCOP"   "FileCopyNumber"                   5;
              "FSCPYS"  "FileNumberOfCopies"               5;
              "ENCRYP"  "Encryption"                       1;
              "FBKGC"   "FileBackgroundColor"              3;
              "ONAME"   "OriginatorName"                  24;
              "OPHONE"  "OriginatorPhoneNumber"           18;
              "FL"      "FileLength"                      12;
              "HL"      "NITFFileHeaderLength"             6;
              "NUMI"    "NumberOfImages"                   3];

    nitf_read_meta!(nif_meta, fields, fid);

    numi = float(join(nitf_meta[end]["value"]));
    nitf_meta, im_lengths = processtopimagesubheadermeta(numi, nitf_meta, fid);
end

"""
Read meta data for NITF 2.1

    nitfparse20(fid::IOStream)
"""
function nitfparse20(fid::IOStream)
    nitf_meta = Dict();
end

"""
Append attributes to metadata type

    nitf_read_meta!(meta::Array, fields::Array, fid)
    :param meta: Array
    :param fields:: Array
    :param fid: IOStream
"""
function nitf_read_meta!(meta::Dict, fields::Array, fid::IOStream)
    # as a an optimization, convert the lengths part of the array to a vector
    data_lengths = vec(fields[:,3]);

    # read all of teh data at once. Extract individual values in the loop.
    data = readbytes(fid, sum(data_lengths));
    data = map(Char, data);

    data_offset = 0;
    for p = 1:size(fields,1)
        start = data_offset + 1;
        stop = data_offset + data_lengths[p];
        meta = [meta; Dict("name"=>fields[p,1],"vname"=>fields[p,2],
                           "value"=>data[start:stop])];
        data_offset = stop;
    end
end

"""
Add the nth image subheader information to the nitf_meta Dict. This Dict
contains the image subheader length and image length for all the images in the
filed

    processtopimagesubheadermeta!(nitf_meta, numi, fid)
    :param numi: number of images
    :param nitf_meta: Dict
    :param fid: IOStream

    :return nitf_meta
"""
function processtopimagesubheadermeta!(nitf_meta, numi, fid)
    nitf_metaISL = [];
    if (numi > 0)
        nitf_meta[end]["vname"] = "ImageAndImageSubheaderLengths";

        fields = ['LISH%03d' 'LengthOfNthImageSubheader'   6;
                  'LI%010d'  'LengthOfNthImage'           10];

        for current = 1:numi
            # setup
            name_t = @sprintf("Image%03d", current);

            # parse the data
            temp = [];
            nitf_read_meta_multi!(temp, fields, fid, current);
            nitf_metaISL = [nitf_metaISL; Dict("name"=>name_t, "vname"=> [name_t "ImageAndSubheaderLengths", "value"=>temp])];
        end

        nitf_meta = [nitf_meta; Dict("name"=>"ImageAndSubHeaderLengths",
                                     "vname"=>fields[p,2],
                                     "value"=>nitf_metaISL)];
    end
end

"""
Append attributes to metadata array

    nitf_read_meta_multi!(meta, fields, fid, index)
"""
function nitf_read_meta_multi!(meta, fields, fid, index)
    # as a an optimization, convert the lengths part of the array to a vector
    data_lengths = vec(fields[:,3]);

    # read all of teh data at once. Extract individual values in the loop.
    data = readbytes(fid, sum(data_lengths));
    data = map(Char, data);

    data_offset = 0;
    for p = 1:size(fields,1)
        start = data_offset + 1;
        stop = data_offset + data_lengths[p];
        name_t = @sprintf(fields[p,1], index);
        name_t1 = @sprintf(fields[p,2], index);
        meta = [meta; Dict("name"=>name_t,"vname"=>name_t1,
                           "value"=>data[start:stop])];
        data_offset = stop;
    end
end


text:
using Vega
using DataFrames

function generate_d_vs(q::Tait, x_min, y_min)
    ids = sort!(collect(vertices(q)))
    vts = [Quon.is_open_vertex(q, v) ? "B" : (Quon.is_genus(q, v) ? "G" : "N") for v in ids]
    xs = [(q.locations[v][1]-x_min)*100 for v in ids]
    ys = [(q.locations[v][2]-y_min)*100 for v in ids]
    degs = [length(trace_vertex(q, v)) for v in ids]
    shifts = zeros(length(ids))
    d_vs = DataFrame(id = ids, v_type = vts, x = xs, y = ys, shift = shifts, degree = degs)
    return d_vs
end

function generate_d_edges(q::Tait)
    hes = sort!(collect(half_edges(q)))
    fs = [face(q, he) for he in hes]
    srcs = [src(q, he) for he in hes]
    dsts = [dst(q, he) for he in hes]
    rs1 = Dict{Int, Float64}()
    rs2 = Dict{Int, Float64}()
    for v in vertices(q)
        hes_v = trace_vertex(q, v)
        n = length(hes_v)
        for i in eachindex(hes_v)
            he = hes_v[i]
            rs1[he] = -(i-1)*2pi/n
            rs2[twin(q, he)] = -(i-1)*2pi/n
        end
    end
    rot_srcs = [rs1[he] for he in hes]
    rot_dsts = [rs2[he] for he in hes]
    d_edges = DataFrame(id = hes, src = srcs, dst = dsts, f = fs, rot_src = rot_srcs, rot_dst = rot_dsts)
    return d_edges
end

function generate_frame(q::Tait)
    x_min, x_max = (0.0, 1.0)
    y_min, y_max = (0.0, 1.0)
    for v in vertices(q)
        loc = q.locations[v]
        loc[1] > x_max && (x_max = loc[1])
        loc[1] < x_min && (x_min = loc[1])
        loc[2] > y_max && (y_max = loc[2])
        loc[2] < y_min && (y_min = loc[2])
    end
    return x_min, x_max, y_min, y_max
end

function plot_vega(q::Tait)
    x_min, x_max, y_min, y_max = generate_frame(q)
    h, w = (x_max-x_min)*100, (y_max-y_min)*100
    println(h, w)
    d_edges = generate_d_edges(q)
    d_vs = generate_d_vs(q, x_min, y_min)
    return @vgplot(
        $schema="https://vega.github.io/schema/vega/v5.json",
        height=h,
        width=w,
        padding=50,
        autosize="none",
        data=[
            {
                name="vertices",
                values = d_vs,
                on=[
                    {
                        modify="whichSymbol",
                        values="newLoc && {x: newLoc.x, y: newLoc.y}",
                        trigger="newLoc"
                    },
                    {
                        modify="whichSymbol",
                        values="Clicked && {shift: whichSymbol.shift + 3.14/whichSymbol.degree}",
                        trigger="Clicked"
                    }
                ],
                transform=[
                    {
                        as="shape",
                        expr="'circle'",
                        type="formula"
                    },
                    {
                        as="color",
                        expr="datum.v_type === 'B' ? '#389826' : (datum.v_type === 'N' ? '#CB3C33' : (datum.v_type === 'G' ? 'white' : 'white'))",
                        type="formula"
                    }
                ]
            },
            {
                name="edges",
                values=d_edges,
                transform=[
                    {
                        key="id",
                        fields=[
                            "dst",
                            "src"
                        ],
                        as=[
                            "t",
                            "s"
                        ],
                        from="vertices",
                        type="lookup"
                    },
                    {
                        as="c1",
                        expr="{'x': datum.s.x + ctRadius * cos(datum.rot_src + datum.s.shift), 'y': datum.s.y + ctRadius * sin(datum.rot_src + datum.s.shift)}",
                        type="formula"
                    },
                    {
                        as="c2",
                        expr="{'x': datum.t.x + ctRadius * cos(datum.rot_dst + datum.t.shift), 'y': datum.t.y + ctRadius * sin(datum.rot_dst + datum.t.shift)}",
                        type="formula"
                    },
                    {
                        as="m",
                        expr="{'x': datum.s.x*1/8 + datum.t.x*1/8 + datum.c1.x*3/8 + datum.c2.x*3/8, 'y': datum.s.y*1/8 + datum.t.y*1/8 + datum.c1.y*3/8 + datum.c2.y*3/8}",
                        type="formula"
                    },
                    {
                        as="mypath",
                        expr="join(['M', toString(datum.s.x), toString(datum.s.y), 'C', toString(datum.c1.x), toString(datum.c1.y), toString(datum.c2.x), toString(datum.c2.y), toString(datum.t.x), toString(datum.t.y)], ' ')",
                        type="formula"
                    }
                ]
            },
            {
                name="fLocs",
                source="edges",
                transform=[
                    {
                        expr="datum.f != 0",
                        type="filter"
                    },
                    {
                        fields=[ "s.x", "s.y", "t.x", "t.y", "c1.x", "c1.y", "c2.x", "c2.y" ],
                        ops=[ "mean", "mean", "mean", "mean", "mean", "mean", "mean", "mean" ],
                        as=[ "sx", "sy", "tx", "ty", "c1x", "c1y", "c2x", "c2y" ],
                        groupby=[ "f" ],
                        type="aggregate"
                    }
                ]
            }
        ],
        signals=[
            {
                name="vSize",
                bind={
                    step=10,
                    max=1000,
                    min=100,
                    input="range"
                },
                value=500
            },
            {
                name="strokeWidth",
                bind={
                    step=0.05,
                    max=3,
                    min=0,
                    input="range"
                },
                value=1.5
            },
            {
                name="edgeWidth",
                bind={
                    step=0.05,
                    max=3,
                    min=0,
                    input="range"
                },
                value=1.5
            },
            {
                name="ctRadius",
                bind={
                    step=5,
                    max=100,
                    min=0,
                    input="range"
                },
                value=50
            },
            {
                name="whichSymbol",
                on=[
                    {
                        events="symbol:mousedown",
                        update="datum"
                    },
                    {
                        events="*:mouseup",
                        update="{}"
                    }
                ],
                value={}
            },
            {
                name="newLoc",
                on=[
                    {
                        events="symbol:mouseout[!event.buttons], window:mouseup",
                        update="false"
                    },
                    {
                        events="symbol:mouseover",
                        update="{x: x(), y: y()}"
                    },
                    {
                        events="[symbol:mousedown, window:mouseup] > window:mousemove!",
                        update="{x: x(), y: y()}"
                    }
                ],
                value=false,
                description="State variable for active node newLoc status."
            },
            {
                name="Clicked",
                on=[
                    {
                        events="symbol:mouseout[!event.buttons], window:mouseup",
                        update="false"
                    },
                    {
                        events="symbol:click",
                        update="true"
                    }
                ],
                value=false,
                description="State variable for active node newLoc status."
            }
        ],
        marks=[
            {
                encode={
                    update={
                        strokeWidth={ signal="edgeWidth" },
                        path={ field="mypath" }
                    },
                    enter={
                        stroke={ value="black" }
                    }
                },
                from={ data="edges" },
                type="path"
            },
            {
                encode={
                    update={
                        stroke={ value="black" },
                        x={ field="x" },
                        strokeWidth={ signal="strokeWidth" },
                        size={ signal="vSize" },
                        y={ field="y" }
                    },
                    enter={
                        shape={ field="shape" },
                        fill={ field="color" }
                    }
                },
                from={ data="vertices" },
                type="symbol"
            },
            {
                encode={
                    update={
                        x={ signal="datum.sx*1/2 + datum.tx*1/2" },
                        text={ field="f" },
                        y={ signal="datum.sy*1/2 + datum.ty*1/2" }
                    },
                    enter={
                        baseline={ value="middle" },
                        align={ value="center" }
                    }
                },
                from={ data="fLocs" },
                type="text"
            },
            {
                encode={
                    update={
                        x={ field="m.x" },
                        y={ field="m.y" }
                    }
                },
                from={ data="edges" },
                type="symbol"
            }
        ]
    )
end


text:
struct BS3ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a32::T
  a41::T
  a42::T
  a43::T
  c1::T2
  c2::T2
  btilde1::T
  btilde2::T
  btilde3::T
  btilde4::T
end

"""
constructBogakiShampine3()

Constructs the tableau object for the Bogakai-Shampine Order 2/3 method.
"""
function BS3ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloats})
  a21 = convert(T,0.5)
  a32 = convert(T,0.75)
  a41 = convert(T,0.2222222222222222)
  a42 = convert(T,0.3333333333333333)
  a43 = convert(T,0.4444444444444444)
  c1 = convert(T2,0.5)
  c2 = convert(T2,0.75)
  # b1 = convert(T,0.2916666666666667)
  # b2 = convert(T,0.25)
  # b3 = convert(T,0.3333333333333333)
  # b4 = convert(T,0.125)
  btilde1 = convert(T,0.06944444444444445)
  btilde2 = convert(T,-0.08333333333333333)
  btilde3 = convert(T,-0.1111111111111111)
  btilde4 = convert(T,0.125)
  BS3ConstantCache(a21,a32,a41,a42,a43,c1,c2,btilde1,btilde2,btilde3,btilde4)
end

"""
constructBogakiShampine3()

Constructs the tableau object for the Bogakai-Shampine Order 2/3 method.
"""
function BS3ConstantCache(T::Type,T2::Type)
  a21 = convert(T,1//2)
  a32 = convert(T,3//4)
  a41 = convert(T,2//9)
  a42 = convert(T,1//3)
  a43 = convert(T,4//9)
  c1 = convert(T2,1//2)
  c2 = convert(T2,3//4)
  # b1 = convert(T,7//24)
  # b2 = convert(T,1//4)
  # b3 = convert(T,1//3)
  # b4 = convert(T,1//8)
  btilde1 = convert(T,5//72)
  btilde2 = convert(T,-1//12)
  btilde3 = convert(T,-1//9)
  btilde4 = convert(T,1//8)
  BS3ConstantCache(a21,a32,a41,a42,a43,c1,c2,btilde1,btilde2,btilde3,btilde4)
end

struct OwrenZen3ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  c1::T2
  c2::T2
  btilde1::T
  btilde2::T
  btilde3::T
  r13::T
  r12::T
  r23::T
  r22::T
  r33::T
  r32::T
end

function OwrenZen3ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloats})
  a21 = convert(T,0.5217391304347826)
  a31 = convert(T,-0.18133333333333335)
  a32 = convert(T,0.9813333333333333)
  a41 = convert(T,0.2152777777777778)
  a42 = convert(T,0.4592013888888889)
  a43 = convert(T,0.3255208333333333)
  c1 = convert(T2,0.5217391304347826)
  c2 = convert(T2,0.8)
  # b1 = convert(T,0.041666666666666664)
  # b2 = convert(T,0.9583333333333334)
  btilde1 = convert(T,-0.1736111111111111)
  btilde2 = convert(T,0.4991319444444444)
  btilde3 = convert(T,-0.3255208333333333)
  r13 = convert(T,0.5694444444444444)
  r12 = convert(T,-1.3541666666666667)
  r23 = convert(T,-0.9184027777777778)
  r22 = convert(T,1.3776041666666667)
  r33 = convert(T,-0.6510416666666666)
  r32 = convert(T,0.9765625)
  OwrenZen3ConstantCache(a21,a31,a32,a41,a42,a43,c1,c2,btilde1,btilde2,btilde3,
                         r13,r12,r23,r22,r33,r32)
end

function OwrenZen3ConstantCache(T, T2)
  a21 = convert(T,12//23)
  a31 = convert(T,-68//375)
  a32 = convert(T,368//375)
  a41 = convert(T,31//144)
  a42 = convert(T,529//1152)
  a43 = convert(T,125//384)
  c1 = convert(T2,12//23)
  c2 = convert(T2,4//5)
  # b1 = convert(T,1//24)
  # b2 = convert(T,23//24)
  btilde1 = convert(T,-25//144)
  btilde2 = convert(T,575//1152)
  btilde3 = convert(T,-125//384)
  r13 = convert(T,41//72)
  r12 = convert(T,-65//48)
  r23 = convert(T,-529//576)
  r22 = convert(T,529//384)
  r33 = convert(T,-125//192)
  r32 = convert(T,125//128)
  OwrenZen3ConstantCache(a21,a31,a32,a41,a42,a43,c1,c2,btilde1,btilde2,btilde3,
                         r13,r12,r23,r22,r33,r32)
end

struct OwrenZen4ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a63::T
  a64::T
  a65::T
  c1::T2
  c2::T2
  c3::T2
  c4::T2
  btilde1::T
  btilde3::T
  btilde4::T
  btilde5::T
  r14::T
  r13::T
  r12::T
  r34::T
  r33::T
  r32::T
  r44::T
  r43::T
  r42::T
  r54::T
  r53::T
  r52::T
  r64::T
  r63::T
  r62::T
end

function OwrenZen4ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloats})
    a21 = convert(T,0.16666666666666666)
    a31 = convert(T,0.03214024835646457)
    a32 = convert(T,0.26515704894083275)
    a41 = convert(T,0.6895990230002036)
    a42 = convert(T,-1.6993690209647874)
    a43 = convert(T,1.6568288214939955)
    a51 = convert(T,-0.09002509947964493)
    a52 = convert(T,0.6817777777777778)
    a53 = convert(T,-0.2402791551882461)
    a54 = convert(T,0.5151931435567799)
    a61 = convert(T,0.08990252172070354)
    a63 = convert(T,0.4360623278236915)
    a64 = convert(T,0.1842858372687918)
    a65 = convert(T,0.2897493131868132)
    c1 = convert(T2,0.16666666666666666)
    c2 = convert(T2,0.2972972972972973)
    c3 = convert(T2,0.6470588235294118)
    c4 = convert(T2,0.8666666666666667)
    # b1 = convert(T,0.27823691460055094)
    # b3 = convert(T,-0.09428374655647383)
    # b4 = convert(T,0.8160468319559229)
    btilde1 = convert(T,0.18833439287984743)
    btilde3 = convert(T,-0.5303460743801653)
    btilde4 = convert(T,0.6317609946871311)
    btilde5 = convert(T,-0.2897493131868132)
    r14 = convert(T,-1.0513495872621479)
    r13 = convert(T,2.922894131082889)
    r12 = convert(T,-2.7816420221000375)
    r34 = convert(T,2.4266369235379472)
    r33 = convert(T,-5.725398502723277)
    r32 = convert(T,3.7348239070090217)
    r44 = convert(T,-0.7705135914137909)
    r43 = convert(T,1.1724555082899983)
    r42 = convert(T,-0.21765607960741548)
    r54 = convert(T,-2.8643157295948325)
    r53 = convert(T,5.149132832816039)
    r52 = convert(T,-1.9950677900343932)
    r64 = convert(T,2.2595419847328246)
    r63 = convert(T,-3.519083969465649)
    r62 = convert(T,1.2595419847328244)
    OwrenZen4ConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                           a61,a63,a64,a65,c1,c2,c3,c4,btilde1,btilde3,btilde4,btilde5,
                           r14,r13,r12,r34,r33,r32,r44,r43,r42,
                           r54,r53,r52,r64,r63,r62)
end

function OwrenZen4ConstantCache(T, T2)
  a21 = convert(T,1//6)
  a31 = convert(T,44//1369)
  a32 = convert(T,363//1369)
  a41 = convert(T,3388//4913)
  a42 = convert(T,-8349//4913)
  a43 = convert(T,8140//4913)
  a51 = convert(T,-36764//408375)
  a52 = convert(T,767//1125)
  a53 = convert(T,-32708//136125)
  a54 = convert(T,210392//408375)
  a61 = convert(T,1697//18876)
  a63 = convert(T,50653//116160)
  a64 = convert(T,299693//1626240)
  a65 = convert(T,3375//11648)
  c1 = convert(T2,1//6)
  c2 = convert(T2,11//37)
  c3 = convert(T2,11//17)
  c4 = convert(T2,13//15)
  # b1 = convert(T,101//363)
  # b3 = convert(T,-1369//14520)
  # b4 = convert(T,11849//14520)
  btilde1 = convert(T,1185//6292)
  btilde3 = convert(T,-4107//7744)
  btilde4 = convert(T,68493//108416)
  btilde5 = convert(T,-3375//11648)
  r14 = convert(T,-866577//824252)
  r13 = convert(T,1806901//618189)
  r12 = convert(T,-104217//37466)
  r34 = convert(T,12308679//5072320)
  r33 = convert(T,-2178079//380424)
  r32 = convert(T,861101//230560)
  r44 = convert(T,-7816583//10144640)
  r43 = convert(T,6244423//5325936)
  r42 = convert(T,-63869//293440)
  r54 = convert(T,-624375//217984)
  r53 = convert(T,982125//190736)
  r52 = convert(T,-1522125//762944)
  r64 = convert(T,296//131)
  r63 = convert(T,-461//131)
  r62 = convert(T,165//131)
  OwrenZen4ConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                         a61,a63,a64,a65,c1,c2,c3,c4,btilde1,btilde3,btilde4,btilde5,
                         r14,r13,r12,r34,r33,r32,r44,r43,r42,
                         r54,r53,r52,r64,r63,r62)
end

struct OwrenZen5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a62::T
  a63::T
  a64::T
  a65::T
  a71::T
  a72::T
  a73::T
  a74::T
  a75::T
  a76::T
  a81::T
  a83::T
  a84::T
  a85::T
  a86::T
  a87::T
  c1::T2
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  btilde1::T
  btilde3::T
  btilde4::T
  btilde5::T
  btilde6::T
  btilde7::T
  r15::T
  r14::T
  r13::T
  r12::T
  r35::T
  r34::T
  r33::T
  r32::T
  r45::T
  r44::T
  r43::T
  r42::T
  r55::T
  r54::T
  r53::T
  r52::T
  r65::T
  r64::T
  r63::T
  r62::T
  r75::T
  r74::T
  r73::T
  r72::T
  r85::T
  r84::T
  r83::T
  r82::T
end

function OwrenZen5ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloats})
    a21 = convert(T,0.16666666666666666)
    a31 = convert(T,0.0625)
    a32 = convert(T,0.1875)
    a41 = convert(T,0.25)
    a42 = convert(T,-0.75)
    a51 = convert(T,-0.75)
    a52 = convert(T,3.75)
    a53 = convert(T,-3.0)
    a54 = convert(T,0.5)
    a61 = convert(T,0.26895043731778423)
    a62 = convert(T,-0.7084548104956269)
    a63 = convert(T,0.8658892128279884)
    a64 = convert(T,0.15462307371928363)
    a65 = convert(T,0.06184922948771345)
    a71 = convert(T,-0.02947695035460993)
    a72 = convert(T,0.18500664893617022)
    a73 = convert(T,0.4802345261121857)
    a74 = convert(T,-0.5337849069148937)
    a75 = convert(T,-0.013090093085106383)
    a76 = convert(T,0.7861107753062541)
    a81 = convert(T,0.08783068783068783)
    a83 = convert(T,0.3006060606060606)
    a84 = convert(T,0.22777777777777777)
    a85 = convert(T,0.027777777777777776)
    a86 = convert(T,0.06218596218596219)
    a87 = convert(T,0.2938217338217338)
    c1 = convert(T2,0.16666666666666666)
    c2 = convert(T2,0.25)
    c3 = convert(T2,0.5)
    c4 = convert(T2,0.5)
    c5 = convert(T2,0.6428571428571429)
    c6 = convert(T2,0.875)
    # b1 = convert(T,-0.1111111111111111)
    # b3 = convert(T,1.2121212121212122)
    # b4 = convert(T,-1.75)
    # b5 = convert(T,-0.08333333333333333)
    # b6 = convert(T,1.7323232323232323)
    btilde1 = convert(T,-0.19894179894179895)
    btilde3 = convert(T,0.9115151515151515)
    btilde4 = convert(T,-1.9777777777777779)
    btilde5 = convert(T,-0.1111111111111111)
    btilde6 = convert(T,1.67013727013727)
    btilde7 = convert(T,-0.2938217338217338)
    r15 = convert(T,1.892063492063492)
    r14 = convert(T,-6.067155067155067)
    r13 = convert(T,7.282458282458283)
    r12 = convert(T,-4.0195360195360195)
    r35 = convert(T,-7.214545454545455)
    r34 = convert(T,20.676923076923078)
    r33 = convert(T,-20.31142191142191)
    r32 = convert(T,7.14965034965035)
    r45 = convert(T,7.866666666666666)
    r44 = convert(T,-18.78205128205128)
    r43 = convert(T,13.508547008547009)
    r42 = convert(T,-2.3653846153846154)
    r55 = convert(T,2)
    r54 = convert(T,-5.294871794871795)
    r53 = convert(T,4.534188034188034)
    r52 = convert(T,-1.2115384615384615)
    r65 = convert(T,-1.4924630924630924)
    r64 = convert(T,1.5785667324128863)
    r63 = convert(T,1.1958838881915805)
    r62 = convert(T,-1.219801565955412)
    r75 = convert(T,-7.051721611721612)
    r74 = convert(T,16.273203719357564)
    r73 = convert(T,-11.978886071193763)
    r72 = convert(T,3.0512256973795435)
    r85 = convert(T,4)
    r84 = convert(T,-8.384615384615385)
    r83 = convert(T,5.769230769230769)
    r82 = convert(T,-1.3846153846153846)
    OwrenZen5ConstantCache(a21,a31,a32,a41,a42,a51,a52,a53,
                           a54,a61,a62,a63,a64,a65,a71,a72,a73,
                           a74,a75,a76,a81,a83,a84,a85,a86,a87,
                           c1,c2,c3,c4,c5,c6,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,
                           r15,r14,r13,r12,r35,r34,r33,r32,r45,r44,
                           r43,r42,r55,r54,r53,r52,r65,r64,r63,
                           r62,r75,r74,r73,r72,r85,r84,r83,r82)
end

function OwrenZen5ConstantCache(T, T2)
  a21 = convert(T,1//6)
  a31 = convert(T,1//16)
  a32 = convert(T,3//16)
  a41 = convert(T,1//4)
  a42 = convert(T,-3//4)
  a43 = convert(T,1)
  a51 = convert(T,-3//4)
  a52 = convert(T,15//4)
  a53 = convert(T,-3)
  a54 = convert(T,1//2)
  a61 = convert(T,369//1372)
  a62 = convert(T,-243//343)
  a63 = convert(T,297//343)
  a64 = convert(T,1485//9604)
  a65 = convert(T,297//4802)
  a71 = convert(T,-133//4512)
  a72 = convert(T,1113//6016)
  a73 = convert(T,7945//16544)
  a74 = convert(T,-12845//24064)
  a75 = convert(T,-315//24064)
  a76 = convert(T,156065//198528)
  a81 = convert(T,83//945)
  a83 = convert(T,248//825)
  a84 = convert(T,41//180)
  a85 = convert(T,1//36)
  a86 = convert(T,2401//38610)
  a87 = convert(T,6016//20475)
  c1 = convert(T2,1//6)
  c2 = convert(T2,1//4)
  c3 = convert(T2,1//2)
  c4 = convert(T2,1//2)
  c5 = convert(T2,9//14)
  c6 = convert(T2,7//8)
  # b1 = convert(T,-1//9)
  # b3 = convert(T,40//33)
  # b4 = convert(T,-7//4)
  # b5 = convert(T,-1//12)
  # b6 = convert(T,343//198)
  btilde1 = convert(T,-188//945)
  btilde3 = convert(T,752//825)
  btilde4 = convert(T,-89//45)
  btilde5 = convert(T,-1//9)
  btilde6 = convert(T,32242//19305)
  btilde7 = convert(T,-6016//20475)
  r15 = convert(T,596//315)
  r14 = convert(T,-4969//819)
  r13 = convert(T,17893//2457)
  r12 = convert(T,-3292//819)
  r35 = convert(T,-1984//275)
  r34 = convert(T,1344//65)
  r33 = convert(T,-43568//2145)
  r32 = convert(T,5112//715)
  r45 = convert(T,118//15)
  r44 = convert(T,-1465//78)
  r43 = convert(T,3161//234)
  r42 = convert(T,-123//52)
  r55 = convert(T,2)
  r54 = convert(T,-413//78)
  r53 = convert(T,1061//234)
  r52 = convert(T,-63//52)
  r65 = convert(T,-9604//6435)
  r64 = convert(T,2401//1521)
  r63 = convert(T,60025//50193)
  r62 = convert(T,-40817//33462)
  r75 = convert(T,-48128//6825)
  r74 = convert(T,96256//5915)
  r73 = convert(T,-637696//53235)
  r72 = convert(T,18048//5915)
  r85 = convert(T,4)
  r84 = convert(T,-109//13)
  r83 = convert(T,75//13)
  r82 = convert(T,-18//13)
  OwrenZen5ConstantCache(a21,a31,a32,a41,a42,a51,a52,a53,
                         a54,a61,a62,a63,a64,a65,a71,a72,a73,
                         a74,a75,a76,a81,a83,a84,a85,a86,a87,
                         c1,c2,c3,c4,c5,c6,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,
                         r15,r14,r13,r12,r35,r34,r33,r32,r45,r44,
                         r43,r42,r55,r54,r53,r52,r65,r64,r63,
                         r62,r75,r74,r73,r72,r85,r84,r83,r82)
end

struct Tsit5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  c1::T2
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a62::T
  a63::T
  a64::T
  a65::T
  a71::T
  a72::T
  a73::T
  a74::T
  a75::T
  a76::T
  btilde1::T
  btilde2::T
  btilde3::T
  btilde4::T
  btilde5::T
  btilde6::T
  btilde7::T
  r11::T
  r12::T
  r13::T
  r14::T
  r22::T
  r23::T
  r24::T
  r32::T
  r33::T
  r34::T
  r42::T
  r43::T
  r44::T
  r52::T
  r53::T
  r54::T
  r62::T
  r63::T
  r64::T
  r72::T
  r73::T
  r74::T
end

function Tsit5ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloats})
  c1 =       convert(T2,0.161)
  c2 =       convert(T2,0.327)
  c3 =       convert(T2,0.9)
  c4 =       convert(T2,0.9800255409045097)
  c5 =       convert(T2,1)
  c6 =       convert(T2,1)
  a21=       convert(T,0.161)
  a31=       convert(T,-0.008480655492356989)
  a32=       convert(T,0.335480655492357)
  a41=       convert(T,2.8971530571054935)
  a42=       convert(T,-6.359448489975075)
  a43=       convert(T,4.3622954328695815)
  a51=       convert(T,5.325864828439257)
  a52=       convert(T,-11.748883564062828)
  a53=       convert(T,7.4955393428898365)
  a54=       convert(T,-0.09249506636175525)
  a61=       convert(T,5.86145544294642)
  a62=       convert(T,-12.92096931784711)
  a63=       convert(T,8.159367898576159)
  a64=       convert(T,-0.071584973281401)
  a65=       convert(T,-0.028269050394068383)
  a71=       convert(T,0.09646076681806523)
  a72=       convert(T,0.01)
  a73=       convert(T,0.4798896504144996)
  a74=       convert(T,1.379008574103742)
  a75=       convert(T,-3.290069515436081)
  a76=       convert(T,2.324710524099774)
  # b1 =       convert(T,0.09468075576583945)
  # b2 =       convert(T,0.009183565540343254)
  # b3 =       convert(T,0.4877705284247616)
  # b4 =       convert(T,1.234297566930479)
  # b5 =       convert(T,-2.7077123499835256)
  # b6 =       convert(T,1.866628418170587)
  # b7 =       convert(T,0.015151515151515152)
  btilde1 =  convert(T,-0.00178001105222577714)
  btilde2 =  convert(T,-0.0008164344596567469)
  btilde3 =  convert(T,0.007880878010261995)
  btilde4 =  convert(T,-0.1447110071732629)
  btilde5 =  convert(T,0.5823571654525552)
  btilde6 =  convert(T,-0.45808210592918697)
  btilde7 =  convert(T,0.015151515151515152)

  r11,r12,r13,r14,r22,r23,r24,r32,r33,r34,r42,r43,r44,r52,r53,r54,r62,r63,r64,r72,r73,r74 = Tsit5Interp(T)
  Tsit5ConstantCache(c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,btilde1,btilde2,btilde3,btilde4,btilde5,btilde6,btilde7,r11,r12,r13,r14,r22,r23,r24,r32,r33,r34,r42,r43,r44,r52,r53,r54,r62,r63,r64,r72,r73,r74)
end

function Tsit5ConstantCache(T::Type,T2::Type)
  c1 =       convert(T2,161//1000)
  c2 =       convert(T2,327//1000)
  c3 =       convert(T2,9//10)
  c4 =       convert(T2,big".9800255409045096857298102862870245954942137979563024768854764293221195950761080302604")
  c5 =       convert(T2,1)
  c6 =       convert(T2,1)
  a21=       convert(T,161//1000)
  a31=       convert(T,big"-.8480655492356988544426874250230774675121177393430391537369234245294192976164141156943e-2")
  a32=       convert(T,big".3354806554923569885444268742502307746751211773934303915373692342452941929761641411569")
  a41=       convert(T,big"2.897153057105493432130432594192938764924887287701866490314866693455023795137503079289")
  a42=       convert(T,big"-6.359448489975074843148159912383825625952700647415626703305928850207288721235210244366")
  a43=       convert(T,big"4.362295432869581411017727318190886861027813359713760212991062156752264926097707165077")
  a51=       convert(T,big"5.325864828439256604428877920840511317836476253097040101202360397727981648835607691791")
  a52=       convert(T,big"-11.74888356406282787774717033978577296188744178259862899288666928009020615663593781589")
  a53=       convert(T,big"7.495539342889836208304604784564358155658679161518186721010132816213648793440552049753")
  a54=       convert(T,big"-.9249506636175524925650207933207191611349983406029535244034750452930469056411389539635e-1")
  a61=       convert(T,big"5.861455442946420028659251486982647890394337666164814434818157239052507339770711679748")
  a62=       convert(T,big"-12.92096931784710929170611868178335939541780751955743459166312250439928519268343184452")
  a63=       convert(T,big"8.159367898576158643180400794539253485181918321135053305748355423955009222648673734986")
  a64=       convert(T,big"-.7158497328140099722453054252582973869127213147363544882721139659546372402303777878835e-1")
  a65=       convert(T,big"-.2826905039406838290900305721271224146717633626879770007617876201276764571291579142206e-1")
  a71=       convert(T,big".9646076681806522951816731316512876333711995238157997181903319145764851595234062815396e-1")
  a72=       convert(T,1//100)
  a73=       convert(T,big".4798896504144995747752495322905965199130404621990332488332634944254542060153074523509")
  a74=       convert(T,big"1.379008574103741893192274821856872770756462643091360525934940067397245698027561293331")
  a75=       convert(T,big"-3.290069515436080679901047585711363850115683290894936158531296799594813811049925401677")
  a76=       convert(T,big"2.324710524099773982415355918398765796109060233222962411944060046314465391054716027841")
  # b1 =       convert(T,big".9468075576583945807478876255758922856117527357724631226139574065785592789071067303271e-1")
  # b2 =       convert(T,big".9183565540343253096776363936645313759813746240984095238905939532922955247253608687270e-2")
  # b3 =       convert(T,big".4877705284247615707855642599631228241516691959761363774365216240304071651579571959813")
  # b4 =       convert(T,big"1.234297566930478985655109673884237654035539930748192848315425833500484878378061439761")
  # b5 =       convert(T,big"-2.707712349983525454881109975059321670689605166938197378763992255714444407154902012702")
  # b6 =       convert(T,big"1.866628418170587035753719399566211498666255505244122593996591602841258328965767580089")
  # b7 =       convert(T,1//66)
  btilde1 =  convert(T,big"-1.780011052225771443378550607539534775944678804333659557637450799792588061629796e-03")
  btilde2 =  convert(T,big"-8.164344596567469032236360633546862401862537590159047610940604670770447527463931e-04")
  btilde3 =  convert(T,big"7.880878010261996010314727672526304238628733777103128603258129604952959142646516e-03")
  btilde4 =  convert(T,big"-1.44711007173262907537165147972635116720922712343167677619514233896760819649515e-01")
  btilde5 =  convert(T,big"5.823571654525552250199376106520421794260781239567387797673045438803694038950012e-01")
  btilde6 =  convert(T,big"-4.580821059291869466616365188325542974428047279788398179474684434732070620889539e-01")
  btilde7 =  convert(T,1//66)

  r11,r12,r13,r14,r22,r23,r24,r32,r33,r34,r42,r43,r44,r52,r53,r54,r62,r63,r64,r72,r73,r74 = Tsit5Interp(T)
  Tsit5ConstantCache(c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,btilde1,btilde2,btilde3,btilde4,btilde5,btilde6,btilde7,r11,r12,r13,r14,r22,r23,r24,r32,r33,r34,r42,r43,r44,r52,r53,r54,r62,r63,r64,r72,r73,r74)
end

"""
Coefficients for the polynomial
bᵢΘ = ri1*Θ + ri2*Θ^2 + ri3*Θ^3 + ...

These are the coefficients of the expanded form of the polynomials from

Runge–Kutta pairs of order 5(4) satisfying only the first column
simplifying assumption

Ch. Tsitouras
"""
function Tsit5Interp(T::Type{<:CompiledFloats})
  r11 = convert(T,1.0)
  r12 = convert(T,-2.763706197274826)
  r13 = convert(T,2.9132554618219126)
  r14 = convert(T,-1.0530884977290216)

  r22 = convert(T,0.13169999999999998)
  r23 = convert(T,-0.2234)
  r24 = convert(T,0.1017)

  r32 = convert(T,3.9302962368947516)
  r33 = convert(T,-5.941033872131505)
  r34 = convert(T,2.490627285651253)

  r42 = convert(T,-12.411077166933676)
  r43 = convert(T,30.33818863028232)
  r44 = convert(T,-16.548102889244902)

  r52 = convert(T,37.50931341651104)
  r53 = convert(T,-88.1789048947664)
  r54 = convert(T,47.37952196281928)

  r62 = convert(T,-27.896526289197286)
  r63 = convert(T,65.09189467479366)
  r64 = convert(T,-34.87065786149661)

  r72 = convert(T,1.5)
  r73 = convert(T,-4)
  r74 = convert(T,2.5)

  return r11,r12,r13,r14,r22,r23,r24,r32,r33,r34,r42,r43,r44,r52,r53,r54,r62,r63,r64,r72,r73,r74
end

"""
Coefficients for the polynomial
bᵢΘ = ri1*Θ + ri2*Θ^2 + ri3*Θ^3 + ...

These are the coefficients of the expanded form of the polynomials from

Runge–Kutta pairs of order 5(4) satisfying only the first column
simplifying assumption

Ch. Tsitouras
"""
function Tsit5Interp(T::Type)
  r11 = convert(T,big"0.999999999999999974283372471559910888475488471328")
  r12 = convert(T,big"-2.763706197274825911336735930481400260916070804192")
  r13 = convert(T,big"2.91325546182191274375068099306808")
  r14 = convert(T,-1.0530884977290216)

  r22 = convert(T,big"0.13169999999999999727")
  r23 = convert(T,big"-0.22339999999999999818")
  r24 = convert(T,0.1017)

  r32 = convert(T,big"3.93029623689475152850687446709813398")
  r33 = convert(T,big"-5.94103387213150473470249202589458001")
  r34 = convert(T,big"2.490627285651252793")

  r42 = convert(T,big"-12.411077166933676983734381540685453484102414134010752")
  r43 = convert(T,big"30.3381886302823215981729903691836576")
  r44 = convert(T,big"-16.54810288924490272")

  r52 = convert(T,big"37.50931341651103919496903965334519631242339792120440212")
  r53 = convert(T,big"-88.1789048947664011014276693541209817")
  r54 = convert(T,big"47.37952196281928122")

  r62 = convert(T,big"-27.896526289197287805948263144598643896")
  r63 = convert(T,big"65.09189467479367152629021928716553658")
  r64 = convert(T,big"-34.87065786149660974")

  r72 = convert(T,1.5)
  r73 = convert(T,-4)
  r74 = convert(T,2.5)

  return r11,r12,r13,r14,r22,r23,r24,r32,r33,r34,r42,r43,r44,r52,r53,r54,r62,r63,r64,r72,r73,r74
end

struct BS5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  c1::T2
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a62::T
  a63::T
  a64::T
  a65::T
  a71::T
  a72::T
  a73::T
  a74::T
  a75::T
  a76::T
  a81::T
  a83::T
  a84::T
  a85::T
  a86::T
  a87::T
  bhat1::T
  bhat3::T
  bhat4::T
  bhat5::T
  bhat6::T
  btilde1::T
  btilde3::T
  btilde4::T
  btilde5::T
  btilde6::T
  btilde7::T
  btilde8::T
  c6::T2
  c7::T2
  c8::T2
  a91::T
  a92::T
  a93::T
  a94::T
  a95::T
  a96::T
  a97::T
  a98::T
  a101::T
  a102::T
  a103::T
  a104::T
  a105::T
  a106::T
  a107::T
  a108::T
  a109::T
  a111::T
  a112::T
  a113::T
  a114::T
  a115::T
  a116::T
  a117::T
  a118::T
  a119::T
  a1110::T
  r016::T
  r015::T
  r014::T
  r013::T
  r012::T
  r036::T
  r035::T
  r034::T
  r033::T
  r032::T
  r046::T
  r045::T
  r044::T
  r043::T
  r042::T
  r056::T
  r055::T
  r054::T
  r053::T
  r052::T
  r066::T
  r065::T
  r064::T
  r063::T
  r062::T
  r076::T
  r075::T
  r074::T
  r073::T
  r072::T
  r086::T
  r085::T
  r084::T
  r083::T
  r082::T
  r096::T
  r095::T
  r094::T
  r093::T
  r106::T
  r105::T
  r104::T
  r103::T
  r102::T
  r116::T
  r115::T
  r114::T
  r113::T
  r112::T
end

"""

An Efficient Runge-Kutta (4,5) Pair by P.Bogacki and L.F.Shampine
 Computers and Mathematics with Applications, Vol. 32, No. 6, 1996, pages 15 to 28
"""
function BS5ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloats})
  c1     =convert(T2,0.16666666666666666)
  c2     =convert(T2,0.2222222222222222)
  c3     =convert(T2,0.42857142857142855)
  c4     =convert(T2,0.6666666666666666)
  c5     =convert(T2,0.75)

  a21    =convert(T,0.16666666666666666)
  a31    =convert(T,0.07407407407407407)
  a32    =convert(T,0.14814814814814814)
  a41    =convert(T,0.13338192419825073)
  a42    =convert(T,-0.47230320699708456)
  a43    =convert(T,0.7674927113702624)
  a51    =convert(T,0.22895622895622897)
  a52    =convert(T,-0.36363636363636365)
  a53    =convert(T,0.2937062937062937)
  a54    =convert(T,0.5076405076405076)
  a61    =convert(T,0.026500355113636364)
  a62    =convert(T,0.23011363636363635)
  a63    =convert(T,0.10772747760052448)
  a64    =convert(T,0.1602190777972028)
  a65    =convert(T,0.225439453125)
  a71    =convert(T,0.18159821692916506)
  a72    =convert(T,-0.38707982536247293)
  a73    =convert(T,0.41288268560074054)
  a74    =convert(T,0.6409913191253497)
  a75    =convert(T,-1.0121439383514517)
  a76    =convert(T,1.1637515420586695)
  a81    =convert(T,0.07279265873015874)
  a83    =convert(T,0.28662437773692473)
  a84    =convert(T,0.19513621794871794)
  a85    =convert(T,0.008638392857142857)
  a86    =convert(T,0.3595655806182122)
  a87    =convert(T,0.07724277210884353)

  bhat1  =convert(T,-0.00234375)
  bhat3  =convert(T,0.0103760754048583)
  bhat4  =convert(T,-0.016490384615384615)
  bhat5  =convert(T,0.018984375)
  bhat6  =convert(T,-0.010526315789473684)
#  btilde1=convert(T,0.07084476451760402)
#  btilde2=convert(T,0)
#  btilde3=convert(T,0.2956730769230769)
#  btilde4=convert(T,0.17965747482208388)
#  btilde5=convert(T,0.029861111111111113)
#  btilde6=convert(T,0.3462886755067825)
#  btilde7=convert(T,0.07176240133870539)
  btilde1=convert(T,-0.0019478942125547064)
  btilde3=convert(T,0.009048699186152193)
  btilde4=convert(T,-0.015478743126634073)
  btilde5=convert(T,0.021222718253968254)
  btilde6=convert(T,-0.013276905111429694)
  btilde7=convert(T,-0.005480370770138146)
  btilde8=convert(T,0.005912495780636172)
  c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110 = BS5Interp(T,T2)
  r016,r015,r014,r013,r012,r036,r035,r034,r033,r032,r046,r045,r044,r043,r042,r056,r055,r054,r053,r052,r066,r065,r064,r063,r062,r076,r075,r074,r073,r072,r086,r085,r084,r083,r082,r096,r095,r094,r093,r106,r105,r104,r103,r102,r116,r115,r114,r113,r112 = BS5Interp_polyweights(T)
  BS5ConstantCache(c1,c2,c3,c4,c5,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,bhat1,bhat3,bhat4,bhat5,bhat6,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,btilde8,c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110,r016,r015,r014,r013,r012,r036,r035,r034,r033,r032,r046,r045,r044,r043,r042,r056,r055,r054,r053,r052,r066,r065,r064,r063,r062,r076,r075,r074,r073,r072,r086,r085,r084,r083,r082,r096,r095,r094,r093,r106,r105,r104,r103,r102,r116,r115,r114,r113,r112)
end

"""

An Efficient Runge-Kutta (4,5) Pair by P.Bogacki and L.F.Shampine
 Computers and Mathematics with Applications, Vol. 32, No. 6, 1996, pages 15 to 28
"""
function BS5ConstantCache(T::Type,T2::Type)
  c1     =convert(T2,1//6)
  c2     =convert(T2,2//9)
  c3     =convert(T2,3//7)
  c4     =convert(T2,2//3)
  c5     =convert(T2,3//4)

  a21    =convert(T,1//6)
  a31    =convert(T,2//27)
  a32    =convert(T,4//27)
  a41    =convert(T,183//1372)
  a42    =convert(T,-162//343)
  a43    =convert(T,1053//1372)
  a51    =convert(T,68//297)
  a52    =convert(T,-4//11)
  a53    =convert(T,42//143)
  a54    =convert(T,1960//3861)
  a61    =convert(T,597//22528)
  a62    =convert(T,81//352)
  a63    =convert(T,63099//585728)
  a64    =convert(T,58653//366080)
  a65    =convert(T,4617//20480)
  a71    =convert(T,174197//959244)
  a72    =convert(T,-30942//79937)
  a73    =convert(T,8152137//19744439)
  a74    =convert(T,666106//1039181)
  a75    =convert(T,-29421//29068)
  a76    =convert(T,482048//414219)
  a81    =convert(T,587//8064)
  a83    =convert(T,4440339//15491840)
  a84    =convert(T,24353//124800)
  a85    =convert(T,387//44800)
  a86    =convert(T,2152//5985)
  a87    =convert(T,7267//94080)

  bhat1  =convert(T,-3//1280)
  bhat3  =convert(T,6561//632320)
  bhat4  =convert(T,-343//20800)
  bhat5  =convert(T,243//12800)
  bhat6  =convert(T,-1//95)
  #  btilde1=convert(T,2479//34992)
  #  btilde2=convert(T,0)
  #  btilde3=convert(T,123//416)
  #  btilde4=convert(T,612941//3411720)
  #  btilde5=convert(T,43//1440)
  #  btilde6=convert(T,2272//6561)
  #  btilde7=convert(T,79937//1113912)
  btilde1=convert(T,-3817//1959552)
  btilde3=convert(T,140181//15491840)
  btilde4=convert(T,-4224731//272937600)
  btilde5=convert(T,8557//403200)
  btilde6=convert(T,-57928//4363065)
  btilde7=convert(T,-23930231//4366535040)
  btilde8=convert(T,3293//556956)
  c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110 = BS5Interp(T,T2)
  r016,r015,r014,r013,r012,r036,r035,r034,r033,r032,r046,r045,r044,r043,r042,r056,r055,r054,r053,r052,r066,r065,r064,r063,r062,r076,r075,r074,r073,r072,r086,r085,r084,r083,r082,r096,r095,r094,r093,r106,r105,r104,r103,r102,r116,r115,r114,r113,r112 = BS5Interp_polyweights(T)
  BS5ConstantCache(c1,c2,c3,c4,c5,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,bhat1,bhat3,bhat4,bhat5,bhat6,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,btilde8,c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110,r016,r015,r014,r013,r012,r036,r035,r034,r033,r032,r046,r045,r044,r043,r042,r056,r055,r054,r053,r052,r066,r065,r064,r063,r062,r076,r075,r074,r073,r072,r086,r085,r084,r083,r082,r096,r095,r094,r093,r106,r105,r104,r103,r102,r116,r115,r114,r113,r112)
end

"""

An Efficient Runge-Kutta (4,5) Pair by P.Bogacki and L.F.Shampine
 Computers and Mathematics with Applications, Vol. 32, No. 6, 1996, pages 15 to 28

Used in the lazy construction of the dense output

k9, k10, k11 are not computed until called in the dense routine
"""
function BS5Interp(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloats})

  c6    = convert(T2,0.5)
  c7    = convert(T2,0.8333333333333334)
  c8    = convert(T2,0.1111111111111111)
  a91   = convert(T,0.07405598958333333)
  a92   = convert(T,0)
  a93   = convert(T,0.28964485093442743)
  a94   = convert(T,0.12839214966168092)
  a95   = convert(T,-0.003779296875)
  a96   = convert(T,0.014230019493177388)
  a97   = convert(T,-0.03379371279761905)
  a98   = convert(T,0.03125)
  a101  = convert(T,-0.06358724036162344)
  a102  = convert(T,0.5742461924818869)
  a103  = convert(T,-0.06365063007249953)
  a104  = convert(T,0.043159777438314964)
  a105  = convert(T,0.8370112883898733)
  a106  = convert(T,-0.34045447246719235)
  a107  = convert(T,0.04926503818334922)
  a108  = convert(T,-0.006882677669165967)
  a109  = convert(T,-0.19577394258960973)
  a111  = convert(T,0.0636090772400987)
  a112  = convert(T,0.01057854182854183)
  a113  = convert(T,0.06600100945670531)
  a114  = convert(T,0.02048391555358402)
  a115  = convert(T,0.003682270330219549)
  a116  = convert(T,0.155258632271002)
  a117  = convert(T,-0.08509702513818027)
  a118  = convert(T,0.1)
  a119  = convert(T,-0.1)
  a1110 = convert(T,-0.12340531043086005)

  return c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110
end

"""

An Efficient Runge-Kutta (4,5) Pair by P.Bogacki and L.F.Shampine
 Computers and Mathematics with Applications, Vol. 32, No. 6, 1996, pages 15 to 28

Used in the lazy construction of the dense output

k9, k10, k11 are not computed until called in the dense routine
"""
function BS5Interp(T::Type,T2::Type)

  c6    = convert(T2,1//2)
  c7    = convert(T2,5//6)
  c8    = convert(T2,1//9)
  a91   = convert(T,455//6144)
  a92   = convert(T,0)
  a93   = convert(T,10256301//35409920)
  a94   = convert(T,2307361//17971200)
  a95   = convert(T,-387//102400)
  a96   = convert(T,73//5130)
  a97   = convert(T,-7267//215040)
  a98   = convert(T,1//32)
  a101  = convert(T,-837888343715//13176988637184)
  a102  = convert(T,30409415//52955362)
  a103  = convert(T,-48321525963//759168069632)
  a104  = convert(T,8530738453321//197654829557760)
  a105  = convert(T,1361640523001//1626788720640)
  a106  = convert(T,-13143060689//38604458898)
  a107  = convert(T,18700221969//379584034816)
  a108  = convert(T,-5831595//847285792)
  a109  = convert(T,-5183640//26477681)
  a111  = convert(T,98719073263//1551965184000)
  a112  = convert(T,1307//123552)
  a113  = convert(T,4632066559387//70181753241600)
  a114  = convert(T,7828594302389//382182512025600)
  a115  = convert(T,40763687//11070259200)
  a116  = convert(T,34872732407//224610586200)
  a117  = convert(T,-2561897//30105600)
  a118  = convert(T,1//10)
  a119  = convert(T,-1//10)
  a1110 = convert(T,-1403317093//11371610250)

  return c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110
end

"""
Coefficients for the polynomial
bᵢΘ = ri1*Θ + ri2*Θ^2 + ri3*Θ^3 + ...

These coefficients are taken from RKSuite

Note that RKSuite has an error: r081 should be 0
and r011 should be 1. This is pretty easy to spot
since the first order interpolation is linear from y₀.
"""
function BS5Interp_polyweights(T::Type{<:CompiledFloats})
  r016   = convert(T,-11.547607240534195)
  r015   = convert(T,36.89579791207878)
  r014   = convert(T,-45.470343197183475)
  r013   = convert(T,27.298136302807084)
  r012   = convert(T,-8.103191118438032)

  r036   = convert(T,-27.245925284262768)
  r035   = convert(T,74.85879078710211)
  r034   = convert(T,-72.2759709558427)
  r033   = convert(T,28.38602193195626)
  r032   = convert(T,-3.4362921012159933)

  r046   = convert(T,-14.307769126529058)
  r045   = convert(T,38.240038148817945)
  r044   = convert(T,-35.469854845413806)
  r043   = convert(T,13.060399314592578)
  r042   = convert(T,-1.3276772735189397)

  r056   = convert(T,-0.7296779223862557)
  r055   = convert(T,1.9817123385873385)
  r054   = convert(T,-1.964240395021645)
  r053   = convert(T,0.8847786781120115)
  r052   = convert(T,-0.16393430643430643)

  r066   = convert(T,19.121088881250603)
  r065   = convert(T,-65.9928405785889)
  r064   = convert(T,77.4690381021765)
  r063   = convert(T,-34.163041154825144)
  r062   = convert(T,3.9253203306051474)

  r076   = convert(T,-14.676126700680273)
  r075   = convert(T,42.17455357142857)
  r074   = convert(T,-42.4062818877551)
  r073   = convert(T,16.83892431972789)
  r072   = convert(T,-1.853826530612245)

  r086   = convert(T,19.453125)
  r085   = convert(T,-54.359375)
  r084   = convert(T,53.671875)
  r083   = convert(T,-21.078125)
  r082   = convert(T,2.3125)

  r096   = convert(T,18.333333333333332)
  r095   = convert(T,-39)
  r094   = convert(T,23)
  r093   = convert(T,-2.3333333333333335)

  r106   = convert(T,-23.400440940191388)
  r105   = convert(T,70.20132282057416)
  r104   = convert(T,-73.55422182095978)
  r103   = convert(T,30.106238940962648)
  r102   = convert(T,-3.3528990003856314)

  r116   = convert(T,35)
  r115   = convert(T,-105)
  r114   = convert(T,117)
  r113   = convert(T,-59)
  r112   = convert(T,12)

  return r016,r015,r014,r013,r012,r036,r035,r034,r033,r032,r046,r045,r044,r043,r042,r056,r055,r054,r053,r052,r066,r065,r064,r063,r062,r076,r075,r074,r073,r072,r086,r085,r084,r083,r082,r096,r095,r094,r093,r106,r105,r104,r103,r102,r116,r115,r114,r113,r112
end

"""
Coefficients for the polynomial
bᵢΘ = ri1*Θ + ri2*Θ^2 + ri3*Θ^3 + ...

These coefficients are taken from RKSuite

Note that RKSuite has an error: r081 should be 0
and r011 should be 1. This is pretty easy to spot
since the first order interpolation is linear from y₀.
"""
function BS5Interp_polyweights(T::Type)
  r016   = convert(T,-12134338393//1050809760)
  r015   = convert(T,12923488183//350269920)
  r014   = convert(T,-2722545893//59875200)
  r013   = convert(T,35856435071//1313512200)
  r012   = convert(T,-3547880131//437837400)

  r036   = convert(T,-33197340367//1218433216)
  r035   = convert(T,65150312289//870309440)
  r034   = convert(T,-27096444225//374902528)
  r033   = convert(T,4323308999//152304152)
  r032   = convert(T,-1046723109//304608304)

  r046   = convert(T,-284800997201//19905339168)
  r045   = convert(T,6343174409579//165877826400)
  r044   = convert(T,-201150852119//5671036800)
  r043   = convert(T,3249645975331//248816739600)
  r042   = convert(T,-55058055073//41469456600)

  r056   = convert(T,-540919//741312)
  r055   = convert(T,85695583//43243200)
  r054   = convert(T,-2903933//1478400)
  r053   = convert(T,3586937//4054050)
  r052   = convert(T,-1772261//10810800)

  r066   = convert(T,7157998304//374350977)
  r065   = convert(T,-41174140576//623918295)
  r064   = convert(T,413114104//5332635)
  r063   = convert(T,-9134977024//267393555)
  r062   = convert(T,2449079168//623918295)

  r076   = convert(T,-138073//9408)
  r075   = convert(T,94471//2240)
  r074   = convert(T,-1329861//31360)
  r073   = convert(T,792103//47040)
  r072   = convert(T,-7267//3920)

  r086   = convert(T,1245//64)
  r085   = convert(T,-3479//64)
  r084   = convert(T,3435//64)
  r083   = convert(T,-1349//64)
  r082   = convert(T,37//16)

  r096   = convert(T,55//3)
  r095   = convert(T,-39)
  r094   = convert(T,23)
  r093   = convert(T,-7//3)

  r106   = convert(T,-1774004627//75810735)
  r105   = convert(T,1774004627//25270245)
  r104   = convert(T,-26477681//359975)
  r103   = convert(T,11411880511//379053675)
  r102   = convert(T,-423642896//126351225)

  r116   = convert(T,35)
  r115   = convert(T,-105)
  r114   = convert(T,117)
  r113   = convert(T,-59)
  r112   = convert(T,12)

  return r016,r015,r014,r013,r012,r036,r035,r034,r033,r032,r046,r045,r044,r043,r042,r056,r055,r054,r053,r052,r066,r065,r064,r063,r062,r076,r075,r074,r073,r072,r086,r085,r084,r083,r082,r096,r095,r094,r093,r106,r105,r104,r103,r102,r116,r115,r114,r113,r112
end

struct DP5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a62::T
  a63::T
  a64::T
  a65::T
  a71::T
  a73::T
  a74::T
  a75::T
  a76::T
  btilde1::T
  btilde3::T
  btilde4::T
  btilde5::T
  btilde6::T
  btilde7::T
  c1::T2
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  d1::T
  d3::T
  d4::T
  d5::T
  d6::T
  d7::T
end

function DP5ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloats})
  a21 = convert(T,0.2)
  a31 = convert(T,0.075)
  a32 = convert(T,0.225)
  a41 = convert(T,0.9777777777777777)
  a42 = convert(T,-3.7333333333333334)
  a43 = convert(T,3.5555555555555554)
  a51 = convert(T,2.9525986892242035)
  a52 = convert(T,-11.595793324188385)
  a53 = convert(T,9.822892851699436)
  a54 = convert(T,-0.2908093278463649)
  a61 = convert(T,2.8462752525252526)
  a62 = convert(T,-10.757575757575758)
  a63 = convert(T,8.906422717743473)
  a64 = convert(T,0.2784090909090909)
  a65 = convert(T,-0.2735313036020583)
  a71 = convert(T,0.09114583333333333)
  a73 = convert(T,0.44923629829290207)
  a74 = convert(T,0.6510416666666666)
  a75 = convert(T,-0.322376179245283)
  a76 = convert(T,0.13095238095238096)
  # b1  = convert(T,0.08991319444444444)
  # b3  = convert(T,0.4534890685834082)
  # b4  = convert(T,0.6140625)
  # b5  = convert(T,-0.2715123820754717)
  # b6  = convert(T,0.08904761904761904)
  # b7  = convert(T,0.025)
  btilde1 = convert(T,-0.0012326388888888888)
  btilde3 = convert(T,0.0042527702905061394)
  btilde4 = convert(T,-0.03697916666666667)
  btilde5 = convert(T,0.05086379716981132)
  btilde6 = convert(T,-0.0419047619047619)
  btilde7 = convert(T,0.025)
  c1 = convert(T2,0.2)
  c2 = convert(T2,0.3)
  c3 = convert(T2,0.8)
  c4 = convert(T2,0.8888888888888888)
  c5 = convert(T2,1)
  c6 = convert(T2,1)
  d1,d3,d4,d5,d6,d7 = DP5_dense_ds(T)
  DP5ConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,c1,c2,c3,c4,c5,c6,d1,d3,d4,d5,d6,d7)
end

function DP5_dense_ds(T::Type{<:CompiledFloats})
  d1  = convert(T,-1.1270175653862835)
  d3  = convert(T,2.675424484351598)
  d4  = convert(T,-5.685526961588504)
  d5  = convert(T,3.5219323679207912)
  d6  = convert(T,-1.7672812570757455)
  d7  = convert(T,2.382468931778144)
  return d1,d3,d4,d5,d6,d7
end

function DP5ConstantCache(T::Type,T2::Type)
  a21 = convert(T,1//5)
  a31 = convert(T,3//40)
  a32 = convert(T,9//40)
  a41 = convert(T,44/45)
  a42 = convert(T,-56//15)
  a43 = convert(T,32//9)
  a51 = convert(T,19372//6561)
  a52 = convert(T,-25360//2187)
  a53 = convert(T,64448//6561)
  a54 = convert(T,-212//729)
  a61 = convert(T,9017//3168)
  a62 = convert(T,-355//33)
  a63 = convert(T,46732//5247)
  a64 = convert(T,49//176)
  a65 = convert(T,-5103//18656)
  a71 = convert(T,35//384)
  a73 = convert(T,500//1113)
  a74 = convert(T,125//192)
  a75 = convert(T,-2187//6784)
  a76 = convert(T,11//84)
  #  b1  = convert(T,5179//57600)
  #  b3  = convert(T,7571//16695)
  #  b4  = convert(T,393//640)
  #  b5  = convert(T,-92097//339200)
  #  b6  = convert(T,187//2100)
  #  b7  = convert(T,1//40)
  btilde1 = convert(T,-71//57600)
  btilde3 = convert(T,71//16695)
  btilde4 = convert(T,-71//1920)
  btilde5 = convert(T,17253//339200)
  btilde6 = convert(T,-22//525)
  btilde7 = convert(T,1//40)
  c1 = convert(T2,1//5)
  c2 = convert(T2,3//10)
  c3 = convert(T2,4//5)
  c4 = convert(T2,8//9)
  c5 = convert(T2,1)
  c6 = convert(T2,1)
  d1,d3,d4,d5,d6,d7 = DP5_dense_ds(T)
  DP5ConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,c1,c2,c3,c4,c5,c6,d1,d3,d4,d5,d6,d7)
end

function DP5_dense_ds(T)
  d1  = convert(T,-12715105075//11282082432)
  d3  = convert(T,87487479700//32700410799)
  d4  = convert(T,-10690763975//1880347072)
  d5  = convert(T,701980252875//199316789632)
  d6  = convert(T,-1453857185//822651844)
  d7  = convert(T,69997945//29380423)
  return d1,d3,d4,d5,d6,d7
end

#=
function DP5_dense_bs(T)
  b1  = convert(T,5179//57600)
  b3  = convert(T,7571//16695)
  b4  = convert(T,393//640)
  b5  = convert(T,-92097//339200)
  b6  = convert(T,187//2100)
  b7  = convert(T,1//40)
  return b1,b3,b4,b5,b6,b7
end
=#

"""
An Optimized Runge-Kutta method for the solution of Orbital Problems
by Z.A. Anastassi and T.E. Simos
Journal of Computational and Applied Mathematics, Volume 175, Issue 1, 1 March 2005, Pages 1 to 9.
"""
struct Anas5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a62::T
  a63::T
  a64::T
  a65::T
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  b1::T
  b3::T
  b4::T
  b5::T
  b6::T
end

function Anas5ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloats})
    a21 = convert(T,0.1)
    a31 = convert(T,-0.2222222222222222)
    a32 = convert(T,0.5555555555555556)
    a41 = convert(T,3.111111111111111)
    a42 = convert(T,-4.444444444444444)
    a43 = convert(T,2.0)
    a51 = convert(T,-1.409625)
    a52 = convert(T,2.1375)
    a53 = convert(T,-0.2295)
    a54 = convert(T,0.401625)
    a61 = convert(T,-4.0)
    a62 = convert(T,5.0)
    a63 = convert(T,0.0)
    a64 = convert(T,0.0)
    a65 = convert(T,0.0)
    c2 = convert(T2,0.1)
    c3 = convert(T2,0.3333333333333333)
    c4 = convert(T2,0.6666666666666667)
    c5 = convert(T2,0.9)
    c6 = convert(T2,1)
    b1 = convert(T,0.1064814814814815)
    b3 = convert(T,0.4632352941176471)
    b4 = convert(T,0.1607142857142857)
    b5 = convert(T,0.3112356053532524)
    b6 = convert(T,-.04166666666666667)
    Anas5ConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,c2,c3,c4,c5,c6,b1,b3,b4,b5,b6)
end

function Anas5ConstantCache(T, T2)
    a21 = convert(T,1//10)
    a31 = convert(T,-2//9)
    a32 = convert(T,5//9)
    a41 = convert(T,28//9)
    a42 = convert(T,-40//9)
    a43 = convert(T,2//1)
    a51 = convert(T,-11277//8000)
    a52 = convert(T,171//80)
    a53 = convert(T,-459//2000)
    a54 = convert(T,3213//8000)
    a61 = convert(T,-4//1)
    a62 = convert(T,5//1)
    a63 = convert(T,0//1)
    a64 = convert(T,0//1)
    a65 = convert(T,0//1)
    c2 = convert(T2,1//10)
    c3 = convert(T2,1//3)
    c4 = convert(T2,2//3)
    c5 = convert(T2,9//10)
    c6 = convert(T2,1//1)
    b1 = convert(T,23//216)
    b3 = convert(T,63//136)
    b4 = convert(T,9//56)
    b5 = convert(T,1000//3213)
    b6 = convert(T,-1//24)
    Anas5ConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,c2,c3,c4,c5,c6,b1,b3,b4,b5,b6)
end

struct MSRK5ConstantCache{T,T1} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  a41::T
  a43::T
  a51::T
  a53::T
  a54::T
  a61::T
  a63::T
  a64::T
  a65::T
  a71::T
  a73::T
  a74::T
  a75::T
  a76::T
  a81::T
  a83::T
  a84::T
  a85::T
  a86::T
  a87::T

  b1::T # == a_{9i} ∀ i=1:8
  b4::T
  b5::T
  b6::T
  b7::T
  b8::T

  c2::T1
  c3::T1
  c4::T1
  c5::T1
  c6::T1
  c7::T1
  c8::T1

end

# Use rational numbers for testing. Define another function that defines the tab using floats.
function MSRK5ConstantCache(T::Type, T1::Type)
  a21 = T(4//45)
  a31 = T(1//30)
  a32 = T(1//10)
  a41 = T(1//20)
  a43 = T(3//20)
  a51 = T(1//2)
  a53 = T(-15//8)
  a54 = T(15//8)
  a61 = T(-11//135)
  a63 = T(23//45)
  a64 = T(-2//27)
  a65 = T(8//45)
  a71 = T(5//108)
  a73 = T(35//72)
  a74 = T(-59//216)
  a75 = T(-25//27)
  a76 = T(3//2)
  a81 = T(31//128)
  a83 = T(-7563//4480)
  a84 = T(233//112)
  a85 = T(3461//2240)
  a86 = T(-765//448)
  a87 = T(153//320)
  b1 = T(29//456)
  b4 = T(11//38)
  b5 = T(2//27)
  b6 = T(11//40)
  b7 = T(4//19)
  b8 = T(224//2565)
  c2 = T1(4//45)
  c3 = T1(2//15)
  c4 = T1(1//5)
  c5 = T1(1//2)
  c6 = T1(8//15)
  c7 = T1(5//6)
  c8 = T1(19//20)

  MSRK5ConstantCache(a21,a31,a32,a41,a43,a51,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,b1,b4,b5,b6,b7,b8,c2,c3,c4,c5,c6,c7,c8)
end


text:
import Voronoi
import Covering
import Draw

using Random
using Printf

#using Algencan

function unpack(x, n)
   r = x[1]

   points = Array{Tuple{Real, Real}}([])

   for i in 1:n
       append!(points, [(x[2i], x[2i + 1])])
   end

   return r, points
end

#function randf(start, finish, n)
#   v = rand(n)
#   return map(x -> start + x*(finish-start), v)
#end

#Em todos os testes, vamos considerar o retangulo de vertices (0,0); (3,0); (3,1); (0,1).

confs = Array{Real}[]

#########
# TESTE 1
#########

push!(confs, [0.3,
1.11, 0.4,
1.44, 0.25,
1.4, 0.77])

#########
# TESTE 2
#########

push!(confs, [0.3,
1.11, 0.4,
1.44, 0.25,
1.4, 0.77,
0.01, 0.01,
2.0, 0.8,
0.2, 0.73,
2.65, 0.51001])

#########
# TESTE 3
#########

push!(confs, [0.3,
1.11, 0.698,
0.887, 0.554,
1.4, 0.77,
1.554, 0.1,
1.8, 0.7,
0.99, 0.177,
1.805, 0.44])

#########
# TESTE 4
#########

push!(confs, [0.2,
1.11, 0.698,
1.05, 0.45,
1.322, 0.66,
1.39, 0.209,
2.99, 0.99,
1.12, 0.22,
1.505, 0.44,
0.01, 0.99,
0.01, 0.01,
2.99, 0.01])

#########
# TESTE 5
#########

push!(confs, [0.2777,
0.01, 0.499,
2.99, 0.499,
2.99, 0.99,
1.5, 0.99,
1.5, 0.01,
0.01, 0.99,
0.01, 0.01,
2.99, 0.01])

#########
# TESTE 6
#########

push!(confs, [0.3,
0.28, 0.8,
0.28, 0.6,
0.28, 0.4,
0.28, 0.2,
0.8, 0.8,
0.8, 0.6,
0.8, 0.4,
0.8, 0.2,
1.32, 0.8,
1.32, 0.6,
1.32, 0.4,
1.32, 0.2,
1.84, 0.8,
1.84, 0.6,
1.84, 0.4,
1.84, 0.2,
2.36, 0.8,
2.36, 0.6,
2.36, 0.4,
2.36, 0.2,
2.84, 0.8,
2.84, 0.6,
2.84, 0.4,
2.84, 0.2])

#########
# TESTE 7
#########

push!(confs, [0.5,
0.0, 0.5,
3.0, 0.5])

#########
# TESTE 8
#########

push!(confs, [0.5,
0.5, 0.5,
2.5, 0.5,
1.5, 0.5])

#########
# TESTE 9
#########

push!(confs, [0.5,
0.0, 1.0,
0.0, 0.0,
3.0, 0.0,
3.0, 1.0])

#########
# TESTE 10
#########

push!(confs, [0.5,
0.0, 1.0,
0.0, 0.0,
3.0, 0.0,
3.0, 1.0,
1.0, 0.5,
2.0, 0.5])

function gradientDescent()
   WIDTH = 3
   HEIGHT = 1

   command = nothing
   for (i, x) in enumerate(confs)
      n = Int((length(x) - 1)/2)
      r, points = unpack(x, n)

      println("=============================")
      println("CASO ", i)
      println("r: ", r)
      println("pontos:")
      for p in points
          println("(", p[1], ", ", p[2], ")")
      end
      println()

      iterations = 1
      fortunes = 1

      time = @elapsed V, f, ∇f, ξ = Voronoi.Optimization.init(points, WIDTH, HEIGHT)

      μ = 0.5
      α = 0.01
      ϵ = 10^-4
      λᵣ = 0.5
      λₛ = 1

      ρ = 0.1

      times = time

      V = Voronoi.Fortune.compute(points, WIDTH, HEIGHT)
      Voronoi.Intersect.intersect(V, Voronoi.Intersect.Rectangle(WIDTH, HEIGHT))
      P = Covering.voronoiDiagramToPartition(V, r)
      #@printf("f = %.10g\t\t|∇f| = %.10g\t\ttime: %s seconds\n", f, ξ, time)
      area, gᵣ, gₛ = Covering.areaAndGradient(P)
      #println("área coberta: ", area)
      @printf("area coberta= %.8g\ng = %.10g\n∂g/∂r = %.8g\n", area, (WIDTH*HEIGHT - area), gᵣ)
      for i in 1:length(gₛ)
          x, y = gₛ[i]
          @printf("∂g/∂c_%d = (%.8g, %.8g)\n", i, x, y)
      end
      @printf("\n")

      Draw.init(WIDTH, HEIGHT)
      Draw.voronoiDiagram(V)
      Draw.coveringPartition(P)
      Draw.commit()

      times += 0#time
      fortunes += 0#forts
      iterations += 1

      if command != "a"
          #println("Press Return for a step or enter \"a\" to animate until the end.")
          command = readline(stdin)
      else
          sleep(0.00001) # make sure stuff is drawn
      end
      println("\n")

      #println("n: ", n, "\t\titer: ", iterations, "\t\tfortunes: ", fortunes, "\t\ttime: ", times)
  end
end

gradientDescent()

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
s = open("input.txt") do file
    readlines(file);
end

t = open("test.txt") do file
    readlines(file);
end

function ParseRules(lines)
    rules = Dict();

    for line in lines
        if length(line) < 2
            break;
        end

        parts = split(line, ": ");
        id = parse(Int, parts[1]);
        options = []

        if occursin("|", parts[2])
            ruleBits = split(parts[2]," | ")
        else
            ruleBits = [parts[2]]
        end

        m = match(r"\"([ab])\"", line);
        if m != nothing
            # This is a base rule
            push!(options, String(m.captures[1]));
        else
            println(ruleBits);
            options = map(x->map(y->parse(Int,y), split(x)), ruleBits);
        end

        rules[id] = options;
    end

    return rules;
end

function SolveRule!(id, rules, solved)
    if haskey(solved, id)
        return;
    end

    # println(id);
    r = rules[id];
    # println(r);
    if length(r) == 1 && typeof(r[1]) == String
        solved[id] = r[1];
        return;
    else
        for ruleOpt in r
            for rule in ruleOpt
                SolveRule!(rule, rules, solved);
            end
        end

        trueForms = map(rule -> map(x->solved[x], rule), r);
        # println(trueForms)
        if length(r) == 2
            part1 = *(trueForms[1]...);
            part2 = *(trueForms[2]...);
            newRule = "("*part1*"|"*part2*")";
        else
            newRule = *(trueForms[1]...);
        end

        solved[id] = newRule;
        println("$(id):$(r) => $(newRule)");
    end
end

function SolveAllRules(rules)
    solved = Dict();
    SolveRule!(0, rules, solved);

    println(length(solved));
    return (Regex("^"*solved[0]*"\$"), solved);
end

rules = ParseRules(s)
(realRule, solved) = SolveAllRules(rules)

messages = s[length(rules)+2:end]
uncorrupt = GetCount(messages, realRule)

# test
# rules_t = ParseRules(t)
# realRule_t = GetFinalRule(rules_t)
# realRule_t = Regex("^"*rules_t[0]*"b\$");
#
# messages_t = t[length(rules_t)+2:end]
# uncorrupt_t = map(x->occursin(realRule_t, x), messages_t)


text:
@symbol_func function K_TB(cur_reactor::AbstractReactor)
  cur_K = 4.027e-2

  cur_K *= K_n(cur_reactor)

  cur_K *= cur_reactor.epsilon

  cur_K /= cur_reactor.max_beta_N

  cur_K *= 1 + cur_reactor.f_D

  cur_K *= 1 + cur_reactor.nu_n

  cur_K *= 1 + cur_reactor.nu_T

  cur_K /= 1 + cur_reactor.nu_n + cur_reactor.nu_T

  cur_K
end


text:
function stress_centrality(g::AbstractGraph,
    vs::AbstractVector=vertices(g))::Vector{Int}

    n_v = nv(g)
    k = length(vs)
    isdir = is_directed(g)

    # Parallel reduction
    stress = @distributed (+) for s in vs
        temp_stress = zeros(Int, n_v)
        if degree(g, s) > 0  # this might be 1?
            state = LightGraphs.dijkstra_shortest_paths(g, s; allpaths=true, trackvertices=true)
            LightGraphs._stress_accumulate_basic!(temp_stress, state, g, s)
        end
        temp_stress
    end
    return stress
end

stress_centrality(g::AbstractGraph, k::Integer) =
    stress_centrality(g, sample(vertices(g), k))


text:
module TestReduce
    using NullableArrays
    using Base.Test

    srand(1)
    f(x) = 5 * x
    f{T<:Number}(x::Nullable{T}) = Nullable(5 * x.value, x.isnull)

    for N in (10, 2050)
        A = rand(N)
        M = rand(Bool, N)
        i = rand(1:N)
        M[i] = true
        j = rand(1:N)
        while j == i
            j = rand(1:N)
        end
        M[j] = false
        X = NullableArray(A)
        Y = NullableArray(A, M)
        B = A[find(x->!x, M)]

        @test isequal(mapreduce(f, +, X), Nullable(mapreduce(f, +, X.values)))
        @test isequal(mapreduce(f, +, Y), Nullable{Float64}())
        v = mapreduce(f, +, Y, skipnull=true)
        @test_approx_eq v.value mapreduce(f, +, B)
        @test !v.isnull

        @test isequal(reduce(+, X), Nullable(reduce(+, X.values)))
        @test isequal(reduce(+, Y), Nullable{Float64}())
        v = reduce(+, Y, skipnull=true)
        @test_approx_eq v.value reduce(+, B)
        @test !v.isnull

        for method in (
            sum,
            prod,
            minimum,
            maximum,
        )
            @test isequal(method(X), Nullable(method(A)))
            @test isequal(method(f, X), Nullable(method(f, A)))
            @test isequal(method(Y), Nullable{Float64}())
            v = method(Y, skipnull=true)
            @test_approx_eq v.value method(B)
            @test !v.isnull
            @test isequal(method(f, Y), Nullable{Float64}())
            v = method(f, Y, skipnull=true)
            @test_approx_eq v.value method(f, B)
            @test !v.isnull
        end

        for method in (
            sumabs,
            sumabs2,
        )
            @test isequal(method(X), Nullable(method(A)))
            @test isequal(method(Y), Nullable{Float64}())
            v = method(Y, skipnull=true)
            @test_approx_eq v.value method(B)
            @test v.isnull == false
        end

        H = rand(Bool, N)
        G = H[find(x->!x, M)]
        U = NullableArray(H)
        V = NullableArray(H, M)

        for op in (
            &,
            |,
        )
            @test isequal(reduce(op, U),
                          Nullable(reduce(op, H)))
            @test isequal(reduce(op, U, skipnull=true),
                          Nullable(reduce(op, H)))
            @test isequal(reduce(op, V, skipnull=true),
                          Nullable(reduce(op, G)))
        end
    end
end


text:
@defcomp totalfactorproductivity begin

	GA		= Variable(index=[time])	#Growth rate of productivity
    AL		= Variable(index=[time])	#Level of total factor productivity

    a0		= Parameter()				#Initial level of total factor productivity
	ga0		= Parameter()				#Initial growth rate for TFP per 5 years
	dela	= Parameter()				#Decline rate of TFP per 5 years

    function run_timestep(p, v, d, t)
		#Define function for GA
        if is_first(t)
            v.GA[t] = p.ga0
        else
            v.GA[t] = v.GA[t - 1] * exp(-p.dela * 5)
        end

		#Define function for AL
		if is_first(t)
			v.AL[t] = p.a0
		else
			v.AL[t] = v.AL[t-1]/(1 - v.GA[t-1])
		end
    end
end