text:
#!/usr/bin/env julia

using Luxor, Colors

function simple_polys()
    gsave()
    sethue(0,0,0)
    vstep = 50
    hstep = 50
    x = 0
    y = 0
    translate(-currentwidth/2 + 50, -currentheight/2 + 50)
    for sides in 3:12
        ngon(x, y, 20, sides, 0, :stroke)
        x += hstep
        if x > currentwidth
            x = 0
            y += vstep
        end
    end
    translate(0, currentheight - 150)
    x = 0
    for sides in 3:12
        ngon(x, y, 20, sides, 0, :fill)
        x += hstep
        if x > currentwidth
            x = 0
            y += vstep
        end
    end
    grestore()
end

function hex_mixtures()
    gsave()
    vstep = 50
    hstep = 50
    for y in -currentheight:vstep:currentheight
        for x in -currentwidth:hstep:currentwidth
            setopacity(0.2)
            setline(rand()*10)
            sethue(backgroundcolors[rand(1:end)])
            gsave()
                 translate(x, y)
                 rotate(rand() * pi)
                 # point version
                 poly(hexagon(rand() * 5, rand() * 5, 15 + rand() * 5), :fill, close=true)
                 setopacity(0.1)
                 poly(hexagon(rand() * 10, rand() * 5, 15 + rand() * 5))
                 setopacity(0.3)
                 randomhue()
                 poly([
                        Point(rand() * -100,rand() * 100),
                        Point(rand() * -100,rand() * 100),
                        Point(rand() * -100,rand() * 100),
                        Point(rand() * -100,rand() * 100)
                     ], :fill)
            grestore()
            sethue(backgroundcolors[rand(1:end)])
            gsave()
                 translate(x,y)
                 ngon(x, y, 50, rand(3:13), 0, :stroke)
                 rotate(rand() * pi)
                 poly(hexagon(rand() * 5, rand() * 5, 15 + rand() * 5), :stroke)
                 ngon(x, y, 50, rand(3:13), 0, :fill)
                 rect(0, 0, hstep-5, vstep-5, :stroke)
                 ngon(x, y, 50, rand(3:13), 0, :none)
            grestore()
        end
    end
    grestore()
end

hexagon(x, y, size) = [Point(x + size * cos(2pi/6 * i), y + size * sin(2pi/6 * i)) for i in 1:6]

function polygon_test(fname)
    global currentwidth = 595 # pts
    global currentheight = 842 # pts
    Drawing(currentwidth, currentheight, fname)

    origin()
    background("antiquewhite")

    global foregroundcolors = diverging_palette(230, 280, 200, s = 0.99, b=0.8)
    global backgroundcolors = diverging_palette(200, 260, 280, s = 0.8, b=0.5)
    setline(2)
    # fill, then clip to heptagon
    setcolor("lightcyan")
    setline(3)
    ngon(0, 0, 270, 7, 0, :fillpreserve) # fill it
    sethue("orange")
    strokepreserve()                     # stroke it
    clip()                               # then use to clip
    hex_mixtures()
    clipreset()
    simple_polys()
    finish()
    println("finished test: output in $(fname)")
end

polygon_test("/tmp/polygon-test.pdf")


text:
"""
    struct GadgetSphere{T} <: AbstractGadgetGeometry
        center::Vector{T}
        radius::T
    end

Defines a sphere by `center` and `radius`.
To be used for [`read_particles_in_geometry`](@ref)
"""
struct GadgetSphere{T} <: AbstractGadgetGeometry
    center::Vector{T}
    radius::T
end

"""
    function get_geometry_center(sphere::GadgetSphere)

Returns the center of `sphere`.
"""
get_geometry_center(sphere::GadgetSphere) = sphere.center

"""
    get_geometry_box_corners(sphere::GadgetSphere)

Returns a tuple with the lower left and upper right corner of a box which contains the `sphere`.
"""
function get_geometry_box_corners(sphere::GadgetSphere)
    sphere.center .- sphere.radius, sphere.center .+ sphere.radius
end

"""
    get_geometry_mask(sphere::GadgetSphere, pos::Matrix{T}) where T

Returns the indices of all particles contained in the `sphere`.
"""
function get_geometry_mask(sphere::GadgetSphere, pos::Matrix{T}) where T
    r² = sphere.radius^2
    mask = @views @. (pos[1,:] - sphere.center[1])^2 + (pos[2,:] - sphere.center[2])^2 + (pos[3,:] - sphere.center[3])^2 ≤ r²

    return findall(mask)
end

ext: jl


text:
# ------------------------------------------------------------------------------------------
# # Character based RNN language model trained on 'The Complete Works of William
# Shakespeare'
# Based on http://karpathy.github.io/2015/05/21/rnn-effectiveness
# ------------------------------------------------------------------------------------------

# Load 'The Complete Works of William Shakespeare'
using Knet
include(Knet.dir("data","gutenberg.jl"))
trn,tst,chars = shakespeare()
map(summary,(trn,tst,chars))

# Print a sample
println(string(chars[trn[1020:1210]]...))

RNNTYPE = :lstm
BATCHSIZE = 256
SEQLENGTH = 100
INPUTSIZE = 168
VOCABSIZE = 84
HIDDENSIZE = 334
NUMLAYERS = 1
DROPOUT = 0.0
LR=0.001
BETA_1=0.9
BETA_2=0.999
EPS=1e-08
EPOCHS = 30;

# Minibatch data
function mb(a)
    N = div(length(a),BATCHSIZE)
    x = reshape(a[1:N*BATCHSIZE],N,BATCHSIZE)' # reshape full data to (B,N) with contiguous rows
    minibatch(x[:,1:N-1], x[:,2:N], SEQLENGTH) # split into (B,T) blocks
end
dtrn,dtst = mb(trn),mb(tst)
map(length, (dtrn,dtst))

# Define model
function initmodel()
    w(d...)=KnetArray(xavier(Float32,d...))
    b(d...)=KnetArray(zeros(Float32,d...))
    r,wr = rnninit(INPUTSIZE,HIDDENSIZE,rnnType=RNNTYPE,numLayers=NUMLAYERS,dropout=DROPOUT)
    wx = w(INPUTSIZE,VOCABSIZE)
    wy = w(VOCABSIZE,HIDDENSIZE)
    by = b(VOCABSIZE,1)
    return r,wr,wx,wy,by
end;

# Given the current character, predict the next character
function predict(ws,xs,hx,cx;pdrop=0)
    r,wr,wx,wy,by = ws
    x = wx[:,xs]                                    # xs=(B,T) x=(X,B,T)
    x = dropout(x,pdrop)
    y,hy,cy = rnnforw(r,wr,x,hx,cx,hy=true,cy=true) # y=(H,B,T) hy=cy=(H,B,L)
    y = dropout(y,pdrop)
    y2 = reshape(y,size(y,1),size(y,2)*size(y,3))   # y2=(H,B*T)
    return wy*y2.+by, hy, cy
end

# Define loss and its gradient
function loss(w,x,y,h;o...)
    py,hy,cy = predict(w,x,h...;o...)
    h[1],h[2] = getval(hy),getval(cy)
    return nll(py,y)
end

lossgradient = gradloss(loss);

function train(model,data,optim)
    hiddens = Any[nothing,nothing]
    Σ,N=0,0
    for (x,y) in data
        grads,loss1 = lossgradient(model,x,y,hiddens;pdrop=DROPOUT)
        update!(model, grads, optim)
        Σ,N=Σ+loss1,N+1
    end
    return Σ/N
end;

function test(model,data)
    hiddens = Any[nothing,nothing]
    Σ,N=0,0
    for (x,y) in data
        Σ,N = Σ+loss(model,x,y,hiddens), N+1
    end
    return Σ/N
end;

# Train model or load from file if exists
using JLD
model=optim=nothing; knetgc()
if !isfile("shakespeare.jld")
    model = initmodel()
    optim = optimizers(model, Adam; lr=LR, beta1=BETA_1, beta2=BETA_2, eps=EPS);    info("Training...")
    @time for epoch in 1:EPOCHS
        @time trnloss = train(model,dtrn,optim) # ~18 seconds
        @time tstloss = test(model,dtst)        # ~0.5 seconds
        println((:epoch, epoch, :trnppl, exp(trnloss), :tstppl, exp(tstloss)))
    end
    save("shakespeare.jld","model",model)
else
    model = load("shakespeare.jld","model")
end
summary(model)

# Sample from trained model
function generate(model,n)
    function sample(y)
        p,r=Array(exp.(y-logsumexp(y))),rand()
        for j=1:length(p); (r -= p[j]) < 0 && return j; end
    end
    h,c = nothing,nothing
    x = findfirst(chars,'\n')
    for i=1:n
        y,h,c = predict(model,[x],h,c)
        x = sample(y)
        print(chars[x])
    end
    println()
end;

generate(model,1000)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using Documenter
using DiffPrivacy

makedocs(
    modules = [DiffPrivacy],
    sitename = "DiffPrivacy"
)

deploydocs(
    repo = "github.com/r0cketr1kky/DiffPrivacy.jl.git",
    deps = Deps.pip("pygments", "mkdocs", "mkdocs-material", "python-markdown-math"),
    make   = () -> run(`mkdocs build`),
    target = "site"
)


text:
"Bohr radius"
const a0 = 0.52917721067

"1 hartree energy = 27.211386245988eV"
const hartree = 27.211386245988


text:
#
# Layout combinators, tiling, and rendering.
#

# Layout tree.

abstract type AbstractBlock end

struct Layout
    blk::AbstractBlock
    args::Vector{Layout}
end

struct LiteralBlock <: AbstractBlock
    atom::Any
    len::Int
end

struct HorizontalBlock <: AbstractBlock
end

struct VerticalBlock <: AbstractBlock
end

struct ChoiceBlock <: AbstractBlock
end

struct PenaltyBlock <: AbstractBlock
    cost::Int
    break_factor::Int
    spill_factor::Int

    PenaltyBlock(; cost::Int=0, break_factor::Int=1, spill_factor::Int=1) =
        new(cost, break_factor, spill_factor)
end

let NO_ARGS = Layout[],
    INDENTS = Layout[Layout(LiteralBlock(" " ^ i, i), NO_ARGS) for i = 1:DEFAULT_LINE_WIDTH]

    global literal, indent

    literal(str::String) =
        Layout(LiteralBlock(str, textwidth(str)), NO_ARGS)

    literal(str::Union{AbstractString,Symbol}) =
        literal(string(str))

    literal(atom, len::Int) =
        Layout(LiteralBlock(atom, len), NO_ARGS)

    indent(i::Int) =
        1 <= i <= length(INDENTS) ? INDENTS[i] : Layout(LiteralBlock(" " ^ i, i), NO_ARGS)
end

const ZERO = literal("")

# Visualization.

show(io::IO, lt::Layout) =
    pprint(io, tile(lt))

# Layout combinators.

(*)(lt::Layout, lts::Layout...) =
    let lts = filter(!isequal(ZERO), collect(Layout, (lt, lts...)))
        if isempty(lts)
            ZERO
        elseif length(lts) == 1
            lts[1]
        else
            Layout(HorizontalBlock(), lts)
        end
    end

(/)(lt::Layout, lts::Layout...) =
    Layout(VerticalBlock(), collect(Layout, (lt, lts...)))

(|)(lt::Layout, lts::Layout...) =
    Layout(ChoiceBlock(), collect(Layout, (lt, lts...)))

penalize(lt::Layout; cost::Int=0, break_factor::Int=1, spill_factor::Int=1) =
    Layout(PenaltyBlock(cost=cost, break_factor=break_factor, spill_factor=spill_factor), Layout[lt])

nobreak(lt::Layout) =
    let lt′ = nobreak(lt.blk, lt.args)
        lt′ isa Layout && lt′.blk == lt.blk && lt′.args == lt.args ? lt : lt′
    end

nobreak(blk::LiteralBlock, args::Vector{Layout}) =
    Layout(blk, args)

function nobreak(blk::HorizontalBlock, args::Vector{Layout})
    args′ = Layout[]
    for arg in args
        arg′ = nobreak(arg)
        arg′ !== nothing || return nothing
        push!(args′, arg′)
    end
    Layout(blk, args′)
end

nobreak(::VerticalBlock, ::Vector{Layout}) = nothing

function nobreak(blk::Union{ChoiceBlock,PenaltyBlock}, args::Vector{Layout})
    args′ = Layout[]
    for arg in args
        arg′ = nobreak(arg)
        arg′ !== nothing || continue
        push!(args′, arg′)
    end
    !isempty(args′) || return nothing
    Layout(blk, args′)
end

