text:
export map, mapvec, map!, map!r, map2!, mapmap, mapmapvec, mapi, mapveci, work, worki
export map2, map3, map4, map5
export mapvec2, mapvec3, mapvec4, mapvec5
export work2, work3, work4, work5
export mapprogress, mapkeys, mapvalues
export share, unshare
export shmap, shmap!, shmap!r, shmap2!, shwork
export pmap, pmap!, pmap!r, pmap2!, pwork
export pmapvec
export lmap, lmap!, lmap!r, lmap2!, lwork
export lmapvec
export hmap, hmap!, hmap!r, hmap2!, hwork
export amap, amap2, amapvec2
export table, ptable, ltable, htable, shtable, tableany, ptableany, ltableany, htableany, shtableany
export sort, sortrev, sortpermrev, uniq, filter, select, reject
export groupdict, groupby
export minelem, maxelem, extremaelem
export isany, areall
export tee
export *
export typed
export call, apply
export fold

typealias Callable Union{Function, Type}

import Base.sort
sort(a, f::Callable; kargs...) = part(a, sortperm(vec(map(a, f)); kargs...))
sort(a, key; kargs...) = part(a, sortperm(vec(extract(a, key)); kargs...))
sortrev(a) = sort(a; rev = true)
sortpermrev(a) = sortperm(a; rev = true)
sortrev(a, f) = sort(a,f; rev = true)

function minelem(a, f::Callable)
    b = map(a,f)
    at(a, findfirst(b .== minimum(b)))
end

function maxelem(a, f::Callable)
    b = map(a,f)
    at(a, findfirst(b .== maximum(b)))
end

function extremaelem(a, f::Callable)
    b = map(a,f)
    mi = minimum(b)
    ma = maximum(b)
    [at(a, findfirst(b .== mi)), at(a, findfirst(b .== ma))]
end

function uniq(a,f = id)
    d = Dict{Any,Int}()
    h = f == id ? a : @p mapvec a f
    for i = len(h):-1:1
        d[at(h,i)] = i
    end
    @p vec d | sort snd | map snd | part a _
end

#######################################
## map, pmap

mapvec(a, f::Callable) = [f(trytoview(a,i)) for i in 1:len(a)]
mapvec2(a, b, f::Callable) = [f(trytoview(a,i),trytoview(b,i)) for i in 1:len(a)]
mapvec3(a, b, c, f::Callable) = [f(trytoview(a,i),trytoview(b,i),trytoview(c,i)) for i in 1:len(a)]
mapvec4(a, b, c, d, f::Callable) = [f(trytoview(a,i),trytoview(b,i),trytoview(c,i),trytoview(d,i)) for i in 1:len(a)]
mapvec5(a, b, c, d, e, f::Callable) = [f(trytoview(a,i),trytoview(b,i),trytoview(c,i),trytoview(d,i),trytoview(e,i)) for i in 1:len(a)]

map2(a, b, f::Callable) = flatten(mapvec2(a,b,f))
map3(a, b, c, f::Callable) = flatten(mapvec3(a,b,c,f))
map4(a, b, c, d, f::Callable) = flatten(mapvec4(a,b,c,d,f))
map5(a, b, c, d, e, f::Callable) = flatten(mapvec5(a,b,c,d,e,f))

mapi(a, f::Callable ) = map2(a, 1:len(a), f)
mapveci(a, f::Callable ) = mapvec2(a, 1:len(a), f)

import Base.map
map(a, f::Callable) = map(unstack(1:len(a)), i->f(at(a,i)))
map(a::AbstractString, f::Callable) = flatten(map(unstack(a),f))
function map{T,N}(a::AbstractArray{T,N}, f::Callable)
    isempty(a) && return []

    r1 = f(fst(a))
    r = arraylike(r1, len(a), a)

    setat!(r,1,r1)
    map_(f,a,r)
    return r
end

@inline function map_{Tin,Nin,Tout,Nout}(f::Callable,a::AbstractArray{Tin,Nin},r::Array{Tout,Nout})
    for i = 2:len(a)
        b = f(at(a,i))
        setat!(r,i,b)
    end
end

function map(a::Dict, f::Callable; kargs...)
    isempty(a) && return a
    makeentry(a::Void) = []
    makeentry(a::Union{Tuple,Pair}) = [a]
    makeentry{T<:Union{Tuple,Pair}}(a::Array{T}) = a
    makeentry(a) = error("FunctionalData: map(::Dict), got entry of type $(typeof(a)), not one of Void, Tuple{Symbol,Any}, Array{Tuple}")

    r = @p vec a | map f | map makeentry | flatten
    d = Dict()
    for x in r
        d[fst(x)] = snd(x)
    end
    d
end

mapkeys(a::Dict, f) = map(a, x -> (f(fst(x)),snd(x)))
mapvalues(a::Dict, f) = map(a, x ->(fst(x),f(snd(x))))
mapmap(a::Dict, f) = [f(v) for (k,v) in a]

function mapmap(a, f)
    isempty(a) && return Any[]
    g(x) = map(x,f)
    map(a, g)
end

function mapmapvec(a, f)
    isempty(a) && return Any[]
    g(x) = mapvec(x,f)
    map(a, g)
end

function map{T<:Real}(a::DenseArray{T,1},f::Callable)
    isempty(a) && return []
    r1 = f(fst(a))
    r = arraylike(r1, len(a), a)
    if isviewable(r)
        rv = view(r,1)
        rv[:] = r1
        next!(rv)
        for i = 2:len(a)
            rv[:] = f(a[i])
            next!(rv)
        end
        r
    else
        setat!(r, 1, r1)
        map_(f,a,r)
        r
    end
end

function map{T<:Real,N}(a::DenseArray{T,N},f::Callable)
    isempty(a) && return []
    v = view(a,1)
    r1 = f(v)
    r = arraylike(r1, len(a), a)
    if isviewable(r)
        rv = view(r,1)
        rv[:] = r1
        next!(v)
        next!(rv)
        for i = 2:len(a)
            rv[:] = f(v)
            next!(v)
            next!(rv)
        end
        r
    else
        setat!(r, 1, r1)
        map_(f,a,r)
        r
    end
end

function map2!{T<:Real,N}(a::DenseArray{T,N}, f1::Callable, f2::Callable)
    isempty(a) && return []
    v = view(a,1)
    r1 = f1(v)
    r = arraylike(r1, len(a), a)
    rv = view(r,1)
    rv[:] = r1
    map2!(view(a, 2:len(a)), view(r, 2:len(a)), f2)
    r
end

function map2!{T<:Real,N,T2<:Real,M}(a::DenseArray{T,N}, r::DenseArray{T2,M}, f::Callable)
    isempty(a) && return []
    v = view(a,1)
    rv = view(r,1)
    for i = 1:len(a)
        f(rv, v)
        next!(v)
        next!(rv)
    end
    r
end

import Base.map!
function map!r{T<:Real,N}(a::DenseArray{T,N},f::Callable)
    isempty(a) && return a
    v = view(a,1)
    for i = 1:len(a)
        v[:] = f(v)
        next!(v)
    end
    a
end

function map!{T<:Real,N}(a::DenseArray{T,N},f::Callable)
    isempty(a) && return a
    v = view(a,1)
    for i = 1:len(a)
        f(v)
        next!(v)
    end
    a
end

work(a,f::Callable) = for i in 1:len(a) f(at(a,i)) end
work(a::Dict,f::Callable) = map(vec(a),(k,v)->(f;nothing))
function work{T<:Real,N}(a::DenseArray{T,N},f::Callable)
    len(a)==0 && return
    v = view(a,1)
    for i = 1:len(a)
        f(v)
        next!(v)
    end
end
lwork(a, f) = (g(x) = (f(x);nothing); lmap(a, g); nothing)
pwork(a, f) = (g(x) = (f(x);nothing); pmap(a, g); nothing)
hwork(a, f) = (g(x) = (f(x);nothing); hmap(a, g); nothing)
shwork(a, f) = (g(x) = (f(x);nothing); shmap(a, g); nothing)
workwork(a, f) = (g(x) = (f(x);nothing); mapmap(a, g); nothing)

work2(a, b, f::Callable) = [(f(at(a,i),at(b,i)); nothing) for i in 1:len(a)]
work3(a, b, c, f::Callable) = [(f(at(a,i),at(b,i),at(c,i)); nothing) for i in 1:len(a)]
work4(a, b, c, d, f::Callable) = [(f(at(a,i),at(b,i),at(c,i),at(d,i)); nothing) for i in 1:len(a)]
work5(a, b, c, d, e_, f::Callable) = [(f(at(a,i),at(b,i),at(c,i),at(d,i),at(e_,i)); nothing) for i in 1:len(a)]
worki(a, f) = mapi(a,(x,i)->(f(x,i);nothing))

share{T<:Real,N}(a::DenseArray{T,N}) = convert(SharedArray, a)
share(a::SharedArray) = a
unshare(a::SharedArray) = sdata(a)

function loopovershared(r::View, a::View, f::Callable)
    v = view(a,1)
    rv = view(r,1)
    for i = 1:len(a)
        rv[:] = f(v)
        next!(v)
        next!(rv)
    end
end

function loopovershared!(a::View, f::Callable)
    v = view(a, 1)
    for i = 1:len(a)
        f(v)
        next!(v)
    end
end

function loopovershared!r(a::View, f::Callable)
    v = view(a, 1)
    for i = 1:len(a)
        v[:] = f(v)
        next!(v)
    end
end

function loopovershared2!(r::View, a::View, f::Callable)
    v = view(a,1)
    rv = view(r,1)
    for i = 1:len(a)
        f(rv, v)
        next!(v)
        next!(rv)
    end
end

function shsetup(a::SharedArray; withfirst = false)
    pids = procs(a)
    if length(pids) > 1
        pids = pids[pids .!= 1]
    end
    ind = (withfirst ? 1 : 2):len(a)
    n = min(length(pids), len(ind))
    inds = partition(ind, n)
    pids, inds, n
end

shmap{T<:Real,N}(a::DenseArray{T,N}, f::Callable) = shmap(share(a), f)
function shmap{T<:Real,N}(a::SharedArray{T,N}, f::Callable)
    pids, inds, n = shsetup(a)

    r1 = f(view(a,1))
    r = sharraylike(r1, len(a))
    rv = view(r,1)
    rv[:] = r1

    @sync for i in 1:n
        @spawnat pids[i] loopovershared(view(r, inds[i]), view(a, inds[i]), f)
    end
    r
end

shmap!{T<:Real,N}(a::DenseArray{T,N}, f::Callable) = shmap!(share(a), f)
function shmap!{T<:Real,N}(a::SharedArray{T,N}, f::Callable)
    pids, inds, n = shsetup(a; withfirst = true)

    @sync for i in 1:n
        @spawnat pids[i] loopovershared!(view(a, inds[i]), f)
    end
    a
end

shmap!r{T<:Real,N}(a::DenseArray{T,N}, f::Callable) = shmap!r(share(a), f)
function shmap!r{T<:Real,N}(a::SharedArray{T,N}, f::Callable)
    pids, inds, n = shsetup(a; withfirst = true)

    @sync for i in 1:n
        @spawnat pids[i] loopovershared!r(view(a, inds[i]), f)
    end
    a
end

shmap2!{T<:Real, N}(a::DenseArray{T,N}, f1::Callable, f2::Callable) = shmap2!(share(a), f1, f2)
function shmap2!{T<:Real, N}(a::SharedArray{T,N}, f1::Callable, f2::Callable)
    r1 = f1(view(a,1))
    r = sharraylike(r1, len(a))
    rv = view(r,1)
    rv[:] = r1
    shmap2!(a, r, f2, withfirst = false)
    r
end

shmap2!{T<:Real,N, T2<:Real, M}(a::DenseArray{T,N}, r::SharedArray{T2,M}, f::Callable) = shmap2!(share(a), r, f)
function shmap2!{T<:Real,N, T2<:Real, M}(a::SharedArray{T,N}, r::SharedArray{T2,M}, f::Callable; withfirst = true)
    pids, inds, n = shsetup(a, withfirst = withfirst)
    @sync for i in 1:n
        @spawnat pids[i] loopovershared2!(view(r, inds[i]), view(a, inds[i]), f)
    end
    r
end

###############################
#  pmap

import Base.pmap
mapper(a, f) = map(a,f)
mappervec(a, f) = mapvec(a,f)
mapper!(a, f) = map!(a,f)
mapper!r(a, f) = map!r(a,f)
mapper2!(a, f1::Callable, f2) = map2!(a,f1,f2)
mapper2!(a, r, f) = map2!(a,r,f)
pmap(a, f::Callable; kargs...) = pmap_internal(mapper, a, f; kargs...)
pmapvec(a, f::Callable; kargs...) = pmap_internal(mappervec, a, f; kargs...)
pmap!(a, f; kargs...) = pmap_internal(mapper!, a, f; kargs...)
pmap!r(a, f; kargs...) = pmap_internal(mapper!r, a, f; kargs...)
pmap2!r(a, f1::Callable, f2::Callable; kargs...) = pmap_internal2!(mapper2!, a, f1, f2; kargs...)
pmap2!r(a, r, f::Callable; kargs...) = pmap_internal2!(mapper2!, a, r, f; kargs...)

