text:
using StaticArrays

# Exports from Vector2D implementation
export Vector2D, Vector2D!, equal, add!, multiply!

"""
Wrapper around StaticArrays.MVector{2,T} for simpler use in
physics applications.
"""
mutable struct Vector2D{T <: Real} <: PhysicalVector
    vec::MVector{2,T}
end

#= ========================================================================================= =#
#  Struct constructors
#= ========================================================================================= =#
"""
    Vector2D(arr::AbstractArray{Real, 1})

Copy constructor from an `AbstractVector`

# Examples

```jldoctest
julia> Vector2D([1., 2])
Vector2D{Float64}([1.0, 2.0])
```
"""
Vector2D(arr::AbstractArray{T,1}) where {T <: Real} = Vector2D{T}( MVector{2,T}(arr[1], arr[2]) )

"""
    Vector2D(vec::Vector2D)

Copy constructor from an `Vector2D`

# Examples

```jldoctest
julia> Vector2D( Vector2D([1., 2]) )
Vector2D{Float64}([1.0, 2.0])
```
"""
Vector2D(vec::Vector2D{T}) where {T <: Real} = Vector2D( vec.vec )

"""
    Vector2D(x::T, y::T) where {T <: Real}

Copy constructor from a pair x and y _explicitly_ of the same type `T <: Real`

# Examples

```jldoctest
julia> Vector2D(1., 2.)
Vector2D{Float64}([1.0, 2.0])
```
"""
Vector2D(x::T, y::T) where {T <: Real} = Vector2D{T}( MVector{2,T}(x,y) )

"""
    Vector2D(x::Real, y::Real)

Copy constructor from a pair x and y of potentially different `Real` types.

# Additional Information
This function makes a call to `Vector2D(::Tuple{S, T})` where `{S <: Real, T <: Real}`.

# Examples

```jldoctest
julia> Vector2D(1, 2.)
Vector2D{Float64}([1.0, 2.0])
```
"""
Vector2D(x::Real, y::Real) = Vector2D((x, y))

"""
    Vector2D(tup::Tuple{S, T}) where {S <: Real, T <: Real}

Copy constructor from a tuple of reals `tup === (x::Real, y::Real)`.

# Additional Information
In the case that the elements are of different types, they are first
passed through `promote(x, y)` and the common-type values are stored
in the output.

# Examples

```jldoctest
julia> Vector2D((1, 2.))
Vector2D{Float64}([1.0, 2.0])
```
"""
function Vector2D(tup::Tuple{S, T}) where {S <: Real, T <: Real}
    vals = promote(tup...)
    return Vector2D{eltype(vals)}( MVector{2, eltype(vals)}(vals) )
end

"""
    Vector2D!(old_vec::Vector2D{P}, new_components::Tuple{S, T}) where {P <: Real, S <: Real, T <: Real}

In-place copy from a tuple of reals `new_components === (x::Real, y::Real)`.

# Additional Information
In the case that the elements are of different types, they are first
passed through `promote(x, y)` and the common-type values are stored
in the output.

# Examples

```jldoctest
julia> A = Vector2D((1, 2.))
Vector2D{Float64}([1.0, 2.0])

julia> Vector2D!(A, (2, 4))

julia> A
Vector2D{Float64}([2.0, 4.0])
```
"""
function Vector2D!( old_vec::Vector2D{P}, new_components::Tuple{S, T} ) where {P <: Real, S <: Real, T <: Real}
    vals = convert( Tuple{eltype(old_vec.vec), eltype(old_vec.vec)}, new_components)
    broadcast!( (x, y) -> y, old_vec.vec, old_vec.vec, vals )
    return nothing
end

#= ========================================================================================= =#
#  Interface requirement definitions.
#= ========================================================================================= =#
"""
    equal(A::Vector2D, B::Vector2D)

Vector2D _equality_ based on component-wise equality.

# Additional Information
Use this function for equality instead of `==` which defaults to object _egality_ since
`Vector2D` is mutable.

# Examples

```jldoctest
julia> A = Vector2D(1., 2)
Vector2D{Float64}([1.0, 2.0])

julia> B = Vector2D(1, 2.)
Vector2D{Float64}([1.0, 2.0])

julia> equal(A, B)
true
```
"""
equal(A::Vector2D, B::Vector2D) = prod(A.vec .== B.vec)

"""
    +(A::Vector2D, B::Vector2D) -> Vector2D

Addition for `Vector2D`. Includes a copy constructor.

```jldoctest
julia> A = Vector2D(1., 2)
Vector2D{Float64}([1.0, 2.0])

julia> B = Vector2D(1, 2.)
Vector2D{Float64}([1.0, 2.0])

julia> A + B
Vector2D{Float64}([2.0, 4.0])
```
"""
Base.:+(A::Vector2D, B::Vector2D) = Vector2D( broadcast(+, A.vec, B.vec) )

"""
    *(λ, A::Vector2D) -> Vector2D

Multiplication by a scalar `λ` for `Vector2D`. Includes a copy constructor.

```jldoctest
julia> A = Vector2D(1., 2)
Vector2D{Float64}([1.0, 2.0])

julia> λ = 2.
2.0

julia> λ * A
Vector2D{Float64}([2.0, 4.0])
```
"""
Base.:*(λ, A::Vector2D) = Vector2D( broadcast(*, A.vec, λ) )

"""
    ⋅(A::Vector2D, B::Vector2D)

Scalar product, i.e. dot product, for `Vector2D`.

```jldoctest
julia> A = Vector2D(1., 2)
Vector2D{Float64}([1.0, 2.0])

julia> B = Vector2D(1, 2.)
Vector2D{Float64}([1.0, 2.0])

julia> A ⋅ B
5.0
```
"""
function ⋅(A::Vector2D, B::Vector2D)
    ans = zero(eltype(A.vec))
    for idx in eachindex(A.vec, B.vec)
        ans += A.vec[idx] * B.vec[idx]
    end
    return ans
end

"""
    add!(A::Vector2D, B::Vector2D) -> nothing

In-place addition for `Vector2D`. The first argument, `A`, is modified.

```jldoctest
julia> A = Vector2D(1., 2)
Vector2D{Float64}([1.0, 2.0])

julia> B = Vector2D(1, 2.)
Vector2D{Float64}([1.0, 2.0])

julia> add!(A, B)

julia> A
Vector2D{Float64}([2.0, 4.0])

julia> B
Vector2D{Float64}([1.0, 2.0])
```
"""
function add!(A::Vector2D, B::Vector2D)
    broadcast!(+, A.vec, A.vec, B.vec)
    return nothing
end

"""
    subtract!(A::Vector2D, B::Vector2D) -> nothing

In-place subtraction for `Vector2D`. The first argument, `A`, is modified.

```jldoctest
julia> A = Vector2D(1., 2)
Vector2D{Float64}([1.0, 2.0])

julia> B = Vector2D(1, 2.)
Vector2D{Float64}([1.0, 2.0])

julia> subtract!(A, B)

julia> A
Vector2D{Float64}([0.0, 0.0])

julia> B
Vector2D{Float64}([1.0, 2.0])
```
"""
function subtract!(A::Vector2D, B::Vector2D)
    broadcast!(-, A.vec, A.vec, B.vec)
    return nothing
end

"""
    multiply!(A::Vector2D, λ) -> nothing

In-place multiplication by a scalar for `Vector2D`. The first argument, `A`, is modified.

```jldoctest
julia> A = Vector2D(1., 2)
Vector2D{Float64}([1.0, 2.0])

julia> λ = 2.
2.0

julia> multiply!(A, λ)

julia> A
Vector2D{Float64}([2.0, 4.0])

julia> λ
2.0
```
"""
function multiply!(A::Vector2D, λ)
    broadcast!(*, A.vec, A.vec, λ)
    return nothing
end

#= ========================================================================================= =#
#  Includes for the Vector2D struct.
#  Each file **must** contain it's own export list.
#= ========================================================================================= =#
include("PeriodicBoundaryConditionHelpers.jl")


text:
using Compat

struct DBM
  filename :: String
  handle :: Ptr{Cvoid}
  DBM(filename::String, flag::String="c") = new(filename, gdbm_open(filename, flag))
end

Base.close(dbm::DBM) = gdbm_close(dbm.handle)

Base.show(io::IO, dbm::DBM) = print(io, "DBM(<$(dbm.filename)>)")

struct Datum
  dptr :: Ptr{UInt8}
  dsize :: Int32
  Datum(str::String) = new(Base.unsafe_convert(Ptr{UInt8}, str), length(str))
  Datum(bytes::Array{UInt8,1}) = new(Base.unsafe_convert(Ptr{UInt8}, bytes), length(bytes))
end

const OPENFLAGS = Dict("r"=>0, "w"=>1, "c"=>2, "n"=>3)

const STOREFLAGS = Dict("r"=>1, "i"=>0)

ext: jl


text:
function no_action(rule, value, first, last, children)
  return Node(rule.name, value, first, last, children, typeof(rule))
end

