text:
const RUN_ON_HOST_IF_NORETURN = Ref(false)

struct FailedInference <: Exception
    f::Any
    kernel_args::Tuple
    kernel_return::Any
    host_args::Tuple
    host_return::Any
end

function Base.showerror(io::IO, err::FailedInference)
    print(io, FailedInference, ": ")
    if err.kernel_return === Union{}
        print(io, "Kernel is inferred to not return (return type is `Union{}`)")
    else
        print(io, "Kernel is inferred to return invalid type: ", err.kernel_return)
    end
    if err.kernel_return === Union{} && err.host_return !== Union{}
        println(io)
        print(io, "Note: on the host, the return type is inferred as ", err.host_return)
    end
    println(io)
    printstyled(io, "HINT"; bold = true, color = :light_black)
    printstyled(
        io,
        ": if this exception is caught as `err``, use `CUDA.code_typed(err)` to",
        " introspect the erronous code.";
        color = :light_black,
    )
end

CUDA.code_typed(err::FailedInference; options...) =
    CUDA.code_typed(err.f, typeof(err.kernel_args); options...)


text:
@testset "Quoted Symbol Syntax" begin
    @test (:a → :b) == (StenoGraphs.@quote_symbols a → b)
    b = :c
    @test (:a → :c) == (StenoGraphs.@quote_symbols a → _(b))
    let err = nothing
        try
            eval(:(StenoGraphs.@quote_symbols a → _(b, c)))
        catch err
        end

        @test err isa Exception
        @test  occursin("Unqote only a single argument. Right: `_(x)`, Wrong: `_(x, y)`.", sprint(showerror, err))
    end
end

@testset "Addition as hcat" begin
    @test StenoGraphs.addition_to_vector(:(a * b + c)) == :([a * b c])
    @test StenoGraphs.addition_to_vector(:(a * (b + c))) == :(a * [b c])
end

@testset "Broadcasting Edges" begin
    @test Edge(:a, [:b :c]) == vec([Edge(:a, :b) Edge(:a, :c)])
    @test Edge([:a :b], :c) == vec([Edge(:a, :c) Edge(:b, :c)])
    @test Edge([:a :b], [:c :d]) == vec([Edge(:a, :c) Edge(:a, :d) Edge(:b, :c) Edge(:b, :d)])
end

@testset "Singleline Arrows" begin
    @test [Edge(:a, :c), Edge(:b, :c)] == @StenoGraph [a b] → c
end

@testset "Multiline Arrows" begin
    @test [Edge(:a, :c), Edge(:b, :c), Edge(:f, :e)] == @StenoGraph begin
        [a b] → c
        e ← f
    end
end

@testset "Multiline Mix Arrow Edge" begin
    @test [Edge(:a, :c), Edge(:b, :c), Edge(:f, :e)] ==
    @StenoGraph begin
        [a b] → c
        Edge(f, e)
    end
end

ext: jl


text:
replace_post(ex, s, s′) = postwalk(x -> x == s ? s′ : x, ex)

function LinearRecEntry(expr::Expr)
    CoeffT = QQ
    _funcs = Symbol[]
    _args = RExpr[]
    function_walk(expr) do f, x
        @assert length(x) == 1
        push!(_funcs, f)
        push!(_args, x[1])
    end
    _funcs = Base.unique(_funcs)
    @assert length(_funcs) > 0 "Not a recurrence: no functions present"
    argsyms = Base.unique(Iterators.flatten(symbols(a) for a in _args))
    @assert length(argsyms) == 1 "More (or less) than one variable in the function arguments, got: $(argsyms)"

    _argsym = argsyms[1]
    @assert !(_argsym in _funcs) "Ambiguous symbol: $(_argsym) is a function and an argument"

    R, argsym = PolynomialRing(CoeffT, string(_argsym))

    s_expr = SymEngine.Basic(expr)
    s_argsym = SymEngine.Basic(_argsym)
    s_args = map(x->(SymEngine.Basic(x) - s_argsym), _args)
    # minarg = minimum(convert(Rational, x) for x in s_args)
    minarg = convert(Int, minimum(s_args))
    if !iszero(minarg)
        _expr = replace_post(expr, _argsym, :($_argsym - $minarg))
        return LinearRecEntry(_expr)
    end
    maxarg = convert(Int, maximum(s_args))

    hom = zero(SymEngine.Basic)
    dicts = [Dict{Symbol,FracElem}() for _ in 1:maxarg + 1]
    for i in 0:maxarg
        for f in _funcs
            fc = SymEngine.SymFunction(f)(s_argsym + i)
            co = SymEngine.coeff(s_expr, fc)
            if !iszero(co)
                hom += co * fc
                # call R twice to deal with integers
                num = numerator(co)
                den = denominator(co)
                cnum = R(R(convert(Expr, num)))
                cden = R(R(convert(Expr, den)))
                dicts[i+1][f] = cnum // cden
            end
        end
    end

    inhom = R(R(convert(Expr, SymEngine.expand(-(s_expr - hom)))))//R(1)

    (coeffs = dicts, inhom = inhom), argsym
end

macro lrs(input)
    @capture(input, begin fields__ end)
    lrs(Vector{Expr}(fields))
end

function entries(ex::Expr...)
    _entries = []
    _args = []
    _syms = Symbol[]
    for x in ex
        if @capture(x, l_ = r_)
            x = :($(unblock(l)) - $(unblock(r)))
            @assert @capture(l, v_(a_))
            push!(_syms, v)
        end
        entry, arg = LinearRecEntry(x)
        push!(_entries, entry)
        push!(_args, arg)
    end
    _args = Base.unique(_args)
    @assert length(_args) == 1
    _syms, _entries, _args[1]
end

function lrs(exprs::Vector{Expr})
    syms, ls, arg = entries(exprs...)
    lrs = LinearRecSystem(arg, syms)
    push!(lrs, ls...)
    lrs
end


text:
## This allows "dekron"ing a PDE operator into two ODE operators
# used in schurfact

# productop is a product of two operator

typealias WrapperOps Union{ConversionWrapper,MultiplicationWrapper,DerivativeWrapper,IntegralWrapper,LaplacianWrapper}

isproductop(a) = iskronop(a)  # all kron ops are product ops

iskronop(::) = false
iskronop(::KroneckerOperator) = true

iskronop(A::Union{WrapperOps,SpaceOperator,ConstantTimesOperator})=iskronop(A.op)
iskronop(A::TimesOperator)=all(iskronop,A.ops)
iskronop{T,DS<:TensorSpace}(::ConstantOperator{T,DS}) = true
iskronop{T,S<:TensorSpace,V<:TensorSpace}(::ZeroOperator{T,S,V}) = true

iskronsumop(::)=false
iskronsumop(A::PlusOperator)=all(a->iskronop(a)||iskronsumop(a),A.ops)
iskronsumop(A::Union{WrapperOps,SpaceOperator,ConstantTimesOperator})=iskronsumop(A.op)

#dekron gives the operators that make up a productop

dekron(K::KroneckerOperator,k)=K.ops[k]

dekron(S::WrapperOps,k)=dekron(S.op,k)
dekron(S::TimesOperator,k)=TimesOperator(Operator{eltype(eltype(S))}[dekron(op,k) for op in S.ops])
dekron(S::SpaceOperator,k)=SpaceOperator(dekron(S.op,k),domainspace(S)[k],rangespace(S)[k])
#TODO: dekron(S::SpaceOperator,k)=SpaceOperator(dekron(S.op,k),domainspace(S)[k],rangespace(S)[k])
dekron(sp::ConstantTimesOperator,k)=k==1?sp.λ*dekron(sp.op,k):dekron(sp.op,k)
dekron{T,DS<:TensorSpace}(C::ConstantOperator{T,DS},k) =
    k==1?ConstantOperator(C.λ,C.space[1]):ConstantOperator(one(C.λ),C.space[k])
dekron{T,S<:TensorSpace,V<:TensorSpace}(C::ZeroOperator{T,S,V},k) =
    ZeroOperator(C.domainspace[k],C.rangespace[k])

dekron(K)=dekron(K,1),dekron(K,2)

sumops(A)=A.ops
sumops(A::Union{SpaceOperator,WrapperOps})=sumops(A.op)
sumops(A::ConstantTimesOperator)=Operator{eltype(A)}[A.λ*op for op in sumops(A.op)]
function dekron(T::Type,A,k,::Colon)
    ret=Array(Operator{T},0)
    if iskronop(A)
        push!(ret,dekron(A,k))
    elseif iskronsumop(A)
        for op in sumops(A)
            da=dekron(op,k,:)
            for a in da
                push!(ret,a)
            end
        end
    else
       error("dekron not defined unless iskronop or iskronsumop" )
    end
    ret
end

dekron(A,k,::Colon)=dekron(eltype(eltype(A)),A,k,:)

dekron(A,::Colon,::Colon)=dekron(A,1,:),dekron(A,2,:)

function reducekronsum!(opsx,opsy)
    @assert length(opsx)==length(opsy)
    k=1
    while k ≤ length(opsx)
        for j=length(opsx):-1:k+1
            if opsx[k]==opsx[j]
               opsy[k]+=opsy[j]
               deleteat!(opsx,j)
               deleteat!(opsy,j)
           end
         end
        k+=1
    end
    opsx,opsy
end

function simplifydekron(A)
    opsx,opsy=dekron(A,:,:)
    map!(simplify,opsx);map!(simplify,opsy)
    reducekronsum!(opsy,opsx)
    reducekronsum!(opsx,opsy)
end


text:
#==============================================================================#
#
#                 some general & simple utility functions
#
#==============================================================================#

"""
    outer(x)

Compute an outer product
"""
outer(x) = x*x'

"""
    get_or_default(container, elem::Symbol, default)

Return `container`.`elem` if it exists, otherwise return `default`
"""
function get_or_default(container, elem::Symbol, default)
    haskey(container, elem) ? getindex(container, elem) : default
end

"""
    static_accessor_HFc(u::SVector, ::Val{T}) where T

Access data stored in the container `u` so that it matches the shapes of H,F,c
and points to the correct points in `u`. `T` is the dimension of the stochastic
process.
"""
function static_accessor_HFc(u::K, ::Val{T}) where {K<:Union{SVector,MVector},T}
    Hidx = SVector{T*T,Int64}(1:T*T)
    Fidx = SVector{T,Int64}((T*T+1):(T*T+T))
    #cidx = SVector{1,Int64}((T*T+T+1):(T*T+T+1))
    reshape(u[Hidx], Size(T,T)), u[Fidx], u[T*T+T+1]
end

@doc raw"""
    standard_guid_prop_time_transf(tt)

Standard time transformation used for guided proposals:
```math
τ: t → t_0 + (t-t_0)\left[2 - \frac{(t-t_0)}{T-t_0}\right],
```
applied to a vector `tt` where $t_0:=$`tt[1]` and $T:=$`tt[end]`
"""
standard_guid_prop_time_transf(tt) = τ(tt[1], tt[end]).(tt)
τ(t₀,T) = (t) ->  t₀ + (t-t₀) * (2-(t-t₀)/(T-t₀))


text:
using PyPlot, ElasticWave
nz = 179; nx = 237; ext= 30;   iflag = 2;
dx = 20. ; dz = 20. ; dt = 2e-3; tmax = 2.0; f0=12.;
vp=4000.*ones(nz,nx);  vs=2000.*ones(nz,nx); rho=2.5e-4*ones(nz,nx);
fidMtx = CreateFidMtx(nz, nx, ext, iflag, vp, vs, rho, dz, dx, dt, f0);
# multi-sources case
isz = 1; isx = 118; ot = 0.0;
flags = vec([false false true true false]);
src = InitSource(isz, isx, nz, nx, ext, iflag, f0, ot, dt, flags);
# =========finite difference modeling==========
root = homedir();
path = join([root, "/Desktop/spt.bin"])
MultiStepForward(path, src, fidMtx, st=3*dt, otype="spt", tmax=tmax)

path1 = join([root, "/Desktop/base.bin"])
lambda = QRbase(path1, path, rk=200, style="mem")
Q = formStateMatrix(path1)
path_Fr = join([root, "/Desktop/Fr.bin"])
lowRankFd(path_Fr, fidMtx, path1, rk=200)
Q = formBaseMatrix(path1, rk)
path_FQ = join([root, "/Desktop/FQ"])
(nz, nx, ext, iflag, dt, nt) = InfoSnapShots(path_FQ)
if iflag == 1
   Nz = nz + 2*ext
elseif iflag == 2
   Nz = nz +   ext
end
Nx = nx + 2*ext
fid = open(path_FQ, "r"); seek(fid, sizeof(Float32)*5)
FQ = reshape(read(fid, Float32, Nz*Nx*10*nt), Nz*Nx*10, nt)
Fr = Q'*FQ
Fr1 = ReadRDFD(path_Fr)
vecnorm(Fr1-Fr)

# ===============test source projection(checked) =================
root = homedir()
path = join([root, "/Desktop/base.bin"])
spj = SrcProj(src, path)
Q = formBaseMatrix(path);
Qt = Q';
Nz = nz+ext; Nx=nx+2*ext; iz=isz; ix=isx+ext;
sc = zeros(size(Qt,1), src.nt);
for it = 1 : src.nt
    ind = (ix-1)*Nz+iz + Nz*Nx*4
    sc[:,it] = sc[:,it] + Qt[:,ind]*src.Txx[it]/2
    ind = (ix-1)*Nz+iz + Nz*Nx*5
    sc[:,it] = sc[:,it] + Qt[:,ind]*src.Txx[it]/2
    ind = (ix-1)*Nz+iz + Nz*Nx*6
    sc[:,it] = sc[:,it] + Qt[:,ind]*src.Tzz[it]/2
    ind = (ix-1)*Nz+iz + Nz*Nx*7
    sc[:,it] = sc[:,it] + Qt[:,ind]*src.Tzz[it]/2
