text:
# Autogenerated wrapper script for ViennaRNA_jll for i686-linux-gnu-cxx11
export AnalyseDists, AnalyseSeqs, Kinfold, Kinfold, RNA2Dfold, RNALalifold, RNALfold, RNAPKplex, RNAaliduplex, RNAalifold, RNAcofold, RNAdistance, RNAdos, RNAduplex, RNAeval, RNAfold, RNAforester, RNAheat, RNAinverse, RNAlocmin, RNAmultifold, RNApaln, RNAparconv, RNApdist, RNAplex, RNAplfold, RNAplot, RNApvmin, RNAsnoop, RNAsubopt, RNAup, b2ct, ct2db, libRNA, popt

using MPFR_jll
using CompilerSupportLibraries_jll
using GSL_jll
JLLWrappers.@generate_wrapper_header("ViennaRNA")
JLLWrappers.@declare_executable_product(AnalyseDists)
JLLWrappers.@declare_executable_product(AnalyseSeqs)
JLLWrappers.@declare_executable_product(Kinfold)
JLLWrappers.@declare_executable_product(Kinfold)
JLLWrappers.@declare_executable_product(RNA2Dfold)
JLLWrappers.@declare_executable_product(RNALalifold)
JLLWrappers.@declare_executable_product(RNALfold)
JLLWrappers.@declare_executable_product(RNAPKplex)
JLLWrappers.@declare_executable_product(RNAaliduplex)
JLLWrappers.@declare_executable_product(RNAalifold)
JLLWrappers.@declare_executable_product(RNAcofold)
JLLWrappers.@declare_executable_product(RNAdistance)
JLLWrappers.@declare_executable_product(RNAdos)
JLLWrappers.@declare_executable_product(RNAduplex)
JLLWrappers.@declare_executable_product(RNAeval)
JLLWrappers.@declare_executable_product(RNAfold)
JLLWrappers.@declare_executable_product(RNAforester)
JLLWrappers.@declare_executable_product(RNAheat)
JLLWrappers.@declare_executable_product(RNAinverse)
JLLWrappers.@declare_executable_product(RNAlocmin)
JLLWrappers.@declare_executable_product(RNAmultifold)
JLLWrappers.@declare_executable_product(RNApaln)
JLLWrappers.@declare_executable_product(RNAparconv)
JLLWrappers.@declare_executable_product(RNApdist)
JLLWrappers.@declare_executable_product(RNAplex)
JLLWrappers.@declare_executable_product(RNAplfold)
JLLWrappers.@declare_executable_product(RNAplot)
JLLWrappers.@declare_executable_product(RNApvmin)
JLLWrappers.@declare_executable_product(RNAsnoop)
JLLWrappers.@declare_executable_product(RNAsubopt)
JLLWrappers.@declare_executable_product(RNAup)
JLLWrappers.@declare_executable_product(b2ct)
JLLWrappers.@declare_executable_product(ct2db)
JLLWrappers.@declare_library_product(libRNA, "libRNA.so")
JLLWrappers.@declare_executable_product(popt)
function __init__()
    JLLWrappers.@generate_init_header(MPFR_jll, CompilerSupportLibraries_jll, GSL_jll)
    JLLWrappers.@init_executable_product(
        AnalyseDists,
        "bin/AnalyseDists",
    )

    JLLWrappers.@init_executable_product(
        AnalyseSeqs,
        "bin/AnalyseSeqs",
    )

    JLLWrappers.@init_executable_product(
        Kinfold,
        "bin/kinwalker",
    )

    JLLWrappers.@init_executable_product(
        Kinfold,
        "bin/Kinfold",
    )

    JLLWrappers.@init_executable_product(
        RNA2Dfold,
        "bin/RNA2Dfold",
    )

    JLLWrappers.@init_executable_product(
        RNALalifold,
        "bin/RNALalifold",
    )

    JLLWrappers.@init_executable_product(
        RNALfold,
        "bin/RNALfold",
    )

    JLLWrappers.@init_executable_product(
        RNAPKplex,
        "bin/RNAPKplex",
    )

    JLLWrappers.@init_executable_product(
        RNAaliduplex,
        "bin/RNAaliduplex",
    )

    JLLWrappers.@init_executable_product(
        RNAalifold,
        "bin/RNAalifold",
    )

    JLLWrappers.@init_executable_product(
        RNAcofold,
        "bin/RNAcofold",
    )

    JLLWrappers.@init_executable_product(
        RNAdistance,
        "bin/RNAdistance",
    )

    JLLWrappers.@init_executable_product(
        RNAdos,
        "bin/RNAdos",
    )

    JLLWrappers.@init_executable_product(
        RNAduplex,
        "bin/RNAduplex",
    )

    JLLWrappers.@init_executable_product(
        RNAeval,
        "bin/RNAeval",
    )

    JLLWrappers.@init_executable_product(
        RNAfold,
        "bin/RNAfold",
    )

    JLLWrappers.@init_executable_product(
        RNAforester,
        "bin/RNAforester",
    )

    JLLWrappers.@init_executable_product(
        RNAheat,
        "bin/RNAheat",
    )

    JLLWrappers.@init_executable_product(
        RNAinverse,
        "bin/RNAinverse",
    )

    JLLWrappers.@init_executable_product(
        RNAlocmin,
        "bin/RNAlocmin",
    )

    JLLWrappers.@init_executable_product(
        RNAmultifold,
        "bin/RNAmultifold",
    )

    JLLWrappers.@init_executable_product(
        RNApaln,
        "bin/RNApaln",
    )

    JLLWrappers.@init_executable_product(
        RNAparconv,
        "bin/RNAparconv",
    )

    JLLWrappers.@init_executable_product(
        RNApdist,
        "bin/RNApdist",
    )

    JLLWrappers.@init_executable_product(
        RNAplex,
        "bin/RNAplex",
    )

    JLLWrappers.@init_executable_product(
        RNAplfold,
        "bin/RNAplfold",
    )

    JLLWrappers.@init_executable_product(
        RNAplot,
        "bin/RNAplot",
    )

    JLLWrappers.@init_executable_product(
        RNApvmin,
        "bin/RNApvmin",
    )

    JLLWrappers.@init_executable_product(
        RNAsnoop,
        "bin/RNAsnoop",
    )

    JLLWrappers.@init_executable_product(
        RNAsubopt,
        "bin/RNAsubopt",
    )

    JLLWrappers.@init_executable_product(
        RNAup,
        "bin/RNAup",
    )

    JLLWrappers.@init_executable_product(
        b2ct,
        "bin/b2ct",
    )

    JLLWrappers.@init_executable_product(
        ct2db,
        "bin/ct2db",
    )

    JLLWrappers.@init_library_product(
        libRNA,
        "lib/libRNA.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        popt,
        "bin/popt",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function __init__()
    initialize_prefs()
end

ext: jl


text:
abstract type AbstractFunJac{J2} end
mutable struct FunJac{F, F2, J, P, M, J2, uType, uType2, Prec, PS} <: AbstractFunJac{J2}
    fun::F
    fun2::F2
    jac::J
    p::P
    mass_matrix::M
    jac_prototype::J2
    prec::Prec
    psetup::PS
    u::uType
    du::uType
    resid::uType2
end
FunJac(fun,jac,p,m,jac_prototype,prec,psetup,u,du) = FunJac(fun,nothing,jac,p,m,jac_prototype,prec,psetup,u,du,nothing)
FunJac(fun,jac,p,m,jac_prototype,prec,psetup,u,du,resid) = FunJac(fun,nothing,jac,p,m,jac_prototype,prec,psetup,u,du,resid)

function cvodefunjac(t::Float64,
                     u::N_Vector,
                     du::N_Vector,
                     funjac::FunJac)
    funjac.u = unsafe_wrap(Vector{Float64}, __N_VGetArrayPointer_Serial(u),length(funjac.u))
    funjac.du = unsafe_wrap(Vector{Float64}, __N_VGetArrayPointer_Serial(du),length(funjac.du))
    _du = funjac.du
    _u = funjac.u
    funjac.fun(_du, _u, funjac.p, t)
    return CV_SUCCESS
end

function cvodefunjac2(t::Float64,
                     u::N_Vector,
                     du::N_Vector,
                     funjac::FunJac)
    funjac.u = unsafe_wrap(Vector{Float64}, __N_VGetArrayPointer_Serial(u),length(funjac.u))
    funjac.du = unsafe_wrap(Vector{Float64}, __N_VGetArrayPointer_Serial(du),length(funjac.du))
    _du = funjac.du
    _u = funjac.u
    funjac.fun2(_du, _u, funjac.p, t)
    return CV_SUCCESS
end

function cvodejac(t::realtype,
                  u::N_Vector,
                  du::N_Vector,
                  J::SUNMatrix,
                  funjac::AbstractFunJac{Nothing},
                  tmp1::N_Vector,
                  tmp2::N_Vector,
                  tmp3::N_Vector)

    funjac.u = unsafe_wrap(Vector{Float64}, __N_VGetArrayPointer_Serial(u),length(funjac.u))
    _u = funjac.u
    funjac.jac(convert(Matrix, J), _u, funjac.p, t)
    return CV_SUCCESS
end

function cvodejac(t::realtype,
                  u::N_Vector,
                  du::N_Vector,
                  _J::SUNMatrix,
                  funjac::AbstractFunJac{<:SparseMatrixCSC},
                  tmp1::N_Vector,
                  tmp2::N_Vector,
                  tmp3::N_Vector)
    jac_prototype = funjac.jac_prototype
    J = convert(SparseMatrixCSC,_J)

    funjac.u = unsafe_wrap(Vector{Float64}, __N_VGetArrayPointer_Serial(u),length(funjac.u))
    _u = funjac.u

    funjac.jac(jac_prototype, _u, funjac.p, t)
    J.nzval .= jac_prototype.nzval
    # Sundials resets the value pointers each time, so reset it too
    @. J.rowval = jac_prototype.rowval - 1
    @. J.colptr = jac_prototype.colptr - 1
    return CV_SUCCESS
end

function idasolfun(t::Float64, u::N_Vector, du::N_Vector, resid::N_Vector, funjac::FunJac)
    funjac.u = unsafe_wrap(Vector{Float64}, __N_VGetArrayPointer_Serial(u),length(funjac.u))
    _u = funjac.u
    funjac.du = unsafe_wrap(Vector{Float64}, __N_VGetArrayPointer_Serial(du),length(funjac.du))
    _du = funjac.du
    funjac.resid = unsafe_wrap(Vector{Float64}, __N_VGetArrayPointer_Serial(resid),length(funjac.resid))
    _resid = funjac.resid
    funjac.fun(_resid, _du, _u, funjac.p, t)
    return IDA_SUCCESS
end

function idajac(t::realtype,
                cj::realtype,
                u::N_Vector,
                du::N_Vector,
                res::N_Vector,
                J::SUNMatrix,
                funjac::AbstractFunJac{Nothing},
                tmp1::N_Vector,
                tmp2::N_Vector,
                tmp3::N_Vector)

    funjac.u = unsafe_wrap(Vector{Float64}, __N_VGetArrayPointer_Serial(u),length(funjac.u))
    _u = funjac.u
    funjac.du = unsafe_wrap(Vector{Float64}, __N_VGetArrayPointer_Serial(du),length(funjac.du))
    _du = funjac.du

    funjac.jac(convert(Matrix, J), _du,
               _u, funjac.p, cj, t)
    return IDA_SUCCESS
end

function idajac(t::realtype,
                cj::realtype,
                u::N_Vector,
                du::N_Vector,
                res::N_Vector,
                _J::SUNMatrix,
                funjac::AbstractFunJac{<:SparseMatrixCSC},
                tmp1::N_Vector,
                tmp2::N_Vector,
                tmp3::N_Vector)

  jac_prototype = funjac.jac_prototype
  J = convert(SparseMatrixCSC,_J)

  funjac.u = unsafe_wrap(Vector{Float64}, __N_VGetArrayPointer_Serial(u),length(funjac.u))
  _u = funjac.u
  funjac.du = unsafe_wrap(Vector{Float64}, __N_VGetArrayPointer_Serial(du),length(funjac.du))
  _du = funjac.du

  funjac.jac(jac_prototype, _du, convert(Vector, _u), funjac.p, cj, t)
  J.nzval .= jac_prototype.nzval
  # Sundials resets the value pointers each time, so reset it too
  @. J.rowval = jac_prototype.rowval - 1
  @. J.colptr = jac_prototype.colptr - 1

  return IDA_SUCCESS
end

function massmat(t::Float64,
                 _M::SUNMatrix,
                 mmf::AbstractFunJac,
                 tmp1::N_Vector,
                 tmp2::N_Vector,
                 tmp3::N_Vector)
  if typeof(mmf.mass_matrix) <: Array
    M = convert(Matrix, _M)
  else
    M = convert(SparseMatrixCSC, _M)
  end
  M .= mmf.mass_matrix

  return IDA_SUCCESS
end

function jactimes(v::N_Vector,
                 Jv::N_Vector,
                 t::Float64,
                 y::N_Vector,
                 fy::N_Vector,
                 fj::AbstractFunJac,
                 tmp::N_Vector)
    DiffEqBase.update_coefficients!(fj.jac_prototype,y,fj.p,t)
    mul!(convert(Vector,Jv),fj.jac_prototype,convert(Vector,v))
    return CV_SUCCESS
end

function idajactimes(
                 t::Float64,
                 y::N_Vector,
                 fy::N_Vector,
                 r::N_Vector,
                 v::N_Vector,
                 Jv::N_Vector,
                 cj::Float64,
                 fj::AbstractFunJac,
                 tmp1::N_Vector,
                 tmp2::N_Vector)
    DiffEqBase.update_coefficients!(fj.jac_prototype,y,fj.p,t)
    mul!(convert(Vector,Jv),fj.jac_prototype,convert(Vector,v))
    return IDA_SUCCESS
end

function precsolve(t::Float64,
                   y::N_Vector,
                   fy::N_Vector,
                   r::N_Vector,
                   z::N_Vector,
                   gamma::Float64,
                   delta::Float64,
                   lr::Int,
                   fj::AbstractFunJac)
    fj.prec(convert(Vector,z),convert(Vector,r),fj.p,t,convert(Vector,y),convert(Vector,fy),gamma,delta,lr)
    return CV_SUCCESS
end

function precsetup(t::Float64,
                   y::N_Vector,
                   fy::N_Vector,
                   jok::Int,
                   jcurPtr::Ref{Int},
                   gamma::Float64,
                   fj::AbstractFunJac)
    fj.psetup(fj.p,t,convert(Vector,y),convert(Vector,fy),jok==1,Base.unsafe_wrap(Vector{Int}, jcurPtr, 1),gamma)
    return CV_SUCCESS
end

function idaprecsolve(t::Float64,
                   y::N_Vector,
                   fy::N_Vector,
                   resid::N_Vector,
                   r::N_Vector,
                   z::N_Vector,
                   gamma::Float64,
                   delta::Float64,
                   lr::Int,
                   fj::AbstractFunJac)
    fj.prec(convert(Vector,z),convert(Vector,r),fj.p,t,convert(Vector,y),convert(Vector,fy),convert(Vector,resid),gamma,delta,lr)
    return IDA_SUCCESS
end

function idaprecsetup(t::Float64,
                      y::N_Vector,
                      fy::N_Vector,
                      rr::N_Vector,
                      gamma::Float64,
                      fj::AbstractFunJac)
    fj.psetup(fj.p,t,convert(Vector,rr),convert(Vector,y),convert(Vector,fy),gamma)
    return IDA_SUCCESS
end


text:
#=
export DirectProductOperation
export isidentity
export directproduct
export ×ˢ

"""
    DirectProductOperation{Ops}

Direct product of symmetry operations.
math```
  g = (g_1, g_2, \\ldots)
```
"""
struct DirectProductOperation{Ops<:Tuple{Vararg{AbstractSymmetryOperation}}}<:AbstractSymmetryOperation
    operations::Ops
    function DirectProductOperation(ops::Vararg{AbstractSymmetryOperation})
        T = typeof(ops)
        return new{T}(ops)
    end
    function DirectProductOperation(ops::T) where {T<:Tuple{Vararg{AbstractSymmetryOperation}}}
        return new{T}(ops)
    end
end

Base.hash(x::P, h::UInt) where {P<:DirectProductOperation} = Base.hash(P, Base.hash(x.operations, h))

Base.:(*)(lhs::P, rhs::P) where {P<:DirectProductOperation} = DirectProductOperation(lhs.operations .* rhs.operations)
Base.:(^)(obj::DirectProductOperation, n::Integer) = DirectProductOperation(obj.operations.^n)
Base.:(==)(lhs::P, rhs::P) where {P<:DirectProductOperation} = all(lhs.operations .== rhs.operations)

Base.inv(obj::DirectProductOperation) = DirectProductOperation(Base.inv.(obj.operations))

function Base.isapprox(lhs::P, rhs::P; atol::Real=0, rtol::Real=Base.rtoldefault(Float64)) where {P<:DirectProductOperation}
    return all(isapprox.(lhs.operations, rhs.operations; atol=atol, rtol=rtol))
end

isidentity(obj::DirectProductOperation) = all(isidentity, obj.operations)

×ˢ(lhs::AbstractSymmetryOperation, rhs::AbstractSymmetryOperation) = DirectProductOperation(lhs, rhs)
×ˢ(lhs::DirectProductOperation, rhs::AbstractSymmetryOperation) = DirectProductOperation(lhs.operations..., rhs)
×ˢ(lhs::AbstractSymmetryOperation, rhs::DirectProductOperation) = DirectProductOperation(lhs, rhs.operations...)
×ˢ(lhs::DirectProductOperation, rhs::DirectProductOperation) = DirectProductOperation(lhs.operations..., rhs.operations...)
function ×ˢ(lhs::AbstractVector{<:AbstractSymmetryOperation}, rhs::AbstractVector{<:AbstractSymmetryOperation})
    return [×ˢ(l, r) for l in lhs, r in rhs]
end

directproduct(lhs::AbstractSymmetryOperation, rhs::AbstractSymmetryOperation) = DirectProductOperation(lhs, rhs)
directproduct(lhs::DirectProductOperation, rhs::AbstractSymmetryOperation) = DirectProductOperation(lhs.operations..., rhs)
directproduct(lhs::AbstractSymmetryOperation, rhs::DirectProductOperation) = DirectProductOperation(lhs, rhs.operations...)
directproduct(lhs::DirectProductOperation, rhs::DirectProductOperation) = DirectProductOperation(lhs.operations..., rhs.operations...)
function directproduct(lhs::AbstractVector{<:AbstractSymmetryOperation}, rhs::AbstractVector{<:AbstractSymmetryOperation})
    return [directproduct(l, r) for l in lhs, r in rhs]
end

function directproduct(::Type{E}, products::Function...) where {E<:DirectProductOperation}
    function product(lhs::E, rhs::E)
        return DirectProductOperation([p(l, r) for (p, l, r) in zip(products, lhs.operations, rhs.operations)]...)
    end
    return product
end
=#


text:
# ---
# title: 494. Target Sum
# id: problem494
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming, Depth-first Search
# link: <https://leetcode.com/problems/target-sum/description/>
# hidden: true
# ---
#
# You are given a list of non-negative integers, a1, a2, ..., an, and a target,
# S. Now you have 2 symbols `+` and `-`. For each integer, you should choose one
# from `+` and `-` as its new symbol.
#
# Find out how many ways to assign symbols to make sum of integers equal to
# target S.
#
# **Example 1:**
#
#
#
#     Input: nums is [1, 1, 1, 1, 1], S is 3.
#     Output: 5
#     Explanation:
#
#     -1+1+1+1+1 = 3
#     +1-1+1+1+1 = 3
#     +1+1-1+1+1 = 3
#     +1+1+1-1+1 = 3
#     +1+1+1+1-1 = 3
#
#     There are 5 ways to assign symbols to make the sum of nums be target 3.
#
#
#
#
# **Constraints:**
#
#   * The length of the given array is positive and will not exceed 20.
#   * The sum of elements in the given array will not exceed 1000.
#   * Your output answer is guaranteed to be fitted in a 32-bit integer.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# ======================================================================================================
# IO
function coefnames(aov::AnovaResult{T, FTest}; kwargs...) where {T <: TableRegressionModel{<: Union{LinearModel, GeneralizedLinearModel}}}
    v = coefnames(aov.model, Val(:anova))
    push!(v, "(Residuals)")
    v
end

coefnames(trm::TableRegressionModel{<: Union{LinearModel, GeneralizedLinearModel}}, anova::Val{:anova}) =
    coefnames(trm.mf, anova)

function coefnames(aov::AnovaResult{T, FTest}) where {T <: TableRegressionModel{<: Union{LinearModel, GeneralizedLinearModel}}}
    v = coefnames(aov.model, Val(:anova))
    push!(v, "(Residuals)")
    v
end

coefnames(trm::TableRegressionModel{<: Union{LinearModel, GeneralizedLinearModel}}, anova::Val{:anova}) =
    vectorize(coefnames(trm.mf.f.rhs, anova))

# anovatable api
function _anovatable(aov::AnovaResult{<: TableRegressionModel{<: LinearModel}, FTest}; kwargs...)
    AnovaTable(hcat(vectorize.((dof(aov), deviance(aov), deviance(aov) ./ dof(aov), teststat(aov), pval(aov)))...),
              ["DOF", "Exp.SS", "Mean Square", "F value","Pr(>|F|)"],
              ["x$i" for i in eachindex(pval(aov))], 5, 4)
end

function _anovatable(aov::AnovaResult{<: TableRegressionModel{<: GeneralizedLinearModel}, FTest}; kwargs...)
    AnovaTable(hcat(vectorize.((dof(aov), deviance(aov), deviance(aov) ./ dof(aov), teststat(aov), pval(aov)))...),
              ["DOF", "ΔDeviance", "Mean ΔDev", "F value","Pr(>|F|)"],
              ["x$i" for i in eachindex(pval(aov))], 5, 4)
end

function _anovatable(aov::AnovaResult{<: TableRegressionModel{<: LinearModel}, LRT}; kwargs...)
    AnovaTable(hcat(vectorize.((dof(aov), deviance(aov), teststat(aov), pval(aov)))...),
              ["DOF", "Res.SS", "χ²", "Pr(>|χ²|)"],
              ["x$i" for i in eachindex(pval(aov))], 4, 3)
end

function _anovatable(aov::AnovaResult{<: TableRegressionModel{<: GeneralizedLinearModel}, LRT}; kwargs...)
    AnovaTable(hcat(vectorize.((dof(aov), deviance(aov), teststat(aov), pval(aov)))...),
              ["DOF", "Deviance", "χ²", "Pr(>|χ²|)"],
              ["x$i" for i in eachindex(pval(aov))], 4, 3)
end

function _anovatable(aov::AnovaResult{<: Tuple, FTest},
                    modeltype1::Type{<: TableRegressionModel{<: LinearModel}},
                    modeltype2::Type{<: TableRegressionModel{<: LinearModel}};
                    kwargs...)

    rs = r2.(aov.model)
    Δrs = _diff(rs)
    AnovaTable(hcat(vectorize.((
                    dof(aov),
                    [NaN, _diff(dof(aov))...],
                    dof_residual(aov),
                    rs,
                    [NaN, Δrs...],
                    deviance(aov),
                    [NaN, _diffn(deviance(aov))...],
                    teststat(aov),
                    pval(aov)
                    ))...),
              ["DOF", "ΔDOF", "Res.DOF", "R²", "ΔR²", "Res.SS", "Exp.SS", "F value", "Pr(>|F|)"],
              ["$i" for i in eachindex(pval(aov))], 9, 8)
end

function _anovatable(aov::AnovaResult{<: Tuple, LRT},
                    modeltype1::Type{<: TableRegressionModel{<: LinearModel}},
                    modeltype2::Type{<: TableRegressionModel{<: LinearModel}};
                    kwargs...)
    rs = r2.(aov.model)
    Δrs = _diff(rs)
    AnovaTable(hcat(vectorize.((
                    dof(aov),
                    [NaN, _diff(dof(aov))...],
                    dof_residual(aov),
                    rs,
                    [NaN, Δrs...],
                    deviance(aov),
                    teststat(aov),
                    pval(aov)
                    ))...),
              ["DOF", "ΔDOF", "Res.DOF", "R²", "ΔR²", "Res.SS", "χ²", "Pr(>|χ²|)"],
              ["$i" for i in eachindex(pval(aov))], 8, 7)
end


text:
PlotlyBase.trace_map(p::SyncPlot, axis) = trace_map(p.plot, axis)
JSON.lower(sp::SyncPlot) = JSON.lower(sp.plot)

PlotlyBase._is3d(p::SyncPlot) = _is3d(p.plot)

# subplot methods on syncplot
Base.hcat(sps::TP...) where {TP<:SyncPlot} = TP(hcat([sp.plot for sp in sps]...))
Base.vcat(sps::TP...) where {TP<:SyncPlot} = TP(vcat([sp.plot for sp in sps]...))
Base.vect(sps::TP...) where {TP<:SyncPlot} = vcat(sps...)
Base.hvcat(rows::Tuple{Vararg{Int}}, sps::TP...) where {TP<:SyncPlot} =
    TP(hvcat(rows, [sp.plot for sp in sps]...))

function PlotlyBase.add_recession_bands!(p::SyncPlot; kwargs...)
    new_shapes = add_recession_bands!(p.plot; kwargs...)
    relayout!(p, shapes=new_shapes)
    new_shapes
end


text:
doc"""
A framework for pretraining neural nets (alternative to random weight initialization).
I expect to implement Deep Belief Net pretraining using Stacked Restricted Boltzmann Machines (RBM)
and Stacked Sparse (Denoising) Autoencoders.
"""
abstract PretrainStrategy

# default
function pretrain(net::NeuralNet, sampler::DataSampler, validationSampler::DataSampler; kwargs...)
  pretrain(DenoisingAutoencoder, net, sampler, validationSampler; kwargs...)
end

doc"Samples from the underlying sampler, but sets y = x.  Used in mapping inputs to themselves in autoencoders."
immutable AutoencoderDataSampler{T<:DataSampler} <: DataSampler
  sampler::T
end
StatsBase.sample(sampler::AutoencoderDataSampler) = (dp = sample(sampler.sampler); DataPoint(dp.x, dp.x))
DataPoints(sampler::AutoencoderDataSampler) = DataPoints([DataPoint(dp.x, dp.x) for dp in DataPoints(sampler.sampler)])

# -----------------------------------------------------------------

immutable DenoisingAutoencoder <: PretrainStrategy end

function pretrain(::Type{DenoisingAutoencoder}, net::NeuralNet, trainSampler::DataSampler, validationSampler::DataSampler;
                    tiedweights::Bool = true,
                    maxiter::Int = 10000,
                    dropout::DropoutStrategy = Dropout(pInput=0.7,pHidden=0.0),  # this is the "denoising" part, which throws out some of the inputs
                    encoderParams::NetParams = NetParams(dropout=dropout),
                    # solverParams::SolverParams = SolverParams(maxiter=maxiter, erroriter=typemax(Int), breakiter=typemax(Int)),  #probably don't set this manually??
                    solverParams::SolverParams = SolverParams(maxiter=maxiter, erroriter=2000, breakiter=2000, stopepochs=40),
                    inputMapping::Mapping = IdentityMapping())

  # lets pre-load the input dataset for simplicity... just need the x vec, since we're trying to map: x --> somthing --> x
  # trainData = DataPoints([DataPoint(dp.x, dp.x) for dp in DataPoints(trainSampler)])
  # println(trainData)
  # trainSampler = SimpleSampler(trainData)

  # set up training data
  trainData = DataPoints(AutoencoderDataSampler(trainSampler))

  # set up validation data
  validationData = DataPoints(AutoencoderDataSampler(validationSampler))

  # for each layer (which is not the output layer), fit the weights/bias as guided by the pretrain strategy
  for layer in net.layers[1:end-1]

    # NOTE: we are mapping input --> hidden --> input, which is why the "center" activation is called hiddenMapping
    #       and the "final" activation is called inputMapping

    # some setup
    hiddenMapping = layer.activation
    inputTransformer = layer === first(net.layers) ? net.inputTransformer : IdentityTransformer()

    # build a neural net which maps: nin -> nout -> nin
    autoencoder = buildNet(layer.nin, layer.nin, [layer.nout];
                            hiddenMapping = hiddenMapping,
                            finalMapping = inputMapping,
                            params = encoderParams,
                            solverParams = solverParams,
                            inputTransformer = inputTransformer)

    # tied weights means w₂ = w₁' ... rebuild the layer with a TransposeView of the first layer's weights
    # l1, l2 = autoencoder.layers
    if tiedweights
      setSecondAutoencoderLayer(autoencoder, autoencoder.layers...)
      # gradientState = ParameterUpdaterState(encoderParams.updater, l2.nin, l2.nout)
      # autoencoder.layers[2] = Layer(l2.nin, l2.nout, l2.activation, l2.p,
      #                               l2.dwState,
      #                               l2.dbState,
      #                               l2.x, TransposeView(l1.w), l2.b, l2.δ,
      #                               l2.Σ, l2.a, l2.r, l2.nextr)
    end

    println("netlayer: $layer  oact: $inputMapping")
    println("autoenc: $autoencoder")

    # solve for the weights and bias... note we're not using stopping criteria... only maxiter
    stats = solve!(autoencoder, SimpleSampler(trainData), SimpleSampler(validationData), true)
    println("  $stats")

    if stats.bestModel == nothing || isnan(stats.bestValidationError) || isinf(stats.bestValidationError)
      warn("Somthing wrong with pretraining model: $stats")
    end

    # save the weights and bias to the layer
    # println("l1: $l1")

    l1 = first(stats.bestModel.layers)
    @assert isa(l1.w, Matrix)
    println(l1)
    layer.w = copy(l1.w)
    layer.b = copy(l1.b)

    # update the inputMapping, so that this layer's activation becomes the next autoencoder's inputMapping
    inputMapping = hiddenMapping

    # feed the data forward to the next layer
    for i in 1:length(trainData)
      x = transform(inputTransformer, trainData[i].x)
      forward!(l1, x, false)
      newx = copy(l1.a)
      trainData[i] = DataPoint(newx, newx)
    end

    # do the same for the validation set
    for i in 1:length(validationData)
      x = transform(inputTransformer, validationData[i].x)
      forward!(l1, x, false)
      newx = copy(l1.a)
      validationData[i] = DataPoint(newx, newx)
    end

    # println(trainData)

  end

  # we're done... net is pretrained now
  return
end

function setSecondAutoencoderLayer(autoencoder::NeuralNet, l1::Layer, l2::Layer)
  autoencoder.layers[2] = Layer(l2.nin, l2.nout, l2.activation, l2.p,
                                l2.dwState,
                                l2.dbState,
                                l2.x, TransposeView(l1.w), l2.b, l2.δ,
                                l2.Σ, l2.a, l2.r, l2.nextr)
end

function setSecondAutoencoderLayer(autoencoder::NeuralNet, l1::NormalizedLayer, l2::NormalizedLayer)
  autoencoder.layers[2] = NormalizedLayer(l2.nin, l2.nout, l2.activation, l2.p,
                                l2.dwState,
                                l2.dbState,
                                l2.dβState,
                                l2.dαState,
                                l2.xvar, l2.w,
                                l2.x, l2.xhat, l2.α, l2.y, l2.δy,
                                l2.b, l2.Σ, l2.a, l2.δΣ,
                                l2.β, l2.r, l2.nextr)
end


text:
export Ball

struct Ball{T <: Real} <: ConvexSet
    center::Vector{T}
    radius::T
end

"""
    Ball()
    Ball(c, r)

A ball with center `c` and radius `r`.
"""
function Box(c::Vector{T}, r::T) where T <: Real
    r ≥ 0 || throw(ArgumentError("radius should be positive"))
    Ball(c, r)
end

"""
    Ball(n)
    Ball(n, r)

A `n` dimensional ball with center `zeros(n)` and radius `r`.
"""
Ball(n::Integer, r::Real=1.0) = Ball(zeros(typeof(r), n), r)

function project!(s::Ball{T}, v::Vector{T}, y::Vector{T}) where T <: Real
    v .= y .- s.center
    distance = sqrt(sum(abs2, v))
    if distance > s.radius
        v .= s.center .+ (s.radius / distance) .* v
    else
        copyto!(v, y)
    end
    v
end

function project(s::Ball{T}, y::Vector{T}) where T <: Real
    project!(s, similar(y), y)
end


text:
using Gadfly, RDatasets

set_default_plot_size(6inch, 3inch)

# see issue #880
plot(dataset("ggplot2", "diamonds"), x="Price", color="Cut",
                   Geom.histogram(bincount=30, density=true))


text:
# This function takes a graph and a set of indices
# and returns an iterator over assignments
# i -> e[1], j -> e[2], k -> e[3], ...
# for all n-edges of the graph, where n is the
# number of indices. The assignments should be
# of type Dict{Symbol,Symbol}.
#
# For an ordinary graph, the 1-edges are the vertices
# and the 2-edges are the edges in the usual sense.
#
# Important: as the order of the indices is not
# defined, this function should assign elements
# to all permutations of the edges, ie. for
# i => 4, j => 5 it should also return i => 5, j => 4
function graph_iterator end

struct BasicGraph{T}
    verts::AbstractVector{T}
    edges::AbstractVector{Pair{T,T}}
end

function graph_iterator(graph::BasicGraph{T}, idcs::Set{Symbol}) where {T}
    if length(idcs) == 1
        i = first(idcs)
        return ( Dict(i => vert) for vert in graph.verts )
    elseif length(idcs) == 2
        i, j = collect(idcs)
        edges_r = ( edge.second => edge.first for edge in graph.edges )
        edges_sym = union(graph.edges, edges_r)

        return ( Dict(i => edge.first, j => edge.second) for edge in edges_sym )
    end

    return nothing
end

struct PeriodicGrid{N}
    size::NTuple{N, Int64}
end

PeriodicGrid(sizes...) = PeriodicGrid(sizes)

function get_neighbour(coords, dim::Int64, amt::Int64, size::NTuple{N,Int64})::Vector{Int64} where N
    ret = collect(coords)
    ret[dim] += amt
    ret[dim] = mod(ret[dim] - 1, size[dim]) + 1    # transform to lie in 1:size[dim]
    ret
end

coord_to_str(coords) = "(" * join(coords, ",") * ")"
grid_coords(grid::PeriodicGrid{N}) where N = Iterators.product([ 1:grid.size[i] for i in 1:N ]...)

function get_edge_iter(i::Symbol, j::Symbol, dim::Int64, grid::PeriodicGrid)
    return Iterators.flatten(( Dict(i => coord_to_str(coords),
                                      j => coord_to_str(get_neighbour(coords, dim, sgn, grid.size)))
                               for coords in grid_coords(grid) )
                               for sgn in [-1,+1])
end

function graph_iterator(grid::PeriodicGrid{N}, idcs::Set{Symbol}) where N
    if length(idcs) == 1
        idx = first(idcs)
        return ( Dict(idx => coord_to_str(coords)) for coords in grid_coords(grid) )
    elseif length(idcs) == 2
        i, j = collect(idcs)
        return Iterators.flatten(get_edge_iter(i, j, dim, grid) for dim in 1:N if grid.size[dim] != 1)
    end

    return nothing
end

# Convert a list of expressions of the form
# G[i], H[j,k] into a dictionary
# of the form Dict(i => G, j => H, k => H)
function get_idx_dict(exprs)::Dict{Symbol,Any}
    ret = Dict{Symbol,Any}()

    for idx_map in exprs
        (isa(idx_map, Expr) && idx_map.head == :ref) || throw("malformed index assignment: $(idx_map)")
        graph_name = idx_map.args[1]
        graph = Base.eval(@__MODULE__, graph_name)
        idcs = idx_map.args[2:end]

        for idx in idcs
            idx isa Symbol || throw("index not a symbol: $(idx)")
            idx in keys(ret) && throw("multiple index assignment: $(idx)")
            ret[idx] = graph
        end
    end

    return ret
end

# Iterate over all valid combinations of indices for a spatial reaction
function iterate_idcs(idcs::Set{Symbol}, idx_dict::AbstractDict)
    isempty(idcs) && return (Dict{Symbol,Any}(),)

    iterators = []
    for graph in unique(values(idx_dict))
        idcs_graph = filter(idx -> idx_dict[idx] == graph, idcs)
        isempty(idcs_graph) && continue
        iterator = graph_iterator(graph, idcs_graph)
        iterator === nothing && throw("graph iteration not implemented for $(length(idcs_graph)) indices")
        push!(iterators, iterator)
    end

    ( merge(dicts...) for dicts in Iterators.product(iterators...) )
end

# Get indices occurring in a reaction
function get_idcs(ex::ExprValues)::Set{Symbol}
    if ex isa Expr
        if ex.head == :ref
            return Set(ex.args[2:end])
        end

        return mapreduce(get_idcs, union, ex.args)
    end

    return Set{Symbol}()
end

# Expand spatial reactants of the form X[i] with the given assignment of
# vertices to indices
function expand_reactant_idcs(ex::ExprValues, idx_map::AbstractDict)
    if ex isa Expr
        if ex.head == :ref
            name = ex.args[1]
            name isa Symbol || throw("only species can be indexed")

            length(ex.args) == 2 || throw("species $(ex.args[1]) can only have one index")

            idx = ex.args[2]
            idx_mapped = idx_map[idx]

            return Symbol(name, "_", idx_mapped)
        end

        return Expr(ex.head, map(arg -> expand_reactant_idcs(arg, idx_map), ex.args)...)
    end

    ex
end

# Expand reaction rates of the form a[i,j] with the given assignment of
# vertices to indices
function expand_rate_idcs(ex::ExprValues, idx_map::AbstractDict)
    if ex isa Expr
        if ex.head == :ref
            name = ex.args[1]
            name isa Symbol || throw("only constant parameters can be indexed")

            idcs = ex.args[2:end]
            idcs_mapped = map(x -> idx_map[x], idcs)
            if length(idcs) == 1
                return Symbol(name, "_", idcs_mapped[1])
            else
                return Symbol(name, "_(", join(idcs_mapped, ","), ")")
            end
        end

        return Expr(ex.head, map(arg -> expand_rate_idcs(arg, idx_map), ex.args)...)
    end

    ex
end

function expand_spatial_reactions!(reactions::AbstractVector{Expr}, rate::ExprValues, sub_line::ExprValues, prod_line::ExprValues, arrow::Symbol,
                                   idx_dict::AbstractDict)
    idcs_sub = get_idcs(sub_line)
    idcs_prod = get_idcs(prod_line)
    idcs = union(idcs_sub, idcs_prod)

    issubset(idcs, keys(idx_dict)) || throw("invalid index in reaction")

    idcs_rates = get_idcs(rate)
    issubset(idcs_rates, idcs) || throw("invalid index in reaction rate")

    for idx_map in iterate_idcs(idcs, idx_dict)
        rate_mapped = expand_rate_idcs(rate, idx_map)
        sub_line_mapped = expand_reactant_idcs(sub_line, idx_map)
        prod_line_mapped = expand_reactant_idcs(prod_line, idx_map)

        r_line = Expr(:call, arrow, sub_line_mapped, prod_line_mapped)

        push!(reactions, Expr(:tuple, rate_mapped, r_line))
    end
end

function get_spatial_reactions(ex::Expr, idx_dict::AbstractDict)
    reactions = Expr[]

    for line in ex.args
        (line isa Expr && line.head == :tuple) || continue

        (rate, r_line) = line.args
        (r_line.head  == :-->) && (r_line = Expr(:call,:→,r_line.args[1],r_line.args[2]))

        expand_spatial_reactions!(reactions, rate, r_line.args[2], r_line.args[3], r_line.args[1], idx_dict)
    end

    Expr(:block, reactions...)
end

function make_spatial_reaction_network(ex::Expr, graphs)
    idx_dict = get_idx_dict(graphs)
    get_spatial_reactions(MacroTools.striplines(ex), idx_dict)
end

macro spatial_reaction_network(ex::Expr, graphs...)
    esc(make_spatial_reaction_network(ex, graphs))
end


text:
push!(LOAD_PATH,"../../src/")
import UNSflow

c = 0.027

u = 10.

AR = 28.14814814814815

data = JLD.load("data.jld", "data")

alphadef = UNSflow.ConstDef(1.)

hdef = UNSflow.ConstDef(0.)

udef = UNSflow.ConstDef(1.)

full_kinem = UNSflow.KinemDef(alphadef, hdef, udef)

pvt = 0.5

geometry = "FlatPlate"

lespcrit = [0.02;]

#surf = UNSflow.ThreeDSurfSimple(AR, full_kinem, geometry, pvt, lespcrit, nspan=10, cref=c, uref=u)

#curfield = UNSflow.ThreeDFieldSimple()

eta = zeros(4)
etad = zeros(4)
#kinem = UNSflow.KinemCantilever(eta, etad)

dtstar = 0.015

nsteps = 40000

t_tot = nsteps*dtstar*c/u

startflag = 1

writeflag = 1

writeInterval = t_tot/8.

delvort = UNSflow.delSpalart(500, 12, 1e-5)
#delvort = UNSflow.delNone()

mat, surf, curfield = UNSflow.stripCantilever(surf, curfield, data, kinem, nsteps, dtstar,startflag, writeflag, writeInterval, delvort)

#UNSflow.makeVortPlots3D()

#UNSflow.makeForcePlots3D()

#UNSflow.makeInfoPlots3D()

UNSflow.cleanWrite()


text:
"""
    legendre(n, x)

Compute the Legendre polynomial of order n at point x.

The Legendre polynomials are a special case of the Jacobi polynomials with
α, β = 0.  The implementation here reflects that.
"""
function legendre(n, x)
    return Jacobi.jacobi(n, 0, 0, x)
end


text:
import Base: push!, pop!, length

"""
    AbstractBit

Inherit from this type to create a custom bit type.
"""
abstract type AbstractBit <: Number end

"""
    AbstractBitstream

Inherit from this type to create a custom bitstream type.

Expected fields:
- `bits::Queue{AbstractBit}`: the underlying bitstream
"""
abstract type AbstractBitstream end

"""
    push!(s::AbstractBitstream, b)

Push a bit(s) `b` onto bitstream `s`.

Fields:
- `s::AbstractBitstream`: the bitstream object
- `b`: the bit(s) to push onto the stream
"""
push!(s::AbstractBitstream, b::AbstractBit) = enqueue!(s.bits, b)
function push!(s::AbstractBitstream, bs::Vector{<:AbstractBit})
    for b in bs
        push!(s, b)
    end
end

"""
    pop!(s::AbstractBitstream)

Pop a bit from bitstream `s`.

Fields:
- `s::AbstractBitstream`: the bitstream object
"""
pop!(s::AbstractBitstream) = dequeue!(s.bits)

"""
    observe(s::AbstractBitstream)

Examine the most recent bit added to the stream without removing it.

Fields:
- `s::AbstractBitstream`: the bitstream object
"""
observe(s::AbstractBitstream) = last(s.bits)

"""
    length(s::AbstractBitstream)

Return the number of bits in `s`.
"""
length(s::AbstractBitstream) = length(s.bits)


text:
using ChaosTools, Test, DelayEmbeddings, DelimitedFiles, Downloads

@testset "exponential-chirp" begin
	file = "exponential-chirp-440-880-5.csv" #generated initially with y = librosa.chirp(440, 880, duration=5.0)
	println("Downloading 2MB file into $(file) in current dir.")
	Downloads.download("https://raw.githubusercontent.com/JuliaDynamics/JuliaDynamics/master/timeseries/$(file)", file)
	y = readdlm(file)[:,1]
	sr = Int64(length(y)/5)
	F0s, time_scale = yin(y, sr, f0_min=200, f0_max=1000, f_step=512, w_len=1024)
	t = range(0,5,length=length(F0s))
	chirp_freqs = 440 .*exp.(log(2)/5 .* t)
	@test all(abs.(F0s .- chirp_freqs) .< 5)
	println("Removing the file.")
	rm(file)
end


text:
# ====================================================================
# Load Different Logger
# ====================================================================
"""
    load_logger_from_config(file_path::AbstractString)::TeeLogger

Create a combined logger from a config file path.

### Note
A combined logger is a `TeeLogger` struct that allows to send the same log message to
all loggers included in the combined logger at once.
"""
function load_logger_from_config(file_path::AbstractString)
    json_string = read(file_path, String)
    json3_dict = JSON3.read(json_string) #, Dict
    return load_logger_from_config(json3_dict)
end

"""
    load_logger_from_config(config::Dict)::TeeLogger

Create a `TeeLoger`, a collection of logger, from a configuration `Dict`.

### Note
-  A `TeeLogger` struct allows to send the same log message  to  all loggers included in the `TeeLoger` at once.
- The configuration `Dict` requires as `"logging"` field, see example for more details.

### Example
```julia
config = Dict(
    ...,
    "logging" => [
        # logger_type => logger_config_dict,
        "SeqLogger" => Dict(...),
        ...
    ],
    ...
)
```

### Returns
`TeeLogger` as defined in `config`
"""
function load_logger_from_config(config::AbstractDict)
    loggers = [
        get_logger(string(logger_type), logger_specs) for (logger_type, logger_specs) in config["logging"]
    ]
    return TeeLogger(loggers...)
end

# ====================================================================
# Load Single Logger
# ====================================================================
"""
    load_seqlogger(logger_config::AbstractDict)::Union{SeqLogger, TransformerLogger}

Return a `SeqLogger` or `TransformerLogger` according to `logger_config`.

### Config Parameters
- `"server_url"` -- required
- `"min_level"` -- required (`"DEBUG", "INFO", "WARN", "ERROR"`)
- `"transformation"` -- optional, default `identity`
- `"api_key"` -- optional, default `""`
- `"batch_size"` -- optional, default `10`

All other config parameters are used as global event properties.

### Example
```julia
log_dict = Dict(
    "server_url" => "http://subdn215:5341/",
    "min_level" => "INFO",
    "batch_size" => 12,
    "App" => "SeqLoggers_Test",
    "Env" => "Test"
)
seq_logger = SeqLoggers.load_seqlogger(log_dict)
```
"""
function load_seqlogger(logger_config::AbstractDict)
    server_url = logger_config["server_url"]
    min_level = logger_config["min_level"] |> get_log_level
    transformation_str = get(logger_config, "transformation", "identity")
    transformation = transformation_str |> get_transformation_function

    # Create a NamedTuple from remaining Config Keyvalue Pairs
    kwarg_keys = filter(
        key -> string(key) ∉ ["server_url", "min_level", "transformation"], keys(logger_config)
    )
    kwarg_keys_names = Tuple(Symbol(key) for key in kwarg_keys)
    kwarg_keys_values = [logger_config[key] for key in kwarg_keys]
    kwarg_keys = NamedTuple{kwarg_keys_names}(kwarg_keys_values)
    return SeqLogger(
        server_url; min_level=min_level, kwarg_keys...
    ) |> transformation
end

"""
    load_consolelogger(logger_config::AbstractDict)::Union{ConsoleLogger, TransformerLogger}

Return a `ConsoleLogger` or `TransformerLogger` according to `logger_config`.

### Config Parameters
- `"min_level"` -- required (`"DEBUG", "INFO", "WARN", "ERROR"`)
- `"transformation"` -- optional, default `identity`

### Example
```julia
logging_config = Dict(
    "min_level" => "ERROR",
    "transformation" => "add_timestamp",
)

seq_logger = SeqLoggers.load_consolelogger(log_dict)
```
"""
function load_consolelogger(logger_config::AbstractDict)
    min_level = logger_config["min_level"] |> get_log_level
    transformation_str = get(logger_config, "transformation", "identity")
    transformation = transformation_str |> get_transformation_function
    return ConsoleLogger(stderr, min_level) |> transformation
end

"""
    load_consolelogger(logger_config::AbstractDict)::AbstractLogger

Return a `MinLevelLogger{FileLogger}` or `TransformerLogger` according to `logger_config`.

### Config Parameters
- `"file_path"` -- required
- `"min_level"` -- required (`"DEBUG", "INFO", "WARN", "ERROR"`)
- `"append"` -- optional, default `true`, append to file if `true`, otherwise truncate file. (See [`LoggingExtras.FileLogger`](@ref) for more information.)
- `"transformation"` -- optional, default `identity`

### Example
```julia
logging_config = Dict(
    "file_path" => "C:/Temp/test.log",
    "min_level" => "ERROR",
    "append" => true,
    "transformation" => "add_timestamp",
)
seq_logger = SeqLoggers.load_filelogger(log_dict)
```
"""
function load_filelogger(logger_config::AbstractDict)
    min_level = logger_config["min_level"] |> get_log_level
    file_path = logger_config["file_path"]
    append = get(logger_config, "append", true)
    transformation_str = get(logger_config, "transformation", "identity")
    transformation = transformation_str |> get_transformation_function
    return MinLevelLogger(FileLogger(file_path; append=append), min_level) |> transformation
end

"""
    load_advanced_filelogger(logger_config::AbstractDict)::AbstractLogger

Return a `DatetimeRotatingFileLogger` or `TransformerLogger` according to `logger_config`.

### Config Parameters
- `"dir_path"` -- required
- `"min_level"` -- required (`"DEBUG", "INFO", "WARN", "ERROR"`)
- `"file_name_pattern"` -- required e.g. `"\\a\\c\\c\\e\\s\\s-YYYY-mm-dd-HH-MM.\\l\\o\\g"`
- `"transformation"` -- optional, default `identity`

### Example
```julia
logging_config = Dict(
    "dir_path" => "C:/Temp",
    "min_level" => "ERROR",
    "file_name_pattern" => "\\a\\c\\c\\e\\s\\s-YYYY-mm-dd-HH-MM.\\l\\o\\g",
    "transformation" => "add_timestamp",
)
seq_logger = SeqLoggers.load_advanced_filelogger(log_dict)
```
"""
function load_advanced_filelogger(logger_config::AbstractDict)
    min_level = logger_config["min_level"] |> get_log_level
    dir_path = logger_config["dir_path"]
    file_name_pattern = logger_config["file_name_pattern"]

    transformation_str = get(logger_config, "transformation", "identity")
    transformation = transformation_str |> get_transformation_function

    return AdvancedFileLogger(
        dir_path,
        file_name_pattern;
        log_format_function=print_standard_format,
        min_level=min_level
    ) |> transformation
end

# ====================================================================
# Logger Type Mapping and Register
# ====================================================================
"""
    get_logger(logger_type::AbstractString, logger_config::AbstractDict)::AbstractLogger

Create logger struct from logger type name and `Dict` with required parameters.

By default, the following logger types are supported:
- `"SeqLogger"` → [`SeqLogger`](@ref)
- `"ConsoleLogger"` → [`ConsoleLogger`](@ref)
- `"FileLogger"` → [`FileLogger`](@ref)

Use [`register_logger!`](@ref) to add custom `AbstractLogger`s.
"""
function get_logger(logger_type::AbstractString, logger_config::AbstractDict)
    logger_constructor = get(LOGGER_TYPE_MAPPING, logger_type, nothing)
    if isnothing(logger_constructor)
        throw(
            ArgumentError(
                "There is no logger corresponding to the key `$logger_type`. " *
                "Available options are $(collect(keys(LOGGER_TYPE_MAPPING))). " *
                "Use `register_logger!` to add new logger types."
            )
        )
    end
    return logger_constructor(logger_config)
end

"""
    register_logger!(logger_type::AbstractString, logger_constructor::Function)

Register a new logger type.

Registering enables the user to use custom `AbstractLogger` struct, defined outside of `SeqLoggers`,
to be used with [`load_logger_from_config`](@ref).
"""
function register_logger!(logger_type::AbstractString, logger_constructor::Function)
    if haskey(LOGGER_TYPE_MAPPING, logger_type)
        @warn "Logger type `$logger_type` already exists and will be overwritten"
    end
    LOGGER_TYPE_MAPPING[logger_type] = logger_constructor
    return nothing
end

const LOGGER_TYPE_MAPPING = Dict(
    "SeqLogger" => load_seqlogger,
    "ConsoleLogger" => load_consolelogger,
    "FileLogger" => load_filelogger,
    "AdvancedFileLogger" => load_advanced_filelogger,
)

# ====================================================================
# Transformation Function and Register
# ====================================================================

"""
    get_transformation_function(key::String)

Convert a string (from config) into a transformation function.

By default, the following transformation functions are supported:
- `"identity"` → [`identity`](@ref): no transformation
- `"add_timestamp"` → [`add_timestamp`](@ref): add timestamp at the beginning of log message

Use [`register_transformation_function!`](@ref) to add custom transformation functions.
"""
function get_transformation_function(key::String)
    return LOGGER_TRANSFORMATION_MAPPING[key]
end

"""
    register_transformation_function!(key::AbstractString, transformation_function::Function)

Register new transformation function.

Registering enables the user to use custom transformation functions, defined outside of `SeqLoggers`,
to be used with [`load_logger_from_config`](@ref).
"""
function register_transformation_function!(key::AbstractString, transformation_function::Function)
    LOGGER_TRANSFORMATION_MAPPING[key] = transformation_function
    return nothing
end

const STANDARD_DATETIME_FORMAT = "yyyy-mm-dd HH:MM:SS"

"""
add_timestamp(logger::AbstractLogger)

Logger transformation function that prepends a timestamp to a logging message.
"""
add_timestamp(logger) = TransformerLogger(logger) do log
    merge(log, (; message="$(Dates.format(now(), STANDARD_DATETIME_FORMAT)): $(log.message)"))
end

const LOGGER_TRANSFORMATION_MAPPING = Dict(
    "identity" => identity,
    "add_timestamp" => add_timestamp,
)

# ====================================================================
# Log Level
# ====================================================================
"""
    get_log_level(key::String)::Logging.LogLevel

Return the `Loggin.LogLevel` corresponding to the input string.
"""
function get_log_level(key::String)
    log_level = get(LOG_LEVEL_MAPPING, uppercase(key), nothing)
    if isnothing(key)
        throw(
            ArgumentError(
                "There is no log level corresponding to the key $key." *
                "Available options are $(collect(keys(LOG_LEVEL_MAPPING)))"
            )
        )
    end
    return log_level
end

const LOG_LEVEL_MAPPING = Dict(
    "INFO" => Logging.Info,
    "INFORMATOIN" => Logging.Info,
    "DEBUG" => Logging.Debug,
    "WARN" => Logging.Warn,
    "WARNING" => Logging.Warn,
    "ERROR" => Logging.Error,
)


text:
doc"""
    VonMises(μ, κ)

The *von Mises distribution* with mean `μ` and concentration `κ` has probability density function

$f(x; \mu, \kappa) = \frac{1}{2 \pi I_0(\kappa)} \exp \left( \kappa \cos (x - \mu) \right)$

```julia
VonMises()       # von Mises distribution with zero mean and unit concentration
VonMises(κ)      # von Mises distribution with zero mean and concentration κ
VonMises(μ, κ)   # von Mises distribution with mean μ and concentration κ
```

External links

* [von Mises distribution on Wikipedia](http://en.wikipedia.org/wiki/Von_Mises_distribution)

"""
immutable VonMises{T<:Real} <: ContinuousUnivariateDistribution
    μ::T      # mean
    κ::T      # concentration
    I0κ::T    # I0(κ), where I0 is the modified Bessel function of order 0

    function (::Type{VonMises{T}}){T}(μ::T, κ::T)
        @check_args(VonMises, κ > zero(κ))
        new{T}(μ, κ, besseli(zero(T), κ))
    end
end

VonMises{T<:Real}(μ::T, κ::T) = VonMises{T}(μ, κ)
VonMises(μ::Real, κ::Real) = VonMises(promote(μ, κ)...)
VonMises(μ::Integer, κ::Integer) = VonMises(Float64(μ), Float64(κ))
VonMises(κ::Real) = VonMises(0.0, κ)
VonMises() = VonMises(0.0, 1.0)

show(io::IO, d::VonMises) = show(io, d, (:μ, :κ))

@distr_support VonMises d.μ - π d.μ + π

#### Conversions

convert{T<:Real}(::Type{VonMises{T}}, μ::Real, κ::Real) = VonMises(T(μ), T(κ))
convert{T<:Real, S<:Real}(::Type{VonMises{T}}, d::VonMises{S}) = VonMises(T(d.μ), T(d.κ))

#### Parameters

params(d::VonMises) = (d.μ, d.κ)
@inline partype{T<:Real}(d::VonMises{T}) = T

#### Statistics

mean(d::VonMises) = d.μ
median(d::VonMises) = d.μ
mode(d::VonMises) = d.μ
var(d::VonMises) = 1 - besseli(1, d.κ) / d.I0κ
# deprecated 12 September 2016
@deprecate circvar(d) var(d)
entropy(d::VonMises) = log(twoπ * d.I0κ) - d.κ * (besseli(1, d.κ) / d.I0κ)

cf(d::VonMises, t::Real) = (besseli(abs(t), d.κ) / d.I0κ) * cis(t * d.μ)

#### Evaluations

pdf(d::VonMises, x::Real) = exp(d.κ * cos(x - d.μ)) / (twoπ * d.I0κ)
logpdf(d::VonMises, x::Real) = d.κ * cos(x - d.μ) - log(d.I0κ) - log2π

cdf(d::VonMises, x::Real) = _vmcdf(d.κ, d.I0κ, x - d.μ, 1e-15)

function _vmcdf(κ::Real, I0κ::Real, x::Real, tol::Real)
    j = 1
    cj = besseli(j, κ) / j
    s = cj * sin(j * x)
    while abs(cj) > tol
        j += 1
        cj = besseli(j, κ) / j
        s += cj * sin(j * x)
    end
    return (x + 2s / I0κ) / twoπ + 1//2
end

#### Sampling

sampler(d::VonMises) = VonMisesSampler(d.μ, d.κ)


text:
include("utils.jl")

using HybridArrays, Random

Random.seed!(42)

@testset "Power manifold" begin

    Ms = Sphere(2)
    Ms1 = PowerManifold(Ms, 5)
    Ms2 = PowerManifold(Ms, 5, 7)
    Mr = Manifolds.Rotations(3)
    Mr1 = PowerManifold(Mr, 5)
    Mr2 = PowerManifold(Mr, 5, 7)

    types_s1 = [Array{Float64,2},
                HybridArray{Tuple{3,StaticArrays.Dynamic()}, Float64, 2}]
    types_s2 = [Array{Float64,3},
                HybridArray{Tuple{3,StaticArrays.Dynamic(),StaticArrays.Dynamic()}, Float64, 3}]

    types_r1 = [Array{Float64,3},
                HybridArray{Tuple{3,3,StaticArrays.Dynamic()}, Float64, 3}]
    types_r2 = [Array{Float64,4},
                HybridArray{Tuple{3,3,StaticArrays.Dynamic(),StaticArrays.Dynamic()}, Float64, 4}]

    retraction_methods = [Manifolds.PowerRetraction(ManifoldsBase.ExponentialRetraction())]
    inverse_retraction_methods = [Manifolds.InversePowerRetraction(ManifoldsBase.LogarithmicInverseRetraction())]

    sphere_dist = Manifolds.uniform_distribution(Ms, @SVector [1.0, 0.0, 0.0])
    power_s1_pt_dist = Manifolds.PowerPointDistribution(Ms1, sphere_dist, randn(Float64, 3, 5))
    power_s2_pt_dist = Manifolds.PowerPointDistribution(Ms2, sphere_dist, randn(Float64, 3, 5, 7))
    sphere_tv_dist = Manifolds.normal_tvector_distribution(Ms, (@MVector [1.0, 0.0, 0.0]), 1.0)
    power_s1_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Ms1), rand(power_s1_pt_dist), sphere_tv_dist)
    power_s2_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Ms2), rand(power_s2_pt_dist), sphere_tv_dist)

    id_rot = @SMatrix [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
    rotations_dist = Manifolds.normal_rotation_distribution(Mr, id_rot, 1.0)
    power_r1_pt_dist = Manifolds.PowerPointDistribution(Mr1, rotations_dist, randn(Float64, 3, 3, 5))
    power_r2_pt_dist = Manifolds.PowerPointDistribution(Mr2, rotations_dist, randn(Float64, 3, 3, 5, 7))
    rotations_tv_dist = Manifolds.normal_tvector_distribution(Mr, MMatrix(id_rot), 1.0)
    power_r1_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Mr1), rand(power_r1_pt_dist), rotations_tv_dist)
    power_r2_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Mr2), rand(power_r2_pt_dist), rotations_tv_dist)

    trim(s::String) = s[1:min(length(s), 20)]

    for T in types_s1
        @testset "Type $(trim(string(T)))..." begin
            pts1 = [convert(T, rand(power_s1_pt_dist)) for _ in 1:3]
            test_manifold(Ms1,
                          pts1;
                          test_reverse_diff = true,
                          test_musical_isomorphisms = true,
                          retraction_methods = retraction_methods,
                          inverse_retraction_methods = inverse_retraction_methods,
                          point_distributions = [power_s1_pt_dist],
                          tvector_distributions = [power_s1_tv_dist],
                          rand_tvector_atol_multiplier = 6.0)
        end
    end
    for T in types_s2
        @testset "Type $(trim(string(T)))..." begin
            pts2 = [convert(T, rand(power_s2_pt_dist)) for _ in 1:3]
            test_manifold(Ms2,
                          pts2;
                          test_reverse_diff = true,
                          test_musical_isomorphisms = true,
                          retraction_methods = retraction_methods,
                          inverse_retraction_methods = inverse_retraction_methods,
                          point_distributions = [power_s2_pt_dist],
                          tvector_distributions = [power_s2_tv_dist],
                          rand_tvector_atol_multiplier = 6.0)
        end
    end

    for T in types_r1
        @testset "Type $(trim(string(T)))..." begin
            pts1 = [convert(T, rand(power_r1_pt_dist)) for _ in 1:3]
            test_manifold(Mr1,
                          pts1;
                          test_reverse_diff = false,
                          test_musical_isomorphisms = true,
                          retraction_methods = retraction_methods,
                          inverse_retraction_methods = inverse_retraction_methods,
                          point_distributions = [power_r1_pt_dist],
                          tvector_distributions = [power_r1_tv_dist],
                          rand_tvector_atol_multiplier = 5.0)
        end
    end
    for T in types_r2
        @testset "Type $(trim(string(T)))..." begin
            pts2 = [convert(T, rand(power_r2_pt_dist)) for _ in 1:3]
            test_manifold(Mr2,
                          pts2;
                          test_reverse_diff = false,
                          test_musical_isomorphisms = true,
                          retraction_methods = retraction_methods,
                          inverse_retraction_methods = inverse_retraction_methods,
                          point_distributions = [power_r2_pt_dist],
                          tvector_distributions = [power_r2_tv_dist],
                          rand_tvector_atol_multiplier = 5.0)
        end
    end