function liftchild(rule, value, first, last, children)  # = default_action for OR-rules
  if length(children) == 0
    return nothing
  elseif length(children) == 1
    return children[1]
  else
    error("More than one child! rule: $rule; children: $children")
  end
end

function childrenarray(rule, value, first, last, children)
  return children
end

function nodevalue(rule, value, first, last, children)
  return value
end


text:
cd(joinpath(Pkg.dir("TWPBVP"), "deps"))
pic = @windows ? "" : "-fPIC"
run(`gfortran -m$WORD_SIZE -ffixed-form $pic -shared -O2 -o libtwpbvpc.so twpbvpc.f`)
# -fdefault-integer-8
# -fdefault-real-8
# -fdefault-integer-8


text:
@compat abstract type OrdinaryDiffEqCache <: DECache end
@compat abstract type OrdinaryDiffEqConstantCache <: OrdinaryDiffEqCache end
@compat abstract type OrdinaryDiffEqMutableCache <: OrdinaryDiffEqCache end
immutable ODEEmptyCache <: OrdinaryDiffEqConstantCache end
immutable ODEChunkCache{CS} <: OrdinaryDiffEqConstantCache end

type CompositeCache{T,F} <: OrdinaryDiffEqCache
  caches::T
  choice_function::F
  current::Int
end

function alg_cache{T}(alg::CompositeAlgorithm,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{T}})
  caches = map((x)->alg_cache(x,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,Val{T}),alg.algs)
  CompositeCache(caches,alg.choice_function,1)
end

alg_cache{F}(alg::OrdinaryDiffEqAlgorithm,prob,callback::F) = ODEEmptyCache()

immutable DiscreteCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  du::rateType
end

u_cache(c::DiscreteCache) = ()
du_cache(c::DiscreteCache) = (c.du)

function alg_cache(alg::Discrete,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  DiscreteCache(u,uprev,rate_prototype)
end

immutable DiscreteConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::Discrete,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = DiscreteConstantCache()

immutable EulerCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
end

function alg_cache(alg::Euler,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  EulerCache(u,uprev,similar(u),zeros(rate_prototype),zeros(rate_prototype))
end

u_cache(c::EulerCache) = ()
du_cache(c::EulerCache) = (c.k,c.fsalfirst)

immutable EulerConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::Euler,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = EulerConstantCache()

immutable MidpointCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  du::rateType
  tmp::uType
  fsalfirst::rateType
end

u_cache(c::MidpointCache) = ()
du_cache(c::MidpointCache) = (c.k,c.du,c.fsalfirst)

immutable MidpointConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::Midpoint,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  tmp = similar(u)
  k = zeros(rate_prototype)
  du = zeros(rate_prototype)
  fsalfirst = zeros(rate_prototype)
  MidpointCache(u,uprev,k,du,tmp,fsalfirst)
end

alg_cache(alg::Midpoint,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = MidpointConstantCache()

immutable RK4Cache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k₂::rateType
  k₃::rateType
  k₄::rateType
  k::rateType
  tmp::uType
end

u_cache(c::RK4Cache) = ()
du_cache(c::RK4Cache) = (c.fsalfirst,c.k₂,c.k₃,c.k₄,c.k)

immutable RK4ConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::RK4,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  k₁ = zeros(rate_prototype)
  k₂ = zeros(rate_prototype)
  k₃ = zeros(rate_prototype)
  k₄ = zeros(rate_prototype)
  k  = zeros(rate_prototype)
  tmp = similar(u)
  RK4Cache(u,uprev,k₁,k₂,k₃,k₄,k,tmp)
end

alg_cache(alg::RK4,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = RK4ConstantCache()

immutable BS3Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::BS3Cache) = (c.atmp,c.utilde)
du_cache(c::BS3Cache) = (c.fsalfirst,c.k2,c.k3,c.k4)

function alg_cache(alg::BS3,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  tab = BS3ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
  k1 = zeros(rate_prototype)
  k2 = zeros(rate_prototype)
  k3 = zeros(rate_prototype)
  k4 = zeros(rate_prototype)
  utilde = similar(u,indices(u))
  atmp = similar(u,uEltypeNoUnits,indices(u))
  tmp = similar(u)
  BS3Cache(u,uprev,k1,k2,k3,k4,utilde,tmp,atmp,tab)
end

alg_cache(alg::BS3,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = BS3ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))

immutable BS5Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  utilde::uArrayType
  uhat::uType
  tmp::uType
  atmp::uEltypeNoUnits
  atmptilde::uEltypeNoUnits
  tab::TabType
end

u_cache(c::BS5Cache) = (c.utilde,c.uhat,c.atmp,c.atmptilde)
du_cache(c::BS5Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8)

function alg_cache(alg::BS5,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  tab = BS5ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
  k1 = zeros(rate_prototype)
  k2 = zeros(rate_prototype)
  k3 = zeros(rate_prototype)
  k4 = zeros(rate_prototype)
  k5 = zeros(rate_prototype)
  k6 = zeros(rate_prototype)
  k7 = zeros(rate_prototype)
  k8 = zeros(rate_prototype)
  utilde = similar(u,indices(u)); uhat = similar(u,indices(u))
  atmp = similar(u,uEltypeNoUnits,indices(u))
  tmp = similar(u)
  atmptilde = similar(u,uEltypeNoUnits,indices(u))
  BS5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,utilde,uhat,tmp,atmp,atmptilde,tab)
end

alg_cache(alg::BS5,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = BS5ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))

immutable Tsit5Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::Tsit5Cache) = (c.utilde,c.atmp)
du_cache(c::Tsit5Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7)

function alg_cache(alg::Tsit5,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  tab = Tsit5ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
  k1 = zeros(rate_prototype)
  k2 = zeros(rate_prototype)
  k3 = zeros(rate_prototype)
  k4 = zeros(rate_prototype)
  k5 = zeros(rate_prototype)
  k6 = zeros(rate_prototype)
  k7 = zeros(rate_prototype)
  utilde = similar(u,indices(u))
  atmp = similar(u,uEltypeNoUnits,indices(u))
  tmp = similar(u)
  atmptilde = similar(u,uEltypeNoUnits,indices(u))
  Tsit5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp,tab)
end

alg_cache(alg::Tsit5,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = Tsit5ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))

immutable DP5Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  dense_tmp3::rateType
  dense_tmp4::rateType
  update::rateType
  bspl::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::DP5Cache) = (c.utilde,c.atmp)
du_cache(c::DP5Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.dense_tmp3,c.dense_tmp4,c.update,c.bspl)

function alg_cache(alg::DP5,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  k1 = zeros(rate_prototype)
  k2 = zeros(rate_prototype)
  k3 = zeros(rate_prototype)
  k4 = zeros(rate_prototype)
  k5 = zeros(rate_prototype)
  k6 = zeros(rate_prototype)
  k7 = zeros(rate_prototype)
  dense_tmp3 = zeros(rate_prototype)
  dense_tmp4 = zeros(rate_prototype)
  update = zeros(rate_prototype)
  bspl = zeros(rate_prototype)
  utilde = similar(u,indices(u))
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits,indices(u))
  tab = DP5ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
  cache = DP5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp,tab)
  cache
end

alg_cache(alg::DP5,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = DP5ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))

immutable DP5ThreadedCache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  dense_tmp3::rateType
  dense_tmp4::rateType
  update::rateType
  bspl::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::DP5ThreadedCache) = (c.utilde,c.atmp)
du_cache(c::DP5ThreadedCache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.dense_tmp3,c.dense_tmp4,c.update,c.bspl)

function alg_cache(alg::DP5Threaded,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  k1 = zeros(rate_prototype)
  k2 = zeros(rate_prototype)
  k3 = zeros(rate_prototype)
  k4 = zeros(rate_prototype)
  k5 = zeros(rate_prototype)
  k6 = zeros(rate_prototype)
  k7 = zeros(rate_prototype)
  dense_tmp3 = zeros(rate_prototype)
  dense_tmp4 = zeros(rate_prototype)
  update = zeros(rate_prototype)
  bspl = zeros(rate_prototype)
  utilde = similar(u,indices(u))
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits,indices(u))
  tab = DP5ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
  DP5ThreadedCache(u,uprev,k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp,tab)
end

immutable Vern6Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::Vern6Cache) = (c.utilde,c.atmp)
du_cache(c::Vern6Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9)

function alg_cache(alg::Vern6,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  tab = Vern6ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
  k1 = zeros(rate_prototype)
  k2 = zeros(rate_prototype); k3 = zeros(rate_prototype); k4 = zeros(rate_prototype);
  k5 = zeros(rate_prototype); k6 = zeros(rate_prototype); k7 = zeros(rate_prototype);
  k8 = zeros(rate_prototype); k9 = zeros(rate_prototype);
  utilde = similar(u,indices(u)); tmp = similar(u); atmp = similar(u,uEltypeNoUnits,indices(u));
  Vern6Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,utilde,tmp,atmp,tab)
end

alg_cache(alg::Vern6,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = Vern6ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))

immutable Vern7Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  utilde::uArrayType
  update::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::Vern7Cache) = (c.utilde,c.update,c.atmp)