function list_layout(items::Vector{Layout};
                     prefix::Union{String,Symbol}="",
                     par::Tuple{String,String}=("(", ")"),
                     sep::String=", ",
                     sep_brk=:end,  # :start, :end, :both, or :none
                     tab::Int=4,
                     nobrk::Int=10)
    !isempty(items) || return literal("$(prefix)$(par[1])$(par[2])")
    head_lt = literal("$(prefix)$(par[1])")
    tail_lt = literal("$(par[2])")
    sepc_lt = literal(sep)
    sepl_lt = literal(sep_brk == :start || sep_brk == :both ? lstrip(sep) : "")
    sepr_lt = literal(sep_brk == :end || sep_brk == :both ? rstrip(sep) : "")
    tab_lt = indent(tab)
    vlt = items[1]
    #hlt = penalize(items[1], break_factor=nobrk)
    hlt = items[1]
    for item in items[2:end]
        vlt = (vlt * sepr_lt) / (sepl_lt * item)
        hlt = hlt * sepc_lt * item
    end
    vlt = (head_lt | (head_lt / tab_lt)) * vlt * tail_lt
    #hlt = penalize(head_lt * hlt * tail_lt, break_factor=nobrk)
    #vlt | hlt
    hlt = nobreak(head_lt * hlt * tail_lt)
    (hlt !== nothing ? vlt | hlt : vlt)::Layout
end

function pair_layout(fst::Layout,
                     snd::Layout;
                     sep::String=" => ",
                     sep_brk=:end,  # :start, :end, :both, or :none
                     tab::Int=4)
    sepc_lt = literal(sep)
    sepl_lt = literal(sep_brk == :start || sep_brk == :both ? lstrip(sep) : "")
    sepr_lt = literal(sep_brk == :end || sep_brk == :both ? rstrip(sep) : "")
    tab_lt = indent(tab)
    ((nobreak(fst) * sepc_lt) | ((fst * sepr_lt) / (tab_lt * sepl_lt))) * snd
end

# Fallback layout.

tile(obj) = literal(repr(obj))

# Layouts for standard types.

tile(p::Pair) =
    pair_layout(tile(p.first), tile(p.second))

tile(t::Tuple) =
    if length(t) == 1
        literal("(") * tile(t[1]) * literal(",)")
    else
        list_layout(Layout[tile(x) for x in t])
    end

tile(v::Vector) =
    list_layout(Layout[tile(x) for x in v], par=("[", "]"))

tile(t::NamedTuple) =
    if length(t) == 1
        ((key, val),) = pairs(t)
        literal("(") * pair_layout(literal(key), tile(val), sep=" = ") * literal(",)")
    else
        list_layout(Layout[pair_layout(literal(key), tile(val), sep=" = ")
                           for (key, val) in pairs(t)])
    end

tile(d::Dict) =
    list_layout(Layout[pair_layout(tile(key), tile(val)) for (key, val) in d],
                prefix=:Dict)

# Layout of the layout tree.

tile(lt::Layout) =
    tile_layout(lt)

tile_layout(lt::Layout, precedence::Int=0) =
    tile_layout(lt.blk, lt.args, precedence)

function tile_layout(blk::LiteralBlock, ::Vector{Layout}, ::Int)
    if blk.atom isa String && blk.len == textwidth(blk.atom)
        if blk.len > 0 && all(ch == ' ' for ch in blk.atom)
            literal("indent($(blk.len))")
        else
            literal("literal($(repr(blk.atom)))")
        end
    else
        literal("literal($(repr(blk.atom)), $(blk.len))")
    end
end

function tile_layout(blk::PenaltyBlock, args::Vector{Layout}, ::Int)
    items = Layout[tile_layout(arg) for arg in args]
    if blk.cost != 0
        push!(items, literal("cost=$(blk.cost)"))
    end
    if blk.break_factor != 1
        push!(items, literal("break_factor=$(blk.break_factor)"))
    end
    if blk.spill_factor != 1
        push!(items, literal("spill_factor=$(blk.spill_factor)"))
    end
    list_layout(items, prefix=:penalize)
end

function tile_layout(blk::Union{HorizontalBlock, VerticalBlock, ChoiceBlock}, args::Vector{Layout}, precedence::Int)
    precedence′ = blk isa ChoiceBlock ? 0 : 1
    items = Layout[]
    first = true
    for arg in args
        item = tile_layout(arg, first ? precedence′ : precedence′ + 1)
        push!(items, item)
        first = false
    end
    par = precedence > precedence′ ? ("(", ")") : ("", "")
    sep = blk isa HorizontalBlock ? " * " :
          blk isa VerticalBlock ? " / " :
          blk isa ChoiceBlock ? " | " : ""
    return list_layout(items, par=par, sep=sep)
end

# Rendering

render(io::IO, lt::Layout, col::Int=0) =
    render(io, lt.blk, lt.args, col)

function render(io::IO, blk::LiteralBlock, ::Vector{Layout}, col::Int)
    print(io, blk.atom)
    col + blk.len
end

function render(io::IO, ::HorizontalBlock, args::Vector{Layout}, col::Int)
    for arg in args
        col = render(io, arg, col)
    end
    col
end

function render(io::IO, ::VerticalBlock, args::Vector{Layout}, col::Int)
    nl = "\n" * " " ^ col
    col′ = col
    first = true
    for arg in args
        if !first
            print(io, nl)
        end
        col′ = render(io, arg, col)
        first = false
    end
    col′
end

render(io::IO, ::Union{ChoiceBlock, PenaltyBlock}, args::Vector{Layout}, col::Int) =
    !isempty(args) ? render(io, args[1], col) : col


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".

#TODO: Use full JPL ephemeris for high precision
function _planet_coords(date::T, num::Integer) where {T<:AbstractFloat}
    rad, long, lat = helio(date, num, true)
    rade, longe, late = helio(date, 3, true)
    sin_lat, cos_lat = sincos(lat)
    sin_long, cos_long = sincos(long)
    sin_late, cos_late = sincos(late)
    sin_longe, cos_longe = sincos(longe)
    x = rad * cos_lat * cos_long - rade * cos_late * cos_longe
    y = rad * cos_lat * sin_long - rade * cos_late * sin_longe
    z = rad * sin_lat - rade * sin_late
    lamda = rad2deg(atan(y,x))
    beta = rad2deg(atan(z, hypot(x,y)))
    ra, dec = euler(lamda, beta, 4)
    return ra, dec
end

"""
    planet_coords(date, num)

### Purpose ###

Find right ascension and declination for the planets when provided a date as input.

### Explanation ###

This function uses the [`helio`](@ref) to get the heliocentric ecliptic coordinates of the
planets at the given date which it then converts these to geocentric ecliptic
coordinates ala "Astronomical Algorithms" by Jean Meeus (1991, p 209).
These are then converted to right ascension and declination using [`euler`](@ref).

The accuracy between the years 1800 and 2050 is better than 1 arcminute for the
terrestial planets, but reaches 10 arcminutes for Saturn. Before 1850 or after 2050
the accuracy can get much worse.

### Arguments ###

* `date`: Can be either a single date or an array of dates. Each element can be
  either a `DateTime` type or Julian Date. It can be a scalar or vector.
* `num`: integer denoting planet number, scalar or vector
  1 = Mercury, 2 = Venus, ... 9 = Pluto. If not in that change, then the
  program will throw an error.

### Output ###

* `ra`: right ascension of planet(J2000), in degrees
* `dec`: declination of the planet(J2000), in degrees

### Example ###

Find the RA, Dec of Venus on 1992 Dec 20

```jldoctest
julia> using AstroLib, Dates

julia> adstring(planet_coords(DateTime(1992,12,20),2))
" 21 05 02.8  -18 51 41"
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
planet_coords(date::Real, num::Integer) = _planet_coords(float(date), num)

function planet_coords(date::AbstractVector{R},
                       num::AbstractVector{<:Integer}) where {R<:Real}
    @assert length(date) == length(num) "date and num arrays should be of the same length"
    typedate = float(R)
    ra_out  = similar(date, typedate)
    dec_out = similar(date, typedate)
    for i in eachindex(date)
        ra_out[i], dec_out[i] = planet_coords(date[i], num[i])
    end
    return ra_out, dec_out
end

planet_coords(date::DateTime, num::Integer) = planet_coords(jdcnv(date), num)


text:
using SimpleNL, Ipopt, JuMP, AmplNLWriter, StatsPlots

include("benchmark_include.jl")

result = []

for (N,nlm,jump,name) in [
    (10000,nlm_luksan_vlcek_501,jump_luksan_vlcek_501,"luksanvlcek"),
    (10000,nlm_ocp,jump_ocp,"hehnandrea")
]
    for i=1:2
        GC.enable(false)
        t1 = @elapsed begin
            m = nlm(;N=N,optimizer=Ipopt.Optimizer,output_file="output/$name-simplenlmodels.out")
            SimpleNL.optimize!(m)
        end
        t2,t3 = parse_ipopt_output("output/$name-simplenlmodels.out")
        GC.enable(true);GC.gc()

        GC.enable(false)
        m = jump(;N=N,optimizer=() -> AmplNLWriter.Optimizer("ipopt",["output_file=output/$name-ampl.out"]))
        JuMP.optimize!(m)
        t4 = m.moi_backend.optimizer.model.ext[:MPBModel].inner.solve_time
        t5,t6 = parse_ipopt_output("output/$name-ampl.out")
        GC.enable(true);GC.gc()

        GC.enable(false)
        t7 = @elapsed begin
            m = jump(;N=N,optimizer=()->Ipopt.Optimizer(output_file="output/$name-jump.out"))
            JuMP.optimize!(m)
        end
        t8,t9 = parse_ipopt_output("output/$name-jump.out")
        GC.enable(true);GC.gc()

        i==2 && push!(result,[[t1,t4,t7],[t2,t5,t8],[t3,t6,t9]])
    end
end

ticklabel = ["SimpleNL","AMPL","JuMP"]
barlabel = ["Model Creation + Solution Time (Wall)" "Solution Time (CPU)" "NLP Function Evaluation Time (CPU)"]

plt1 = groupedbar([result[1][1] result[1][2] result[1][3]], title="Luksan & Vlcek (1999)", bar_position = :dodge, legend=:topleft,
                  bar_width=0.7, xticks=(1:3, ticklabel), label=barlabel, ylabel="Time (sec)", framestyle=:box, ylim = (0,1.3*maximum(vcat(result[1]...))),dpi = 300)
plt2 = groupedbar([result[2][1] result[2][2] result[2][3]], title="Hehn & Andrea (2011, modified)", bar_position = :dodge, legend=:topleft,
                  bar_width=0.7, xticks=(1:3, ticklabel), label=barlabel, ylabel="Time (sec)", framestyle=:box, ylim = (0,1.3*maximum(vcat(result[2]...))),dpi = 300)
savefig(plt1,"output/luksanvlcek.png")
savefig(plt2,"output/hehnandrea.png")


text:
# redis client in julia
# wraps the stable hiredis C library

module HiRedis

using Docile
using Logging

Logging.configure(level=WARNING)
@docstrings

const REDIS_ERR = -1
const REDIS_OK = 0

const REDIS_REPLY_STRING = 1
const REDIS_REPLY_ARRAY = 2
const REDIS_REPLY_INTEGER = 3
const REDIS_REPLY_NIL = 4
const REDIS_REPLY_STATUS = 5
const REDIS_REPLY_ERROR = 6

global redisContext = 0
global pipelinedCommandCount = 0

type RedisReadTask
    rtype::Int32
    elements::Int32
    idx::Int32
    obj::Ptr{Void}
    parent::Ptr{RedisReadTask}
    privdata::Ptr{Void}
end

function create_string(task::Ptr{RedisReadTask}, str::Ptr{Uint8}, len::Uint)
    # not implemented
    ret::Ptr{Void} = 0
    ret
end

function create_array(task::Ptr{RedisReadTask}, len::Int32)
    # not implemented
    ret::Ptr{Void} = 0
    ret
end

function create_integer(task::Ptr{RedisReadTask}, int::Int)
    # not implemented
    ret::Ptr{Void} = 0
    ret
end

function create_nil(task::Ptr{RedisReadTask})
    # not implemented
    ret::Ptr{Void} = 0
    ret
end

function free_object(obj::Ptr{Void})
    # not implemented
    ret::Void = 0
    ret
end

const create_string_c = cfunction(create_string, Ptr{Void}, (Ptr{RedisReadTask}, Ptr{Uint8}, Uint))

const create_array_c = cfunction(create_array, Ptr{Void}, (Ptr{RedisReadTask}, Int32))

const create_integer_c = cfunction(create_integer, Ptr{Void}, (Ptr{RedisReadTask}, Int))

const create_nil_c = cfunction(create_nil, Ptr{Void}, (Ptr{RedisReadTask},))

const free_object_c = cfunction(free_object, Void, (Ptr{Void},))

type RedisReplyObjectFunctions
    create_string_c
    create_array_c
    create_integer_c
    create_nil_c
    free_object_c
end

type RedisReader
    err::Int32
    errstr::Ptr{Uint8}
    buf::Ptr{Uint8}
    pos::Uint
    len::Uint
    maxbuf::Uint
    rstack::Array{RedisReadTask,1}
    ridx::Int32
    reply::Ptr{Void}
    fn::Ptr{RedisReplyObjectFunctions}
    privdata::Ptr{Void}
end

type RedisContext
    err::Int32
    errstr::Ptr{Uint8}
    fd::Int32
    flags::Int32
    obuf::Ptr{Uint8}
    reader::Ptr{RedisReader}
end

type RedisReply
    rtype::Int32                  # REDIS_REPLY_*
    integer::Uint64               # The integer when type is REDIS_REPLY_INTEGER
    len::Int32                    # Length of string
    str::Ptr{Uint8}               # Used for both REDIS_REPLY_ERROR and REDIS_REPLY_STRING
    elements::Uint                # number of elements, for REDIS_REPLY_ARRAY
    element::Ptr{Ptr{RedisReply}} # elements vector for REDIS_REPLY_ARRAY
end

function start_session(host::ASCIIString = "127.0.0.1", port::Int = 6379)
    global redisContext = ccall((:redisConnect, "libhiredis"), Ptr{RedisContext}, (Ptr{Uint8}, Int32), host, port)
end

function end_session()
    if redisContext != 0 # isdefined(:redisContext)
        ccall((:redisFree, "libhiredis"), Void, (Ptr{RedisContext},), redisContext::Ptr{RedisContext})
    end
end

@doc "Free memory allocated to objects returned from hiredis" ->
function free_reply_object(redisReply)
    ccall((:freeReplyObject, "libhiredis"), Void, (Ptr{RedisReply},), redisReply)
end

@doc """
Appends commands to an output buffer. Pipelining is sending a batch of commands
to redis to be processed in bulk. It cuts down the number of network requests.
""" ->
function pipeline_command(command::ASCIIString)
    if redisContext == 0 # !isdefined(:redisContext)
        start_session()
    end
    debug(string("RedisClient.pipeline_command: ", command))
    global pipelinedCommandCount += 1
    ccall((:redisAppendCommand, "libhiredis"), Int32, (Ptr{RedisContext}, Ptr{Uint8}), redisContext::Ptr{RedisContext}, command)
end

@doc """
In a blocking context, this function first checks if there are unconsumed
replies to return and returns one if so. Otherwise, it flushes the output
buffer to the socket and reads until it has a reply.
""" ->
function call_get_reply(redisReply::Array{Ptr{RedisReply},1})
    ccall((:redisGetReply, "libhiredis"), Int32, (Ptr{RedisContext}, Ptr{Ptr{RedisReply}}), redisContext::Ptr{RedisContext}, redisReply)
end

@doc """
Calls call_get_reply until the pipelinedCommandCount is 0 or an error
is returned. Adds the results from each reply to an Array, then returns the
Array.
""" ->
function get_reply()
    redisReply = Array(Ptr{RedisReply},1)
    results = Any[]
    while pipelinedCommandCount::Int > 0 && call_get_reply(redisReply) == REDIS_OK
        push!(results, get_result(redisReply[1]))
        global pipelinedCommandCount -= 1
    end
    results
end

@doc """
Converts the reply object from hiredis into a String, int, or Array
as appropriate the the reply type.
""" ->
function get_result(redisReply::Ptr{RedisReply})
    r = unsafe_load(redisReply)
    if r.rtype == REDIS_REPLY_ERROR
        error(bytestring(r.str))
    end
    ret::Any = nothing
    if r.rtype == REDIS_REPLY_STRING
        ret = bytestring(r.str)
    elseif r.rtype == REDIS_REPLY_INTEGER
        ret = int(r.integer)
    elseif r.rtype == REDIS_REPLY_ARRAY
        n = int(r.elements)
        results = ASCIIString[]
        replies = pointer_to_array(r.element, n)
        for i in 1:n
            ri = unsafe_load(replies[i])
            push!(results, bytestring(ri.str))
        end
        ret = results
    end
    free_reply_object(redisReply)
    ret
end

@doc "Pipelines a block of ordinary blocking calls." ->
macro pipeline(expr::Expr)
    Expr(:block, map(x ->
        begin
            if x.args[1] in names(HiRedis)
                args = copy(x.args)
                push!(args, Expr(:kw, :pipeline, true))
                Expr(x.head, args...)
            else
                x
            end
        end, filter(x -> typeof(x) == Expr, expr.args))...)
end

@doc "Issues a blocking command to hiredis." ->
function do_command(command::ASCIIString)
    if redisContext == 0 # !isdefined(:redisContext)
#         error("redisContext not defined. Please call RedisClient.start_session.")
        start_session()
    end
    debug(string("RedisClient.do_command: ", command))
    redisReply = ccall((:redisvCommand, "libhiredis"), Ptr{RedisReply}, (Ptr{RedisContext}, Ptr{Uint8}), redisContext::Ptr{RedisContext}, command)
    get_result(redisReply)
end

@doc "Issues a blocking command to hiredis, accepting command arguments as an Array." ->
function do_command{S<:Any}(argv::Array{S,1})
    if redisContext == 0 # !isdefined(:redisContext)
        start_session()
    end
    redisReply = ccall((:redisCommandArgv, "libhiredis"), Ptr{RedisReply}, (Ptr{RedisContext}, Int32, Ptr{Ptr{Uint8}}, Ptr{Uint}), redisContext::Ptr{RedisContext}, length(argv), argv, C_NULL)
    get_result(redisReply)
end

@doc "Switches between blocking and pipelined command execution according to flag." ->
function docommand(cmd::ASCIIString, pipeline::Bool)
    (pipeline || (pipelinedCommandCount::Int > 0)) ? pipeline_command(cmd) : do_command(cmd)
end

# @doc "Set the string value of a key." Dict{Any,Any}(
#     :params => Dict(
#         :key           => "The key to set",
#         :value         => "The value to set",
#         :ex            => "Set the specified expire time, in seconds",
#         :px            => "Set the specified expire time, in milliseconds",
#         :pipeline      => "A flag to indicate that the command should be pipelined"
#     )) ->
@doc md"""
Set the string value of a key.