end


text:
__precompile__()

module AutoUrban

using AutomotiveDrivingModels
using AutoViz
using EzXML
using Printf
using LinearAlgebra

include("simulation/main.jl")
include("drivermodels/main.jl")
include("render/main.jl")
include("roadway/main.jl")

end


text:
using DifferentialEquations
using DataInterpolations
using Plots
using ForwardDiff

#===============================================================#
abstract type DynamicalSystem end

struct LorenzSystem <: DynamicalSystem end
struct LotkaVolterra <: DynamicalSystem end
struct NSECylinder <: DynamicalSystem end
struct HelicopterData <: DynamicalSystem end

function datagen(data::LorenzSystem)
    function noise(u)
        for i = 1:length(u)
            for j = 1:length(u[i])
                u[i][j] += 0.1*randn()
            end
        end
        return u
    end

    function lorenz!(du,u,p,t)
        du[1] = p[1]*(u[2]-u[1])
        du[2] = u[1]*(p[2]-u[3]) - u[2]
        du[3] = u[1]*u[2] - p[3]*u[3]
    end

    u₀ = [1,0,0]
    p = [10,12,8/3] # Change d_lorenz
    tspan = (0,100)
    t = 0.01:0.01:100
    prob = ODEProblem(lorenz!,u₀,tspan,p)
    sol = solve(prob,Tsit5(),saveat=t)
    # noisy = noise(sol.u)
    fig = plot(sol,vars=(1,2,3))
    display(fig)
    return sol