du_cache(c::Vern7Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9,c.k10)

function alg_cache(alg::Vern7,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  tab = Vern7ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
  k1 = zeros(rate_prototype); k2 = zeros(rate_prototype); k3 = zeros(rate_prototype); k4 = zeros(rate_prototype);
  k5 = zeros(rate_prototype); k6 = zeros(rate_prototype); k7 = zeros(rate_prototype); k8 = zeros(rate_prototype);
  k9 = zeros(rate_prototype); k10 = zeros(rate_prototype); utilde = similar(u,indices(u)); update = similar(u,indices(u))
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits,indices(u))
  Vern7Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,utilde,update,tmp,atmp,tab)
end

alg_cache(alg::Vern7,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = Vern7ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))

immutable Vern8Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  k11::rateType
  k12::rateType
  k13::rateType
  utilde::uArrayType
  update::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::Vern8Cache) = (c.utilde,c.update,c.atmp)
du_cache(c::Vern8Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9,c.k10,c.k11,c.k12,c.k13)

function alg_cache(alg::Vern8,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  tab = Vern8ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
  k1 = zeros(rate_prototype); k2 = zeros(rate_prototype); k3 = zeros(rate_prototype);
  k4 = zeros(rate_prototype);
  k5 = zeros(rate_prototype); k6 = zeros(rate_prototype); k7 = zeros(rate_prototype);
  k8 = zeros(rate_prototype); tmp = similar(u)
  k9 = zeros(rate_prototype); k10 = zeros(rate_prototype); k11 = zeros(rate_prototype);
  k12 = zeros(rate_prototype); k13 = zeros(rate_prototype)
  utilde = similar(u,indices(u)); update = similar(u,indices(u));
  atmp = similar(u,uEltypeNoUnits,indices(u))
  Vern8Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,utilde,update,tmp,atmp,tab)
end

alg_cache(alg::Vern8,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = Vern8ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))

immutable Vern9Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  k11::rateType
  k12::rateType
  k13::rateType
  k14::rateType
  k15::rateType
  k16::rateType
  utilde::uArrayType
  update::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::Vern9Cache) = (c.utilde,c.update,c.atmp)
du_cache(c::Vern9Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9,c.k10,c.k11,c.k12,c.k13,c.k14,c.k15,c.k16)

function alg_cache(alg::Vern9,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  tab = Vern9ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
  k1 = zeros(rate_prototype); k2 = zeros(rate_prototype);k3 = zeros(rate_prototype);
  k4 = zeros(rate_prototype);
  k5 = zeros(rate_prototype); k6 = zeros(rate_prototype);k7 = zeros(rate_prototype);
  k8 = zeros(rate_prototype);
  k9 = zeros(rate_prototype); k10 = zeros(rate_prototype); k11 = zeros(rate_prototype);
  k12 = zeros(rate_prototype); update = similar(u,indices(u))
  k13 = zeros(rate_prototype); k14 = zeros(rate_prototype); k15 = zeros(rate_prototype);
  k16 =zeros(rate_prototype);
  utilde = similar(u,indices(u)); tmp = similar(u); atmp = similar(u,uEltypeNoUnits,indices(u));
  Vern9Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,utilde,update,tmp,atmp,tab)
end

alg_cache(alg::Vern9,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = Vern9ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))

immutable TanYam7Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  k::rateType
  tab::TabType
end

u_cache(c::TanYam7Cache) = (c.utilde,c.atmp)
du_cache(c::TanYam7Cache) = (c.fsalfirst,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9,c.k10,c.k)

function alg_cache(alg::TanYam7,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  tab = TanYam7ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
  k1 = zeros(rate_prototype); k2 = zeros(rate_prototype) ; k3 = zeros(rate_prototype); k4 = zeros(rate_prototype)
  k5 = zeros(rate_prototype); k6 = zeros(rate_prototype) ; k7 = zeros(rate_prototype); k8 = zeros(rate_prototype)
  k9 = zeros(rate_prototype); k10= zeros(rate_prototype) ;
  utilde = similar(u,indices(u)); tmp = similar(u); atmp = similar(u,uEltypeNoUnits,indices(u)); k = zeros(rate_prototype)
  TanYam7Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,utilde,tmp,atmp,k,tab)
end

alg_cache(alg::TanYam7,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = TanYam7ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))

immutable DP8Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  k11::rateType
  k12::rateType
  k13::rateType
  k14::rateType
  k15::rateType
  k16::rateType
  kupdate::rateType
  update::uArrayType
  udiff::rateType
  bspl::rateType
  dense_tmp3::rateType
  dense_tmp4::rateType
  dense_tmp5::rateType
  dense_tmp6::rateType
  dense_tmp7::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  atmp2::uEltypeNoUnits
  tab::TabType
end

u_cache(c::DP8Cache) = (c.update,c.utilde,c.atmp,c.atmp2)
du_cache(c::DP8Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9,c.k10,c.k11,c.k12,c.k13,c.k14,c.k15,c.k16,c.kupdate,c.udiff,c.bspl,c.dense_tmp3,c.dense_tmp4,c.dense_tmp5,c.dense_tmp6,c.dense_tmp7)

function alg_cache(alg::DP8,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  k1 = zeros(rate_prototype); k2  = zeros(rate_prototype); k3  = zeros(rate_prototype);  k4 = zeros(rate_prototype)
  k5 = zeros(rate_prototype); k6  = zeros(rate_prototype); k7  = zeros(rate_prototype);  k8 = zeros(rate_prototype)
  k9 = zeros(rate_prototype); k10 = zeros(rate_prototype); k11 = zeros(rate_prototype); k12 = zeros(rate_prototype)
  kupdate = zeros(rate_prototype); utilde = similar(u,indices(u));
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits,indices(u)); atmp2 = similar(u,uEltypeNoUnits,indices(u)); update = similar(u,indices(u))
  k13 = zeros(rate_prototype)
  k14 = zeros(rate_prototype)
  k15 = zeros(rate_prototype)
  k16 = zeros(rate_prototype)
  udiff = zeros(rate_prototype)
  bspl = zeros(rate_prototype)
  # dense_tmp1 = udiff
  # dense_tmp2 = bspl
  dense_tmp3 = zeros(rate_prototype)
  dense_tmp4 = zeros(rate_prototype)
  dense_tmp5 = zeros(rate_prototype)
  dense_tmp6 = zeros(rate_prototype)
  dense_tmp7 = zeros(rate_prototype)
  tab = DP8ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
  DP8Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,kupdate,
           update,udiff,bspl,dense_tmp3,dense_tmp4,dense_tmp5,dense_tmp6,dense_tmp7,
           utilde,tmp,atmp,atmp2,tab)
end

alg_cache(alg::DP8,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = DP8ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))

immutable TsitPap8Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  k11::rateType
  k12::rateType
  k13::rateType
  utilde::uArrayType
  update::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  k::rateType
  tab::TabType
end

u_cache(c::TsitPap8Cache) = (c.utilde,c.atmp,c.update)
du_cache(c::TsitPap8Cache) = (c.fsalfirst,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9,c.k10,c.k11,c.k12,c.k13,c.k)

function alg_cache(alg::TsitPap8,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  tab = TsitPap8ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
  k1 = zeros(rate_prototype); k2 = zeros(rate_prototype); k3 = zeros(rate_prototype); k4 = zeros(rate_prototype)
  k5 = zeros(rate_prototype); k6 = zeros(rate_prototype); k7 = zeros(rate_prototype); k8 = zeros(rate_prototype)
  k9 = zeros(rate_prototype); k10 = zeros(rate_prototype); k11 = zeros(rate_prototype); k12 = zeros(rate_prototype)
  k13 = zeros(rate_prototype); update = similar(u,indices(u)); utilde = similar(u,indices(u)); k = zeros(rate_prototype)
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits,indices(u));
  TsitPap8Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,utilde,update,tmp,atmp,k,tab)
end

alg_cache(alg::TsitPap8,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = TsitPap8ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))