workerpool(pids) =  VERSION < v"0.5-" ? pids : WorkerPool(pids)

function pmapsetup(a; pids = workers())
    if length(pids) > 1
        pids = pids[pids .!= 1]
    end
    n = min(length(pids)*10, len(a))
    inds = partition(1:len(a), n)
    pids, inds, n
end

function pmapparts(a, inds, n)
    if isa(a, DenseArray) && eltype(a)<:Real
        parts = [view(a, inds[i]) for i in 1:n]
    else
        parts = [part(a, inds[i]) for i in 1:n]
    end
end

function pmap_exec(g, a; nworkers = typemax(Int), vec = false, kargs...)
    pids, inds, n = pmapsetup(a; kargs...)
    parts = pmapparts(a, inds, n)
    pids = workerpool(take(pids, nworkers))
    if VERSION < v"0.5-"
        r = Base.pmap(x->(yield();g(x)), parts, pids = pids)
    else
        r = Base.pmap(pids, x->(yield();g(x)), parts)
    end
    for x in r
        isa(x,RemoteException) && rethrow(x)
    end
    flatten(r)
end

function pmap_internal(mapf::Callable, a, f::Callable; pids = workers(), kargs...)
    g(a) = mapf(a,f)
    pmap_exec(g, a; pids = pids, kargs...)
end

function pmap_internal2!(mapf::Callable, a, f1::Callable, f2::Callable; pids = workers(), kargs...)
    g(a) = mapf(a,f1,f2)
    pmap_exec(g, a; pids = pids, kargs...)
end

function pmap_internal2!(mapf::Callable, a, r, f::Callable; nworkers = typemax(Int), kargs...)
    g(a) = mapf(fst(a), snd(a), f)
    pids, inds, n = pmapsetup(a; kargs...)
    partsa = pmapparts(a, inds, n)
    partsr = pmapparts(r, inds, n)
    parts = zip(partsa, partsr)
    pids = workerpool(take(pids,nworkers))
    @show pids
    r = Base.pmap(g, parts, pids = pids)
    flatten(r)
end

export localworkers
localworkers(pid = myid()) = (r = sort(procs(pid)); len(r) > 1 ? r[2:end] : r)

export hostpids
hostpids() = @p map unstack(workers()) procs | map sort | unique | map x->x[1]==1 && len(x)>1 ? x[2] : x[1]

export hmap, hmap!, hmap!r, hmap2!r
hmap(a, f; kargs...) = pmap_internal(mapper, a, f; pids = hostpids(), kargs...)
hmap!(a, f; kargs...) = pmap_internal(mapper!, a, f; pids = hostpids(), kargs...)
hmap!r(a, f; kargs...) = pmap_internal(mapper!r, a, f; pids = hostpids(), kargs...)
hmap2!r(a, f1::Callable, f2::Callable; kargs...) = pmap_internal2!(mapper2!, a, f1, f2; pids = hostpids(), kargs...)
hmap2!r(a, r, f::Callable; kargs...) = pmap_internal2!(mapper2!, a, r, f; pids = hostpids(), kargs...)

lmap(a, f; kargs...) = pmap_internal(mapper, a, f; pids = localworkers(), kargs...)
lmapvec(a, f; kargs...) = pmap_internal(mappervec, a, f; pids = localworkers(), kargs...)
lmap!(a, f; kargs...) = pmap_internal(mapper!, a, f; pids = localworkers(), kargs...)
lmap!r(a, f; kargs...) = pmap_internal(mapper!r, a, f; pids = localworkers(), kargs...)
lmap2!r(a, f1::Callable, f2::Callable; kargs...) = pmap_internal2!(mapper2!, a, f1, f2; pids = localworkers(), kargs...)
lmap2!r(a, r, f::Callable; kargs...) = pmap_internal2!(mapper2!, a, r, f; pids = localworkers(), kargs...)

amapvec(a,f; kargs...) = amap(a,f; n = 10, mapper = mapvec)
function amap(a,f; n = 10, mapper = map)
    n = min(len(a),n)
    a = @p partition a n
    r = Array{Any}(n)
    @sync for i in 1:n
        @async r[i] = mapper(a[i], f)
    end
    flatten(r)
end
amap2(a,b,f; kargs...) = amap(czip(unstack(a),unstack(b)), x->f(fst(x),snd(x)); kargs...)
amapvec2(a,b,f; kargs...) = amapvec(czip(unstack(a),unstack(b)), x->f(fst(x),snd(x)); kargs...)

table(a...; kargs...) = table_internal(map, a[end], a[1:end-1]...; flat = true, kargs...)
ptable(a...; kargs...) = table_internal(pmap, a[end], a[1:end-1]...; flat = true, kargs...)
ltable(a...; kargs...) = table_internal(lmap, a[end], a[1:end-1]...; flat = true, kargs...)
htable(a...; kargs...) = table_internal(hmap, a[end], a[1:end-1]...; flat = true, kargs...)

table(f::Callable, a...; kargs...) = table_internal(map, f, a...; flat = true, kargs...)
ptable(f::Callable, a...; kargs...) = table_internal(pmap, f, a...; flat = true, kargs...)
ltable(f::Callable, a...; kargs...) = table_internal(lmap, f, a...; flat = true, kargs...)
htable(f::Callable, a...; kargs...) = table_internal(hmap, f, a...; flat = true, kargs...)

tableany(a...; kargs...) = table_internal(map, a[end], a[1:end-1]...; flat = false, kargs...)
ptableany(a...; kargs...) = table_internal(pmap, a[end], a[1:end-1]...; flat = false, kargs...)
ltableany(a...; kargs...) = table_internal(lmap, a[end], a[1:end-1]...; flat = false, kargs...)
htableany(a...; kargs...) = table_internal(hmap, a[end], a[1:end-1]...; flat = false, kargs...)

tableany(f::Callable, a...; kargs...) = table_internal(map, f, a...; flat = false, kargs...)
ptableany(f::Callable, a...; kargs...) = table_internal(pmap, f, a...; flat = false, kargs...)
ltableany(f::Callable, a...; kargs...) = table_internal(lmap, f, a...; flat = false, kargs...)
htableany(f::Callable, a...; kargs...) = table_internal(hmap, f, a...; flat = false, kargs...)

function table_internal(mapf, f, args...; flat = true, kargs...)
    a = [isa(x,Range) ? collect(x) : x for x in args]
    s = @p col [len(x) for x in a]
    if any(s.==0)
        error("FunctionalData.xtable: empty arguments, lengths of the $(len(args)) arguments are $(vec(s))")
    end
    S = tuple(s...)
    getarg(sub) = [at(a[i],sub[i]) for i in 1:length(a)]
    args = [getarg(ind2sub(S,x)) for x in 1:prod(s)]
    g(x) = f(x...)
    r = mapf(args,g; kargs...)

    if flat
        if length(fst(r))==1
            news = s
        else
            s1 = siz(fst(r))
            s1 = s1[end] == 1 ? s1[1:end-1] : s1
            news = vcat(s1,s)
        end
        @p flatten r | reshape news
    else
        @p reshape r s
    end
end

areall(a, f::Callable) = !any(a, not*f)

function isany(a, f::Callable)
    for i in 1:len(a)
        f(at(a,i)) && return true
    end
    return false
end

tee(a,f) = (f(a);a)

import Base.*
*(f::Callable, g::Callable) = (a...) -> f(g(a...))

typed{N}(a::Array{Any,N}) = isempty(a) ? [] : convert(Array{typeof(a[1]),N}, a)
typed(a) = a

import Base.filter
filter(r::Regex, f2::Callable) = error()
filter(f1::Callable, f2::Callable) = error()
function filter(a, f::Callable)
    ind = vec(typed(map(a,f)))
    isempty(ind) ? [] : part(a, ind)
end

import Base.select
select(a, f::Callable) = filter(a, f)
reject(a, f::Callable) = select(a, not*f)

function groupdict(a,f::Function = id)
    d = Dict()
    inds = @p map a f
    for i in 1:len(a)
        ind = @p at inds i
        d[ind] = push!(get(d,ind,[]), at(a,i))
    end
    mapvalues(d,flatten)
end

function groupby(a,f = id)
    r = groupdict(a,f)
    ks = collect(keys(r))
    try
        ks = sort(ks)
    end
    values(r,ks)
end

import Base.apply
apply(f::Callable, f2::Callable) = error("undefined")
apply(f::Callable, args...) = f(args...)
apply(args, f::Callable) = f(args)

function fold(a, f::Callable)
    if isempty(a)
        return a
    end
    if len(a) == 1
        return fst(a)
    end
    r = f(fst(a), snd(a))
    for i = 3:len(a)
        r = f(r,at(a,i))
    end
    r
end


text:
xnodes(ϕ) = repeat(reshape(ϕ.grid.xC, ϕ.grid.Nx, 1, 1), 1, ϕ.grid.Ny, ϕ.grid.Nz)
ynodes(ϕ) = repeat(reshape(ϕ.grid.yC, 1, ϕ.grid.Ny, 1), ϕ.grid.Nx, 1, ϕ.grid.Nz)
znodes(ϕ) = repeat(reshape(ϕ.grid.zC, 1, 1, ϕ.grid.Nz), ϕ.grid.Nx, ϕ.grid.Ny, 1)

xnodes(ϕ::FaceFieldX) = repeat(reshape(ϕ.grid.xF[1:end-1], ϕ.grid.Nx, 1, 1), 1, ϕ.grid.Ny, ϕ.grid.Nz)
ynodes(ϕ::FaceFieldY) = repeat(reshape(ϕ.grid.yF[1:end-1], 1, ϕ.grid.Ny, 1), ϕ.grid.Nx, 1, ϕ.grid.Nz)
znodes(ϕ::FaceFieldZ) = repeat(reshape(ϕ.grid.zF[2:end],   1, 1, ϕ.grid.Nz), ϕ.grid.Nx, ϕ.grid.Ny, 1)

zerofunk(args...) = 0

function set_ic!(model; u=zerofunk, v=zerofunk, w=zerofunk, T=zerofunk, S=zerofunk)
    model.velocities.u.data .= u.(xnodes(model.velocities.u), ynodes(model.velocities.u), znodes(model.velocities.u))
    model.velocities.v.data .= v.(xnodes(model.velocities.v), ynodes(model.velocities.v), znodes(model.velocities.v))
    model.velocities.w.data .= w.(xnodes(model.velocities.w), ynodes(model.velocities.w), znodes(model.velocities.w))
    model.tracers.T.data    .= T.(xnodes(model.tracers.T),    ynodes(model.tracers.T),    znodes(model.tracers.T))
    model.tracers.S.data    .= S.(xnodes(model.tracers.S),    ynodes(model.tracers.S),    znodes(model.tracers.S))
    return nothing
end

function u_relative_error(model, u)
    u_num = model.velocities.u

    u_ans = FaceFieldX(u.(
        xnodes(u_num), ynodes(u_num), znodes(u_num), model.clock.time),
        model.grid)

    return mean(
        (u_num.data .- u_ans.data).^2 ) / mean(u_ans.data.^2)
end

function w_relative_error(model, w)
    w_num = model.velocities.w

    w_ans = FaceFieldZ(
        w.(xnodes(w_num), ynodes(w_num), znodes(w_num), model.clock.time),
        model.grid)

    return mean(
        (w_num.data .- w_ans.data).^2 ) / mean(w_ans.data.^2)
end

function T_relative_error(model, T)
    T_num = model.tracers.T

    T_ans = FaceFieldZ(
        T.(xnodes(T_num), ynodes(T_num), znodes(T_num), model.clock.time),
        model.grid)

    return mean(
        (T_num.data .- T_ans.data).^2 ) / mean(T_ans.data.^2)
end

function test_diffusion_simple(fld)
    Nx, Ny, Nz = 1, 1, 16
    Lx, Ly, Lz = 1, 1, 1
    κ = 1
    eos = LinearEquationOfState(βS=0, βT=0)

    model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=κ, κ=κ, eos=eos)

    if fld ∈ (:u, :v, :w)
        field = getfield(model.velocities, fld)
    else
        field = getfield(model.tracers, fld)
    end

    value = π
    field.data .= value

    Δt = 0.01 # time-step much less than diffusion time-scale
    Nt = 10
    time_step!(model, Nt, Δt)

    !any(@. !isapprox(value, field.data))
end

function test_diffusion_budget(field_name)
    Nx, Ny, Nz = 1, 1, 16
    Lx, Ly, Lz = 1, 1, 1
    κ = 1
    eos = LinearEquationOfState(βS=0, βT=0)

    model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=κ, κ=κ, eos=eos)

    if field_name ∈ (:u, :v, :w)
        field = getfield(model.velocities, field_name)
    else
        field = getfield(model.tracers, field_name)
    end

    half_Nz = floor(Int, Nz/2)
    @views @. field.data[:, :, 1:half_Nz]  = -1
    @views @. field.data[:, :, half_Nz:end] =  1

    mean_init = mean(field.data)
    τκ = Lz^2 / κ # diffusion time-scale
    Δt = 0.0001 * τκ # time-step much less than diffusion time-scale
    Nt = 100

    time_step!(model, Nt, Δt)
    isapprox(mean_init, mean(field.data))