Params:

* key - The key to set
* value - The value to set
* ex - Set the specified expire time, in seconds
* px - Set the specified expire time, in milliseconds
* pipeline - A flag to indicate that the command should be pipelined
""" ->
function kvset(key::ASCIIString, value::Any; ex::Int=0, px::Int=0, pipeline::Bool=false)
    cmd = string("SET ", key, " ", value)
    if ex > 0
        cmd = string(cmd, " EX ", ex)
    elseif px > 0
        cmd = string(cmd, " PX ", ex)
    end
    docommand(cmd, pipeline)
end

@doc "Get the value of a key." ->
function kvget(key::ASCIIString; pipeline::Bool=false)
    docommand(string("GET ", key), pipeline)
end

@doc "Increment the integer value of a key by one." ->
function incr(key::ASCIIString; pipeline::Bool=false)
    docommand(string("INCR ", key), pipeline)
end

@doc "Increment the integer value of a key by the given amount." ->
function incrby(key::ASCIIString, by::Int; pipeline::Bool=false)
    docommand(string("INCRBY ", key, " ", by), pipeline)
end

@doc "Delete a key." ->
function del(key::ASCIIString)
    docommand(string("DEL ", key))
end

@doc "Determine if a key exists." ->
function exists(key::ASCIIString)
    docommand(string("EXISTS ", key))
end

@doc "Find all keys matching the given pattern." ->
function getkeys(pattern::ASCIIString)
    docommand(string("KEYS ", pattern))
end

@doc "Return a serialized version of the value stored at the specified key." ->
function rdump(key::ASCIIString)
    docommand(string("DUMP ", key))
end

@doc "Determine the type stored at key." ->
function rtype(key::ASCIIString)
    docommand(string("TYPE ", key))
end

@doc "Set the string value of a hash field." ->
function hset(key::ASCIIString, attr_name::ASCIIString, attr_value::Any; pipeline::Bool=false)
    #TODO do_command(["HSET %s %s %s", key, attr_name, string(attr_value)])
    docommand(string("HSET ", key, " ", attr_name, " ", attr_value), pipeline)
end

@doc "Get the value of a hash field." ->
function hget(key::ASCIIString, attr_name::ASCIIString; pipeline::Bool=false)
    docommand(string("HGET ", key, " ", attr_name), pipeline)
end

@doc """
Set multiple hash fields to multiple values. A variable number of arguments
follow key in `field` `value` format, e.g.:
    `hmset("myhash", "field1", "Hello", "field2", "World")`
