text:
"""`Jn, J0, T, Ti, Tn = jacobian(qin, DH::DH, tool=eye(4))`
Calculates the jacobian in base (J0) and tool frame (Jn) as well as the forward kinematics (T),
given the joint angles and the DH-parameters"""
function jacobian(q, DH::DH, tool=eye(4))
    n_joints    = size(q,1)
    dhpar       = DH.dhpar
    # q           = qin + DH.offset # Moved to dh2Tn
    trans       = dh2Tn(DH,q,tool)

    # ## jacobn (tool)
    Jn = zeros(6,n_joints)
    U = tool #Alternatively transformation to TCP <-------------- obs
    for j=n_joints:-1:1
        U = trans[:,:,j] * U
        # revolute axis
        d = [-U[1,1]*U[2,4]+U[2,1]*U[1,4];
        -U[1,2]*U[2,4]+U[2,2]*U[1,4];
        -U[1,3]*U[2,4]+U[2,3]*U[1,4]]
        delta = U[3,1:3]	# nz oz az

        Jn[:,j] = [d; delta]
    end

    ## fkine
    Ti          = zeros(4,4,n_joints+1)
    Ti[:,:,1]   = trans[:,:,1]
    T           = trans[:,:,1]
    for i=2:n_joints
        T           = T * trans[:,:,i]
        Ti[:,:,i] = Ti[:,:, i-1] * trans[:,:,i]
    end
    Ti[:,:,end] = Ti[:,:,end-1]*tool
    T           = T*tool

    # jacob0 (base)
    R = T[1:3,1:3]
    J0 = [R zeros(3,3); zeros(3,3) R] * Jn

    return Jn, J0, T, Ti, trans
end

"""
`jacobianPOE(q, xi)` Returns The jacobian in 1) the base frame, 2) the tool frame. It can most likely be rewritten to be faster.
"""
function jacobianPOE(q, xi)
    #Page 117 in Murray 94
    n_joints  = size(q,1)
    Jss       = zeros(6,n_joints)
    adint     = eye(6)
    T         = eye(4)
    for i = 1:n_joints
        Jss[:,i] = adint*xi[:,i]
        Ti      = expξ(xi[:,i],q[i])
        T       = T*Ti
        adint   = adint*ad(Ti)
    end
    T = T * expξ(xi[:,end-1],1)*expξ(xi[:,end],1) # T is now the full FK

    # Now we do a bit messing around to get the jacobian we are used to from extCtrl
    Jbb = adi(T)*Jss
    Jsb = [T[1:3,1:3] zeros(3,3); zeros(3,3) T[1:3,1:3]] * Jbb

    return Jsb, Jbb, T
end

function jacobianPOEb(q, xi)
    #Page 117 in Murray 94
    n_joints  = size(q,1)
    Jbb       = zeros(6,n_joints)
    T         = expξ2(xi[:,end-1],1)*expξ2(xi[:,end],1)
    adint     = adi(T)
    Ti        = eye(4)
    for i = n_joints:-1:1
        expξ2!(Ti,xi[:,i],q[i])
        adint   = adint*adi(Ti)
        Jbb[:,i]= adint*xi[:,i]
        T       = Ti*T
    end

    # Now we do a bit messing around to get the jacobian we are used to from extCtrl
    #Jbb = adi(T)*Jss
    Jsb = [T[1:3,1:3] zeros(3,3); zeros(3,3) T[1:3,1:3]] * Jbb

    return Jsb, Jbb, T
end

"""
`jacobianPOE(q, xi, T)` Returns The jacobian in 1) the base frame, 2) the tool frame and returns the FK. It must be handed a precomputed end transform `Tend = expξ2(xi[:,end-1],1)*expξ2(xi[:,end],1)`
"""
@fastmath function jacobianPOEikine(q, xi,T)
    #Page 117 in Murray 94
    n_joints  = size(q,1)
    Jbb       = zeros(6,n_joints)
    adint     = adi(T)
    Ti        = eye(4)
    @inbounds for i = n_joints:-1:1
        expξ2!(Ti,xi[:,i],q[i])
        adint   = adint*adi(Ti)
        Jbb[:,i]= adint*xi[:,i]
        T       = Ti*T
    end

    # Now we do a bit messing around to get the jacobian we are used to from extCtrl
    Jsb = [T[1:3,1:3] zeros(3,3); zeros(3,3) T[1:3,1:3]] * Jbb

    return Jsb, Jbb, T
end

"""Computes a set of local transformation matrices given the DH-parameters of a robot. Can be sent an optional joint-angle vector and a tool"""
function dh2Tn{P}(DH, q::VecOrMat{P}, tool=eye(4))
    dhpar    = DH.dhpar
    n_joints = size(dhpar,1)
    Tn       = zeros(P,4,4,n_joints+1)
    q = q .+ DH.offset
    for j=1:n_joints
        Tn[:,:,j] = [    cos(q[j])   -sin(q[j])*cos(dhpar[j,1])    sin(q[j])*sin(dhpar[j,1])    dhpar[j,2]*cos(q[j]);
        sin(q[j])   cos(q[j])*cos(dhpar[j,1])     -cos(q[j])*sin(dhpar[j,1])   dhpar[j,2]*sin(q[j]);
        0                 sin(dhpar[j,1])                     cos(dhpar[j,1])            dhpar[j,4];
        0                 0                                   0                                  1 ]
    end
    Tn[:,:,end] = tool
    return Tn
end

"""Computes a set of local transformation matrices given the DH-parameters of a robot. Can be sent an optional joint-angle vector and a tool"""
function dh2Tn(DH, q=zeros(size(DH.dhpar,1)), tool=eye(4))
    dhpar    = DH.dhpar
    n_joints = size(dhpar,1)
    Tn       = zeros(4,4,n_joints+1)
    q = q .+ DH.offset
    for j=1:n_joints
        Tn[:,:,j] = [    cos(q[j])   -sin(q[j])*cos(dhpar[j,1])    sin(q[j])*sin(dhpar[j,1])    dhpar[j,2]*cos(q[j]);
        sin(q[j])   cos(q[j])*cos(dhpar[j,1])     -cos(q[j])*sin(dhpar[j,1])   dhpar[j,2]*sin(q[j]);
        0                 sin(dhpar[j,1])                     cos(dhpar[j,1])            dhpar[j,4];
        0                 0                                   0                                  1 ]
    end
    Tn[:,:,end] = tool
    return Tn
end

"""`fkinePOE(xi0,q)` Forward kinematics using POE"""
function fkinePOE(xi0,q)
    Tfull = expξ(xi0[:,end-1],1)*expξ(xi0[:,end],1) # We first handle the tool transform, which does not include a joint variable
    for j = size(xi0,2)-2:-1:1
        ξ   = xi0[:,j]
        Tfull = expξ(ξ,q[j])*Tfull
    end
    return Tfull
end

"""`fkineLPOE(Tn0,xi,q)` Forward kinematics using LPOE"""
function fkineLPOE(Tn0,xi,q)
    T = eye(4)
    n = size(xi,2)-1
    for j = 1:n
        T = T*Tn0[:,:,j]*expξ(xi[:,j],q[j]) # turn all joints
    end
    T = T*Tn0[:,:,end]
    return T
end