end

function test_diffusion_cosine(fld)
    Nx, Ny, Nz = 1, 1, 128
    Lx, Ly, Lz = 1, 1, π/2
    κ = 1
    eos = LinearEquationOfState(βS=0, βT=0)

    model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=κ, κ=κ, eos=eos)

    if fld == :w
        throw("There are no boundary condition tests for w yet.")
    elseif fld ∈ (:u, :v)
        field = getfield(model.velocities, fld)
    else
        field = getfield(model.tracers, fld)
    end

    zC = model.grid.zC
    m = 2
    @views @. field.data[1, 1, :] = cos.(m*zC)

    diffusing_cosine(κ, m, z, t) = exp(-κ*m^2*t) * cos(m*z)

    τκ = Lz^2 / κ # diffusion time-scale
    Δt = 1e-6 * τκ # time-step much less than diffusion time-scale
    Nt = 100

    time_step!(model, Nt, Δt)

    field_numerical = dropdims(field.data, dims=(1, 2))

    !any(@. !isapprox(field_numerical, diffusing_cosine(κ, m, zC, model.clock.time), atol=1e-6, rtol=1e-6))
end

function internal_wave_test(; N=128, Nt=10)

    # Internal wave parameters
     ν = κ = 1e-9
     L = 2π
    z₀ = -L/3
     δ = L/20
    a₀ = 1e-3
     m = 16
     k = 1
     f = 0.2
     ℕ = 1.0
     σ = sqrt( (ℕ^2*k^2 + f^2*m^2) / (k^2 + m^2) )

    # Numerical parameters
     N = 128
    Δt = 0.01 * 1/σ

    cᵍ = m * σ / (k^2 + m^2) * (f^2/σ^2 - 1)
     U = a₀ * k * σ   / (σ^2 - f^2)
     V = a₀ * k * f   / (σ^2 - f^2)
     W = a₀ * m * σ   / (σ^2 - ℕ^2)
     Θ = a₀ * m * ℕ^2 / (σ^2 - ℕ^2)

    a(x, y, z, t) = exp( -(z - cᵍ*t - z₀)^2 / (2*δ)^2 )

    u(x, y, z, t) =           a(x, y, z, t) * U * cos(k*x + m*z - σ*t)
    v(x, y, z, t) =           a(x, y, z, t) * V * sin(k*x + m*z - σ*t)
    w(x, y, z, t) =           a(x, y, z, t) * W * cos(k*x + m*z - σ*t)
    T(x, y, z, t) = ℕ^2 * z + a(x, y, z, t) * Θ * sin(k*x + m*z - σ*t)

    u₀(x, y, z) = u(x, y, z, 0)
    v₀(x, y, z) = v(x, y, z, 0)
    w₀(x, y, z) = w(x, y, z, 0)
    T₀(x, y, z) = T(x, y, z, 0)

    # Create a model where temperature = buoyancy.
    model = Model(N=(N, 1, N), L=(L, L, L), ν=ν, κ=κ,
                    eos=LinearEquationOfState(βT=1.),
                    constants=PlanetaryConstants(f=f, g=1.))

    set_ic!(model, u=u₀, v=v₀, w=w₀, T=T₀)

    time_step!(model, Nt, Δt)

    # Tolerance was found by trial and error...
    u_relative_error(model, u) < 1e-4
end

function passive_tracer_advection_test(; N=128, κ=1e-12, Nt=100)
    L, U, V = 1.0, 0.5, 0.8
    δ, x₀, y₀ = L/15, L/2, L/2

    Δt = 0.05 * L/N / sqrt(U^2 + V^2)

    T(x, y, z, t) = exp( -((x - U*t - x₀)^2 + (y - V*t - y₀)^2) / (2*δ^2) )
    u₀(x, y, z) = U
    v₀(x, y, z) = V
    T₀(x, y, z) = T(x, y, z, 0)

    model = Model(N=(N, N, 2), L=(L, L, L), ν=κ, κ=κ)

    set_ic!(model, u=u₀, v=v₀, T=T₀)
    time_step!(model, Nt, Δt)

    # Error tolerance is a bit arbitrary
    return T_relative_error(model, T) < 1e-4
end

ext: jl


text:
# Parly derived from latex_symbols.jl, which is a part of Julia
# License is MIT: http://julialang.org/license

const greek_letters =
    ("Alpha"   => "A",
     "Beta"    => "B",
     "Gamma"   => "G",
     "Delta"   => "D",
     "Epsilon" => "E",
     "Zeta"    => "Z",
     "Eta"     => "H",
     "Theta"   => "J",
     "Iota"    => "I",
     "Kappa"   => "K",
     "Lambda"  => "L",
     "Mu"      => "M",
     "Nu"      => "N",
     "Xi"      => "X",
     "Omicron" => "U",
     "Pi"      => "P",
     "Rho"     => "R",
     "Sigma"   => "S",
     "Tau"     => "T",
     "Upsilon" => "Y",
     "Phi"     => "F",
     "Chi"     => "C",
     "Psi"     => "W",
     "Omega"   => "O",
     "alpha"   => "a",
     "beta"    => "b",
     "gamma"   => "g",
     "delta"   => "d",
     "epsilon" => "e",
     "zeta"    => "z",
     "eta"     => "h",
     "theta"   => "j",
     "iota"    => "i",
     "kappa"   => "k",
     "lambda"  => "l",
     "mu"      => "m",
     "nu"      => "n",
     "xi"      => "x",
     "omicron" => "u",
     "pi"      => "p",
     "rho"     => "r",
     "sigma"   => "s",
     "tau"     => "t",
     "upsilon" => "y",
     "phi"     => "f",
     "chi"     => "c",
     "psi"     => "w",
     "omega"   => "o",
)

const var_greek =
    ("varTheta"   => "J",
     "nabla"      => "n",
     "partial"    => "d", # partial differential
     "varepsilon" => "e",
     "varsigma"   => "s",
     "vartheta"   => "j",
     "varkappa"   => "k",
     "varphi"     => "f",
     "varrho"     => "r",
     "varpi"      => "p"
     )

const digits = (
    "zero"  => "0",
    "one"   => "1",
    "two"   => "2",
    "three" => "3",
    "four"  => "4",
    "five"  => "5",
    "six"   => "6",
    "seven" => "7",
    "eight" => "8",
    "nine"  => "9"
    )

const remove_lead_char = "AEt"
const remove_prefix = ("APL", "Elz", "Elx", "El", "textascii", "text")

const replace_lead_char = "Bm"
const replace_prefix =
    (("Bbb",       "d",   "bb"),	# double-struck or blackboard
     ("mbfsans",   "sb",  "bsans"),     # bold sans-serif
     ("mbfscr",    "cb",  "bscr"),	# bold cursive script
     ("mbffrak",   "fb",  "bfrak"),     # bold fraktur
     ("mbfitsans", "sib", "bisans"),    # bold italic sans-serif
     ("mbfit",     "ib",  "bi"),        # bold italic
     ("mbf",       "b",   "bf"),        # bold
     ("mfrak",     "f",   "frak"),      # fraktur
     ("mitsans",   "si",  "isans"),     # italic sans-serif
     ("mitBbb",    "di",  "bbi"),       # italic double-struck (or blackboard)
     ("mit",       "i",   "it"),        # italic
     ("msans",     "s",   "sans"),      # sans-serif
     ("mscr",      "c",   "scr"),       # cursive script
     ("mtt",       "t",   "tt")         # teletype (monospaced)
)

const remove_name = ("Elxsqcup", "Elxuplus", "ElOr", "textTheta", "Elzbar")

const replace_name = (
    "textasciiacute"  => "textacute",
    "textasciibreve"  => "textbreve",
    "textasciimacron" => "highminus",
    "textphi"         => "ltphi",
    "Eulerconst"      => "eulermascheroni",
    "Hermaphrodite"   => "hermaphrodite",
    "Planckconst"     => "planck",
    )

const manual = [
    "cbrt"        => "\u221B", # synonym of \cuberoot
    "mars"        => "♂",      # synonym of \male
    "pprime"      => "″",      # synonym of \dprime
    "ppprime"     => "‴",      # synonym of \trprime
    "pppprime"    => "⁗",      # synonym of \qprime
    "backpprime"  => "‶",      # synonym of \backdprime
    "backppprime" => "‷",      # synonym of \backtrprime
    "emptyset"    => "∅",      # synonym of \varnothing
    "llbracket"   => "⟦",      # synonym of \lBrack
    "rrbracket"   => "⟧",      # synonym of \rBrack
    "xor"         => "⊻",      # synonym of \veebar
    "iff"         => "⟺",
    "implies"     => "⟹",
    "impliedby"   => "⟸",
    "to"          => "→",
    "euler"       => "ℯ",

    # Misc. Math and Physics
    "del"         => "∇",      # synonym of \nabla (combining character)
    "sout"        => "\u0336", # synonym of \Elzbar (from ulem package)
    "strike"      => "\u0336", # synonym of \Elzbar
    "zbar"        => "\u0336", # synonym of \Elzbar

    # Avoid getting "incorrect" synonym
    "imath"       => "ı",
    "jmath"       => "ȷ",
    "i_imath"     => "\U1d6a4",     # mathematical italic small dotless i
    "i_jmath"     => "\U1d6a5",     # mathematical italic small dotless j
    "hbar"        => "\u0127",      # ħ synonym of \Elzxh
    "AA"          => "\u00c5",      # Å
    "Upsilon"     => "\u03a5",      # Υ
    "setminus"    => "\u2216",      # ∖ synonym of \smallsetminus
    "ddot{i}"     => "\u00cf",      # is ddot{\imath} in unicode.xml
    "bigsetminus" => "\u29f5",      # add to allow access to standard setminus
    "circlearrowleft"  => "\u21ba", # ↺ synonym of acwopencirclearrow
    "circlearrowright" => "\u21bb", # ↻ synonym of cwopencirclearrow

    "ohm"         => "Ω",
    "leq"         => "≤",  # synonym of le
    "geq"         => "≥",  # synonym of ge
    "bbsemi"      => "⨟",
    "ith"         => "ℎ",  # mathematical italic small h (planck constant)
    "tricolon"    => "⁝",  # tricolon
    "join"        => "⨝",  # synonym of Join
]

# Vulgar fractions
const fractions = [
    "1/4"  => "¼", # vulgar fraction one quarter
    "1/2"  => "½", # vulgar fraction one half
    "3/4"  => "¾", # vulgar fraction three quarters
    "1/7"  => "⅐",# vulgar fraction one seventh
    "1/9"  => "⅑", # vulgar fraction one ninth
    "1/10" => "⅒", # vulgar fraction one tenth
    "1/3"  => "⅓", # vulgar fraction one third
    "2/3"  => "⅔", # vulgar fraction two thirds
    "1/5"  => "⅕", # vulgar fraction one fifth
    "2/5"  => "⅖", # vulgar fraction two fifths
    "3/5"  => "⅗", # vulgar fraction three fifths
    "4/5"  => "⅘", # vulgar fraction four fifths
    "1/6"  => "⅙", # vulgar fraction one sixth
    "5/6"  => "⅚", # vulgar fraction five sixths
    "1/8"  => "⅛", # vulgar fraction one eigth
    "3/8"  => "⅜", # vulgar fraction three eigths
    "5/8"  => "⅝", # vulgar fraction five eigths
    "7/8"  => "⅞", # vulgar fraction seventh eigths
    "1/"   => "⅟", # fraction numerator one
    "0/3"  => "↉", # vulgar fraction zero thirds
    "1/4"  => "¼", # vulgar fraction one quarter
]

const superscripts = [
    "^0" => "⁰",
    "^1" => "¹",
    "^2" => "²",
    "^3" => "³",
    "^4" => "⁴",
    "^5" => "⁵",
    "^6" => "⁶",
    "^7" => "⁷",
    "^8" => "⁸",
    "^9" => "⁹",
    "^+" => "⁺",
    "^-" => "⁻",
    "^=" => "⁼",
    "^(" => "⁽",
    "^)" => "⁾",
    "^a" => "ᵃ",
    "^b" => "ᵇ",
    "^c" => "ᶜ",
    "^d" => "ᵈ",
    "^e" => "ᵉ",
    "^f" => "ᶠ",
    "^g" => "ᵍ",
    "^h" => "ʰ",
    "^i" => "ⁱ",
    "^j" => "ʲ",
    "^k" => "ᵏ",
    "^l" => "ˡ",
    "^m" => "ᵐ",
    "^n" => "ⁿ",
    "^o" => "ᵒ",
    "^p" => "ᵖ",
    "^r" => "ʳ",
    "^s" => "ˢ",
    "^t" => "ᵗ",
    "^u" => "ᵘ",
    "^v" => "ᵛ",
    "^w" => "ʷ",
    "^x" => "ˣ",
    "^y" => "ʸ",
    "^z" => "ᶻ",
    "^A" => "ᴬ",
    "^B" => "ᴮ",
    "^D" => "ᴰ",
    "^E" => "ᴱ",
    "^G" => "ᴳ",
    "^H" => "ᴴ",
    "^I" => "ᴵ",
    "^J" => "ᴶ",
    "^K" => "ᴷ",
    "^L" => "ᴸ",
    "^M" => "ᴹ",
    "^N" => "ᴺ",
    "^O" => "ᴼ",
    "^P" => "ᴾ",
    "^R" => "ᴿ",
    "^T" => "ᵀ",
    "^U" => "ᵁ",
    "^V" => "ⱽ",
    "^W" => "ᵂ",
    "^alpha" => "ᵅ",
    "^beta" => "ᵝ",
    "^gamma" => "ᵞ",
    "^delta" => "ᵟ",
    "^epsilon" => "ᵋ",
    "^theta" => "ᶿ",
    "^iota" => "ᶥ",
    "^phi" => "ᵠ",
    "^chi" => "ᵡ",
    "^Phi" => "ᶲ",

    "^uparrow" => "ꜛ",
    "^downarrow" => "ꜜ",
    "^!" => "ꜝ",
]