""" ->
function hmset(key::ASCIIString, argv::Any...; pipeline::Bool=false)
    cmd = string("HMSET ", key)
    for arg in argv
        cmd = string(cmd, " ", string(arg))
    end
    docommand(cmd, pipeline)
end

@doc """
Set multiple hash fields to multiple values. Fields and values are provided
as an Array.
""" ->
function hmset{S<:Any}(key::ASCIIString, attrs::Array{S,1}; pipeline::Bool=false)
    cmd = string("HMSET ", key)
    for attr in attrs
        cmd = string(cmd, " ", string(attr))
    end
    docommand(cmd, pipeline)
end

@doc """
Set multiple hash fields to multiple values. Fields and values are provided
as a Dict.
""" ->
function hmset(key::ASCIIString, attrs::Dict{ASCIIString,Any}; pipeline::Bool=false)
    cmd = string("HMSET ", key)
    for (field, val) in attrs
        cmd = string(cmd, " ", field, " ", string(val))
    end
    docommand(cmd, pipeline)
end

@doc """
Get the values of all the given hash fields. Multiple fields are provided
as additional arguments.
""" ->
function hmget(key::ASCIIString, argv::ASCIIString...; pipeline::Bool=false)
    cmd = string("HMGET ", key)
    for arg in argv
        cmd = string(cmd, " ", arg)
    end
    docommand(cmd, pipeline)
end

@doc """
Get the values of all the given hash fields. Multiple fields are provided
as an Array.
""" ->
function hmget(key::ASCIIString, fields::Array{ASCIIString,1}; pipeline::Bool=false)
    cmd = string("HMGET ", key)
    for field in fields
        cmd = string(cmd, " ", field)
    end
    docommand(cmd, pipeline)
end

@doc "Get all the fields and values in a hash." ->
function hgetall(key::ASCIIString)
    reply::Array{Any,1} = do_command(string("HGETALL ", key))
    n = length(reply)
    dict = Dict{ASCIIString,Any}()
    if n > 1 && mod(n, 2) == 0
        for i = 1:2:n
            dict[reply[i]] = reply[i + 1]
        end
    end
    dict
end

@doc """
Delete one or more hash fields. Multiple fields are provided
as additional arguments.
""" ->
function hdel(key::ASCIIString, argv::ASCIIString...; pipeline::Bool=false)
    cmd = string("HDEL ", key)
    for arg in argv
        cmd = string(cmd, " ", arg)
    end
    docommand(cmd, pipeline)
end

@doc """
Delete one or more hash fields. Multiple fields are provided
as an Array.
""" ->
function hdel(key::ASCIIString, fields::Array{ASCIIString,1}; pipeline::Bool=false)
    cmd = string("HDEL ", key)
    for field in fields
        cmd = string(cmd, " ", field)
    end
    docommand(cmd, pipeline)
end

@doc "Determine if a hash field exists." ->
function hexists(key::ASCIIString, field::ASCIIString; pipeline::Bool=false)
    docommand(string("HEXISTS ", key, " ", field), pipeline)
end

@doc "Get all the fields in a hash." ->
function hkeys(key::ASCIIString; pipeline::Bool=false)
    docommand(string("HKEYS ", key), pipeline)
end

@doc "Get all the values in a hash." ->
function hvals(key::ASCIIString; pipeline::Bool=false)
    docommand(string("HVALS ", key), pipeline)
end

@doc "Get the number of fields in a hash." ->
function hlen(key::ASCIIString; pipeline::Bool=false)
    docommand(string("HLEN ", key), pipeline)
end

@doc "Increment the integer value of a hash field by the given number." ->
function hincrby(key::ASCIIString, field::ASCIIString, increment::Int; pipeline::Bool=false)
    docommand(string("HINCRBY ", key, " ", field, " ", increment), pipeline)
end

@doc "Increment the float value of a hash field by the given number." ->
function hincrby(key::ASCIIString, field::ASCIIString, increment::Float64; pipeline::Bool=false)
    docommand(string("HINCRBYFLOAT ", key, " ", field, " ", increment), pipeline)
end

@doc "Add one or more members to a set." ->
function sadd(key::ASCIIString, argv::Any...; pipeline::Bool=false)
    cmd = string("SADD ", key)
    for arg in argv
        if isa(arg, Array) || isa(arg, Range)
            for member in arg
                cmd = string(cmd, " ", string(member))
            end
        else
            cmd = string(cmd, " ", string(arg))
        end
    end
    docommand(cmd, pipeline)
end

@doc "Get all the members in a set." ->
function smembers(key::ASCIIString; pipeline::Bool=false)
    docommand(string("SMEMBERS ", key), pipeline)
end

@doc "Determine if a given value is a member of a set." ->
function sismember(key::ASCIIString, member::Any; pipeline::Bool=false)
    docommand(string("SISMEMBER ", key, " ", string(member)), pipeline)
end

@doc "Get the number of members in a set." ->
function scard(key::ASCIIString; pipeline::Bool=false)
    docommand(string("SCARD ", key), pipeline)
end

@doc """
Remove one or more members from a set. Members are provided
as additional arguments.
""" ->
function srem(key::ASCIIString, argv::Any...; pipeline::Bool=false)
    cmd = string("SREM ", key)
    for arg in argv
        cmd = string(cmd, " ", string(arg))
    end
    docommand(cmd, pipeline)
end

@doc """
Remove one or more members from a set. Members are provided
as an Array.
""" ->
function srem{S<:Any}(key::ASCIIString, members::Array{S,1}; pipeline::Bool=false)
    cmd = string("SREM ", key)
    for member in members
        cmd = string(cmd, " ", string(member))
    end
    docommand(cmd, pipeline)
end

@doc "Subtract multiple sets. Multiple sets are provided as additional arguments." ->
function sdiff(key::ASCIIString, argv::ASCIIString...; pipeline::Bool=false)
    cmd = string("SDIFF ", key)
    for arg in argv
        cmd = string(cmd, " ", arg)
    end
    docommand(cmd, pipeline)
end

@doc "Subtract multiple sets. Multiple sets are provided as an Array." ->
function sdiff(keys::Array{ASCIIString,1}; pipeline::Bool=false)
    cmd = "SDIFF"
    for key in keys
        cmd = string(cmd, " ", key)
    end
    docommand(cmd, pipeline)
end

@doc "Intersect multiple sets. Multiple sets are provided as additional arguments." ->
function sinter(key::ASCIIString, argv::ASCIIString...; pipeline::Bool=false)
    cmd = string("SINTER ", key)
    for arg in argv
        cmd = string(cmd, " ", arg)
    end
    docommand(cmd, pipeline)
end

@doc "Intersect multiple sets. Multiple sets are provided as an Array." ->
function sinter(keys::Array{ASCIIString,1}; pipeline::Bool=false)
    cmd = "SINTER"
    for key in keys
        cmd = string(cmd, " ", key)
    end
    docommand(cmd, pipeline)
end

@doc "Add multiple sets. Multiple sets are provided as additional arguments." ->
function sunion(key::ASCIIString, argv::ASCIIString...; pipeline::Bool=false)
    cmd = string("SUNION ", key)
    for arg in argv
        cmd = string(cmd, " ", arg)
    end
    docommand(cmd, pipeline)
end

@doc "Add multiple sets. Multiple sets are provided as an Array." ->
function sunion(keys::Array{ASCIIString,1}; pipeline::Bool=false)
    cmd = "SUNION"
    for key in keys
        cmd = string(cmd, " ", key)
    end
    docommand(cmd, pipeline)
end

@doc "Change the selected database for the current connection." ->
function selectdb(db::Int)
    do_command(string("SELECT ", db))
end

@doc "Remove all keys from the current database." ->
function flushdb()
    do_command("FLUSHDB")
end

@doc "Remove all keys from all databases." ->
function flushall()
    do_command("FLUSHALL")
end

export start_session, end_session,                                                  # session
    kvset, kvget, incr, incrby, del, exists, getkeys, rdump, rtype,                 # key-value
    hset, hget, hmset, hmget, hgetall, hdel, hexists, hkeys, hvals, hlen, hincrby,  # hash sets
    sadd, smembers, sismember, scard, srem, sdiff, sinter, sunion,                  # sets
    selectdb, flushdb, flushall,                                                    # management
    @pipeline, get_reply,                                                           # pipelining
    do_command, pipeline_command                                                    # generic

end # module


text:
# # Illustration of Jacobi Fields
#
# This tutorial illustrates the usage of Jacobi Fields within
# `Manopt.jl`.
# For this tutorial you should be familiar with the basic terminology on a
# manifold like the exponential and logarithmic map as well as
# [shortest geodesic](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.shortest_geodesic-Tuple{Manifold,Any,Any})s.
#
# We first initialize the manifold
exportFolder = joinpath(@__DIR__,"..","..","docs","src","assets","images","tutorials") #src
using Manopt, Manifolds
# and we define some colors from [Paul Tol](https://personal.sron.nl/~pault/)
using Colors
black = RGBA{Float64}(colorant"#000000")
TolVibrantOrange = RGBA{Float64}(colorant"#EE7733")
TolVibrantCyan = RGBA{Float64}(colorant"#33BBEE")
TolVibrantTeal = RGBA{Float64}(colorant"#009988")
nothing #hide
# Assume we have two points on the equator of the [Sphere](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html) $\mathcal M = \mathbb S^2$
M = Sphere(2)
p,q = [ [1.,0.,0.], [0.,1.,0.]]
# their connecting [shortest geodesic](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.shortest_geodesic-Tuple{Manifold,Any,Any}) (sampled at `100` points)
geodesicCurve = shortest_geodesic(M,p,q,[0:0.1:1.0...]);
nothing #hide
# looks as follows using the [`asymptote_export_S2_signals`](@ref) export
asymptote_export_S2_signals(exportFolder*"/jacobiGeodesic.asy"; #src
    curves=[geodesicCurve], points = [ [p,q] ], #src
    colors=Dict(:curves => [black], :points => [TolVibrantOrange]), #src
    dotSize = 3.5, lineWidth = 0.75, cameraPosition = (1.,1.,.5) #src
)#src
render_asymptote(exportFolder*"/jacobiGeodesic.asy"; render = 2) #src
#md #
#md # ```julia
#md # asymptote_export_S2_signals("jacobiGeodesic.asy";
#md #     render = asyResolution,
#md #     curves=[geodesicCurve], points = [ [x,y] ],
#md #     colors=Dict(:curves => [black], :points => [TolVibrantOrange]),
#md #     dotSize = 3.5, lineWidth = 0.75, cameraPosition = (1.,1.,.5)
#md # )
#md # render_asymptote("jacobiGeodesic.asy"; render = 2)
#md # ```
#md #
#md # ![A geodesic connecting two points on the equator](../assets/images/tutorials/jacobiGeodesic.png)
#
# where $x$ is on the left. Then this tutorial solves the following task:
#
# Given a direction $X_p∈ T_x\mathcal M$, for example
X = [0.,0.4,0.5]
# we move the start point $x$ into, how does any point on the geodesic move?
#
# Or mathematically: Compute $D_p g(t; p,q)$ for some fixed $t∈[0,1]$
# and a given direction $X_p$.
# Of course two cases are quite easy: For $t=0$ we are in $x$ and how $x$ “moves”
# is already known, so $D_x g(0;p,q) = X$. On the other side, for $t=1$,
# $g(1; p,q) = q$ which is fixed, so $D_p g(1; p,q)$ is the zero tangent vector
# (in $T_q\mathcal M$).
#
# For all other cases we employ a [`jacobi_field`](@ref), which is a (tangent)
# vector field along the [shortest geodesic](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.shortest_geodesic-Tuple{Manifold,Any,Any}) given as follows: The _geodesic variation_
# $\Gamma_{g,X}(s,t)$ is defined for some $\varepsilon > 0$ as
# ````math
# \Gamma_{g,X}(s,t):=\exp{\gamma_{p,X}(s)}[t\log_{g(s;p,X)}p],\qquad s∈(-\varepsilon,\varepsilon),\ t∈[0,1].
# ````
# Intuitively we make a small step $s$ into direction $\xi$ using the geodesic
# $g(\cdot; p,X)$ and from $r=g(s; p,X)$ we follow (in $t$) the geodesic
# $g(\cdot; r,q)$. The corresponding Jacobi field~\(J_{g,X}\)
# along~\(g(\cdot; p,q)\) is given
#
# ````math
# J_{g,X}(t):=\frac{D}{\partial s}\Gamma_{g,X}(s,t)\Bigl\rvert_{s=0}$
# ````
#
# which is an ODE and we know the boundary conditions $J_{g,X}(0)=X$ and
# $J_{g,X}(t) = 0$. In symmetric spaces we can compute the solution, since the
# system of ODEs decouples, see for example [do Carmo](#doCarmo1992),
# Chapter 4.2. Within `Manopt.jl` this is implemented as
# [`jacobi_field`](@ref)`(M,p,q,t,X[,β])`, where the optional parameter (function)
# `β` specifies, which Jacobi field we want to evaluate and the one used here is
# the default.
#
# We can hence evaluate that on the points on the geodesic at
T = [0:0.1:1.0...]
nothing #hide
# namely
r = shortest_geodesic(M,p,q,T)
nothing #hide
# the geodesic moves as
W = jacobi_field.(Ref(M), Ref(p), Ref(q), T, Ref(X) )
# which can also be called using [`differential_geodesic_startpoint`](@ref).
# We can add to the image above by creating extended tangent vectors
# the include their base points
V = [ Tuple([a,b]) for (a,b) ∈ zip(r,W) ]
# and add that as one further set to the Asymptote export.
asymptote_export_S2_signals(exportFolder*"/jacobiGeodesicdifferential_geodesic_startpoint.asy"; #src
    curves=[geodesicCurve], points = [ [p,q], r], tVectors = [V], #src
    colors=Dict( #src
        :curves => [black], #src
        :points => [TolVibrantOrange,TolVibrantCyan], #src
        :tvectors => [TolVibrantCyan] #src
    ), #src
    dotSizes = [3.5,2.], lineWidth = 0.75, cameraPosition = (1.,1.,.5) #src
) #src
render_asymptote(exportFolder*"/jacobiGeodesicdifferential_geodesic_startpoint.asy"; render = 2) #src
#md #
#md # ```julia
#md # asymptote_export_S2_signals("jacobiGeodesicdifferential_geodesic_startpoint.asy";
#md #     render = asyResolution,
#md #     curves=[geodesicCurve], points = [ [x,y], Z], tVectors = [Vx],
#md #     colors=Dict(
#md #         :curves => [black],
#md #         :points => [TolVibrantOrange,TolVibrantCyan],
#md #         :tvectors => [TolVibrantCyan]
#md #     ),
#md #     dotSizes = [3.5,2.], lineWidth = 0.75, cameraPosition = (1.,1.,.5)
#md # )
#md # render_asymptote("jacobiGeodesicdifferential_geodesic_startpoint.asy"; render = 2)
#md # ```
#
#md # ![A Jacobi field for $D_xg(t,x,y)[\eta]$](../assets/images/tutorials/jacobiGeodesicdifferential_geodesic_startpoint.png)
#
# If we further move the end point, too, we can derive that Differential in direction
Xq = [0.2,0.,-0.5]
W2 = differential_geodesic_endpoint.(Ref(M),Ref(p),Ref(q),T,Ref(Xq))
V2 = [ Tuple([a,b]) for (a,b) ∈ zip(r,W2) ]
# and we can combine both keeping the base point
V3 = [ Tuple([a,b]) for (a,b) ∈ zip(r,W2+W) ]
asymptote_export_S2_signals(exportFolder*"/jacobiGeodesicResult.asy"; #src
    curves=[geodesicCurve], points = [ [p,q], r], tVectors = [V,V2,V3], #src
    colors=Dict( #src
        :curves => [black], #src
        :points => [TolVibrantOrange,TolVibrantCyan], #src
        :tvectors => [TolVibrantCyan,TolVibrantCyan,TolVibrantTeal] #src
    ), #src
    dotSizes = [3.5,2.], lineWidth = 0.75, cameraPosition = (1.,1.,0.) #src
) #src
render_asymptote(exportFolder*"/jacobiGeodesicResult.asy"; render = 2) #src
#md # ```julia
#md # asymptote_export_S2_signals("jacobiGeodesicResult.asy";
#md #    render = asyResolution,
#md #    curves=[geodesicCurve], points = [ [x,y], Z], tVectors = [Vx,Vy,Vb],
#md #    colors=Dict(
#md #        :curves => [black],
#md #        :points => [TolVibrantOrange,TolVibrantCyan],
#md #        :tvectors => [TolVibrantCyan,TolVibrantCyan,TolVibrantTeal]
#md #   ),
#md #   dotSizes = [3.5,2.], lineWidth = 0.75, cameraPosition = (1.,1.,0.)
#md # )
#md # render_asymptote("jacobiGeodesicResult.asy"; render = 2)
#md # ```
#md #
#md # ![A Jacobi field for the effect of two differentials (blue) in sum (teal)](../assets/images/tutorials/jacobiGeodesicResult.png)
#
# ## Literature
#
# ```@raw html
# <ul><li id="doCarmo1992">[<a>doCarmo1992</a>] do Carmo, M. P.:
#    <emph>Riemannian Geometry</emph>, Mathematics: Theory & Applications,
#    Birkhäuser Basel, 1992, ISBN: 0-8176-3490-8</li>
# <li id="BergmannGousenbourger2018">[<a>BergmannGousenbourger2018</a>]
#   Bergmann, R.; Gousenbourger, P.-Y.: <emph>A variational model for data
#   fitting on manifolds by minimizing the acceleration of a Bézier curve</emph>,
#   Frontiers in Applied Mathematics and Statistics, 2018.
#   doi: <a href="https://dx.doi.org/10.3389/fams.2018.00059">10.3389/fams.2018.00059</a></li>
# </ul>
# ```


text:
using Test
using ResourcePools

@testset "PooledAbstractArray" begin
    a = [1 2; 3 4]

    is_disposed = false
    r = PooledAbstractArray(a, (r) -> is_disposed = true)
    @test ref_count(r) == 1
    @test resource(r) === a
    @test r == a

    retain!(r)
    @test ref_count(r) == 2

    release!(r)
    @test ref_count(r) == 1

    @assert !is_disposed "is_disposed is true"
    release!(r)
    @test ref_count(r) == 0
    @test is_disposed

    @test_throws ArgumentError release!(r)
    @test_throws ArgumentError retain!(r)
end


text:
using RedPitayaDAQServer
using PyPlot

rp = RedPitaya("rp-f04972.local")
#rp = RedPitaya("rp-f07083.local")
#rp = RedPitaya("192.168.178.46")

dec = 8
modulus = 4800
base_frequency = 125000000
periods_per_step = 10
samples_per_period = div(modulus, dec)#*periods_per_step
periods_per_frame = 100 # about 0.5 s frame length
frame_period = dec*samples_per_period*periods_per_frame / base_frequency
slow_dac_periods_per_frame = div(periods_per_frame, periods_per_step)

@info samples_per_period, periods_per_frame, slow_dac_periods_per_frame

decimation(rp, dec)
samplesPerPeriod(rp, samples_per_period)
periodsPerFrame(rp, periods_per_frame)
passPDMToFastDAC(rp, true)
slowDACPeriodsPerFrame(rp, slow_dac_periods_per_frame)

numSlowDACChan(rp, 2)
lutA = collect(range(0,0.3,length=slow_dac_periods_per_frame))
lutB = collect(ones(slow_dac_periods_per_frame))

lutEnableDACA = ones(Bool, slow_dac_periods_per_frame)
lutEnableDACA[1:2:end] .= false
lutEnableDACB = ones(Bool, slow_dac_periods_per_frame)
lutEnableDACB[2:2:end] .= false

modeDAC(rp, "STANDARD")
frequencyDAC(rp,1,1, base_frequency / modulus)

freq = frequencyDAC(rp,1,1)
println(" frequency = $(freq)")
jumpSharpnessDAC(rp, 1, 0.01) # controls the sharpness of the jump for the square
signalTypeDAC(rp, 1 , "SQUARE")
phaseDAC(rp, 1, 1, 0.0 ) # Phase has to be given in between 0 and 1
amplitudeDACNext(rp, 1, 1, 0.1)

ramWriterMode(rp, "TRIGGERED")
triggerMode(rp, "INTERNAL")
masterTrigger(rp, false)
startADC(rp)
masterTrigger(rp, true)

sleep(0.1)

enableDACLUT(rp, collect( cat(lutEnableDACA,lutEnableDACB,dims=2)' ) )

lut = cat(lutA,lutB*0.1,dims=2)'
setSlowDACLUT(rp, collect(lut))

currFr = enableSlowDAC(rp, true, 1, frame_period, 0.5)
uCurrentPeriod = readData(rp, currFr, 1)

sleep(0.5)

lut = cat(lutA,lutB*0.4,dims=2)'
setSlowDACLUT(rp, collect(lut))

currFr = enableSlowDAC(rp, true, 1, frame_period, 0.5)
uCurrentPeriod2 = readData(rp, currFr, 1)

figure(1)
clf()
subplot(2,1,1)
plot(vec(uCurrentPeriod[:,1,:,:]))
plot(vec(uCurrentPeriod[:,2,:,:]))
legend(("Rx1", "Rx2"))
subplot(2,1,2)
plot(vec(uCurrentPeriod2[:,1,:,:]))
plot(vec(uCurrentPeriod2[:,2,:,:]))
legend(("Rx1", "Rx2"))

stopADC(rp)

max_stars_repo_name: fredrikekre/RedPitayaDAQServer

max_issues_repo_name: fuerpy/RedPitayaDAQServer

max_forks_repo_name: fuerpy/RedPitayaDAQServer


text:
import MathOptInterface
using Test

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOF = MOI.FileFormats.MOF

const TEST_MOF_FILE = "test.mof.json"

@test sprint(show, MOF.Model()) == "A MathOptFormat Model"

include("nonlinear.jl")

struct UnsupportedSet <: MOI.AbstractSet end
struct UnsupportedFunction <: MOI.AbstractFunction end

function test_model_equality(model_string, variables, constraints; suffix="")
    model = MOF.Model(validate = true)
    MOIU.loadfromstring!(model, model_string)
    MOI.write_to_file(model, TEST_MOF_FILE * suffix)
    model_2 = MOF.Model()
    MOI.read_from_file(model_2, TEST_MOF_FILE * suffix)
    MOIU.test_models_equal(model, model_2, variables, constraints)
    MOF.validate(TEST_MOF_FILE * suffix)
end

@testset "Error handling: read_from_file" begin
    failing_models_dir = joinpath(@__DIR__, "failing_models")

    @testset "Non-empty model" begin
        model = MOF.Model(warn=true)
        MOI.add_variable(model)
        @test !MOI.is_empty(model)
        exception = ErrorException(
            "Cannot read model from file as destination model is not empty.")
        @test_throws exception MOI.read_from_file(
            model, joinpath(@__DIR__, "empty_model.mof.json"))
        options = MOF.get_options(model)
        @test options.warn
        MOI.empty!(model)
        @test MOI.is_empty(model)
        MOI.read_from_file(
            model, joinpath(@__DIR__, "empty_model.mof.json"))
        options2 = MOF.get_options(model)
        @test options2.warn
    end

    @testset "$(filename)" for filename in filter(
        f -> endswith(f, ".mof.json"), readdir(failing_models_dir))
        @test_throws Exception MOI.read_from_file(MOF.Model(),
            joinpath(failing_models_dir, filename))
    end
end
@testset "Names" begin
    @testset "Blank variable name" begin
        model = MOF.Model()
        variable_index = MOI.add_variable(model)
        @test_throws Exception MOF.moi_to_object(variable_index, model)
        MOI.FileFormats.create_unique_names(model, warn=true)
        @test MOF.moi_to_object(variable_index, model) ==
            MOF.OrderedObject("name" => "x1")
    end
    @testset "Duplicate variable name" begin
        model = MOF.Model()
        x = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x")
        y = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), y, "x")
        @test MOF.moi_to_object(x, model) == MOF.OrderedObject("name" => "x")
        @test MOF.moi_to_object(y, model) == MOF.OrderedObject("name" => "x")
        MOI.FileFormats.create_unique_names(model, warn=true)
        @test MOF.moi_to_object(x, model) == MOF.OrderedObject("name" => "x")
        @test MOF.moi_to_object(y, model) == MOF.OrderedObject("name" => "x_1")
    end
    @testset "Blank constraint name" begin
        model = MOF.Model()
        x = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x")
        c = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.ZeroOne())
        name_map = Dict(x => "x")
        MOI.FileFormats.create_unique_names(model, warn=true)
        @test MOF.moi_to_object(c, model, name_map)["name"] == "c1"
    end
    @testset "Duplicate constraint name" begin
        model = MOF.Model()
        x = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x")
        c1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.LessThan(1.0))
        c2 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
        MOI.set(model, MOI.ConstraintName(), c1, "c")
        MOI.set(model, MOI.ConstraintName(), c2, "c")
        name_map = Dict(x => "x")
        @test MOF.moi_to_object(c1, model, name_map)["name"] == "c"
        @test MOF.moi_to_object(c2, model, name_map)["name"] == "c"
        MOI.FileFormats.create_unique_names(model, warn=true)
        @test MOF.moi_to_object(c1, model, name_map)["name"] == "c_1"
        @test MOF.moi_to_object(c2, model, name_map)["name"] == "c"
    end
end
@testset "round trips" begin
    @testset "Empty model" begin
        model = MOF.Model(validate = true)
        MOI.write_to_file(model, TEST_MOF_FILE)
        model_2 = MOF.Model(validate = true)
        MOI.read_from_file(model_2, TEST_MOF_FILE)
        MOIU.test_models_equal(model, model_2, String[], String[])
    end
    @testset "FEASIBILITY_SENSE" begin
        model = MOF.Model(validate = true)
        x = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x")
        MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
        MOI.write_to_file(model, TEST_MOF_FILE)
        model_2 = MOF.Model(validate = true)
        MOI.read_from_file(model_2, TEST_MOF_FILE)
        MOIU.test_models_equal(model, model_2, ["x"], String[])
    end
    @testset "Empty function term" begin
        model = MOF.Model(validate = true)
        x = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x")
        c = MOI.add_constraint(model,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0),
            MOI.GreaterThan(1.0)
        )
        MOI.set(model, MOI.ConstraintName(), c, "c")
        MOI.write_to_file(model, TEST_MOF_FILE)
        model_2 = MOF.Model(validate = true)
        MOI.read_from_file(model_2, TEST_MOF_FILE)
        MOIU.test_models_equal(model, model_2, ["x"], ["c"])
    end
    @testset "min objective" begin
        test_model_equality("""
            variables: x
            minobjective: x
        """, ["x"], String[])
    end
    @testset "max objective" begin
        test_model_equality("""
            variables: x
            maxobjective: x
        """, ["x"], String[], suffix=".gz")
    end
    @testset "min scalaraffine" begin
        test_model_equality("""
            variables: x
            minobjective: 1.2x + 0.5
        """, ["x"], String[])
    end
    @testset "max scalaraffine" begin
        test_model_equality("""
            variables: x
            maxobjective: 1.2x + 0.5
        """, ["x"], String[], suffix=".gz")
    end
    @testset "singlevariable-in-lower" begin
        test_model_equality("""
            variables: x
            minobjective: 1.2x + 0.5
            c1: x >= 1.0
        """, ["x"], ["c1"])
    end
    @testset "singlevariable-in-upper" begin
        test_model_equality("""
            variables: x
            maxobjective: 1.2x + 0.5
            c1: x <= 1.0
        """, ["x"], ["c1"], suffix=".gz")
    end
    @testset "singlevariable-in-interval" begin
        test_model_equality("""
            variables: x
            minobjective: 1.2x + 0.5
            c1: x in Interval(1.0, 2.0)
        """, ["x"], ["c1"])
    end
    @testset "singlevariable-in-equalto" begin
        test_model_equality("""
            variables: x
            minobjective: 1.2x + 0.5
            c1: x == 1.0
        """, ["x"], ["c1"])
    end
    @testset "singlevariable-in-zeroone" begin
        test_model_equality("""
            variables: x
            minobjective: 1.2x + 0.5
            c1: x in ZeroOne()
        """, ["x"], ["c1"])
    end
    @testset "singlevariable-in-integer" begin
        test_model_equality("""
            variables: x
            minobjective: 1.2x + 0.5
            c1: x in Integer()
        """, ["x"], ["c1"])
    end
    @testset "singlevariable-in-Semicontinuous" begin
        test_model_equality("""
            variables: x
            minobjective: 1.2x + 0.5
            c1: x in Semicontinuous(1.0, 2.0)
        """, ["x"], ["c1"])
    end
    @testset "singlevariable-in-Semiinteger" begin
        test_model_equality("""
            variables: x
            minobjective: 1.2x + 0.5
            c1: x in Semiinteger(1.0, 2.0)
        """, ["x"], ["c1"])
    end
    @testset "scalarquadratic-objective" begin
        test_model_equality("""
            variables: x
            minobjective: 1.0*x*x + -2.0x + 1.0
        """, ["x"], String[])
    end
    @testset "SOS1" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in SOS1([1.0, 2.0, 3.0])
        """, ["x", "y", "z"], ["c1"])
    end
    @testset "SOS2" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in SOS2([1.0, 2.0, 3.0])
        """, ["x", "y", "z"], ["c1"])
    end
    @testset "Reals" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in Reals(3)
        """, ["x", "y", "z"], ["c1"])
    end
    @testset "Zeros" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in Zeros(3)
        """, ["x", "y", "z"], ["c1"])
    end
    @testset "Nonnegatives" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in Nonnegatives(3)
        """, ["x", "y", "z"], ["c1"])
    end
    @testset "Nonpositives" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in Nonpositives(3)
        """, ["x", "y", "z"], ["c1"])
    end
    @testset "PowerCone" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in PowerCone(2.0)
        """, ["x", "y", "z"], ["c1"])
    end
    @testset "DualPowerCone" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in DualPowerCone(0.5)
        """, ["x", "y", "z"], ["c1"])
    end
    @testset "GeometricMeanCone" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in GeometricMeanCone(3)
        """, ["x", "y", "z"], ["c1"])
    end
    @testset "vectoraffine-in-zeros" begin
        test_model_equality("""
            variables: x, y
            minobjective: x
            c1: [1.0x + -3.0, 2.0y + -4.0] in Zeros(2)
        """, ["x", "y"], ["c1"])
    end
    @testset "vectorquadratic-in-nonnegatives" begin
        test_model_equality("""
            variables: x, y
            minobjective: x
            c1: [1.0*x*x + -2.0x + 1.0, 2.0y + -4.0] in Nonnegatives(2)
        """, ["x", "y"], ["c1"])
    end
    @testset "ExponentialCone" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in ExponentialCone()
        """, ["x", "y", "z"], ["c1"])
    end
    @testset "DualExponentialCone" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in DualExponentialCone()
        """, ["x", "y", "z"], ["c1"])
    end
    @testset "SecondOrderCone" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in SecondOrderCone(3)
        """, ["x", "y", "z"], ["c1"])
    end
    @testset "RotatedSecondOrderCone" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in RotatedSecondOrderCone(3)
        """, ["x", "y", "z"], ["c1"])
    end
    @testset "PositiveSemidefiniteConeTriangle" begin
        test_model_equality("""
            variables: x1, x2, x3
            minobjective: x1
            c1: [x1, x2, x3] in PositiveSemidefiniteConeTriangle(2)
        """, ["x1", "x2", "x3"], ["c1"])
    end
    @testset "PositiveSemidefiniteConeSquare" begin
        test_model_equality("""
            variables: x1, x2, x3, x4
            minobjective: x1
            c1: [x1, x2, x3, x4] in PositiveSemidefiniteConeSquare(2)
        """, ["x1", "x2", "x3", "x4"], ["c1"])
    end
    @testset "LogDetConeTriangle" begin
        test_model_equality("""
            variables: t, u, x1, x2, x3
            minobjective: x1
            c1: [t, u, x1, x2, x3] in LogDetConeTriangle(2)
        """, ["t", "u", "x1", "x2", "x3"], ["c1"])
    end
    @testset "LogDetConeSquare" begin
        test_model_equality("""
            variables: t, u, x1, x2, x3, x4
            minobjective: x1
            c1: [t, u, x1, x2, x3, x4] in LogDetConeSquare(2)
        """, ["t", "u", "x1", "x2", "x3", "x4"], ["c1"])
    end
    @testset "RootDetConeTriangle" begin
        test_model_equality("""
            variables: t, x1, x2, x3
            minobjective: x1
            c1: [t, x1, x2, x3] in RootDetConeTriangle(2)
        """, ["t", "x1", "x2", "x3"], ["c1"])
    end
    @testset "RootDetConeSquare" begin
        test_model_equality("""
            variables: t, x1, x2, x3, x4
            minobjective: x1
            c1: [t, x1, x2, x3, x4] in RootDetConeSquare(2)
        """, ["t", "x1", "x2", "x3", "x4"], ["c1"])
    end
    @testset "IndicatorSet" begin
        test_model_equality("""
            variables: x, y
            minobjective: x
            c1: [x, y] in IndicatorSet{ACTIVATE_ON_ONE}(GreaterThan(1.0))
            c2: x >= 0.0
        """, ["x", "y"], ["c1", "c2"])

        test_model_equality("""
            variables: x, y
            minobjective: x
            c1: [x, y] in IndicatorSet{ACTIVATE_ON_ZERO}(GreaterThan(1.0))
            c2: x >= 0.0
        """, ["x", "y"], ["c1", "c2"])
    end
    @testset "NormOneCone" begin
        test_model_equality("""
            variables: x, y
            minobjective: x
            c1: [x, y] in NormOneCone(2)
            c2: x >= 0.0
        """, ["x", "y"], ["c1", "c2"])
    end
    @testset "NormInfinityCone" begin
        test_model_equality("""
            variables: x, y
            minobjective: x
            c1: [x, y] in NormInfinityCone(2)
            c2: x >= 0.0
        """, ["x", "y"], ["c1", "c2"])
    end
    @testset "RelativeEntropyCone" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in RelativeEntropyCone(3)
            c2: x >= 0.0
        """, ["x", "y", "z"], ["c1", "c2"])
    end
    @testset "NormSpectralCone" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in NormSpectralCone(1, 2)
        """, ["x", "y", "z"], ["c1"])
    end
    @testset "NormNuclearCone" begin
        test_model_equality("""
            variables: x, y, z
            minobjective: x
            c1: [x, y, z] in NormNuclearCone(1, 2)
        """, ["x", "y", "z"], ["c1"])
    end
    # Clean up
    sleep(1.0)  # allow time for unlink to happen
    rm(TEST_MOF_FILE, force=true)
    rm(TEST_MOF_FILE * ".gz", force=true)