immutable ExplicitRKCache{uType,rateType,uEltypeNoUnits,ksEltype,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  utilde::rateType
  uEEst::rateType
  atmp::uEltypeNoUnits
  fsalfirst::ksEltype
  fsallast::ksEltype
  kk::Vector{ksEltype}
  tab::TabType
end

u_cache(c::ExplicitRKCache) = (c.utilde,c.atmp,c.update)
du_cache(c::ExplicitRKCache) = (c.uEEst,c.kk...)

function alg_cache(alg::ExplicitRK,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  kk = Vector{typeof(rate_prototype)}(0)
  for i = 1:alg.tableau.stages
    push!(kk,zeros(rate_prototype))
  end
  fsalfirst = kk[1]
  if isfsal(alg.tableau)
    fsallast = kk[end]
  else
    fsallast = zeros(rate_prototype)
  end
  utilde = zeros(rate_prototype)
  tmp = similar(u)
  atmp = similar(u,uEltypeNoUnits,indices(u))
  uEEst = zeros(rate_prototype)
  tab = ExplicitRKConstantCache(alg.tableau,rate_prototype)
  ExplicitRKCache(u,uprev,tmp,utilde,uEEst,atmp,fsalfirst,fsallast,kk,tab)
end

immutable ExplicitRKConstantCache{MType,VType,KType} <: OrdinaryDiffEqConstantCache
  A::MType
  c::VType
  α::VType
  αEEst::VType
  stages::Int
  kk::KType
end

function ExplicitRKConstantCache(tableau,rate_prototype)
  @unpack A,c,α,αEEst,stages = tableau
  A = A' # Transpose A to column major looping
  kk = Array{typeof(rate_prototype)}(stages) # Not ks since that's for integrator.opts.dense
  ExplicitRKConstantCache(A,c,α,αEEst,stages,kk)
end

alg_cache(alg::ExplicitRK,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = ExplicitRKConstantCache(alg.tableau,rate_prototype)

immutable Feagin10Cache{uType,uEltypeNoUnits,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  k11::rateType
  k12::rateType
  k13::rateType
  k14::rateType
  k15::rateType
  k16::rateType
  k17::rateType
  tmp::uType
  atmp::uEltypeNoUnits
  k::rateType
  tab::TabType
end

u_cache(c::Feagin10Cache) = (c.atmp,)
du_cache(c::Feagin10Cache) = (c.fsalfirst,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9,c.k10,c.k11,c.k12,c.k13,c.k14,c.k15,c.k16,c.k17,c.k)

function alg_cache(alg::Feagin10,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  tab = Feagin10ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
  k1 = zeros(rate_prototype); k2 = zeros(rate_prototype); k3 = zeros(rate_prototype); k4 = zeros(rate_prototype); k5 = zeros(rate_prototype)
  k6 = zeros(rate_prototype); k7 = zeros(rate_prototype); k8 = zeros(rate_prototype); k9 = zeros(rate_prototype); k10 = zeros(rate_prototype)
  k11 = zeros(rate_prototype); k12 = zeros(rate_prototype); k13 = zeros(rate_prototype); k14 = zeros(rate_prototype)
  k15 = zeros(rate_prototype); k16 = zeros(rate_prototype); k17 = zeros(rate_prototype)
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits,indices(u)); k = zeros(rate_prototype)

  Feagin10Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,k17,tmp,atmp,k,tab)
end

alg_cache(alg::Feagin10,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = Feagin10ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))

immutable Feagin12Cache{uType,uEltypeNoUnits,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  k11::rateType
  k12::rateType
  k13::rateType
  k14::rateType
  k15::rateType
  k16::rateType
  k17::rateType
  k18::rateType
  k19::rateType
  k20::rateType
  k21::rateType
  k22::rateType
  k23::rateType
  k24::rateType
  k25::rateType
  tmp::uType
  atmp::uEltypeNoUnits
  k::rateType
  tab::TabType
end

u_cache(c::Feagin12Cache) = (c.atmp,)
du_cache(c::Feagin12Cache) = (c.fsalfirst,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9,c.k10,c.k11,c.k12,c.k13,c.k14,c.k15,c.k16,c.k17,c.k18,c.k19,c.k20,c.k21,c.k22,c.k23,c.k24,c.k25,c.k)

function alg_cache(alg::Feagin12,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  tab = Feagin12ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
  k1 = zeros(rate_prototype); k2 = zeros(rate_prototype); k3 = zeros(rate_prototype); k4 = zeros(rate_prototype); k5 = zeros(rate_prototype)
  k6 = zeros(rate_prototype); k7 = zeros(rate_prototype); k8 = zeros(rate_prototype); k9 = zeros(rate_prototype); k10 = zeros(rate_prototype)
  k11 = zeros(rate_prototype); k12 = zeros(rate_prototype); k13 = zeros(rate_prototype); k14 = zeros(rate_prototype)
  k15 = zeros(rate_prototype); k16 = zeros(rate_prototype); k17 = zeros(rate_prototype); k18 = zeros(rate_prototype)
  k19 = zeros(rate_prototype); k20 = zeros(rate_prototype); k21 = zeros(rate_prototype); k22 = zeros(rate_prototype)
  k23 = zeros(rate_prototype); k24 = zeros(rate_prototype); k25 = zeros(rate_prototype)
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits,indices(u)); k = zeros(rate_prototype)

  Feagin12Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,k17,k18,k19,k20,k21,k22,k23,k24,k25,tmp,atmp,k,tab)
end

alg_cache(alg::Feagin12,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = Feagin12ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))

immutable Feagin14Cache{uType,uEltypeNoUnits,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  k9::rateType
  k10::rateType
  k11::rateType
  k12::rateType
  k13::rateType
  k14::rateType
  k15::rateType
  k16::rateType
  k17::rateType
  k18::rateType
  k19::rateType
  k20::rateType
  k21::rateType
  k22::rateType
  k23::rateType
  k24::rateType
  k25::rateType
  k26::rateType
  k27::rateType
  k28::rateType
  k29::rateType
  k30::rateType
  k31::rateType
  k32::rateType
  k33::rateType
  k34::rateType
  k35::rateType
  tmp::uType
  atmp::uEltypeNoUnits
  k::rateType
  tab::TabType
end

u_cache(c::Feagin14Cache) = (c.atmp,)
du_cache(c::Feagin14Cache) = (c.fsalfirst,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8,c.k9,c.k10,c.k11,c.k12,c.k13,c.k14,c.k15,c.k16,c.k17,c.k18,c.k19,c.k20,c.k21,c.k22,c.k23,c.k24,c.k25,c.k26,c.k27,c.k28,c.k29,c.k30,c.k31,c.k32,c.k33,c.k34,c.k35,c.k)

function alg_cache(alg::Feagin14,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  tab = Feagin14ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
  k1 = zeros(rate_prototype); k2 = zeros(rate_prototype); k3 = zeros(rate_prototype); k4 = zeros(rate_prototype); k5 = zeros(rate_prototype)
  k6 = zeros(rate_prototype); k7 = zeros(rate_prototype); k8 = zeros(rate_prototype); k9 = zeros(rate_prototype); k10 = zeros(rate_prototype)
  k11 = zeros(rate_prototype); k12 = zeros(rate_prototype); k13 = zeros(rate_prototype); k14 = zeros(rate_prototype)
  k15 = zeros(rate_prototype); k16 = zeros(rate_prototype); k17 = zeros(rate_prototype); k18 = zeros(rate_prototype)
  k19 = zeros(rate_prototype); k20 = zeros(rate_prototype); k21 = zeros(rate_prototype); k22 = zeros(rate_prototype)
  k23 = zeros(rate_prototype); k24 = zeros(rate_prototype); k25 = zeros(rate_prototype)
  k26 = zeros(rate_prototype); k27 = zeros(rate_prototype); k28 = zeros(rate_prototype)
  k29 = zeros(rate_prototype); k30 = zeros(rate_prototype); k31 = zeros(rate_prototype); k32 = zeros(rate_prototype)
  k33 = zeros(rate_prototype); k34 = zeros(rate_prototype); k35 = zeros(rate_prototype)
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits,indices(u)); k = zeros(rate_prototype)

  Feagin14Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,
                k17,k18,k19,k20,k21,k22,k23,k24,k25,k26,k27,k28,k29,k30,
                k31,k32,k33,k34,k35,tmp,atmp,k,tab)
end

alg_cache(alg::Feagin14,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}}) = Feagin14ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))