const subscripts = [
    "_0" => "₀",
    "_1" => "₁",
    "_2" => "₂",
    "_3" => "₃",
    "_4" => "₄",
    "_5" => "₅",
    "_6" => "₆",
    "_7" => "₇",
    "_8" => "₈",
    "_9" => "₉",
    "_+" => "₊",
    "_-" => "₋",
    "_=" => "₌",
    "_(" => "₍",
    "_)" => "₎",
    "_a" => "ₐ",
    "_e" => "ₑ",
    "_h" => "ₕ",
    "_i" => "ᵢ",
    "_j" => "ⱼ",
    "_k" => "ₖ",
    "_l" => "ₗ",
    "_m" => "ₘ",
    "_n" => "ₙ",
    "_o" => "ₒ",
    "_p" => "ₚ",
    "_r" => "ᵣ",
    "_s" => "ₛ",
    "_t" => "ₜ",
    "_u" => "ᵤ",
    "_v" => "ᵥ",
    "_x" => "ₓ",
    "_schwa" => "ₔ",
    "_beta" => "ᵦ",
    "_gamma" => "ᵧ",
    "_rho" => "ᵨ",
    "_phi" => "ᵩ",
    "_chi" => "ᵪ"
]

const mansym = [manual, fractions, superscripts, subscripts]
const mantyp = ["manual", "fractions", "superscripts", "subscripts"]


text:
# Overloadable skinner for column names. This provides a generic method for naming the
# fields of a specific tuple instance. This data type provides an iterable interface to
# the field names
immutable DataView{N}
	tuple::NTuple{N}
end

# Generate the column name skins for the column vector store. Each column name-index pair is
# stored in a pair of constant functions, essentially forming a look up table, or
# enumeration. The functions are typed by the singleton types encoding the object id of the
# column vector tuple store, and either the symbol of the column name, or the index of the
# column.
function _DataView(view::DataView, columnname::Symbol, columnindex::Int64)
	eval(:(DataView(::Type{Val{$(view)}}, ::Type{Val{$(QuoteNode(columnname))}}) = $(columnindex)))
	eval(:(DataView(::Type{Val{$(view)}}, ::Type{Val{$(columnindex)}}) = $(QuoteNode(columnname))))
end

# Courtesy method to ensure the names are valid
Base.isvalid{N}(V::DataView{N}) =
	length(unique(getcolumnnames(V))) == N &&
	all([(getcolumnindex(V, getcolumnname(V, c)) == c)::Bool for c = 1:N])

# Iterable methods
Base.start(::DataView) = 1
Base.next(V::DataView, r::Int64) = (V[r], r + 1)
Base.done(V::DataView, r::Int64) = r > length(V)
Base.eltype(::Type{DataView}) = Symbol
Base.length{N}(V::DataView{N}) = N

# Collection methods
Base.isempty(D::DataView) = length(D) < 1
Base.endof(D::DataView) = length(D)

# Indexable methods
Base.setindex!(V::DataView, n::Symbol, c::Int64) = _DataView(V, n, c)
Base.getindex(V::DataView, c::Int64) = DataView(Val{V}, Val{c})
Base.getindex(V::DataView, n::Symbol) = DataView(Val{V}, Val{n})

# Array methods
Base.ndims(::DataView) = 1
Base.size(V::DataView) = (length(V))
Base.linearindexing(::Type{DataView}) = Base.LinearFast()


text:
# Goal: Test the correctness of state initialization
using ProbNumDiffEq
using OrdinaryDiffEq
using LinearAlgebra
using Test

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems;
importodeproblems();
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_fitzhughnagumo, prob_ode_pleiades

d = 2
q = 6
D = d * (q + 1)

a, b = 1.1, -0.5
f(u, p, t) = [a * u[1], b * u[2]]
u0 = [0.1, 1.0]
tspan = (0.0, 5.0)
t0, T = tspan
prob = ODEProblem(f, u0, tspan)
p = prob.p

# True Solutions and derivatives
u(t) = [a^0 * u0[1] * exp(a * t), u0[2] * exp(b * t)]
du(t) = [a^1 * u0[1] * exp(a * t), b * u0[2] * exp(b * t)]
ddu(t) = [a^2 * u0[1] * exp(a * t), (b)^2 * u0[2] * exp(b * t)]
dddu(t) = [a^3 * u0[1] * exp(a * t), (b)^3 * u0[2] * exp(b * t)]
ddddu(t) = [a^4 * u0[1] * exp(a * t), (b)^4 * u0[2] * exp(b * t)]
dddddu(t) = [a^5 * u0[1] * exp(a * t), (b)^5 * u0[2] * exp(b * t)]
ddddddu(t) = [a^6 * u0[1] * exp(a * t), (b)^6 * u0[2] * exp(b * t)]
true_init_states = [u(t0); du(t0); ddu(t0); dddu(t0); ddddu(t0); dddddu(t0); ddddddu(t0)]

@testset "Taylormode initialization" begin
    @testset "OOP" begin
        dfs = ProbNumDiffEq.taylormode_get_derivatives(
            prob.u0,
            prob.f,
            prob.p,
            prob.tspan[1],
            q,
        )
        @test length(dfs) == q + 1
        @test true_init_states ≈ vcat(dfs...)
    end

    @testset "IIP" begin
        f!(du, u, p, t) = (du .= f(u, p, t))
        prob = ODEProblem(f!, u0, tspan)

        dfs = ProbNumDiffEq.taylormode_get_derivatives(
            prob.u0,
            prob.f,
            prob.p,
            prob.tspan[1],
            q,
        )
        @test length(dfs) == q + 1
        @test true_init_states ≈ vcat(dfs...)
    end
end

@testset "Compare TaylorModeInit and ClassicSolverInit" begin
    # This has not worked before with the Taylor-Mode init!
    # The high dimensions made the runtimes explode

    prob = prob_ode_fitzhughnagumo
    d = length(prob.u0)

    integ1 = init(prob, EK0(order=8))
    tm_init = integ1.cache.x.μ
    Proj1 = integ1.cache.Proj

    @testset "Order $o" for o in (1, 2, 3, 4, 5)
        integ2 = init(prob, EK0(order=o, initialization=ClassicSolverInit(init_on_du=true)))
        rk_init = integ2.cache.x.μ
        Proj2 = integ2.cache.Proj

        # @info "how do things look?" tm_init rk_init
        # error()

        # These are fit via the initial values + autodiff, and should be good
        @test Proj1(0) * tm_init ≈ Proj2(0) * rk_init
        @test Proj1(1) * tm_init ≈ Proj2(1) * rk_init
        (o > 1) && @test Proj1(2) * tm_init ≈ Proj2(2) * rk_init

        # # Test if the others are correct, up to order 5
        (5 > o > 2) && @test Proj1(3) * tm_init ≈ Proj2(3) * rk_init rtol = 1e-1
        (5 > o > 3) && @test Proj1(4) * tm_init ≈ Proj2(4) * rk_init rtol = 1e-1
        # (o == 5) && @test Proj1(3) * tm_init ≈ Proj2(3) * rk_init rtol=6e-1
        # (o == 5) && @test Proj1(4) * tm_init ≈ Proj2(4) * rk_init rtol=2e-2
        # (o == 5) && @test Proj1(5) * tm_init ≈ Proj2(5) * rk_init rtol=8e-1
        # 8e-1 is only because of macOS, on linux this test should be much tighter

        # Test if the covariance covers the true error
        for i in 3:o
            _rk = Proj2(i) * rk_init
            _tm = Proj1(i) * tm_init
            err = _rk .- _tm
            C = ProbNumDiffEq.X_A_Xt(integ2.cache.x.Σ, Proj2(i))
            @assert isdiag(C)
            whitened_err = err ./ sqrt.(diag(C.mat))
            @test all(abs.(whitened_err) .< 4e-1)
        end
    end
end

@testset "ClassicSolverInit enables solving Pleiades with order 5" begin
    # This has not worked before with the Taylor-Mode init, it simply took too long!
    prob = prob_ode_pleiades
    sol1 = solve(prob, Vern9(), abstol=1e-10, reltol=1e-10, save_everystep=false)
    sol2 = solve(
        prob,
        EK1(
            order=5,
            smooth=false,
            initialization=ClassicSolverInit(alg=Rosenbrock32(), init_on_du=true),
        ),
        abstol=1e-6,
        reltol=1e-6,
        save_everystep=false,
    )
    @test sol1.u[end] ≈ sol2.u[end] rtol = 1e-5
end


text:
struct AnimatedGIF{T, A<:AbstractArray} <: AbstractArray{T, 3}
    data::A
    fps::Int

    # NOTE: the default fps is chosen purely by experience and may be changed in the future
    function AnimatedGIF{T, A}(data::A; fps::Real=min(10, size(data, 3)÷2)) where {T, A<:AbstractArray}
        if fps < 1
            fps = 1
            @warn "FPS should be larger than 1" FPS=fps
        end

        new{T, A}(data, ceil(Int, fps))
    end
end
AnimatedGIF(data::AbstractArray; kwargs...) = AnimatedGIF{eltype(data), typeof(data)}(data; kwargs...)

# Display gif requires ImageMagick v1.2.0, which requires Julia 1.3
Base.showable(::MIME"image/gif", agif::AnimatedGIF) = true
function Base.show(io::IO, ::MIME"image/gif", agif::AnimatedGIF)
    FileIO.save(_format_stream(format"GIF", io), agif.data; fps=agif.fps)
end

Base.size(A::AnimatedGIF) = size(A.data)
Base.axes(A::AnimatedGIF) = axes(A.data)
Base.@propagate_inbounds Base.getindex(A::AnimatedGIF, inds::Vararg{Int}) = getindex(A.data, inds...)

"""
    gif(img; kwargs...)
    gif(frames; kwargs...)

Convert 3D `img` or 2D frame list `frames` to animated gif array.

!!! compat "ImageMagick 1.2"
    ImageMagick at least 1.2 is required to generate gif image. For ImageMagick 0.y.z versions,
    ImageMagick at least 0.7.7 also works.

# Arguments

- `img::AbstractArray{T, 3}`: the last dimension is the time axis.
- `frames`: vector of 2d arrays.

# Parameters

- `fps::Int`: frame per second.

# Examples:

You can use this to view a 3D image:

```julia
using ImageShow, TestImages
ImageShow.gif(testimage("mri-stack"))
```

Or use it to dynamically check how different parameters change the function behavior:

```julia
using TestImages, ImageShow, ImageTransformations
img = testimage("cameraman")
ImageShow.gif([imrotate(img, θ, axes(img)) for θ in -π/4:π/16:π/4]; fps=3)
```

See also `mosaic`, provided by `MosaicViews`/`ImageCore`, for a 2d alternative of `gif`.
"""
gif(img::AbstractArray{<:Colorant, 3}; kwargs...) = AnimatedGIF(img; kwargs...)
gif(img::AbstractArray{<:Real, 3}; kwargs...) = AnimatedGIF(of_eltype(Gray, img); kwargs...)
function gif(frames::AbstractVector{<:AbstractMatrix}; kwargs...)
    FT = eltype(eltype(frames))
    fillvalue = zero(FT) # require ColorVectorSpace for RGB and Gray

    # FIXME: sym_paddedviews+StackView makes it unable to infer the types
    # TODO: sym_paddedviews(zero(FT), frames...) would just break whatever the fancy lazy
    #       feature `frames` has and thus slows the code
    sz = size(first(frames))
    frames = all(A->sz == size(A), frames) ? frames : sym_paddedviews(zero(FT), frames...)
    stacked_frames = StackView{FT}(frames, Val(3))
    gif(stacked_frames; kwargs...)
end

# Because any object can be callable (not just `Function`), we add this method to distinguish
# the `f` function version and throw a menaingful message.
# TODO: Maybe we could just support `gif(framestack...)`?
gif(::AbstractMatrix, ::AbstractMatrix...) = throw(ArgumentError("Do you mean `gif([A, Bs...])"))