end


text:
module TestWithTimeout

using Terminators
using Test
using ..Utils: exec

function test_no_timeout()
    ok = Ref(false)
    Terminators.withtimeout(10) do
        ok[] = true
    end
    @test ok[]
end

function test_timeout()
    proc = Terminators.withtimeout(300) do
        exec("""
        using Terminators
        Terminators.withtimeout(0.1) do
            while true
                $(VERSION ≥ v"1.4" ? "GC.safepoint()" : "yield()")
            end
        end
        """)
    end
    @test !success(proc)
    @test occursin("Timeout", proc.stderr)
    @test occursin("terminating process", proc.stderr)
end

end  # module


text:
ALL_TASK_OFFSETS[("linux-x86_64", v"1.6.1")] = Dict(
    :END => 312,
    :_isexception => 58,
    :_state => 56,
    :bufsz => 288,
    :copy_stack => 296,
    :ctx => 80,
    :donenotify => 24,
    :eh => 72,
    :excstack => 64,
    :gcstack => 304,
    :logstate => 40,
    :next => 0,
    :prio => 62,
    :queue => 8,
    :result => 32,
    :sizeof_ctx => 200,
    :start => 48,
    :sticky => 57,
    :stkbuf => 280,
    :tid => 60,
    :tls => 16,
    :tls_base_context => 6656,
)