end
vecnorm(spj.sc-sc)

# ========= test more MOR =======================
using ElasticWave
nz = 179; nx = 237; ext= 30;   iflag = 1;
dx = 20. ; dz = 20. ; dt = 2e-3; tmax = 2.0; f0=12.;
vp=4000.*ones(nz,nx);  vs=2000.*ones(nz,nx); rho=2.5e-4*ones(nz,nx);
fidMtx = CreateFidMtx(nz, nx, ext, iflag, vp, vs, rho, dz, dx, dt, f0);
# multi-sources case
isz = 1; isx = 118; ot = 0.0;
flags = vec([true false false false false]);
src = InitSource(isz, isx, nz, nx, ext, iflag, f0, ot, dt, flags);
# =========finite difference modeling==========
root = homedir();
path = join([root, "/Desktop/spt.bin"])
MultiStepForward(path, src, fidMtx, st=dt, otype="spt", tmax=tmax)

rk = 200
pathe = join([root, "/Desktop/base.bin"])
lambda = QRbase(pathe, path, rk=rk, style="mem")

pathFr = join([root, "/Desktop/Fr.bin"])
lowRankFd(pathFr, fidMtx, pathe;)

spj = SrcProj(src, pathe, rk=rk)
pathc = join([root, "/Desktop/wc.bin"])
Fr = ReadRDFD(pathFr);

MultiStepForward_MOR(pathc, pathFr, spj, tmax=tmax)
path1 = join([root, "/Desktop/spt1.bin"])
ModSpts2Spts(path1, pathe, pathc, rk=rk)

path1 = join([root, "/Desktop/spt1.bin"])
MultiStepForward(path1, src, fidMtx, st=dt, otype="spt", tmax=tmax)

Nz = nz+ext; Nx = nx+2*ext; len = Nz*Nx;
Vx = tmp[1:len] + tmp[len+1:2*len]
Vx = reshape(Vx, Nz, Nx)

spt = ReadSnapShot(pathe, 50);
vx = reshape(spt.Vxx+spt.Vxz, 239, 297);
SeisPlot(vx)

spt = ReadSnapShot(pathe, 120);
txx = reshape(spt.Txxx+spt.Txxz, 239, 297);
SeisPlot(txx)

# path_spt2vx = join([root, "/Desktop/spt2vx.bin"])
# ExtractCpt(path_spt2vx, path, itype="spt")
# path_wfd2vx = join([root, "/Desktop/wfd2vx.bin"])
# ExtractCpt(path_wfd2vx, path1, itype="wfd")
#
# fid  = open(path_spt2vx, "r")
# fid1 = open(path_wfd2vx, "r")
# tmp = 0.0; tmp1 = 0.0
# nz=read(fid, Int32); nx=read(fid, Int32); ext=read(fid, Int32); iflag=read(fid, Int32); dt=read(fid, Float32);
# sptSize = sizeof(Float32)*nz*nx
# nt = floor(Int64, (filesize(fid)-sizeof(Float32)*5)/sptSize)
#     it = 450
#     pos = (it-1)*sptSize +  sizeof(Float32)*5
#     seek(fid,  pos)
#     seek(fid1, pos)
#     vx = read(fid, Float32, nz*nx)
#     vx1 = read(fid1, Float32, nz*nx)
#     vx = reshape(vx, Int64(nz), Int64(nx))
#     figure()
#     vx1 = reshape(vx1, Int64(nz), Int64(nx))

max_stars_repo_name: WneleiGao/ElasticWave


text:
using Test, MDInclude

@testset "MDInclude" begin
    @mdinclude("tests.md")
end


text:
using xim
using Test

xim_haircut_testsum(1,2)
xim_basics_testsum(3,4)

@testset "xim.jl" begin
    # Write your own tests here.

xim_haircut_testsum(1,2)

#println("hi")

end

max_stars_repo_name: NelsonMCLopes/xim

max_issues_repo_name: NelsonMCLopes/xim

max_forks_repo_name: NelsonMCLopes/xim


text:
# generic operations on associative collections

abstract Associative{K,V}

const secret_table_token = :__c782dbf1cf4d6a2e5e3865d7e95634f2e09b5902__

haskey(d::Associative, k) = in(k,keys(d))

function in(p::(Any,Any), a::Associative)
    v = get(a,p[1],secret_table_token)
    !is(v, secret_table_token) && (v == p[2])
end

function summary(t::Associative)
    n = length(t)
    string(typeof(t), " with ", n, (n==1 ? " entry" : " entries"))
end

function show{K,V}(io::IO, t::Associative{K,V})
    if isempty(t)
        print(io, typeof(t),"()")
    else
        if K === Any && V === Any
            delims = ['{','}']
        else
            delims = ['[',']']
        end
        print(io, delims[1])
        first = true
        for (k, v) in t
            first || print(io, ',')
            first = false
            show(io, k)
            print(io, "=>")
            show(io, v)
        end
        print(io, delims[2])
    end
end

function _truncate_at_width_or_chars(str, width, chars="", truncmark="…")
    truncwidth = strwidth(truncmark)
    (width <= 0 || width < truncwidth) && return ""

    wid = truncidx = lastidx = 0
    idx = start(str)
    while !done(str, idx)
        lastidx = idx
        c, idx = next(str, idx)
        wid += charwidth(c)
        wid >= width - truncwidth && truncidx == 0 && (truncidx = lastidx)
        (wid >= width || c in chars) && break
    end

    lastidx != 0 && str[lastidx] in chars && (lastidx = prevind(str, lastidx))
    truncidx == 0 && (truncidx = lastidx)
    if lastidx < sizeof(str)
        return bytestring(SubString(str, 1, truncidx) * truncmark)
    else
        return bytestring(str)
    end
end

showdict(t::Associative; kw...) = showdict(STDOUT, t; kw...)
function showdict{K,V}(io::IO, t::Associative{K,V}; limit::Bool = false,
                       sz=(s = tty_size(); (s[1]-3, s[2])))
    rows, cols = sz
    print(io, summary(t))
    isempty(t) && return
    print(io, ":")

    if limit
        rows < 2   && (print(io, " …"); return)
        cols < 12  && (cols = 12) # Minimum widths of 2 for key, 4 for value
        cols -= 6 # Subtract the widths of prefix "  " separator " => "
        rows -= 2 # Subtract the summary and final ⋮ continuation lines

        # determine max key width to align the output, caching the strings
        ks = Array(String, min(rows, length(t)))
        keylen = 0
        for (i, k) in enumerate(keys(t))
            i > rows && break
            ks[i] = sprint(show, k)
            keylen = clamp(length(ks[i]), keylen, div(cols, 3))
        end
    end

    for (i, (k, v)) in enumerate(t)
        print(io, "\n  ")
        limit && i > rows && (print(io, rpad("⋮", keylen), " => ⋮"); break)

        if limit
            key = rpad(_truncate_at_width_or_chars(ks[i], keylen, "\r\n"), keylen)
        else
            key = sprint(show, k)
        end
        print(io, key)
        print(io, " => ")

        val = sprint(show, v)
        if limit
            val = _truncate_at_width_or_chars(val, cols - keylen, "\r\n")
        end
        print(io, val)
    end
end

immutable KeyIterator{T<:Associative}
    dict::T
end
immutable ValueIterator{T<:Associative}
    dict::T
end

summary{T<:Union(KeyIterator,ValueIterator)}(iter::T) =
    string(T.name, " for a ", summary(iter.dict))

show(io::IO, iter::Union(KeyIterator,ValueIterator)) = show(io, collect(iter))

showkv(iter::Union(KeyIterator,ValueIterator); kw...) = showkv(STDOUT, iter; kw...)
function showkv{T<:Union(KeyIterator,ValueIterator)}(io::IO, iter::T; limit::Bool = false,
                                                     sz=(s = tty_size(); (s[1]-3, s[2])))
    rows, cols = sz
    print(io, summary(iter))
    isempty(iter) && return
    print(io, ". ", T<:KeyIterator ? "Keys" : "Values", ":")
    if limit
        rows < 2 && (print(io, " …"); return)
        cols < 4 && (cols = 4)
        cols -= 2 # For prefix "  "
        rows -= 2 # For summary and final ⋮ continuation lines
    end

    for (i, v) in enumerate(iter)
        print(io, "\n  ")
        limit && i >= rows && (print(io, "⋮"); break)

        str = sprint(show, v)
        limit && (str = _truncate_at_width_or_chars(str, cols, "\r\n"))
        print(io, str)
    end
end

length(v::Union(KeyIterator,ValueIterator)) = length(v.dict)
isempty(v::Union(KeyIterator,ValueIterator)) = isempty(v.dict)
eltype(v::KeyIterator) = eltype(v.dict)[1]
eltype(v::ValueIterator) = eltype(v.dict)[2]

start(v::Union(KeyIterator,ValueIterator)) = start(v.dict)
done(v::Union(KeyIterator,ValueIterator), state) = done(v.dict, state)

function next(v::KeyIterator, state)
    n = next(v.dict, state)
    n[1][1], n[2]
end

function next(v::ValueIterator, state)
    n = next(v.dict, state)
    n[1][2], n[2]
end

in(k, v::KeyIterator) = !is(get(v.dict, k, secret_table_token),
                            secret_table_token)

keys(a::Associative) = KeyIterator(a)
values(a::Associative) = ValueIterator(a)

function copy(a::Associative)
    b = similar(a)
    for (k,v) in a
        b[k] = v
    end
    return b
end

function merge!(d::Associative, others::Associative...)
    for other in others
        for (k,v) in other
            d[k] = v
        end
    end
    return d
end
merge(d::Associative, others::Associative...) = merge!(copy(d), others...)

function filter!(f::Function, d::Associative)
    for (k,v) in d
        if !f(k,v)
            delete!(d,k)
        end
    end
    return d
end
filter(f::Function, d::Associative) = filter!(f,copy(d))

eltype{K,V}(a::Associative{K,V}) = (K,V)

function isequal(l::Associative, r::Associative)
    if isa(l,ObjectIdDict) != isa(r,ObjectIdDict)
        return false
    end
    if length(l) != length(r) return false end
    for (key, value) in l
        if !isequal(value, get(r, key, secret_table_token))
            return false
        end
    end
    true
end

function ==(l::Associative, r::Associative)
    if isa(l,ObjectIdDict) != isa(r,ObjectIdDict)
        return false
    end
    if length(l) != length(r) return false end
    for (key, value) in l
        if value != get(r, key, secret_table_token)
            return false
        end
    end
    true
end

# some support functions

_tablesz(x::Integer) = x < 16 ? 16 : one(x)<<((sizeof(x)<<3)-leading_zeros(x-1))

function getindex(t::Associative, key)
    v = get(t, key, secret_table_token)
    if is(v, secret_table_token)
        throw(KeyError(key))
    end
    return v
end

# t[k1,k2,ks...] is syntactic sugar for t[(k1,k2,ks...)].  (Note
# that we need to avoid dispatch loops if setindex!(t,v,k) is not defined.)
getindex(t::Associative, k1, k2, ks...) = getindex(t, tuple(k1,k2,ks...))
setindex!(t::Associative, v, k1, k2, ks...) = setindex!(t, v, tuple(k1,k2,ks...))

push!(t::Associative, key, v) = setindex!(t, v, key)

# hashing objects by identity

type ObjectIdDict <: Associative{Any,Any}
    ht::Array{Any,1}
    ObjectIdDict() = new(cell(32))

    function ObjectIdDict(itr)
        d = ObjectIdDict()
        for (k,v) in itr
            d[k] = v
        end
        d
    end

    function ObjectIdDict(o::ObjectIdDict)
        N = length(o.ht)
        ht = cell(N)
        ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, Uint),
              ht, o.ht, N*sizeof(Ptr))
        new(ht)
    end
end

similar(d::ObjectIdDict) = ObjectIdDict()

function setindex!(t::ObjectIdDict, v::ANY, k::ANY)
    t.ht = ccall(:jl_eqtable_put, Array{Any,1}, (Any, Any, Any), t.ht, k, v)
    return t
end

get(t::ObjectIdDict, key::ANY, default::ANY) =
    ccall(:jl_eqtable_get, Any, (Any, Any, Any), t.ht, key, default)

pop!(t::ObjectIdDict, key::ANY, default::ANY) =
    ccall(:jl_eqtable_pop, Any, (Any, Any, Any), t.ht, key, default)

function pop!(t::ObjectIdDict, key::ANY)
    val = pop!(t, key, secret_table_token)
    !is(val,secret_table_token) ? val : throw(KeyError(key))
end

function delete!(t::ObjectIdDict, key::ANY)
    ccall(:jl_eqtable_pop, Any, (Any, Any), t.ht, key)
    t
end

empty!(t::ObjectIdDict) = (t.ht = cell(length(t.ht)); t)

_oidd_nextind(a, i) = int(ccall(:jl_eqtable_nextind, Csize_t, (Any, Csize_t), a, i))

start(t::ObjectIdDict) = _oidd_nextind(t.ht, 0)
done(t::ObjectIdDict, i) = (i == -1)
next(t::ObjectIdDict, i) = ((t.ht[i+1],t.ht[i+2]), _oidd_nextind(t.ht, i+2))

function length(d::ObjectIdDict)
    n = 0
    for pair in d
        n+=1
    end
    n
end

copy(o::ObjectIdDict) = ObjectIdDict(o)

# dict