"""
    gif(f, Xs; kwargs...)
    gif(f, Xs, Ys...; kwargs...)

A lazy version of `gif([f(X) for X in Xs]; kwargs...)` that allocates memory only when needed.

# Parameters

- `fps::Int`: frame per second.

# Examples

Rotate the image and see how things going:

```julia
using TestImages, ImageShow, ImageTransformations
img = testimage("cameraman")
ImageShow.gif(-π/4:π/16:π/4]; fps=3) do θ
    imrotate(img, θ, axes(img))
end
```

The following example is less meaningful, but it shows how multiple arguments are passed:

```julia
sizes = 16:4:64
values = range(0, stop=1, length=length(sizes))
gif = ImageShow.gif(values, sizes) do v, x
    fill(RGB(v, v, v), ntuple(_->x, 2)...)
end
```
"""
gif(f, arg1, args...; kwargs...) = gif(mappedarray(f, arg1, args...); kwargs...)
# MappedArrays are not efficient here https://github.com/JuliaArrays/MappedArrays.jl/issues/46
gif(frames::AbstractMappedArray; kwargs...) = gif(collect(frames); kwargs...)


text:
#=
   From Allen Downey
   https://twitter.com/AllenDowney/status/1063460117029535746
   """
   Today's Bayesian problem of the week: Suppose we visit a wild animal preserve where we
   know that the only animals are lions and tigers and bears, but we don't know how
   many of each there are.

   During the tour, we see 3 lions, 2 tigers, and 1 bear. Assuming that every animal had an equal
   chance to appear in our sample, estimate the prevalence of each species.

   What is the probability that the next animal we see is a bear?
   """

   Also see: https://towardsdatascience.com/estimating-probabilities-with-bayesian-modeling-in-python-7144be007815

   Cf ~/blog/lions_tigers_and_bears.blog
      ~/webppl/lions_tigers_and_bears.wppl

    From https://towardsdatascience.com/estimating-probabilities-with-bayesian-modeling-in-python-7144be007815
    """
    Species: lions    Prevalence: 44.44%.
    Species: tigers   Prevalence: 33.33%.
    Species: bears    Prevalence: 22.22%.
    """

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function lions_tigers_and_bears()
    lion = 1
    tiger = 2
    bear = 3

    # Prior
    # probLion  ~ Beta(1,1)
    # probTiger ~ Beta(1,1)
    # probBear  ~ Beta(1,1)

    # The Dirichlet distribution ensures that the sum of probabilities is 1
    # i.e. we don't have to ensure this via some specific constraint.
    v = [3, 2, 1] # We saw three lions, two tigers, and one bear
    # v = [2,2,2]
    x ~ Dirichlet(v)

    # The probabilities to calculate ("aliased" for simplicity)
    probLion  = x[1]
    probTiger = x[2]
    probBear  = x[3]

    n = 7
    o = tzeros(n)
    for i in 1:n
        o[i] ~ Categorical([probLion,probTiger,probBear])
    end

    # It shouldn't matter in what order we see the different animals.
    true ~ Dirac(o[1] == lion)
    true ~ Dirac(o[2] == lion)
    true ~ Dirac(o[3] == lion)
    true ~ Dirac(o[4] == tiger)
    true ~ Dirac(o[5] == tiger)
    true ~ Dirac(o[6] == bear)

    # return o[7] == lion
    # return o[7] == tiger
    return o[7] == bear

end

model = lions_tigers_and_bears()
num_chains = 4

# chains = sample(model, Prior(), 10_000)

# chains = sample(model, MH(), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, MH(), 10_000)

# chains = sample(model, PG(15), MCMCThreads(), 1_000, num_chains)
# chains = sample(model, PG(15), 1_000)

chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, SMC(1000), 10_000)

# chains = sample(model, IS(), 1_000)

#
display(chains)
show_var_dist_pct(chains,:len,1000)

genq = generated_quantities(model, chains)
show_var_dist_pct(genq,1000)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module RandUtils

using Compat

import Compat: UTF8String, view

const SEED = MersenneTwister(1)
const DEFAULT_ELTYPE = typeof(rand())

samerand(args...) = rand(deepcopy(SEED), args...)
samesprand(args...) = sprand(deepcopy(SEED), args...)

if VERSION >= v"0.5.0-dev+3807"
    samesprandbool(args...) = sprand(deepcopy(SEED), Bool, args...)
else
    samesprandbool(args...) = sprandbool(deepcopy(SEED), args...)
end

randvec(T, n) = samerand(T, n)
randvec(n) = samerand(n)

randmat(T, n) = samerand(T, n, n)
randmat(n) = samerand(n, n)

export samerand, samesprand, samesprandbool, randvec, randmat

end


text:
########### state #############
function CSRound!(CS::CarSt, spre::Float64, vpre::Float64)
    CS.s = round(CS.s, spre)
    CS.v = round(CS.v, vpre)
    return CS
end

function CSRound(CS::CarSt, spre::Float64, vpre::Float64)
    newCS = deepcopy(CS)
    CSRound!(newCS, spre, vpre)
    return newCS
end

function SSRound(SS::Sts, spre::Float64, vpre::Float64)
    newSS = deepcopy(SS) #
    CSRound!(newSS.Ego, spre, vpre)
    CSRound!(newSS.Other, spre, vpre)
    return newSS
end

# helper function to handel the situation when calculated s and v go out of range
function sv_boundry!(DP::DrivePOMDP, CS::CarSt)
    CS.s = min(CS.s, DP.sV[end])
    CS.s = max(CS.s, DP.sV[1])
    CS.v = min(CS.v, DP.vV[end])
    CS.v = max(CS.v, DP.vV[1])
    return nothing
end

POMDPs.states(DP::DrivePOMDP) = DP.SSpace # the state space in this model
POMDPs.n_states(DP::DrivePOMDP) = length(DP.SSpace)
function POMDPs.stateindex(DP::DrivePOMDP, Ss::Sts)
      idx = findfirst(x->x==Ss, DP.SSpace)
      if idx == nothing
          Ss_R = SSRound(Ss, DP.Δs, DP.Δv)
          sv_boundry!(DP, Ss_R.Ego)
          sv_boundry!(DP, Ss_R.Other)
          idx_R = findfirst(x->x==Ss_R, DP.SSpace)
          if idx_R == nothing
              error("state not found: $Ss_R")
          else
              return idx_R
          end
      else
          return idx
      end
end


text:
include("razzle-dazzle-root-beer.jl")
include("rosettacode-amb.jl")
include("sudoku.jl")
include("emirp.jl")
include("misc.jl")


text:
# `ql.jl` requires Turing, Distributions, Plots, and StatsPlots
using Turing, Distributions, Plots, StatsPlots

# definitions of type alias
Action = Int64
Reward = Real
Probability = Real
QValue = Real

# definition of the environment
struct Bandit
    p::Array{Probability, 1}
end

function step(bandit::Bandit, action::Action)
    p = bandit.p[action]
    if p >= rand()
        return 1.
    end
    return 0.
end

# definition of the model
mutable struct Agent
    q::Array{QValue, 1}
    alpha::Real
    beta::Real
end

function from_env(alpha::Real, beta::Real, env::Bandit)
    n = length(env.p)
    q = zeros(n)
    return Agent(q, alpha, beta)
end

function softmax(agent::Agent)
    qmax = maximum(agent.q)
    qexp = exp.(agent.beta * (agent.q .- qmax))
    return qexp / sum(qexp)
end

function choose_action(agent::Agent)
    p = softmax(agent)
    return rand(Categorical(p))
end

function update(agent::Agent, action::Action, reward::Reward)
    agent.q[action] += agent.alpha * (reward - agent.q[action])
end

# definition of steps in the simulation
function step(agent::Agent, env::Bandit)
    action = choose_action(agent)
    reward = step(env, action)
    update(agent, action, reward)
    return action, reward
end

function run(agent::Agent, env::Bandit, trial::Int64)
    result = Array{Tuple{Action, Reward}, 1}(undef, trial)
    for i in 1:trial
        ar = step(agent, env)
        result[i] = ar
    end
    return result
end

# run simulation and fit q-learning model to the simulated data
bandit = Bandit([0.1, 0.5])
agent = from_env(0.1, 2., bandit)
result = run(agent, bandit, 1000)

@model QLearning(env::Bandit, result::Array{Tuple{Action, Reward}}) = begin
    T = length(result)
    alpha ~ Beta(1, 1)
    beta ~ Gamma(1., 100.)
    agent = from_env(alpha, beta, env)

    for t in 1:(T - 1)
        a, r = result[t]
        p = softmax(agent)
        result[t][1] ~ Categorical(p)
        update(agent, a, r)
    end
end

chain = sample(QLearning(bandit, result), NUTS(), 500)
plot(chain)


text:
function genLPM(lammbda::AbstractVector, T::Int, X::AbstractMatrix)
    # generating the Linear Probability Model
    # lammbda is a vector (n x 1)
    # X is the matrix (T x n)

    y = BitArray(undef, T)
    for i = 1:T
        y[i] = dot(X[i,:], lammbda) > randn()
    end

    return y

end

genMC(PIt::Array{Float64,3}, T::Int; Tdrop::Int = 0, S0::Int = 1) = genEndoMC(PIt, T; Tdrop = Tdrop, S0 = S0)

function genEndoMC(PIt::Array{Float64,3}, T::Int; Tdrop::Int = 0, S0::Int = 1)
    S  = [S0]; # fix starting point
    for i = 2:T + Tdrop
        cumu_p = cumsum(PIt[S[i - 1], :,i])
        push!(S, sum(rand() .> cumu_p) + 1)
    end
    return S[Tdrop + 1:end]
end

simMC(PI::Array{Float64,2}, T::Int; Tdrop::Int = 0, S0::Int = 1) = genMC(PI, T, Tdrop = Tdrop, S0 = S0)

function genMC(PI::Array{Float64,2}, T::Int; Tdrop::Int = 0, S0::Int = 1)
    S  = [S0] # fix starting point
    cumPI = cumsum(PI;dims = 2)
    for i = 2:T
        push!(S, sum(rand() .> cumPI[S[i - 1], :]) + 1);
    end
    return S[Tdrop + 1:end]
end

max_stars_repo_name: zymbuzz/myJuliaUtils

max_issues_repo_name: zymbuzz/myJuliaUtils

max_forks_repo_name: zymbuzz/myJuliaUtils


text:
#
# Generate initial-velocities
#

import Random

"""

```
initial_velocities(N;T=298.15,R=8.3145,m=1.0,seed=123)
```

Generates a vector of initial velocities for `N` particles of
mass `m`, in two dimensions according to the Maxwell-Boltzmann distribution. By
default the result will be output in ms⁻¹, according to the dimensions
of the gas constant `R`.

## Example

```
v = initial_velocities(10_000)
```

Use a random seed:

```
v = initial_velocities(10_000,100,0.5,seed=-1)
```

"""
function initial_velocities(N;T=298.15,R=8.3145,m=1.0,seed=123)
  # Use a different seed if desired
  seed > 0 ?  Random.seed!(seed) : Random.seed!()
  # Create velocity vector
  v = [ initial_velocity(T,R,m) for i in 1:N ]
  return v
end

# initial_velocity for one particle
function initial_velocity(T,R,m)
  # Inverse of cumulative distribution of p(v)dv
  vmod = sqrt( -(2*R*T/m) * log(1-rand()) )
  θ = 2π*rand()
  return [ vmod*cos(θ), vmod*sin(θ) ]
end

# From the data structure
initial_velocities(data::Data;T=298.15,R=8.3145,m=1.0,seed=123) =
   initial_velocities(data.N,T=T,R=R,m=m,seed=seed)


text:
const confpath = joinpath(dirname(dirname(@__FILE__)), "deps", "config.json")
if isfile(confpath)
    const confd = JSON.parsefile(confpath)
    if !haskey(confd, "dbsock")
        error("set `dbsock` key in $(confpath) to have a valid ",
            "ZeroMQ connection string.")
    elseif !haskey(confd, "dsn")
        error("set `dsn` key in $(confpath) to the name of a valid ODBC DSN.")
    end

    if haskey(confd, "password") && !haskey(confd, "username")
        error("set `username` key in $(confpath) if providing a `password` key.")
    end

    if !haskey(confd, "username")
        confd["username"] = ""
    end
    if !haskey(confd, "password")
        confd["password"] = ""
    end
else
    error("configuration file not found at $(confpath).")
end


text:
using GLVisualize, GeometryTypes, Reactive

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0f0,1f0,360))
end

description = """
Visualizing an animated 3D vector field.
"""

N = 7
# generate some rotations
function rotation_func(t)
    t = (t == 0f0 ? 0.01f0 : t)
    [Vec3f0(sin(x/t), cos(y/(t/2f0)), sqrt(t+z^2))*3f0 for x=1:N, y=1:N, z=1:N]
end

# us Reactive.map to transform the timesignal signal into the arrow flow
flow = map(rotation_func, timesignal)

# create a visualisation
s = step(linspace(0.,1.,N))
vis = visualize(flow, scale = Vec3f0(0.1))
_view(vis, window)

if !isdefined(:runtests)
    @async renderloop(window)
end
# vis.children[].vertexarray.program
# dot(Vec3f0(0, 0, 1), Vec3f0(0, 0, 1))
#
# _view(visualize(
#     (Pyramid(Point3f0(0, 0, -0.5), 1f0, 1f0), [Point3f0(0)]),
#     rotation = [Vec3f0(1, 0, 0)]
# ))