text:
export Toeplitz, levinson, levinson!

struct Toeplitz{T<:AbstractFloat,S} <: AbstractMatrix{T}
    n::Int
    l::Int
    u::Int
    t::S # vector of length l + u + 1 < 2n
    f::Vector{T}
    b::Vector{T}
    function Toeplitz(n, l, u, t::AbstractVector)
        T = float(eltype(t))
        S = Vector{T}
        l <= u + 1 && -l + u + 1 <= 2n || throw(ArgumentError("l + u + 1 < 2n"))
        length(t) >= -l + u + 1 || throw(ArgumentError("vector insufficient size"))
        new{T,S}(n, l, u, t, Vector{T}(undef,0), Vector{T}(undef,0))
    end
end

function Base.getindex(A::Toeplitz{T}, i, j) where T
    A.l <= j - i <= A.u ? getindex(A.t, j-i-A.l+1) : zero(T)
end
Base.size(A::Toeplitz) = (A.n, A.n)

function levinson!(A::Toeplitz{T}) where T
    n = A.n
    if length(A.f) == 0
        resize!(A.f, n*(n+1)÷2)
        resize!(A.b, n*(n+1)÷2)
        fill!(A.f, zero(T))
        fill!(A.b, zero(T))

        toeplitz_forward(A, 1)[1] = toeplitz_backward(A, 1)[1] = inv(float(A[1,1]))
        for i = 2:n
            fi = toeplitz_forward(A, i)
            fi1 = toeplitz_forward(A, i-1)
            bi = toeplitz_backward(A, i)
            bi1 = toeplitz_backward(A, i-1)
            efi = sum(A[i,j] * fi1[j] for j = 1:i-1)
            ebi = sum(A[1,j+1] * bi1[j] for j = 1:i-1)
            ebf = inv(one(T) - efi*ebi)
            fi[1] = fi1[1] * ebf
            fi[2:i-1] = fi1[2:i-1] * ebf - bi1[1:i-2] * efi * ebf
            fi[i] = -bi1[i-1] * efi * ebf
            bi[1] = -fi1[1] * ebi * ebf
            bi[2:i-1] = bi1[1:i-2] * ebf - fi1[2:i-1] * ebi * ebf
            bi[i] = bi1[i-1] * ebf
        end
    end
    A
end

toeplitz_forward(A::Toeplitz, i) = view(A.f, i*(i-1)÷2+1:i*(i+1)÷2)
toeplitz_backward(A::Toeplitz, i) = view(A.b, i*(i-1)÷2+1:i*(i+1)÷2)

function levinson(A::Toeplitz{T}, B::AbstractArray{T,N}) where {T,N}
    N <= 2 || throw(ArgumentError("rhs must be vector or matrix"))
    n = size(B,1)
    n == A.n || throw(DimensionMismatch("rhs not correct size"))
    r = N == 2 ? size(B, 2) : 1
    x = N == 1 ? Vector{T}(undef, n) : Matrix{T}(undef, n, n)
    for k = 1:r
        x[1,k] = A[1,1] \ B[1,k]
        for i = 2:n
            bi = toeplitz_backward(A, i)
            exi = sum(A[i,j] * x[j,k] for j = 1:i-1)
            x[1:i,k] += bi * (B[i,k] - exi)
        end
    end
    x
end


text:
module OpenTelemetryUber

using Reexport

@reexport using OpenTelemetryAPI
@reexport using OpenTelemetrySDK
@reexport using OpenTelemetryExporterOtlpProtoGrpc
@reexport using OpenTelemetryExporterPrometheus
@reexport using OpenTelemetryProto

# instrumentations
@reexport using OpenTelemetryInstrumentationBase
@reexport using OpenTelemetryInstrumentationDownloads
@reexport using OpenTelemetryInstrumentationDistributed
@reexport using OpenTelemetryInstrumentationHTTP
@reexport using OpenTelemetryInstrumentationGenie

function init(;
    meter_provider = global_meter_provider(),
    tracer_provider = global_tracer_provider(),
)
    global_meter_provider!(meter_provider)
    global_tracer_provider!(tracer_provider)
    OpenTelemetryInstrumentationBase.init(;
        meter_provider = meter_provider,
        tracer_provider = tracer_provider,
    )

    OpenTelemetryInstrumentationDownloads.init(;
        meter_provider = meter_provider,
        tracer_provider = tracer_provider,
    )
end

end # module


text:
"""
    constraint_equations(f::ODEFunction)

Returns the constraint equations from an ODEFunction f used in DifferentialEquations.jl.

f must be in Mass Matrix form meaning: M ẋ = f(x).
"""
function constraint_equations(f::ODEFunction)
    M = f.mass_matrix
    if M == true*I
        error("There are no constraints in the system!")
    end
    len_M = size(M, 1)
    g_idx = findall([M[i,i] for i in 1:len_M] .== 0)
    g(x) = (dx = similar(x);
            f(dx, x, nothing, 0.0);
            dx[g_idx])
end

"""
    AmbientForcing(f::ODEFunction, z, tau, Frand)
Draws random initial condition form the constraint manifold M.

This consists of the following steps:
    1. Draw a random vector Frand from the ambient space
    2. Projetion back onto the tangetial space TxM by using PN

# Arguments
    f: An ODEFunction in mass matrix form
    initial_con: A point which lays on the manifold M, could e.g. be a fixed point
    dist_args: Arguments of the distribution function dist
    dist: Distribution function e.g. Uniform, Normal....
    tau: Integration time
"""
function ambient_forcing(f::ODEFunction, z, tau, Frand)
    g = constraint_equations(f)
    prob = ODEProblem(ambient_forcing_ODE, z, (0.0, tau), (g, Frand))
    sol = solve(prob, Tsit5())
    return sol[end]
end

"""
    Proj_N(A)

Calculates the orthogonal projection matrix on a subspace N.

A is a matrix containg the basis vectors of N as columns.
The matrix (A^T*A)^-1 recovers the norm.
"""
Proj_N(A) = A * inv(transpose(A) * A) * transpose(A)

"""
    ambient_forcing_ODE(u, p, t)
The manifold preserving ODE version of a vector h.

ż = Pn(ker(Jg)) * h
# Arguments
    `p[1]`: g the function defining the manifold
    `p[2]` h constant vector from ambient space
    `u0`: Inital condition
    `t`: Integration time
"""
function ambient_forcing_ODE(u, p, t)
    g, h = p
    Jacg = jacobian(g, u)
    N = nullspace(Jacg)
    du = Proj_N(N) * h # Projection back onto the Manifold
end


text:
# This file includes functions for managing opentdb.com tokens. A token is a unique string
# which the server uses to identify an individual user. The opentdb.com API will use your
# unique token to ensure that you don't get a trivia question you've already seen.
#
# By default, this package will request a token for you and store it in `data/token.txt`.
# Any time you call `Token()` again, that saved token will be reused. You can get a brand
# new token with `request_token()` or you can simply delete the `data/token.txt` file if
# you want to completely start over.

default_token_file() = joinpath(@__DIR__, "..", "data", "token.txt")

"""
Basic struct wrapping a single string holding an opentdb.com API token
"""
struct Token
    value::String
end