"""
`function ikinePOE(xi,T,q0; maxiter=100, λ = 1e0, tol = 1e-12, verbose = false, adaptive = true)`
Iterative inverse kinematics
`λ`sets the start value of the regularizing parameter (kind of like inverse step size). Bigger value means slower but more robust convergence, if adaptive is set to true, λ is adapted to speed up convergence.
"""
function ikinePOE(xi,T,q0; maxiter=100, λ = 1e0, tol = 1e-12, verbose = false, adaptive = true)
    q = q0
    err = Array(Float64,maxiter)
    Tend = expξ2(xi[:,end-1],1)*expξ2(xi[:,end],1)
    J,V,Tc = jacobianPOEikine(q,xi,Tend)
    err[1] = norm(twistcoords(logT(T*trinv(Tc))))
    newjac = true
    iter = 2
    while iter <= maxiter

        dq = (J'J + λ*I)\J'*twistcoords(logT(T*trinv(Tc)))
        qc = q + dq
        J,V,Tc = jacobianPOEikine(qc,xi,Tend)
        err[iter] = norm(twistcoords(logT(T*trinv(Tc))))
        if err[iter] > err[iter-1]
            if adaptive
                λ *= 8
            end
            err[iter] = err[iter-1]
            newjac = false
        else
            q = qc
            if adaptive
                λ /= 2
            end
            newjac = true
        end
        verbose && println("Error: ",round(err[iter],7), " λ: ", round(λ,10), " Norm dq: ", round(norm(dq),7))
        if err[iter] < tol || norm(dq) < tol
            err = err[1:iter]
            break
        end

        iter += 1
    end
    return q, err
end

"""
Run this function with a string representing the robot you want the kinematic functions for, e.g. `get_kinematic_functions("YuMi")`
Currently supports YuMi and ABB IRB7600
returns `fkine(q), ikine(T,q0), jacobian(q)`
TODO: implement YuMi_left, Yumi_right Tbase*ikinePOE
"""
function get_kinematic_functions(robot)
    robot = lowercase(robot)
    if robot == "yumi" || robot == "frida"
        dh = DHYuMi()
    elseif robot == "7600" || robot == "irb7600"
        dh = DH7600()
    end
    Tn0 = dh2Tn(dh)
    xi = DH2twistsPOE(dh)
    xiL = DH2twistsLPOE(dh)
    Z = zeros(3,3)
    if robot == "yumileft"
        baseAnglesLeft  = [-0.63 , 0.95 , -0.18]
        Rbase           = rpy2R(baseAnglesLeft,"xyz")
        Tbase           = eye(4)
        Tbase[1:3,1:3]  = Rbase
        fkine = (q) -> Tbase*fkinePOE(xi,q)
        jacobian = (q) -> [Rbase Z;Z Rbase]*jacobianPOE(q,xi)
        ikine = (T,q0, maxiter=100, λ = 1e0, tol = 1e-12, verbose = false) -> ikinePOE(xi,trinv(Tbase)*T,q0,maxiter=maxiter, λ = λ, tol = tol, verbose = verbose)
    elseif robot == "yumiright"
        baseAnglesLeft  = [0.63 , 0.95 , 0.18]
        Rbase           = rpy2R(baseAnglesLeft,"xyz")
        Tbase           = eye(4)
        Tbase[1:3,1:3]  = Rbase
        fkine = (q) -> Tbase*fkinePOE(xi,q)
        jacobian = (q) -> [Rbase Z;Z Rbase]*jacobianPOE(q,xi)
        ikine = (T,q0, maxiter=100, λ = 1e0, tol = 1e-12, verbose = false) -> ikinePOE(xi,trinv(Tbase)*T,q0,maxiter=maxiter, λ = λ, tol = tol, verbose = verbose)
    else
        fkine = (q) -> fkinePOE(xi,q)
        jacobian = (q) -> jacobianPOE(q,xi)
        ikine = (T,q0, maxiter=100, λ = 1e0, tol = 1e-12, verbose = false) -> ikinePOE(xi,T,q0,maxiter=maxiter, λ = λ, tol = tol, verbose = verbose)
    end

    return fkine, ikine, jacobian

end

function testikine()
    dh = DH7600();
    q0 = randn(6)
    xi = DH2twistsPOE(dh)
    T0 = fkinePOE(xi,q0)
    Tt = deepcopy(T0)
    Tt[1:3,4] += 0.1*[1,1,1]
    Tt[1:3,1:3] *= expω(1*π/180*randn(3))
    Profile.clear()
    gc()
    @time q,err = ikinePOE(xi,Tt,q0,verbose=true)

    # pp = semilogy(err)
    # display(pp)
end

function testJacobian()
    include("../POEutils.jl")
    dh = DH7600();
    q = [1,1,0,0,0,0];
    q = randn(6)
    AAA,BBB,T0,Ti0,Tn0 = jacobian(zeros(6),dh, eye(4));
    Jn,J0,Tjac,Tijac,Tnjac = jacobian(q,dh, eye(4));

    xiL = DH2twistsLPOE(Tn0)
    xi = DH2twistsPOE(Tn0)
    Jsb, Jbb = jacobianPOE(q,xi)

    display(round(J0,3))
    display(round(Jsb,3))
end


text:
abstract type RosenbrockTableau{T,T2} end
struct RosenbrockFixedTableau{T,T2}<:RosenbrockTableau{T,T2}
    a::Array{T,2}
    C::Array{T,2}
    b::Array{T,1}
    gamma::T2
    d::Array{T,1}
    c::Array{T2,1}
end

struct RosenbrockAdaptiveTableau{T,T2}<:RosenbrockTableau{T,T2}
    a::Array{T,2}
    C::Array{T,2}
    b::Array{T,1}
    btilde::Array{T,1}
    gamma::T2
    d::Array{T,1}
    c::Array{T2,1}
end

"""
    @_bitarray2boolarray RosenbrockTableau(tab.a.!=0,...)

Transform BitArray (in the form of `xs.!=0` ) into 1D-Array of Bools by
`[i for i in xs.!=0]` to satisfy the type constraint of RosenbrockTableau
"""
macro _bitarray2boolarray(expr)
    args=[:([i for i in $arg]) for arg in expr.args[2:end]]
    args[end-2]=:(tab.gamma!=0)
    esc(:($(expr.args[1])($(args...))))
end

"""
    _masktab(tab)

Convert normal tableau into a dummy tableau consisting of Bools. We use dummy tableaus
where we only care about whether values in the tableau are zeros.
"""
_masktab(tab::RosenbrockFixedTableau)=@_bitarray2boolarray RosenbrockFixedTableau(tab.a.!=0,tab.C.!=0,tab.b.!=0,tab.gamma!=0,tab.d.!=0,tab.c.!=0)
_masktab(tab::RosenbrockAdaptiveTableau)=@_bitarray2boolarray RosenbrockAdaptiveTableau(tab.a.!=0,tab.C.!=0,tab.b.!=0,tab.btilde.!=0,tab.gamma!=0,tab.d.!=0,tab.c.!=0)

"""
    _common_nonzero_vals(tab::RosenbrockTableau)

Return the common nonzero symbols in the tableau. Typical return value:
`[[:a21,:a31,:a32],[:C21,:C31,:C32],[:b1,:b2,:b3],:gamma,[:d1,:d2,:d3],[:c1,:c2,:c3]]`
"""
function _common_nonzero_vals(tab::RosenbrockTableau)
    nzvals=[]
    push!(nzvals,[Symbol(:a,ind[1],ind[2]) for ind in findall(!iszero,tab.a)])
    push!(nzvals,[Symbol(:C,ind[1],ind[2]) for ind in findall(!iszero,tab.C)])
    push!(nzvals,[Symbol(:b,ind) for ind in findall(!iszero,tab.b)])
    push!(nzvals,:gamma)
    push!(nzvals,[Symbol(:d,ind) for ind in findall(!iszero,tab.d)])
    push!(nzvals,[Symbol(:c,ind) for ind in findall(!iszero,tab.c)])
    nzvals
end

"""
    _nonzero_vals(tab::RosenbrockFixedTableau)

Return all the nonzero symbols in the tableau. Typical return value:
`[:a21,:a31,:a32,:C21,:C31,:C32,:b1,:b2,:b3,:gamma,:d1,:d2,:d3,:c1,:c2,:c3]`
"""
function _nonzero_vals(tab::RosenbrockFixedTableau)
    nzvals=_common_nonzero_vals(tab)
    vcat(nzvals...)
end

"""
    _nonzero_vals(tab::RosenbrockAdaptiveTableau)

Typical return value:
`[:a21,:a31,:a32,:C21,:C31,:C32,:b1,:b2,:b3,:btilde1,:btilde2,:btilde3,:gamma,:d1,:d2,:d3,:c1,:c2,:c3]`
"""
function _nonzero_vals(tab::RosenbrockAdaptiveTableau)
    nzvals=_common_nonzero_vals(tab)
    push!(nzvals,[Symbol(:btilde,ind) for ind in findall(!iszero,tab.btilde)])
    vcat(nzvals...)
end

"""
    _push_assigns!(valexprs,inds,name,type)

Insert a series of field statements like `[:(c2::T2),:(c3::T2)]` into the array `valexprs`.

# Arguments
- `valexpr::Array{Expr,1}`: the array to be inserted
- `inds`: an iterator that gives indices
- `name::Symbol`: the prefix name of the values
- `type::Symbol`: type in the statements
"""
function _push_assigns!(valexprs,inds,name,type::Symbol)
    for ind in inds
        push!(valexprs,:($(Symbol(name,"$(Tuple(ind)...)"))::$type))
    end
end

"""
    gen_tableau_struct(tab::RosenbrockTableau,tabstructname::Symbol)

Generate the tableau struct expression from a given tableau emulating those in
`tableaus/rosenbrock_tableaus.jl`. The statements of `aij`,`Cij` and `ci` are generated
according to the nonzero values of `tab.a`,`tab.C` and `tab.c` while others are generated
from their indices. One may choose to pass in a dummy tableau with type `<:RosenbrockTalbeau{Bool,Bool}`
to fully control the tableau struct.
"""
function gen_tableau_struct(tab::RosenbrockTableau,tabstructname::Symbol)
    valexprs=Array{Expr,1}()
    _push_assigns!(valexprs,findall(!iszero,tab.a),:a,:T)
    _push_assigns!(valexprs,findall(!iszero,tab.C),:C,:T)
    _push_assigns!(valexprs,eachindex(tab.b),:b,:T)
    if typeof(tab)<:RosenbrockAdaptiveTableau
        _push_assigns!(valexprs,eachindex(tab.btilde),:btilde,:T)
    end
    push!(valexprs,:(gamma::T2))
    _push_assigns!(valexprs,eachindex(tab.d),:d,:T)
    _push_assigns!(valexprs,findall(!iszero,tab.c),:c,:T2)
    quote struct $tabstructname{T,T2}
        $(valexprs...)
        end
    end
end

"""
    gen_tableau(tab::RosenbrockTableau,tabstructexpr::Expr,tabname::Symbol)

Generate the tableau function expression emulating those in `tableaus/rosenbrock_tableaus.jl`.
It takes in the tableau struct expression (generated by gen_tableau_sturct(...) or written by hand)
to make sure the actual values of the tableau are organized in the right order.
"""
function gen_tableau(tab::RosenbrockTableau,tabstructexpr::Expr,tabname::Symbol)
    @capture(tabstructexpr, struct T_ fields__ end) || error("incorrect tableau expression")
    tabstructname = namify(T)
    valsym2tabdict=Dict("a"=>tab.a,"C"=>tab.C,"gamma"=>tab.gamma,"c"=>tab.c,"d"=>tab.d,"b"=>tab.b)
    if typeof(tab)<:RosenbrockAdaptiveTableau
        valsym2tabdict["btilde"]=tab.btilde
    end
    pattern=r"^([a-zA-Z]+)([1-9]{0,2})$"
    assignexprs = Expr[]
    valsyms = Symbol[]
    for field in fields
        if @capture(field, valsym_Symbol::valtype_)
            push!(valsyms, valsym)
            m = match(pattern, String(valsym))
            val = valsym2tabdict[m[1]][(parse(Int, i) for i in m[2])...]
            push!(assignexprs, :($valsym = convert($valtype, $val)))
        end
    end
    quote function $tabname(T, T2)
            $(assignexprs...)
            $tabstructname($(valsyms...))
        end
    end
end

"""
    gen_cache_struct(tab::RosenbrockTableau,cachename::Symbol,constcachename::Symbol)

Generate cache struct expression emulating those in `caches/rosenbrock_caches.jl`.
The length of k1,k2,... in the mutable cache struct is determined by the length of `tab.b`
because in the end of Rosenbrock's method, we have: `y_{n+1}=y_n+ki*bi`.
"""
function gen_cache_struct(tab::RosenbrockTableau,cachename::Symbol,constcachename::Symbol)
    kstype=[:($(Symbol(:k,i))::rateType) for i in 1:length(tab.b)]
    constcacheexpr=quote struct $constcachename{TF,UF,Tab,JType,WType,F} <: OrdinaryDiffEqConstantCache
            tf::TF
            uf::UF
            tab::Tab
            J::JType
            W::WType
            linsolve::F
        end
    end
    cacheexpr=quote
        @cache mutable struct $cachename{uType,rateType,uNoUnitsType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
            u::uType
            uprev::uType
            du::rateType
            du1::rateType
            du2::rateType
            $(kstype...)
            fsalfirst::rateType
            fsallast::rateType
            dT::rateType
            J::JType
            W::WType
            tmp::rateType
            atmp::uNoUnitsType
            tab::TabType
            tf::TFType
            uf::UFType
            linsolve_tmp::rateType
            linsolve::F
            jac_config::JCType
            grad_config::GCType
        end
    end
    constcacheexpr,cacheexpr
end

"""
    gen_algcache(cacheexpr::Expr,constcachename::Symbol,algname::Symbol,tabname::Symbol)

Generate expressions for `alg_cache(...)` emulating those in `caches/rosenbrock_caches.jl`.
"""
function gen_algcache(cacheexpr::Expr,constcachename::Symbol,algname::Symbol,tabname::Symbol)
    @capture(cacheexpr, @cache mutable struct T_ fields__ end) || error("incorrect cache expression")
    cachename = namify(T)
    ksinit = Expr[]
    valsyms = Symbol[]
    for field in fields
        if @capture(field, valsym_Symbol::valtype_)
            push!(valsyms, valsym)

            if match(r"^k[1-9]+$", String(valsym)) !== nothing
                push!(ksinit, :($valsym = zero(rate_prototype)))
            end
        end
    end

    quote
        function alg_cache(alg::$algname,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
            tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
            uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
            J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(false))
            linsolve = alg.linsolve(Val{:init},uf,u)
            $constcachename(tf,uf,$tabname(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)),J,W,linsolve)
        end
        function alg_cache(alg::$algname,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
            du = zero(rate_prototype)
            du1 = zero(rate_prototype)
            du2 = zero(rate_prototype)
            $(ksinit...)
            fsalfirst = zero(rate_prototype)
            fsallast = zero(rate_prototype)
            dT = zero(rate_prototype)
            J,W = build_J_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits,Val(true))
            tmp = zero(rate_prototype)
            atmp = similar(u, uEltypeNoUnits)
            tab = $tabname(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

            tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
            uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
            linsolve_tmp = zero(rate_prototype)
            linsolve = alg.linsolve(Val{:init},uf,u)
            grad_config = build_grad_config(alg,f,tf,du1,t)
            jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
            $cachename($(valsyms...))
        end
    end
end

"""
    gen_initialize(cachename::Symbol,constcachename::Symbol)

Generate expressions for `initialize!(...)` in `perform_step/rosenbrock_perform_step.jl`.
It only generates a default version of `initialize!` which support 3rd-order Hermite interpolation.
"""
function gen_initialize(cachename::Symbol,constcachename::Symbol)
    quote
        function initialize!(integrator, cache::$constcachename)
            integrator.kshortsize = 2
            integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
            integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t)
            integrator.destats.nf += 1

            # Avoid undefined entries if k is an array of arrays
            integrator.fsallast = zero(integrator.fsalfirst)
            integrator.k[1] = integrator.fsalfirst
            integrator.k[2] = integrator.fsallast
          end

          function initialize!(integrator, cache::$cachename)
            integrator.kshortsize = 2
            @unpack fsalfirst,fsallast = cache
            integrator.fsalfirst = fsalfirst
            integrator.fsallast = fsallast
            resize!(integrator.k, integrator.kshortsize)
            integrator.k .= [fsalfirst,fsallast]
            integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t)
            integrator.destats.nf += 1
          end
    end
end

"""
    gen_constant_perform_step(tabmask::RosenbrockTableau{Bool,Bool},cachename::Symbol,n_normalstep::Int,specialstepexpr=:nothing)

Generate non-inplace version of `perform_step!` expression emulating those in `perform_step/rosenbrock_perform_step.jl`.
The `perform_step!` function calculates `k1,k2,k3,...` defined by `(-W)ki=f(u+aij*kj,t+ci*dt)+di*dt*dT+Cij*kj*dt, i=1,2,...,n_normalstep`
and then gives the result by `y_{n+1}=y_n+ki*bi`. Terms with 0s (according to tabmask) are skipped in the expressions.
Special steps can be added before calculating `y_{n+1}`. The non-inplace `perform_step!` assumes the mass_matrix==I.
"""
function gen_constant_perform_step(tabmask::RosenbrockTableau{Bool,Bool},cachename::Symbol,n_normalstep::Int,specialstepexpr=:nothing)
    unpacktabexpr=:(@unpack ()=cache.tab)
    unpacktabexpr.args[3].args[1].args=_nonzero_vals(tabmask)
    dtCijexprs=[:($(Symbol(:dtC,Cind[1],Cind[2]))=$(Symbol(:C,Cind[1],Cind[2]))/dt) for Cind in findall(!iszero,tabmask.C)]
    dtdiexprs=[:($(Symbol(:dtd,dind))=dt*$(Symbol(:d,dind))) for dind in eachindex(tabmask.d)]
    iterexprs=[]
    for i in 1:n_normalstep
        aijkj=[:($(Symbol(:a,i+1,j))*$(Symbol(:k,j))) for j in findall(!iszero,tabmask.a[i+1,:])]
        Cijkj=[:($(Symbol(:dtC,i+1,j))*$(Symbol(:k,j))) for j in findall(!iszero,tabmask.C[i+1,:])]
        push!(iterexprs,
        quote
            $(Symbol(:k,i)) = _reshape(W\-_vec(linsolve_tmp), axes(uprev))
            integrator.destats.nsolve += 1
            u=+(uprev,$(aijkj...))
            du = f(u, p, t+$(Symbol(:c,i+1))*dt)
            integrator.destats.nf += 1
            if mass_matrix == I
                linsolve_tmp=+(du,$(Symbol(:dtd,i+1))*dT,$(Cijkj...))
            else
                linsolve_tmp=du+$(Symbol(:dtd,i+1))*dT+mass_matrix*(+($(Cijkj...)))
            end
        end)
    end
    push!(iterexprs,specialstepexpr)
    n=length(tabmask.b)
    biki=[:($(Symbol(:b,i))*$(Symbol(:k,i))) for i in 1:n]
    push!(iterexprs,
    quote
        $(Symbol(:k,n))=_reshape(W\-_vec(linsolve_tmp), axes(uprev))
        integrator.destats.nsolve += 1
        u=+(uprev,$(biki...))
    end)

    adaptiveexpr=[]
    if typeof(tabmask)<:RosenbrockAdaptiveTableau
        btildeiki=[:($(Symbol(:btilde,i))*$(Symbol(:k,i))) for i in findall(!iszero,tabmask.btilde)]
        push!(adaptiveexpr,quote
            if integrator.opts.adaptive
                utilde =  +($(btildeiki...))
                atmp = calculate_residuals(utilde, uprev, u, integrator.opts.abstol,
                                        integrator.opts.reltol,integrator.opts.internalnorm,t)
                integrator.EEst = integrator.opts.internalnorm(atmp,t)
            end
        end)
    end
    quote
        @muladd function perform_step!(integrator, cache::$cachename, repeat_step=false)
            @unpack t,dt,uprev,u,f,p = integrator
            @unpack tf,uf = cache
            $unpacktabexpr

            $(dtCijexprs...)
            $(dtdiexprs...)
            dtgamma = dt*gamma

            mass_matrix = integrator.f.mass_matrix

            # Time derivative
            tf.u = uprev
            dT = ForwardDiff.derivative(tf, t)

            W = calc_W(integrator, cache, dtgamma, repeat_step, true)
            linsolve_tmp = integrator.fsalfirst + dtd1*dT #calc_rosenbrock_differentiation!

            $(iterexprs...)

            integrator.fsallast = f(u, p, t + dt)
            integrator.destats.nf += 1

            integrator.k[1] = integrator.fsalfirst
            integrator.k[2] = integrator.fsallast
            integrator.u = u

            $(adaptiveexpr...)
        end
    end