text:
@testset "790.domino-and-tromino-tiling.jl" begin
    let
        res = [1, 2, 5, 11, 24, 53, 117, 258, 569, 1255]
        for i in 1:10
            @test num_tilings(i) == res[i]
        end
    end
end


text:
const FluxOptimizer = Union{
    Descent, Momentum, Nesterov, RMSProp,
    ADAM, AdaMax, ADAGrad, ADADelta, RADAM,
    AMSGrad, NADAM
}

abstract type OptimParam end
struct LR <: OptimParam end
struct Mom <: OptimParam end

optimparamfield(::O, op::Type{P}) where {O, P<:OptimParam} = error(
    "No parameter $(P) on optimizer $(O)."
)
optimparamfield(::FluxOptimizer, ::Type{LR}) = :eta

function setoptimparam!(opt, optimparam::Type{<:OptimParam}, x)
    setfield!(opt, optimparamfield(opt, optimparam), x)
end
getoptimparam(opt, optimparam::Type{<:OptimParam}) = getfield(opt, optimparamfield(opt, optimparam))

# TODO: add momentum, weight decay and other...
getoptimparams(opt) = Dict(LR => getoptimparam(opt, LR))


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Qt5Tools_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Qt5Tools")
JLLWrappers.@generate_main_file("Qt5Tools", UUID("a9c6e4b1-b2fb-56d5-96a9-25f276f13840"))
end  # module Qt5Tools_jll


text:
module GenericBlas

import BlasIntrinsics

using Base.Threads

#https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects-1

#create function/type alias
#SSE
double2   = BlasIntrinsics.double2
float4    = BlasIntrinsics.float4
#AVX/AVX2
double4   = BlasIntrinsics.double4
float8    = BlasIntrinsics.float8
#AVX512
double8   = BlasIntrinsics.double8
float16   = BlasIntrinsics.float16

mm_mul    = BlasIntrinsics.mm_mul
mm_add    = BlasIntrinsics.mm_add
mm_sub    = BlasIntrinsics.mm_sub
mm_mad    = BlasIntrinsics.mm_mad

#collection of generic blas operations (e.g., on vector registers etc.)

abstract type GenericMultiblasFunctor{T1<:Any, T2<:Any} end

mutable struct ConvertC2R{DstTp <: Any, SrcScalarTp <: AbstractFloat} <: GenericMultiblasFunctor{DstTp, SrcScalarTp}
  #scalars:
  #vectors
  dst::Vector{DstTp}
  src::Vector{Complex{SrcScalarTp}}

  ConvertC2R{DstTp, SrcScalarTp}(y::Vector{DstTp}, x::Vector{Complex{SrcScalarTp}}) where {DstTp <: Any} where {SrcScalarTp <: AbstractFloat} = new(y, x)

end # struct ConvertC2R

# Helper methods
#function register_t1(f::GenericMultiblasFunctor{T1<:Any, T2<:Any}) where T1 where T2; return T1; end
#function register_t2(f::GenericMultiblasFunctor{T1<:Any, T2<:Any}) where T1 where T2; return T2; end

@inline function (f::ConvertC2R)(i::Int64)
             dst_t = Float32 #register_t1(f)
             if     dst_t == double4 || dst_t == float4
               c=[real(f.src[2i-1]), imag(f.src[2i-1]), real(f.src[2i]), imag(f.src[2i]) ]
               f.dst[i] = dst_t(ntuple(i->c[i],4))
             elseif dst_t == double8 || dst_t == float8
               c=[real(f.src[4i-3]), imag(f.src[4i-3]), real(f.src[4i-2]), imag(f.src[4i-2]), real(f.src[4i-1]), imag(f.src[4i-1]), real(f.src[4i]), imag(f.src[4i]) ]
               f.dst[i] = dst_t(ntuple(i->c[i],8))
             elseif dst_t ==float16
               c=[real(f.src[8i-7]), imag(f.src[8i-7]),
                  real(f.src[8i-6]), imag(f.src[8i-6]),
                  real(f.src[8i-5]), imag(f.src[8i-5]),
                  real(f.src[8i-4]), imag(f.src[8i-4]),
                  real(f.src[8i-3]), imag(f.src[8i-3]),
                  real(f.src[8i-2]), imag(f.src[8i-2]),
                  real(f.src[8i-1]), imag(f.src[8i-1]),
                  real(f.src[8i  ]), imag(f.src[8i  ]) ]
               f.dst[i] = dst_t(ntuple(i->c[i],16))
             elseif dst_t <: AbstractFloat
               f.dst[2i-1] = real(src[i])
               f.dst[2i  ] = real(src[i])
             else
               println("Register type is not supported.")
             end

end

#mutable struct ConvertR2C{DstScalarTp <: AbstractFloat, SrcTp <: Any} <: GenericMultiblasFunctor{DstScalarTp, SrcTp}
  #scalars:
  #vectors
#  dst::Vector{Complex{DstScalarTp}}
#  src::Vector{SrcTp}

#  ConvertC2R{DstScalarTp, SrcTp}(y_::Vector{Complex{DstScalarTp}}, x_::Vector{SrcTp}) where {DstScalarTp <: AbstractFloat} where {SrcTp <: Any} = new(y_, x_)

#end #

#@inline function (f::ConvertR2C)(i::Int64)
#          f.dst[i] = f.src[2i-1] + f.src[2i]*im;
#end

mutable struct cpy{DstTp <: AbstractArray, SrcTp <: AbstractArray} <: GenericMultiblasFunctor{DstTp, SrcTp}
  #scalars:
  #vectors
  dst::Vector{DstTp}
  src::Vector{SrcTp}

  cpy{DstTp, SrcTp}(y_::DstTp, x_::SrcTp) where {DstTp <: AbstractArray} where {SrcTp <: AbstractArray}  = new(y_, x_)

end #

@inline function (f::cpy)(i::Int64)
          f.dst[i] = f.src[i]
end

#@inline function cpy(y::AbstractArray, x::AbstractArray)
#  if pointer_from_objref(y) == pointer_from_objref(x); return; end
#@threads for i in 1:length(x);  y[i] = x[i]; end
#end #cpy

mutable struct axpy{ScalarTp <: AbstractFloat, VectorTp <: Any} <: GenericMultiblasFunctor{ScalarTp, VectorTp}
  #scalars:
  a::ScalarTp
  #vectors
  x::Vector{VectorTp}
  y::Vector{VectorTp}

  axpy{ScalarTp, VectorTp}(a_::ScalarTp, x_::Vector{VectorTp}, y_::Vector{VectorTp}) where {ScalarTp <: AbstractFloat} where {VectorTp <: Any}  = new(a_, x_, y_)

end #

@inline function (f::axpy)(i::Int64)
          f.y[i] = mm_mad(f.a, f.x[i], f.y[i])
end

mutable struct xpay{VectorTp <: Any, ScalarTp <: AbstractFloat} <: GenericMultiblasFunctor{VectorTp, ScalarTp}

  #vectors
  x::Vector{VectorTp}
  y::Vector{VectorTp}
  #scalars:
  a::ScalarTp

  xpay{VectorTp, ScalarTp}(x_::Vector{VectorTp}, y_::Vector{VectorTp}, a_::ScalarTp) where {VectorTp <: Any} where {ScalarTp <: AbstractFloat}  = new(x_, y_, a_)

end #

@inline function (f::xpay)(i::Int64)
          f.y[i] = mm_mad(f.a, f.y[i], f.x[i])
end

mutable struct xmy{VectorTp <: Any} <: GenericMultiblasFunctor{VectorTp, VectorTp}

  #vectors
  x::Vector{VectorTp}
  y::Vector{VectorTp}

  xmy{VectorTp}(x_::Vector{VectorTp}, y_::Vector{VectorTp}) where {VectorTp <: Any}  = new(x_, y_)

end #

@inline function (f::xmy)(i::Int64)
          f.y[i] = mm_sub(f.x[i], f.y[i])
end

mutable struct xpy{VectorTp <: Any} <: GenericMultiblasFunctor{VectorTp, VectorTp}

  #vectors
  x::Vector{VectorTp}
  y::Vector{VectorTp}

  xpy{VectorTp}(x_::Vector{VectorTp}, y_::Vector{VectorTp}) where {VectorTp <: Any}  = new(x_, y_)

end #

@inline function (f::xpy)(i::Int64)
          f.y[i] = mm_add(f.x[i], f.y[i])
end

#First performs the operation y[i] += a*x[i]
#Second performs the operator x[i] -= a*z[i]

mutable struct axpyXmaz{ScalarTp <: AbstractFloat, VectorTp <: Any} <: GenericMultiblasFunctor{ScalarTp, VectorTp}
  #scalar
  a::ScalarTp
  #vectors
  x::Vector{VectorTp}
  y::Vector{VectorTp}
  z::Vector{VectorTp}

  axpyXmaz{ScalarTp, VectorTp}(a_::ScalarTp, x_::Vector{VectorTp}, y_::Vector{VectorTp}, z_::Vector{VectorTp}) where {ScalarTp <: AbstractFloat} where {VectorTp <: Any}  = new(a_, x_, y_, z_)

end #

@inline function (f::axpyXmaz)(i::Int64)
          f.y[i] = mm_mad(+f.a, f.x[i], f.y[i])
          f.x[i] = mm_mad(-f.a, f.z[i], f.x[i])
end

#First performs the operation x[i] = x[i] + a*p[i]
#Second performs the operator p[i] = u[i] + b*p[i]

mutable struct axpyZpbx{ScalarTp <: AbstractFloat, VectorTp <: Any} <: GenericMultiblasFunctor{ScalarTp, VectorTp}
  #scalar
  a::ScalarTp
  b::ScalarTp
  #vectors
  x::Vector{VectorTp}
  p::Vector{VectorTp}
  u::Vector{VectorTp}

  axpyZpbx{ScalarTp, VectorTp}(a_::ScalarTp, b_::ScalarTp, x_::Vector{VectorTp}, p_::Vector{VectorTp}, u_::Vector{VectorTp}) where {ScalarTp <: AbstractFloat} where {VectorTp <: Any}  = new(a_, b_, x_, p_, u_)

end #

@inline function (f::axpyZpbx)(i::Int64)
          f.x[i] = mm_mad(f.a, f.p[i], f.x[i])
          f.p[i] = mm_mad(f.b, f.p[i], f.u[i])
end

mutable struct norm2{VectorTp <: Any} <: GenericMultiblasFunctor{VectorTp, VectorTp}

  #vectors
  x::Vector{VectorTp}

  norm2{VectorTp}(x_::Vector{VectorTp}) where {VectorTp <: Any}  = new(x_)

end #

@inline function (f::norm2)(i::Int64)
          return mm_mul(f.x[i], f.x[i])
end

@inline function plus(x::AbstractFloat, y::AbstractFloat)
          return mm_add(x, y)
end

end #GenericBlas

max_stars_repo_name: alexstrel/solver-tuner

max_issues_repo_name: alexstrel/solver-tuner

max_forks_repo_name: alexstrel/solver-tuner