"""
Construct a Token. By default, this will use an existing saved token in
`data/token.txt`. If no token exists, one will be requested and saved.
"""
Token() = get_or_request_token()

"""
Request a brand new token.
"""
function request_token()
    json = opentdb_request("https://opentdb.com/api_token.php?command=request")
    Token(json["token"])
end

"""
Save the given token using the givne file name.
"""
function save(filename::AbstractString, token::Token)
    open(filename, write=true) do file
        write(file, token.value)
    end
end

"""
Construct a Token. By default, this will use an existing saved token in
`filename`. If no token exists, one will be requested and saved.
"""
function get_or_request_token(filename=default_token_file())
    if isfile(filename)
        token = Token(open(r -> strip(read(r, String)), filename))
    else
        token = request_token()
        mkpath(dirname(filename))
        save(filename, token)
    end
    token
end

"""
Tell the opentdb.com API to reset the given token. This deletes the history of that token
on the server, so when using this token again you may see questions that you have already
seen.
"""
function reset!(t::Token)
    opentdb_request("https://opentdb.com/api_token.php?command=reset&token=" * t.value)
end


text:
## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 1 on the triangle *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TRIANGLE_G1N1 = (
    (SVector(0.333333333333333, 0.333333333333333), 0.500000000000000),)
## 0 negative weights, 0 points outside of the triangle,  total sum of the
## weights is 0.5

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 2 on the triangle *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TRIANGLE_G2N3 = (
    (SVector(0.166666666666667, 0.166666666666667), 0.166666666666667),
    (SVector(0.166666666666667, 0.666666666666667), 0.166666666666667),
    (SVector(0.666666666666667, 0.166666666666667), 0.166666666666667))
## 0 negative weights, 0 points outside of the triangle,  total sum of the
## weights is 0.5

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 3 on the triangle *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TRIANGLE_G3N4 = (
    (SVector(0.333333333333333, 0.333333333333333), -0.281250000000000),
    (SVector(0.200000000000000, 0.200000000000000), 0.260416666666667),
    (SVector(0.200000000000000, 0.600000000000000), 0.260416666666667),
    (SVector(0.600000000000000, 0.200000000000000), 0.260416666666667))
## 1 negative weight, 0 points outside of the triangle,  total sum of the
## weights is 0.5

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 4 on the triangle *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TRIANGLE_G4N6 = (
    (SVector(0.445948490915965, 0.445948490915965), 0.111690794839005),
    (SVector(0.445948490915965, 0.108103018168070), 0.111690794839005),
    (SVector(0.108103018168070, 0.445948490915965), 0.111690794839005),
    (SVector(0.091576213509771, 0.091576213509771), 0.054975871827661),
    (SVector(0.091576213509771, 0.816847572980459), 0.054975871827661),
    (SVector(0.816847572980459, 0.091576213509771), 0.054975871827661))
## 0 negative weights, 0 points outside of the triangle,  total sum of the
## weights is 0.5

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 5 on the triangle *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TRIANGLE_G5N7 = (
    (SVector(0.333333333333333, 0.333333333333333), 0.112500000000000),
    (SVector(0.470142064105115, 0.470142064105115), 0.066197076394253),
    (SVector(0.470142064105115, 0.059715871789770), 0.066197076394253),
    (SVector(0.059715871789770, 0.470142064105115), 0.066197076394253),
    (SVector(0.101286507323456, 0.101286507323456), 0.062969590272414),
    (SVector(0.101286507323456, 0.797426985353087), 0.062969590272414),
    (SVector(0.797426985353087, 0.101286507323456), 0.062969590272414))
## 0 negative weights, 0 points outside of the triangle,  total sum of the
## weights is 0.5

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 6 on the triangle *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TRIANGLE_G6N12 = (
    (SVector(0.249286745170910, 0.249286745170910), 0.058393137863189),
    (SVector(0.249286745170910, 0.501426509658179), 0.058393137863189),
    (SVector(0.501426509658179, 0.249286745170910), 0.058393137863189),
    (SVector(0.063089014491502, 0.063089014491502), 0.025422453185104),
    (SVector(0.063089014491502, 0.873821971016996), 0.025422453185104),
    (SVector(0.873821971016996, 0.063089014491502), 0.025422453185104),
    (SVector(0.310352451033785, 0.636502499121399), 0.041425537809187),
    (SVector(0.636502499121399, 0.053145049844816), 0.041425537809187),
    (SVector(0.053145049844816, 0.310352451033785), 0.041425537809187),
    (SVector(0.310352451033785, 0.053145049844816), 0.041425537809187),
    (SVector(0.636502499121399, 0.310352451033785), 0.041425537809187),
    (SVector(0.053145049844816, 0.636502499121399), 0.041425537809187))
## 0 negative weights, 0 points outside of the triangle,  total sum of the
## weights is 0.5

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 7 on the triangle *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TRIANGLE_G7N13 = (
    (SVector(0.333333333333333, 0.333333333333333), -0.074785022233841),
    (SVector(0.260345966079040, 0.260345966079040), 0.087807628716604),
    (SVector(0.260345966079040, 0.479308067841920), 0.087807628716604),
    (SVector(0.479308067841920, 0.260345966079040), 0.087807628716604),
    (SVector(0.065130102902216, 0.065130102902216), 0.026673617804419),
    (SVector(0.065130102902216, 0.869739794195568), 0.026673617804419),
    (SVector(0.869739794195568, 0.065130102902216), 0.026673617804419),
    (SVector(0.312865496004874, 0.638444188569810), 0.038556880445128),
    (SVector(0.638444188569810, 0.048690315425316), 0.038556880445128),
    (SVector(0.048690315425316, 0.312865496004874), 0.038556880445128),
    (SVector(0.312865496004874, 0.048690315425316), 0.038556880445128),
    (SVector(0.638444188569810, 0.312865496004874), 0.038556880445128),
    (SVector(0.048690315425316, 0.638444188569810), 0.038556880445128))
## 1 negative weight, 0 points outside of the triangle,  total sum of the
## weights is 0.5

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 8 on the triangle *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TRIANGLE_G8N16 = (
    (SVector(0.333333333333333, 0.333333333333333), 0.072157803838894),
    (SVector(0.459292588292723, 0.459292588292723), 0.047545817133643),
    (SVector(0.459292588292723, 0.081414823414554), 0.047545817133643),
    (SVector(0.081414823414554, 0.459292588292723), 0.047545817133643),
    (SVector(0.170569307751760, 0.170569307751760), 0.051608685267359),
    (SVector(0.170569307751760, 0.658861384496480), 0.051608685267359),
    (SVector(0.658861384496480, 0.170569307751760), 0.051608685267359),
    (SVector(0.050547228317031, 0.050547228317031), 0.016229248811599),
    (SVector(0.050547228317031, 0.898905543365938), 0.016229248811599),
    (SVector(0.898905543365938, 0.050547228317031), 0.016229248811599),
    (SVector(0.263112829634638, 0.728492392955404), 0.013615157087217),
    (SVector(0.728492392955404, 0.008394777409958), 0.013615157087217),
    (SVector(0.008394777409958, 0.263112829634638), 0.013615157087217),
    (SVector(0.263112829634638, 0.008394777409958), 0.013615157087217),
    (SVector(0.728492392955404, 0.263112829634638), 0.013615157087217),
    (SVector(0.008394777409958, 0.728492392955404), 0.013615157087217))
## 0 negative weights, 0 points outside of the triangle,  total sum of the
## weights is 0.5

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 9 on the triangle *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TRIANGLE_G9N19 = (
    (SVector(0.3333333333330, 0.3333333333333), 0.04856789814140),
    (SVector(0.4896825191990, 0.4896825191990), 0.01566735011355),
    (SVector(0.4896825191990, 0.0206349616025), 0.01566735011355),
    (SVector(0.0206349616025, 0.4896825191990), 0.01566735011355),
    (SVector(0.4370895914930, 0.4370895914930), 0.03891377050240),
    (SVector(0.4370895914930, 0.1258208170140), 0.03891377050240),
    (SVector(0.1258208170140, 0.4370895914930), 0.03891377050240),
    (SVector(0.1882035356190, 0.1882035356190), 0.03982386946360),
    (SVector(0.1882035356190, 0.6235929287620), 0.03982386946360),
    (SVector(0.6235929287620, 0.1882035356190), 0.03982386946360),
    (SVector(0.0447295133945, 0.0447295133945), 0.01278883782935),
    (SVector(0.0447295133945, 0.9105409732110), 0.01278883782935),
    (SVector(0.9105409732110, 0.0447295133945), 0.01278883782935),
    (SVector(0.2219629891610, 0.7411985987840), 0.02164176968865),
    (SVector(0.7411985987840, 0.0368384120547), 0.02164176968865),
    (SVector(0.0368384120547, 0.2219629891610), 0.02164176968865),
    (SVector(0.2219629891610, 0.0368384120547), 0.02164176968865),
    (SVector(0.7411985987840, 0.2219629891610), 0.02164176968865),
    (SVector(0.0368384120547, 0.7411985987840), 0.02164176968865))
## 0 negative weights, 0 points outside of the triangle,  total sum of the
## weights is 0.5

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 10 on the triangle *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TRIANGLE_G10N25 = (
    (SVector(0.3333333333330, 0.3333333333330), 0.04540899519140),
    (SVector(0.4855776333840, 0.4855776333840), 0.01836297887825),
    (SVector(0.4855776333840, 0.0288447332327), 0.01836297887825),
    (SVector(0.0288447332327, 0.4855776333840), 0.01836297887825),
    (SVector(0.1094815754850, 0.1094815754850), 0.02266052971775),
    (SVector(0.1094815754850, 0.7810368490300), 0.02266052971775),
    (SVector(0.7810368490300, 0.1094815754850), 0.02266052971775),
    (SVector(0.3079398387640, 0.5503529418210), 0.03637895842270),
    (SVector(0.5503529418210, 0.1417072194150), 0.03637895842270),
    (SVector(0.1417072194150, 0.3079398387640), 0.03637895842270),
    (SVector(0.3079398387640, 0.1417072194150), 0.03637895842270),
    (SVector(0.5503529418210, 0.3079398387640), 0.03637895842270),
    (SVector(0.1417072194150, 0.5503529418210), 0.03637895842270),
    (SVector(0.2466725606400, 0.7283239045970), 0.01416362126555),
    (SVector(0.7283239045970, 0.0250035347627), 0.01416362126555),
    (SVector(0.0250035347627, 0.2466725606400), 0.01416362126555),
    (SVector(0.2466725606400, 0.0250035347627), 0.01416362126555),
    (SVector(0.7283239045970, 0.2466725606400), 0.01416362126555),
    (SVector(0.0250035347627, 0.7283239045970), 0.01416362126555),
    (SVector(0.0668032510122, 0.9236559335870), 0.00471083348185),
    (SVector(0.9236559335870, 0.0095408154003), 0.00471083348185),
    (SVector(0.0095408154003, 0.0668032510122), 0.00471083348185),
    (SVector(0.0668032510122, 0.0095408154003), 0.00471083348185),
    (SVector(0.9236559335870, 0.0668032510122), 0.00471083348185),
    (SVector(0.0095408154003, 0.9236559335870), 0.00471083348185))
## 0 negative weights, 0 points outside of the triangle,  total sum of the
## weights is 0.5

## -----------------------------------------------------------------------------
#*! Quadrature rule for an interpolation of order 12 on the triangle *#
#* 'Higher-order Finite Elements', P.Solin, K.Segeth and I. Dolezel *#

const TRIANGLE_G12N33 = (
    (SVector(0.4882173897740, 0.4882173897740), 0.01286553322025),
    (SVector(0.4882173897740, 0.0235652204524), 0.01286553322025),
    (SVector(0.0235652204524, 0.4882173897740), 0.01286553322025),
    (SVector(0.4397243922940, 0.4397243922940), 0.02184627226900),
    (SVector(0.4397243922940, 0.1205512154110), 0.02184627226900),
    (SVector(0.1205512154110, 0.4397243922940), 0.02184627226900),
    (SVector(0.2712103850120, 0.2712103850120), 0.03142911210895),
    (SVector(0.2712103850120, 0.4575792299760), 0.03142911210895),
    (SVector(0.4575792299760, 0.2712103850120), 0.03142911210895),
    (SVector(0.1275761455420, 0.1275761455420), 0.01739805646535),
    (SVector(0.1275761455420, 0.7448477089170), 0.01739805646535),
    (SVector(0.7448477089170, 0.1275761455420), 0.01739805646535),
    (SVector(0.0213173504532, 0.0213173504532), 0.00308313052578),
    (SVector(0.0213173504532, 0.9573652990940), 0.00308313052578),
    (SVector(0.9573652990940, 0.0213173504532), 0.00308313052578),
    (SVector(0.2757132696860, 0.6089432357800), 0.02018577888320),
    (SVector(0.6089432357800, 0.1153434945350), 0.02018577888320),
    (SVector(0.1153434945350, 0.2757132696860), 0.02018577888320),
    (SVector(0.2757132696860, 0.1153434945350), 0.02018577888320),
    (SVector(0.6089432357800, 0.2757132696860), 0.02018577888320),
    (SVector(0.1153434945350, 0.6089432357800), 0.02018577888320),
    (SVector(0.2813255809900, 0.6958360867880), 0.01117838660115),
    (SVector(0.6958360867880, 0.0228383322223), 0.01117838660115),
    (SVector(0.0228383322223, 0.2813255809900), 0.01117838660115),
    (SVector(0.2813255809900, 0.0228383322223), 0.01117838660115),
    (SVector(0.6958360867880, 0.2813255809900), 0.01117838660115),
    (SVector(0.0228383322223, 0.6958360867880), 0.01117838660115),
    (SVector(0.1162519159080, 0.8580140335440), 0.00865811555435),
    (SVector(0.8580140335440, 0.0257340505483), 0.00865811555435),
    (SVector(0.0257340505483, 0.1162519159080), 0.00865811555435),
    (SVector(0.1162519159080, 0.0257340505483), 0.00865811555435),
    (SVector(0.8580140335440, 0.1162519159080), 0.00865811555435),
    (SVector(0.0257340505483, 0.8580140335440), 0.00865811555435))