type Dict{K,V} <: Associative{K,V}
    slots::Array{Uint8,1}
    keys::Array{K,1}
    vals::Array{V,1}
    ndel::Int
    count::Int
    deleter::Function

    function Dict()
        n = 16
        new(zeros(Uint8,n), Array(K,n), Array(V,n), 0, 0, identity)
    end
    function Dict(ks, vs)
        # TODO: eventually replace with a call to Dict(zip(ks,vs))
        n = min(length(ks), length(vs))
        h = Dict{K,V}()
        for i=1:n
            h[ks[i]] = vs[i]
        end
        return h
    end
    function Dict(kv)
        h = Dict{K,V}()
        for (k,v) in kv
            h[k] = v
        end
        return h
    end
end
Dict() = Dict{Any,Any}()

Dict{K,V}(ks::AbstractArray{K}, vs::AbstractArray{V}) = Dict{K,V}(ks,vs)
Dict(ks, vs) = Dict{Any,Any}(ks, vs)

# conversion between Dict types
function convert{K,V}(::Type{Dict{K,V}},d::Dict)
    h = Dict{K,V}()
    for (k,v) in d
        ck = convert(K,k)
        if !haskey(h,ck)
            h[ck] = convert(V,v)
        else
            error("key collision during dictionary conversion")
        end
    end
    return h
end
convert{K,V}(::Type{Dict{K,V}},d::Dict{K,V}) = d

# syntax entry points
Dict{K,V}(ks::(K...), vs::(V...)) = Dict{K  ,V  }(ks, vs)
Dict{K  }(ks::(K...), vs::Tuple ) = Dict{K  ,Any}(ks, vs)
Dict{V  }(ks::Tuple , vs::(V...)) = Dict{Any,V  }(ks, vs)

Dict{K,V}(kv::AbstractArray{(K,V)}) = Dict{K,V}(kv)
Dict{K,V}(kv::Associative{K,V}) = Dict{K,V}(kv)

similar{K,V}(d::Dict{K,V}) = (K=>V)[]

function serialize(s, t::Dict)
    serialize_type(s, typeof(t))
    write(s, int32(length(t)))
    for (k,v) in t
        serialize(s, k)
        serialize(s, v)
    end
end

function deserialize{K,V}(s, T::Type{Dict{K,V}})
    n = read(s, Int32)
    t = T(); sizehint(t, n)
    for i = 1:n
        k = deserialize(s)
        v = deserialize(s)
        t[k] = v
    end
    return t
end

hashindex(key, sz) = (int(hash(key)) & (sz-1)) + 1

isslotempty(h::Dict, i::Int) = h.slots[i] == 0x0
isslotfilled(h::Dict, i::Int) = h.slots[i] == 0x1
isslotmissing(h::Dict, i::Int) = h.slots[i] == 0x2

function rehash{K,V}(h::Dict{K,V}, newsz)
    olds = h.slots
    oldk = h.keys
    oldv = h.vals
    sz = length(olds)
    newsz = _tablesz(newsz)
    if h.count == 0
        resize!(h.slots, newsz)
        fill!(h.slots, 0)
        resize!(h.keys, newsz)
        resize!(h.vals, newsz)
        h.ndel = 0
        return h
    end

    slots = zeros(Uint8,newsz)
    keys = Array(K, newsz)
    vals = Array(V, newsz)
    count0 = h.count
    count = 0

    for i = 1:sz
        if olds[i] == 0x1
            k = oldk[i]
            v = oldv[i]
            index = hashindex(k, newsz)
            while slots[index] != 0
                index = (index & (newsz-1)) + 1
            end
            slots[index] = 0x1
            keys[index] = k
            vals[index] = v
            count += 1

            if h.count != count0
                # if items are removed by finalizers, retry
                return rehash(h, newsz)
            end
        end
    end

    h.slots = slots
    h.keys = keys
    h.vals = vals
    h.count = count
    h.ndel = 0

    return h
end

function sizehint(d::Dict, newsz)
    oldsz = length(d.slots)
    if newsz <= oldsz
        # todo: shrink
        # be careful: rehash() assumes everything fits. it was only designed
        # for growing.
        return d
    end
    # grow at least 25%
    newsz = max(newsz, (oldsz*5)>>2)
    rehash(d, newsz)
end

function empty!{K,V}(h::Dict{K,V})
    fill!(h.slots, 0x0)
    sz = length(h.slots)
    h.keys = Array(K, sz)
    h.vals = Array(V, sz)
    h.ndel = 0
    h.count = 0
    return h
end

# get the index where a key is stored, or -1 if not present
function ht_keyindex{K,V}(h::Dict{K,V}, key)
    sz = length(h.keys)
    iter = 0
    maxprobe = max(16, sz>>6)
    index = hashindex(key, sz)
    keys = h.keys

    while true
        if isslotempty(h,index)
            break
        end
        if !isslotmissing(h,index) && isequal(key,keys[index])
            return index
        end

        index = (index & (sz-1)) + 1
        iter+=1
        iter > maxprobe && break
    end

    return -1
end

# get the index where a key is stored, or -pos if not present
# and the key would be inserted at pos
# This version is for use by setindex! and get!
function ht_keyindex2{K,V}(h::Dict{K,V}, key)
    sz = length(h.keys)
    iter = 0
    maxprobe = max(16, sz>>6)
    index = hashindex(key, sz)
    avail = 0
    keys = h.keys

    while true
        if isslotempty(h,index)
            avail < 0 && return avail
            return -index
        end

        if isslotmissing(h,index)
            if avail == 0
                # found an available slot, but need to keep scanning
                # in case "key" already exists in a later collided slot.
                avail = -index
            end
        elseif isequal(key, keys[index])
            return index
        end

        index = (index & (sz-1)) + 1
        iter+=1
        iter > maxprobe && break
    end

    avail < 0 && return avail

    rehash(h, h.count > 64000 ? sz*2 : sz*4)

    return ht_keyindex2(h, key)
end

function _setindex!(h::Dict, v, key, index)
    h.slots[index] = 0x1
    h.keys[index] = key
    h.vals[index] = v
    h.count += 1

    sz = length(h.keys)
    # Rehash now if necessary
    if h.ndel >= ((3*sz)>>2) || h.count*3 > sz*2
        # > 3/4 deleted or > 2/3 full
        rehash(h, h.count > 64000 ? h.count*2 : h.count*4)
    end
end

function setindex!{K,V}(h::Dict{K,V}, v0, key0)
    key = convert(K,key0)
    if !isequal(key,key0)
        error(key0, " is not a valid key for type ", K)
    end
    v = convert(V,  v0)

    index = ht_keyindex2(h, key)

    if index > 0
        h.keys[index] = key
        h.vals[index] = v
    else
        _setindex!(h, v, key, -index)
    end

    return h
end

function get!{K,V}(h::Dict{K,V}, key0, default)
    key = convert(K,key0)
    if !isequal(key,key0)
        error(key0, " is not a valid key for type ", K)
    end

    index = ht_keyindex2(h, key)

    index > 0 && return h.vals[index]

    v = convert(V,  default)
    _setindex!(h, v, key, -index)
    return v
end

function get!{K,V}(default::Function, h::Dict{K,V}, key0)
    key = convert(K,key0)
    if !isequal(key,key0)
        error(key0, " is not a valid key for type ", K)
    end

    index = ht_keyindex2(h, key)

    index > 0 && return h.vals[index]

    v = convert(V,  default())
    _setindex!(h, v, key, -index)
    return v
end

# NOTE: this macro is specific to Dict, not Associative, and should
#       therefore not be exported as-is: it's for internal use only.
macro get!(h, key0, default)
    quote
        K, V = eltype($(esc(h)))
        key = convert(K, $(esc(key0)))
        isequal(key, $(esc(key0))) || error($(esc(key0)), " is not a valid key for type ", K)
        idx = ht_keyindex2($(esc(h)), key)
        if idx < 0
            idx = -idx
            v = convert(V, $(esc(default)))
            _setindex!($(esc(h)), v, key, idx)
        else
            @inbounds v = $(esc(h)).vals[idx]
        end
        v
    end
end

function getindex{K,V}(h::Dict{K,V}, key)
    index = ht_keyindex(h, key)
    return (index<0) ? throw(KeyError(key)) : h.vals[index]::V
end

function get{K,V}(h::Dict{K,V}, key, deflt)
    index = ht_keyindex(h, key)
    return (index<0) ? deflt : h.vals[index]::V
end

function get{K,V}(deflt::Function, h::Dict{K,V}, key)
    index = ht_keyindex(h, key)
    return (index<0) ? deflt() : h.vals[index]::V
end

haskey(h::Dict, key) = (ht_keyindex(h, key) >= 0)
in{T<:Dict}(key, v::KeyIterator{T}) = (ht_keyindex(v.dict, key) >= 0)

function getkey{K,V}(h::Dict{K,V}, key, deflt)
    index = ht_keyindex(h, key)
    return (index<0) ? deflt : h.keys[index]::K
end

function _pop!(h::Dict, index)
    val = h.vals[index]
    _delete!(h, index)
    return val
end

function pop!(h::Dict, key)
    index = ht_keyindex(h, key)
    index > 0 ? _pop!(h, index) : throw(KeyError(key))
end

function pop!(h::Dict, key, default)
    index = ht_keyindex(h, key)
    index > 0 ? _pop!(h, index) : default
end

function _delete!(h::Dict, index)
    h.slots[index] = 0x2
    ccall(:jl_arrayunset, Void, (Any, Uint), h.keys, index-1)
    ccall(:jl_arrayunset, Void, (Any, Uint), h.vals, index-1)
    h.ndel += 1
    h.count -= 1
    h
end

function delete!(h::Dict, key)
    index = ht_keyindex(h, key)
    if index > 0; _delete!(h, index); end
    h
end

function skip_deleted(h::Dict, i)
    L = length(h.slots)
    while i<=L && !isslotfilled(h,i)
        i += 1
    end
    return i
end

start(t::Dict) = skip_deleted(t, 1)
done(t::Dict, i) = done(t.vals, i)
next(t::Dict, i) = ((t.keys[i],t.vals[i]), skip_deleted(t,i+1))

isempty(t::Dict) = (t.count == 0)
length(t::Dict) = t.count

next{T<:Dict}(v::KeyIterator{T}, i) = (v.dict.keys[i], skip_deleted(v.dict,i+1))
next{T<:Dict}(v::ValueIterator{T}, i) = (v.dict.vals[i], skip_deleted(v.dict,i+1))

# weak key dictionaries

function weak_key_delete!(t::Dict, k)
    # when a weak key is finalized, remove from dictionary if it is still there
    wk = getkey(t, k, secret_table_token)
    if !is(wk,secret_table_token) && is(wk.value, k)
        delete!(t, k)
    end
end

function add_weak_key(t::Dict, k, v)
    if is(t.deleter, identity)
        t.deleter = x->weak_key_delete!(t, x)
    end
    t[WeakRef(k)] = v
    # TODO: it might be better to avoid the finalizer, allow
    # wiped WeakRefs to remain in the table, and delete them as
    # they are discovered by getindex and setindex!.
    finalizer(k, t.deleter)
    return t
end

function weak_value_delete!(t::Dict, k, v)
    # when a weak value is finalized, remove from dictionary if it is still there
    wv = get(t, k, secret_table_token)
    if !is(wv,secret_table_token) && is(wv.value, v)
        delete!(t, k)
    end
end

function add_weak_value(t::Dict, k, v)
    t[k] = WeakRef(v)
    finalizer(v, x->weak_value_delete!(t, k, x))
    return t
end

type WeakKeyDict{K,V} <: Associative{K,V}
    ht::Dict{Any,V}

    WeakKeyDict() = new((Any=>V)[])
end
WeakKeyDict() = WeakKeyDict{Any,Any}()

setindex!{K}(wkh::WeakKeyDict{K}, v, key) = add_weak_key(wkh.ht, convert(K,key), v)

function getkey{K}(wkh::WeakKeyDict{K}, kk, deflt)
    k = getkey(wkh.ht, kk, secret_table_token)
    if is(k, secret_table_token)
        return deflt
    end
    return k.value::K
end

get{K}(wkh::WeakKeyDict{K}, key, def) = get(wkh.ht, key, def)
get{K}(def::Function, wkh::WeakKeyDict{K}, key) = get(def, wkh.ht, key)
get!{K}(wkh::WeakKeyDict{K}, key, def) = get!(wkh.ht, key, def)
get!{K}(def::Function, wkh::WeakKeyDict{K}, key) = get!(def, wkh.ht, key)
pop!{K}(wkh::WeakKeyDict{K}, key) = pop!(wkh.ht, key)
pop!{K}(wkh::WeakKeyDict{K}, key, def) = pop!(wkh.ht, key, def)
delete!{K}(wkh::WeakKeyDict{K}, key) = delete!(wkh.ht, key)
empty!(wkh::WeakKeyDict)  = (empty!(wkh.ht); wkh)
haskey{K}(wkh::WeakKeyDict{K}, key) = haskey(wkh.ht, key)
getindex{K}(wkh::WeakKeyDict{K}, key) = getindex(wkh.ht, key)
isempty(wkh::WeakKeyDict) = isempty(wkh.ht)

start(t::WeakKeyDict) = start(t.ht)
done(t::WeakKeyDict, i) = done(t.ht, i)
function next{K}(t::WeakKeyDict{K}, i)
    kv, i = next(t.ht, i)
    ((kv[1].value::K,kv[2]), i)
end
length(t::WeakKeyDict) = length(t.ht)

max_stars_repo_name: jey/julia

max_issues_repo_name: Aerlinger/julia

max_forks_repo_name: Aerlinger/julia


text:
using Mathy
using Test

include("../docs/utils.jl")

@testset "Mathy.jl" begin
    @test (mathy_makedocs(); true)
end


text:
module IOstrings

export IOstrings_outpath_name

using Dates

include("IOstrings_outpath_name.jl")

end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


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