type Rosenbrock23Cache{uType,uArrayType,rateType,vecuType,JType,TabType,TFType,UFType,F,JCType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k₁::rateType
  k₂::rateType
  k₃::rateType
  du1::rateType
  du2::rateType
  f₁::rateType
  vectmp::vecuType
  vectmp2::vecuType
  vectmp3::vecuType
  fsalfirst::rateType
  fsallast::rateType
  dT::uArrayType
  J::JType
  W::JType
  tmp::uArrayType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::vecuType
  factorization::F
  jac_config::JCType
end

u_cache(c::Rosenbrock23Cache) = (c.dT,c.tmp)
du_cache(c::Rosenbrock23Cache) = (c.k₁,c.k₂,c.k₃,c.du1,c.du2,c.f₁,c.fsalfirst,c.fsallast,c.linsolve_tmp)
jac_cache(c::Rosenbrock23Cache) = (c.J,c.W)
vecu_cache(c::Rosenbrock23Cache) = (c.vectmp,c.vectmp2,c.vectmp3)

type Rosenbrock32Cache{uType,uArrayType,rateType,vecuType,JType,TabType,TFType,UFType,F,JCType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k₁::rateType
  k₂::rateType
  k₃::rateType
  du1::rateType
  du2::rateType
  f₁::rateType
  vectmp::vecuType
  vectmp2::vecuType
  vectmp3::vecuType
  fsalfirst::rateType
  fsallast::rateType
  dT::uArrayType
  J::JType
  W::JType
  tmp::uArrayType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::vecuType
  factorization::F
  jac_config::JCType
end

u_cache(c::Rosenbrock32Cache) = (c.dT,c.tmp)
du_cache(c::Rosenbrock32Cache) = (c.k₁,c.k₂,c.k₃,c.du1,c.du2,c.f₁,c.fsalfirst,c.fsallast,c.linsolve_tmp)
jac_cache(c::Rosenbrock32Cache) = (c.J,c.W)
vecu_cache(c::Rosenbrock32Cache) = (c.vectmp,c.vectmp2,c.vectmp3)

function alg_cache(alg::Rosenbrock23,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  k₁ = zeros(rate_prototype)
  k₂ = zeros(rate_prototype)
  k₃ = zeros(rate_prototype)
  du1 = zeros(rate_prototype)
  du2 = zeros(rate_prototype)
  # f₀ = similar(u) fsalfirst
  f₁ = zeros(rate_prototype)
  vectmp = vec(similar(u,indices(u)))
  vectmp2 = vec(similar(u,indices(u)))
  vectmp3 = vec(similar(u,indices(u)))
  fsalfirst = zeros(rate_prototype)
  fsallast = zeros(rate_prototype)
  dT = similar(u,indices(u))
  J = zeros(uEltypeNoUnits,length(u),length(u)) # uEltype?
  W = similar(J);
  tmp = similar(u,indices(u))
  tab = Rosenbrock23ConstantCache(uEltypeNoUnits,identity,identity)
  vf = VectorF(f,size(u))
  vfr = VectorFReturn(f,size(u))
  tf = TimeGradientWrapper(vf,uprev,du2)
  uf = UJacobianWrapper(vfr,t)
  linsolve_tmp = vec(similar(u,indices(u)))
  jac_config = ForwardDiff.JacobianConfig{determine_chunksize(u,alg)}(vec(du1),vec(uprev))
  Rosenbrock23Cache(u,uprev,k₁,k₂,k₃,du1,du2,f₁,vectmp,vectmp2,vectmp3,fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,alg.factorization,jac_config)
end

function alg_cache(alg::Rosenbrock32,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  k₁ = zeros(rate_prototype)
  k₂ = zeros(rate_prototype)
  k₃ = zeros(rate_prototype)
  du1 = zeros(rate_prototype)
  du2 = zeros(rate_prototype)
  # f₀ = similar(u) fsalfirst
  f₁ = zeros(rate_prototype)
  vectmp = vec(similar(u,indices(u)))
  vectmp2 = vec(similar(u,indices(u)))
  vectmp3 = vec(similar(u,indices(u)))
  fsalfirst = zeros(rate_prototype)
  fsallast = zeros(rate_prototype)
  dT = similar(u,indices(u))
  J = zeros(uEltypeNoUnits,length(u),length(u)) # uEltype?
  W = similar(J); tmp = similar(u,indices(u))
  tab = Rosenbrock32ConstantCache(uEltypeNoUnits,identity,identity)
  vf = VectorF(f,size(u))
  vfr = VectorFReturn(f,size(u))
  tf = TimeGradientWrapper(vf,uprev,du2)
  uf = UJacobianWrapper(vfr,t)
  linsolve_tmp = vec(similar(u,indices(u)))
  jac_config = ForwardDiff.JacobianConfig{determine_chunksize(u,alg)}(vec(du1),vec(uprev))
  Rosenbrock32Cache(u,uprev,k₁,k₂,k₃,du1,du2,f₁,vectmp,vectmp2,vectmp3,fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,alg.factorization,jac_config)
end

immutable Rosenbrock23ConstantCache{T,TF,UF} <: OrdinaryDiffEqConstantCache
  c₃₂::T
  d::T
  tf::TF
  uf::UF
end

function Rosenbrock23ConstantCache(T::Type,tf,uf)
  c₃₂ = T(6 + sqrt(2))
  d = T(1/(2+sqrt(2)))
  Rosenbrock23ConstantCache(c₃₂,d,tf,uf)
end

function alg_cache(alg::Rosenbrock23,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}})
  tf = TimeDerivativeWrapper(f,u)
  uf = UDerivativeWrapper(f,t)
  Rosenbrock23ConstantCache(uEltypeNoUnits,tf,uf)
end

immutable Rosenbrock32ConstantCache{T,TF,UF} <: OrdinaryDiffEqConstantCache
  c₃₂::T
  d::T
  tf::TF
  uf::UF
end

function Rosenbrock32ConstantCache(T::Type,tf,uf)
  c₃₂ = T(6 + sqrt(2))
  d = T(1/(2+sqrt(2)))
  Rosenbrock32ConstantCache(c₃₂,d,tf,uf)
end

function alg_cache(alg::Rosenbrock32,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}})
  tf = TimeDerivativeWrapper(f,u)
  uf = UDerivativeWrapper(f,t)
  Rosenbrock32ConstantCache(uEltypeNoUnits,tf,uf)
end

type ImplicitEulerCache{uType,uArrayType,vecuType,DiffCacheType,rateType,rhsType,adfType,CS} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  uhold::vecuType
  dual_cache::DiffCacheType
  u_old::uArrayType
  tmp::uType
  k::rateType
  fsalfirst::rateType
  rhs::rhsType
  adf::adfType
end

u_cache(c::ImplicitEulerCache)    = (c.uprev2,c.u_old)
du_cache(c::ImplicitEulerCache)   = (c.k,c.fsalfirst)
vecu_cache(c::ImplicitEulerCache) = (c.uhold,)
dual_cache(c::ImplicitEulerCache) = (c.dual_cache,)

function alg_cache(alg::ImplicitEuler,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  u_old = similar(u,indices(u)); k = zeros(rate_prototype)
  dual_cache = DiffCache(u,Val{determine_chunksize(u,alg)})
  uhold = vec(u) # this makes uhold the same values as integrator.u
  rhs = RHS_IE(f,u_old,t,t,dual_cache,size(u),eachindex(u))
  fsalfirst = zeros(rate_prototype)
  if alg_autodiff(alg)
    adf = autodiff_setup(rhs,uhold,alg)
  else
    adf = nothing
  end
  tmp = u_old
  ImplicitEulerCache{typeof(u),typeof(u_old),typeof(uhold),typeof(dual_cache),typeof(k),
                     typeof(rhs),typeof(adf),determine_chunksize(u,alg)}(
                     u,uprev,uprev2,uhold,dual_cache,u_old,tmp,k,fsalfirst,rhs,adf)
end

immutable ImplicitEulerConstantCache{vecuType,rhsType,adfType,CS} <: OrdinaryDiffEqConstantCache
  uhold::vecuType
  u_old::vecuType
  rhs::rhsType
  adf::adfType
end

function alg_cache(alg::ImplicitEuler,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}})
  uhold = Vector{typeof(u)}(1)
  u_old = Vector{typeof(u)}(1)
  rhs = RHS_IE_Scalar(f,u_old,t,t)
  if alg_autodiff(alg)
    adf = autodiff_setup(rhs,uhold,alg)
  else
    adf = nothing
  end
  ImplicitEulerConstantCache{typeof(uhold),typeof(rhs),typeof(adf),1}(uhold,u_old,rhs,adf)
end

type TrapezoidCache{uType,uArrayType,vecuType,DiffCacheType,rateType,rhsType,adfType,CS} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  uhold::vecuType
  u_old::uArrayType
  fsalfirst::rateType
  dual_cache::DiffCacheType
  tmp::uType
  k::rateType
  rhs::rhsType
  adf::adfType
end

u_cache(c::TrapezoidCache)    = (c.uprev2,c.u_old)
du_cache(c::TrapezoidCache)   = (c.k,c.fsalfirst)
vecu_cache(c::TrapezoidCache) = (c.uhold,)
dual_cache(c::TrapezoidCache) = (c.dual_cache,)

function alg_cache(alg::Trapezoid,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
  u_old = similar(u,indices(u)); k = zeros(rate_prototype)
  uhold = vec(u); fsalfirst = zeros(rate_prototype)
  dual_cache = DiffCache(u,Val{determine_chunksize(u,alg)})
  rhs = RHS_Trap(f,u_old,fsalfirst,t,t,size(u),dual_cache,eachindex(u))
  if alg_autodiff(alg)
    adf = autodiff_setup(rhs,uhold,alg)
  else
    adf = nothing
  end
  tmp = u_old
  TrapezoidCache{typeof(u),typeof(u_old),typeof(uhold),typeof(dual_cache),typeof(k),
    typeof(rhs),typeof(adf),determine_chunksize(u,alg)}(
    u,uprev,uprev2,uhold,u_old,fsalfirst,dual_cache,tmp,k,rhs,adf)
end

immutable TrapezoidConstantCache{vecuType,rhsType,adfType,CS} <: OrdinaryDiffEqConstantCache
  uhold::vecuType
  u_old::vecuType
  rhs::rhsType
  adf::adfType
end

function alg_cache(alg::Trapezoid,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}})
  uhold = Vector{typeof(u)}(1)
  u_old = Vector{typeof(u)}(1)
  rhs = RHS_Trap_Scalar(f,u_old,rate_prototype,t,t)
  if alg_autodiff(alg)
    adf = autodiff_setup(rhs,uhold,alg)
  else
    adf = nothing
  end
  TrapezoidConstantCache{typeof(uhold),typeof(rhs),typeof(adf),1}(uhold,u_old,rhs,adf)