text:
###############################################################################
#
# Leaf parameters container
#
###############################################################################
"""
    mutable struct Leaf{FT}

Struct to store leaf information.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef mutable struct Leaf{FT<:AbstractFloat}
    # Temperature related
    "Temperature `[K]`"
    T    ::FT = 298.15
    "Old Temperature `[K]`, if not T, run leaf_temperature_dependence!"
    T_old::FT = 0

    # Photosynthesis system
    "Rate constant for thermal dissipation"
    Kd       ::FT = 0.85
    "Rate constant for fluorescence (const)"
    Kf       ::FT = 0.05
    "Reversible NPQ rate constant (initially zero)"
    Kr       ::FT = 0
    "Sustained NPQ rate constant (for seasonal changes, default is zero)"
    Ks       ::FT = 0
    "Rate constant for photochemistry (all reaction centers open)"
    Kp       ::FT = 4
    "Maximal Kp"
    Kp_max   ::FT = 4
    "max PSII yield (Kr=0, all RC open)"
    maxPSII  ::FT = Kp/(Kp+Kf+Kd)
    "Fraction of absorbed light used by PSII ETR"
    PSII_frac::FT = 0.5

    # CO₂ pressures
    "Leaf internal CO₂ partial pressure `[Pa]`"
    p_i  ::FT = 10
    "Leaf surface CO₂ partial pressure `[Pa]`"
    p_s  ::FT = 40
    "Saturation H₂O vapor pressure `[Pa]`"
    p_sat::FT = saturation_vapor_pressure(T)
    "Leaf diffusive conductance to CO₂ `[mol m⁻² s⁻¹]`"
    g_bc ::FT = 3 / 1.35
    "Leaf diffusive conductance to CO₂ `[mol m⁻² s⁻¹]`"
    g_lc ::FT = 0.01

    # Photosynthesis related
    "RubisCO limited photosynthetic rate `[μmol m⁻² s⁻¹]`"
    Ac     ::FT = 0
    "Light limited photosynthetic rate `[μmol m⁻² s⁻¹]`"
    Aj     ::FT = 0
    "Gross photosynthetic rate `[μmol m⁻² s⁻¹]`"
    Ag     ::FT = 0
    "Net photosynthetic rate `[μmol m⁻² s⁻¹]`"
    An     ::FT = 0
    "Product limited photosynthetic rate `[μmol m⁻² s⁻¹]`"
    Ap     ::FT = 0
    "Electron transport `[μmol m⁻² s⁻¹]`"
    J      ::FT = 0
    "Potential Electron Transport Rate `[μmol m⁻² s⁻¹]`"
    J_pot  ::FT = 0
    "Maximal electron transport rate `[μmol m⁻² s⁻¹]`"
    Jmax   ::FT = 120
    "Maximal electron transport rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Jmax25 ::FT = 120
    "RubisCO coefficient Kc `[Pa]`"
    Kc     ::FT = 0
    "RubisCO coefficient Ko `[Pa]`"
    Ko     ::FT = 0
    "PEP coefficient Ko `[Pa]`"
    Kpep   ::FT = 0
    "Michaelis-Menten's coefficient `[Pa]`"
    Km     ::FT = 0
    "Respiration rate `[μmol m⁻² s⁻¹]`"
    Rd     ::FT = 1
    "Respiration rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Rd25   ::FT = 1
    "Maximal carboxylation rate `[μmol m⁻² s⁻¹]`"
    Vcmax  ::FT = 60
    "Maximal carboxylation rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Vcmax25::FT = 60
    "Maximal PEP carboxylation rate `[μmol m⁻² s⁻¹]`"
    Vpmax  ::FT = 120
    "Maximal PEP carboxylation rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Vpmax25::FT = 120
    "CO₂ compensation point with the absence of Rd `[Pa]`"
    Γ_star ::FT = 0

    # Well watered condition values (to use with β function over PS)
    "Well watered maximal electron transport rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Jmax25WW ::FT = Jmax25
    "Well watered respiration rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Rd25WW   ::FT = Rd25
    "Well watered maximal carboxylation rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Vcmax25WW::FT = Vcmax25
    "Well watered maximal PEP carboxylation rate at 298.15 K `[μmol m⁻² s⁻¹]`"
    Vpmax25WW::FT = Vpmax25

    # Fluorescence related
    "Total efficiency, incl. photorespiration `[mol CO₂ mol⁻¹ e-]`"
    e2c::FT = 1 / 6
    "dark adapted yield (`Kp=0`)"
    Fm ::FT = 0
    "light adapted yield (`Kp=0`)"
    Fm′::FT = 0
    "dark-adapted fluorescence yield (`Kp=max`)"
    Fo ::FT = 0
    "light-adapted fluorescence yield in the dark (`Kp=max`)"
    Fo′::FT = 0
    "Actual electron transport rate `[μmol m⁻² s⁻¹]`"
    Ja ::FT = 0
    "Non-Photochemical quenching "
    NPQ::FT = 0
    "Photochemical quenching"
    qQ ::FT = 0
    "energy quenching"
    qE ::FT = 0
    "PSII yield"
    φ  ::FT = 0
    "Steady-state (light-adapted) yield (aka Fs)"
    ϕs ::FT = 0

    # Environment related
    "Absorbed photosynthetic active radiation `[μmol m⁻² s⁻¹]`"
    APAR::FT = 100
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using CSV
using DataFrames
using DelimitedFiles
using Glob
using Plots; pgfplotsx()

# Create plots
function plot_gmres_orth()

	# Plot labels
	labels = Dict(
		# kernels
		"ser" => "SparseArrays", "par" =>  "MtSpMV.jl", "mkl" => "",
		# formats
		"csr" => "", "csc" => "",
		# orth
		"mgs" => ", MGS", "cgs" => ", CGS", "dgks" => ", DGKS",
		# precon
		"none" => "", "jacobi" => "", "ilu0" => "",
		"parjacobi" => "", "parilu0" => "",
	)

   	# Plot linestyles
   	linestyles = [:solid, :dash, :dot]

	# Benchmark folders
	folders = ["1e6_100-samples", "1e6_1000-samples", "1e7_100-samples"]

	# Benchmark subfolders
	subfolders = ["par-ser", "ser-par", "par-par"]

	# Iterate over matrix size number of samples folders
	for (ifolder, folder) in enumerate(folders)

		# Get folders with csv files
		csvfolder_parser = joinpath(folder, subfolders[1])
		csvfolder_parser = joinpath(csvfolder_parser, readdir(csvfolder_parser)[1])
		csvfolder_serpar = joinpath(folder, subfolders[2])
		csvfolder_serpar = joinpath(csvfolder_serpar, readdir(csvfolder_serpar)[1])
		csvfolder_parpar = joinpath(folder, subfolders[3])
		csvfolder_parpar = joinpath(csvfolder_parpar, readdir(csvfolder_parpar)[1])
		csvfolders = [csvfolder_parser, csvfolder_serpar, csvfolder_parpar]

	   	# Get result for serial sparse serial dense (orth=mgs)
	   	file = joinpath(csvfolder_serpar, "csv", "runtimes_ser_csr_mgs_none.csv")
	   	df = DataFrame!(CSV.File(file, header=true))
	   	y = convert(Vector, df[!, 5])
	   	y_ser_ser = similar(y)
	   	y_ser_ser .= y[1]

	   	# Iterate over csv folders
	   	for (icsvfolder, csvfolder) in enumerate(csvfolders)

	   		# Create plots
	   		plt1 = plot(ylabel="Runtime in seconds", xlabel="Number of threads",
	   					legend=:topright, framestyle=:box)
	   		plt2 = plot(ylabel="Parallel speedup", xlabel="Number of threads",
	   					legend=:none, framestyle=:box)

			# Get all csv files with runtimes in sorted order
			files = glob(string("runtimes*"), joinpath(csvfolder, "csv"))
			if length(files) == 3
				files = [files[3], files[1], files[2]]
			end

			# Add results from all files to plots
			for (ifile, file) in enumerate(files)

				# Get benchmark results
			    df = DataFrame!(CSV.File(file, header=true))
			    x = convert(Vector, df[!, 1])
			    y = convert(Vector, df[!, 5])

			    # Create label from filename
			    infos = file[findlast('/', file) + 1 : findlast('.', file) - 1]
			    infos = infos[findfirst('_', infos) + 1 : length(infos)]
			    infos = split(infos, "_")
			    label = ""
			    for info in infos
					label = string(label, labels[info])
				end

				# Add runtimes to first plot
				plot!(plt1, x, y, label=label, linestyles=linestyles[ifile])

				# Add speedups to second plot
				y_speedups = y_ser_ser ./ y
				plot!(plt2, x, y_speedups, linestyles=linestyles[ifile])

			end # for

			# Plot side by side and save
			width = 530
			height = width / 2
			plot(plt1, plt2, layout=(1,2), size=(width,height))
			filename = string("gmres_orth_", folder, "_", subfolders[icsvfolder])
			savefig(string(filename, ".png"))
			savefig(string(filename, ".tex"))

		end # for

	end # for

end # function

max_stars_repo_name: jkrch/JuliaParGMRES-Plots

max_issues_repo_name: jkrch/JuliaParGMRES-Plots

max_forks_repo_name: jkrch/JuliaParGMRES-Plots


text:
#!/usr/bin/env julia
using RobotOS
@rosimport geometry_msgs.msg: Point, Pose, Pose2D, PoseStamped, Vector3, Twist
@rosimport std_srvs.srv: Empty, SetBool
@rosimport nav_msgs.srv.GetPlan
@rosimport gazebo_msgs.msg: ModelState, ModelStates
@rosimport gazebo_msgs.srv: SetModelState, GetModelState, GetWorldProperties, GetLinkState, SetLinkState

rostypegen()
using geometry_msgs.msg
using std_srvs.srv
using nav_msgs.srv
using nav_msgs.srv.GetPlan
using gazebo_msgs.srv
using gazebo_msgs.msg

function loop(get_link_state, pub)
    loop_rate = Rate(100.0)
    while ! is_shutdown()
        linkName = "base_footprint"  # TODO make this a parameter

        # Get the current position of the Gazebo model
        gs = GetLinkStateRequest()
        gs.link_name = linkName
        gs_r = get_link_state(gs)
        if !gs_r.success
            error(string(" calling /gazebo/get_link_state service: ", gs_r.status_message))
        end

        npt = gs_r.link_state.pose
        publish(pub, npt)
        rossleep(loop_rate)
    end
end

function main()
    init_node("rosjl_link_position")

    # Set up service to get Gazebo link state
    const get_link_state = ServiceProxy("/gazebo/get_link_state",GetLinkState)
    println("Waiting for 'gazebo/get_link_state' service...")
    wait_for_service("gazebo/get_link_state")

    # TODO use getparam to get robot name
    pub = Publisher{Pose}("/hmmwv/base_footprint_link_pose", queue_size=10)

    loop(get_link_state, pub)
end

if ! isinteractive()
    main()
end


text:
module ComplexBC

import JuMP
import MathOptInterface
import BranchAndBound
import PowerModels
import Mosek, MosekTools
import Ipopt
import LinearAlgebra: norm, rank

const MOI = MathOptInterface
const BB = BranchAndBound
const PM = PowerModels

const MOSEK_OPTIMIZER = JuMP.optimizer_with_attributes(Mosek.Optimizer, "QUIET" => true)
const MAX_NODE_ALLOWED = 10000

include("utils.jl")
include("data_structures.jl")
include("pm_exts.jl")
include("key_steps.jl")
include("bb_exts.jl")
include("initialize.jl")

function main(file::String)::Nothing
    data = PM.parse_file(file)
    pm = PM.instantiate_model(data, NodeWRMPowerModel, PM.build_opf)
    optimizer = JuMP.optimizer_with_attributes(Mosek.Optimizer, "QUIET" => true)
    # optimizer = optimizer_with_attributes(SCS.Optimizer, "max_iters" => 10000, "verbose" => 0)
    JuMP.set_optimizer(pm.model, optimizer)

    tree, node = initialize(pm, file)

    @time BB.run(tree)

    if isempty(node.auxiliary_data["rank1"])
        println("No rank-1 solution found")
    else
        (bound, id) = findmin(node.auxiliary_data["rank1"])
        println("Best bound obtained at $(id), bound value $(bound)")
    end
end

end # module

max_stars_repo_name: kibaekkim/ComplexBC

max_issues_repo_name: kibaekkim/ComplexBC

max_forks_repo_name: kibaekkim/ComplexBC


text:
"""
A stochastic flexibility design problem featuring a rounding strategy. Inspired
on an example problem from Pulsipher JL, Zavala VM. A scalable stochastic
programming approach for the design of flexible systems. Computers & Chemical
Engineering. 2019 Sep 2;128:69-76.
"""

using InfiniteOpt, JuMP, Clp, Distributions

# Set the covariance matrix for the uncertain parameters
θ_nom = [0.; 60.; 10.]
covar = [80. 0 0; 0 80. 0; 0 0 120.]

# Set the dimensions
n_z = 3; n_θ = 3; n_d = 3

# Set the problem parameters
c = ones(n_d) / sqrt(n_d)
c_max = 5
U = 10000
num_samples = 100

# Initialize the model
m = InfiniteModel(Clp.Optimizer, seed = true) # use seed to get same samples
set_silent(m)

# Set the uncertainty parameters
@infinite_parameter(m, θ[i = 1:n_θ] in MvNormal(θ_nom, covar))

# Initialize the variables
@infinite_variable(m, 0 <= y(θ) <= 1)
# @infinite_variable(m, y(θ), Bin)
@infinite_variable(m, z[1:n_z](θ))
@hold_variable(m, d[1:n_d] >= 0)

# Set objective function
@objective(m, Max, expect(1 - y, θ, num_supports = num_samples))

# Set the line capacity constraints
@constraint(m, f1, -z[1] - 35 - d[1] <= y * U)
@constraint(m, f2, z[1] - 35 - d[1] <= y * U)
@constraint(m, f3, -z[2] - 50 - d[2] <= y * U)
@constraint(m, f4, z[1] - 50 - d[2] <= y * U)

# Set the generator capacity constraints
@constraint(m, f5, -z[3] <= y * U)
@constraint(m, f6, z[3] - 100 - d[3] <= y * U)

# Set the node balance constraints
@constraint(m, h1, z[1] - θ[1] == 0)
@constraint(m, h2, -z[1] -z[2] + z[3] - θ[2] == 0)
@constraint(m, h3, z[2] - θ[3] == 0)

# Enforce the minimum SF
@constraint(m, max_cost, sum(c[i] * d[i] for i = 1:n_d) <= c_max)

# Solve and and obtain results
optimize!(m)
if has_values(m)
    opt_y = value(y)
    opt_d = value.(d)
    opt_obj = objective_value(m)
end

# Estimate the value of SF
SF = 1 - sum(opt_y .>= 1e-8) / num_samples

# Print the results
print("------------------RESULTS------------------\n")
print("Optimal Objective:     ", opt_obj, "\n")
print("Optimal Cost:          ", sum(c[i] * opt_d[i] for i = 1:n_d), "\n")
print("Maximum Cost:          ", c_max, "\n")
print("Predicted SF:          ", 100 * SF, "%\n")
print("Optimal Design Values: ", opt_d, "\n\n")


text:
# AUTO GENERATED FILE - DO NOT EDIT

export column

"""
    column(;kwargs...)
    column(children::Any;kwargs...)
    column(children_maker::Function;kwargs...)