end

"""
    gen_perform_step(tabmask::RosenbrockTableau{Bool,Bool},cachename::Symbol,n_normalstep::Int,specialstepexpr=:nothing)

Generate inplace version of `perform_step!` expression emulating those in `perform_step/rosenbrock_perform_step.jl`.
The inplace `perform_step!` produces the same result as the non-inplace version except that it treats the mass_matrix appropriately.
"""
function gen_perform_step(tabmask::RosenbrockTableau{Bool,Bool},cachename::Symbol,n_normalstep::Int,specialstepexpr=:nothing)
    unpacktabexpr=:(@unpack ()=cache.tab)
    unpacktabexpr.args[3].args[1].args=_nonzero_vals(tabmask)
    dtCij=[:($(Symbol(:dtC,"$(Cind[1])$(Cind[2])"))=$(Symbol(:C,"$(Cind[1])$(Cind[2])"))/dt) for Cind in findall(!iszero,tabmask.C)]
    dtdi=[:($(Symbol(:dtd,dind[1]))=dt*$(Symbol(:d,dind[1]))) for dind in eachindex(tabmask.d)]
    iterexprs=[]
    for i in 1:n_normalstep
        ki=Symbol(:k,i)
        dtdj=Symbol(:dtd,i+1)
        aijkj=[:($(Symbol(:a,i+1,j))*$(Symbol(:k,j))) for j in findall(!iszero,tabmask.a[i+1,:])]
        dtCijkj=[:($(Symbol(:dtC,i+1,j))*$(Symbol(:k,j))) for j in findall(!iszero,tabmask.C[i+1,:])]
        repeatstepexpr=[]
        if i==1
            repeatstepexpr=[:(!repeat_step)]
        end
        push!(iterexprs,quote
            cache.linsolve(vec($ki), W, vec(linsolve_tmp), $(repeatstepexpr...))
            @.. $ki = -$ki
            integrator.destats.nsolve += 1
            @.. u = +(uprev,$(aijkj...))
            f( du,  u, p, t+$(Symbol(:c,i+1))*dt)
            integrator.destats.nf += 1
            if mass_matrix == I
                @.. linsolve_tmp = +(du,$dtdj*dT,$(dtCijkj...))
            else
                @.. du1 = +($(dtCijkj...))
                mul!(du2,mass_matrix,du1)
                @.. linsolve_tmp = du + $dtdj*dT + du2
            end
        end)
    end
    push!(iterexprs,specialstepexpr)
    n=length(tabmask.b)
    ks=[Symbol(:k,i) for i in 1:n]
    klast=Symbol(:k,n)
    biki=[:($(Symbol(:b,i))*$(Symbol(:k,i))) for i in 1:n]
    push!(iterexprs,quote
        cache.linsolve(vec($klast), W, vec(linsolve_tmp))
        @.. $klast = -$klast
        integrator.destats.nsolve += 1
        @.. u = +(uprev,$(biki...))
    end)

    adaptiveexpr=[]
    if typeof(tabmask)<:RosenbrockAdaptiveTableau
        btildeiki=[:($(Symbol(:btilde,i))*$(Symbol(:k,i))) for i in findall(!iszero,tabmask.btilde)]
        push!(adaptiveexpr,quote
            utilde=du
            if integrator.opts.adaptive
                @.. utilde = +($(btildeiki...))
                calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol,
                                    integrator.opts.reltol,integrator.opts.internalnorm,t)
                integrator.EEst = integrator.opts.internalnorm(atmp,t)
            end
        end)
    end
    quote
        @muladd function perform_step!(integrator, cache::$cachename, repeat_step=false)
            @unpack t,dt,uprev,u,f,p = integrator
            @unpack du,du1,du2,fsallast,dT,J,W,uf,tf,$(ks...),linsolve_tmp,jac_config,atmp = cache
            $unpacktabexpr

            # Assignments
            sizeu  = size(u)
            uidx = eachindex(integrator.uprev)
            mass_matrix = integrator.f.mass_matrix

            # Precalculations
            $(dtCij...)
            $(dtdi...)
            dtgamma = dt*gamma
            calc_rosenbrock_differentiation!(integrator, cache, dtd1, dtgamma, repeat_step, true)

            $(iterexprs...)

            f( fsallast,  u, p, t + dt)
            integrator.destats.nf += 1

            $(adaptiveexpr...)
        end
    end
end

"""
    RosenbrockW6S4OSTableau()

Rahunanthan, A., & Stanescu, D. (2010). High-order W-methods.
Journal of computational and applied mathematics, 233(8), 1798-1811.
"""
function RosenbrockW6S4OSTableau()
    a=[0                  0                  0                  0                  0;
       0.5812383407115008 0                  0                  0                  0;
       0.9039624413714670 1.8615191555345010 0                  0                  0;
       2.0765797196750000 0.1884255381414796 1.8701589674910320 0                  0;
       4.4355506384843120 5.4571817986101890 4.6163507880689300 3.1181119524023610 0;
       10.791701698483260 -10.05691522584131 14.995644854284190 5.2743399543909430 1.4297308712611900]
    C=[0                  0                  0                  0                  0;
       -2.661294105131369 0                  0                  0                  0;
       -3.128450202373838 0.0000000000000000 0                  0                  0;
       -6.920335474535658 -1.202675288266817 -9.733561811413620 0                  0;
       -28.09530629102695 20.371262954793770 -41.04375275302869 -19.66373175620895 0;
       9.7998186780974000 11.935792886603180 3.6738749290132010 14.807828541095500 0.8318583998690680]
    b=[6.4562170746532350,-4.853141317768053,9.7653183340692600,2.0810841772787230,0.6603936866352417,0.6000000000000000]
    gamma=0.2500000000000000
    d=[0.2500000000000000,0.0836691184292894,0.0544718623516351,-0.3402289722355864,0.0337651588339529,-0.0903074267618540]
    c=[0                 ,0.1453095851778752,0.3817422770256738,0.6367813704374599,0.7560744496323561,0.9271047239875670]
    RosenbrockFixedTableau(a,C,b,gamma,d,c)
end

"""
    @RosenbrockW6S4OS(part)

Generate code for the RosenbrockW6S4OS method.
`part` should be one of `:tableau`, `:cache`, `:init`, `:performstep`.
`@RosenbrockW6S4OS(:tableau)` should be placed in `tableaus/rosenbrock_tableaus.jl`.
`@RosenbrockW6S4OS(:cache)` should be placed in `caches/rosenbrock_caches.jl`.
`@RosenbrockW6S4OS(:init)` and `@RosenbrockW6S4OS(:performstep)` should be
placed in `perform_step/rosenbrock_perform_step.jl`.
"""
macro RosenbrockW6S4OS(part)
    tab=RosenbrockW6S4OSTableau()
    tabmask=_masktab(tab)
    algname=:RosenbrockW6S4OS
    tabname=:RosenbrockW6S4OSTableau
    tabstructname=:RosenbrockW6STableau
    cachename=:RosenbrockW6SCache
    constcachename=:RosenbrockW6SConstantCache
    n_normalstep=length(tab.b)-1
    if part.value==:tableau
        #println("Generating const cache")
        tabstructexpr=gen_tableau_struct(tabmask,tabstructname)
        tabexpr=gen_tableau(tab,tabstructexpr,tabname)
        return esc(quote $([tabstructexpr,tabexpr]...) end)
    elseif part.value==:cache
        #println("Generating cache")
        constcacheexpr,cacheexpr=gen_cache_struct(tabmask,cachename,constcachename)
        algcacheexpr=gen_algcache(cacheexpr,constcachename,algname,tabname)
        return esc(quote $([constcacheexpr,cacheexpr,algcacheexpr]...) end)
    elseif part.value==:init
        #println("Generating initialize")
        return esc(gen_initialize(cachename,constcachename))
    elseif part.value==:performstep
        #println("Generating perform_step")
        constperformstepexpr=gen_constant_perform_step(tabmask,constcachename,n_normalstep)
        performstepexpr=gen_perform_step(tabmask,cachename,n_normalstep)
        return esc(quote $([constperformstepexpr,performstepexpr]...) end)
    else
        throw(ArgumentError("Unknown parameter!"))
        nothing
    end
end

"""
    Ros4dummyTableau()

Generate a dummy tableau for ROS4 methods. It can be considered as performing elementwise OR to the masks
of those sepecific tableaus: `Ros4dummyTableau()==_masktab(RosShamp4Tableau()) OR _masktab(Veldd4Tableau()) OR ...`
ROS4 methods have the property of a4j==a3j so a is a 3*3 matrix instead of a 4*4 matrix and c is a 1*3 vector instead of a 1*4 vector.
"""
function Ros4dummyTableau()#create a tabmask for all ROS4 methods where false->0,true->non-0
    a=[false false false;
       true  false false;
       true  true  false]
    C=[false false false false;
       true  false false false;
       true  true  false false;
       true  true  true  false]
    b=[true,true,true,true]
    btilde=[true,true,true,true]
    gamma=true
    c=[false,true,true]
    d=[true,true,true,true]
    RosenbrockAdaptiveTableau(a,C,b,btilde,gamma,d,c)
end