## 0 negative weights, 0 points outside of the triangle,  total sum of the
## weights is 0.5

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
"""
 Tests for SIE and SME methods
"""

import Statistics # for mean values of trajectories
import LinearAlgebra # for the normn
using StochasticDiffEq
using Test
using Random
using DiffEqDevTools
#using DiffEqGPU

function prob_func(prob, i, repeat)
    remake(prob,seed=seeds[i])
end

"""
 Test Scalar SDEs (oop)
"""

@info "Scalar oop noise"

# PC exercise 14.2.2
numtraj = Int(7e3)
u₀ = 0.1
f(u,p,t) = p[1]*u
g(u,p,t) = p[2]*u
dts = 1 .//2 .^(6:-1:3)
tspan = (0.0,1.0)
p = [3//2,1//100]

h1(z) = z
#analytical_sol(t) = E(f(X(t))) =

seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

prob = SDEProblem(f,g,u₀,tspan,p)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h1(sol[end]),false),
        prob_func = prob_func
        )

sim = test_convergence(dts,ensemble_prob,SIEA(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=u₀.*exp(1.0*(p[1]))
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.35
println("SIEA:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,SMEA(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=u₀.*exp(1.0*(p[1]))
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.34
println("SMEA:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,SIEB(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=u₀.*exp(1.0*(p[1]))
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.35
println("SIEB:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,SMEB(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=u₀.*exp(1.0*(p[1]))
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.34
println("SMEB:", sim.𝒪est[:weak_final])

"""
 Test Scalar SDEs (iip)
"""

@info "Scalar iip noise"

u₀ = [0.1]
f1!(du,u,p,t) =  (du[1] = p[1]*u[1])
g1!(du,u,p,t) =  (du[1] = p[2]*u[1])
dts = 1 .//2 .^(6:-1:3)
tspan = (0.0,1.0)
p = [3//2,1//100]

prob = SDEProblem(f1!,g1!,u₀,tspan,p)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h1(sol[end][1]),false),
        prob_func = prob_func
        )

numtraj = Int(7e3)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,SIEA(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=u₀.*exp(1.0*(p[1]))
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.35
println("SIEA:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,SMEA(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=u₀.*exp(1.0*(p[1]))
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.34
println("SMEA:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,SIEB(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=u₀.*exp(1.0*(p[1]))
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.35
println("SIEB:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,SMEB(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=u₀.*exp(1.0*(p[1]))
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.35 # order is 2.34
println("SMEB:", sim.𝒪est[:weak_final])

"""
 Test Diagonal noise SDEs (iip), SIAM Journal on Numerical Analysis, 47 (2009), pp. 1713–1738
"""

@info "Diagonal noise"

u₀ = [0.1,0.1]
function f3!(du,u,p,t)
  du[1] = 3//2*u[1]
  du[2] = 3//2*u[2]
end
function g3!(du,u,p,t)
  du[1] = 1//10*u[1]
  du[2] = 1//10*u[2]
end
dts = 1 .//2 .^(5:-1:1)
tspan = (0.0,1.0)

h3(z) = z^2 # == 1//10**exp(3//2*t) if h3(z) = z and  == 1//100**exp(301//100*t) if h3(z) = z^2 )

prob = SDEProblem(f3!,g3!,u₀,tspan)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h3(sol[end][1]),false),
        prob_func = prob_func
        )

numtraj = Int(5e4)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,SIEA(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("SIEA:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,SMEA(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("SMEA:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,SIEB(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("SIEB:", sim.𝒪est[:weak_final])

sim = test_convergence(dts,ensemble_prob,SMEB(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
)
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("SMEB:", sim.𝒪est[:weak_final])


text:
using Statistics
using Oceananigans.Architectures

using OffsetArrays: OffsetArray
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid

import Base: getindex, size, show, minimum, maximum
import Statistics: mean

import Oceananigans.Fields: AbstractField, AbstractDataField, AbstractReducedField, Field, ReducedField, minimum, maximum, mean, location, short_show, KernelComputedField
import Oceananigans.Grids: new_data
import Oceananigans.BoundaryConditions: FieldBoundaryConditions

struct CubedSphereFaces{E, F}
    faces :: F
end

function CubedSphereFaces(faces::F) where F
    E = typeof(faces[1])
    return CubedSphereFaces{E, F}(faces)
end

@inline Base.getindex(f::CubedSphereFaces, i::Int) = @inbounds f.faces[i]

#####
##### Dispatch the world / insane type unions
#####

const CubedSphereData = CubedSphereFaces{<:OffsetArray}

# Some dispatch foo to make a type union for CubedSphereFaceField...
#
# Conformal cubed sphere grid wrapped in ImmersedBoundaryGrid:
const ImmersedConformalCubedSphereFaceGrid = ImmersedBoundaryGrid{FT, TX, TY, TZ, <:ConformalCubedSphereFaceGrid} where {FT, TX, TY, TZ}

# CubedSphereFaceField:
const NonImmersedCubedSphereFaceField = AbstractField{X, Y, Z, A, <:ConformalCubedSphereFaceGrid} where {X, Y, Z, A}
const ImmersedCubedSphereFaceField    = AbstractField{X, Y, Z, A, <:ImmersedConformalCubedSphereFaceGrid} where {X, Y, Z, A}

const CubedSphereFaceField = Union{NonImmersedCubedSphereFaceField{X, Y, Z, A},
                                      ImmersedCubedSphereFaceField{X, Y, Z, A}} where {X, Y, Z, A}

const NonImmersedCubedSphereAbstractReducedFaceField = AbstractReducedField{X, Y, Z, A, D, <:ConformalCubedSphereFaceGrid} where {X, Y, Z, A, D}
const ImmersedCubedSphereAbstractReducedFaceField = AbstractReducedField{X, Y, Z, A, D, <:ImmersedConformalCubedSphereFaceGrid} where {X, Y, Z, A, D}

const CubedSphereAbstractReducedFaceField = Union{NonImmersedCubedSphereAbstractReducedFaceField{X, Y, Z, A},
                                                     ImmersedCubedSphereAbstractReducedFaceField{X, Y, Z, A}} where {X, Y, Z, A}

# CubedSphereField

# Flavors of CubedSphereField
const CubedSphereField                     = Field{X, Y, Z, A, <:CubedSphereData} where {X, Y, Z, A}
const CubedSphereReducedField              = ReducedField{X, Y, Z, A, <:CubedSphereData} where {X, Y, Z, A}
const CubedSphereAbstractField             = AbstractField{X, Y, Z, A, <:ConformalCubedSphereGrid} where {X, Y, Z, A}
const CubedSphereAbstractDataField         = AbstractDataField{X, Y, Z, A, <:ConformalCubedSphereGrid} where {X, Y, Z, A}

const AbstractCubedSphereField{X, Y, Z, A} = Union{    CubedSphereAbstractField{X, Y, Z, A},
                                                   CubedSphereAbstractDataField{X, Y, Z, A},
                                                        CubedSphereReducedField{X, Y, Z, A},
                                                               CubedSphereField{X, Y, Z, A}} where {X, Y, Z, A}

#####
##### new data
#####

function new_data(FT, arch::CPU, grid::ConformalCubedSphereGrid, (X, Y, Z))
    faces = Tuple(new_data(FT, arch, face_grid, (X, Y, Z)) for face_grid in grid.faces)
    return CubedSphereFaces{typeof(faces[1]), typeof(faces)}(faces)
end

function new_data(FT, arch::GPU, grid::ConformalCubedSphereGrid, (X, Y, Z))
    faces = Tuple(new_data(FT, arch, face_grid, (X, Y, Z)) for face_grid in grid.faces)
    return CubedSphereFaces{typeof(faces[1]), typeof(faces)}(faces)
end

#####
##### FieldBoundaryConditions
#####

function FieldBoundaryConditions(grid::ConformalCubedSphereGrid, (X, Y, Z); user_defined_bcs...)

    faces = Tuple(
        inject_cubed_sphere_exchange_boundary_conditions(
            FieldBoundaryConditions(face_grid, (X, Y, Z); user_defined_bcs...),
            face_index,
            grid.face_connectivity
        )
        for (face_index, face_grid) in enumerate(grid.faces)
    )

    return CubedSphereFaces{typeof(faces[1]), typeof(faces)}(faces)
end

#####
##### Utils
#####

function Base.show(io::IO, field::AbstractCubedSphereField)
    LX, LY, LZ = location(field)
    arch = architecture(field)
    A = typeof(arch)
    return print(io, "$(typeof(field).name.wrapper) at ($LX, $LY, $LZ)\n",
          "├── architecture: $A\n",
          "└── grid: $(short_show(field.grid))")
end

@inline function interior(field::AbstractCubedSphereField)
    faces = Tuple(interior(face_field_face) for face_field in faces(field))
    return CubedSphereFaces(faces)
end

Base.size(data::CubedSphereData) = (size(data.faces[1])..., length(data.faces))

@inline function get_face(field::CubedSphereField, face_index)
    X, Y, Z = location(field)

    # Should we define a new lower-level constructor for Field that doesn't call validate_field_data?
    return Field{X, Y, Z}(get_face(field.data, face_index),
                          field.architecture,
                          get_face(field.grid, face_index),
                          get_face(field.boundary_conditions, face_index))
end

@inline function get_face(reduced_field::CubedSphereReducedField, face_index)
    X, Y, Z = location(reduced_field)

    return ReducedField{X, Y, Z}(get_face(reduced_field.data, face_index),
                                 reduced_field.architecture,
                                 get_face(reduced_field.grid, face_index),
                                 reduced_field.dims,
                                 get_face(reduced_field.boundary_conditions, face_index))
end

@inline function get_face(kernel_field::CubedSphereAbstractDataField, face_index)
    X, Y, Z = location(kernel_field)

    return Field{X, Y, Z}(get_face(kernel_field.data, face_index),
                                 kernel_field.architecture,
                                 get_face(kernel_field.grid, face_index),
                                 get_face(kernel_field.boundary_conditions, face_index))
end

faces(field::AbstractCubedSphereField) = Tuple(get_face(field, face_index) for face_index in 1:length(field.data.faces))

minimum(field::AbstractCubedSphereField; dims=:) = minimum(minimum(face_field; dims) for face_field in faces(field))
maximum(field::AbstractCubedSphereField; dims=:) = maximum(maximum(face_field; dims) for face_field in faces(field))
mean(field::AbstractCubedSphereField; dims=:) = mean(mean(face_field; dims) for face_field in faces(field))

minimum(f::Function, field::AbstractCubedSphereField; dims=:) = minimum(minimum(f, face_field; dims) for face_field in faces(field))
maximum(f::Function, field::AbstractCubedSphereField; dims=:) = maximum(maximum(f, face_field; dims) for face_field in faces(field))
mean(f::Function, field::AbstractCubedSphereField; dims=:) = mean(mean(f, face_field; dims) for face_field in faces(field))

λnodes(field::CubedSphereFaceField{LX, LY, LZ}) where {LX, LY, LZ} = λnodes(LX(), LY(), LZ(), field.grid)
φnodes(field::CubedSphereFaceField{LX, LY, LZ}) where {LX, LY, LZ} = φnodes(LX(), LY(), LZ(), field.grid)


text:
@enum RunTask begin
    visualize
    simulate
    test
end

runTask = simulate

"""
    run(path; task::RunTask=simulate)

Load model and run a simulation task of file specified by `path`. This function sets the global
variable `Modia3D.runTask` as defined by `task` and executes `include(path)`. Afterwards
`Modia3D.runTask` is reset to `simulate`. Other supported simulation tasks are `visualize` for
visualization and `test` for simulation with result checks.
"""
function run(path; task::RunTask=simulate)
    global runTask = task
    include(path)
    global runTask = simulate
    return nothing
end