A Column component.
Row Column
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional): The children of the element
- `style` (Dict; optional): The inline styles
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `className` (String; optional): The class of the element
- `columnSizes` (Array of Strings; optional): The size of the column with the display size, sm-4, lg-16 ...
- `offsetSizes` (Array of Strings; optional): The size of the offset with the display size, lg-2 ...
"""
function column(; kwargs...)
        available_props = Symbol[:children, :style, :id, :className, :columnSizes, :offsetSizes]
        wild_props = Symbol[]
        return Component("column", "Column", "dash_carbon_components", available_props, wild_props; kwargs...)
end

column(children::Any; kwargs...) = column(;kwargs..., children = children)
column(children_maker::Function; kwargs...) = column(children_maker(); kwargs...)

max_stars_repo_name: Matheus-Rangel/dash-carbon-components

max_issues_repo_name: Matheus-Rangel/dash-carbon-components

max_forks_repo_name: Matheus-Rangel/dash-carbon-components


text:
using StaticArrays
using Distributions
using Random
using Bridge
using Statistics, LinearAlgebra
using GaussianDistributions
using BridgeSDEInference
const SDE = BridgeSDEInference

## Helper functions and directories

SRC_DIR = joinpath(Base.source_dir(), "..", "..", "src")
OUT_DIR = joinpath(Base.source_dir(), "..", "..", "output")
mkpath(OUT_DIR)

DIR = "auxiliary"
include(joinpath(SRC_DIR, DIR, "data_simulation_fns.jl"))
include(joinpath(SRC_DIR, DIR, "utility_functions.jl"))

### Import example SDE LotkaVolterraDiffusion

using BridgeSDEInference: LotkaVolterraDiffusion, LotkaVolterraDiffusionAux

#=
Stochastic Lotka-Volterra System (SDE) with
parameters (α, β, γ, δ, σ1, σ2) with drift

     b(x) = [α*x[1] - β*x[1]*x[2],
             δ*x[1]*x[2] - γ*x[2]]

and volatility

    σ(x) = [ σ1 0
             0  σ2 ]

# parameters: see src/examples/lotka_volterra.jl
=#

### Generate data with parameters

Random.seed!(4)
θˣ = [1.5, 1.0, 3.0, 1.0, 0.2, 0.2] #  (α, β, γ, δ, σ1, σ2)
Pˣ = LotkaVolterraDiffusion(θˣ...)

x0 = ℝ{2}(2.5, 2.) # starting point
w0 = ℝ{2}(0.0, 0.0) # starting point of driving Wiener process

dt, T =  1/5000, 20.0 # time grid
tt = 0.0:dt:T

XX, _ = simulate_segment(w0, x0, Pˣ, tt) # simulate process

### Subsampling observations

skip = 2000

# Set observation scheme: partial observations of first coordinate
# with observation operator L and Gaussian noise
Σdiagel = 1.0
Σ = @SMatrix[Σdiagel]
L = @SMatrix[1.0 0.0]

# Observe with noise at sample times
obs_time, obs_vals = XX.tt[1:skip:end], [rand(Gaussian(L*x, Σ)) for x in XX.yy[1:skip:end]]

### Setup MCMC scheme to infer parameters given observations

θ_init = copy(θˣ)
Pˣ = LotkaVolterraDiffusion(θ_init...)

# Setup: For observation intervals define linearizations of the diffusion to be used
# in generating MCMC proposals
P̃ = [LotkaVolterraDiffusionAux(θ_init..., t₀, u, T, v) for (t₀, T, u, v)
     in zip(obs_time[1:end-1], obs_time[2:end], obs_vals[1:end-1], obs_vals[2:end])]

# Setup: Partial observations: register observation scheme of each observation
model_setup = DiffusionSetup(Pˣ, P̃, PartObs())
set_observations!(model_setup, [L for _ in P̃], [Σ for _ in P̃], obs_vals, obs_time)
set_imputation_grid!(model_setup, 1/1000)

# Setup: Priors
set_x0_prior!(model_setup,
              GsnStartingPt(x0, @SMatrix [20.0 0.0;
                                          0.0 20.0;]),
              x0)

# Initialize sampler
set_auxiliary!(model_setup; skip_for_save=10^0,
               adaptive_prop=NoAdaptation())
initialise!(eltype(x0), model_setup, Vern7(), false, NoChangePt(100))

# Parameters for adaptively finding good Metropolis stepsize
readj = (100, 0.001, 0.001, 999.9, 0.234, 50)

# Finalize setup:
# Define MCMC steps for 1.) path imputation, 2.) conjugate updates for
# linear parameters α, β, γ, δ and 3.) - 8.) random walk steps
# for parameters α, β, γ, δ, σ1, σ2 (in that order)
mcmc_setup = MCMCSetup(
      Imputation(NoBlocking(), 0.99, Vern7()),
      ParamUpdate(ConjugateUpdt(), [1,2,3,4], θ_init, nothing,
                  MvNormal(fill(0.0, 4), diagm(0=>fill(1000.0, 4))),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, [1,2,3,4]))),
      ParamUpdate(MetropolisHastingsUpdt(), 1, θ_init,
                  UniformRandomWalk(0.1, true), ImproperPrior(),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 1)), readj),
      ParamUpdate(MetropolisHastingsUpdt(), 2, θ_init,
                  UniformRandomWalk(0.1, true), ImproperPrior(),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 2)), readj),
      ParamUpdate(MetropolisHastingsUpdt(), 3, θ_init,
                  UniformRandomWalk(0.1, true), ImproperPrior(),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 3)), readj),
      ParamUpdate(MetropolisHastingsUpdt(), 4, θ_init,
                  UniformRandomWalk(0.1, true), ImproperPrior(),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 4)), readj),
      ParamUpdate(MetropolisHastingsUpdt(), 5, θ_init,
                  UniformRandomWalk(0.1, true), ImproperPosPrior(),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 5)), readj),
      ParamUpdate(MetropolisHastingsUpdt(), 6, θ_init,
                  UniformRandomWalk(0.1, true), ImproperPosPrior(),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P̃, 6)), readj))

# Set up schedule for MCMC: Impute, and estimate parameters α, β, and δ (Steps 1, 3, 4, 6)
schedule = MCMCSchedule(1*10^3, [[1],[3,4,6]], # which (groups of) steps in which order
                        (save=1*10^3, verbose=10^2, warm_up=100,
                         readjust=(x->x%100==0), fuse=(x->false)))

### Run MCMC

Random.seed!(4)
state, history = mcmc(mcmc_setup, schedule, model_setup)

# most importantly: parameter samples
θs = history.θ_chain

# Compare posterior mean with true parameters, for example
@show [mean(θs) θˣ]

### Post processing: plot results

include(joinpath(SRC_DIR, DIR, "plotting_fns.jl"))
plot_chains(history; truth=θˣ)

plot_paths(state, history, schedule; obs=(times=obs_time[2:end],
                     vals=[[v[1] for v in obs_vals[2:end]]], indices=[1]))


text:
using SourceWalk, Test

@testset "SourceWalk" begin

include("patch.jl")

end


text:
# implement Tables.jl's row access interface
# https://tables.juliadata.org/stable/#Implementing-the-Interface-(i.e.-becoming-a-Tables.jl-source)-1

import Tables
# import Tables: istable, rowaccess, columnaccess, rows, columnnames, getcolumn, AbstractRow
import Base: eltype, length, iterate

# a table structure local to TableScraper
struct Table
    rows
    columnnames
end

struct TableRow <: Tables.AbstractRow
    row::Int
    source::Table
end

Tables.istable(::Table)=true
Tables.rowaccess(::Table)=true
Tables.columnaccess(::Table)=true
Tables.columnnames(t::Table)=t.columnnames
Tables.rows(t::Table)=t

Base.eltype(::Table) = TableRow
Base.length(t::Table) = length(t.rows)
Base.iterate(t::Table, st = 1) = st > length(t) ? nothing : (TableRow(st, t), st+1)

function Tables.getcolumn(t::TableRow, ::Type, col::Int, nm::Symbol)
     tbl = getfield(t, :source)
    row = tbl.rows[t.row]
    row[col]
end

function Tables.getcolumn(t::TableRow, i::Int)
    tbl = getfield(t, :source)
    row_num = getfield(t, :row)
    row = tbl.rows[row_num]
    row[i]
end

function Tables.getcolumn(t::TableRow, nm::Symbol)
    tbl = getfield(t, :source)
    row_num = getfield(t, :row)
    row = tbl.rows[row_num]
    col = indexin([string(nm)], tbl.columnnames)[1]
    row[col]
end

Tables.columnnames(t::TableRow) = getfield(t, :source).columnnames


text:
function dataTransf(x::AbstractVecOrMat, transfId::Vector{Int})
    y = similar(x)
    n_loss = 0
    for i = 1:size(x, 2)
        (y[:,i], n_loss_temp) = dataTransf(x[:,i], transfId[i])
        n_loss = max(n_loss, n_loss_temp)
    end
    return (y = y, n_loss = n_loss)
end

function dataTransf(x::AbstractVecOrMat, transfId::Int)
#   Categories:
#   ------------------------------------------------------------------------
#   1. Levels (Y_{t})
#   2. First Difference (Y_{t}-Y_{t-1})
#   3. Seasonal difference in quarterly data ((Y_{t})-(Y_{t-4}))
#   4. Seasonal difference in monthly data ((Y_{t})-(Y_{t-12}))
#   5. Log-levels (Ln(Y_{t}))
#   6. Log-First Difference (Ln(Y_{t})-Ln(Y_{t-1}))
#   7. Seasonal log difference in quarterly data (Ln(Y_{t})-Ln(Y_{t-4}))
#   8. Seasonal log difference in monthly data (Ln(Y_{t})-Ln(Y_{t-12}))
#   9. Detrending Ln(Y_{t}) by HP filter using quarterly data
#   10. Detrending Ln(Y_{t}) by HP filter using monthly data
#   11. Seasonal difference of a detrended Ln(Y_{t}) by HP filter using quarterly data
#   12. Seasonal difference of a detrended Ln(Y_{t}) by HP filter using monthly data
#   13. Seasonal difference of a detrended Ln(Y_{t}) by removing a linear trend
#   14. Detrended Y_{t} by removing a linear trend (for already SA data)
#   15. Detrended Ln(Y_{t}) by removing a linear trend (for already SA data)
#  =========================================================================
    y = similar(x)

    if transfId == 1
        y = copy(x)
        n_loss = 0
    elseif transfId == 2
        y[2:end,:] = x[2:end,:] - x[1:end - 1,:]
        n_loss = 1
    elseif transfId == 3
        y[5:end,:] = x[5:end,:] - x[1:end - 4,:]
        n_loss = 4
    elseif transfId == 4
        y[13:end,:] = x[13:end,:] - x[1:end - 12,:]
        n_loss = 12
    elseif transfId == 5
        y = log.(x)
        n_loss = 0
    elseif transfId == 6
        y[2:end,:] = log.(x[2:end,:]) - log.(x[1:end - 1,:])
        n_loss = 1
    elseif transfId == 7
        y[5:end,:] = log.(x[5:end,:]) - log.(x[1:end - 4,:])
        n_loss = 4
    elseif transfId == 8
        y[13:end,:] = log.(x[13:end,:]) - log.(x[1:end - 12,:])
        n_loss = 12
        # case 7
        #     yt=hpfilter(log(x),1600);
        #     y=log(x)-yt;
        #     n_loss=0;
        # case 8
        #     yt=hpfilter(log(x),14400);
        #     y=log(x)-yt;
        #     n_loss=0;
        # case 9
        #     yt=hpfilter(log(x),1600);
        #     yt=log(x)-yt;
        #     y(5:end,:)=yt(5:end,:)-yt(1:end-4,:);
        #     n_loss=4;
        # case 10
        #     yt=hpfilter(log(x),14400);
        #     yt=log(x)-yt;
        #     y(13:end,:)=yt(13:end,:)-yt(1:end-12,:);
        #     n_loss=12;
    elseif transfId == 13
        yt = log.(x)
        yt = yt[13:end,:] - yt[1:end - 12,:]
        y[14:end,:] = yt[2:end,:] - yt[1:end - 1,:]
        n_loss = 13
    elseif transfId == 14
        y = detrend(x, power=1)
        n_loss = 0
    elseif transfId == 15
        y = detrend(log.(x), power=1)
        n_loss = 0
    else
        error("Please specify a valid transformation index")
    end

    return (y = y, n_loss = n_loss)
end

function detrend(y::AbstractVecOrMat;power::Int = 1)
    T = size(y, 1)
    X = similar(y, T, power + 1)
    X[:,1] = ones(T)
    if power > 0
        for i = 1:power
            X[:,i + 1] = collect(1:T).^i
        end
    end
    b = X \ y
    return y - X * b
end