end

get_chunksize(cache::DECache) = error("This cache does not have a chunksize.")
get_chunksize{uType,DiffCacheType,rateType,CS}(cache::ImplicitEulerCache{uType,DiffCacheType,rateType,CS}) = CS
get_chunksize{uType,DiffCacheType,rateType,CS}(cache::TrapezoidCache{uType,DiffCacheType,rateType,CS}) = CS
get_chunksize{CS}(cache::ODEChunkCache{CS}) = CS


text:
fname = "savetuple.csv"
csvfile = open(fname,"w")
# writing headers:
write(csvfile, "ColName A, ColName B, ColName C\n")
for i = 1:10
	tup(i) = tuple(rand(Float64,3)...)
	write(csvfile, join(tup(i),","), "\n")
end
close(csvfile)


text:
module PFCNonUniform
export make_advect_1D!

function make_advect_1D!(Δx_::Vector{Float64})
    ξ = 2.0
    Δx = copy(Δx_)

    min_Δx = minimum(Δx)/maximum(Δx)
    ξ = (1.0 + min_Δx)*(1.0 + 2*min_Δx)/(3.0 + (min_Δx - 1.0/min_Δx)^2)

    f_tmp::Vector{Float64} = similar(Δx)

    function ϵ⁺(f::Float64, g::Float64)::Float64
        if f < g
            return min(g - f, ξ*f)
        else
            return max(g - f, -ξ*(1 - f))
        end
    end

    function ϵ⁻(f::Float64, g::Float64)::Float64
        if f < g
            return max(f - g, -ξ*f)
        else
            return min(f - g, ξ*(1 - f))
        end
    end

    Φ⁺(α, f, Δx) = α*(f[2] + (Δx[2] - α)/(Δx[3] + Δx[2] + Δx[1])*(
                    ϵ⁺(f[2], f[3])/(Δx[3] + Δx[2])*(Δx[2] + Δx[1] - α) +
                    ϵ⁻(f[2], f[1])/(Δx[2] + Δx[1])*(Δx[3] + α)))

    Φ⁻(α, f, Δx) = α*(f[2] - (Δx[2] + α)/(Δx[3] + Δx[2] + Δx[1])*(
                    ϵ⁺(f[2], f[3])/(Δx[3] + Δx[2])*(Δx[1] - α) +
                    ϵ⁻(f[2], f[1])/(Δx[2] + Δx[1])*(Δx[2] + Δx[3] + α)))

    function advect_1D!(f, α::Float64)
        f_tmp[:] = copy(f)
        if α > 0
            Φ = Φ⁺(α, [f[end],f[1],f[2]], [Δx[end],Δx[1],Δx[2]])
            f_tmp[1] = f_tmp[1] - Φ/Δx[1]
            f_tmp[2] = f_tmp[2] + Φ/Δx[2]
            for i in 2:length(Δx)-1
                Φ = Φ⁺(α, f[i-1:i+1], Δx[i-1:i+1])
                f_tmp[i] = f_tmp[i] - Φ/Δx[i]
                f_tmp[i+1] = f_tmp[i+1] + Φ/Δx[i+1]
            end
            Φ = Φ⁺(α, [f[end-1],f[end],f[1]], [Δx[end-1],Δx[end],Δx[1]])
            f_tmp[end] = f_tmp[end] - Φ/Δx[end]
            f_tmp[1] = f_tmp[1] + Φ/Δx[1]
        else
            Φ = Φ⁻(α, [f[end],f[1],f[2]], [Δx[end],Δx[1],Δx[2]])
            f_tmp[1] = f_tmp[1] + Φ/Δx[1]
            f_tmp[end] = f_tmp[end] - Φ/Δx[end]
            for i in 2:length(Δx)-1
                Φ = Φ⁻(α, f[i-1:i+1], Δx[i-1:i+1])
                f_tmp[i] = f_tmp[i] + Φ/Δx[i]
                f_tmp[i-1] = f_tmp[i-1] - Φ/Δx[i-1]
            end
            Φ = Φ⁻(α, [f[end-1],f[end],f[1]], [Δx[end-1],Δx[end],Δx[1]])
            f_tmp[end] = f_tmp[end] + Φ/Δx[end]
            f_tmp[end-1] = f_tmp[end-1] - Φ/Δx[end-1]
        end
        f[:] = copy(f_tmp)
    end

    return advect_1D!
end

end


text:
import Pkg

using Pkg.TOML
 const MLJ_VERSION =
     VersionNumber(TOML.parsefile(joinpath(dirname(@__DIR__),
                                           "Project.toml"))["version"])


text:
"""
The Parsing module handles the translation from `.proto` file text into
NamedTuples which contain the parsed content.
The `parse_proto` method maps either text or IO to a nested NamedTuple.
"""
module Parsing

using CombinedParsers
using CombinedParsers.Regexp: whitespace, whitespace_newline

export parse_proto

# a period
period = CharIn('.')

# optional whitespace
ws = Optional(whitespace)

# equals sign with optional surrounding whitespace
_equals = ws * "=" * ws

# handle parsers separated by commas and whitespace
function commasep(x, delim=','*ws)
    join(Repeat(x), delim, infix=:prefix) do (x1, xrest)
        ret = Any[]
        push!(ret, x1)
        for xi in xrest
            push!(ret, xi[2])
        end
        ret
    end
end

# letters and digits
letter = CharIn('A':'Z', 'a':'z')
decimalDigit = CharIn('0':'9')
octalDigit = CharIn('0':'7')
hexDigit = CharIn('0':'9', 'A':'F', 'a':'f')
letterExtended = Either{Any}(letter, decimalDigit, CharIn('_'))

ident = !(letter * Repeat(letterExtended))
fullIdent = !(ident * Repeat(period * ident))
messageName = ident
enumName = ident
fieldName = ident
oneofName = ident
mapName = ident
serviceName = ident
rpcName = ident
messageType = !(Optional(period) * Repeat(ident * period) * messageName)
enumType = messageType

# integer literals
decimalLit = CharIn('1':'9') * Repeat(decimalDigit)
octalLit = CharIn('0') * Repeat(octalDigit)
hexLit = CharIn('0') * CharIn("xX") * Repeat1(hexDigit)
intLit = !Either{Any}(decimalLit, octalLit, hexLit)

# floating point literals
decimals = Repeat1(decimalDigit)
exponent = CharIn("eE") * Optional(CharIn("+-")) * decimals
floatLit = !Either{Any}(
    decimals * period * Optional(decimals) * Optional(exponent),
    decimals * exponent,
    period * decimals * Optional(exponent)
)

# boolean
boolLit = !Either{Any}("true", "false")

# string literals
# hexEscape = '\\' * CharIn("xX") * hexDigit * hexDigit
# octEscape = '\\' * Repeat(3, 3, octalDigit)
# charEscape = CharIn("\a\b\f\n\r\t\v\\\'\"")
# charValue = hexEscape | octEscape | charEscape | CharIn("\0\n\\")
charValue = AnyChar()
_quote = CharIn("'", '"')
strLit = Either{Any}(
    Sequence(2, "'", !Repeat(charValue), "'"),
    Sequence(2, '"', !Repeat(charValue), '"')
)

# empty statement
emptyStatement = ";";

# comments
blockComment = Sequence(2,
    "/*", !Repeat(AnyChar()), "*/"
)
lineComment = Sequence(2,
    "//", !Repeat(CharNotIn("\n\r")), "\n"
)
lineComments = map(lineComment * Repeat(Optional(whitespace_newline) * lineComment)) do (c1, cRest)
    # make one big comment from the multiple line comments
    comms = cat([c1], [ci[2] for ci in cRest], dims=1)
    join(comms, "\n")
end

# this grabs a comment block with any whitespace before or after.
# we'll use this in most places we want to grab comments.
filler = Sequence(2,
    Optional(whitespace_newline),
    Optional(Either{Any}(blockComment, lineComments)),
    Optional(whitespace_newline)
)

# constant
constant = Either{Any}(
	fullIdent,
	Optional(CharIn("+-")) * intLit,
	Optional(CharIn("+-")) * floatLit,
	strLit,
	boolLit
);

# syntax to define proto version
# example:   syntax="proto3";
syntax = !("syntax" * _equals * _quote * "proto3" * _quote * ";");

# import statement
# example: import public "other.proto";
_import = Sequence(
	"import", whitespace,
	:mod=>Optional(Sequence(1, Either("weak", "public"), whitespace)),
	:path=>strLit, ";"
);

# package statement
# example: package foo.bar;
_package = Sequence(
	"package", whitespace, :package=>fullIdent, ";"
);

# option statement
# example:  option java_package = "com.example.foo";
optionName = !Sequence(
    Either{Any}(ident, "(" * fullIdent * ")"),
    Repeat(period * ident)
)
_option = Sequence(
    :comments=>filler,
    "option", whitespace,
    :key=>optionName, _equals, :value=>constant, ";"
)