"""
    AbstractDataProvider

The root type for all data provider. A data provider should implement the following interfaces:

* [`get_batch_size`](@ref)
* [`provide_data`](@ref)
* [`provide_label`](@ref)

As well as the Julia iterator interface (see
[the Julia manual](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-iteration-1)).
Normally this involves defining:

* `Base.eltype(provider) -> AbstractDataBatch`
* `Base.iterate(provider[, state]) -> (AbstractDataBatch, AbstractDataProvider)`
"""
abstract type AbstractDataProvider end

"""
    get_batch_size(provider) -> Int

# Arguments:
* `provider::AbstractDataProvider`: the data provider.

Returns the mini-batch size of the provided data. All the provided data should have the same mini-batch size (i.e. the last dimension).
"""
get_batch_size

"""
    provide_data(provider) -> Vector{Tuple{Base.Symbol, Tuple}}

# Arguments:
* `provider::AbstractDataProvider`: the data provider.

Returns a vector of (name, shape) pairs describing the names of the data it provides, and the corresponding shapes.

"""
provide_data

"""
    provide_label(provider) -> Vector{Tuple{Base.Symbol, Tuple}}

# Arguments:
* `provider::AbstractDataProvider`: the data provider.

Returns a vector of (name, shape) pairs describing the names of the labels it provides, and the corresponding shapes.
"""
provide_label

"""
    AbstractDataProviderState

   Base type for data provider states.
"""
abstract type AbstractDataProviderState end

"""
    AbstractDataBatch

   Base type for a data mini-batch. It should implement the following interfaces:

* [`count_samples`](@ref)
* [`get_data`](@ref)
* [`get_label`](@ref)

The following utility functions will be automatically defined:

* [`get`](@ref)
* [`load_data!`](@ref)
* [`load_label!`](@ref)
"""
abstract type AbstractDataBatch end

"""
    count_samples(provider, batch) -> Int

# Arguments:
* `batch::AbstractDataBatch`: the data batch object.

Returns the number of samples in this batch. This number should be greater than 0, but less than or equal to the batch size. This is used to indicate at the end of the data set, there might not be enough samples for a whole mini-batch.

"""
count_samples

"""
    get_data(provider, batch) -> Vector{NDArray}

# Arguments:
* `provider::AbstractDataProvider`: the data provider.
* `batch::AbstractDataBatch`: the data batch object.

Returns a vector of data in this batch, should be in the same order as declared in `provide_data() <AbstractDataProvider.provide_data>`.

The last dimension of each `NDArray` should always match the batch_size, even when `count_samples` returns a value less than the batch size. In this case,      the data provider is free to pad the remaining contents with any value.
"""
get_data

"""
    get_label(provider, batch) -> Vector{NDArray}

# Arguments:
* `provider::AbstractDataProvider`: the data provider.
* `batch::AbstractDataBatch`: the data batch object.

Returns a vector of labels in this batch. Similar to [`get_data`](@ref).
"""
get_label

"""
    DataBatch

A basic subclass of `AbstractDataBatch`, that implement the interface by
accessing member fields.
"""
mutable struct DataBatch{T,S,N,M} <: AbstractDataBatch
  data  :: Vector{NDArray{T,N}}
  label :: Vector{NDArray{S,M}}
  count :: Int
end

count_samples(batch::DataBatch) = batch.count

get_data(::Provider, batch::DataBatch) where {Provider<:AbstractDataProvider} =
  batch.data

get_label(::Provider, batch::DataBatch) where {Provider<:AbstractDataProvider} =
  batch.label

"""
    SlicedNDArray

A alias type of `Tuple{UnitRange{Int},NDArray}`.
"""
const SlicedNDArray = Tuple{UnitRange{Int},<:NDArray}

function _load_general!(provider :: AbstractDataProvider, batch :: AbstractDataBatch,
                        targets :: Vector{<:Vector{<:SlicedNDArray}}, loader::Function)
  data = loader(provider, batch)
  for (d_src, d_targets) in zip(data, targets)
    for (slice_idx, d_dst) in d_targets
      copy!(d_dst, slice(d_src, slice_idx))
    end
  end
end

"""
    load_data!(provider, batch, targets)

# Arguments:
* `provider::AbstractDataProvider`: the data provider.
* `batch::AbstractDataBatch`: the data batch object.
* `targets::Vector{Vector{SlicedNDArray}}`: the targets to load data into.

The targets is a list of the same length as number of data provided by this provider.
Each element in the list is a list of `SlicedNDArray`. This list described a
spliting scheme of this data batch into different slices, each slice is specified by
a slice-ndarray pair, where *slice* specify the range of samples in the mini-batch
that should be loaded into the corresponding *ndarray*.

This utility function is used in data parallelization, where a mini-batch is splited
and computed on several different devices.
"""
function load_data!(provider :: AbstractDataProvider, batch :: AbstractDataBatch,
                    targets :: Vector{<:Vector{<:SlicedNDArray}})
  _load_general!(provider, batch, targets, get_data)
end

"""
    load_label!(provider, batch, targets)

* `provider::AbstractDataProvider provider`: the data provider.
* `batch::AbstractDataBatch batch`: the data batch object.
* `targets::Vector{Vector{SlicedNDArray}}`: the targets to load label into.

The same as [`load_data!`](@ref), except that this is for loading labels.
"""
function load_label!(provider :: AbstractDataProvider, batch :: AbstractDataBatch,
                     targets :: Vector{<:Vector{<:SlicedNDArray}})
  _load_general!(provider, batch, targets, get_label)
end

function load_data!(provider :: AbstractDataProvider, batch :: AbstractDataBatch,
                    targets :: Vector{<:NDArray})
  for (src, dst) in zip(get_data(provider, batch), targets)
    copy!(dst, src)
  end
end
function load_label!(provider :: AbstractDataProvider, batch :: AbstractDataBatch,
                     targets :: Vector{<:NDArray})
  for (src, dst) in zip(get_label(provider, batch), targets)
    copy!(dst, src)
  end
end

import Base.get
"""
    get(provider, batch, name) -> NDArray

* `provider::AbstractDataProvider`: the data provider.
* `batch::AbstractDataBatch`: the data batch object.
* `name::Symbol`: the name of the data to get, should be one of the names
  provided in either `provide_data() <AbstractDataProvider.provide_data>`
  or `provide_label() <AbstractDataProvider.provide_label>`.

Returns the corresponding data array corresponding to that name.
"""
function get(provider::AbstractDataProvider, batch::AbstractDataBatch, name::Symbol)
  for (idx, (k, s)) in enumerate(provide_data(provider))
    if name == k
      return get_data(provider, batch)[idx]
    end
  end
  for (idx, (k, s)) in enumerate(provide_label(provider))
    if name == k
      return get_label(provider, batch)[idx]
    end
  end
  error("$name is not provided by this data provider")
end

"""
    eachbatch(provider::AbstractDataProvider)

Allows you to perform operations on data every epoch. This is especially useful
when you need to perform real-time augmentation of the data.

# Arguments:
* `provider`: an instance of the custom DataProvider type. You must return this
instance after modifying its fields.

"""
eachbatch(provider::AbstractDataProvider) = provider

"""
    ArrayDataProvider

A convenient tool to iterate `NDArray` or Julia `Array`.

    ArrayDataProvider(data[, label]; batch_size, shuffle, data_padding, label_padding)

Construct a data provider from `NDArray` or Julia Arrays.

# Arguments:
* `data`: the data, could be
  * a `NDArray`, or a Julia Array. This is equivalent to `:data => data`.
  * a name-data pair, like `:mydata => array`, where `:mydata` is the name of the data
  * and `array` is an `NDArray` or a Julia Array.
  * a list of name-data pairs.

* `label`: the same as the `data` parameter. When this argument is omitted, the constructed provider will provide no labels.
* `batch_size::Int`: the batch size, default is 0, which means treating the whole array as a single mini-batch.
* `shuffle::Bool`: turn on if the data should be shuffled at every epoch.
* `data_padding::Real`: when the mini-batch goes beyond the dataset boundary, there might
  be less samples to include than a mini-batch. This value specify a scalar to pad the
  contents of all the missing data points.
* `label_padding::Real`: the same as `data_padding`, except for the labels.

TODO: remove `data_padding` and `label_padding`, and implement rollover that copies
the last or first several training samples to feed the padding.
"""
mutable struct ArrayDataProvider{T,N} <: AbstractDataProvider
  data_arrays   :: Vector{Array{T,N}}
  data_names    :: Vector{Symbol}
  label_arrays
  label_names   :: Vector{Symbol}
  batch_size    :: Int
  sample_count  :: Int
  shuffle       :: Bool
  data_padding  :: MX_float
  label_padding :: MX_float

  data_batch
  label_batch
end

# Julia's type system is sometimes very frustrating. You cannot specify a function
# with argument Vector{Pair} to expect to be matched when calling with the parameter
# [:foo => zeros(2,3), :bar => zeros(3)] because the type inference gives very specific
# results, about the parametric type in the Pair{T1,T2} type, thus does not match the
# generic Pair type. In general, Int <: Number but Vector{Int} <: Vector{Number} is not
# true. So let us just use Any here...
function ArrayDataProvider(data; batch_size::Int = 0, shuffle::Bool = false,
                           data_padding::Real = 0, label_padding::Real = 0)
  ArrayDataProvider(data, [], batch_size = batch_size, shuffle = shuffle,
                    data_padding = data_padding, label_padding = label_padding)
end

function ArrayDataProvider(data, label; batch_size::Int = 0, shuffle::Bool = false,
                           data_padding::Real = 0, label_padding::Real = 0)
  asarr(arr :: Array{T}) where {T} = convert(Array{MX_float}, arr)
  asarr(arr :: NDArray) = copy(arr)

  if isa(data, Union{NDArray, Array}) && eltype(data) <: Real
    data_names  = [:data]
    data_arrays = Array{MX_float}[asarr(data)]
  elseif isa(data, Pair)
    @assert isa(data.first, Base.Symbol) && isa(data.second, Union{NDArray, Array})
    data_names  = [data.first]
    data_arrays = Array{MX_float}[asarr(data.second)]
  elseif isa(data, Vector) || isa(data, Tuple)
    map(data) do d
      @assert isa(d, Pair) && isa(d.first, Base.Symbol) && isa(d.second, Union{NDArray, Array})
    end
    data_names  = Base.Symbol[d.first for d in data]
    data_arrays = Array{MX_float}[asarr(d.second) for d in data]
  else
    error("Invalid data argument type")
  end

  if isa(label, Union{NDArray, Array}) && eltype(label) <: Real
    label_names  = [:softmax_label]
    label_arrays = Array{MX_float}[asarr(label)]
  elseif isa(label, Pair)
    @assert isa(label.first, Base.Symbol) && isa(label.second, Union{NDArray, Array})
    label_names  = [label.first]
    label_arrays = Array{MX_float}[asarr(label.second)]
  elseif isa(label, Vector) || isa(label, Tuple)
    map(label) do d
      @assert isa(d, Pair) && isa(d.first, Base.Symbol) && isa(d.second, Union{NDArray, Array})
    end
    label_names  = Base.Symbol[d.first for d in label]
    label_arrays = Array{MX_float}[asarr(d.second) for d in label]
  else
    error("Invalid label argument type")
  end

  @assert length(data_arrays) > 0
  sample_count = size(data_arrays[1])[end]
  for i = 1:length(data_names)
    @assert(size(data_arrays[i])[end] == sample_count,
            "Number of samples in  $(data_names[i]) is mismatch with $(data_names[1])")
  end
  for i = 1:length(label_names)
    @assert(size(label_arrays[i])[end] == sample_count,
            "Number of samples in  $(label_names[i]) is mismatch with $(data_names[1])")
  end

  if batch_size == 0
    batch_size = sample_count
  end
  @assert 0 < batch_size <= sample_count

  function gen_batch_nds(arrs :: Vector{Array{MX_float}}, bsize :: Int)
    map(arrs) do arr
      shape = size(arr)
      empty(shape[1:end-1]..., bsize)
    end
  end

  data_batch  = gen_batch_nds(data_arrays, batch_size)
  label_batch = gen_batch_nds(label_arrays, batch_size)

  # reshape data and labels into 2D tensors, so that it is easier to work with them
  data_arrays = map(data_arrays) do arr
    reshape(arr, prod(size(arr)[1:end-1]), size(arr)[end])
  end
  label_arrays = map(label_arrays) do arr
    reshape(arr, prod(size(arr)[1:end-1]), size(arr)[end])
  end

  ArrayDataProvider(data_arrays, data_names, label_arrays, label_names, batch_size,
                    sample_count, shuffle, MX_float(data_padding), MX_float(label_padding),
                    data_batch, label_batch)
end

provide_data(provider::ArrayDataProvider) =
  collect(zip(provider.data_names, map(size, provider.data_batch)))

provide_label(provider::ArrayDataProvider) =
  collect(zip(provider.label_names, map(size, provider.label_batch)))

get_batch_size(provider::ArrayDataProvider) = provider.batch_size

struct ArrayDataProviderState <: AbstractDataProviderState
  curr_idx :: Int
end

Base.eltype(provider :: ArrayDataProvider) = ArrayDataProviderState

struct ArrayDataBatch <: AbstractDataBatch
  idx :: UnitRange{Int}
end

function _start(provider::ArrayDataProvider)
  if provider.shuffle
    # re-shuffle all data
    idx_perm = randperm(provider.sample_count)
    provider.data_arrays = map(x->x[:,idx_perm], provider.data_arrays)
    provider.label_arrays = map(x->x[:,idx_perm], provider.label_arrays)
  end

  return ArrayDataProviderState(1)
end

function Base.iterate(provider::ArrayDataProvider,
                      state::ArrayDataProviderState = _start(provider))
  (state.curr_idx > provider.sample_count) && return nothing
  idx = state.curr_idx:Base.min(state.curr_idx+provider.batch_size-1, provider.sample_count)
  return (ArrayDataBatch(idx), ArrayDataProviderState(idx.stop+1))
end

function count_samples(provider :: ArrayDataProvider, batch :: ArrayDataBatch)
  return length(batch.idx)
end

function get_data(provider :: ArrayDataProvider, batch :: ArrayDataBatch)
  for (src, dst) in zip(provider.data_arrays, provider.data_batch)
    copy_ignore_shape!(dst[1:length(batch.idx)], src[:, batch.idx])
    if length(batch.idx) < provider.batch_size
      dst[length(batch.idx)+1:provider.batch_size] = provider.data_padding
    end
  end
  return provider.data_batch
end
function get_label(provider :: ArrayDataProvider, batch :: ArrayDataBatch)
  for (src, dst) in zip(provider.label_arrays, provider.label_batch)
    copy_ignore_shape!(dst[1:length(batch.idx)], src[:, batch.idx])
    if length(batch.idx) < provider.batch_size
      dst[length(batch.idx)+1:provider.batch_size] = provider.label_padding
    end
  end
  return provider.label_batch
end

"""
    MXDataProvider

A data provider that wrap built-in data iterators from libmxnet. See below for
a list of built-in data iterators.
"""
mutable struct MXDataProvider <: AbstractDataProvider
  handle     :: MX_DataIterHandle
  data_shape :: Vector{Tuple{Symbol,Tuple}}
  label_shape:: Vector{Tuple{Symbol,Tuple}}
  batch_size :: Int

  # those two a auxiliary variables to help avoid calling reset
  # but still pre-fetch first batch to get shape information
  first_epoch:: Bool
  first_batch:: Bool
end

function _reset_data_iter(handle :: MX_DataIterHandle)
  @mxcall(:MXDataIterBeforeFirst, (MX_handle,), handle)
end
function _iter_next(handle :: MX_DataIterHandle)
  ref_ret = Ref{Cint}(0)
  @mxcall(:MXDataIterNext, (MX_handle, Ref{Cint}), handle, ref_ret)
  return Bool(ref_ret[])
end
function _get_data(handle :: MX_DataIterHandle)
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXDataIterGetData, (MX_handle, Ref{MX_handle}), handle, ref_hdr)
  return NDArray(MX_NDArrayHandle(ref_hdr[]), false)
end
function _get_label(handle :: MX_DataIterHandle)
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXDataIterGetLabel, (MX_handle, Ref{MX_handle}), handle, ref_hdr)
  return NDArray(MX_NDArrayHandle(ref_hdr[]), false)
end

function MXDataProvider(handle     :: MX_DataIterHandle;
                        data_name  :: Symbol = :data,
                        label_name :: Union{Symbol,Nothing} = :softmax_label,
                        kwargs...) # for convenience, we ignore the rest keyword arguments
  # init iterator, load the first batch and get shapes
  @assert(_iter_next(handle), "Failed to load the first batch in MXDataProvider")
  data_shape = Tuple{Base.Symbol, Tuple}[(data_name, size(_get_data(handle)))]
  if !isa(label_name, Nothing)
    label_shape = Tuple{Base.Symbol, Tuple}[(label_name::Base.Symbol, size(_get_label(handle)))]
  else
    label_shape = Tuple{Base.Symbol, Tuple}[]
  end

  MXDataProvider(handle, data_shape, label_shape, data_shape[1][2][end], true, true)
end

provide_data(provider::MXDataProvider) = provider.data_shape
provide_label(provider::MXDataProvider) = provider.label_shape
get_batch_size(provider::MXDataProvider) = provider.batch_size

mutable struct MXDataProviderState <: AbstractDataProviderState
  has_next :: Bool
end
struct MXDataBatch <: AbstractDataBatch
end

Base.eltype(::MXDataProvider) = MXDataBatch

function _start(provider::MXDataProvider)
  if !provider.first_epoch
    _reset_data_iter(provider.handle)
  else
    provider.first_epoch = false
  end

  return MXDataProviderState(true)
end

function _done(provider::MXDataProvider, state::MXDataProviderState)
  if provider.first_batch
    state.has_next = true
    provider.first_batch = false
  else
    state.has_next = _iter_next(provider.handle)
  end
  return !state.has_next
end

function Base.iterate(provider::MXDataProvider, state::MXDataProviderState = _start(provider))
  _done(provider, state) && return nothing
  MXDataBatch(), state
end

function get_data(provider :: MXDataProvider, batch :: MXDataBatch)
  return NDArray[_get_data(provider.handle)]
end
function get_label(provider :: MXDataProvider, batch :: MXDataBatch)
  return NDArray[_get_label(provider.handle)]
end
function count_samples(provider :: MXDataProvider, batch :: MXDataBatch)
  ref_pad = Ref{Cint}(0)
  @mxcall(:MXDataIterGetPadNum, (MX_handle, Ref{Cint}), provider.handle, ref_pad)
  return provider.batch_size - Int(ref_pad[])
end

function _get_iter_creators()
  n_ref = Ref{MX_uint}(0)
  h_ref = Ref{Ptr{MX_handle}}(0)
  @mxcall(:MXListDataIters, (Ref{MX_uint}, Ref{Ptr{MX_handle}}), n_ref, h_ref)

  return unsafe_wrap(Array, h_ref[], n_ref[])
end

function _get_iter_name(hdr :: MX_handle)
  ref_name      = Ref{char_p}(0)
  ref_desc      = Ref{char_p}(0)
  ref_narg      = Ref{MX_uint}(0)
  ref_arg_names = Ref{char_pp}(0)
  ref_arg_types = Ref{char_pp}(0)
  ref_arg_descs = Ref{char_pp}(0)

  @mxcall(:MXDataIterGetIterInfo,
          (MX_handle, Ref{char_p}, Ref{char_p}, Ref{MX_uint}, Ref{char_pp}, Ref{char_pp}, Ref{char_pp}),
          hdr, ref_name, ref_desc, ref_narg, ref_arg_names, ref_arg_types, ref_arg_descs)

  return Symbol(unsafe_string(ref_name[]))
end

const _iter_creator_cache = Dict{Symbol,MX_handle}()
function _populate_iter_creator_cache!()
  empty!(_iter_creator_cache)
  h_creators = _get_iter_creators()
  for handle in h_creators
    name = _get_iter_name(handle)
    _iter_creator_cache[name] = handle
  end
end

_get_iter_creator(name::Symbol) = _iter_creator_cache[name]

function _define_data_iter_creator(hdr :: MX_handle)
  ref_name      = Ref{char_p}(0)
  ref_desc      = Ref{char_p}(0)
  ref_narg      = Ref{MX_uint}(0)
  ref_arg_names = Ref{char_pp}(0)
  ref_arg_types = Ref{char_pp}(0)
  ref_arg_descs = Ref{char_pp}(0)

  @mxcall(:MXDataIterGetIterInfo,
          (MX_handle, Ref{char_p}, Ref{char_p}, Ref{MX_uint}, Ref{char_pp}, Ref{char_pp}, Ref{char_pp}),
          hdr, ref_name, ref_desc, ref_narg, ref_arg_names, ref_arg_types, ref_arg_descs)

  iter_name = Symbol(unsafe_string(ref_name[]))

  isprovider =  endswith(string(iter_name), "Iter")
  signature = _format_signature(Int(ref_narg[]), ref_arg_names)
  f_desc = "    " * string(iter_name) * "(" *signature * ")\n\n"
  if isprovider
    f_desc *= "Can also be called with the alias `$(string(iter_name)[1:end-4] * "Provider")`.\n"
  end
  f_desc *= unsafe_string(ref_desc[]) * "\n\n"
  f_desc *= "# Arguments:\n"
  f_desc *= "* `data_name::Symbol`: keyword argument, default `:data`. The name of the data.\n"
  f_desc *= "* `label_name::Symbol`: keyword argument, default `:softmax_label`. " *
            "The name of the label. Could be `nothing` if no label is presented in this dataset.\n\n"
  f_desc *= _format_docstring(Int(ref_narg[]), ref_arg_names, ref_arg_types, ref_arg_descs) * "\n"
  f_desc *= "Returns the constructed `MXDataProvider`."

  if isprovider
    alias_name = Symbol(string(iter_name)[1:end-4] * "Provider")
  else
    alias_name = nothing
  end

  defun = quote
    @doc $f_desc
    function $iter_name(; kwargs...)
      arg_keys = String[string(k) for (k,v) in kwargs]
      arg_vals = String[dump_mx_param(v) for (k,v) in kwargs]
      ref_hdr  = Ref{MX_handle}(0)

      local hdr = _get_iter_creator($(QuoteNode(iter_name)))
      @mxcall(:MXDataIterCreateIter, (MX_handle, MX_uint, char_pp, char_pp, Ref{MX_handle}),
              hdr, length(arg_keys), arg_keys, arg_vals, ref_hdr)

      return MXDataProvider(MX_DataIterHandle(ref_hdr[]); kwargs...)
    end
    $(isprovider ? :(const $alias_name = $iter_name) : :())

  end
  defun
end

macro _import_io_iterators()
  creators = _get_iter_creators()
  defs = Expr[]
  for handle in creators
    push!(defs, _define_data_iter_creator(handle))
  end
  esc(quote
    $(defs...)
  end)
end

@_import_io_iterators()

max_stars_repo_name: hsali/incubator-mxnet

max_issues_repo_name: hsali/incubator-mxnet

max_forks_repo_name: hsali/incubator-mxnet


text:
export MRIGARKExplicit
export MRIGARKERK33aSandu, MRIGARKERK45aSandu

"""
    MRIParam(p, γs, Rs, ts, Δts)

Construct a type for passing the data around for the `MRIGARKExplicit` explicit
time stepper to follow on methods. `p` is the original user defined ODE
parameters, `γs` and `Rs` are the MRI parameters and stage values, respectively.
`ts` and `Δts` are the stage time and stage time step.
"""
struct MRIParam{P, T, AT, N, M}
    p::P
    γs::NTuple{M, SArray{NTuple{1, N}, T, 1, N}}
    Rs::NTuple{N, AT}
    ts::T
    Δts::T
    function MRIParam(
        p::P,
        γs::NTuple{M},
        Rs::NTuple{N, AT},
        ts,
        Δts,
    ) where {P, M, N, AT}
        T = eltype(γs[1])
        new{P, T, AT, N, M}(p, γs, Rs, ts, Δts)
    end
end

# We overload get property to access the original param
function Base.getproperty(mriparam::MRIParam, s::Symbol)
    if s === :p
        p = getfield(mriparam, :p)
        return p isa MRIParam ? p.p : p
    else
        getfield(mriparam, s)
    end
end

"""
    MRIGARKExplicit(f!, fastsolver, Γs, γ̂s, Q, Δt, t0)

Construct an explicit MultiRate Infinitesimal General-structure Additive
Runge--Kutta (MRI-GARK) scheme to solve

```math
    \\dot{y} = f(y, t) + g(y, t)
```

where `f` is the slow tendency function and `g` is the fast tendency function;
see Sandu (2019).

The fast tendency is integrated using the `fastsolver` and the slow tendency
using the MRI-GARK scheme. Namely, at each stage the scheme solves

```math
               v(T_i) &= Y_i \\\\
             \\dot{v} &= f(v, t) + \\sum_{j=1}^{i} \\bar{γ}_{ij}(t) R_j \\\\
    \\bar{γ}_{ijk}(t) &= \\sum_{k=0}^{NΓ-1} γ_{ijk} τ(t)^k / Δc_s \\\\
                 τ(t) &= (t - t_s) / Δt \\\\
              Y_{i+1} &= v(T_i + c_s * Δt)
```

where ``Y_1 = y_n`` and ``y_{n+1} = Y_{Nstages+1}``.

Here ``R_j = g(Y_j, t_0 + c_j * Δt)`` is the tendency for stage ``j``,
``γ_{ijk}`` are the GARK coupling coefficients,
``NΓ`` is the number of sets of GARK coupling coefficients there are
``Δc_s = \\sum_{j=1}^{Nstages} γ_{sj1} = c_{s+1} - c_s`` is the scaling
increment between stage times. The ODE for ``v(t)`` is solved using the
`fastsolver`.  Note that this form of the scheme is based on Definition 2.2 of
Sandu (2019), but ODE for ``v(t)`` is written to go from ``t_s`` to
``T_i + c_s * Δt`` as opposed to ``0`` to ``1``.

Currently only [`LowStorageRungeKutta2N`](@ref) schemes are supported for
`fastsolver`

The coefficients defined by `γ̂s` can be used for an embedded scheme (only the
last stage is different).

The available concrete implementations are:

  - [`MRIGARKERK33aSandu`](@ref)
  - [`MRIGARKERK45aSandu`](@ref)

### References
 - [Sandu2019](@cite)
"""
mutable struct MRIGARKExplicit{T, RT, AT, Nstages, NΓ, FS, Nstages_sq} <:
               AbstractODESolver
    "time step"
    dt::RT
    "time"
    t::RT
    "elapsed time steps"
    steps::Int
    "rhs function"
    slowrhs!::Any
    "Storage for RHS during the `MRIGARKExplicit` update"
    Rstages::NTuple{Nstages, AT}
    "RK coefficient matrices for coupling coefficients"
    Γs::NTuple{NΓ, SArray{NTuple{2, Nstages}, RT, 2, Nstages_sq}}
    "RK coefficient matrices for embedded scheme"
    γ̂s::NTuple{NΓ, SArray{NTuple{1, Nstages}, RT, 1, Nstages}}
    "RK coefficient vector C (time scaling)"
    Δc::SArray{NTuple{1, Nstages}, RT, 1, Nstages}
    "fast solver"
    fastsolver::FS

    function MRIGARKExplicit(
        slowrhs!,
        fastsolver,
        Γs,
        γ̂s,
        Q::AT,
        dt,
        t0,
    ) where {AT <: AbstractArray}
        NΓ = length(Γs)
        Nstages = size(Γs[1], 1)
        T = eltype(Q)
        RT = real(T)

        # Compute the Δc coefficients
        Δc = sum(Γs[1], dims = 2)[:]

        # Scale in the Δc to the Γ and γ̂, and convert to real type
        Γs = ntuple(k -> RT.(Γs[k] ./ Δc), NΓ)
        γ̂s = ntuple(k -> RT.(γ̂s[k] / Δc[Nstages]), NΓ)

        # Convert to real type
        Δc = RT.(Δc)

        # create storage for the stage values
        Rstages = ntuple(i -> similar(Q), Nstages)

        FS = typeof(fastsolver)
        new{T, RT, AT, Nstages, NΓ, FS, Nstages^2}(
            RT(dt),
            RT(t0),
            0,
            slowrhs!,
            Rstages,
            Γs,
            γ̂s,
            Δc,
            fastsolver,
        )
    end
end

function dostep!(Q, mrigark::MRIGARKExplicit, param, time::Real)
    dt = mrigark.dt
    fast = mrigark.fastsolver

    Rs = mrigark.Rstages
    Δc = mrigark.Δc
    Nstages = length(Δc)
    slowrhs! = mrigark.slowrhs!
    Γs = mrigark.Γs
    NΓ = length(Γs)

    ts = time
    groupsize = 256
    for s in 1:Nstages
        # Stage dt
        dts = Δc[s] * dt

        p = param isa MRIParam ? param.p : param
        slowrhs!(Rs[s], Q, p, ts, increment = false)
        if param isa MRIParam
            # fraction of the step slower stage increment we are on
            τ = (ts - param.ts) / param.Δts
            event = Event(array_device(Q))
            event = mri_update_rate!(array_device(Q), groupsize)(
                realview(Rs[s]),
                τ,
                param.γs,
                param.Rs;
                ndrange = length(realview(Rs[s])),
                dependencies = (event,),
            )
            wait(array_device(Q), event)
        end

        γs = ntuple(k -> ntuple(j -> Γs[k][s, j], s), NΓ)
        mriparam = MRIParam(param, γs, realview.(Rs[1:s]), ts, dts)
        updatetime!(mrigark.fastsolver, ts)
        solve!(Q, mrigark.fastsolver, mriparam; timeend = ts + dts)

        # update time
        ts += dts
    end
end

@kernel function mri_update_rate!(dQ, τ, γs, Rs)
    i = @index(Global, Linear)
    @inbounds begin
        NΓ = length(γs)
        Ns = length(γs[1])
        dqi = dQ[i]

        for s in 1:Ns
            ri = Rs[s][i]
            sc = γs[NΓ][s]
            for k in (NΓ - 1):-1:1
                sc = sc * τ + γs[k][s]
            end
            dqi += sc * ri
        end

        dQ[i] = dqi
    end
end