end

#===============================================================#

function denoise!(sol)

end

#===============================================================#

function munge(A::Vector{Vector{Float64}})
    l1 = length(A)
    l2 = length(A[1])
    M = zeros(l2,l1)
    for  i = 1:l1
        M[:,i] = A[i]
    end
    M
end

#===============================================================#

abstract type AbstractDiff  end

struct AnalyticalDeriv <: AbstractDiff end
struct FiniteDiff <: AbstractDiff end
struct TotalVariationalDerivativative <: AbstractDiff end

function differentiate(sol,o::TotalVariationalDerivativative)
    # interp = LagrangeInterpolation(munge(sol.u),sol.t)
    interp = QuadraticInterpolation(munge(sol.u),sol.t)
    y_reconstruct = interp.(sol.t)
    v = [DataInterpolations.derivative(interp,k) for k in sol.t]
    temp = munge(v)
    display(plot(temp[1,:],temp[2,:],temp[3,:]))
    return v
end

function differentiate(sol,o::FiniteDiff)
    q = sol.u[2:end] .- sol.u[1:end-1]
    q = q./(sol.t[2] - sol.t[1])
    w = munge(q)
    fig = plot(w[1,:],w[2,:],w[3,:])
    display(fig)
    return w
end

function differentiate(sol,ds::LorenzSystem,o::AnalyticalDeriv)
    function lorenz(u) # Modify when changing parameters
        u1 = 10*(u[2] - u[1])
        u2 = u[1]*(12 - u[3]) - u[2]
        u3 = u[1]*u[2] - 8/3 * u[3]
        return [u1 u2 u3]
    end
    deriv = []
    for i = 1:length(sol.u)
        append!(deriv,lorenz(sol.u[i]))
    end
    deriv = reshape(deriv,(3,length(sol.u)))
    fig = plot(deriv[1,:],deriv[2,:],deriv[3,:])
    display(fig)
    return deriv
end

#===============================================================#
function n2_terms(p::Int)
    return Int(p*(p+1)/2)
end

function quadratic(X::Array)
    """
        Assume data is passed in the form x,y,z
        Solve this issue for an array of arbitrary size.
    """
    s = size(X) # 3,1001
    p = s[1]
    s = s[2]
    n = n2_terms(p) # 6
    C = []
    for i=1:p
        for j=i:p
            append!(C,X[i,:].*X[j,:])
        end
    end
    reshape(C,(n,s))
end

function cubic(X::Array)
    """
        Assume data is passed in the form x,y,z
        Solve this issue for an array of arbitrary size.
    """
    s = size(X) # 3,1001
    p = s[1]
    s = s[2]
    n = n2_terms(p) # 6
    C = []
    for i=1:p
        for j=i:p
            for k = j:p
                append!(C,X[i,:].*X[j,:] .* X[k,:])
            end
        end
    end
    # return C
    reshape(C,(10,s))
end

function linear(X::Array)
    """
        Assuming that the data passed is linear.
    """
    return X
end

function basis(X)
    s = size(X)
    ntsteps = s[2]
    nparams = s[1]
    # Since we consider unit, linear and quadratic terms.
    n = 1 + 3 + n2_terms(nparams) + 10 # You are better than this!
    θ = zeros(n,ntsteps)
    θ[1,:] = ones(ntsteps)
    θ[2:4,:] = linear(X)
    θ[5:10,:] = quadratic(X)
    θ[11:end,:] = cubic(X)
    return θ
end

# X = rand(3,1000)
# cubic(X)
# basis(X)
#===============================================================#
abstract type AbstractOptimizer end

struct LSTSQ <: AbstractOptimizer
    abstol::Float64
    LSTSQ() = new(1e-16)
end

struct STLSQ <: AbstractOptimizer
    abstol::Float64
    λ::Float64
    STLSQ(thres::Float64) = new(1e-16,thres)
end

##TODO : SR3

function _optimize(θ::Matrix{T},v::Vector{Vector{T}},opt::LSTSQ) where T
    θ = permutedims(θ,(2,1))
    v = permutedims(munge(v))
    return θ\v
end

function _optimize(θ::Matrix{T},v::Matrix{T},opt::LSTSQ) where T
    θ = permutedims(θ,(2,1))
    v = permutedims(v)
    return θ\v
end

function _optimize(θ::Matrix{T},v::Matrix{T},opt::STLSQ) where T
    iter = 0
    abstol = 1e-8
    θ = permutedims(θ,(2,1))
    v = permutedims(v)
    maxiter = maximum(collect(size(θ)))
    convlimit = abstol # Is the limit on the difference betweem two iterations.
    conv = Inf

    λ = opt.λ
    ξ = θ\v
    smallnums = abs.(ξ) .< λ
    bignums = @. !smallnums[:,1]
    x = similar(ξ)

    for i = 1:maxiter
        smallnums .= abs.(ξ) .< λ
        bignums .= @. !smallnums[:,1]
        y = similar(ξ) # Iteration level least square estimate.
        ξ[smallnums] .= 0
        # θ[:,bignums] * ξ[bignums,i] = v[:,i]
        for i = 1:size(v,2)
            bignums .= @. !smallnums[:,i]
            ξ[bignums,i] .= θ[:,bignums] \ v[:,i]
        end
        conv = LinearAlgebra.norm2(y - ξ)
        display(conv)
    end
    ξ[smallnums] .= 0
    return ξ
end

function _optimize(θ::Matrix{T},v::Vector{Vector{T}},opt::STLSQ) where T
    _optimize(θ,munge(v),opt)
end

#===============================================================#

struct Basis
    unit(x) = 1
    linear(x) = x
    square(x) = x^2
    q_mixed(x,y) = x*y
    cubic(x) = x^3
    c_mix1(x,y) = x
end

function _remake(ξ,ds::LorenzSystem)

    function make(du,u,p,t)

    end
end

#===============================================================#

# Main function

# Trail 1 : TVD Derivative functions
function trail1()
    obj_lorenz = LorenzSystem()
    sol = datagen(obj_lorenz)
    denoise!(sol)
    diff = TotalVariationalDerivativative()
    v = differentiate(sol,diff)
    X = munge(sol.u)
    θ = basis(X)
    opt = LSTSQ()
    ξ = _optimize(θ,v,opt)
end

# Trail 2 : Finite Diff Derivative
function trail2()
    obj_lorenz = LorenzSystem()
    sol = datagen(obj_lorenz)
    denoise!(sol)
    diff = FiniteDiff()
    w = differentiate(sol,diff)
    X = munge(sol.u)[:,2:1001]
    θ = basis(X)
    opt = LSTSQ()
    ξ = _optimize(θ,w,opt)
end

# Trail 3 : AnalyticalDeriv and the STLSQ optimizer.
function trail3()
    obj_lorenz = LorenzSystem()
    sol = datagen(obj_lorenz)
    denoise!(sol)
    # diff = AnalyticalDeriv()
    # d = Matrix{Float64}(differentiate(sol,obj_lorenz,diff))
    diff = TotalVariationalDerivativative()
    d = differentiate(sol,diff)
    X = munge(sol.u)
    θ = basis(X)
    display(LinearAlgebra.cond(θ))
    opt = STLSQ(0.05)
    ξ = _optimize(θ,d,opt)
    display(ξ)
    return θ,ξ
end

# Conclusion : Optimizer works!!
θ,ξ = trail3()

ξ

LinearAlgebra.norm2(θ)

#TODO: Reimplement the construction of the library function
#   - Test based implmentation of the library function with unit tests along every step of the way.

#TODO:  Automate the remake function
#   - List of basis  functions
#   - Every velocity vector is a linear combination of all these functions taken together.
#   - Use lambda functions
#   - This needs to be brute forced for the moment. Think of a better approach

#TODO: Implement the sparse regressor SR3
#   - This is a wish.
#   - The LSTSQ might be failing, but it is hard to determine the need for a new regressor without fixing the problem of the condition number.

# LinearAlgebra.cond(θ)
# function lremade!(du,u,p,t)
#     du[1] = -10*u[1] + 10*u[2]
#     du[2] = 15.53 - 15.86*u[1] + 11.26*u[2] + 0.84*u[3]
#     du[3] = -15.63 + 6.99*u[1] + 3.63*u[2] - 3.799*u[3]
# end
#
# p = [10,12,8/3]
# u₀ = [1,0,0]
# tspan = (0,100)
# t = 0:0.01:100
# remade =  ODEProblem(lremade!,u₀,tspan,p)
# sol1 = solve(remade,Tsit5(),saveat=t)
# plot(sol1,vars=(1,2,3))
# savefig("./examples/Learned.png")
#
# obj = LorenzSystem()
# sol = datagen(obj)
# savefig("./examples/Original.png")


text:
# globalization.jl
# this file contains the methods for applying differnent globalizations
# to Newtons method

#------------------------------------------------------------------------------
# Inexact-Newton-Krylov
#------------------------------------------------------------------------------
@doc """
### NonlinearSolvers.updateKrylov

  This function up dates the relative toleranced for the linear solve when doing
  inexact Newton-Krylov.

  Inputs/Outputs:
    newton_data: the NewtonData object

"""->
function updateKrylov(newton_data::NewtonData)

  # this function provides the opportunity to update the krylov tolerances.
  # Currently we don't change the tolerances, because leaving the relative
  # tolerance in place seems to work ok
  if newton_data.use_inexact_nk
    norm_i = newton_data.res_norm_i
    norm_i_1 = newton_data.res_norm_i_1
    gamma = newton_data.krylov_gamma
    #reltol = newton_data.ls.reltol*(norm_i/norm_i_1)^gamma
    reltol = newton_data.ls.reltol
    setTolerances(newton_data.ls, reltol, -1, -1, -1)
    #println(BSTDOUT, "updating krylov reltol to ", reltol)
  end

  return nothing
end

#------------------------------------------------------------------------------
# Psuedo-Transient Continuation (aka. Implicit Euler)
#------------------------------------------------------------------------------

"""
  Reinitailize the underlying data.  This function should be called every
  time newtonInner() is entered
"""
function clearEulerConstants(ls::LinearSolver)

  lo = getInnerLO(ls.lo, NewtonLinearObject)
  clearEulerConstants(lo)

  if !(typeof(ls.pc) <: PCNone)
    pc = getInnerPC(ls.pc, NewtonLinearObject)
    clearEulerConstants(pc)
  end

  return nothing
end

function clearEulerConstants(lo::NewtonLinearObject)

  lo.idata.res_norm_i = 0
  lo.idata.res_norm_i_1 = 0

  return nothing
end

"""
  Record the most recent residual norm

  **Inputs**

   * ls: a [`LinearSolver`](@ref) object
   * res_norm_i: must be a real number
"""
function recordEulerResidual(ls::LinearSolver, res_norm_i)

  lo = getInnerLO(ls.lo, NewtonLinearObject)
  lo.idata.res_norm_i = res_norm_i

  if !(typeof(ls.pc) <: PCNone)
    pc = getInnerPC(ls.pc, NewtonLinearObject)
    pc.idata.res_norm_i = res_norm_i
  end
end

"""
  Get the most recent residual norm and the residual norm the last time
  the PC was recalculated.  This function should be used to inspect the
  values only.  [`useEulerConstants`](@ref) should be used when updating the
  diagonal term in the Jacobian.

  **Inputs**

   * pc or lo: a PC or LO object

  **Outputs**

   * res_norm_i_1: second most recent residual norm
   * res_norm_i: most recent residual norm
"""
function getEulerConstants(pc::AbstractPC)

  pc2 = getInnerPC(pc, NewtonLinearObject)
  return pc2.idata.res_norm_i_1, pc2.idata.res_norm_i
end

function getEulerConstants(lo::AbstractLO)

  lo = getInnerLO(lo, NewtonLinearObject)
  return lo.idata.res_norm_i_1, lo.idata.res_norm_i
end

"""
  Similar to [`getEulerConstants`](@ref), this function both returns the
  Euler constants and shifts them so that the next relative residual will
  be correct.
"""
function useEulerConstants(pc::AbstractPC)
  pc2 = getInnerPC(pc, NewtonLinearObject)
  t1 = pc2.idata.res_norm_i_1
  t2 = pc2.idata.res_norm_i
  pc2.idata.res_norm_i_1 = t2
  return t1, t2
end

function useEulerConstants(lo::AbstractLO)
  lo2 = getInnerLO(lo, NewtonLinearObject)
  t1 = lo2.idata.res_norm_i_1
  t2 = lo2.idata.res_norm_i
  lo2.idata.res_norm_i_1 = t2
  return t1, t2
end

"""
  Check if Euler globalization is initialized
"""
function isEulerInitialized(lo::AbstractLO)
  lo2 = getInnerLO(lo, NewtonLinearObject)
  return !(lo2.idata.res_norm_i == 0 && lo2.idata.res_norm_i_1 == 0) &&
          lo2.idata.use_implicit_euler
end

function isEulerInitialized(pc::AbstractPC)
  pc2 = getInnerPC(pc, NewtonLinearObject)
  return !(pc2.idata.res_norm_i == 0 && pc2.idata.res_norm_i_1 == 0) &&
           pc2.idata.use_implicit_euler
end

function isEulerInitialized(pc::PCNone)
  return false  # can't globalize what doesn't exist
end

@doc """
### NonlinearSolver.calcTauVec

  This function calculates the spatially varying vector for the timestep.

  The timestep varies according to tau/(1 + |J|^(1/d) ), where `|J|` is the
  determinant of the mapping jacobian and `d` is the dimensionality (2D ord
  3D).

  **Inputs**

  *  mesh
  *  opts
  *  tau:  timestep factor

  **Inputs/Outputs**

   * tau_vec: vector (of length numDof) populated with tau*spatially varying
             factor

  Aliasing restrictions:  none
"""->
function calcTauVec(mesh, opts, tau, tau_vec)
# calculate the spatially varying pseudo-timestep
  #TODO: make tau_vec = 1/tau_vec, so we don't have to do fp division when
  #      applying it
  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      for k=1:mesh.numDofPerNode
        dof = mesh.dofs[k, j, i]
        tau_vec[dof] = tau/(1 + sqrt(real(mesh.jac[j, i])))  #TODO: should be |J|^1/dim
#        tau_vec[dof] = tau
      end
    end
  end

  if mesh.myrank == 0
    println(BSTDOUT, "average tau value = ", mean(tau_vec))
  end

  return nothing

end

"""
  Reinitilize an `ImplicitEulerData` object (for example, when calling
  `newtonInner` a second time).

  Uses opts["euler_tau"] to reset the pseudo-time step

  **Inputs**

   * mesh
   * opts
   * idata: `ImplicitEulerData`
"""
function reinitImplicitEuler(mesh, opts, idata::ImplicitEulerData)

  idata.res_norm_i = 0
  idata.res_norm_i_1 = 0
  idata.tau_l = opts["euler_tau"]
  calcTauVec(mesh, opts, idata.tau_l, idata.tau_vec)

  return nothing
end

@doc """
### NonlinearSolver.updateEuler

  This function updates the tau paraemter for the implicit Euler
  globalization.  tau_vec is also updated

  Inputs/Outputs:
    lo: any kind of Newton PC or LO object

"""->
function updateEuler(lo::NewtonLinearObject)
  # updates the tau parameter for the Implicit Euler globalization
  # norm_i is the residual step norm, norm_i_1 is the previous residual norm

  if !lo.idata.use_implicit_euler
    return nothing
  end

  tau_l_old = lo.idata.tau_l
  res_norm_i_1, res_norm_i = useEulerConstants(lo)

  # on the first Jacobian calculate, don't update
  if res_norm_i_1 == 0
    return nothing
  end

  tau_update = res_norm_i_1/res_norm_i

  # update tau
  lo.idata.tau_l = lo.idata.tau_l * tau_update

  if lo.myrank == 0
    println(BSTDOUT, "tau_update factor = ", tau_update)
  end

  scale!(lo.idata.tau_vec, tau_update)

  return nothing
end

@doc """
### NonlinearSolvers.applyEuler

  This function updates the jacobian matrix with the term from the implicit
  Euler globalization.  The term is eqn.M/tau_vec.  Methods are available for
  dense, sparse, Petsc jacobians, as well as jacobian-vector products.

 ** Inputs**
    mesh
    sbp
    eqn
    opts

  **Inputs/Outputs**
    lo:  a [`NewtonHasMat`](@ref) t contaiing tau_vec and the Jacobian matrix

"""->
function applyEuler(mesh, sbp, eqn, opts, lo::NewtonHasMat)
# this allocations memory every time
# should there be a reusable array for this?
# maybe something in lo?
# for explicitly stored jacobian only

  if !isEulerInitialized(lo)
    return nothing
  end

  if mesh.myrank == 0
    println(BSTDOUT, "applying Implicit Euler globalization")
    println(BSTDOUT, "average tau value = ", mean(lo.idata.tau_vec))
  end

  lo2 = getBaseObject(lo)
#  println("euler globalization tau = ", lo.tau_l)
  # create the indices

  val = [1/lo.idata.tau_l]
  idx = PetscInt[0]
  idy = PetscInt[0]
  #TODO: replace this with MatDiagonalSet
  for i=1:mesh.numDof
    idx[1] = i + mesh.dof_offset
    idy[1] = i + mesh.dof_offset
    val[1] = -eqn.M[i]/lo.idata.tau_vec[i]
    set_values1!(lo2.A, idx, idy, val, ADD_VALUES)
#    PetscMatSetValues(jac, idx, idy, val, ADD_VALUES)
  end

  return nothing
end

"""
  Updates a jacobian-vector product with the effects of the globalization

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * vec: vector that the jacobian is being multiplied by
   * lo: a linaer operator or preconditioner, usually a matrix-free one

  **Inputs/Outputs**

   * b: result vector, updated
"""
function applyEuler(mesh, sbp, eqn, opts, vec::AbstractArray,
                    lo::NewtonLinearObject, b::AbstractArray)
# apply the diagonal update term to the jacobian vector product

  if !isEulerInitialized(lo)
    return nothing
  end

  for i=1:mesh.numDof
    b[i] -= eqn.M[i]*(1/lo.idata.tau_vec[i])*vec[i]
  end

  return nothing
end

# Globalize by adding a diagonal term to the Jacobian.  Rather crude.
#------------------------------------------------------------------------------

function addDiagonal(mesh, sbp, eqn, jac)
# add the mass matrix to the jacobian

  for i=1:mesh.numDof
    idx = PetscInt[i-1]
    idy = PetscInt[i-1]
    vals = [100*eqn.M[i]]

#     println("adding ", vals, " to jacobian entry ", i, ",", i)
    PetscMatSetValues(jac, idx, idy, vals, ADD_VALUES)
  end

  return nothing

end

#------------------------------------------------------------------------------


text:
@testset "Gabor" begin
    v1 = rand(3)
    v2 = rand(3)
    ell = abs(rand())
    p = abs(rand())
    k = GaborKernel(; ell=ell, p=p)
    @test k.ell ≈ ell atol = 1e-5
    @test k.p ≈ p atol = 1e-5

    k_manual = exp(-sqeuclidean(v1, v2) / (2 * k.ell^2)) * cospi(euclidean(v1, v2) / k.p)
    @test k(v1, v2) ≈ k_manual atol = 1e-5

    lhs_manual = transform(SqExponentialKernel(), 1 / k.ell)(v1, v2)
    rhs_manual = transform(CosineKernel(), 1 / k.p)(v1, v2)
    @test k(v1, v2) ≈ lhs_manual * rhs_manual atol = 1e-5

    k = GaborKernel()
    @test k.ell ≈ 1.0 atol = 1e-5
    @test k.p ≈ 1.0 atol = 1e-5
    @test repr(k) == "Gabor Kernel (ell = 1.0, p = 1.0)"

    test_interface(k, Vector{Float64})

    test_ADs(x -> GaborKernel(; ell=x[1], p=x[2]), [ell, p]; ADs=[:Zygote])

    # Tests are also failing randomly for ForwardDiff and ReverseDiff but randomly
end


text:
"""
Mirrors C struct `mongoc_query_flags_t`.

These flags correspond to the MongoDB wire protocol.
They may be bitwise or’d together. They may modify how a query is performed in the MongoDB server.

From: http://mongoc.org/libmongoc/current/mongoc_query_flags_t.html
"""
primitive type QueryFlags sizeof(Cint) * 8 end

Base.convert(::Type{T}, t::QueryFlags) where {T<:Number} = reinterpret(Cint, t)
Base.convert(::Type{QueryFlags}, n::T) where {T<:Number} = reinterpret(QueryFlags, n)
QueryFlags(u::Cint) = convert(QueryFlags, u)
QueryFlags(i::Number) = QueryFlags(Cint(i))
Cint(flags::QueryFlags) = convert(Cint, flags)
Base.:(|)(flag1::QueryFlags, flag2::QueryFlags) = QueryFlags( Cint(flag1) | Cint(flag2) )
Base.:(&)(flag1::QueryFlags, flag2::QueryFlags) = QueryFlags( Cint(flag1) & Cint(flag2) )
Base.show(io::IO, flags::QueryFlags) = print(io, "QueryFlags($(Cint(flags)))")

"Specify no query flags."
const QUERY_FLAG_NONE              = QueryFlags(0)

"Cursor will not be closed when the last data is retrieved. You can resume this cursor later."
const QUERY_FLAG_TAILABLE_CURSOR   = QueryFlags(1 << 1)

"Allow query of replica set secondaries."
const QUERY_FLAG_SLAVE_OK          = QueryFlags(1 << 2)

"Used internally by MongoDB."
const QUERY_FLAG_OPLOG_REPLAY      = QueryFlags(1 << 3)

"The server normally times out an idle cursor after an inactivity period (10 minutes). This prevents that."
const QUERY_FLAG_NO_CURSOR_TIMEOUT = QueryFlags(1 << 4)

"Use with MONGOC_QUERY_TAILABLE_CURSOR. Block rather than returning no data. After a period, time out."
const QUERY_FLAG_AWAIT_DATA        = QueryFlags(1 << 5)

"Stream the data down full blast in multiple “reply” packets. Faster when you are pulling down a lot of data and you know you want to retrieve it all."
const QUERY_FLAG_EXHAUST           = QueryFlags(1 << 6)

"Get partial results from mongos if some shards are down (instead of throwing an error)."
const QUERY_FLAG_PARTIAL           = QueryFlags(1 << 7)

"`URI` is a wrapper for C struct `mongoc_uri_t`."
mutable struct URI
    uri::String
    handle::Ptr{Cvoid}

    function URI(uri_string::String)
        err = BSONError()
        handle = mongoc_uri_new_with_error(uri_string, err)
        if handle == C_NULL
            error("Failed to parse URI $uri_string. Error Message: $(err)")
        end
        new_uri = new(uri_string, handle)
        finalizer(destroy!, new_uri)
        return new_uri
    end
end

"`Client` is a wrapper for C struct `mongoc_client_t`."
mutable struct Client
    uri::String
    handle::Ptr{Cvoid}

    function Client(uri::URI)
        client_handle = mongoc_client_new_from_uri(uri.handle)
        if client_handle == C_NULL
            error("Failed connecting to URI $uri.")
        end
        client = new(uri.uri, client_handle)
        finalizer(destroy!, client)
        return client
    end
end

abstract type AbstractDatabase end
abstract type AbstractCollection end

"`Database` is a wrapper for C struct `mongoc_database_t`."
mutable struct Database <: AbstractDatabase
    client::Client
    name::String
    handle::Ptr{Cvoid}

    function Database(client::Client, db_name::String)
        db = new(client, db_name, mongoc_client_get_database(client.handle, db_name))
        finalizer(destroy!, db)
        return db
    end
end

"`Collection` is a wrapper for C struct `mongoc_collection_t`."
mutable struct Collection <: AbstractCollection
    database::Database
    name::String
    handle::Ptr{Cvoid}

    function Collection(database::Database, collection_name::String)
        collection_handle = mongoc_database_get_collection(database.handle, collection_name)
        if collection_handle == C_NULL
            error("Failed creating collection $collection_name on db $(database.name).")
        end
        collection = new(database, collection_name, collection_handle)
        finalizer(destroy!, collection)
        return collection
    end
end

const CursorSource = Union{Client, AbstractDatabase, AbstractCollection}

"`Cursor` is a wrapper for C struct `mongoc_cursor_t`."
mutable struct Cursor{T<:CursorSource}
    source::T
    handle::Ptr{Cvoid}

    function Cursor(source::T, handle::Ptr{Cvoid}) where {T<:CursorSource}
        cursor = new{T}(source, handle)
        finalizer(destroy!, cursor)
        return cursor
    end
end

mutable struct BulkOperation
    collection::Collection
    handle::Ptr{Cvoid}
    executed::Bool

    function BulkOperation(collection::Collection; options::Union{Nothing, BSON}=nothing)
        options_handle = options == nothing ? C_NULL : options.handle
        handle = mongoc_collection_create_bulk_operation_with_opts(collection.handle, options_handle)
        if handle == C_NULL
            error("Failed to create a new bulk operation.")
        end
        bulk_operation = new(collection, handle, false)
        finalizer(destroy!, bulk_operation)
        return bulk_operation
    end
end

struct InsertOneResult
    reply::BSON
    inserted_oid::Union{Nothing, String}
end

struct BulkOperationResult
    reply::BSON
    server_id::UInt32
    inserted_oids::Vector{Union{Nothing, String}}
end

#
# Basic functions for types
#

function destroy!(uri::URI)
    if uri.handle != C_NULL
        mongoc_uri_destroy(uri.handle)
        uri.handle = C_NULL
    end
    nothing
end

function destroy!(client::Client)
    if client.handle != C_NULL
        mongoc_client_destroy(client.handle)
        client.handle = C_NULL
    end
    nothing
end

function destroy!(database::Database)
    if database.handle != C_NULL
        mongoc_database_destroy(database.handle)
        database.handle = C_NULL
    end
    nothing
end

function destroy!(collection::Collection)
    if collection.handle != C_NULL
        mongoc_collection_destroy(collection.handle)
        collection.handle = C_NULL
    end
    nothing
end

function destroy!(cursor::Cursor)
    if cursor.handle != C_NULL
        mongoc_cursor_destroy(cursor.handle)
        cursor.handle = C_NULL
    end
    nothing
end

function destroy!(bulk_operation::BulkOperation)
    if bulk_operation.handle != C_NULL
        mongoc_bulk_operation_destroy(bulk_operation.handle)
        bulk_operation.handle = C_NULL
        bulk_operation.executed = true
    end
    nothing
end

#
# Session Types
#

mutable struct SessionOptions
    handle::Ptr{Cvoid}

    function SessionOptions(; casual_consistency::Bool=true)
        session_options_handle = mongoc_session_opts_new()
        if session_options_handle == C_NULL
            error("Couldn't create SessionOptions.")
        end

        session_options = new(session_options_handle)
        finalizer(destroy!, session_options)
        set_casual_consistency!(session_options, casual_consistency)
        return session_options
    end
end

mutable struct Session
    client::Client
    options::SessionOptions
    handle::Ptr{Cvoid}

    function Session(client::Client; options::SessionOptions=SessionOptions())
        err = BSONError()
        session_handle = mongoc_client_start_session(client.handle, options.handle, err)
        if session_handle == C_NULL
            error("$err")
        end
        session = new(client, options, session_handle)
        finalizer(destroy!, session)
        return session
    end
end

struct DatabaseSession <: AbstractDatabase
    database::Database
    session::Session
end

struct CollectionSession <: AbstractCollection
    database_session::DatabaseSession
    collection::Collection
end

#=
struct BulkOperationSession
    collection_session::CollectionSession
    bulk_operation::BulkOperation

    function BulkOperationSession(collection_session::CollectionSession, options::Union{Nothing, BSON})
        options_with_session = _join(options, get_session(collection_session))
        bulk_operation = BulkOperation(collection_session.collection, options=options_with_session)
        return BulkOperationSession(collection_session, bulk_operation)
    end
end
=#

function destroy!(session_options::SessionOptions)
    if session_options.handle != C_NULL
        mongoc_session_opts_destroy(session_options.handle)
        session_options.handle = C_NULL
    end
    nothing
end

function destroy!(session::Session)
    if session.handle != C_NULL
        mongoc_client_session_destroy(session.handle)
        session.handle = C_NULL
    end
    nothing
end


text:
using GLMakie
function plot_chain(t, x, FULL=false; skip=0, color=:black, title="")
    r = eachindex(t)[1:1+skip:end]
    fig3 = Figure(resolution=(2000,2000), title=title)
    e = 6
    for i in 1:e^2
        u = CartesianIndices((e,e))[i]
        if u[1] == u[2]
            FULL && GLMakie.lines(fig3[u[1],u[2]], t, getindex.(x, u[1]),  color=color)
            !FULL && GLMakie.scatter(fig3[u[1],u[2]], t[r], getindex.(x, u[1])[r], markersize=0.5,  color=color)
        elseif u[1] < u[2]
            FULL && GLMakie.lines(fig3[u[1],u[2]], getindex.(x, u[1]),  getindex.(x, u[2]), linewidth=0.5, color=(color,1.0))
            !FULL && GLMakie.scatter(fig3[u[1],u[2]], getindex.(x, u[1])[r],  getindex.(x, u[2])[r], linewidth=0.5, color=(color,1.0), markersize= 0.5, strokewidth=0)
        end
    end
    fig3
end
function plot_chain!(fig3, t, x, FULL=false; color=:black, skip=0)
    r = eachindex(t)[1:1+skip:end]
    e = 6
    for i in 1:e^2
        u = CartesianIndices((e,e))[i]
        if u[1] == u[2]
            FULL && GLMakie.lines!(fig3[u[1],u[2]], t, getindex.(x, u[1]), color=color)
            !FULL && GLMakie.scatter!(fig3[u[1],u[2]], t[r], getindex.(x, u[1])[r], markersize=0.5, color=color)
        elseif u[1] < u[2]
            FULL && GLMakie.lines!(fig3[u[1],u[2]], getindex.(x, u[1]),  getindex.(x, u[2]), linewidth=0.5, color=(color,1.0))
            !FULL && GLMakie.scatter!(fig3[u[1],u[2]], getindex.(x, u[1])[r],  getindex.(x, u[2])[r], linewidth=0.5, color=(color,1.0), markersize= 0.5, strokewidth=0)
        end
    end
    fig3
end


text:
# Automatically generated using Clang.jl

using ....VectorEngine: VEPtr         # added to autogenerated code

const veo_args = Cvoid
const VEDAargs = Ptr{veo_args}
const veo_ptr = UInt64
const VEDAdevice = Cint
const VEDAidx = UInt32
const VEDAoffset = UInt64
const VEDAstream = Cint
const VEDAfunction = UInt64
const VEDAhost_function = Ptr{Cvoid}
const VEDAstream_callback = Ptr{Cvoid}
const VEDAdeviceptr = VEPtr{Int8}      # modified from autogenerated code
const VEDAcontext = Ptr{__VEDAcontext}
const VEDAmodule = Ptr{__VEDAmodule}

@cenum VEDAresult_enum::UInt32 begin
    VEDA_SUCCESS = 0
    VEDA_ERROR_ALREADY_INITIALIZED = 1
    VEDA_ERROR_CANNOT_CREATE_CONTEXT = 2
    VEDA_ERROR_CANNOT_CREATE_DEVICE = 3
    VEDA_ERROR_CANNOT_CREATE_STREAM = 4
    VEDA_ERROR_CANNOT_POP_CONTEXT = 5
    VEDA_ERROR_CANT_READ_SENSOR_FILE = 6
    VEDA_ERROR_FUNCTION_NOT_FOUND = 7
    VEDA_ERROR_INITIALIZING_DEVICE = 8
    VEDA_ERROR_INVALID_ARGS = 9
    VEDA_ERROR_INVALID_CONTEXT = 10
    VEDA_ERROR_INVALID_COREIDX = 11
    VEDA_ERROR_INVALID_DEVICE = 12
    VEDA_ERROR_INVALID_MODULE = 13
    VEDA_ERROR_INVALID_REQID = 14
    VEDA_ERROR_INVALID_STREAM = 15
    VEDA_ERROR_INVALID_VALUE = 16
    VEDA_ERROR_MODULE_NOT_FOUND = 17
    VEDA_ERROR_NOT_IMPLEMENTED = 18
    VEDA_ERROR_NOT_INITIALIZED = 19
    VEDA_ERROR_NO_DEVICES_FOUND = 20
    VEDA_ERROR_NO_SENSOR_FILE = 21
    VEDA_ERROR_OFFSETTED_VPTR_NOT_ALLOWED = 22
    VEDA_ERROR_OUT_OF_BOUNDS = 23
    VEDA_ERROR_OUT_OF_MEMORY = 24
    VEDA_ERROR_SHUTTING_DOWN = 25
    VEDA_ERROR_UNABLE_TO_DETECT_DEVICES = 26
    VEDA_ERROR_UNKNOWN_CONTEXT = 27
    VEDA_ERROR_UNKNOWN_DEVICE = 28
    VEDA_ERROR_UNKNOWN_KERNEL = 29
    VEDA_ERROR_UNKNOWN_MODULE = 30
    VEDA_ERROR_UNKNOWN_PPTR = 31
    VEDA_ERROR_UNKNOWN_ATTRIBUTE = 32
    VEDA_ERROR_UNKNOWN_STREAM = 33
    VEDA_ERROR_UNKNOWN_VPTR = 34
    VEDA_ERROR_VEDA_LD_LIBRARY_PATH_NOT_DEFINED = 35
    VEDA_ERROR_VEO_COMMAND_ERROR = 36
    VEDA_ERROR_VEO_COMMAND_EXCEPTION = 37
    VEDA_ERROR_VEO_COMMAND_UNFINISHED = 38
    VEDA_ERROR_VEO_COMMAND_UNKNOWN_ERROR = 39
    VEDA_ERROR_VEO_STATE_BLOCKED = 40
    VEDA_ERROR_VEO_STATE_RUNNING = 41
    VEDA_ERROR_VEO_STATE_SYSCALL = 42
    VEDA_ERROR_VEO_STATE_UNKNOWN = 43
    VEDA_ERROR_VPTR_ALREADY_ALLOCATED = 44
    VEDA_ERROR_TODO = 45
    VEDA_ERROR_ASSERT = 46
    VEDA_ERROR_CONTEXT_IS_DESTROYED = 47
    VEDA_ERROR_UNKNOWN = 48
end

@cenum VEDAdevice_attribute_enum::UInt32 begin
    VEDA_DEVICE_ATTRIBUTE_CLOCK_RATE = 0
    VEDA_DEVICE_ATTRIBUTE_CLOCK_BASE = 1
    VEDA_DEVICE_ATTRIBUTE_MULTIPROCESSOR_COUNT = 2
    VEDA_DEVICE_ATTRIBUTE_MEMORY_CLOCK_RATE = 3
    VEDA_DEVICE_ATTRIBUTE_L1D_CACHE_SIZE = 4
    VEDA_DEVICE_ATTRIBUTE_L1I_CACHE_SIZE = 5
    VEDA_DEVICE_ATTRIBUTE_L2_CACHE_SIZE = 6
    VEDA_DEVICE_ATTRIBUTE_LLC_CACHE_SIZE = 7
    VEDA_DEVICE_ATTRIBUTE_ABI_VERSION = 8
    VEDA_DEVICE_ATTRIBUTE_SINGLE_TO_DOUBLE_PRECISION_PERF_RATIO = 9
    VEDA_DEVICE_ATTRIBUTE_FIREWARE_VERSION = 10
    VEDA_DEVICE_ATTRIBUTE_MODEL = 11
end

@cenum VEDAargs_intent_enum::UInt32 begin
    VEDA_ARGS_INTENT_IN = 0
    VEDA_ARGS_INTENT_INOUT = 1
    VEDA_ARGS_INTENT_OUT = 2
end

@cenum VEDAcontext_mode_enum::UInt32 begin
    VEDA_CONTEXT_MODE_OMP = 0
    VEDA_CONTEXT_MODE_SCALAR = 1
end

const VEDAresult = VEDAresult_enum
const VEDAdevice_attribute = VEDAdevice_attribute_enum
const VEDAargs_intent = VEDAargs_intent_enum
const VEDAcontext_mode = VEDAcontext_mode_enum


text:
using Test
using MPI
using ClimateMachine.Mesh.Topologies

function main()
    MPI.Init()

    comm = MPI.COMM_WORLD
    crank = MPI.Comm_rank(comm)
    csize = MPI.Comm_size(comm)

    @assert csize == 2

    topology = StackedBrickTopology(
        comm,
        (1:4, 5:8, 9:12),
        periodicity = (false, true, false),
        boundary = ((1, 2), (3, 4), (5, 6)),
        connectivity = :face,
    )

    elems = topology.elems
    realelems = topology.realelems
    ghostelems = topology.ghostelems
    sendelems = topology.sendelems
    elemtocoord = topology.elemtocoord
    elemtoelem = topology.elemtoelem
    elemtoface = topology.elemtoface
    elemtoordr = topology.elemtoordr
    elemtobndy = topology.elemtobndy
    nabrtorank = topology.nabrtorank
    nabrtorecv = topology.nabrtorecv
    nabrtosend = topology.nabrtosend

    globalelemtoface = [
        1 1 1 2 2 2 2 2 2 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
        1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2
        4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
        3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
        5 6 6 5 6 6 5 6 6 5 6 6 5 6 6 5 6 6 5 6 6 5 6 6 5 6 6
        5 5 6 5 5 6 5 5 6 5 5 6 5 5 6 5 5 6 5 5 6 5 5 6 5 5 6
    ]

    globalelemtoordr = ones(Int, size(globalelemtoface))

    globalelemtocoord = Array{Int}(undef, 3, 8, 27)

    globalelemtocoord[:, :, 1] =
        [1 2 1 2 1 2 1 2; 5 5 6 6 5 5 6 6; 9 9 9 9 10 10 10 10]
    globalelemtocoord[:, :, 2] =
        [1 2 1 2 1 2 1 2; 5 5 6 6 5 5 6 6; 10 10 10 10 11 11 11 11]
    globalelemtocoord[:, :, 3] =
        [1 2 1 2 1 2 1 2; 5 5 6 6 5 5 6 6; 11 11 11 11 12 12 12 12]
    globalelemtocoord[:, :, 4] =
        [2 3 2 3 2 3 2 3; 5 5 6 6 5 5 6 6; 9 9 9 9 10 10 10 10]
    globalelemtocoord[:, :, 5] =
        [2 3 2 3 2 3 2 3; 5 5 6 6 5 5 6 6; 10 10 10 10 11 11 11 11]
    globalelemtocoord[:, :, 6] =
        [2 3 2 3 2 3 2 3; 5 5 6 6 5 5 6 6; 11 11 11 11 12 12 12 12]
    globalelemtocoord[:, :, 7] =
        [2 3 2 3 2 3 2 3; 6 6 7 7 6 6 7 7; 9 9 9 9 10 10 10 10]
    globalelemtocoord[:, :, 8] =
        [2 3 2 3 2 3 2 3; 6 6 7 7 6 6 7 7; 10 10 10 10 11 11 11 11]
    globalelemtocoord[:, :, 9] =
        [2 3 2 3 2 3 2 3; 6 6 7 7 6 6 7 7; 11 11 11 11 12 12 12 12]
    globalelemtocoord[:, :, 10] =
        [1 2 1 2 1 2 1 2; 6 6 7 7 6 6 7 7; 9 9 9 9 10 10 10 10]
    globalelemtocoord[:, :, 11] =
        [1 2 1 2 1 2 1 2; 6 6 7 7 6 6 7 7; 10 10 10 10 11 11 11 11]
    globalelemtocoord[:, :, 12] =
        [1 2 1 2 1 2 1 2; 6 6 7 7 6 6 7 7; 11 11 11 11 12 12 12 12]
    globalelemtocoord[:, :, 13] =
        [1 2 1 2 1 2 1 2; 7 7 8 8 7 7 8 8; 9 9 9 9 10 10 10 10]
    globalelemtocoord[:, :, 14] =
        [1 2 1 2 1 2 1 2; 7 7 8 8 7 7 8 8; 10 10 10 10 11 11 11 11]
    globalelemtocoord[:, :, 15] =
        [1 2 1 2 1 2 1 2; 7 7 8 8 7 7 8 8; 11 11 11 11 12 12 12 12]
    globalelemtocoord[:, :, 16] =
        [2 3 2 3 2 3 2 3; 7 7 8 8 7 7 8 8; 9 9 9 9 10 10 10 10]
    globalelemtocoord[:, :, 17] =
        [2 3 2 3 2 3 2 3; 7 7 8 8 7 7 8 8; 10 10 10 10 11 11 11 11]
    globalelemtocoord[:, :, 18] =
        [2 3 2 3 2 3 2 3; 7 7 8 8 7 7 8 8; 11 11 11 11 12 12 12 12]
    globalelemtocoord[:, :, 19] =
        [3 4 3 4 3 4 3 4; 7 7 8 8 7 7 8 8; 9 9 9 9 10 10 10 10]
    globalelemtocoord[:, :, 20] =
        [3 4 3 4 3 4 3 4; 7 7 8 8 7 7 8 8; 10 10 10 10 11 11 11 11]
    globalelemtocoord[:, :, 21] =
        [3 4 3 4 3 4 3 4; 7 7 8 8 7 7 8 8; 11 11 11 11 12 12 12 12]
    globalelemtocoord[:, :, 22] =
        [3 4 3 4 3 4 3 4; 6 6 7 7 6 6 7 7; 9 9 9 9 10 10 10 10]
    globalelemtocoord[:, :, 23] =
        [3 4 3 4 3 4 3 4; 6 6 7 7 6 6 7 7; 10 10 10 10 11 11 11 11]
    globalelemtocoord[:, :, 24] =
        [3 4 3 4 3 4 3 4; 6 6 7 7 6 6 7 7; 11 11 11 11 12 12 12 12]
    globalelemtocoord[:, :, 25] =
        [3 4 3 4 3 4 3 4; 5 5 6 6 5 5 6 6; 9 9 9 9 10 10 10 10]
    globalelemtocoord[:, :, 26] =
        [3 4 3 4 3 4 3 4; 5 5 6 6 5 5 6 6; 10 10 10 10 11 11 11 11]
    globalelemtocoord[:, :, 27] =
        [3 4 3 4 3 4 3 4; 5 5 6 6 5 5 6 6; 11 11 11 11 12 12 12 12]

    globalelemtobndy = [
        1 1 1 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 2 2 2 2 2 2 2
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
        5 0 0 5 0 0 5 0 0 5 0 0 5 0 0 5 0 0 5 0 0 5 0 0 5 0 0
        0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6 0 0 6
    ]

    if crank == 0
        nrealelem = 12
        globalelems = [
            1,
            2,
            3, # 1
            4,
            5,
            6, # 2
            7,
            8,
            9, # 3
            10,
            11,
            12, # 4
            13,
            14,
            15, # 5
            16,
            17,
            18, # 6
            22,
            23,
            24, # 8
            25,
            26,
            27,
        ] # 9

        elemtoelem_expect = [
            1 2 3 1 2 3 10 11 12 4 5 6 7 8 9 16 17 18 19 20 21 22 23 24
            4 5 6 22 23 24 19 20 21 7 8 9 13 14 15 16 17 18 1 2 3 4 5 6
            13 14 15 16 17 18 4 5 6 1 2 3 13 14 15 16 17 18 19 20 21 22 23 24
            10 11 12 7 8 9 16 17 18 13 14 15 13 14 15 16 17 18 19 20 21 22 23 24
            1 1 2 2 4 5 3 7 8 4 10 11 5 14 15 6 17 18 7 20 21 8 23 24
            2 3 1 5 6 2 8 9 3 11 12 4 13 14 5 16 17 6 19 20 7 22 23 8
        ]

        nabrtorank_expect = [1]
        nabrtorecv_expect = UnitRange{Int}[1:12]
        nabrtosend_expect = UnitRange{Int}[1:12]
    elseif crank == 1
        nrealelem = 15

        globalelems = [
            13,
            14,
            15, # 5
            16,
            17,
            18, # 6
            19,
            20,
            21, # 7
            22,
            23,
            24, # 8
            25,
            26,
            27, # 9
            1,
            2,
            3, # 1
            4,
            5,
            6, # 2
            7,
            8,
            9, # 3
            10,
            11,
            12,
        ] # 4

        elemtoelem_expect = [
            1 2 3 1 2 3 4 5 6 22 23 24 19 20 21 4 5 6 19 20 21 22 23 24 7 8 9
            4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 16 17 18 19 20 21 22 23 24 25 26 27
            25 26 27 22 23 24 10 11 12 13 14 15 7 8 9 16 17 18 19 20 21 22 23 24 25 26 27
            16 17 18 19 20 21 13 14 15 7 8 9 10 11 12 16 17 18 19 20 21 22 23 24 25 26 27
            1 1 2 2 4 5 3 7 8 4 10 11 5 13 14 6 17 18 7 20 21 8 23 24 9 26 27
            2 3 1 5 6 2 8 9 3 11 12 4 14 15 5 16 17 6 19 20 7 22 23 8 25 26 9
        ]

        nabrtorank_expect = [0]
        nabrtorecv_expect = UnitRange{Int}[1:12]
        nabrtosend_expect = UnitRange{Int}[1:12]
    end

    @test elems == 1:length(globalelems)
    @test realelems == 1:nrealelem
    @test ghostelems == (nrealelem + 1):length(globalelems)

    @test elemtocoord == globalelemtocoord[:, :, globalelems]
    @test elemtoface[:, realelems] ==
          globalelemtoface[:, globalelems[realelems]]
    @test elemtoelem == elemtoelem_expect
    @test elemtobndy == globalelemtobndy[:, globalelems]
    @test elemtoordr == ones(eltype(elemtoordr), size(elemtoordr))
    @test nabrtorank == nabrtorank_expect
    @test nabrtorecv == nabrtorecv_expect
    @test nabrtosend == nabrtosend_expect

    @test collect(realelems) ==
          sort(union(topology.exteriorelems, topology.interiorelems))
    @test unique(sort(sendelems)) == topology.exteriorelems
    @test length(intersect(topology.exteriorelems, topology.interiorelems)) == 0
end

main()