# fields
fieldType = !Either{Any}(
    "double", "float", "int32", "int64", "uint32", "uint64",
    "sint32", "sint64", "fixed32", "fixed64", "sfixed32", "sfixed64",
    "bool", "string", "bytes", messageType, enumType
)
# get the intLit as an Int64
fieldNumber = map(i -> parse(Int64, i), intLit)
fieldOption = Sequence(:name=>optionName, _equals, :val=>constant)
fieldOptions = fieldOption * Repeat("," * ws * fieldOption)
fieldMaybeOptions = Sequence(2,
    ws,
    Optional(Sequence(2, '[', fieldOptions, ']')),
    ';'
)
field = Sequence(
    Optional("repeated" * whitespace),
    :type=>fieldType, whitespace,
    :name=>fieldName, _equals, :num=>fieldNumber,
    fieldMaybeOptions
)

# oneofs
oneofField = Sequence(
    :type=>fieldType, whitespace,
    :name=>fieldName, _equals, :num=>fieldNumber,
    fieldMaybeOptions
)
oneofFieldWithComments = Sequence(
    :comments=>filler,
    :field=>Either{Any}(_option, oneofField, emptyStatement)
)
oneof = Sequence(
    :comments=>filler,
    "oneof", whitespace, :oneofName=>oneofName, ws, '{',
    :fields=>Repeat(oneofFieldWithComments),
    filler, '}'
)

# maps
keyType = Either(
    "int32", "int64', uint32", "unit64", "sint32", "sint64",
    "fixed32", "fixed64", "sfixed32", "sfixed64", "bool", "string"
)
mapField = Sequence(
    "map", ws, "<",
    :keyType=>keyType,
    ",", ws,
    :valueType=>fieldType,
    ">", ws,
    :name=>mapName, _equals, :num=>fieldNumber,
    fieldMaybeOptions
)

# reserved statements
# example: reserved 2, 15, 9 to 11;
# example: reserved "foo", "bar";
fieldNames = commasep(fieldName)
range = Either{Any}(
    !Sequence(intLit, " to ", intLit | "max"),
    intLit
)
ranges = commasep(range)
reserved = Sequence(
    "reserved", whitespace, :reserved=>Either{Any}(ranges, fieldNames)
)

# enums
# example:
#	enum XX {
#		option allow_alias = true;
#		UNKNOWN = 0;
#	}
enumValueOption = optionName * _equals * constant
enumField = Sequence(
    :fieldname=>ident, _equals,
    :fieldnum=>map(i -> parse(Int64, i), !Sequence(Optional('-'), intLit)),
    Optional(Sequence('[', enumValueOption,
            Repeat(Sequence(',', ws, enumValueOption)), ']')), ';'
)
enumFieldWithComments = Sequence(
    :comments=>filler,
    :field=>Either{Any}(_option, enumField, emptyStatement)
)
enum = Sequence(
    :comments=>filler,
    "enum", whitespace, :enumName=>enumName, ws, '{',
    :fields=>Repeat(enumFieldWithComments),
    filler,
    '}'
)

# message definitions
innerFieldWithComments = Sequence(
    :comments=>filler,
    :field=>Either{Any}(field, _option, oneof, mapField, reserved, emptyStatement)
)
innerMessage = Sequence(
    :comments=>filler,
    "message", whitespace, :messageName=>messageName, ws,
    '{',
    :fields=>Repeat(innerFieldWithComments),
    filler,
    '}'
)
messageFieldWithComments = Either{Any}(
    innerMessage,
    enum,
    innerFieldWithComments
)
# Sequence(
#     :comments=>filler,
#     :field=>Either{Any}(field, enum, _option, oneof, mapField, reserved, emptyStatement, innerMessage)
# )
message = Sequence(
    :comments=>filler,
    "message", whitespace, :messageName=>messageName, ws,
    '{',
    :fields=>Repeat(messageFieldWithComments),
    filler,
    '}'
)

# service definitions

# returns true when matching "stream "
isstream = map(x -> x=="stream ", Optional("stream "))

rpc = Sequence(
    :comments=>filler,
    "rpc", whitespace,
    :rpcName=>rpcName, ws,
    '(', :inputIsStream=>isstream, :inputType=>messageType, ')',
    whitespace, "returns", whitespace,
    '(', :outputIsStream=>isstream, :outputType=>messageType, ')',
    ws, Either{Any}(
        Sequence('{', :extra=>Optional(_option | emptyStatement), "};"),
        ';'
    )
)

service = Sequence(
    :comments=>filler,
    "service", whitespace, :serviceName=>serviceName, ws, '{',
    :rpcs=>Repeat(Sequence(2, whitespace_newline, Either{Any}(_option, rpc, emptyStatement))),
    whitespace_newline, '}'
)

# a complete proto file
topLevelDef = Either{Any}(message, enum, service)
proto = Sequence(
    :comment=>filler,
    syntax,
    :statements=>Repeat(Sequence(2,
        whitespace_newline,
        Either{Any}(_import, _package, _option, topLevelDef, emptyStatement)
    )),
    filler,
    :unparsed=>!Repeat(AnyChar()) # anything we missed
)

function parse_proto(text::AbstractString)
    proto(text)
end

function parse_proto(io::IO)
    text = read(io, String)
    parse_proto(text)
end

end # module


text:
using StatsModels
using RegressionFormulae
using Test

include("dummymod.jl")

dat = (; y=zeros(3), a=1:3, b=11:13, c=21:23, d=31:33, e=["u", "i", "o"])

@testset "error checking" begin
    @test_throws ArgumentError (term(:b),) ^ term(:a)
    @test_throws ArgumentError (term(:b),) ^ term(2.5)
    @test_throws ArgumentError (term(:b),) ^ term(-2)

    sch = schema(dat)
    @test_throws ArgumentError apply_schema(@formula(y ~ ^(a, b, c)), sch, RegressionModel)
    @test_throws ArgumentError apply_schema(@formula(y ~ (a+b)^2.5), sch, RegressionModel)
end

@testset "powers of sums" begin
    m = fit(DummyMod, @formula(y ~ (a + b + c + d)^3), dat)
    @test coefnames(m) == ["(Intercept)", "a", "b", "c", "d",
                           "a & b", "a & c", "a & d", "b & c", "b & d", "c & d",
                           "a & b & c", "a & b & d", "a & c & d", "b & c & d"]

    m = fit(DummyMod, @formula(y ~ (1 + a + b + c + d)^3), dat)
    @test_broken coefnames(m) == ["(Intercept)", "a", "b", "c", "d",
                                  "a & b", "a & c", "a & d", "b & c", "b & d", "c & d",
                                  "a & b & c", "a & b & d", "a & c & d", "b & c & d"]

    m = fit(DummyMod, @formula(y ~ (a + b + e)^2), dat)
    @test coefnames(m) == ["(Intercept)", "a", "b", "e: o", "e: u",
                           "a & b", "a & e: o", "a & e: u", "b & e: o", "b & e: u"]
end

@testset "embedded interactions" begin
    m = fit(DummyMod, @formula(y ~ (a + b + c * d)^3), dat)
    cn = coefnames(m)
    @test_broken !("a & c & c & d" in cn)
    @test_broken cn == unique(cn)
end


text:
"""
This package provides an MIT license, dependency-free implementation of
Algorithm 1 of "Sinkhorn Divergences for Unbalanced Optimal Transport" [SFVTP19
(http://arxiv.org/abs/1910.12958)], which allows calculation of the optimal
transport and Sinkhorn divergence between two positive measures (with possibly
different total mass), where mass creation and destruction is penalized by one
of several possible φ-divergences.
"""
module UnbalancedOptimalTransport

using LinearAlgebra: norm

export DiscreteMeasure,
    OT!, optimal_coupling!, sinkhorn_divergence!, unbalanced_sinkhorn!, cost_matrix

struct DiscreteMeasure{P,LP,S,T}
    density::P
    log_density::LP
    set::S
    dual_potential::Vector{T}
    cache::Vector{T}
end

"""
    DiscreteMeasure(density, [log_density], set) -> DiscreteMeasure

Construct a `DiscreteMeasure` object for use in [`unbalanced_sinkhorn!`](@ref)
and related functions.

* `density` should be strictly positive; zero elements should instead be removed
  from `set`
* `log_density` should be equal to `log.(density)` and can be omitted
  (in which case its computed automatically)
* `set` is a collection so that `density[i]` is the probability of the element
  `set[i]` occurring (where `i ∈ eachindex(density, set)`).

"""
function DiscreteMeasure(density::P, log_density::LP, set::S) where {P,LP,S}
    T = eltype(density)
    n = length(density)
    n == length(log_density) ||
    throw(ArgumentError("`density`, `log_density`, and `set` (if supplied) must have equal length"))
    set !== nothing &&
    length(set) != n &&
    throw(ArgumentError("`density`, `log_density`, and `set` must have equal length"))
    dual_potential = zeros(T, n)
    cache = similar(dual_potential)
    DiscreteMeasure{P,LP,S,T}(density, log_density, set, dual_potential, cache)