"""
    MRIGARKERK33aSandu(f!, fastsolver, Q; dt, t0 = 0, δ = -1 // 2)

The 3rd order, 3 stage scheme from Sandu (2019). The parameter `δ` defaults to
the value suggested by Sandu, but can be varied.
"""
function MRIGARKERK33aSandu(slowrhs!, fastsolver, Q; dt, t0 = 0, δ = -1 // 2)
    T = eltype(Q)
    RT = real(T)
    #! format: off
    Γ0 = [
                1 // 3           0 // 1          0 // 1
        (-6δ - 7) // 12  (6δ + 11) // 12         0 // 1
                0 // 1    (6δ - 5) // 12  (3 - 2δ) // 4
    ]
    γ̂0 = [     1 // 12          -1 // 3          7 // 12]

    Γ1 = [
               0 // 1          0 // 1   0 // 1
        (2δ + 1) // 2  -(2δ + 1) // 2   0 // 1
               1 // 2  -(2δ + 1) // 2   δ // 1
    ]
    γ̂1 = [     0 // 1          0 // 1   0 // 1]
    #! format: on
    MRIGARKExplicit(slowrhs!, fastsolver, (Γ0, Γ1), (γ̂0, γ̂1), Q, dt, t0)
end

"""
    MRIGARKERK45aSandu(f!, fastsolver, Q; dt, t0 = 0)

The 4th order, 5 stage scheme from Sandu (2019).
"""
function MRIGARKERK45aSandu(slowrhs!, fastsolver, Q; dt, t0 = 0)
    T = eltype(Q)
    RT = real(T)
    #! format: off
    Γ0 = [
                  1 // 5                 0 // 1                 0 // 1               0 // 1             0 // 1
                -53 // 16              281 // 80                0 // 1               0 // 1             0 // 1
         -36562993 // 71394880    34903117 // 17848720  -88770499 // 71394880        0 // 1             0 // 1
          -7631593 // 71394880  -166232021 // 35697440    6068517 // 1519040   8644289 // 8924360       0 // 1
            277061 // 303808       -209323 // 1139280    -1360217 // 1139280   -148789 // 56964    147889 // 45120
    ]
    γ̂0 = [-1482837 // 759520        175781 // 71205       -790577 // 1139280     -6379 // 56964        47 // 96]
    Γ1 = [
               0 // 1                0 // 1               0 // 1               0 // 1             0 // 1
             503 // 80            -503 // 80              0 // 1               0 // 1             0 // 1
        -1365537 // 35697440   4963773 // 7139488  -1465833 // 2231090         0 // 1             0 // 1
        66974357 // 35697440  21445367 // 7139488        -3 // 1        -8388609 // 4462180       0 // 1
          -18227 // 7520             2 // 1               1 // 1               5 // 1        -41933 // 7520
    ]
    γ̂1 = [  6213 // 1880         -6213 // 1880            0 // 1               0 // 1             0 // 1]
    #! format: on
    MRIGARKExplicit(slowrhs!, fastsolver, (Γ0, Γ1), (γ̂0, γ̂1), Q, dt, t0)
end

max_forks_repo_name: jm-c/CLIMA


text:
# will be removed once merged to Tokenize
# https://github.com/JuliaLang/Tokenize.jl/pull/153
import Tokenize.Tokens.Token

Meta.parse(t::T) where {T <: Union{Token, Array{Token}}} = Meta.parse(untokenize(t))

"""
    UndefToken

An abstract type that shows if the an evaluated `Token` is defined.
"""
abstract type UndefToken end

"""
    tisdefined(t)
    tisdefined(modul, t)

Check if a `Token` that represent a `Symbol` is defined. Throws error if the `Token` is representing an `Expr`.

# Examples
```julia
julia> t = collect(tokenize("Int64"))

julia> tisdefined(t)
true
```
"""
function tisdefined(modul, t::T) where {T <: Union{Token, Array{Token}}}
    pt = Meta.parse(t)
    return isdefined(modul,pt)
end

tisdefined(t::T) where {T <: Union{Token, Array{Token}}} = tisdefined(@__MODULE__, t)

tisdefined(modul, t::T) where {T <: Union{Symbol, Expr}} = isdefined(modul, t)
tisdefined(t::T) where {T <: Union{Symbol, Expr}} = isdefined(@__MODULE__, t)

"""
    tevalfast(t::T)
    tevalfast(t::T, check_isdefined::Bool=false)
    tevalfast(modul = @__MODULE__, t::T, check_isdefined::Bool=false)

Parses and evaluates a `Token` that represents a `Symbol` or `Expr` in Julia. For `Symbol` It dispatches on the fast method based on the parsed Token.

If you set `check_isdefined` to `true`, and `t` is not defined in the scope it returns `UndefToken` instead of throwing an error.

# Examples
```julia
julia> t = collect(tokenize("Int64"))

julia> tevalfast(t)
Int64
```
"""
function tevalfast(modul::Module, t::T, check_isdefined::Bool = false) where {T <: Union{Token, Array{Token}}}
    pt = Meta.parse(t)
    if check_isdefined && !(tisdefined(modul,pt))
        return UndefToken
    end
    return evalfast(modul,pt)
end

tevalfast(t::T, check_isdefined::Bool = false) where {T <: Union{Token, Array{Token}}} = tevalfast(@__MODULE__, t, check_isdefined)

"""
    evalfast(x)
    evalfast(modul, x)

Evaluates `x` fast. Uses `getfield` if `x` is a `Symbol`, and uses `eval` if it is an `Expr`
# Examples
```julia
julia> evalfast(:Int64)

julia> evalfast(:([5]))

```
"""
evalfast(modul::Module, x::Expr)= modul.eval(x)
evalfast(modul::Module, x::Symbol)= getfield(modul,x)
evalfast(x::Expr)= Core.eval(@__MODULE__, x)
evalfast(x::Symbol)= getfield(@__MODULE__,x)

"""
    tgetfield(t::T)
    tgetfield(t::T, check_isdefined::Bool=false)
    tgetfield(modul = @__MODULE__, t::T, check_isdefined::Bool=false)

See [`tevalfast`](@ref) for fast Token evaluation.

Parses and evaluates a Token that represents a `Symbol` in Julia. For `Symbol` it is similar to eval, but much faster.

If you set `check_isdefined` to `true`, and `t` is not defined in the scope it returns `UndefToken` instead of throwing an error.

# Examples
```julia
julia> t = collect(tokenize("Int64"))

julia> tgetfield(t)
Int64
```
"""
function tgetfield(modul::Module, t::T, check_isdefined::Bool = false) where {T <: Union{Token, Array{Token}}}
    pt = Meta.parse(t)
    if check_isdefined && !(tisdefined(modul,pt))
        return UndefToken
    end
    return getfield(modul,pt)
end

tgetfield(t::T, check_isdefined::Bool = false) where {T <: Union{Token, Array{Token}}} = tgetfield(@__MODULE__, t, check_isdefined)

"""
    teval(t::T, check_isdefined::Bool=false)
    teval(t::T, check_isdefined::Bool=false)
    teval(modul = @__MODULE__, t::T, check_isdefined::Bool=false)

See [`tevalfast`](@ref) for fast Token evaluation.

Parses and evaluates a `Token`.

If you set `check_isdefined` to `true`, and `t` is not defined in the scope it returns `UndefToken` instead of throwing an error.

# Examples
```julia
julia> t = collect(tokenize("Int64"))

julia> teval(t)
Int64
```
"""
function teval(modul::Module, t::T, check_isdefined::Bool = false) where {T <: Union{Token, Array{Token}}}
    pt = Meta.parse(t)
    if check_isdefined && !(tisdefined(modul,pt))
        return UndefToken
    end
    return Core.eval(modul, pt)
end
teval(t::T, check_isdefined::Bool = false) where {T <: Union{Token, Array{Token}}} = teval(@__MODULE__, t, check_isdefined)

"""
    ttypeof(t)
    ttypeof(t::T, check_isdefined::Bool = false)

Returns the type of an evaluated Token

If you set `check_isdefined` to `true`, and `t` is not defined in the scope it returns `UndefToken` instead of throwing an error.

# Examples
```julia
julia> t = collect(tokenize("Int64"))

julia> ttypeof(t)
DataType
```
"""
ttypeof(t::T, check_isdefined::Bool = false) where {T <: Union{Token, Array{Token}}} = typeof(tevalfast(t, check_isdefined))

"""
    tisa(t, T::Type)
    tisa(t::T, Tspecified::Type, check_isdefined::Bool = false)

Compares the specified type with the type of an evaluated Token

If you set `check_isdefined` to `true`, and `t` is not defined in the scope it returns `UndefToken` instead of throwing an error.

# Examples
```julia
julia> t = collect(tokenize("Int64"))

julia> tisa(t, DataType)
true
```
"""
tisa(t::T, Tspecified::Type, check_isdefined::Bool = false) where {T <: Union{Token, Array{Token}}} = isa(tevalfast(t, check_isdefined), Tspecified)


text:
function TableauHSPARKSymmetricProjection(name, q::Tableau{T}, p::Tableau{T}, d=[]; R∞=1) where {T}

    @assert q.s == p.s

    o = min(q.o, p.o)

    a_q = q.a
    a_p = p.a

    α_q = zeros(T, q.s, 2)
    α_q[:,1] .= 0.5

    α_p = zeros(T, p.s, 2)
    α_p[:,1] .= 0.5

    a_q̃ = vcat(zero(q.b)', q.b')
    a_p̃ = vcat(zero(p.b)', p.b')

    α_q̃ = [[0.0  0.0]
           [0.5  0.5]]
    α_p̃ = [[0.0  0.0]
           [0.5  0.5]]

    b_q = q.b
    b_p = p.b
    β_q = [0.5, 0.5]
    β_p = [0.5, 0.5]

    c_q = q.c
    c_p = p.c
    c_λ = [ 0.0, 1.0]
    d_λ = [ 0.5, 0.5*R∞]

    ω_λ = reshape(T[1  1  0], (1,3))
    δ_λ = reshape(T[-1   R∞], (1,2))

    if length(d) == 0
        return TableauHSPARKprimary(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ)
    else
        @assert length(d) == q.s == p.s

        return TableauHSPARKprimary(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ, d)
    end

end

"Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and symmetric projection."
function TableauHSPARKLobattoIIIAIIIBpSymmetric(s)
    TableauHSPARKSymmetricProjection(Symbol("HSPARKLobattoIIIAIIIB($s)pSymmetric"), TableauLobattoIIIA(s), TableauLobattoIIIB(s), get_lobatto_nullvector(s))
end

"Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and symmetric projection."
function TableauHSPARKLobattoIIIBIIIApSymmetric(s)
    TableauHSPARKSymmetricProjection(Symbol("HSPARKLobattoIIIBIIIA($s)pSymmetric"), TableauLobattoIIIB(s), TableauLobattoIIIA(s), get_lobatto_nullvector(s))
end

"Tableau for Gauss-Legendre method with s stages and symplectic projection."
function TableauHSPARKGLRKpSymmetric(s)
    glrk = TableauGauss(s)
    TableauHSPARKSymmetricProjection(Symbol("HSPARKGLRK($s)"), glrk, glrk)
end


text:
"""
    MLFlowClient

[MLFlowClient](https://github.com/JuliaAI.jl) is a [Julia](https://julialang.org/) package for working with [MLFlow](https://mlflow.org/) using the REST [API v2.0](https://www.mlflow.org/docs/latest/rest-api.html).

`MLFlowClient` allows you to create and manage `MLFlow` experiments, runs, and log metrics and artifacts. If you are not familiar with `MLFlow` and its concepts, please refer to [MLFlow documentation](https://mlflow.org/docs/latest/index.html).

# Limitations

- no authentication support.
- when storing artifacts, the assumption is that MLFlow and this library run on the same server. Artifacts are stored using plain filesystem operations. Therefore, `/mlruns` or the specified `artifact_location` must be accessible to both the MLFlow server (read), and this library (write).
"""
module MLFlowClient

using Dates
using UUIDs
using HTTP
using URIs
using JSON
using ShowCases
using FilePathsBase: AbstractPath

include("types.jl")
export
    MLFlow,
    MLFlowExperiment,
    MLFlowRunStatus,
    MLFlowRunInfo,
    get_run_id,
    MLFlowRunData,
    get_params,
    MLFlowRunDataMetric,
    MLFlowRun,
    get_info,
    get_data,
    MLFlowArtifactFileInfo,
    MLFlowArtifactDirInfo,
    get_path,
    get_size

include("utils.jl")
export
    generatefilterfromparams

include("experiments.jl")
export
    createexperiment,
    getexperiment,
    getorcreateexperiment,
    deleteexperiment,
    listexperiments

include("runs.jl")
export
    createrun,
    getrun,
    updaterun,
    deleterun,
    searchruns,
    logparam,
    logmetric,
    logartifact,
    listartifacts

end

max_stars_repo_name: JuliaAI/MLFlowClient

max_issues_repo_name: JuliaAI/MLFlowClient

max_forks_repo_name: JuliaAI/MLFlowClient


text:
using Random
using StatsBase
using Distributions
using StatsPlots

include("helper_functions.jl")

begin
    λ1 = 1
    λ2 = 1.5
    α = 1
    β = 1
    μ(α) = 1 / (1.5 * (α))
    ν(α) = 1 / (0.7 * (α))
    parameters = (λ1, λ2, α, μ, ν)
    he_parameters = (λ1, λ2, α, μ, ν, β)
    # opponent_parameters = (λ2, λ1, α, μ, ν)
    # opponent_he_parameters = (λ2, λ1, α, μ, ν, β)
    # qs = 0:0.001:1
end

Hyperexponential(parameters, weights) = MixtureModel(Exponential, parameters, weights)
RTDist(p, q, λ1, λ2, α, μ, ν) = Hyperexponential([(1 / (ν(α) - λ1 * p - λ2 * q)), (1 / (μ(α) - λ1 * (1 - p)))], [p, 1 - p])
RT = RTDist(0.57, 0.57, 1, 1, 1, μ, ν)
RT2 = RTDist(0.57, 0.57, 0.5, 0.5, 2, μ, ν)
density(log.(rand(RT, 1000000)), lw = 4, label = "Theoretical Response Distribution")
density!(log.(rand(RT2, 1000000)), lw = 4, label = "Theoretical Response Distribution")

JSs = 0.1:0.01:1
AR_JS_ratio = 1
ARs = AR_JS_ratio ./ JSs

begin
RTs = RTDist.(0.57, 0.57, ARs, ARs, JSs, μ, ν)
plot_mean = plot(JSs,
    mean.(RTs),
    title = "Mean Response Time by Job Size",
    xlabel = "Job Size",
    ylabel = "Mean Response Time",
    # label = "Theoretical Mean",
    size = [400, 300],
    legend = :none)

plot_variance = plot(JSs,
    var.(RTs),
    title = "Variance by Job Size",
    xlabel = "Job Size",
    ylabel = "Variance of Response Time",
    # label = "Theoretical Variance",
    size = [400, 300],
    legend = :none)

mean_var_plot = plot(plot_mean, plot_variance, size = [800, 300], margin = 12pt)
end

savefig("outputs/theoretical_mean_var.pdf")


text:
module SpglibTests
using Spglib
using FactCheck

facts("Get Symmetry of simple lattices") do
  cell = [0 0.5 0.5; 0.5 0 0.5; 0.5 0.5 0]
  positions = transpose([0 0 0; 0.25 0.25 0.25])
  facts("diamond") do
    types = [:Si, :Si]
    symops = symmetry_operations(cell, positions, types)
    @fact length(symops) --> 48
    symbol = international_symbol(cell, positions, types)
    @fact symbol --> "Fd-3m"
    symbol = schoenflies_symbol(cell, positions, types)
    @fact symbol --> "Oh^7"
  end

  facts("zincblende") do
    types = [:Zn, :S]
    symops = symmetry_operations(cell, positions, types)
    @fact length(symops) --> 24
    symbol = international_symbol(cell, positions, types)
    @fact symbol --> "F-43m"
    symbol = schoenflies_symbol(cell, positions, types)
    @fact symbol --> "Td^2"
  end
end
exitstatus()
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Testing

@testset "2d poisson problem with known analytical solution" begin
    # from FENiCS tutorial, u(x,y) = 1 + x² + 2y² on [0x1]×[0,1]
    # and u₀(x,y) = 1 + x² + 2y², f(x,y) = -6

    mesh_file = Pkg.dir("JuliaFEM")*"/test/testdata/primitives.med"
    mesh = aster_read_mesh(mesh_file, "UNITSQUARE_6X4")

    field = Problem(Heat, "unit square, 6x4 triangular mesh", 1)
    field.elements = create_elements(mesh, "UNITSQUARE")
    field.properties.formulation = "2D"
    update!(field, "thermal conductivity", 1.0)
    update!(field, "heat source", -6.0)

    bc = Problem(Dirichlet, "u₀(x,y) = 1 + x² + 2y²", 1, "temperature")
    #bc.properties.order = 2
    #bc.properties.dual_basis = true
    bc.properties.variational = false
    bc.elements = create_elements(mesh, "FACE1", "FACE2", "FACE3", "FACE4")
    function u0(element, ip, time)
        x, y = element("geometry", ip, time)
        return 1 + x^2 + 2*y^2
    end
    update!(bc, "temperature 1", u0)

    solver = LinearSolver(field, bc)
    solver()

    T_fem = Float64[]
    T_acc = Float64[]
    for (nid, X) in field("geometry")
        push!(T_fem, field("temperature", X)[1])
        push!(T_acc, 1.0 + X[1]^2 + 2*X[2]^2)
    end

    @test maximum(abs(T_fem-T_acc)) < 1.0e-12

    # gradient of field is
    gradT(X) = [2*X[1] 4*X[2]]
    X = [0.5, 0.5]
    gradT1 = gradT(X)
    gradT2 = field("temperature", X, solver.time, Val{:Grad})
    info("gradT1 = $gradT1, gradT2 = $gradT2")
    # [1.1666666666666625 1.5000000000000018] quite big difference ..?
    @test isapprox(gradT1, gradT2; rtol=25.0e-2)
end


text:
# # Tilted bottom boundary layer example
#
# This example simulates a two-dimensional oceanic bottom boundary layer
# in a domain that's tilted with respect to gravity. We simulate the perturbation
# away from a constant along-slope (y-direction) velocity constant density stratification.
# This perturbation develops into a turbulent bottom boundary layer due to momentum
# loss at the bottom boundary modeled with a quadratic drag law.
#
# This example illustrates
#
#   * changing the direction of gravitational acceleration in the buoyancy model;
#   * changing the axis of rotation for Coriolis forces.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, Plots"
# ```
#
# ## The domain
#
# We create a ``400 × 100`` meter ``x, z`` grid with ``128 × 32`` cells
# and finer resolution near the bottom,

using Oceananigans

Lx = 400 # m
Lz = 100 # m
Nx = 128
Nz = 64

## Creates a grid with near-constant spacing `refinement * Lz / Nz`
## near the bottom:
refinement = 1.8 # controls spacing near surface (higher means finer spaced)
stretching = 10  # controls rate of stretching at bottom

## "Warped" height coordinate
h(k) = (Nz + 1 - k) / Nz

## Linear near-surface generator
ζ(k) = 1 + (h(k) - 1) / refinement

## Bottom-intensified stretching function
Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching))

## Generating function
z_faces(k) = - Lz * (ζ(k) * Σ(k) - 1)

grid = RectilinearGrid(topology = (Periodic, Flat, Bounded),
                       size = (Nx, Nz),
                       x = (0, Lx),
                       z = z_faces,
                       halo = (3, 3))

# Let's make sure the grid spacing is both finer and near-uniform at the bottom,

using Plots

plot(grid.Δzᵃᵃᶜ[1:Nz], grid.zᵃᵃᶜ[1:Nz],
     marker = :circle,
     ylabel = "Depth (m)",
     xlabel = "Vertical spacing (m)",
     legend = nothing)

# ## Tilting the domain
#
# We use a domain that's tilted with respect to gravity by

θ = 3 # degrees

# so that ``x`` is the along-slope direction, ``z`` is the across-sloce direction that
# is perpendicular to the bottom, and the unit vector anti-aligned with gravity is

ĝ = (sind(θ), 0, cosd(θ))

# Changing the vertical direction impacts both the `gravity_unit_vector`
# for `Buoyancy` as well as the `rotation_axis` for Coriolis forces,

buoyancy = Buoyancy(model=BuoyancyTracer(), gravity_unit_vector=ĝ)
coriolis = ConstantCartesianCoriolis(f=1e-4, rotation_axis=ĝ)

# where we have used a constant Coriolis parameter ``f = 10⁻⁴ \rm{s}⁻¹``.
# The tilting also affects the kind of density stratified flows we can model.
# In particular, a constant density stratification in the tilted
# coordinate system

@inline constant_stratification(x, y, z, t, p) = p.N² * (x * p.ĝ[1] + z * p.ĝ[3])

# is _not_ periodic in ``x``. Thus we cannot explicitly model a constant stratification
# on an ``x``-periodic grid such as the one used here. Instead, we simulate periodic
# _perturbations_ away from the constant density stratification by imposing
# a constant stratification as a `BackgroundField`,

B_field = BackgroundField(constant_stratification, parameters=(; ĝ, N² = 1e-5))

# where ``N² = 10⁻⁵ \rm{s}⁻¹`` is the background buoyancy gradient.

# ## Bottom drag
#
# We impose bottom drag that follows Monin-Obukhov theory.
# We include the background flow in the drag calculation,
# which is the only effect the background flow enters the problem,

V∞ = 0.1 # m s⁻¹
z₀ = 0.1 # m (roughness length)
κ = 0.4 # von Karman constant
z₁ = znodes(Center, grid)[1] # Closest grid center to the bottom
cᴰ = (κ / log(z₁ / z₀))^2 # Drag coefficient

@inline drag_u(x, y, t, u, v, p) = - p.cᴰ * √(u^2 + (v + p.V∞)^2) * u
@inline drag_v(x, y, t, u, v, p) = - p.cᴰ * √(u^2 + (v + p.V∞)^2) * (v + p.V∞)

drag_bc_u = FluxBoundaryCondition(drag_u, field_dependencies=(:u, :v), parameters=(; cᴰ, V∞))
drag_bc_v = FluxBoundaryCondition(drag_v, field_dependencies=(:u, :v), parameters=(; cᴰ, V∞))

u_bcs = FieldBoundaryConditions(bottom = drag_bc_u)
v_bcs = FieldBoundaryConditions(bottom = drag_bc_v)

# ## Create the `NonhydrostaticModel`
#
# We are now ready to create the model. We create a `NonhydrostaticModel` with an
# `UpwindBiasedFifthOrder` advection scheme, a `RungeKutta3` timestepper,
# and a constant viscosity and diffusivity.

ν = 1e-4 # m² s⁻¹, small-ish
closure = ScalarDiffusivity(ν=ν, κ=ν)

model = NonhydrostaticModel(; grid, buoyancy, coriolis, closure,
                            timestepper = :RungeKutta3,
                            advection = UpwindBiasedFifthOrder(),
                            tracers = :b,
                            boundary_conditions = (u=u_bcs, v=v_bcs),
                            background_fields = (; b=B_field))

# ## Create and run a simulation
#
# We are now ready to create the simulation. We begin by setting the initial time step
# conservatively, based on the smallest grid size of our domain and set-up a

using Oceananigans.Units
using Oceananigans.Grids: min_Δz

simulation = Simulation(model, Δt = 0.5 * min_Δz(grid) / V∞, stop_time = 2days)

# We use `TimeStepWizard` to adapt our time-step and print a progress message,

using Printf

wizard = TimeStepWizard(max_change=1.1, cfl=0.7)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(4))

start_time = time_ns() # so we can print the total elapsed wall time

progress_message(sim) =
    @printf("Iteration: %04d, time: %s, Δt: %s, max|w|: %.1e m s⁻¹, wall time: %s\n",
            iteration(sim), prettytime(time(sim)),
            prettytime(sim.Δt), maximum(abs, sim.model.velocities.w),
            prettytime((time_ns() - start_time) * 1e-9))

simulation.callbacks[:progress] = Callback(progress_message, IterationInterval(100))

# ## Add outputs to the simulation
#
# We add outputs to our model using the `NetCDFOutputWriter`,

u, v, w = model.velocities
b = model.tracers.b
B = model.background_fields.tracers.b

total_b = b + B
total_v = v + V∞
ω_y = ∂z(u) - ∂x(w)

outputs = (; u, total_v, w, total_b, ω_y)

simulation.output_writers[:fields] = NetCDFOutputWriter(model, outputs;
                                                        filename = joinpath(@__DIR__, "tilted_bottom_boundary_layer.nc"),
                                                        schedule = TimeInterval(20minutes),
                                                        overwrite_existing = true)

# Now we just run it!

run!(simulation)

# ## Visualize the results
#
# First we load the required package to load NetCDF output files and define the coordinates for
# plotting using existing objects:

using NCDatasets

xω, zω = xnodes(ω_y), znodes(ω_y)
xv, zv = xnodes(total_v), znodes(total_v)

kwargs = (xlabel = "along-slope distance",
          zlabel = "across-slope distance",
          fill = true,
          levels = 20,
          linewidth = 0,
          color = :balance,
          aspect = :equal,
          colorbar = true,
          xlim = (0, Lx),
          zlim = (0, Lz))

# Read in the simulation's `output_writer` for the two-dimensional fields and then create an
# animation showing the ``y``-component of vorticity.

ds = NCDataset(simulation.output_writers[:fields].filepath, "r")

anim = @animate for (frame, time) in enumerate(ds["time"])
    @info "Plotting frame $frame of $(length(ds["time"]))..."

    ω_y = ds["ω_y"][:, 1, :, frame]
    total_v = ds["total_v"][:, 1, :, frame]

    ω_max = maximum(abs, ω_y)
    v_max = maximum(abs, total_v)

    plot_ω = contour(xω, zω, ω_y', clim = (-0.015, +0.015),
                     title = "y-vorticity, ω_y at t = " * prettytime(time);
                     kwargs...)

    plot_v = contour(xv, zv, total_v', clim = (-V∞, +V∞),
                     title = "Total along-slope velocity at t = " * prettytime(time);
                     kwargs...)

    plot(plot_ω, plot_v, layout = (2, 1), size = (800, 440))
end

mp4(anim, "tilted_bottom_boundary_layer.mp4", fps=12)

# It's always good practice to close the NetCDF files when we are done.

close(ds)


text:
using LaTeXStrings

import OceanTurbulenceParameterEstimation.TurbulenceClosureParameters: closure_with_parameters

parameter_guide = Dict(:Cᴰ => (name = "Dissipation parameter (TKE equation)", latex = L"C^D",
          default = 2.9079, bounds = (0.0, 10.0)), :Cᴸᵇ => (name = "Mixing length parameter", latex = L"C^{\ell}_b",
          default = 1.1612, bounds = (0.0, 10.0)), :Cᵟu => (name = "Ratio of mixing length to grid spacing", latex = L"C^{\delta}u",
          default = 0.5, bounds = (0.0, 3.0)), :Cᵟc => (name = "Ratio of mixing length to grid spacing", latex = L"C^{\delta}c",
          default = 0.5, bounds = (0.0, 3.0)), :Cᵟe => (name = "Ratio of mixing length to grid spacing", latex = L"C^{\delta}e",
          default = 0.5, bounds = (0.0, 3.0)), :Cᵂu★ => (name = "Mixing length parameter", latex = L"C^W_{u\star}",
          default = 3.6188, bounds = (0.0, 20.0)), :CᵂwΔ => (name = "Mixing length parameter", latex = L"C^Ww\Delta",
          default = 1.3052, bounds = (0.0, 5.0)), :CᴷRiʷ => (name = "Stability function parameter", latex = L"C^KRi^w",
          default = 0.7213, bounds = (0.0, 5.0)), :CᴷRiᶜ => (name = "Stability function parameter", latex = L"C^KRi^c",
          default = 0.7588, bounds = (-1.5, 4.0)), :Cᴷu⁻ => (name = "Velocity diffusivity LB", latex = L"C^Ku^-",
          default = 0.1513, bounds = (0.0, 2.0)), :Cᴷuʳ => (name = "Velocity diffusivity (UB-LB)/LB", latex = L"C^Ku^r",
          default = 3.8493, bounds = (0.0, 50.0)), :Cᴷc⁻ => (name = "Velocity diffusivity LB", latex = L"C^Kc^-",
          default = 0.3977, bounds = (0.0, 5.0)), :Cᴷcʳ => (name = "Velocity diffusivity (UB-LB)/LB", latex = L"C^Kc^r",
          default = 3.4601, bounds = (0.0, 50.0)), :Cᴷe⁻ => (name = "Velocity diffusivity LB", latex = L"C^Ke^-",
          default = 0.1334, bounds = (0.0, 3.0)), :Cᴷeʳ => (name = "Velocity diffusivity (UB-LB)/LB", latex = L"C^Ke^r",
          default = 8.1806, bounds = (0.0, 50.0)), :Cᴬu => (name = "Convective adjustment velocity parameter", latex = L"C^A_U",
          default = 0.0057, bounds = (0.0, 0.2)), :Cᴬc => (name = "Convective adjustment tracer parameter", latex = L"C^A_C",
          default = 0.6706, bounds = (0.0, 2.0)), :Cᴬe => (name = "Convective adjustment TKE parameter", latex = L"C^A_E",
          default = 0.2717, bounds = (0.0, 2.0)),
)

bounds(name) = parameter_guide[name].bounds
default(name) = parameter_guide[name].default

function named_tuple_map(names, f)
     names = Tuple(names)
     return NamedTuple{names}(f.(names))
end

"""
    ParameterSet()

Parameter set containing the names `names` of parameters, and a
NamedTuple `settings` mapping names of "background" parameters
to their fixed values to be maintained throughout the calibration.
"""
struct ParameterSet{N,S}
        names :: N
     settings :: S
end

"""
    ParameterSet(names::Set; nullify = Set())

Construct a `ParameterSet` containing all of the information necessary
to build a closure with the specified default parameters and settings,
given a Set `names` of the parameter names to be tuned, and a Set `nullify`
of parameters to be set to zero.
"""
function ParameterSet(names::Set; nullify = Set())
     zero_set = named_tuple_map(nullify, name -> 0.0)
     bkgd_set = named_tuple_map(keys(parameter_guide), name -> default(name))
     settings = merge(bkgd_set, zero_set) # order matters: `zero_set` overrides `bkgd_set`
     return ParameterSet(Tuple(names), settings)
end

names(ps::ParameterSet) = ps.names

###
### Define some convenient parameter sets based on the present CATKE formulation
###

required_params = Set([:Cᵟu, :Cᵟc, :Cᵟe, :Cᴰ, :Cᴸᵇ, :Cᵂu★, :CᵂwΔ, :Cᴷu⁻, :Cᴷc⁻, :Cᴷe⁻])
ri_depen_params = Set([:CᴷRiʷ, :CᴷRiᶜ, :Cᴷuʳ, :Cᴷcʳ, :Cᴷeʳ])
conv_adj_params = Set([:Cᴬu, :Cᴬc, :Cᴬe])

CATKEParametersRiDependent = ParameterSet(union(required_params, ri_depen_params); nullify = conv_adj_params)
CATKEParametersRiIndependent = ParameterSet(union(required_params); nullify = union(conv_adj_params, ri_depen_params))
CATKEParametersRiDependentConvectiveAdjustment = ParameterSet(union(required_params, conv_adj_params, ri_depen_params))
CATKEParametersRiIndependentConvectiveAdjustment = ParameterSet(union(required_params, conv_adj_params); nullify = ri_depen_params)


text:
#=

# 2D Lid-driven cavity (AD-Newton)
([source code](SOURCE_URL))

This example solves the lid-driven cavity problem where one seeks
a velocity ``\mathbf{u}`` and pressure ``\mathbf{p}`` of the incompressible Navier--Stokes problem
```math
\begin{aligned}
- \mu \Delta \mathbf{u} + (\mathbf{u} \cdot \nabla) \mathbf{u} + \nabla p & = 0\\
\mathrm{div}(u) & = 0
\end{aligned}
```
where ``\mathbf{u} = (1,0)`` along the top boundary of a square domain.

This examples highlights the use of automatic differentation to obtain Newton derivatives of nonlinear PDEOperators. The user can
switch between a Newton obtained by automatic differentation (ADnewton = true) and a 'manual' Newton scheme (ADnewton = false). The runtime
of the manual newton is slightly faster, but requires much more code input from the user that is more prone to errors.

=#

using PDELib
using Printf

function fe_liddrivencavity_autonewton(; verbosity = 0, Plotter = nothing, ADnewton = true, nref=2)

    function boundary_data_top!(result)
        result[1] = 1;
        result[2] = 0;
    end

    ## grid
    xgrid = uniform_refine(grid_unitsquare(Triangle2D), nref);

    ## problem parameters
    viscosity = 1e-2
    maxIterations = 50  # termination criterion 1 for nonlinear mode
    maxResidual = 1e-12 # termination criterion 2 for nonlinear mode
    broken_p = false

    ## choose one of these (inf-sup stable) finite element type pairs
    #FETypes = [H1P2{2,2}, H1P1{1}] # Taylor--Hood
    #FETypes = [H1P2B{2,2}, H1P1{1}]; broken_p = true # P2-bubble
    #FETypes = [H1CR{2}, H1P0{1}] # Crouzeix--Raviart
    #FETypes = [H1MINI{2,2}, H1P1{1}] # MINI element on triangles only
    FETypes = [H1BR{2}, H1P0{1}]; broken_p = true # Bernardi--Raugel

    #####################################################################################
    #####################################################################################

    ## negotiate data functions to the package
    user_function_bnd = DataFunction(boundary_data_top!, [2,2]; name = "u_bnd", dependencies = "", quadorder = 0)

    ## load linear Stokes problem prototype and assign data
    ## we are adding the nonlinar convection term ourself below
    ## to discuss the details
    StokesProblem = IncompressibleNavierStokesProblem(2; viscosity = viscosity, nonlinear = false)
    add_boundarydata!(StokesProblem, 1, [1,2,4], HomogeneousDirichletBoundary)
    add_boundarydata!(StokesProblem, 1, [3], BestapproxDirichletBoundary; data = user_function_bnd)

    ## store matrix of Laplace operator for nonlinear solver
    StokesProblem.LHSOperators[1,1][1].store_operator = true

    ## add Newton for convection term
    if ADnewton
        ## AUTOMATIC DIFFERENTATION
        ## requries kernel function for NonlinearForm action (u * grad) u
        function ugradu_kernel_AD(result, input)
            ## input = [u, grad(u)]
            ## compute (u * grad) u = grad(u)*u
            for j = 1 : 2
                result[j] = 0.0
                for k = 1 : 2
                    result[j] += input[k]*input[2 + (j-1)*2+k]
                end
            end
            return nothing
        end
        action_kernel = ActionKernel(ugradu_kernel_AD, [2,6]; dependencies = "", quadorder = 1)

        ## generate and add nonlinear PDEOperator (modifications to RHS are taken care of automatically)
        NLConvectionOperator = GenerateNonlinearForm("(u * grad) u  * v", [Identity, Gradient], [1,1], Identity, action_kernel; ADnewton = true)
        add_operator!(StokesProblem, [1,1], NLConvectionOperator)
    else
        ## MANUAL DIFFERENTATION
        ## uses the following kernel function for the linearised convection operator
        function ugradu_kernel_nonAD(result, input_current, input_ansatz)
            ## input_current = [current, grad(current)]
            ## input_ansatz = [ansatz, grad(ansatz)]
            ## compute (current * grad) ansatz + (current * grad) ansatz
            for j = 1 : 2
                result[j] = 0.0
                for k = 1 : 2
                    result[j] += input_current[k]*input_ansatz[2 + (j-1)*2+k]
                    result[j] += input_ansatz[k]*input_current[2 + (j-1)*2+k]
                end
            end
            return nothing
        end
        ## and a similar kernel function for the right-hand side
        function ugradu_kernel_rhs(result, input_current)
            ## input_current = [current, grad(current)]
            ## input_ansatz = [ansatz, grad(ansatz)]
            ## compute (current * grad) current
            for j = 1 : 2
                result[j] = 0.0
                for k = 1 : 2
                    result[j] += input_current[k]*input_current[2 + (j-1)*2+k]
                end
            end
            return nothing
        end
        newton_action_kernel = NLActionKernel(ugradu_kernel_nonAD, [2,6]; dependencies = "", quadorder = 1)
        action_kernel_rhs = ActionKernel(ugradu_kernel_rhs, [2,6]; dependencies = "", quadorder = 1)

        ## generate and add nonlinear PDEOperator (modifications to RHS are taken care of by optional arguments)
        NLConvectionOperator = GenerateNonlinearForm("(u * grad) u  * v", [Identity, Gradient], [1,1], Identity, newton_action_kernel; ADnewton = false, action_kernel_rhs = action_kernel_rhs)
        add_operator!(StokesProblem, [1,1], NLConvectionOperator)
    end

    ## generate FESpaces
    FESpaceVelocity = FESpace{FETypes[1]}(xgrid)
    FESpacePressure = FESpace{FETypes[2]}(xgrid; broken = broken_p)
    Solution = FEVector{Float64}("Stokes velocity",FESpaceVelocity)
    append!(Solution,"Stokes pressure",FESpacePressure)

    ## show configuration and solve Stokes problem
    Base.show(StokesProblem)
    GradientRobustMultiPhysics.solve!(Solution, StokesProblem; verbosity = verbosity, maxIterations = maxIterations, maxResidual = maxResidual)

    ## plot
    GradientRobustMultiPhysics.plot(Solution, [1,2], [Identity, Identity]; Plotter = Plotter, verbosity = verbosity, use_subplots = true)
    true
end


text:
module Mypackage

# Write your package code here.
greet()=print("Hello world")

include("myfunction.jl")
export(my_f)

end