"""
    RosShamp4Tableau()

L. F. Shampine, Implementation of Rosenbrock Methods,
ACM Transactions on Mathematical Software (TOMS), 8: 2, 93-113.
doi:10.1145/355993.355994
"""
function RosShamp4Tableau()
    a=[0      0     0;
       2      0     0;
       48//25 6//25 0]
    C=[ 0         0        0    0;
       -8         0        0    0;
        372//25   12//5    0    0;
       -112//125 -54//125 -2//5 0]
    b=[19//9,1//2,25//108,125//108]
    btilde=[17//54,7//36,0,125//108]
    gamma=1//2
    c=[0,1,3//5]
    d=[1//2,-3//2,242//100,116//1000]#2.42,0.116
    RosenbrockAdaptiveTableau(a,C,b,btilde,gamma,d,c)
end

"""
    Veldd4Tableau()

van Veldhuizen, D-stability and Kaps-Rentrop-methods,
M. Computing (1984) 32: 229. doi:10.1007/BF02243574
"""
function Veldd4Tableau()
    a=[0                 0                 0;
       2                 0                 0;
       4.812234362695436 4.578146956747842 0]
    C=[ 0                  0                  0                 0;
       -5.333333333333331  0                  0                 0;
        6.100529678848254  1.804736797378427  0                 0;
       -2.540515456634749 -9.443746328915205 -1.988471753215993 0]
    b=[4.289339254654537,5.036098482851414,0.6085736420673917,1.355958941201148]
    btilde=[2.175672787531755,2.950911222575741,-0.7859744544887430,-1.355958941201148]
    gamma=0.2257081148225682
    c=[0,0.4514162296451364,0.8755928946018455]
    d=[0.2257081148225682,-0.04599403502680582,0.5177590504944076,-0.03805623938054428]
    RosenbrockAdaptiveTableau(a,C,b,btilde,gamma,d,c)
end

"""
    Velds4Tableau()

van Veldhuizen, D-stability and Kaps-Rentrop-methods,
M. Computing (1984) 32: 229. doi:10.1007/BF02243574
"""
function Velds4Tableau()
    a=[0    0    0;
       2    0    0;
       7//4 1//4 0]
    C=[ 0     0    0 0;
       -8     0    0 0;
       -8    -1    0 0;
        1//2 -1//2 2 0]
    b=[4//3,2//3,-4//3,4//3]
    btilde=[-1//3,-1//3,0,-4//3]
    gamma=1//2
    c=[0,1,1//2]
    d=[1//2,-3//2,-3//4,1//4]
    RosenbrockAdaptiveTableau(a,C,b,btilde,gamma,d,c)
end

"""
    GRK4TTableau()

Kaps, P. & Rentrop, Generalized Runge-Kutta methods of order four
with stepsize control for stiff ordinary differential equations.
P. Numer. Math. (1979) 33: 55. doi:10.1007/BF01396495
"""
function GRK4TTableau()
    a=[0                 0                 0;
       2                 0                 0;
       4.524708207373116 4.163528788597648 0]
    C=[ 0                  0                   0                 0;
       -5.071675338776316  0                   0                 0;
        6.020152728650786  0.1597506846727117  0                 0;
       -1.856343618686113 -8.505380858179826  -2.084075136023187 0]
    b=[3.957503746640777,4.624892388363313,0.6174772638750108,1.282612945269037]
    btilde=[2.302155402932996,3.073634485392623,-0.8732808018045032,-1.282612945269037]
    gamma=0.231
    c=[0,0.462,0.8802083333333334]
    d=[0.2310000000000000,-0.03962966775244303,0.5507789395789127,-0.05535098457052764]
    RosenbrockAdaptiveTableau(a,C,b,btilde,gamma,d,c)
end

"""
    GRK4ATableau()

Kaps, P. & Rentrop, Generalized Runge-Kutta methods of order four
with stepsize control for stiff ordinary differential equations.
P. Numer. Math. (1979) 33: 55. doi:10.1007/BF01396495
"""
function GRK4ATableau()
    a=[0                 0                  0;
       1.108860759493671 0                  0;
       2.377085261983360 0.1850114988899692 0]
    C=[ 0                 0                  0                 0;
       -4.920188402397641 0                  0                 0;
        1.055588686048583 3.351817267668938  0                 0;
        3.846869007049313 3.427109241268180 -2.162408848753263 0]
    b=[1.845683240405840,0.1369796894360503,0.7129097783291559,0.6329113924050632]
    btilde=[0.04831870177201765,-0.6471108651049505,0.2186876660500240,-0.6329113924050632]
    gamma=0.395
    c=[0,0.438,0.87]
    d=[0.395,-0.3726723954840920,0.06629196544571492,0.4340946962568634]
    RosenbrockAdaptiveTableau(a,C,b,btilde,gamma,d,c)
end

"""
    Ros4LSTableau()

E. Hairer, G. Wanner, Solving ordinary differential equations II,
stiff and differential-algebraic problems. Computational mathematics (2nd revised ed.), Springer (1996)
"""
function Ros4LSTableau()
    a=[0                 0                  0;
       2                 0                  0;
       1.867943637803922 0.2344449711399156 0]
    C=[ 0                  0                   0                  0;
       -7.137615036412310  0                   0                  0;
        2.580708087951457  0.6515950076447975  0                  0;
       -2.137148994382534 -0.3214669691237626 -0.6949742501781779 0]
    b=[2.255570073418735,0.2870493262186792,0.4353179431840180,1.093502252409163]
    btilde=[-0.2815431932141155,-0.07276199124938920,-0.1082196201495311,-1.093502252409163]
    gamma=0.5728200000000000
    c=[0,1.145640000000000,0.6552168638155900]
    d=[0.5728200000000000,-1.769193891319233,0.7592633437920482,-0.1049021087100450]
    RosenbrockAdaptiveTableau(a,C,b,btilde,gamma,d,c)
end

"""
    @Rosenbrock4(part)

Generate code for the Rosenbrock4 methods: RosShamp4, Veldd4, Velds4, GRK4A, GRK4T, Ros4LStab.
`part` should be one of `:tableau`, `:cache`, `:performstep`.
`@Rosenbrock4(:tableau)` should be placed in `tableaus/rosenbrock_tableaus.jl`.
`@Rosenbrock4(:cache)` should be placed in `caches/rosenbrock_caches.jl`.
`@Rosenbrock4(:performstep)` should be placed in `perform_step/rosenbrock_perform_step.jl`.
The `initialize!` function for Rosenbrock4 methods is already included in `rosenbrock_perform_step.jl`.
The special property of ROS4 methods that a4j==a3j requires a special step in `perform_step!` that
calculates `linsolve_tmp` from the previous `du` which reduces a function call.
"""
macro Rosenbrock4(part)
    tabmask=Ros4dummyTableau()#_masktab(tab)
    cachename=:Rosenbrock4Cache
    constcachename=:Rosenbrock4ConstantCache
    RosShamp4tabname=:RosShamp4Tableau
    Veldd4tabname=:Veldd4Tableau
    Velds4tabname=:Velds4Tableau
    GRK4Ttabname=:GRK4TTableau
    GRK4Atabname=:GRK4ATableau
    Ros4LStabname=:Ros4LStabTableau
    n_normalstep=2 #for the third step a4j=a3j which reduced one function call
    if part.value==:tableau
        #println("Generating tableau for Rosenbrock4")
        tabstructexpr=gen_tableau_struct(tabmask,:Ros4Tableau)
        tabexprs=Array{Expr,1}()
        push!(tabexprs,tabstructexpr)
        push!(tabexprs,gen_tableau(RosShamp4Tableau(),tabstructexpr,RosShamp4tabname))
        push!(tabexprs,gen_tableau(Veldd4Tableau(),tabstructexpr,Veldd4tabname))
        push!(tabexprs,gen_tableau(Velds4Tableau(),tabstructexpr,Velds4tabname))
        push!(tabexprs,gen_tableau(GRK4TTableau(),tabstructexpr,GRK4Ttabname))
        push!(tabexprs,gen_tableau(GRK4ATableau(),tabstructexpr,GRK4Atabname))
        push!(tabexprs,gen_tableau(Ros4LSTableau(),tabstructexpr,Ros4LStabname))
        return esc(quote $(tabexprs...) end)
    elseif part.value==:cache
        #println("Generating cache for Rosenbrock4")
        constcacheexpr,cacheexpr=gen_cache_struct(tabmask,cachename,constcachename)
        cacheexprs=Array{Expr,1}([constcacheexpr,cacheexpr])
        push!(cacheexprs,gen_algcache(cacheexpr,constcachename,:RosShamp4,RosShamp4tabname))
        push!(cacheexprs,gen_algcache(cacheexpr,constcachename,:Veldd4,Veldd4tabname))
        push!(cacheexprs,gen_algcache(cacheexpr,constcachename,:Velds4,Velds4tabname))
        push!(cacheexprs,gen_algcache(cacheexpr,constcachename,:GRK4T,GRK4Ttabname))
        push!(cacheexprs,gen_algcache(cacheexpr,constcachename,:GRK4A,GRK4Atabname))
        push!(cacheexprs,gen_algcache(cacheexpr,constcachename,:Ros4LStab,Ros4LStabname))
        return esc(quote $(cacheexprs...) end)
    elseif part.value==:performstep
        #println("Generating perform_step for Rosenbrock4")
        specialstepconst=quote
            k3 = _reshape(W\-_vec(linsolve_tmp), axes(uprev))
            integrator.destats.nsolve += 1
            #u = uprev  + a31*k1 + a32*k2 #a4j=a3j
            #du = f(u, p, t+c3*dt) #reduced function call
            if mass_matrix == I
                linsolve_tmp =  du + dtd4*dT + dtC41*k1 + dtC42*k2 + dtC43*k3
            else
                linsolve_tmp = du + dtd4*dT + mass_matrix * (dtC41*k1 + dtC42*k2 + dtC43*k3)
            end
        end
        specialstep=quote
            cache.linsolve(vec(k3), W, vec(linsolve_tmp))
            @.. k3 = -k3
            integrator.destats.nsolve += 1
            #@.. u = uprev + a31*k1 + a32*k2 #a4j=a3j
            #f( du,  u, p, t+c3*dt) #reduced function call
            if mass_matrix == I
                @.. linsolve_tmp = du + dtd4*dT + dtC41*k1 + dtC42*k2 + dtC43*k3
            else
                @.. du1 = dtC41*k1 + dtC42*k2 + dtC43*k3
                mul!(du2,mass_matrix,du1)
                @.. linsolve_tmp = du + dtd4*dT + du2
            end
        end
        constperformstepexpr=gen_constant_perform_step(tabmask,constcachename,n_normalstep,specialstepconst)
        performstepexpr=gen_perform_step(tabmask,cachename,n_normalstep,specialstep)
        return esc(quote $([constperformstepexpr,performstepexpr]...) end)
    else
        throw(ArgumentError("Unknown parameter!"))
        nothing
    end
end

#ROS34PW methods (Rang and Angermann, 2005)
"""
    Ros34dummyTableau()

Generate a dummy tableau for ROS34W methods proposed by Rang and Angermann. This type of methods has 4 steps.
"""
function Ros34dummyTableau()
    a=[false false false false;
       true  false false false;
       true  true  false false;
       true  true  true  false]
    C=[false false false false;
       true  false false false;
       true  true  false false;
       true  true  true  false]
    b=[true,true,true,true]
    btilde=[true,true,true,true]
    gamma=true
    c=[false,true,true,true]
    d=[true,true,true,true]
    RosenbrockAdaptiveTableau(a,C,b,btilde,gamma,d,c)
end

"""
    _transformtab(Alpha,Gamma,B,Bhat)

Transform the tableau from values in the paper into values used in OrdinaryDiffEq according to p112 in Hairer and Wanner.

E. Hairer, G. Wanner, Solving ordinary differential equations II, stiff and
differential-algebraic problems. Computational mathematics (2nd revised ed.), Springer (1996)
"""
function _transformtab(Alpha,Gamma,B,Bhat)
    invGamma=inv(Gamma)
    a=Alpha*invGamma
    C=diagm(0=>diag(invGamma))-invGamma
    b=[(transpose(B)*invGamma)...]# [2Darray...]=>1Darray
    btilde=[(transpose(B-Bhat)*invGamma)...]
    gamma=Gamma[1,1]#Gamma11==Gamma22==...==Gammass
    d=[sum(Gamma,dims=2)...]#di=sum_j Gamma_ij
    c=[sum(Alpha,dims=2)...]#ci=sum_j Alpha_ij
    (a,C,b,btilde,d,c)
end

"""
    ROS34PW1aTableau()

L-Stable, Rosenbrock-W method with order 3 and 4 inner steps

Rang, J., & Angermann, L. (2005). New Rosenbrock W-methods of order 3 for partial
differential algebraic equations of index 1. BIT Numerical Mathematics, 45(4), 761-787.
"""
function ROS34PW1aTableau()
    gamma=4.358665215084590e-1
    Alpha=[0                 0                    0   0;
           2.218787467653286 0                    0   0;
           0                 0                    0   0; # can reduce one function call with specialized perform_step
           1.208587690772214 7.511610241919324e-2 0.5 0]
    Gamma=[ gamma                 0                    0     0;
           -2.218787467653286     gamma                0     0;
           -9.461966143940745e-2 -7.913526735718213e-3 gamma 0;
           -1.870323744195384    -9.624340112825115e-2 2.726301276675511e-1 gamma]
    B=[3.285609536316354e-1,-5.785609536316354e-1,0.25,1]
    Bhat=[-0.25,0,0.25,1]#B-Bhat[3:4]==[0,0]
    a,C,b,btilde,d,c=_transformtab(Alpha,Gamma,B,Bhat)
    RosenbrockAdaptiveTableau(a,C,b,btilde,gamma,d,c)
end

"""
    ROS34PW1bTableau()

L-Stable, Rosenbrock-W method with order 3 and 4 inner steps

Rang, J., & Angermann, L. (2005). New Rosenbrock W-methods of order 3 for partial
differential algebraic equations of index 1. BIT Numerical Mathematics, 45(4), 761-787.
"""
function ROS34PW1bTableau()
    gamma=4.358665215084590e-1
    Alpha=[0                 0 0   0;
           2.218787467653286 0 0   0;
           2.218787467653286 0 0   0; # can reduce one function call with specialized perform_step
           1.453923375357884 0 0.1 0]
    Gamma=[ gamma              0                    0     0;
           -2.218787467653286  gamma                0     0;
           -2.848610224639349 -5.267530183845237e-2 gamma 0;
           -1.128167857898393 -1.677546870499461e-1 5.452602553351021e-2 gamma]
    B=[5.495647928937977e-1,-5.507258170857301e-1,0.25,7.511610241919324e-1]
    Bhat=[-1.161024191932427e-3,0,0.25,7.511610241919324e-1]#B-Bhat[3:4]==[0,0]
    a,C,b,btilde,d,c=_transformtab(Alpha,Gamma,B,Bhat)
    RosenbrockAdaptiveTableau(a,C,b,btilde,gamma,d,c)
end

"""
    ROS34PW2Tableau()

A stiffy accurate Rosenbrock-W method with order 3 and 4 inner steps whose
embedded method is strongly A-stable with Rinf~=0.48

Rang, J., & Angermann, L. (2005). New Rosenbrock W-methods of order 3 for partial
differential algebraic equations of index 1. BIT Numerical Mathematics, 45(4), 761-787.
"""
function ROS34PW2Tableau()
    gamma=4.3586652150845900e-1
    Alpha=[0                      0                     0 0;
           8.7173304301691801e-1  0                     0 0;
           8.4457060015369423e-1 -1.1299064236484185e-1 0 0;
           0                      0                     1 0]
    Gamma=[ gamma                  0                     0     0;
           -8.7173304301691801e-1  gamma                 0     0;
           -9.0338057013044082e-1  5.4180672388095326e-2 gamma 0;
            2.4212380706095346e-1 -1.2232505839045147    5.4526025533510214e-1 gamma]
    B=[2.4212380706095346e-1,-1.2232505839045147,1.5452602553351020,4.3586652150845900e-1]
    Bhat=[3.7810903145819369e-1,-9.6042292212423178e-2,0.5,2.1793326075422950e-1]
    a,C,b,btilde,d,c=_transformtab(Alpha,Gamma,B,Bhat)
    RosenbrockAdaptiveTableau(a,C,b,btilde,gamma,d,c)
end

"""
    ROS34PW3Tableau()

an A-stable (Rinf~=0.63), Rosenbrock-W method with order 4 and 4 inner steps.

Rang, J., & Angermann, L. (2005). New Rosenbrock W-methods of order 3 for partial
differential algebraic equations of index 1. BIT Numerical Mathematics, 45(4), 761-787.
"""
function ROS34PW3Tableau()#4th order
    gamma=1.0685790213016289
    Alpha=[0                      0                     0 0;
           2.5155456020628817     0                     0 0;
           5.0777280103144085e-1  0.75                  0 0;
           1.3959081404277204e-1 -3.3111001065419338e-1 8.2040559712714178e-1 0]
    Gamma=[ gamma                  0                      0     0;
           -2.5155456020628817     gamma                  0     0;
           -8.7991339217106512e-1 -9.6014187766190695e-1  gamma 0;
           -4.1731389379448741e-1  4.1091047035857703e-1 -1.3558873204765276 gamma]
    B=[2.2047681286931747e-1,2.7828278331185935e-3,7.1844787635140066e-3,7.6955588053404989e-1]
    Bhat=[3.1300297285209688e-1,-2.8946895245112692e-1,9.7646597959903003e-1,0]
    a,C,b,btilde,d,c=_transformtab(Alpha,Gamma,B,Bhat)
    RosenbrockAdaptiveTableau(a,C,b,btilde,gamma,d,c)
end

"""
    @ROS34PW(part)

Generate code for the ROS34PW methods: ROS34PW1a, ROS34PW1b, ROS34PW2, ROS34PW3.
`part` should be one of `:tableau`, `:cache`, `:init`, `:performstep`.
`@ROS34PW(:tableau)` should be placed in `tableaus/rosenbrock_tableaus.jl`.
`@ROS34PW(:cache)` should be placed in `caches/rosenbrock_caches.jl`.
`@ROS34PW(:init)` and `@ROS34PW(:performstep)` should be placed in
`perform_step/rosenbrock_perform_step.jl`.
"""
macro ROS34PW(part)
    tabmask=Ros34dummyTableau()
    cachename=:ROS34PWCache
    constcachename=:ROS34PWConstantCache
    ROS34PW1atabname=:ROS34PW1aTableau
    ROS34PW1btabname=:ROS34PW1bTableau
    ROS34PW2tabname=:ROS34PW2Tableau
    ROS34PW3tabname=:ROS34PW3Tableau
    n_normalstep=length(tabmask.b)-1
    if part.value==:tableau
        tabstructexpr=gen_tableau_struct(tabmask,:Ros34Tableau)
        tabexprs=Array{Expr,1}([tabstructexpr])
        push!(tabexprs,gen_tableau(ROS34PW1aTableau(),tabstructexpr,ROS34PW1atabname))
        push!(tabexprs,gen_tableau(ROS34PW1bTableau(),tabstructexpr,ROS34PW1btabname))
        push!(tabexprs,gen_tableau(ROS34PW2Tableau(),tabstructexpr,ROS34PW2tabname))
        push!(tabexprs,gen_tableau(ROS34PW3Tableau(),tabstructexpr,ROS34PW3tabname))
        return esc(quote $(tabexprs...) end)
    elseif part.value==:cache
        constcacheexpr,cacheexpr=gen_cache_struct(tabmask,cachename,constcachename)
        cacheexprs=Array{Expr,1}([constcacheexpr,cacheexpr])
        push!(cacheexprs,gen_algcache(cacheexpr,constcachename,:ROS34PW1a,ROS34PW1atabname))
        push!(cacheexprs,gen_algcache(cacheexpr,constcachename,:ROS34PW1b,ROS34PW1btabname))
        push!(cacheexprs,gen_algcache(cacheexpr,constcachename,:ROS34PW2,ROS34PW2tabname))
        push!(cacheexprs,gen_algcache(cacheexpr,constcachename,:ROS34PW3,ROS34PW3tabname))
        return esc(quote $(cacheexprs...) end)
    elseif part.value==:init
        return esc(gen_initialize(cachename,constcachename))
    elseif part.value==:performstep
        performstepexprs=Array{Expr,1}()
        push!(performstepexprs,gen_constant_perform_step(tabmask,constcachename,n_normalstep))
        push!(performstepexprs,gen_perform_step(tabmask,cachename,n_normalstep))
        return esc(quote $(performstepexprs...) end)
    else
        throw(ArgumentError("Unknown parameter!"))
        nothing
    end
end

#=========================================================================================
# How to add a new method
1. `OrdinaryDiffEq.jl`: export <Algorithm_name>
2. `alg_utils.jl`: alg_order(alg::<Algorithm_name>)=<Algorithm_order>
    if the method is a W-method, add isWmethod(alg::<Algorithm_name>) = true as well
3. `algorithms.jl`: add algorithm struct
4. `generic_rosenbrock.jl`:
    a. write dummy tableau function (or generate from actual tableau using _masktab) for generating
       table sturct, cache struct and perform_step
    b. write tableau function. When only `Alpha, Gamma, B, Bhat` are given, use _transformtab
    c. write macro with :tableau, :cache, :init and :performstep
    d. put the macros in the right places.

# How to refactor methods into generic ones
RUN CONVERGENCE TESTS BETWEEN ANY OF THE TWO STEPS!
1. write tableau function and macro definition in this file
2. replace the tableau function (usually named with `XXXConstCache()`) using `gen_tableau()`
   and the original tableau struct expression in `tableaus/rosenbrock_tableaus.jl`
3. replace the `perform_step!` methods in `perform_step/rosenbrock_perform_step.jl` using
   `gen_perform_step()` and `gen_constant_perform_step()`
4. replace cache struct and `alg_cache` in `caches/rosenbrock_caches.jl` using `gen_cache_struct()`
   and `gen_algcache()`
5. If the method only have 3rd-order Hermite interpolation, you can replace `initialize!()`
   in `perform_step/rosenbrock_perform_step.jl` with `gen_initialize()`
DONE

# How to debug
Use macroexpand like `macroexpand(OrdinaryDiffEq,:(@ROS34PW(:performstep)))` and check the
generated codes.

`Revise.jl` is not compatible with macros. One may want to manually re-eval files that use
the macro like `@eval OrdinaryDiffEq include(...)`

# If you want to refactor Rosenbrock methods ...
You need to change respective places in this file.
1. `perform_step/rosenbrock_perform_step.jl` -> `gen_perform_step()`, `gen_constant_perform_step()`,
    `gen_initialize()` and special step expressions in macro definitions
2. `caches/rosenbrock_caches.jl` ->  `gen_algcache()`, `gen_cache_struct()`
3. `tableaus/rosenbrock_tableaus.jl` -> `gen_tableau_struct()` and `gen_tableau()`
=========================================================================================#

ext: jl


text:
"""
    NoncentralChisq(ν, λ)

The *noncentral chi-squared distribution* with `ν` degrees of freedom and noncentrality parameter `λ` has the probability density function

```math
f(x; \\nu, \\lambda) = \\frac{1}{2} e^{-(x + \\lambda)/2} \\left( \\frac{x}{\\lambda} \\right)^{\\nu/4-1/2} I_{\\nu/2-1}(\\sqrt{\\lambda x}), \\quad x > 0
```

It is the distribution of the sum of squares of `ν` independent [`Normal`](@ref) variates with individual means ``\\mu_i`` and

```math
\\lambda = \\sum_{i=1}^\\nu \\mu_i^2
```

```julia
NoncentralChisq(ν, λ)     # Noncentral chi-squared distribution with ν degrees of freedom and noncentrality parameter λ

params(d)    # Get the parameters, i.e. (ν, λ)
```

External links

* [Noncentral chi-squared distribution on Wikipedia](https://en.wikipedia.org/wiki/Noncentral_chi-squared_distribution)
"""
struct NoncentralChisq{T<:Real} <: ContinuousUnivariateDistribution
    ν::T
    λ::T
    NoncentralChisq{T}(ν::T, λ::T) where {T <: Real} = new{T}(ν, λ)
end

function NoncentralChisq(ν::T, λ::T) where {T <: Real}
    @check_args(NoncentralChisq, ν > zero(ν) && λ >= zero(λ))
    return NoncentralChisq{T}(ν, λ)
end

NoncentralChisq(ν::T, λ::T, ::NoArgCheck) where {T<:Real} = NoncentralChisq{T}(ν, λ)
NoncentralChisq(ν::Real, λ::Real) = NoncentralChisq(promote(ν, λ)...)
NoncentralChisq(ν::Integer, λ::Integer) = NoncentralChisq(float(ν), float(λ))

@distr_support NoncentralChisq 0.0 Inf

#### Conversions

function convert(::Type{NoncentralChisq{T}}, ν::S, λ::S) where {T <: Real, S <: Real}
    NoncentralChisq(T(ν), T(λ))
end
function convert(::Type{NoncentralChisq{T}}, d::NoncentralChisq{S}) where {T <: Real, S <: Real}
    NoncentralChisq(T(d.ν), T(d.λ), NoArgCheck())
end

### Parameters

params(d::NoncentralChisq) = (d.ν, d.λ)
@inline partype(d::NoncentralChisq{T}) where {T<:Real} = T

### Statistics

mean(d::NoncentralChisq) = d.ν + d.λ
var(d::NoncentralChisq) = 2(d.ν + 2d.λ)
skewness(d::NoncentralChisq) = 2sqrt2*(d.ν + 3d.λ)/sqrt(d.ν + 2d.λ)^3
kurtosis(d::NoncentralChisq) = 12(d.ν + 4d.λ)/(d.ν + 2d.λ)^2

function mgf(d::NoncentralChisq, t::Real)
    exp(d.λ * t/(1 - 2t))*(1 - 2t)^(-d.ν/2)
end

function cf(d::NoncentralChisq, t::Real)
    cis(d.λ * t/(1 - 2im*t))*(1 - 2im*t)^(-d.ν/2)
end

### Evaluation & Sampling

@_delegate_statsfuns NoncentralChisq nchisq ν λ

# TODO: remove RFunctions dependency
@rand_rdist(NoncentralChisq)
rand(d::NoncentralChisq) = StatsFuns.RFunctions.nchisqrand(d.ν, d.λ)


text:
"""
    ms_to_mol(G,T,P,R,K₀)
    ms_to_mol(G,T,P)

Conversion of a conductance `G` from ``m\\ s^{-1}`` to ``mol\\ m^{-2}\\ s^{-1}``.

# Arguments

- `G` (``m\\ s^{-1}``): conductance
- `T` (°C): air temperature
- `P` (kPa): air pressure
- `R` (``J\\ mol^{-1}\\ K^{-1}``): universal gas constant.
- `K₀` (°C): absolute zero

# See also

[`mol_to_ms`](@ref) for the inverse process.
"""
function ms_to_mol(G,T,P,R,K₀)
    G * f_ms_to_mol(T,P,R,K₀)
end

function ms_to_mol(G,T,P)
    constants = Constants()
    ms_to_mol(G,T,P,constants.R,constants.K₀)
end

"""
    ms_to_mol(G,T,P,R,K₀)
    ms_to_mol(G,T,P)

Conversion of a conductance `G` from ``mol\\ m^{-2}\\ s^{-1}`` to ``m\\ s^{-1}``.

# Arguments

- `G` (``m\\ s^{-1}``): conductance
- `T` (°C): air temperature
- `P` (kPa): air pressure
- `R` (``J\\ mol^{-1}\\ K^{-1}``): universal gas constant.
- `K₀` (°C): absolute zero

# See also

[`ms_to_mol`](@ref) for the inverse process.
"""
function mol_to_ms(G,T,P,R,K₀)
    G / f_ms_to_mol(T,P,R,K₀)
end

function mol_to_ms(G,T,P)
    constants = Constants()
    mol_to_ms(G,T,P,constants.R,constants.K₀)
end

"""
Conversion factor between conductance in ``m\\ s^{-1}`` to ``mol\\ m^{-2}\\ s^{-1}``.

# Arguments

- `T` (°C): air temperature
- `P` (kPa): air pressure
- `R` (``J\\ mol^{-1}\\ K^{-1}``): universal gas constant.
- `K₀` (°C): absolute zero
"""
function f_ms_to_mol(T,P,R,K₀)
    (P * 1000) / (R * (T - K₀))
end

"""
    rh_from_vpd(VPD,eₛ)

Conversion between VPD and rh.

# Examples

```julia
eₛ = e_sat(Tₐ)
rh_from_vpd(1.5,eₛ)
```
"""
function rh_from_vpd(VPD,eₛ)
    one(VPD) - VPD / eₛ
end

"""
    rh_from_e(VPD,eₛ)

Conversion between e (kPa) and rh (0-1).

# Examples

```julia
rh_from_e(1.5,25.0)
```
"""
function rh_from_e(e,Tₐ)
    eₛ = e_sat(Tₐ)
    min(one(e), e/eₛ)
end

"""
    gbh_to_gbw(gbh, Gbₕ_to_Gbₕ₂ₒ = Constants().Gbₕ_to_Gbₕ₂ₒ)
    gbw_to_gbh(gbh, Gbₕ_to_Gbₕ₂ₒ = Constants().Gbₕ_to_Gbₕ₂ₒ)

Boundary layer conductance for water vapor from boundary layer conductance for heat.

# Arguments

- `gbh` (m s-1): boundary layer conductance for heat under mixed convection.
- `Gbₕ_to_Gbₕ₂ₒ`: conversion factor.

# Note

Gbₕ is the sum of free and forced convection. See [`gbₕ_free`](@ref) and [`gbₕ_forced`](@ref).
"""
function gbh_to_gbw(gbh, Gbₕ_to_Gbₕ₂ₒ = Constants().Gbₕ_to_Gbₕ₂ₒ)
    gbh * Gbₕ_to_Gbₕ₂ₒ
end

function gbw_to_gbh(gbh, Gbₕ_to_Gbₕ₂ₒ = Constants().Gbₕ_to_Gbₕ₂ₒ)
    gbh / Gbₕ_to_Gbₕ₂ₒ
end

"""
    gsc_to_gsw(Gₛ, Gsc_to_Gsw = Constants().Gsc_to_Gsw)

Conversion of a stomatal conductance for CO₂ into stomatal conductance for H₂O.
"""
function gsc_to_gsw(Gₛ, Gsc_to_Gsw = Constants().Gsc_to_Gsw)
    Gₛ * Gsc_to_Gsw
end

"""
    gsw_to_gsc(Gₛ, Gsc_to_Gsw = Constants().Gsc_to_Gsw)

Conversion of a stomatal conductance for H₂O into stomatal conductance for CO₂.
"""
function gsw_to_gsc(Gₛ, Gsc_to_Gsw = Constants().Gsc_to_Gsw)
    Gₛ / Gsc_to_Gsw
end


text:
using AutoHashEquals: @auto_hash_equals
using Bijections: Bijection, inverse
using DataStructures: DefaultDict
using Dates

# This is the identifier with which every tea file must start. It is used to ensure that
# the file is a valid tea-file, and also that the endianness matches.
const MAGIC_VALUE = Int64(0x0d0e0a0402080500)

abstract type AbstractSection end

struct TeaFileMetadata
    item_start::Int64  # Byte-index of the start of the item area.
    item_end::Int64  # Byte-index for the end of the item area, or 0 for EOF.
    sections::Vector{AbstractSection}

    function TeaFileMetadata(item_start, item_end, sections)
        metadata = new(item_start, item_end, sections)
        _verify_metadata(metadata)
        return metadata
    end
end

function Base.hash(x::TeaFileMetadata, h::UInt)
    return hash(x.item_start, hash(x.item_end, hash(x.sections, hash(:TeaFileMetadata, h))))
end

function Base.:(==)(a::TeaFileMetadata, b::TeaFileMetadata)
    return (
        a.item_start == b.item_start &&
        a.item_end == b.item_end &&
        a.sections == b.sections
    )
end

const _FIELD_DATA_TYPE_TO_ID = Bijection(Dict{DataType,Int32}(
    Int8 => 1,
    Int16 => 2,
    Int32 => 3,
    Int64 => 4,
    UInt8 => 5,
    UInt16 => 6,
    UInt32 => 7,
    UInt64 => 8,
    Float32 => 9,
    Float64 => 10,
    # We do not support .NET decimals, but this reserves 0x200.
    # User-defined types must have identifiers at or above 0x1000
))

field_type_id(T::Type) = _FIELD_DATA_TYPE_TO_ID[T]

struct Field
    type_id::Int32
    offset::Int32  # Byte offset inside the item
    name::String
end

field_type(field::Field) = inverse(_FIELD_DATA_TYPE_TO_ID, field.type_id)

struct NameValue{T <: Union{Int32,Float64,String,Base.UUID}}
    name::String
    value::T
end

const _NAME_VALUE_TYPE_TO_KIND = Bijection(Dict{DataType, Int32}(
    Int32 => 1,
    Float64 => 2,
    String => 3,
    Base.UUID => 4,
))

kind(::NameValue{T}) where {T} = _NAME_VALUE_TYPE_TO_KIND[T]
name_value_type(kind::Int32) = inverse(_NAME_VALUE_TYPE_TO_KIND, kind)

@auto_hash_equals struct ItemSection <: AbstractSection
    item_size::Int32  # Size of each item in bytes.
    item_name::String
    fields::Vector{Field}
end

@auto_hash_equals struct TimeSection <: AbstractSection
    epoch::Int64  # Number of days from 0001-01-01 to the origin.
    ticks_per_day::Int64
    time_field_offsets::Vector{Int32}
end

struct ContentDescriptionSection <: AbstractSection
    description::String
end

@auto_hash_equals struct NameValueSection <: AbstractSection
    name_values::Vector{NameValue}
end

const _SECTION_TYPE_TO_ID = Bijection(Dict{DataType, Int32}(
    ItemSection => 0x0a,
    TimeSection => 0x40,
    ContentDescriptionSection => 0x80,
    NameValueSection => 0x81,
))

section_id(::T) where {T <: AbstractSection} = _SECTION_TYPE_TO_ID[T]
function section_type(id::Int32)::(Type{T} where {T <: AbstractSection})
    return inverse(_SECTION_TYPE_TO_ID, id)
end

"""Get a section of the specified type."""
function get_section(metadata::TeaFileMetadata, type::Type)
    return only(filter(section -> isa(section, type), metadata.sections))
end

"""
Get all Field instances which are marked as representing times in the metadata.
"""
function get_time_fields(metadata::TeaFileMetadata)::Vector{Field}
    # There should never be more than one time section, and this will throw if there is not
    # a time section. Same for the item section.
    time_section = get_section(metadata, TimeSection)
    item_section = get_section(metadata, ItemSection)

    result = Field[]
    field, iter_field_state = iterate(item_section.fields)
    for time_field_offset in time_section.time_field_offsets
        while time_field_offset > field.offset
            field, iter_field_state = iterate(item_section.fields, iter_field_state)
        end
        if field.offset != time_field_offset
            throw(ArgumentError("Inconsistent metadata; could not find time field."))
        end

        push!(result, field)
    end

    return result
end

"""Get the time field which is the primary index for the tea file."""
get_primary_time_field(metadata::TeaFileMetadata) = first(get_time_fields(metadata))

"""Division that will throw if `y` does not fit into `x` an integer number of times."""
function _duration_div(x::Dates.FixedPeriod, y::T)::Int where {T <: Dates.FixedPeriod}
    # TODO use Dates.divexact?
    value, remainder = divrem(convert(T, x).value, y.value)
    if remainder != 0
        throw(ArgumentError("$y % $x != 0, division not defined."))
    end
    return value
end

"""
Get an epoch in Tea-format given it specified as a datetime.

This is the number of days since 0001-01-01 as an Int64.
"""
tea_epoch(epoch::DateTime)::Int64 = Day(epoch - DateTime(1, 1, 1)).value

"""
Given a tick duration, convert it to an (integer) number of ticks per day.
"""
function ticks_per_day(tick_duration::Dates.FixedPeriod)::Int64
    return _duration_div(Day(1), tick_duration)
end

"""
Return true iff the specfied time section implies that Int64 time fields can safely be cast
to Julia `DateTime` objects.
"""
function is_julia_time_compatible(time_section::TimeSection)::Bool
    julia_epoch_utc = DateTime(0, 1, 1) - Millisecond(Dates.DATETIMEEPOCH)
    if time_section.epoch != tea_epoch(julia_epoch_utc)
        @debug "Epoch wrong: expected $julia_epoch, got $(time_section.epoch)"
        return false
    end
    if time_section.ticks_per_day != ticks_per_day(Millisecond(1))
        return false
    end

    return true
end

"""
Return true iff the specified `Item` has a compatible memory layout with that specified
in `item_section`.

Note that we *do not* inspect field names, just types and offsets.

For fields that are of Julia `DateTime` type, we verify that this corresponds to
    * an Int64 field
    * a field marked as a "time" field
    * that the epoch & ticks_per_day match that used in the Julia `DateTime` type.
"""
function is_item_compatible(Item::Type, metadata::TeaFileMetadata)::Bool
    if !isbitstype(Item)
        # Only types with a well defined layout in memory can be compatible.
        return false
    end

    item_section::ItemSection = get_section(metadata, ItemSection)

    if sizeof(Item) != item_section.item_size
        return false
    end

    if fieldcount(Item) != length(item_section.fields)
        return false
    end

    have_done_time_checks = false
    for i in 1:fieldcount(Item)
        field = item_section.fields[i]
        if fieldtype(Item, i) == DateTime
            # DateTime fields are special, and require more extensive checks to be made
            # (see the docstring).
            if field_type(field) != Int64
                return false
            end

            if !have_done_time_checks
                # We now check the epoch and ticks-per-day for compatibility with Julia
                # date times. Since this is not field-specific, we ensure that we only do
                # this check once.
                time_section::TimeSection = get_section(metadata, TimeSection)
                if !is_julia_time_compatible(time_section)
                    return false
                end
                have_done_time_checks = true
            end

        elseif field_type(field) != fieldtype(Item, i)
            return false
        end

        if field.offset != fieldoffset(Item, i)
            return false
        end
    end

    return true
end

"""
The minimum item_start compatible with the given `metadata`.
"""
function minimum_item_start(sections::Vector{AbstractSection})::Int64
    # We need to determine item_start. We do so automatically by calculating the size of
    # the metadata on disk, and then rounding up to a multiple of 8 bytes.
    # Remeber that each section additionally has fields denoting the section type and the
    # offset of the next section.
    metadata_size = (
        sizeof(Int64) * 4  # Top-level metadata
        + length(sections) * 2 * sizeof(Int32)  # Section metadata
        + (isempty(sections) ? 0 : sum(_tea_size, sections))  # Section payloads
    )
    return round(Int, metadata_size / 8, RoundUp) * 8
end

"""Throw if `metadata` is invalid."""
function _verify_metadata(metadata::TeaFileMetadata)::Nothing
    if metadata.item_start % 8 != 0
        throw(ArgumentError(
            "Invalid item_start: $(metadata.item_start) is not a multiple of 8."
        ))
    end

    # Verify that item_start is large enough to be after the whole metadata section.
    min_item_start = minimum_item_start(metadata.sections)
    if metadata.item_start < min_item_start
        throw(ArgumentError(
            "item_start is too small: $(metadata.item_start) < $min_item_start, "
        ))
    end

    if metadata.item_end != 0 && metadata.item_end < metadata.item_start
        throw(ArgumentError(
            "Invalid item_end: $(metadata.item_end) is before item_start."
        ))
    end

    # We should have at most one section of any type.
    type_to_section = Dict{Type, AbstractSection}()
    type_to_count = DefaultDict{Type, Int}(0)
    for section in metadata.sections
        type_to_count[typeof(section)] += 1
        type_to_section[typeof(section)] = section
    end
    for (type, count) in type_to_count
        if count > 1
            throw(ArgumentError(
                "Found $count sections of type $type -- there should be at most one."
            ))
        end
    end

    # At this point we could have zero or one time section.
    if haskey(type_to_section, TimeSection)
        # If the time section is present, we will do some tests
        time_section::TimeSection = type_to_section[TimeSection]

        # We *must* have an item section if we have a time section.
        if !haskey(type_to_section, ItemSection)
            throw(ArgumentError("There is a TimeSection, but no ItemSection."))
        end
        item_section = type_to_section[ItemSection]

        unclaimed_field_offsets = Set((field.offset for field in item_section.fields))
        # Verify that time field offsets match up to known fields, and that each field is
        # referred to exactly once.
        for time_field_offset in time_section.time_field_offsets
            if !(time_field_offset in unclaimed_field_offsets)
                throw(ArgumentError(
                    "Time field offset $time_field_offset does not match to a field."
                ))
            end

            # We have claimed this field.
            delete!(unclaimed_field_offsets, time_field_offset)
        end
    end

    if haskey(type_to_section, ItemSection)
        item_section::ItemSection = type_to_section[ItemSection]

        offsets = [field.offset for field in item_section.fields]
        # `item_size` takes us to the start of the next item.
        push!(offsets, item_section.item_size)
        actual_field_sizes = diff(offsets)

        minimum_field_sizes = [sizeof(field_type(field)) for field in item_section.fields]

        for (i, (actual, minimum)) in enumerate(
                zip(actual_field_sizes, minimum_field_sizes))
            # This verifies that all item offsets are strictly increasing, and each field
            # has at least enough space on disk to store the type.
            if actual < minimum
                throw(ArgumentError(
                    "offsets imply field $i has size $actual, "
                    * "but it should be at least $minimum"
                ))
            end
        end
    end

    return nothing
end


text:
using Test
using ASE
using JuLIP.Testing

h0("Testing ASE.jl")

@testset "ASE" begin
   @testset "Atoms"  begin include("testase.jl"); end
   @testset "Calculators" begin include("test_calculators.jl"); end
   @testset "JuLIP vs ASE" begin include("test_asevsjulip.jl"); end
end


text:
@testset "LogEvent" begin

@testset "Timestamp" begin
    time_in_ms = round(Int64, time() * 1000)

    event = LogEvent("Foo", time_in_ms)
    @test CloudWatchLogs.timestamp(event) == time_in_ms

    dt = DateTime(Dates.UTM(time_in_ms + Dates.UNIXEPOCH))
    event = LogEvent("Foo", dt)
    @test CloudWatchLogs.timestamp(event) == time_in_ms

    zdt = ZonedDateTime(dt, tz"UTC")
    event = LogEvent("Foo", zdt)
    @test CloudWatchLogs.timestamp(event) == time_in_ms

    event = LogEvent("Foo")
    one_hour = Dates.value(Millisecond(Hour(1)))
    @test time_in_ms <= CloudWatchLogs.timestamp(event) <= time_in_ms + one_hour
end

@testset "Message" begin
    @test CloudWatchLogs.message(LogEvent("Foo")) == "Foo"
end

@testset "Bad construction" begin
    @test_throws ArgumentError LogEvent("")
    @test_throws ArgumentError LogEvent("Foo", -45)
    @test_warn getlogger("CloudWatchLogs") "Log Event message cannot be more than" begin
        LogEvent("A" ^ (MAX_EVENT_SIZE - 25))
    end
end

end


text:
# Author: AJ Tyler
# Date: 24/01/22
module CommonNeighbourAnalysis # This module implements the common neighbour analysis algorithm
export CNA # This function returns the total and atomic CNA profiles
export adjacencyGraph # This function returns a graphical representation of a configuration.
using Graphs # Use Graphs package to represent cluster structure and LongestPaths to calculate longest bond path

"""
This function implements the Common Neighbour Analysis (CNA) algorithm and retuns the total and atomic CNA profile.

Inputs:
configuration: (Matrix{Float64}) x,y,z positions of the cluster's atoms.
N: (Int64) Number of atoms in the cluster
rCut: (Float64) Cutoff radius in units of equilibrium bond length.
B: (Bool) Whether configurations were found under a magnetic field of 0.3 au or not.
EBL: (Float64) Equilibrium bond length (angstroms) for MP2 data for Neon2.

Outputs:
totalProfile: (Dictionary) Maps the bond type triplet identifier (i,j,k) to their frequency for the entire cluster
atomicProfile: (Dictionary) Maps the bond type triplet identifier (i,j,k) to their frequency for each atom
"""
function CNA(configuration,N,rCut,B,EBL)
	bondGraph = adjacencyGraph(configuration,N,rCut,B,EBL) # Compute graph representation of configuration.
	# Compute triplet identifiers for each pair of bonded atoms (neighbours)
	totalProfile = Dict{String,Int}() # Create Dictionary to store triplet frequencies
	atomicProfile = [Dict{String,Int}() for i in 1:N] # Create vector of dictionaries to store triplet frequencies
	for atom1 in 1:N # For all atoms
		# neighbourhood1: The subgraph which only contains the bonds between the neighbours of atom1
		# map1: map1[i] is the vertex number in bondgraph that vertex i in neighbourhood1 corresponds to
		neighbourhood1, map1 = induced_subgraph(bondGraph, neighbors(bondGraph,atom1))
		for atom2 in map1 # For each of atom1's neighbours
			if atom2 > atom1 # Check not double counting bonds
				# neighbourhood2: The subgraph which only contains the bonds between the neighbours of atom2
				# map2: map2[i] is the vertex number in bondgraph that vertex i in neighbourhood2 corresponds to
				neighbourhood2, map2 = induced_subgraph(bondGraph, neighbors(bondGraph,atom2))
				commonNeighbourhood = bondGraph[intersect(map1,map2)] # Obtain the subgraph with the common neighbours of atom1 and atom2
				i = size(commonNeighbourhood,1) # Compute the number of common neighbours
				bondsToProcess = Set(edges(commonNeighbourhood))
				j = length(bondsToProcess) # Compute the number of bonds between the common neighbours
				kMax = 0 # Initialise longest bond path length
				while (!isempty(bondsToProcess)) # While still bonds left to process
					#thisCluster = Graph(i) # Useful for debugging to see longest path construction
					atomsToProcess = Set{Int64}() # Initialise list of atoms left to process in current longest path construction
					k = 1 # Initialise current longest path length
					nextBond = pop!(bondsToProcess) # Remove a bond to process and start processing it
					#add_edge!(thisCluster,nextBond) # Useful for debugging to see longest path construction
					# Add both atoms of the bond to be processed
					push!(atomsToProcess,nextBond.src)
					push!(atomsToProcess,nextBond.dst)
					while (!isempty(atomsToProcess)) # While have atoms to process
						atom = pop!(atomsToProcess) # Remove an atom to process and start processing it
						for bond in bondsToProcess # For all the bonds left to process
							if (bond.src == atom) # If the source of the bond was the atom to process
								push!(atomsToProcess,bond.dst) # Add it's bonded atom to be processed
							elseif (bond.dst == atom) # If the source of the bond was the atom to process
								push!(atomsToProcess,bond.src) # Add it's bonded atom to be processed
							else # If atom being processed is not in the current bond
								continue # Skip to next bond
							end
							#add_edge!(thisCluster,pop!(bondsToProcess,bond)) Useful for debugging to see longest path construction
							pop!(bondsToProcess) # Finished processing the bond
							k += 1 # Increment the longest bond path length
						end
					end
					kMax = max(k,kMax)
				end

				key = "($i,$j,$(Int(kMax)))" # Create triplet identifier key for dictionary
				push!(totalProfile, key => get!(totalProfile,key,0)+1) # Increments total CNA triplet frequency by 1
				push!(atomicProfile[atom1], key => get!(atomicProfile[atom1],key,0)+1) # Increments atom1 triplet frequency by 1
				push!(atomicProfile[atom2], key => get!(atomicProfile[atom2],key,0)+1) # Increments atom2 triplet frequency by 1
			end
		end
	end
	return totalProfile, atomicProfile # Return the total and atomic CNA profiles
end

"""
This function constructs an undirected graph representation of a configuration. The nodes are the atoms which have an edge between them if
they are bonded. Atoms are bonded if the distance between them is less than rCut*EBL.

Inputs:
configuration: (Matrix{Float64}) x,y,z positions of the cluster's atoms.
N: (Int64) Number of atoms in the cluster
rCut: (Float64) Cutoff radius in units of equilibrium bond length.
B: (Bool) Whether configurations were found under a magnetic field of 0.3 au or not.
EBL: (Float64) Equilibrium bond length (angstroms) for MP2 data for Neon2.

Output:
bondGraph: (SimpleGraph) Graph representation of the configuration.
"""
function adjacencyGraph(configuration,N,rCut,B,EBL)
	# Create graph representation of bonded atoms, where the atoms are vertices are the vertices and the edges are bonds
	bondGraph = SimpleGraph(N) # Constructs a SimpleGraph object with N atoms(vertices)
	# For all pairs of atoms
	for j in 1:N-1
		for i in j+1:N
			diff = configuration[j,:]-configuration[i,:] # calculate vector between atoms
			norm = sqrt(sum(diff.*diff)) # Compute length of bond
			if B # If strong magnetic field present
				theta = acos(diff[3]/norm) # Compute angle bond makes with B field (z-axis) in radians
				rCutScaled = EBL(theta)*rCut # Compute rCut is terms of equilibrium bond length in the theta direction
			else # If no magnetic field
				rCutScaled = EBL*rCut # Use standard equilibrium Bond Length.
			end
			 # If the distance between atoms is less the cut-off radius
			if (norm < rCutScaled)
				add_edge!(bondGraph,i,j) # Add bond (edge) to graph
			end
		end
	end

	return bondGraph # Return graph
end

"""
This function computes the equilibrium bond length in angstroms. If a 0.3 au magentic field is present, then bond angle with the B field
direction is used in a harmonic function to output the altered equlibrium bond length due to magnetic effects.

Input:
theta: (Float64) Angle the bond makes with the magnetic field in radians.

Output:
EBL: (Float64) The equilibrium bond length at angle theta.
"""
function EBL(theta)
	if theta > pi/2 # If theta is greater than pi/2 radians
		theta = pi-theta # Change theta to be within [0,pi/2] radians for interpolation purposes
	end
	return 0.188*cos(2*theta)+0.01073*cos(4*theta)-0.00269*cos(8*theta)+3.06621 # Return prefitted harmonic function
end

end # End of module


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

testtables = [
  (
    "basic",
    (col1=Int64[1,2,3,4,5,6,7,8,9,10],),
    NamedTuple(),
    NamedTuple(),
    nothing
  ),
  (
    "missing values",
    (col1=Union{Int64, Missing}[1,2,3,4,5,6,7,8,9,missing],),
    NamedTuple(),
    NamedTuple(),
    nothing
  ),
  (
    "primitive types",
    (
      col1=[missing, missing, missing, missing],
      col2=Union{UInt8, Missing}[0, 1, 2, missing],
      col3=Union{UInt16, Missing}[0, 1, 2, missing],
      col4=Union{UInt32, Missing}[0, 1, 2, missing],
      col5=Union{UInt64, Missing}[0, 1, 2, missing],
      col6=Union{Int8, Missing}[0, 1, 2, missing],
      col7=Union{Int16, Missing}[0, 1, 2, missing],
      col8=Union{Int32, Missing}[0, 1, 2, missing],
      col9=Union{Int64, Missing}[0, 1, 2, missing],
      col10=Union{Float16, Missing}[0, 1, 2, missing],
      col11=Union{Float32, Missing}[0, 1, 2, missing],
      col12=Union{Float64, Missing}[0, 1, 2, missing],
      col13=[true, false, true, missing],
    ),
    NamedTuple(),
    NamedTuple(),
    nothing
  ),
  (
    "arrow date/time types",
    (
      col14=[zero(Arrow.Decimal{Int32(2), Int32(2), Int128}), zero(Arrow.Decimal{Int32(2), Int32(2), Int128}), zero(Arrow.Decimal{Int32(2), Int32(2), Int128}), missing],
      col15=[zero(Arrow.Date{Arrow.Meta.DateUnit.DAY, Int32}), zero(Arrow.Date{Arrow.Meta.DateUnit.DAY, Int32}), zero(Arrow.Date{Arrow.Meta.DateUnit.DAY, Int32}), missing],
      col16=[zero(Arrow.Time{Arrow.Meta.TimeUnit.SECOND, Int32}), zero(Arrow.Time{Arrow.Meta.TimeUnit.SECOND, Int32}), zero(Arrow.Time{Arrow.Meta.TimeUnit.SECOND, Int32}), missing],
      col17=[zero(Arrow.Timestamp{Arrow.Meta.TimeUnit.SECOND, nothing}), zero(Arrow.Timestamp{Arrow.Meta.TimeUnit.SECOND, nothing}), zero(Arrow.Timestamp{Arrow.Meta.TimeUnit.SECOND, nothing}), missing],
      col18=[zero(Arrow.Interval{Arrow.Meta.IntervalUnit.YEAR_MONTH, Int32}), zero(Arrow.Interval{Arrow.Meta.IntervalUnit.YEAR_MONTH, Int32}), zero(Arrow.Interval{Arrow.Meta.IntervalUnit.YEAR_MONTH, Int32}), missing],
      col19=[zero(Arrow.Duration{Arrow.Meta.TimeUnit.SECOND}), zero(Arrow.Duration{Arrow.Meta.TimeUnit.SECOND}), zero(Arrow.Duration{Arrow.Meta.TimeUnit.SECOND}), missing],
    ),
    NamedTuple(),
    (convert=false,),
    nothing
  ),
  (
    "list types",
    (
      col1=Union{String, Missing}["hey", "there", "sailor", missing],
      col2=Union{Vector{UInt8}, Missing}[b"hey", b"there", b"sailor", missing],
      col3=Union{Vector{Int64}, Missing}[Int64[1], Int64[2], Int64[3], missing],
      col4=Union{NTuple{2, Vector{Int64}},Missing}[(Int64[1], Int64[2]), missing, missing, (Int64[3], Int64[4])],
      col5=Union{NTuple{2, UInt8}, Missing}[(0x01, 0x02), (0x03, 0x04), missing, (0x05, 0x06)],
      col6=NamedTuple{(:a, :b), Tuple{Int64, String}}[(a=Int64(1), b="hey"), (a=Int64(2), b="there"), (a=Int64(3), b="sailor"), (a=Int64(4), b="jo-bob")],
    ),
    NamedTuple(),
    NamedTuple(),
    nothing
  ),
  (
    "unions",
    (
      col1=Arrow.DenseUnionVector( Union{Int64, Float64, Missing}[1, 2.0, 3, 4.0, missing]),
      col2=Arrow.SparseUnionVector(Union{Int64, Float64, Missing}[1, 2.0, 3, 4.0, missing]),
    ),
    NamedTuple(),
    NamedTuple(),
    nothing
  ),
  (
    "dict encodings",
    (
      col1=Arrow.DictEncode(Int64[4, 5, 6]),
    ),
    NamedTuple(),
    NamedTuple(),
    function (tt)
      col1 = copy(tt.col1)
      @test typeof(col1) == PooledVector{Int64, Int8, Vector{Int8}}
    end
  ),
  (
    "more dict encodings",
    (
      col1=Arrow.DictEncode(NamedTuple{(:a, :b), Tuple{Int64, Union{String, Missing}}}[(a=Int64(1), b=missing), (a=Int64(1), b=missing), (a=Int64(3), b="sailor"), (a=Int64(4), b="jo-bob")]),
    ),
    NamedTuple(),
    NamedTuple(),
    nothing
  ),
  (
    "PooledArray",
    (
      col1=PooledArray([4,5,6,6]),
    ),
    NamedTuple(),
    NamedTuple(),
    nothing
  ),
  (
    "auto-converting types",
    (
      col1=[Date(2001, 1, 2), Date(2010, 10, 10), Date(2020, 12, 1)],
      col2=[Time(1, 1, 2), Time(13, 10, 10), Time(22, 12, 1)],
      col3=[DateTime(2001, 1, 2), DateTime(2010, 10, 10), DateTime(2020, 12, 1)],
      col4=[ZonedDateTime(2001, 1, 2, TimeZone("America/Denver")), ZonedDateTime(2010, 10, 10, TimeZone("America/Denver")), ZonedDateTime(2020, 12, 1, TimeZone("America/Denver"))]
    ),
    NamedTuple(),
    NamedTuple(),
    nothing
  ),
  (
    "Map",
    (
      col1=[Dict(Int32(1) => Float32(3.14)), missing],
    ),
    NamedTuple(),
    NamedTuple(),
    nothing
  ),
  (
    "non-standard types",
    (
      col1=[:hey, :there, :sailor],
      col2=['a', 'b', 'c'],
      col3=Arrow.DictEncode(['a', 'a', 'b'])
    ),
    NamedTuple(),
    NamedTuple(),
    nothing
  ),
  (
    "large lists",
    (
      col1=Union{String, Missing}["hey", "there", "sailor", missing],
      col2=Union{Vector{UInt8}, Missing}[b"hey", b"there", b"sailor", missing],
      col3=Union{Vector{Int64}, Missing}[Int64[1], Int64[2], Int64[3], missing],
      col4=Union{NTuple{2, Vector{Int64}},Missing}[(Int64[1], Int64[2]), missing, missing, (Int64[3], Int64[4])],
      col5=Union{NTuple{2, UInt8}, Missing}[(0x01, 0x02), (0x03, 0x04), missing, (0x05, 0x06)],
      col6=NamedTuple{(:a, :b), Tuple{Int64, String}}[(a=Int64(1), b="hey"), (a=Int64(2), b="there"), (a=Int64(3), b="sailor"), (a=Int64(4), b="jo-bob")],
    ),
    (largelists=true,),
    NamedTuple(),
    nothing
  ),
  (
    "dictencode keyword",
    (
      col1=Int64[1,2,3,4],
      col2=Union{String, Missing}["hey", "there", "sailor", missing],
      col3=Arrow.DictEncode(NamedTuple{(:a, :b), Tuple{Int64, Union{String, Missing}}}[(a=Int64(1), b=missing), (a=Int64(1), b=missing), (a=Int64(3), b="sailor"), (a=Int64(4), b="jo-bob")]),
      col4=[:a, :b, :c, missing],
      col5=[Date(2020, 1, 1) for x = 1:4]
    ),
    (dictencode=true,),
    NamedTuple(),
    nothing
  ),
  (
    "nesteddictencode keyword",
    (
      col1=NamedTuple{(:a, :b), Tuple{Int64, Union{Missing, NamedTuple{(:c,), Tuple{String}}}}}[(a=Int64(1), b=missing), (a=Int64(1), b=missing), (a=Int64(3), b=(c="sailor",)), (a=Int64(4), b=(c="jo-bob",))],
    ),
    (dictencode=true, dictencodenested=true,),
    NamedTuple(),
    nothing
  ),
  (
    "Julia unions",
    (
      col1=Union{Int, String}[1, "hey", 2, "ho"],
      col2=Union{Char, NamedTuple{(:a,), Tuple{Symbol}}}['a', (a=:hey,), 'b', (a=:ho,)],
    ),
    (denseunions=false,),
    NamedTuple(),
    nothing
  ),
  (
    "Decimal256",
    (
      col1=[zero(Arrow.Decimal{Int32(2), Int32(2), Arrow.Int256}), zero(Arrow.Decimal{Int32(2), Int32(2), Arrow.Int256}), zero(Arrow.Decimal{Int32(2), Int32(2), Arrow.Int256}), missing],
    ),
    NamedTuple(),
    (convert=false,),
    nothing
  ),
];

function testtable(nm, t, writekw, readkw, extratests)
  println("testing: $nm")
  io = IOBuffer()
  Arrow.write(io, t; writekw...)
  seekstart(io)
  tt = Arrow.Table(io; readkw...)
  @test length(tt) == length(t)
  @test all(isequal.(values(t), values(tt)))
  extratests !== nothing && extratests(tt)
  seekstart(io)
  str = Arrow.Stream(io; readkw...)
  tt = first(str)
  @test length(tt) == length(t)
  @test all(isequal.(values(t), values(tt)))
  # compressed
  io = IOBuffer()
  Arrow.write(io, t; compress=((:lz4, :zstd)[rand(1:2)]), writekw...)
  seekstart(io)
  tt = Arrow.Table(io; readkw...)
  @test length(tt) == length(t)
  @test all(isequal.(values(t), values(tt)))
  extratests !== nothing && extratests(tt)
  seekstart(io)
  str = Arrow.Stream(io; readkw...)
  tt = first(str)
  @test length(tt) == length(t)
  @test all(isequal.(values(t), values(tt)))
  # file
  io = IOBuffer()
  Arrow.write(io, t; file=true, writekw...)
  seekstart(io)
  tt = Arrow.Table(io; readkw...)
  @test length(tt) == length(t)
  @test all(isequal.(values(t), values(tt)))
  extratests !== nothing && extratests(tt)
  seekstart(io)
  str = Arrow.Stream(io; readkw...)
  tt = first(str)
  @test length(tt) == length(t)
  @test all(isequal.(values(t), values(tt)))
  return
end


text:
function getoffset(a::TopLevel{D}) where {D}
    return oneunit(Address{D}) - Address(dim_min(addresses(a)))
end

"""
    mappables(a::TopLevel{D}, address::Address{D})

Return a `Vector{Path{Component}}` of paths to mappable components at `address`.
"""
function mappables(toplevel::TopLevel, ruleset::RuleSet, address::Address)
    return [
        path for
        path in walk_children(toplevel, address) if ismappable(ruleset, toplevel[path])
    ]
end

"""
    check(c::AbstractComponent)

Check a component for ports that are not connected to any link. Only applies
to ports visible to the level of hierarchy of `c`. That is, children of `c` will
not be checked.

Returns a `Vector{PortPath}` of unused ports.
"""
function check(component::AbstractComponent)
    # Gather all visible ports of `c` and all connected visible ports of `c`.
    # The unused ports are just the set diference.
    all_ports = visible_ports(component)
    used_ports = connected_ports(component)

    unused_ports = setdiff(all_ports, used_ports)

    # Debug printing
    if length(unused_ports) > 0
        @debug begin
            n_unused = length(unused_ports)
            str = "Component $(component.name) has $(n_unused) unused ports. $unused_ports"
        end
    end
    return unused_ports
end

################################################################################
# visible_ports

function visible_ports(component::AbstractComponent)
    # Get ports defined at this level of hierarchy
    paths = portpaths(component)
    # Collect all external ports of children.
    for (name, child) in component.children
        for port in portpaths(child)
            path = catpath(name, port)
            push!(paths, path)
        end
    end
    return paths
end

@doc """
    visible_ports(component::AbstractComponent)

Return `Vector{PortPath}` of the ports of `component` and the ports of the
children of `component`.
""" visible_ports

################################################################################
# Connectedness methods.
################################################################################

function connectedlink(toplevel::TopLevel, portpath::Path{Port}, portclass::PortClass)
    # Extract the port type of the top level definition.
    port = toplevel[portpath]

    # If the port direction differs from "class", then we must go up one level
    # in the component hierarchy to get the attached link.
    up_one = port.class == portclass

    split_amount = up_one ? 2 : 1
    cpath, ppath = splitpath(portpath, split_amount)

    # Get the component type and look up the port name in the "portlink"
    # dictionary to get its connected link.
    component = toplevel[cpath]
    link = component.portlink[ppath]

    # Prefix the name of the component to the link's name to get the full path
    # relative to the toplevel.
    return catpath(cpath, Path{Link}(link.name))
end

function connectedports(toplevel::TopLevel, linkpath::Path{Link}, class::PortClass)
    link = toplevel[linkpath]
    # Check directionality
    if class == Input
        ports = sources(link)
    elseif class == Output
        ports = dests(link)
    else
        throw(DomainError())
    end

    # Strip the name of the link to just get the prefix.
    linkprefix = striplast(linkpath)
    return catpath.(Ref(linkprefix), ports)

    # If the link is global, its sources and destinations are already full
    # address paths, just return these arrays directly.
    #
    # Otherwise, need to prefix the paths local to the link with the path to
    # the component where the links are defined.

    # isgloballink(linkpath) && (return ports)
    # linkprefix = striplast(linkpath)
    # return [catpath(linkprefix, i) for i in ports]
end

################################################################################
# isconnected
################################################################################

# Default to false
function isconnected(t::TopLevel, a::Path, b::Path)
    @error "Undefined connection"
    return false
end

function isconnected(toplevel::TopLevel, portpath::Path{Port}, linkpath::Path{Link})

    # Get the output link connected to the port
    actuallink = connectedlink(toplevel, portpath, Output)
    if linkpath != actuallink
        @error """
            $portpath is not connected to $linkpath.

            It is connected to $actuallink.
            """
        return false
    end
    actualports = connectedports(toplevel, linkpath, Input)
    if portpath ∉ actualports
        @error """
            $linkpath is not connected to $portpath.

            Actual sink ports: $actualports.
            """
    end
    return true
end

function isconnected(toplevel::TopLevel, linkpath::Path{Link}, portpath::Path{Port})
    # Get the output link connected to the port
    actuallink = connectedlink(toplevel, portpath, Input)
    if linkpath != actuallink
        @error """
            $portpath is not connected to $linkpath.

            It is connected to $actuallink.
            """
        return false
    end
    actualports = connectedports(toplevel, linkpath, Output)
    if portpath ∉ actualports
        @error """
            $linkpath is not connected to $portpath.

            Actual sink ports: $actualports.
            """
    end
    return true
end

function isconnected(
    toplevel::TopLevel,
    portpath::Path{Port},
    #cpath::Path{Component}
    componentpath::Path{Component},
)

    # Check if the port is an input of the component by collecting all of the
    # inputs ports of the component and checking for path equality.
    component = toplevel[componentpath]
    component_inputs = catpath.(Ref(componentpath), portpaths(component, (Input,)))

    if portpath ∉ component_inputs
        @error "$portpath is not an input of $componentpath."
        return false
    end
    return true
end

function isconnected(
    toplevel::TopLevel, componentpath::Path{Component}, portpath::Path{Port}
)

    # Do the same trick as above - get all the output ports of the component
    # by path and check if the given portpath is in them.
    component = toplevel[componentpath]
    component_outputs = catpath.(Ref(componentpath), portpaths(component, (Output,)))

    if portpath ∉ component_outputs
        @error "$portpath is not an input of $componentpath."
        return false
    end
    return true
end

@doc """
    isconnected(toplevel, a::AbstractPath, b::AbstractPath)

Return `true` if architectural component referenced by path `a` is
architecturally connected to that referenced by path `b`.

The order of the arguments is important for directed components. For example,
if `a` references a port that is a source for link `b` in `toplevel`, then

```julia
julia> isconnected(toplevel, a, b)
true

julia> isconnected(toplevel, b, a)
false
```

If one of `a` or `b` is of type `ComponentPath`, then only ports are considered
connected.
""" isconnected


text:
# fbp_ramp.jl

export fbp_ramp
export ramp_flat
export ramp_arc

using MIRT: SinoFanFlat, SinoFanArc

"""
    h, nn = fbp_ramp(how, n, ds, dsd)
'ramp-like' filters for parallel-beam and fan-beam FBP reconstruction.
This sampled band-limited approach avoids the aliasing that would be
caused by sampling the ramp directly in the frequency domain.

in
- `how::Symbol`
  `:arc` (3rd generation CT) or `:flat` (for parallel too)
- `n::Int` : # of samples (must be even)
- `ds::RealU` : sample spacing (in distance units, e.g., cm)
- `dsd::RealU` : source-to-detector distance, for `:arc` case

out
- `h::AbstractMatrix{<:Number}` : samples of band-limited ramp filter

"""
fbp_ramp(::SinoPar, n::Int, ds::RealU, dsd::RealU) = ramp_flat(n, ds)
fbp_ramp(::SinoFanFlat, n::Int, ds::RealU, dsd::RealU) = ramp_flat(n, ds)
fbp_ramp(::SinoFanArc, n::Int, ds::RealU, dsd::RealU) = ramp_arc(n, ds, dsd)

function ramp_arc(n::Int, ds::RealU, dsd::RealU)
    isodd(n) && throw("n must be even")

    if n/2 * ds / dsd > 0.9 * π/2
        println("angle is $(rad2deg(n/2 * ds / dsd)) degrees: too large")
        throw("physically impossible arc geometry")
    end

    nn = -(n÷2):(n÷2-1)
    h = zeros(size(nn))
    h[nn .== 0] .= 1 / (4 * abs2(ds))
    odd = isodd.(nn)
    h[odd] = @.( -1 / abs2((pi * dsd * sin(nn[odd] * ds / dsd))))

    return h, nn
end

function ramp_flat(n::Int, ds::RealU)

    isodd(n) && throw("n must be even")

    nn = -(n÷2):(n÷2-1)
    h = zeros(size(nn))
    h[n÷2+1] = 1 / 4
    odd = isodd.(nn)
    h[odd] .= -1 ./ abs2.(π * nn[odd])
    h ./= abs2(ds)

    return h, nn
end


text:
"""
    `status_stopping_to_stats(stp::AbstractStopping)`

Return the status in GenericExecutionStats from a Stopping.
"""
function status_stopping_to_stats(stp::AbstractStopping)
  stp_status = status(stp)
  convert = Dict([
    (:Optimal, :first_order),
    (:SubProblemFailure, :unknown),
    (:SubOptimal, :acceptable),
    (:Unbounded, :unbounded),
    (:UnboundedPb, :unbounded),
    (:Stalled, :stalled),
    (:IterationLimit, :max_iter),
    (:TimeLimit, :max_time),
    (:EvaluationLimit, :max_eval),
    (:ResourcesOfMainProblemExhausted, :max_eval),
    (:Infeasible, :infeasible),
    (:DomainError, :exception),
    (:StopByUser, :user),
    (:Exception, :exception),
    (:Unknown, :unknown),
  ])
  return convert[stp_status]
end

"""
    `stopping_to_stats(stp::NLPStopping)`

Initialize a GenericStats from Stopping
"""
function stopping_to_stats(stp::NLPStopping)
  nlp = stp.pb
  cx = stp.current_state.cx
  return GenericExecutionStats(
    status_stopping_to_stats(stp),
    stp.pb,
    solution = stp.current_state.x,
    objective = stp.current_state.fx,
    primal_feas = get_ncon(nlp) == 0 ? zero(eltype(cx)) :
                  maximum(max.(cx - get_ucon(nlp), get_lcon(nlp) - cx, 0)),
    dual_feas = norm(stp.current_state.current_score, Inf), # stp.current_state.res
    multipliers = stp.current_state.lambda,
    iter = stp.meta.nb_of_stop,
    elapsed_time = stp.current_state.current_time - stp.meta.start_time,
  )
end

function stats_status_to_meta!(stp::AbstractStopping, stats::GenericExecutionStats)
  return stats_status_to_meta!(stp, stats.status)
end

function stats_status_to_meta!(stp::AbstractStopping, status::Symbol)
  #Update the meta boolean with the output message
  if status == :first_order
    stp.meta.optimal = true
  end
  if status == :acceptable
    stp.meta.suboptimal = true
  end
  if status == :infeasible
    stp.meta.infeasible = true
  end
  if status == :small_step
    stp.meta.stalled = true
  end
  if status == :max_eval
    stp.meta.resources = true
  end
  if status == :max_iter
    stp.meta.iteration_limit = true
  end
  if status == :max_time
    stp.meta.tired = true
  end
  if status ∈ [:neg_pred, :not_desc]
    stp.meta.fail_sub_pb = true
  end
  if status == :unbounded
    stp.meta.unbounded = true
  end
  if status == :user
    stp.meta.stopbyuser = true
  end
  if status ∈ [:stalled, :small_residual, :small_step]
    stp.meta.stalled = true
  end
  if status == :exception
    stp.meta.exception = true
  end #available ≥ 0.2.6

  return stp
end


text:
using BinomialSynapses
using CUDA
using Distributions: Exponential

using BenchmarkTools
using Test

const RUN_BENCHMARKS = false # optional intermediate benchmarks

@testset "BinomialSynapses.jl" begin
    include("models.jl")
    include("propagate.jl")
    include("timestep.jl")
    include("emission.jl")
    include("likelihood.jl")
    include("jitter.jl")
    include("resample.jl")
    include("filter.jl")
    include("statistics.jl")
    include("simulate.jl")

    include("_integration.jl")
end


text:
export
    FoldAssigner,
    RandomFoldAssigner,
    OrderedFoldAssigner,

    assign_folds

abstract type FoldAssigner end
assign_folds(fa::FoldAssigner, nitems::Int, nfolds::Int) = error("get_fold_assignment not implemented for $a")

"""
    RandomFoldAssigner
Assign folds using a random permutation
"""
type RandomFoldAssigner <: FoldAssigner
    use_all_samples::Bool # true - use all samples and assign extra ones to the first folds
                          # false - ensure all folds are of the same size
    RandomFoldAssigner(use_all_samples::Bool=true) = new(use_all_samples)
end
function assign_folds(fa::RandomFoldAssigner, nitems::Int, nfolds::Int)
    assignment = zeros(Int, nitems)

    N = fa.use_all_samples ? nitems : nfolds*div(nitems, nfolds)

    p = randperm(nitems)
    fold = 1
    for i in 1 : N
        assignment[p[i]] = fold
        fold += 1
        if fold > nfolds
            fold = 1
        end
    end

    assignment
end

"""
    OrderedFoldAssigner
Assigns the first values to fold 1, the next block to fold 2, etc.
"""
type OrderedFoldAssigner <: FoldAssigner
    use_all_samples::Bool # true - use all samples and assign extra ones to the first folds
                          # false - ensure all folds are of the same size
    OrderedFoldAssigner(use_all_samples::Bool=true) = new(use_all_samples)
end
function assign_folds(fa::OrderedFoldAssigner, nitems::Int, nfolds::Int)
    assignment = zeros(Int, nitems)

    fold_count = div(nitems, nfolds)
    n_folds_with_extra = rem(nitems, nfolds)

    i = 0
    for fold in 1 : nfolds
        this_fold_count = fold_count
        if fa.use_all_samples && fold ≤ n_folds_with_extra
            this_fold_count += 1
        end

        for j in 1 : this_fold_count
            assignment[i+=1] = fold
        end
    end

    assignment
end