end

DiscreteMeasure(density, set = nothing) = DiscreteMeasure(density, log.(density), set)

Base.eltype(::DiscreteMeasure{P,LP,S,T}) where {P,LP,S,T} = T
Base.length(a::DiscreteMeasure) = length(a.density)

include("divergences.jl")
include("sinkhorn.jl")
include("optimized_methods.jl")
include("utilities.jl")

end


text:
export narrow_class_group

doc"""
***
    power_reduce2(A::NfOrdIdl, e::fmpz) -> NfOrdIdl, FacElem{nf_elem}
> Computes $B$ and $\alpha$ in factored form, such that $\alpha B = A^e$
"""
function power_reduce2(A::NfOrdIdl, e::fmpz)
  A_orig = deepcopy(A)
  e_orig = e

  O = order(A)
  K= nf(O)
  if norm(A) > abs(discriminant(O))
    A, a = reduce_ideal2(A)
    @hassert :PID_Test 1 a*A_orig == A
    # A_old * inv(a) = A_new
    #so a^e A_old^e = A_new^e
    al = FacElem(Dict(a=>-e))
  else
    al = FacElem(Dict(K(1) => fmpz(1)))
  end

  #we have A_orig^e = (A*a)^e = A^e*a^e = A^e*al and A is now small

  if e < 0
    B = inv(A)
    A = numerator(B)
    al *= FacElem(Dict(K(denominator(B)) => fmpz(e)))
    e = -e
  end
  # A^e = A^(e/2)^2 A or A^(e/2)^2
  # al * A^old^(e/2) = A_new
  if e>1
    C, cl = power_reduce2(A, div(e, 2))
    @hassert :PID_Test 1 C*evaluate(cl) == A^Int(div(e, 2))

    C2 = C^2
    al = al*cl^2
    if norm(C2) > abs(discriminant(O))
      C2, a = reduce_ideal2(C2)
      al *= inv(a)
    end

    if isodd(e)
      A = C2*A
      if norm(A) > abs(discriminant(O))
        A, a = reduce_ideal2(A)
        al *= inv(a)
      end
    else
      A = C2
    end
    return A, al
  else
    @assert e==1
    return A, al
  end
end

function ==(A::NfOrdIdl, B::NfOrdFracIdl)
  C = simplify(B)
  if C.den != 1 ||
     C.num != A
    return false
  else
    return true
  end
end
==(A::NfOrdFracIdl, B::NfOrdIdl) = B==A

doc"""
***
    reduce_ideal2(A::FacElem{NfOrdIdl}) -> NfOrdIdl, FacElem{nf_elem}
> Computes $B$ and $\alpha$ in factored form, such that $\alpha B = A$.
"""
function reduce_ideal2(I::FacElem{NfOrdIdl, NfOrdIdlSet})
  O = order(first(keys(I.fac)))
  K = nf(O)
  fst = true
  a = FacElem(Dict(K(1) => fmpz(1)))
  A = ideal(O, 1)

  for (k,v) = I.fac
    @assert order(k) == O
    if iszero(v)
      continue
    end
    if fst
      A, a = power_reduce2(k, v)
      @hassert :PID_Test 1 (v>0? k^Int(v) : inv(k)^Int(-v)) == A*evaluate(a)
      fst = false
    else
      B, b = power_reduce2(k, v)
      @hassert :PID_Test (v>0? k^Int(v) : inv(k)^Int(-v)) == B*evaluate(b)
      A = A*B
      a = a*b
      if norm(A) > abs(discriminant(O))
        A, c = reduce_ideal2(A)
        a = a*FacElem(Dict(K(c) => fmpz(-1)))
      end
    end
  end
  @hassert :PID_Test 1 A*evaluate(a) == evaluate(I)
  return A, a
end

mutable struct MapNarrowClassGrp{T} <: Map{T, NfOrdIdlSet, HeckeMap, MapNarrowClassGrp}
  header::MapHeader

  function MapNarrowClassGrp{T}() where {T}
    return new{T}()
  end
end

function show(io::IO, mC::MapNarrowClassGrp)
  println(io, "NarrowClassGroup map of $(codomain(mC))")
end

#TODO: make it for AnticNumberField and for NfOrd. I that
#      case also return OrdElem

doc"""
***
    elements_with_all_signs(L::NfOrd) -> Array{nf_elem, 1}
> Finds elements $x_i$ in the number field s.th the elements
> totally positive at all but the $i$-th real embedding, but
> negative at the $i$-th.
"""
function elements_with_all_signs(L::NfOrd)
  r1, r2 = signature(L)
  K = nf(L)

  S = DiagonalGroup([2 for i=1:r1])

  inf_places = infinite_places(K)

  function logS(x)
    return S([x[P] > 0 ? 0 : 1 for P in inf_places])
  end

  s = typeof(S[1])[]
  g = nf_elem[]
  u, mu = sub(S, s)
  b = 10
  cnt = 0
  while true
    a = rand(L, b)
    t = logS(signs(K(a)))
    if !haspreimage(mu, t)[1]
      push!(s, t)
      push!(g, K(a))
      u, mu = sub(S, s)
      if order(u) == order(S)
        break
      end
    else
      cnt += 1
      if cnt > 100
        b *= 2
        cnt = 0
      end
    end
  end
  X = DiagonalGroup([2 for i=1:r1])
  hX = GrpAbFinGenMap(X, X, vcat([x.coeff for x=s]))
  r = nf_elem[]
  for i=1:r1
    y = haspreimage(hX, X[i])[2]
    push!(r, prod([g[i]^Int(y[i]) for i=1:r1]))
  end
  return r
end

doc"""
***
    narrow_class_group(L::NfOrd) -> GrpAbFinGen, Map
> Compute the narrow (or strict) class group of $L$, ie. the
> group of invertable ideals modulo the group of totally positive elements.
> In case the field has no real embedding, this is just the class group.
"""
function narrow_class_group(L::NfOrd)
  C, mC = class_group(L)

  r1, r2 = signature(L)
  if r1 == 0
    return C, mC
  end

  U, mU = unit_group_fac_elem(L)

  K = nf(L)

  inf_places = infinite_places(K)

  # 1 -> K^*/K^+U -> Cl^+ -> Cl -> 1
  # K^*/K^+ = C_2^{r_1}

  gensC = [mC(C[i]) for i=1:ngens(C)]

  S = DiagonalGroup([2 for i=1:r1])

  function logS(x)
    return S([x[P] > 0 ? 0 : 1 for P in inf_places])
  end
  s, ms = quo(S, [logS(signs(mU(U[i]))) for i=1:ngens(U)])
  s, mms = snf(s)
  ms = mms*inv(ms)

  # we want the extension 1-> s -> X -> C -> 1
  # generators are gensC supplemented by gens(s)

  RC = rels(C)
  R = hcat(zero_matrix(FlintZZ, rows(RC), ngens(s)), RC)
  for i=1:ngens(C)
    A, al = power_reduce2(gensC[i], order(C[i]))
    be = principal_gen_fac_elem(A)
    x = ms(logS(signs(al*be)))
    for j=1:ngens(s)
      R[i,j] = x.coeff[1, j]
    end
  end
  R = vcat(R, hcat(rels(s), zero_matrix(FlintZZ, nrels(s), ngens(C))))

  X = AbelianGroup(R)

  el = elements_with_all_signs(L)
  sign_gen = [preimage(ms, s[i]) for i=1:ngens(s)]
  sg = [ prod([el[i]^Int(x[i]) for i=1:length(el)]) for x=sign_gen]

  function exp(a::GrpAbFinGenElem)
    return L(prod([sg[i]^Int(a.coeff[1,i]) for i=1:length(sg)])) * image(mC, C([a.coeff[1, i+length(sg)] for i=1:ngens(C)]))
  end

  function log(A::NfOrdIdl)
    a = preimage(mC, A)
    B = FacElem(Dict(A => fmpz(-1))) * FacElem(gensC, [a.coeff[1,i] for i=1:ngens(C)])
    A, c = reduce_ideal2(B)
    b = c*principal_gen_fac_elem(A)
    return X(hcat(ms(logS(signs(b))).coeff, a.coeff))
  end

  mp = MapNarrowClassGrp{typeof(X)}()
  mp.header = MapHeader(X, NfOrdIdlSet(L), exp, log)
  return X, mp
end


text:
# v0.6

function combsort!(x::Array)::Array
    gap, swaps = length(x), true
    while gap > 1 || swaps
        gap = floor(Int, gap / 1.25)
        i, swaps = 0, false
        while i + gap < length(x)
            if x[i+1] > x[i+1+gap]
                x[i+1], x[i+1+gap] = x[i+1+gap], x[i+1]
                swaps = true
            end
            i += 1
        end
    end
    return x
end

x = randn(100)
@show x combsort!(x)
@assert issorted(x)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode