text:
struct CompileInfo
    invcheckon::Ref{Bool}
end
CompileInfo() = CompileInfo(Ref(true))

function compile_body(m::Module, body::AbstractVector, info)
    out = []
    for ex in body
        ex_ = compile_ex(m, ex, info)
        ex_ !== nothing && push!(out, ex_)
    end
    return out
end

deleteindex!(d::AbstractDict, index) = delete!(d, index)

@inline function map_func(x::Symbol)
    if x == :+=
        PlusEq, false
    elseif x == :.+=
        PlusEq, true
    elseif x == :-=
        MinusEq, false
    elseif x == :.-=
        MinusEq, true
    elseif x == :*=
        MulEq, false
    elseif x == :.*=
        MulEq, true
    elseif x == :/=
        DivEq, false
    elseif x == :./=
        DivEq, true
    elseif x == :⊻=
        XorEq, false
    elseif x == :.⊻=
        XorEq, true
    else
        error("`$x` can not be mapped to a reversible function.")
    end
end

# e.g. map `x += sin(z)` => `PlusEq(sin)(x, z)`.
function to_standard_format(ex::Expr)
    head::Symbol = ex.head
    F, isbcast = map_func(ex.head)
    a, b = ex.args
    if !isbcast
        @smatch b begin
            :($f($(args...); $(kwargs...))) => :($F($f)($a, $(args...); $(kwargs...)))
            :($f($(args...))) => :($F($f)($a, $(args...)))
            :($x || $y) => :($F($logical_or)($a, $x, $y))
            :($x && $y) => :($F($logical_and)($a, $x, $y))
            _ => :($F(identity)($a, $b))
        end
    else
        @smatch b begin
            :($f.($(args...); $(kwargs...))) => :($F($f).($a, $(args...); $(kwargs...)))
            :($f.($(args...))) => :($F($f).($a, $(args...)))
            :($f($(args...); $(kwargs...))) => :($F($(removedot(f))).($a, $(args...); $(kwargs...)))
            :($f($(args...))) => :($F($(removedot(f))).($a, $(args...)))
            _ => :($F(identity).($a, $b))
        end
    end
end
logical_or(a, b) = a || b
logical_and(a, b) = a && b

"""
    compile_ex(m::Module, ex, info)

Compile a NiLang statement to a regular julia statement.
"""
function compile_ex(m::Module, ex, info)
    @smatch ex begin
        :($a += $b) || :($a .+= $b) ||
        :($a -= $b) || :($a .-= $b) ||
        :($a *= $b) || :($a .*= $b) ||
        :($a /= $b) || :($a ./= $b) ||
        :($a ⊻= $b) || :($a .⊻= $b) => compile_ex(m, to_standard_format(ex), info)
        :(($t1=>$t2)($x)) => assign_ex(x, :(convert($t2, $x)), info.invcheckon[])
        :(($t1=>$t2).($x)) => assign_ex(x, :(convert.($t2, $x)), info.invcheckon[])

        # multi args expanded in preprocessing
        # general
        :($x ↔ $y) => begin
            e1 = isemptyvar(x)
            e2 = isemptyvar(y)
            if e1 && e2
                nothing
            elseif e1 && !e2
                _push_value(x, _pop_value(y), info.invcheckon[])
            elseif !e1 && e2
                _push_value(y, _pop_value(x), info.invcheckon[])
            else
                tmp = gensym("temp")
                Expr(:block, :($tmp = $y), assign_ex(y, x, info.invcheckon[]), assign_ex(x, tmp, info.invcheckon[]))
            end
        end

        # stack
        :($s[end] → $x) => begin
            if info.invcheckon[]
                y = gensym("result")
                Expr(:block, :($y=$loaddata($x, $pop!($s))), _invcheck(y, x), assign_ex(x, y, info.invcheckon[]))
            else
                y = gensym("result")
                Expr(:block, :($y=$loaddata($x, $pop!($s))), assign_ex(x, y, info.invcheckon[]))
            end
        end
        :($s[end+1] ← $x) => :($push!($s, $_copy($x)))
        # dict
        :($x[$index] ← $tp) => begin
            assign_expr = :($x[$index] = $tp)
            if info.invcheckon[]
                Expr(:block, _assert_nokey(x, index), assign_expr)
            else
                assign_expr
            end
        end
        :($x[$index] → $tp) => begin
            delete_expr = :($(deleteindex!)($x, $index))
            if info.invcheckon[]
                Expr(:block, _invcheck(:($x[$index]), tp), delete_expr)
            else
                delete_expr
            end
        end
        # general
        :($x ← $tp) => :($x = $tp)
        :($x → $tp) => begin
            if info.invcheckon[]
                _invcheck(x, tp)
            end
        end

        :($f($(args...))) => begin
            assignback_ex(ex, info.invcheckon[])
        end
        :($f.($(allargs...))) => begin
            args, kwargs = seperate_kwargs(allargs)
            symres = gensym("results")
            ex = :($symres = $unzipped_broadcast($kwargs, $f, $(args...)))
            Expr(:block, ex, assign_vars(args, symres, info.invcheckon[]).args...)
        end
        Expr(:if, _...) => compile_if(m, copy(ex), info)
        :(while ($pre, $post); $(body...); end) => begin
            whilestatement(pre, post, compile_body(m, body, info), info)
        end
        :(for $i=$range; $(body...); end) => begin
            forstatement(i, range, compile_body(m, body, info), info, nothing)
        end
        :(@simd $line for $i=$range; $(body...); end) => begin
            forstatement(i, range, compile_body(m, body, info), info, Symbol("@simd")=>line)
        end
        :(@threads $line for $i=$range; $(body...); end) => begin
            forstatement(i, range, compile_body(m, body, info), info, Symbol("@threads")=>line)
        end
        :(@avx $line for $i=$range; $(body...); end) => begin
            forstatement(i, range, compile_body(m, body, info), info, Symbol("@avx")=>line)
        end
        :(begin $(body...) end) => begin
            Expr(:block, compile_body(m, body, info)...)
        end
        :(@safe $line $subex) => subex
        :(@inbounds $line $subex) => Expr(:macrocall, Symbol("@inbounds"), line, compile_ex(m, subex, info))
        :(@invcheckoff $line $subex) => begin
            state = info.invcheckon[]
            info.invcheckon[] = false
            ex = compile_ex(m, subex, info)
            info.invcheckon[] = state
            ex
        end
        :(@cuda $line $(args...)) => begin
            fcall = @smatch args[end] begin
                :($f($(args...))) => Expr(:call,
                    Expr(:->,
                        :(args...),
                        Expr(:block,
                            :($f(args...)),
                            nothing
                        )
                    ),
                    args...
                )
                _ => error("expect a function after @cuda, got $(args[end])")
            end
            Expr(:macrocall, Symbol("@cuda"), line, args[1:end-1]..., fcall)
        end
        :(@launchkernel $line $device $thread $ndrange $f($(args...))) => begin
            res = gensym("results")
            Expr(:block,
                :($res = $f($device, $thread)($(args...); ndrange=$ndrange)),
                :(wait($res))
            )
        end
        :(nothing) => ex
        ::Nothing => ex
        ::LineNumberNode => ex
        _ => error("statement not supported: `$ex`")
    end
end

function compile_if(m::Module, ex, info)
    pres = []
    posts = []
    ex = analyse_if(m, ex, info, pres, posts)
    Expr(:block, pres..., ex, posts...)
end

function analyse_if(m::Module, ex, info, pres, posts)
    var = gensym("branch")
    if ex.head == :if
        pre, post = ex.args[1].args
        ex.args[1] = var
    elseif ex.head == :elseif
        pre, post = ex.args[1].args[2].args
        ex.args[1].args[2] = var
    end
    push!(pres, :($var = $pre))
    if info.invcheckon[]
        push!(posts, _invcheck(var, post))
    end
    ex.args[2] = Expr(:block, compile_body(m, ex.args[2].args, info)...)
    if length(ex.args) == 3
        if ex.args[3].head == :elseif
            ex.args[3] = analyse_if(m, ex.args[3], info, pres, posts)
        elseif ex.args[3].head == :block
            ex.args[3] = Expr(:block, compile_body(m, ex.args[3].args, info)...)
        end
    end
    ex
end

function whilestatement(precond, postcond, body, info)
    ex = Expr(:block,
        Expr(:while,
            precond,
            Expr(:block, body...),
        ),
    )
    if info.invcheckon[]
        pushfirst!(ex.args, _invcheck(postcond, false))
        push!(ex.args[end].args[end].args,
            _invcheck(postcond, true)
        )
    end
    ex
end

function forstatement(i, range, body, info, mcr)
    assigns, checkers = compile_range(range)
    exf = Expr(:for, :($i=$range), Expr(:block, body...))
    if !(mcr isa Nothing)
        exf = Expr(:macrocall, mcr.first, mcr.second, exf)
    end
    if info.invcheckon[]
        Expr(:block, assigns..., exf, checkers...)
    else
        exf
    end
end

_pop_value(x) = @smatch x begin
    :($s[end]) => :($pop!($s))
    :($s[$ind]) => :($pop!($s, $ind))  # dict (notice pop over vector elements is not allowed.)
    :($x::$T) => :($(_pop_value(x))::$T)
    :(($(args...)),) => Expr(:tuple, _pop_value.(args)...)
    _ => x
end

_push_value(x, val, invcheck) = @smatch x begin
    :($s[end+1]) => :($push!($s, $val))
    :($s[$arg]::∅) => begin
        ex = :($s[$arg] = $val)
        if invcheck
            Expr(:block, _assert_nokey(s, arg), ex)
        else
            ex
        end
    end
    _ => assign_ex(x, val, invcheck)
end

function _assert_nokey(x, index)
    str = "dictionary `$x` already has key `$index`"
    Expr(:if, :(haskey($x, $index)), :(throw(InvertibilityError($str))))
end

_copy(x) = copy(x)
_copy(x::Tuple) = copy.(x)

export @code_julia

"""
    @code_julia ex

Get the interpreted expression of `ex`.

```julia
julia> @code_julia x += exp(3.0)
quote
    var"##results#267" = ((PlusEq)(exp))(x, 3.0)
    x = var"##results#267"[1]
    try
        (NiLangCore.deanc)(3.0, var"##results#267"[2])
    catch e
        @warn "deallocate fail: `3.0 → var\"##results#267\"[2]`"
        throw(e)
    end
end

julia> @code_julia @invcheckoff x += exp(3.0)
quote
    var"##results#257" = ((PlusEq)(exp))(x, 3.0)
    x = var"##results#257"[1]
end
```
"""
macro code_julia(ex)
    QuoteNode(compile_ex(__module__, ex, CompileInfo()))
end

compile_ex(m::Module, ex) = compile_ex(m, ex, CompileInfo())

export @i

"""
    @i function fname(args..., kwargs...) ... end
    @i struct sname ... end

Define a reversible function/type.

```jldoctest; setup=:(using NiLangCore)
julia> @i function test(out!, x)
           out! += identity(x)
       end

julia> test(0.2, 0.8)
(1.0, 0.8)
```
See `test/compiler.jl` for more examples.
"""
macro i(ex)
    ex = gen_ifunc(__module__, ex)
    ex.args[1] = :(Base.@__doc__ $(ex.args[1]))
    esc(ex)
end

# generate the reversed function
function gen_ifunc(m::Module, ex)
    mc, fname, args, ts, body = precom(m, ex)
    fname = _replace_opmx(fname)

    # implementations
    ftype = get_ftype(fname)

    head = :($fname($(args...)) where {$(ts...)})
    dfname = dual_fname(fname)
    dftype = get_ftype(dfname)
    fdef1 = Expr(:function, head, Expr(:block, compile_body(m, body, CompileInfo())..., functionfoot(args)))
    dualhead = :($dfname($(args...)) where {$(ts...)})
    fdef2 = Expr(:function, dualhead, Expr(:block, compile_body(m, dual_body(m, body), CompileInfo())..., functionfoot(args)))
    if mc !== nothing
        fdef1 = Expr(:macrocall, mc[1], mc[2], fdef1)
        fdef2 = Expr(:macrocall, mc[1], mc[2], fdef2)
    end
    #ex = :(Base.@__doc__ $fdef1; if $ftype != $dftype; $fdef2; end)
    ex = Expr(:block, fdef1,
        Expr(:if, :($ftype != $dftype), fdef2),
        )
end

export nilang_ir
"""
    nilang_ir(ex; reversed::Bool=false)

Get the NiLang reversible IR from the function expression `ex`,
return the reversed function if `reversed` is `true`.

This IR is not directly executable on Julia, please use
`macroexpand(Main, :(@i function .... end))` to get the julia expression of a reversible function.

```jldoctest; setup=:(using NiLangCore)
julia> ex = :(@inline function f(x!::T, y) where T
                @routine begin
                    anc ← zero(T)
                    anc += identity(x!)
                end
                x! += y * anc
                ~@routine
           end);

julia> NiLangCore.nilang_ir(Main, ex) |> NiLangCore.rmlines
:(@inline function f(x!::T, y) where T
          begin
              anc ← zero(T)
              anc += identity(x!)
          end
          x! += y * anc
          begin
              anc -= identity(x!)
              anc → zero(T)
          end
      end)

julia> NiLangCore.nilang_ir(Main, ex; reversed=true) |> NiLangCore.rmlines
:(@inline function (~f)(x!::T, y) where T
          begin
              anc ← zero(T)
              anc += identity(x!)
          end
          x! -= y * anc
          begin
              anc -= identity(x!)
              anc → zero(T)
          end
      end)
```
"""
function nilang_ir(m::Module, ex; reversed::Bool=false)
    mc, fname, args, ts, body = precom(m, ex)
    fname = _replace_opmx(fname)

    # implementations
    if reversed
        dfname = :(~$fname) # use fake head for readability
        head = :($dfname($(args...)) where {$(ts...)})
        body = dual_body(m, body)
    else
        head = :($fname($(args...)) where {$(ts...)})
    end
    fdef = Expr(:function, head, Expr(:block, body...))
    if mc !== nothing
        fdef = Expr(:macrocall, mc[1], mc[2], fdef)
    end
    fdef
end

# seperate and return `args` and `kwargs`
@inline function seperate_kwargs(args)
    if length(args) > 0 && args[1] isa Expr && args[1].head == :parameters
        args = args[2:end], args[1]
    else
        args, Expr(:parameters)
    end
end

# add a `return` statement to the end of the function body.
function functionfoot(args)
    args = get_argname.(seperate_kwargs(args)[1])
    if length(args) == 1
        if args[1] isa Expr && args[1].head == :(...)
            args[1].args[1]
        else
            args[1]
        end
    else
        :(($(args...),))
    end
end

# to provide the eye candy for defining `x += f(args...)` like functions
_replace_opmx(ex) = @smatch ex begin
    :(:+=($f)) => :($(gensym())::PlusEq{typeof($f)})
    :(:-=($f)) => :($(gensym())::MinusEq{typeof($f)})
    :(:*=($f)) => :($(gensym())::MulEq{typeof($f)})
    :(:/=($f)) => :($(gensym())::DivEq{typeof($f)})
    :(:⊻=($f)) => :($(gensym())::XorEq{typeof($f)})
    _ => ex
end

export @instr

"""
    @instr ex

Execute a reversible instruction.
"""
macro instr(ex)
    ex = precom_ex(__module__, ex, NiLangCore.PreInfo())
    #variable_analysis_ex(ex, SymbolTable())
    esc(Expr(:block, NiLangCore.compile_ex(__module__, ex, CompileInfo()), nothing))
end

# the range of for statement
compile_range(range) = @smatch range begin
    :($start:$step:$stop) => begin
        start_, step_, stop_ = gensym("start"), gensym("step"), gensym("stop")
        Any[:($start_ = $start),
        :($step_ = $step),
        :($stop_ = $stop)],
        Any[_invcheck(start_, start),
        _invcheck(step_, step),
        _invcheck(stop_, stop)]
    end
    :($start:$stop) => begin
        start_, stop_ = gensym("start"), gensym("stop")
        Any[:($start_ = $start),
        :($stop_ = $stop)],
        Any[_invcheck(start_, start),
        _invcheck(stop_, stop)]
    end
    :($list) => begin
        list_ = gensym("iterable")
        Any[:($list_ = deepcopy($list))],
        Any[_invcheck(list_, list)]
    end
end

"""
    get_ftype(fname)

Return the function type, e.g.

* `obj::ABC` => `ABC`
* `f` => `typeof(f)`
"""
function get_ftype(fname)
    @smatch fname begin
        :($x::$tp) => tp
        _ => :($NiLangCore._typeof($fname))
    end
end


text:
"""
    statustext(::Int) -> String

`String` representation of a HTTP status code.

## Examples
```julia
julia> statustext(200)
"OK"

julia> statustext(404)
"Not Found"
```
"""
statustext(status) = Base.get(STATUS_MESSAGES, status, "Unknown Code")

const STATUS_MESSAGES = (()->begin
    v = fill("Unknown Code", 530)
    v[100] = "Continue"
    v[101] = "Switching Protocols"
    v[102] = "Processing"                            # RFC 2518 => obsoleted by RFC 4918
    v[103] = "Early Hints"
    v[200] = "OK"
    v[201] = "Created"
    v[202] = "Accepted"
    v[203] = "Non-Authoritative Information"
    v[204] = "No Content"
    v[205] = "Reset Content"
    v[206] = "Partial Content"
    v[207] = "Multi-Status"                          # RFC 4918
    v[208] = "Already Reported"                      # RFC5842
    v[226] = "IM Used"                               # RFC3229
    v[300] = "Multiple Choices"
    v[301] = "Moved Permanently"
    v[302] = "Moved Temporarily"
    v[303] = "See Other"
    v[304] = "Not Modified"
    v[305] = "Use Proxy"
    v[307] = "Temporary Redirect"
    v[308] = "Permanent Redirect"                    # RFC7238
    v[400] = "Bad Request"
    v[401] = "Unauthorized"
    v[402] = "Payment Required"
    v[403] = "Forbidden"
    v[404] = "Not Found"
    v[405] = "Method Not Allowed"
    v[406] = "Not Acceptable"
    v[407] = "Proxy Authentication Required"
    v[408] = "Request Time-out"
    v[409] = "Conflict"
    v[410] = "Gone"
    v[411] = "Length Required"
    v[412] = "Precondition Failed"
    v[413] = "Request Entity Too Large"
    v[414] = "Request-URI Too Large"
    v[415] = "Unsupported Media Type"
    v[416] = "Requested Range Not Satisfiable"
    v[417] = "Expectation Failed"
    v[418] = "I'm a teapot"                        # RFC 2324
    v[421] = "Misdirected Request"                 # RFC 7540
    v[422] = "Unprocessable Entity"                # RFC 4918
    v[423] = "Locked"                              # RFC 4918
    v[424] = "Failed Dependency"                   # RFC 4918
    v[425] = "Unordered Collection"                # RFC 4918
    v[426] = "Upgrade Required"                    # RFC 2817
    v[428] = "Precondition Required"               # RFC 6585
    v[429] = "Too Many Requests"                   # RFC 6585
    v[431] = "Request Header Fields Too Large"     # RFC 6585
    v[440] = "Login Timeout"
    v[444] = "nginx error: No Response"
    v[451] = "Unavailable For Legal Reasons"       # RFC7725
    v[495] = "nginx error: SSL Certificate Error"
    v[496] = "nginx error: SSL Certificate Required"
    v[497] = "nginx error: HTTP -> HTTPS"
    v[499] = "nginx error or Antivirus intercepted request or ArcGIS error"
    v[500] = "Internal Server Error"
    v[501] = "Not Implemented"
    v[502] = "Bad Gateway"
    v[503] = "Service Unavailable"
    v[504] = "Gateway Time-out"
    v[505] = "HTTP Version Not Supported"
    v[506] = "Variant Also Negotiates"             # RFC 2295
    v[507] = "Insufficient Storage"                # RFC 4918
    v[508] = "Loop Detected"                       # RFC5842
    v[509] = "Bandwidth Limit Exceeded"
    v[510] = "Not Extended"                        # RFC 2774
    v[511] = "Network Authentication Required"     # RFC 6585
    v[520] = "CloudFlare Server Error: Unknown"
    v[521] = "CloudFlare Server Error: Connection Refused"
    v[522] = "CloudFlare Server Error: Connection Timeout"
    v[523] = "CloudFlare Server Error: Origin Server Unreachable"
    v[524] = "CloudFlare Server Error: Connection Timeout"
    v[525] = "CloudFlare Server Error: Connection Failed"
    v[526] = "CloudFlare Server Error: Invalid SSL Ceritificate"
    v[527] = "CloudFlare Server Error: Railgun Error"
    v[530] = "Site Frozen"
    return v
end)()

ext: jl


text:
using ExponentialAction, ChainRulesCore, ChainRulesTestUtils, LinearAlgebra, Test
using ExponentialAction: _opnormInf

@testset "Utilities" begin
    @testset "_opnormInf" begin
        A = randn(ComplexF64, 10)
        @test _opnormInf(A) ≈ _opnormInf(reshape(A, 10, 1))
        @test _opnormInf(A) ≈ norm(A, Inf)
        A = randn(ComplexF64, 10, 10)
        @test _opnormInf(A) ≈ opnorm(A, Inf)
    end

    @testset "_opnormInf non-differentiable" begin
        test_rrule(_opnormInf, randn(ComplexF64, 10) ⊢ NoTangent(); atol=1e-6)
        test_rrule(_opnormInf, randn(ComplexF64, 10, 10) ⊢ NoTangent(); atol=1e-6)
    end

    @testset "opnormest1" begin
        A = randn(ComplexF64, 10, 10)
        @test ExponentialAction.opnormest1(A) ≈ opnorm(A, 1)
    end
end


text:
a = [1, 2]
A = NonlinearBandits.action_matrix(a, 2)
@test A == [1.0 0.0; 0.0 1.0]

d = 5
num_actions = 2
layer_sizes = [64, 64]
encoder = NeuralEncoder(d, num_actions, layer_sizes)
X = rand(5, 2)
@test size(encoder(X)) == (layer_sizes[end], 2)

y = rand(1, 2)
Ypred = cpu(encoder.nn(gpu(X)))
ls = (y[1, 1] - Ypred[a[1], 1])^2 + (y[1, 2] - Ypred[a[2], 2])^2
ls /= 2
@test encoder.loss(gpu(X), gpu(A), gpu(y)) ≈ ls

# Create bandit trajectory
d, n = 2, 5000
mf = (x -> 3 * x[1] + 5 * x[2], x -> 4 * x[1])
X = rand(d, n)
a = rand(1:2, n)
r = zeros(1, n)
for i in 1:n
    r[1, i] = mf[a[i]](X[:, i]) + rand(Normal(0, 0.1))
end

enc = NeuralEncoder(d, length(mf), [32, 32])
fit!(enc, X, a, r, 10; verbose=false)
Ypred = cpu(enc.nn(gpu(X)))
@test mae(Ypred[1:1, :], mapslices(mf[1], X; dims=1)) < 0.1
@test mae(Ypred[2:2, :], mapslices(mf[2], X; dims=1)) < 0.1

# Test updating NeuralLinear policy appears to work
limits = repeat([-1.0 1.0], d, 1)
csampler = UniformContexts(limits)
rsampler = GaussianRewards(mf)
metrics = (FunctionalRegret(mf),)

# Test NeuralLinear appears to be learning the optimal policy.
batch_size = 1000
num_batches = 3
initial_batches = retrain_freq = 1
nl = NeuralLinear(d, length(mf), layer_sizes, initial_batches, retrain_freq, 20)
metrics = (FunctionalRegret(mf),)
driver = StandardDriver(csampler, nl, rsampler, metrics)
run!(num_batches, batch_size, driver)

regret = metrics[1].regret[(end - batch_size + 1):end]
@test sum(regret .== 0.0) / length(regret) >= 0.95


text:
"""
This file intends to genrate Julia methods according C++ function:

void atg_abs_out(tensor *out__, tensor out, tensor self) {
  ...
}

->

export abs_out
function abs_out(out::Tensor, self::Tensor)
    outputs = Int[0]
    ccall((:atg_abs_out, libtorch_capi),
          Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
          pointer(outputs), out.pointer, self.pointer)
    return tensor_from_ptr(Ptr{Cvoid}(outputs[1]))
end
"""

const PROJECT_DIR = (@__DIR__) |> dirname |> dirname
const CPP_API_FILE = joinpath(PROJECT_DIR, "csrc", "torch_api_generated.cpp.h")
const JUL_API_FILE = joinpath(PROJECT_DIR, "src", "thc", "thc.jl")

const FUNC_SIG_REG = r"(\w+)\s+(\*?atg_\w+)\((.+)\)\s*{"

const JULIA_KEYWORDS = Set(["function", "end"])

const C_TYPE_MAP = Dict(
    "void"     => "Cvoid",
    "tensor*"  => "Ptr{Cvoid}",
    "tensor"   => "Ptr{Cvoid}",
    "scalar"   => "Ptr{Cvoid}",
    "int"      => "Cint",
    "int*"     => "Ptr{Cvoid}",
    "int64_t"  => "Clonglong",
    "int64_t*" => "Ptr{Cvoid}",
    "double"   => "Cdouble",
)

const J_TYPE_MAP = Dict(
    "void"     => "Any",
    "tensor*"  => "Array{Tensor{T,N}}",
    "tensor"   => "Tensor",
    "scalar"   => "TorchNumber",
    "int"      => "Int",
    "int*"     => "Array{Int}",
    "int64_t"  => "Int64",
    "int64_t*" => "Array{Int64}",
    "double"   => "Float64",
)

struct APIFunction
    cpp_signature::String
    func_name::String
    return_type::String
    output_count::Int
    args::Vector{Pair{String, String}} # name -> type

    function APIFunction(m::RegexMatch, n::Int)
        csig = strip(m.match, [' ', '{', '\n'])
        args = parse_args(m[3])
        ret_type = m[1]
        fname = m[2]
        if fname[1] == '*'
            fname = fname[6:end]
            ret_type *= '*'
        else
            fname = fname[5:end]
        end
        new(csig, fname, ret_type, n, args)
    end
end

function parse_args(args::AbstractString)
    arg_list = strip.(split(args, ','))
    arg_pairs = map(arg_list) do arg
        info = strip.(split(arg)) # type, name
        if info[2][1] == '*'
            info[1] *= '*'
            info[2] = info[2][2:end]
        end
        info[2] in JULIA_KEYWORDS && (info[2] *= "_")
        info[2] => info[1]
    end
    arg_pairs
end

function julia_source(f::APIFunction)
    if  length(f.args) < 1
        @warn "E1: can't translate function [$(f.cpp_signature)], ignored."
        return "# $(f.func_name) ignored"
    end

    for arg in f.args
        if !haskey(C_TYPE_MAP, arg.second)
            @warn "E2: can't translate function [$(f.cpp_signature)], ignored."
            return "# $(f.func_name) ignored"
        end
    end

    lines = [""]
    # in-place op: pow_ -> pow!, pow_1 -> pow1!, ...
    jl_fname = f.func_name
    sufix_m = match(r"(\w+[^_])_(\d*)$", jl_fname)
    sufix_m != nothing && (jl_fname = "$(sufix_m[1])$(sufix_m[2])!")

    if Symbol(jl_fname) in names(Base)
        if !in(jl_fname, ["div"])
            push!(lines, "import Base.$(jl_fname)")
        end
    end

    push!(lines, doc(f, jl_fname)) # docs

    start = f.args[1].first == "out__" ? 2 : 1
    para_type = any(x -> x.second == "tensor*", f.args[start:end]) ?
        " where {T,N}" : ""
    ccall_ret = start == 1 ? "Ptr{Int}" : "Cvoid"

    push!(lines, "function $(jl_fname)($(julia_args(f)))$(para_type)")
    push!(lines, julia_locals(f))
    push!(lines, "    __cret = ccall((:atg_$(f.func_name), libtorch_capi),")
    push!(lines, "                 $(ccall_ret), ($(ccall_args(f))),")
    push!(lines, "                 $(ccall_julia_args(f)))")
    push!(lines, return_statement(f))
    push!(lines, "end")

    return join(lines, "\n")
end

function doc(f::APIFunction, jl_fname::AbstractString)
    cpp_sig = replace(f.cpp_signature, "_" => "\\\\_")
    lines = ["\n"]
    push!(lines, "\"\"\"")
    push!(lines, "    $(jl_fname)($(julia_args(f)))")
    push!(lines, "")
    push!(lines, " Wrapper of C++ function $(cpp_sig)")
    push!(lines, "\"\"\"")
    join(lines, "\n")
end

function julia_args(f::APIFunction)
    args = []
    start = f.args[1].first == "out__" ? 2 : 1
    for i in start:length(f.args)
        p = f.args[i]
        if endswith(p.first, "_len") && endswith(f.args[i-1].first, "_data")
            nothing
        else
            push!(args, "$(p.first)::$(J_TYPE_MAP[p.second])")
        end
    end
    join(args, ", ")
end

function julia_locals(f::APIFunction)
    lines = []
    for i in 1:length(f.args)
        p = f.args[i]
        if endswith(p.first, "_len") && endswith(f.args[i-1].first, "_data")
            push!(lines, "    $(p.first) = length($(f.args[i-1].first))")
        elseif p.second == "scalar"
            push!(lines, "    $(p.first)_s_ = Scalar($(p.first))")
        elseif p.second == "tensor*"
            if p.first == "out__"
                output_init = join(repeat(["0"], f.output_count), ", ")
                push!(lines, "    outputs__ = Int[$(output_init)]")
            else
                push!(lines, "    $(p.first)_ta_ = map(x->x.pointer, $(p.first))")
            end
        end
    end
    join(lines, "\n")
end

function ccall_args(f::APIFunction)
    length(f.args) == 1 && return C_TYPE_MAP[f.args[1].second] * ","
    args = map(p -> C_TYPE_MAP[p.second], f.args)
    join(args, ", ")
end

function ccall_julia_args(f::APIFunction)
    args = map(f.args) do p
        p.second == "tensor*" && p.first == "out__" && return "outputs__"
        p.second == "tensor*" && return "$(p.first)_ta_"
        p.second == "tensor"  && return "$(p.first).pointer"
        p.second == "scalar"  && return "$(p.first)_s_.pointer"
        return p.first
    end
    join(args, ", ")
end

function return_statement(f::APIFunction)
    if match(r"_\d*$", f.func_name) != nothing
        return "    return self"
    elseif f.return_type == "void" && f.args[1].first == "out__"
        lines = []
        for i in 1:f.output_count
            push!(lines,
                  "    __o_$(i) = tensor_from_ptr(Ptr{Cvoid}(outputs__[$(i)]))")
        end
        push!(lines,
              "    return " * join(map(x-> "__o_$x", 1:f.output_count), ", "))
        return join(lines, "\n")
    elseif f.return_type == "tensor*"
        lines = []
        push!(lines, "    ptrs__, i__ = Int[], 1")
        push!(lines, "    while true")
        push!(lines, "        ptr__ = unsafe_load(__cret, i__)")
        push!(lines, "        ptr__ == 0 && break")
        push!(lines, "        push!(ptrs__, ptr__)")
        push!(lines, "        i__ += 1")
        push!(lines, "    end")
        push!(lines, "    ccall(:free, Cvoid, (Ptr{Cvoid},), __cret)")
        push!(lines, "    return map(x -> tensor_from_ptr(Ptr{Nothing}(x)), ptrs__)")
        return join(lines, "\n")
    end
    return ""
end

function main()
    count = 0
    source_lines = readlines(CPP_API_FILE)
    output = open(JUL_API_FILE, "w")
    func_match = nothing
    output_count = 0

    write(output, "# !!! THIS FILE IS AUTO-GENERATED, PLEASE DO NOT MODIFY. !!!\n\n")
    write(output, "module ThC\n") # module start
    write(output, """
@static if Sys.islinux()
    using LibTorchCAPI_jll
elseif Sys.isapple()
    const libtorch_capi = :libtorch_capi
end
""")
    write(output, "using ..ThArrays: Tensor, Scalar, TorchNumber, tensor_from_ptr\n")

    for line in source_lines
        m = match(FUNC_SIG_REG, line)
        if m != nothing # start of a function
            if func_match != nothing # deal with the previous function
                f = APIFunction(func_match, output_count)
                write(output, julia_source(f))
                output_count = 0
                count += 1
            end
            func_match = m
        end

        if func_match != nothing # in a function
            if match(r"out__\[\d+\]\s*=\s*new", line) != nothing
                output_count += 1
            elseif match(r"out__\[\D+\]\s*=\s*new", line) != nothing
                output_count = 1
            end
        end
    end

    if func_match != nothing # the last function
        f = APIFunction(func_match, output_count)
        write(output, julia_source(f))
        count += 1
    end

    write(output, "\n")
    write(output, "include(\"thc-opt.jl\")\n")
    write(output, "\n")
    write(output, "end\n") # module end

    close(output)

    @info "$(count) methods generated!\n"
end

main()


text:
# to test vtk-files
OVERWRITE_CHECKSUMS = false
checksums_file = joinpath(dirname(@__FILE__), "checksums.sha1")
checksum_list = read(checksums_file, String)
if OVERWRITE_CHECKSUMS
    csio = open(checksums_file, "w")
else
    csio = open(checksums_file, "r")
end

@testset "Grid, DofHandler, vtk" begin
    for (celltype, dim) in ((Line,                   1),
                            (QuadraticLine,          1),
                            (Quadrilateral,          2),
                            (QuadraticQuadrilateral, 2),
                            (Triangle,               2),
                            (QuadraticTriangle,      2),
                            (Hexahedron,             3),
                            (Tetrahedron,            3))

        # create test grid, do some operations on it and then test
        # the resulting sha1 of the stored vtk file
        # after manually checking the exported vtk
        nels = ntuple(x->5, dim)
        right = Vec{dim, Float64}(ntuple(x->1.5, dim))
        left = -right
        grid = generate_grid(celltype, nels, left, right)

        transform!(grid, x-> 2x)

        radius = 2*1.5
        addcellset!(grid, "cell-1", [1,])
        addcellset!(grid, "middle-cells", x -> norm(x) < radius)
        addfaceset!(grid, "middle-faceset", x -> norm(x) < radius)
        addfaceset!(grid, "right-faceset", getfaceset(grid, "right"))
        addnodeset!(grid, "middle-nodes", x -> norm(x) < radius)

        gridfilename = "grid-$(Ferrite.celltypes[celltype])"
        vtk_grid(gridfilename, grid) do vtk
            vtk_cellset(vtk, grid, "cell-1")
            vtk_cellset(vtk, grid, "middle-cells")
            vtk_nodeset(vtk, grid, "middle-nodes")
        end

        # test the sha of the file
        sha = bytes2hex(open(SHA.sha1, gridfilename*".vtu"))
        if OVERWRITE_CHECKSUMS
            write(csio, sha, "\n")
        else
            @test chomp(readline(csio)) == sha
            rm(gridfilename*".vtu")
        end

        # Create a DofHandler, add some things, write to file and
        # then check the resulting sha
        dofhandler = DofHandler(grid)
        push!(dofhandler, :temperature, 1)
        push!(dofhandler, :displacement, 3)
        close!(dofhandler)
        ch = ConstraintHandler(dofhandler)
        dbc = Dirichlet(:temperature, union(getfaceset(grid, "left"), getfaceset(grid, "right-faceset")), (x,t)->1)
        add!(ch, dbc)
        dbc = Dirichlet(:temperature, getfaceset(grid, "middle-faceset"), (x,t)->4)
        add!(ch, dbc)
        for d in 1:dim
            dbc = Dirichlet(:displacement, union(getfaceset(grid, "left")), (x,t) -> d, d)
            add!(ch, dbc)
        end
        close!(ch)
        update!(ch, 0.0)
        Random.seed!(1234)
        u = rand(ndofs(dofhandler))
        apply!(u, ch)

        dofhandlerfilename = "dofhandler-$(Ferrite.celltypes[celltype])"
        vtk_grid(dofhandlerfilename, dofhandler) do vtk
            vtk_point_data(vtk, ch)
            vtk_point_data(vtk, dofhandler, u)
        end

        # test the sha of the file
        sha = bytes2hex(open(SHA.sha1, dofhandlerfilename*".vtu"))
        if OVERWRITE_CHECKSUMS
            write(csio, sha, "\n")
        else
            @test chomp(readline(csio)) == sha
            rm(dofhandlerfilename*".vtu")
        end

    end

end # of testset

if OVERWRITE_CHECKSUMS
    close(csio)
end

# right = Vec{2, Float64}(ntuple(x->1.5, dim))
# left = -right

@testset "Grid utils" begin

    grid = Ferrite.generate_grid(QuadraticQuadrilateral, (1, 1), Vec((0.,0.)), Vec((1.,1.)))

    addcellset!(grid, "cell_set", [1]);
    node_set = Set(1:getnnodes(grid))
    addnodeset!(grid, "node_set", node_set)

    @test getnodesets(grid) == Dict("node_set" => node_set)

    @test getnodes(grid, [1]) == [getnodes(grid, 1)] # untested

    @test length(getnodes(grid, "node_set")) == 9

    @test collect(getcoordinates(getnodes(grid, 5)).data) ≈ [0.5, 0.5]

    @test getcells(grid, "cell_set") == [getcells(grid, 1)]

    f(x) = Tensor{1,1,Float64}((1 + x[1]^2 + 2x[2]^2, ))

    values = compute_vertex_values(grid, f)
    @test f([0.0, 0.0]) == values[1]
    @test f([0.5, 0.5]) == values[5]
    @test f([1.0, 1.0]) == values[9]

    @test compute_vertex_values(grid, collect(1:9), f) == values

    # Can we test this in a better way? The set makes the order random.
    @test length(compute_vertex_values(grid, "node_set", f)) == 9

    # CellIterator on a grid without DofHandler
    grid = generate_grid(Triangle, (4,4))
    n = 0
    ci = CellIterator(grid)
    @test length(ci) == getncells(grid)
    for c in ci
        getcoordinates(c)
        getnodes(c)
        n += cellid(c)
    end
    @test n == div(getncells(grid)*(getncells(grid) + 1), 2)
end


text:
function serial_pi_estimation(n_samples)
    n_inside = 0;

    x = zeros(2);

    for _ in 1:n_samples
        @. x = rand();
        if(norm(x,2)<1)
            n_inside+=1;
        end
    end

    pi_estimate = 4 * n_inside/n_samples;

    return pi_estimate

end

function pi_estimation(n_samples)
    n_inside = 0;

    x_inside = zeros(Bool, n_samples);

    x = zeros(2);
    Threads.@threads for j in 1:n_samples
        @. x = rand();
        if(norm(x,2)<1)
            x_inside[j] = 1;
        end
    end

    n_inside = sum(x_inside);

    pi_estimate = 4 * n_inside/n_samples;

    return pi_estimate

end

function atomic_pi_estimation(n_samples)

    n_inside = Atomic{Int}(0);

    x = zeros(2);

    Threads.@threads for _ in 1:n_samples
        @. x = rand();
        if(norm(x,2)<1)
            atomic_add!(n_inside,1);
        end
    end

    pi_estimate = 4 * n_inside[]/n_samples;

    return pi_estimate

end


text:
@testset "Det" begin
    Mod.@polyvar x y

    @test det([x 1; y 1]) == x - y
    @test det([x 1; x 1]) == 0
    @test det([x 1 1 1; x y 1 2; 0 0 0 1; x 0 y 0]) == -x*y^2 + 2*x*y - x
end


text:
export SiteLabel, BitSiteLabel, Bit, Spin, Half
export up, down

##############
# Site Label
##############

"""
    SiteLabel

Site Labels, type of each element, like Spin, Bit, etc.
"""
abstract type SiteLabel <: AbstractTag end

"""
    BitSiteLabel <: SiteLabel

Binary Site, like Spin, Bit, etc.
"""
abstract type BitSiteLabel <: SiteLabel end

#################################################################
# Interfaces
import Base: eltype, is

"""
    eltype(BitSiteLabel)

get tag datatype.
"""
function eltype(::Type{T}) where {T <: BitSiteLabel} end

"""
    up(BitSiteLabel) -> tag

up tag for this label. e.g. `1` for `Bit`, `0.5` for `Half`.
"""
function up(::Type{T}) where {T <: BitSiteLabel} end

"""
    down(BitSiteLabel) -> tag

down tag for this label. e.g. `0` for `Bit`, `-0.5` for `Half`.
"""
function down(::Type{T}) where {T <: BitSiteLabel} end

"""
    is(x, BitSiteLabel) -> Bool

check if x is of given BitSiteLabel
"""
is(::Type{T}, x::Number) where {T <: BitSiteLabel} =
    x == up(T) ? true : x == down(T) ? true : false
is(::Type{T}, x::NTuple{N, D}) where {T <: BitSiteLabel, D <: Number, N} =
    all(is(T, each) for each in x)

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

"""
    Spin <: BitSiteLabel

Binary State, Spin. Has two state: -1, 1
"""
abstract type Spin <: BitSiteLabel end

eltype(::Type{Spin}) = Float32
up(::Type{Spin}) = convert(eltype(Spin), 1.0f0)
down(::Type{Spin}) = convert(eltype(Spin), -1.0f0)

"""
    Bit <: BitSiteLabel

Binary State, Bit. Has two state: 0, 1
"""
abstract type Bit <: BitSiteLabel end

eltype(::Type{Bit}) = Float32
up(::Type{Bit}) = convert(eltype(Bit), 1.0f0)
down(::Type{Bit}) = convert(eltype(Bit), 0.0f0)

"""
    Half <: BitSiteLabel

Binary State, Half. Has two state: -0.5, 0.5
"""
abstract type Half <: BitSiteLabel end

eltype(::Type{Half}) = Float32
up(::Type{Half}) = convert(eltype(Half), 0.5f0)
down(::Type{Half}) = convert(eltype(Half), -0.5f0)


text:
#functionality in this file implements writing .vtu files
#complient to https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf
import Base64
#import CodecZlib
##
function vtk_write1(file_name::String,mesh::Mesh,data=[],binary=true,compressed=false)
   open(file_name*".vtu","w") do fid
      if compressed
         compressor=" compressor=\"vtkZLibDataCompressor\""
      else
         compressor=""
      end
      write(fid,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" header_type=\"UInt64\"$compressor>\n")
      write(fid,"\t<UnstructuredGrid>\n")
      write(fid,"\t\t<Piece NumberOfPoints=\"$(size(mesh.points,2))\" NumberOfCells=\"$(length(mesh.tetrahedra))\">\n")
      write(fid,"\t\t\t<Points>\n")
      #write(fid,"\t\t\t\t<DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">")
      #for pnt in mesh.points
      #   write(fid,string(pnt)*" ")
      #end
      #write(fid,"</DataArray>")
      write_data_array(fid,Dict("Coordinates"=>mesh.points),binary,compressed)
      write(fid,"\t\t\t</Points>\n")
      write(fid,"\t\t\t<Cells>\n")
      # write(fid,"\t\t\t\t<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">")
      # for tet in mesh.tetrahedra
      #    for pnt in tet
      #       write(fid,string(Int32(pnt-1))*" ")
      #    end
      # end
      # write(fid,"</DataArray>\n")
      N_points=size(mesh.points,2)
      #connectivioty must be Int32. Paraview is serious about this in its latest
      #versions! https://discourse.paraview.org/t/cannot-open-vtu-files-with-paraview-5-8/3759
      connectivity=Array{Int32}(undef,length(mesh.tetrahedra)*4)
      idx=1
      for tet in mesh.tetrahedra
         for pnt_idx in tet
            connectivity[idx]=Int32(pnt_idx-1)
            idx+=1
         end
      end

      write_data_array(fid,Dict("connectivity"=>connectivity),binary,compressed)

      write(fid,"\t\t\t\t<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">")
      off=0
      for tet in mesh.tetrahedra
            off+=length(tet)
            write(fid,string(off)*" ")
      end
      write(fid,"</DataArray>\n")
      write(fid,"\t\t\t\t<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">")
      for tet in mesh.tetrahedra
         for pnt in tet
            write(fid,"10 ")
         end
      end
      write(fid,"</DataArray>\n")
      write(fid,"\t\t\t</Cells>\n")

      write(fid,"\t\t\t<PointData>\n")
      write_data_array(fid,data,binary,compressed)
      write(fid,"\t\t\t</PointData>\n")
      write(fid,"\t\t</Piece>\n")
      write(fid,"\t</UnstructuredGrid>\n")
      write(fid,"</VTKFile>\n")
   end
   return nothing
end
##
function vtk_write0(file_name::String,mesh::Mesh,data=[],binary=true,compressed=false)
   open(file_name*".vtu","w") do fid
      if compressed
         compressor=" compressor=\"vtkZLibDataCompressor\""
      else
         compressor=""
      end
      write(fid,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" header_type=\"UInt64\"$compressor>\n")
      write(fid,"\t<UnstructuredGrid>\n")
      write(fid,"\t\t<Piece NumberOfPoints=\"$(size(mesh.points,2))\" NumberOfCells=\"$(length(mesh.tetrahedra))\">\n")
      write(fid,"\t\t\t<Points>\n")
      # write(fid,"\t\t\t\t<DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">")
      # for pnt in mesh.points
      #    write(fid,string(pnt)*" ")
      # end
      # write(fid,"</DataArray>")
      write_data_array(fid,Dict("Coordinates"=>mesh.points),binary,compressed)
      write(fid,"\t\t\t</Points>\n")
      write(fid,"\t\t\t<Cells>\n")
      write(fid,"\t\t\t\t<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">")
      for tet in mesh.tetrahedra
         for pnt in tet
            write(fid,string(Int32(pnt-1))*" ")
         end
      end
      write(fid,"</DataArray>\n")
      write(fid,"\t\t\t\t<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">")
      off=0
      for tet in mesh.tetrahedra
            off+=length(tet)
            write(fid,string(off)*" ")
      end
      write(fid,"</DataArray>\n")
      write(fid,"\t\t\t\t<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">")
      for tet in mesh.tetrahedra
         for pnt in tet
            write(fid,"10 ")
         end
      end
      write(fid,"</DataArray>\n")
      write(fid,"\t\t\t</Cells>\n")

      write(fid,"\t\t\t<CellData>\n")
      write_data_array(fid,data,binary,compressed)
      write(fid,"\t\t\t</CellData>\n")
      write(fid,"\t\t</Piece>\n")
      write(fid,"\t</UnstructuredGrid>\n")
      write(fid,"</VTKFile>\n")
   end
   return nothing
end

function vtk_write2(file_name::String,mesh::Mesh,data=[],binary=true,compressed=false)
   open(file_name*".vtu","w") do fid
      if compressed
         compressor=" compressor=\"vtkZLibDataCompressor\""
      else
         compressor=""
      end
      write(fid,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" header_type=\"UInt64\"$compressor>\n")
      write(fid,"\t<UnstructuredGrid>\n")
      write(fid,"\t\t<Piece NumberOfPoints=\"$(size(mesh.points,2)+length(mesh.lines))\" NumberOfCells=\"$(length(mesh.tetrahedra))\">\n")
      write(fid,"\t\t\t<Points>\n")
      #write(fid,"\t\t\t\t<DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">")
      # for pnt in mesh.points
      #    write(fid,string(pnt)*" ")
      # end
      # for line in mesh.lines
      #    pnts=sum(mesh.points[:,line],dims=2)/2
      #    for pnt in pnts
      #       write(fid,string(pnt)*" ")
      #    end
      # end
      #write(fid,"</DataArray>")
      N_points=size(mesh.points,2)
      points=Array{Float64}(undef,3,N_points+length(mesh.lines))
      points[:,1:N_points]=mesh.points
      for (idx,line) in enumerate(mesh.lines)
         points[:,N_points+idx]=sum(mesh.points[:,line],dims=2)/2
      end
      write_data_array(fid,Dict("coordinates"=>points),binary,compressed)
      write(fid,"\t\t\t</Points>\n")
      write(fid,"\t\t\t<Cells>\n")
      write(fid,"\t\t\t\t<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">")
      Npoints=size(mesh.points,2)
      for tet in mesh.tetrahedra
         for pnt in tet
            write(fid,string(Int32(pnt-1))*" ")
         end
         pnt=get_line_idx(mesh,[tet[1]; tet[2]])+Npoints
         write(fid,string(Int32(pnt-1))*" ")
         pnt=get_line_idx(mesh,[tet[2]; tet[3]])+Npoints
         write(fid,string(Int32(pnt-1))*" ")
         pnt=get_line_idx(mesh,[tet[3]; tet[1]])+Npoints
         write(fid,string(Int32(pnt-1))*" ")
         pnt=get_line_idx(mesh,[tet[1]; tet[4]])+Npoints
         write(fid,string(Int32(pnt-1))*" ")
         pnt=get_line_idx(mesh,[tet[2]; tet[4]])+Npoints
         write(fid,string(Int32(pnt-1))*" ")
         pnt=get_line_idx(mesh,[tet[3]; tet[4]])+Npoints
         write(fid,string(Int32(pnt-1))*" ")
      end
      write(fid,"</DataArray>\n")
      write(fid,"\t\t\t\t<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">")
      off=0
      for tet in mesh.tetrahedra
            off+=10#length(tet)
            write(fid,string(off)*" ")
      end
      write(fid,"</DataArray>\n")
      write(fid,"\t\t\t\t<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">")
      for tet in mesh.tetrahedra
         #for pnt in tet
            write(fid,"24 ")
         #end
      end
      write(fid,"</DataArray>\n")
      write(fid,"\t\t\t</Cells>\n")

      write(fid,"\t\t\t<PointData>\n")
      write_data_array(fid,data,binary,compressed)
      write(fid,"\t\t\t</PointData>\n")
      write(fid,"\t\t</Piece>\n")
      write(fid,"\t</UnstructuredGrid>\n")
      write(fid,"</VTKFile>\n")
   end
   return nothing
end
##
function vtk_write_tri(file_name::String,mesh::Mesh,data=[],binary=true,compressed=false)
   open(file_name*".vtu","w") do fid
      if compressed
         compressor=" compressor=\"vtkZLibDataCompressor\""
      else
         compressor=""
      end
      write(fid,"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\" header_type=\"UInt64\"$compressor>\n")
      write(fid,"\t<UnstructuredGrid>\n")
      write(fid,"\t\t<Piece NumberOfPoints=\"$(size(mesh.points,2))\" NumberOfCells=\"$(length(mesh.triangles))\">\n")
      write(fid,"\t\t\t<Points>\n")
      # write(fid,"\t\t\t\t<DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">")
      # for pnt in mesh.points
      #    write(fid,string(pnt)*" ")
      # end
      # write(fid,"</DataArray>")
      write_data_array(fid,Dict("Coordinates"=>mesh.points),binary,compressed)
      write(fid,"\t\t\t</Points>\n")
      write(fid,"\t\t\t<Cells>\n")
      write(fid,"\t\t\t\t<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">")
      for tri in mesh.triangles
         for pnt in tri
            write(fid,string(Int32(pnt-1))*" ")
         end
      end
      write(fid,"</DataArray>\n")
      write(fid,"\t\t\t\t<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">")
      off=0
      for tri in mesh.triangles
            off+=length(tri)
            write(fid,string(off)*" ")
      end
      write(fid,"</DataArray>\n")
      write(fid,"\t\t\t\t<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">")
      for tri in mesh.triangles
         for pnt in tri
            write(fid,"5 ")
         end
      end
      write(fid,"</DataArray>\n")
      write(fid,"\t\t\t</Cells>\n")

      write(fid,"\t\t\t<CellData>\n")
      write_data_array(fid,data,binary,compressed)
      write(fid,"\t\t\t</CellData>\n")
      write(fid,"\t\t</Piece>\n")
      write(fid,"\t</UnstructuredGrid>\n")
      write(fid,"</VTKFile>\n")
   end
   return nothing
end

##

"""
    vtk_write(file_name, mesh, data)

Write vtk-file containing the datavalues `data` given on the usntructured grid `mesh`.

# Arguments
- `file_name::String`: name given to the written files. The name will be preceeded by a specifier. See Notes below.
- `mesh::Mesh`: mesh associated with the data.
- `data::Dict`: Dictionary containing the data.

# Notes
The routine automatically sorts the data according to its type and writes it in up to three diffrent files.
Data that is constant on a tetrahedron goes into `"\$(filename)_const.vtu"`,
data that is linearly interpolated on a tetrahedron goes into `"\$(filename)_lin.vtu"`,
and data that is quadratically interpolated on a tetrahedron goes into `"\$(filename)_quad.vtu"`
"""
function vtk_write(file_name::String,mesh::Mesh,data::Dict,binary=false,compressed=false) #TODO: write into one file
   data0=Dict()
   data1=Dict()
   data2=Dict()
   data_tri=Dict()
   for (key,val) in data
      if maximum(size(val))==length(mesh.tetrahedra)
         data0[key]=val
      elseif maximum(size(val))==size(mesh.points,2) #TODO: check for rare case of same number of points and tets
         data1[key]=val
      elseif maximum(size(val))==size(mesh.points,2)+length(mesh.lines)
         data2[key]=val
      elseif maximum(size(val))==length(mesh.triangles)
         data_tri[key]=val
      else
         println("WARNING: data length of '$key' does not fit to any mesh dimension.")
      end
   end

   if length(data0)!=0
      vtk_write0(file_name*"_const",mesh,data0,binary,compressed)
   end
   if length(data1)!=0
      vtk_write1(file_name*"_lin",mesh,data1,binary,compressed)
   end
   if length(data2)!=0
      vtk_write2(file_name*"_quad",mesh,data2,binary,compressed)
   end
   if length(data_tri)!=0
      vtk_write_tri(file_name*"_tri",mesh,data_tri,binary,compressed)
   end
   return
end
## bloch utilities
# function bloch_expand(mesh::Mesh,sol,b=:b)
#     naxis=mesh.dos.naxis
#     nxsector=mesh.dos.nxsector
#     DOS=mesh.dos.DOS
#     v=zeros(ComplexF64,naxis+nxsector*DOS)
#     v[1:naxis]=sol.v[1:naxis]
#     B=sol.params[b]
#     for s=0:DOS-1
#         v[naxis+1+s*nxsector:naxis+(s+1)*nxsector]=sol.v[naxis+1:naxis+nxsector].*exp(+2.0im*pi/DOS*B*s)
#     end
#     return v
# end
# function bloch_expand(mesh::Mesh,vec::Array,B::Real=0)
#     naxis=mesh.dos.naxis
#     nxsector=mesh.dos.nxsector
#     DOS=mesh.dos.DOS
#     v=zeros(ComplexF64,naxis+nxsector*DOS)
#     v[1:naxis]=sol.v[1:naxis]
#     for s=0:DOS-1
#         v[naxis+1+s*nxsector:naxis+(s+1)*nxsector]=vec[naxis+1:naxis+nxsector].*exp(+2.0im*pi/DOS*B*s)
#     end
#     return v
# end

## helper functions for writing vtk data

function write_data_array(fid,data,binary=true,compressed=false)
   if binary
      if compressed
         iob = IOBuffer()
         #io.append=true
         #iob = Base64.Base64EncodePipe(iob);
      else
         iob = Base64.Base64EncodePipe(fid);
      end
      format="\"binary\""
   else
      format="\"ascii\""
   end

   for dataname in keys(data)
      if size(data[dataname],2)>1 #write vector output #TODO: Sanity checks for three components
         datatype="$(typeof(data[dataname][1]))"
         bytesize=sizeof(data[dataname][1])
         write(fid,"\t\t\t\t<DataArray type=\"$datatype\" NumberOfComponents=\"3\" Name=\"$dataname\" format=$format>")
         write(fid,"\n\t\t\t\t\t")
         if binary
            len=length(data[dataname])
            if !compressed
               num_of_bytes=UInt64(length(data[dataname])*bytesize)
               write(iob, num_of_bytes)
            end
            for vec in data[dataname]
               for num in vec
                  write(iob, num)#write number to base64 coded stream
               end
            end
         else
            for vec in data[dataname]
               write(fid,string(vec)*" ")
            end
         end
      else
         datatype="$(typeof(data[dataname][1]))"
         bytesize=sizeof(data[dataname][1])
         #datatype="Float64"
         write(fid,"\t\t\t\t<DataArray type=\"$datatype\" Name=\"$dataname\" format=$format>")
         write(fid,"\n\t\t\t\t\t")
         if binary
               #num_of_bytes=UInt64(length(data[dataname])*bytesize)
            if !compressed
               num_of_bytes=UInt64(length(data[dataname])*bytesize)
               write(iob, num_of_bytes)
            end
            for datum in data[dataname]
               idxx=0
               #println("$idxx : $dataname")
               for (idxc,num) in enumerate(datum)
                  idxx+=1
                  #println("$idxx : $num")
                  write(iob, num)#write number to base64 coded stream
               end
            end
         else
            for datum in data[dataname]
                  write(fid,string(datum)*" ")
            end
         end
      end
      if binary
         #close(iob) #close io stream
         #zip compress data
         if compressed
            block_zcompress!(fid,take!(iob))
         end
         #write(fid,String(take!(io))) #write encoded data to file
      end
      write(fid,"\n\t\t\t\t</DataArray>\n")
   end
   if compressed
      close(iob)
   end
end

function block_zcompress!(fid,data,block_length=2^10)
    iob64_encode = Base64.Base64EncodePipe(fid);
    len=length(data)
    number_of_blocks=len÷block_length
    length_of_last_partial_block=len%block_length
    if length_of_last_partial_block!=0
        number_of_blocks+=1
    end
    compressed_block_size = Array{UInt64}(undef,number_of_blocks)
    compressed_data=[]
    for idx=1:number_of_blocks
        if idx==number_of_blocks
            block=data[((idx-1)*block_length+1):end]
        else
            block=data[((idx-1)*block_length+1):(idx*block_length)]
        end
        compressed_block = CodecZlib.transcode(CodecZlib.GzipCompressor, block)
        compressed_block_size[idx]=sizeof(compressed_block)
        append!(compressed_data,[compressed_block])
    end

    #write compressed binary entry
    #io = IOBuffer()
    #header
    write(iob64_encode, UInt64(number_of_blocks))
    write(iob64_encode, UInt64(block_length))
    write(iob64_encode, UInt64(length_of_last_partial_block))
    for idx=1:number_of_blocks
      write(iob64_encode, UInt64(compressed_block_size[idx]))
    end
    #data
    for idx=1:number_of_blocks
        for num in compressed_data[idx]
            write(iob64_encode, num)
        end
    end

    return nothing
end


text:
# This file was generated, do not modify it. # hide
bar(v) = v[end-1:end]
combine(df, :a => foo, :b => bar)


text:
abstract type AbstractSymTenArray{T,N} <: AbstractArray{SymTen{T},N} end

Base.size(v::AbstractSymTenArray) =
    size(xxvec(v))

@inline Base.@propagate_inbounds function Base.getindex(v::AbstractSymTenArray{T,N},i::Int) where {T,N}
    xx = xxvec(v)
    xy = xyvec(v)
    xz = xzvec(v)
    yy = yyvec(v)
    yz = yzvec(v)
    zz = zzvec(v)
    xxv = xx[i]
    xyv = xy[i]
    xzv = xz[i]
    yyv = yy[i]
    yzv = yz[i]
    zzv = zz[i]
    return SymTen{T}(xxv,xyv,xzv,yyv,yzv,zzv)
end

@inline Base.@propagate_inbounds function Base.getindex(v::AbstractSymTenArray{T,N},I::Vararg{Int,N}) where {T,N}
    xx = xxvec(v)
    xy = xyvec(v)
    xz = xzvec(v)
    yy = yyvec(v)
    yz = yzvec(v)
    zz = zzvec(v)
    xxv = xx[I...]
    xyv = xy[I...]
    xzv = xz[I...]
    yyv = yy[I...]
    yzv = yz[I...]
    zzv = zz[I...]
    return SymTen{T}(xxv,xyv,xzv,yyv,yzv,zzv)
end

@inline Base.@propagate_inbounds function Base.setindex!(v::AbstractSymTenArray,ten,i::Int)
    xx = xxvec(v)
    xy = xyvec(v)
    xz = xzvec(v)
    yy = yyvec(v)
    yz = yzvec(v)
    zz = zzvec(v)
    setindex!(xx,ten.xx,i)
    setindex!(xy,ten.xy,i)
    setindex!(xz,ten.xz,i)
    setindex!(yy,ten.yy,i)
    setindex!(yz,ten.yz,i)
    setindex!(zz,ten.zz,i)
    return v
end

@inline Base.@propagate_inbounds function Base.setindex!(v::AbstractSymTenArray{T,N},ten,I::Vararg{Int,N}) where {T,N}
    xx = xxvec(v)
    xy = xyvec(v)
    xz = xzvec(v)
    yy = yyvec(v)
    yz = yzvec(v)
    zz = zzvec(v)
    setindex!(xx,ten.xx,I...)
    setindex!(xy,ten.xy,I...)
    setindex!(xz,ten.xz,I...)
    setindex!(yy,ten.yy,I...)
    setindex!(yz,ten.yz,I...)
    setindex!(zz,ten.zz,I...)
    return v
end

struct SymTenArray{T,N,A<:AbstractArray{T,N},B<:AbstractArray{T,N},C<:AbstractArray{T,N},D<:AbstractArray{T,N},E<:AbstractArray{T,N},F<:AbstractArray{T,N}} <: AbstractSymTenArray{T,N}
    xx::A
    xy::B
    xz::C
    yy::D
    yz::E
    zz::F
end

@generated function Base.IndexStyle(::Type{<:SymTenArray{T,N,A,B,C,D,E,F}}) where {T,N,A,B,C,D,E,F}
    indexstyle = all(x->(IndexStyle(x)===IndexLinear()),(A,B,C,D,E,F)) ? IndexLinear() : IndexCartesian()
    return :($indexstyle)
end

@inline SymTenArray{T}(dims::Vararg{Int,N}) where {T,N} = SymTenArray(zeros(T,dims...),zeros(T,dims...),zeros(T,dims...),zeros(T,dims...),zeros(T,dims...),zeros(T,dims...))
@inline SymTenArray(dims::Vararg{Int,N}) where {N} = SymTenArray{Float64}(dims...)

@inline xxvec(v::SymTenArray) =
    v.xx
@inline xyvec(v::SymTenArray) =
    v.xy
@inline xzvec(v::SymTenArray) =
    v.xz
@inline yyvec(v::SymTenArray) =
    v.yy
@inline yzvec(v::SymTenArray) =
    v.yz
@inline zzvec(v::SymTenArray) =
    v.zz

Base.similar(a::SymTenArray) = SymTenArray(similar(a.xx),similar(a.xy),similar(a.xz),similar(a.yy),similar(a.yz),similar(a.zz))

Base.read!(io::NTuple{6,A},a::SymTenArray) where{A<:Union{<:IO,<:AbstractString}} = (read!(io[1],a.xx); read!(io[2],a.xy); read!(io[3],a.xz); read!(io[4],a.yy); read!(io[5],a.yz); read!(io[6],a.zz))
Base.write(io::NTuple{6,A},a::SymTenArray) where{A<:Union{<:IO,<:AbstractString}} = (write(io[1],a.xx); write(io[2],a.xy); write(io[3],a.xz); write(io[4],a.yy); write(io[5],a.yz); write(io[6],a.zz))

abstract type AbstractSymTrTenArray{T,N} <: AbstractArray{SymTen{T},N} end

@inline Base.@propagate_inbounds function Base.getindex(v::AbstractSymTrTenArray{T,N},i::Int) where {T,N}
    xx = xxvec(v)
    xy = xyvec(v)
    xz = xzvec(v)
    yy = yyvec(v)
    yz = yzvec(v)
    xxv = xx[i]
    xyv = xy[i]
    xzv = xz[i]
    yyv = yy[i]
    yzv = yz[i]
    return SymTen{T}(xxv,xyv,xzv,yyv,yzv,-(xxv+yyv))
end

@inline Base.@propagate_inbounds function Base.getindex(v::AbstractSymTrTenArray{T,N},I::Vararg{Int,N}) where {T,N}
    xx = xxvec(v)
    xy = xyvec(v)
    xz = xzvec(v)
    yy = yyvec(v)
    yz = yzvec(v)
    xxv = xx[I...]
    xyv = xy[I...]
    xzv = xz[I...]
    yyv = yy[I...]
    yzv = yz[I...]
    return SymTen{T}(xxv,xyv,xzv,yyv,yzv,-(xxv+yyv))
end

@inline Base.@propagate_inbounds function Base.setindex!(v::AbstractSymTrTenArray,ten,i::Int)
    xx = xxvec(v)
    xy = xyvec(v)
    xz = xzvec(v)
    yy = yyvec(v)
    yz = yzvec(v)
    ten = traceless(ten)
    setindex!(xx,ten.xx,i)
    setindex!(xy,ten.xy,i)
    setindex!(xz,ten.xz,i)
    setindex!(yy,ten.yy,i)
    setindex!(yz,ten.yz,i)
    return v
end

@inline Base.@propagate_inbounds function Base.setindex!(v::AbstractSymTrTenArray{T,N},vec,I::Vararg{Int,N}) where {T,N}
    xx = xxvec(v)
    xy = xyvec(v)
    xz = xzvec(v)
    yy = yyvec(v)
    yz = yzvec(v)
    ten = traceless(vec)
    setindex!(xx,ten.xx,I...)
    setindex!(xy,ten.xy,I...)
    setindex!(xz,ten.xz,I...)
    setindex!(yy,ten.yy,I...)
    setindex!(yz,ten.yz,I...)
    return v
end

Base.size(v::AbstractSymTrTenArray) =
    size(xxvec(v))

struct SymTrTenArray{T,N,A<:AbstractArray{T,N},B<:AbstractArray{T,N},C<:AbstractArray{T,N},D<:AbstractArray{T,N},E<:AbstractArray{T,N}} <: AbstractSymTrTenArray{T,N}
    xx::A
    xy::B
    xz::C
    yy::D
    yz::E
end

@generated function Base.IndexStyle(::Type{<:SymTrTenArray{T,N,A,B,C,D,E}}) where {T,N,A,B,C,D,E}
    indexstyle = all(x->(IndexStyle(x)===IndexLinear()),(A,B,C,D,E)) ? IndexLinear() : IndexCartesian()
    return :($indexstyle)
end

@inline SymTrTenArray{T}(dims::Vararg{Int,N}) where {T,N} = SymTrTenArray(zeros(T,dims...),zeros(T,dims...),zeros(T,dims...),zeros(T,dims...),zeros(T,dims...))
@inline SymTrTenArray(dims::Vararg{Int,N}) where {N} = SymTrTenArray{Float64}(dims...)

@inline xxvec(v::SymTrTenArray) =
    v.xx
@inline xyvec(v::SymTrTenArray) =
    v.xy
@inline xzvec(v::SymTrTenArray) =
    v.xz
@inline yyvec(v::SymTrTenArray) =
    v.yy
@inline yzvec(v::SymTrTenArray) =
    v.yz

Base.similar(a::SymTrTenArray) = SymTrTenArray(similar(a.xx),similar(a.xy),similar(a.xz),similar(a.yy),similar(a.yz))

Base.read!(io::NTuple{5,A},a::SymTrTenArray) where{A<:Union{<:IO,<:AbstractString}} = (read!(io[1],a.xx); read!(io[2],a.xy); read!(io[3],a.xz); read!(io[4],a.yy); read!(io[5],a.yz))
Base.write(io::NTuple{5,A},a::SymTrTenArray) where{A<:Union{<:IO,<:AbstractString}} = (write(io[1],a.xx); write(io[2],a.xy); write(io[3],a.xz); write(io[4],a.yy); write(io[5],a.yz))

abstract type AbstractAntiSymTenArray{T,N} <: AbstractArray{AntiSymTen{T},N} end

@inline Base.@propagate_inbounds function Base.getindex(v::AbstractAntiSymTenArray{T,N},i::Int) where {T,N}
    xy = xyvec(v)
    xz = xzvec(v)
    yz = yzvec(v)
    xyv = xy[i]
    xzv = xz[i]
    yzv = yz[i]
    return AntiSymTen{T}(xyv,xzv,yzv)
end

@inline Base.@propagate_inbounds function Base.getindex(v::AbstractAntiSymTenArray{T,N},I::Vararg{Int,N}) where {T,N}
    xy = xyvec(v)
    xz = xzvec(v)
    yz = yzvec(v)
    xyv = xy[I...]
    xzv = xz[I...]
    yzv = yz[I...]
    return AntiSymTen{T}(xyv,xzv,yzv)
end

@inline Base.@propagate_inbounds function Base.setindex!(v::AbstractAntiSymTenArray,ten,i::Int)
    xy = xyvec(v)
    xz = xzvec(v)
    yz = yzvec(v)
    setindex!(xy,ten.xy,i)
    setindex!(xz,ten.xz,i)
    setindex!(yz,ten.yz,i)
    return v
end

@inline Base.@propagate_inbounds function Base.setindex!(v::AbstractAntiSymTenArray{T,N},vec,I::Vararg{Int,N}) where {T,N}
    xy = xyvec(v)
    xz = xzvec(v)
    yz = yzvec(v)
    setindex!(xy,ten.xy,I...)
    setindex!(xz,ten.xz,I...)
    setindex!(yz,ten.yz,I...)
    return v
end

Base.size(v::AbstractAntiSymTenArray) =
    size(xyvec(v))

struct AntiSymTenArray{T,N,A<:AbstractArray{T,N},B<:AbstractArray{T,N},C<:AbstractArray{T,N}} <: AbstractAntiSymTenArray{T,N}
    xy::A
    xz::B
    yz::C
end

@generated function Base.IndexStyle(::Type{<:AntiSymTenArray{T,N,A,B,C}}) where {T,N,A,B,C}
    indexstyle = all(x->(IndexStyle(x)===IndexLinear()),(A,B,C)) ? IndexLinear() : IndexCartesian()
    return :($indexstyle)
end

@inline AntiSymTenArray{T}(dims::Vararg{Int,N}) where {T,N} = AntiSymTenArray(zeros(T,dims...),zeros(T,dims...),zeros(T,dims...))
@inline AntiSymTenArray(dims::Vararg{Int,N}) where {N} = AntiSymTenArray{Float64}(dims...)

@inline xyvec(v::AntiSymTenArray) =
    v.xy
@inline xzvec(v::AntiSymTenArray) =
    v.xz
@inline yzvec(v::AntiSymTenArray) =
    v.yz

Base.similar(a::AntiSymTenArray) = AntiSymTenArray(similar(a.xy),similar(a.xz),similar(a.yz))

Base.read!(io::NTuple{3,A},a::AntiSymTenArray) where{A<:Union{<:IO,<:AbstractString}} = (read!(io[1],a.xy); read!(io[2],a.xz); read!(io[3],a.yz))
Base.write(io::NTuple{3,A},a::AntiSymTenArray) where{A<:Union{<:IO,<:AbstractString}} = (write(io[1],a.xy); write(io[2],a.xz); write(io[3],a.yz))


text:
"""
Gene struct.

# Fields
````
id :: String
name :: Union{String, Nothing}
notes :: Dict{String, Vector{String}}
annotation :: Dict{String, Union{Vector{String}, String}}
````
"""
mutable struct Gene
    id::String
    name::Maybe{String}
    notes::Notes
    annotations::Annotations

    Gene(id::String = ""; name = nothing, notes = Notes(), annotations = Annotations()) =
        new(id, name, notes, annotations)
end


text:
export Gcomb, Gls, Gkl, Gls2, Gpref
"""
combine different G's
"""
mutable struct Gcomb
    Gs::Array
end

function (G::Gcomb)(U)
    res = G.Gs[1](U)
    for k=2:length(G.Gs)
        res += G.Gs[k](U)
    end
    return res
end
function Base.show(io::IO, G::Gcomb)
  print(io, G.Gs[1])
  for k=2:length(G.Gs)
      print(io, " + $(G.Gs[k])")
  end
end

"""
Least-Squares Terminal Cost
"""
mutable struct Gls
    rho0
    rho1
    rho0x::AbstractVector # = rho0(X0), stored for efficiency
    rho1x::AbstractVector # = rho1(X0), stored for efficiency
    mu::Real      # = penalty parameter
end

function Base.show(io::IO, G::Gls)
  print(io, "$(G.mu) ⋅ Gls(U)")
end

function (G::Gls)(U::AbstractArray{R}) where R <: Real

    d   = size(U,1)-4
    return G.mu*R(0.5)* ( G.rho0x./ exp.(-U[d+1,:]) - G.rho1(U[1:d,:]) ).^2 .* (exp.(-U[d+1,:])./G.rho0x)

end

function getDeltaG(G::Gls,U::AbstractArray)
    (d,nex) = size(U)
    d      -= 4
    detDy = exp.(U[d+1,:])
    return G.mu.*(G.rho0x ./detDy - G.rho1(U[1:d,:]))
end

"""
KL Divergence Terminal Cost
"""
mutable struct Gkl
    rho0
    rho1
    rho0x::Vector # = rho0(X0), stored for efficiency
    rho1x::Vector # = rho1(X0), stored for efficiency
    mu::Real      # = penalty parameter
end
function (G::Gkl)(U)
    (d,nex) = size(U)
    d -= 4
    return G.mu .* (log.(G.rho0x) - U[d+1,:] - log.(G.rho1(U[1:d,:])))
end
function Base.show(io::IO, G::Gkl)
  print(io, "$(G.mu) ⋅ Gkl(U)")
end

function getDeltaG(G::Gkl,U::AbstractArray{R})  where R <: Real
    (d,nex) = size(U)
    d      -= 4
    return G.mu.*(R(1.0) .+ log.(G.rho0x) - U[d+1,:] - log.(G.rho1(U[1:d,:])))
end

"""
 Preference Terminal Cost
"""
mutable struct Gpref
    Pref::Function  # preference function
    rho0x::Vector   #
    mu::Real        # = penalty parameter
end
function (G::Gpref)(U)
    (d,nex) = size(U)
    d -= 4
    return G.mu .* G.Pref(U[1:d,:])
end

function getDeltaG(G::Gpref,U::AbstractArray{R})  where R <: Real
    (d,nex) = size(U)
    d      -= 4
    return G.mu.*(G.Pref(U[1:d,:]))
end

max_stars_repo_name: EmoryMLIP/MFGnet

max_issues_repo_name: EmoryMLIP/MFGnet

max_forks_repo_name: EmoryMLIP/MFGnet


text:
cd(@__DIR__); include("setups/grid23x22.jl")
gr(dpi = 200)

##
x = 242

## (a)
l = 11 # vertical
plt = grid_eigenvector_plot(l)
savefig(plt, "../figs/Grid$(Nx)x$(Ny)_eigenvector_l$(l-1).png")

## (d)
plt = grid_NGWFvector_plot(l, x, D, 𝚽)
savefig(plt, "../figs/Grid$(Nx)x$(Ny)_NGWF_DAG_sig02_l$(l-1)_x$(x).png")

## (b)
l = 16  # horizontal
plt = grid_eigenvector_plot(l)
savefig(plt, "../figs/Grid$(Nx)x$(Ny)_eigenvector_l$(l-1).png")

## (e)
plt = grid_NGWFvector_plot(l, x, D, 𝚽)
savefig(plt, "../figs/Grid$(Nx)x$(Ny)_NGWF_DAG_sig02_l$(l-1)_x$(x).png")

## (c)
l = 30  # mix
plt = grid_eigenvector_plot(l)
savefig(plt, "../figs/Grid$(Nx)x$(Ny)_eigenvector_l$(l-1).png")

## (f)
plt = grid_NGWFvector_plot(l, x, D, 𝚽)
savefig(plt, "../figs/Grid$(Nx)x$(Ny)_NGWF_DAG_sig02_l$(l-1)_x$(x).png")


text:
module Monitor
    #
    # just a wee thing to share will web apps
    # to make progress bars, etc.

    using Observables
    using UUIDs

    export Progress

    struct Progress
        uuid   :: UUID
        phase  :: String
        thread :: Int
        count  :: Int
        step   :: Int
        size   :: Int
    end

    # run_obs = Observable( Progress("", 0, 0 ))

end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
function map_to_circuit(a, ::Int, ::Union{ChebyshevSparse, Product})
	return a
end

function map_to_circuit(a, b::Real, ::ChebyshevTower)
	return a * b
end

function calculate_diff_evalue(feature_map_circuit::AbstractBlock, cost::AbstractBlock, var::AbstractBlock, N::Int, theta, x, pc; mapping::FeatureMap)
	diff_evalue = 0.0
	for j in 1:N
		diff_evalue_pos = real.(expect(cost, zero_state(N) => new_circuit(feature_map_circuit, var, x, theta, N, pc, j, π/2; mapping)))

		diff_evalue_neg = real.(expect(cost, zero_state(N) => new_circuit(feature_map_circuit, var, x, theta, N, pc, j, -π/2; mapping)))

		diff_evalue += (diff_evalue_pos - diff_evalue_neg) * map_to_circuit(gradient(_x -> phi(_x, pc, mapping), x)[1], j, mapping)
	end
	diff_evalue /= 2
	return diff_evalue
end


text:
function mergesort(v)
    if length(v) <= 1
        return v
    end
    return merge(mergesort(v[1:end÷2]), mergesort(v[end÷2+1:end]))
end

function merge(v1, v2)
    if isempty(v1)
        return v2

    elseif isempty(v2)
        return v1

    elseif v2[1] < v1[1]
        return [v2[1]; merge(v1, v2[2:end])]

    else # v1[1] <= v2[1]
        return [v1[1]; merge(v1[2:end], v2)]
    end
end

v = [144, 89, 4, 9, 95, 12, 86, 25]
println("unsorted array = ", v)

sorted_v = mergesort(v)

println("merge sorted array = ", sorted_v)

max_stars_repo_name: AllAlgorithms/julia

max_issues_repo_name: AllAlgorithms/julia

max_forks_repo_name: AllAlgorithms/julia


text:
using Revise
using StatsBase
using Statistics
using OnlineStats
using BenchmarkTools
using WeightedOnlineStats

x = (collect(reverse((1:10_000_000)*0.0000001)))
w = (collect((1:10_000_000)*0.0000001))
x1 = x[1:5_000_000]
x2 = (x[5_000_001:end])
w1 = w[1:5_000_000]
w2 = (w[5_000_001:end])

vw = WeightedVariance()
vw1 = WeightedVariance()
vw2 = WeightedVariance()

fit!(vw, x, w)
fit!(vw1, x1, w1)
fit!(vw2, x2, w2)

@code_warntype merge!(vw1, vw2)

@benchmark merge($vw1, $vw2)
#=
BenchmarkTools.Trial.1.1:
ll. 133-145 (w/o 135-138)
  memory estimate:  544 bytes
  allocs estimate:  7
  --------------
  minimum time:     3.140 μs (0.00% GC)
  median time:      3.257 μs (0.00% GC)
  mean time:        3.476 μs (1.52% GC)
  maximum time:     272.601 μs (97.37% GC)
  --------------
  samples:          10000
  evals/sample:     8

BenchmarkTools.Trial.1.2:
ll. 133-145 (w/o 139-141)
  memory estimate:  544 bytes
  allocs estimate:  7
  --------------
  minimum time:     3.141 μs (0.00% GC)
  median time:      3.206 μs (0.00% GC)
  mean time:        3.380 μs (1.54% GC)
  maximum time:     268.521 μs (97.16% GC)
  --------------
  samples:          10000
  evals/sample:     8

BenchmarkTools.Trial.2:
ll. 151-161
  memory estimate:  544 bytes
  allocs estimate:  7
  --------------
  minimum time:     3.136 μs (0.00% GC)
  median time:      3.229 μs (0.00% GC)
  mean time:        3.347 μs (1.53% GC)
  maximum time:     262.741 μs (97.26% GC)
  --------------
  samples:          10000
  evals/sample:     8
=#
x = (collect(reverse((1:10_000_000)*0.0000001)))
w = (collect((1:10_000_000)*0.0000001))
x1 = x[1:5_000_000]
x2 = BigFloat.(x[5_000_001:end])
w1 = w[1:5_000_000]
w2 = BigFloat.(w[5_000_001:end])

vw = WeightedVariance()
vw1 = WeightedVariance()
vw2 = WeightedVariance(BigFloat)

fit!(vw, x, w)
fit!(vw1, x1, w1)
fit!(vw2, x2, w2)

@code_warntype merge(vw1, vw2)
merge(vw1, vw2)
merge(vw2, vw1)
@benchmark merge($vw1, $vw2)
@benchmark merge($vw2, $vw1)
#=
BenchmarkTools.Trial.1.1:
ll. 133-145 (w/o 135-138)
  memory estimate:  2.28 KiB
  allocs estimate:  39
  --------------
  minimum time:     4.955 μs (0.00% GC)
  median time:      5.060 μs (0.00% GC)
  mean time:        5.603 μs (3.95% GC)
  maximum time:     346.025 μs (94.78% GC)
  --------------
  samples:          10000
  evals/sample:     7

BenchmarkTools.Trial.1.2:
ll. 133-145 (w/o 139-141)
  memory estimate:  2.06 KiB
  allocs estimate:  35
  --------------
  minimum time:     4.779 μs (0.00% GC)
  median time:      4.861 μs (0.00% GC)
  mean time:        5.366 μs (3.63% GC)
  maximum time:     348.297 μs (96.43% GC)
  --------------
  samples:          10000
  evals/sample:     7

BenchmarkTools.Trial.2:
ll. 151-161
  memory estimate:  2.50 KiB
  allocs estimate:  43
  --------------
  minimum time:     5.398 μs (0.00% GC)
  median time:      5.515 μs (0.00% GC)
  mean time:        6.042 μs (3.80% GC)
  maximum time:     426.032 μs (96.56% GC)
  --------------
  samples:          10000
  evals/sample:     6
=#

using Traceur

x = rand(100)
x2 = rand(100, 2)
w = rand(100)

@trace fit!(WeightedSum(), x, w)
@trace fit!(WeightedSum(Float32), x, w)

@trace fit!(WeightedMean(), x, w)
@trace fit!(WeightedMean(Float32), x, w)

@trace fit!(WeightedVariance(), x, w)
@trace fit!(WeightedVariance(Float32), x, w)

@trace fit!(WeightedCovMatrix(), x2, w)
@trace fit!(WeightedCovMatrix(Float32), x2, w)


text:
module NonlinearSolve

using Reexport
using UnPack: @unpack
using FiniteDiff, ForwardDiff
using ForwardDiff: Dual
using Setfield
using StaticArrays
using RecursiveArrayTools
using LinearAlgebra
import ArrayInterface
import IterativeSolvers
import RecursiveFactorization

@reexport using SciMLBase

abstract type AbstractNonlinearProblem{uType,isinplace} end
abstract type AbstractNonlinearSolveAlgorithm end
abstract type AbstractBracketingAlgorithm <: AbstractNonlinearSolveAlgorithm end
abstract type AbstractNewtonAlgorithm{CS,AD} <: AbstractNonlinearSolveAlgorithm end
abstract type AbstractNonlinearSolver end
abstract type AbstractImmutableNonlinearSolver <: AbstractNonlinearSolver end

include("utils.jl")
include("jacobian.jl")
include("types.jl")
include("solve.jl")
include("bisection.jl")
include("falsi.jl")
include("raphson.jl")
include("scalar.jl")

# DiffEq styled algorithms
export Bisection, Falsi, NewtonRaphson

end # module


text:
using POMDPs, StatsBase, POMDPToolbox, QMDP, RLInterface, AutomotiveDrivingModels, AutoViz, SARSOP, Images
using ArgParse
s = ArgParseSettings()
@add_arg_table s begin
    "--port"
        help = "set port for zmq communication"
        arg_type = Int
        default=9393
end
parsed_args = parse_args(ARGS, s)

include("occluded_crosswalk_env.jl")
include("pomdp_types.jl")
include("spaces.jl")
include("transition.jl")
include("observation.jl")
include("belief.jl")
include("decomposition.jl")
include("adm_helpers.jl")
include("render_helpers.jl")

function add_noise(o::Vector{Float64})
    return o + 1e-3*randn()
end

function POMDPs.initialobs(pomdp::SingleOCPOMDP, s::SingleOCState, rng::AbstractRNG)
    o = gen(DDNNode(:o), pomdp, s, SingleOCAction(0.), s, rng)
    return o
end

pomdp = SingleOCPOMDP()
pomdp.p_birth = 0.3
pomdp.pos_res = 0.5
pomdp.vel_res = 0.5
pomdp.pos_obs_noise = 0.5
pomdp.vel_obs_noise = 0.5

ip = "127.0.0.1"
port = parsed_args["port"]

RLInterface.run_env_server(ip, port, pomdp)

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
using GeoPhyInv
using SparseArrays
using StatsBase
using LinearAlgebra
using Random
using LinearAlgebra
using Test
using ForwardDiff
using Calculus

nx=21
nz=21
nt=4
nznx=nz*nx
mgrid=[range(-div(nz,2), step=1.0, length=nz), range(-div(nx,2), step=1.0, length=nx)]
tgrid=range(0.0,step=0.5, length=nt)
@info "Grids are all set."

Qv=abs.(randn(nz,nx))
η=abs.(randn(nz,nx))
k=abs.(randn(nz,nx))
σ=abs.(randn(nz,nx))
p=randn(nz,nx,nt)
@info "Medium parameters allocated."

ageom=AGeom(mgrid, SSrcs(1), Srcs(1), Recs(30))
update!(ageom, SSrcs(), [0,0], 5, [0,2π])
update!(ageom, Recs(), [0,0], 5, [0,2π])
ACQ=SparseMatrixCSC(ageom[1],mgrid);
@info "ACQ will be used to project ψ onto receivers."

paE=PoissonExpt(p, tgrid, mgrid, Qv, k, η, σ, ACQ)
GeoPhyInv.mod!(paE)

data=paE[:data]
@info string("The dimensions of data are (nt,nr)=",size(data))


text:
#test new formulation of algorithm, for adaptive sde solving

using ProjectedGPE, DifferentialEquations

using PyPlot, Interact

#import ProjectedGPE.Lgp!
#definition of the C-region
M = 30

#system params
g = 0.1
ω = 1.0
γ = 0.1
t0 = 1/ω
μ = 25.0

#Time grid
ti = 0.              #initial time
tf = 80*t0         #final time
Nt = 50             #size of time vector
t  = collect(linspace(ti,tf,Nt))

dt = 0.01*t0    #initial step size [ - should have dt ≪ 2π/μ]

function makexktrans2(basis,M,ω=1.0)
#n-field to x-space
#x,wx,Tx = nfieldtrans(basis,M,nx,ω)
#2M rule but with no variable changes - suitable for doing 2field integrals on a 4 field grid
#actually any linear operator may be computed spectrally this way
x,wx = gausshermite(2M) #rule order set by terms in eq
wx = wx.*exp.(x.^2)/√ω
Tx = eigmat("Hermite",M,x/√ω,ω) #take grid to spectral
#n-field to k-space
k,wk = gausshermite(2M) #set rule order by terms in eq
wk = wk.*exp.(k.^2)*√ω
Tk = eigmat("Hermite",M,k*√ω,1/ω)
Tk = Tk*diagm((-im).^(0:M-1)) #take grid to spectral
#Transform from x --> k via auxiliary states.
Txk = Tk*Tx'
    return x,wx,Tx,k,wk,Tk,Txk
end

x,wx,Tx,k,wk,Tk,Txk = makexktrans2("Hermite",M)

#M modes on 4-field grid, gives 2M points with correct weight for 4-field product
x4,wx4,Tx4 = nfieldtrans("Hermite",M,4)

#create transform from 2-field grid (at 2M-rule precision) to 4-field grid
Tx24 = Tx4*Tx'

#PGPE kinetic term
#out of place
function kinetic(ψ)
    return Txk'*(0.5*wk.*k.^2.*(Txk*(wx.*ψ)))
end

#in place (mutating)
function kinetic!(ψ,dψ)
    dψ .= Txk'*(0.5*wk.*k.^2.*(Txk*(wx.*ψ)))
end

#initialize fields
#ψ0=randn(size(x))+im*randn(size(x)) #create new random initial state
c0=zeros(M)+im*zeros(M) #zero pad to length 2M
c0[25]=sqrt(100) # initial state
ψ0=Tx*c0
ϕ = Tx24*(wx.*ψ0)

function nlin(ψ)
    #Tx*(P.*Tx'*(wx.*abs.(ψ).^2.*ψ))  #right weight??
    ϕ = Tx24*(wx.*ψ)
    return Tx24'*(wx4.*abs.(ϕ).^2.*ϕ)
end

#Time evolution

#dPGPE in reservoir frame (μ is energy zero)

#out of place
function Lgp(t,ψ)
    #return -im*(1-im*γ)*((g*abs.(ψ).^2 .+ 0.25.*x.^2.*exp(-.5*x.^2) .- μ).*ψ .+ kinetic(ψ))
    ψ = -im*(1-im*γ)*((g*abs.(ψ).^2 .- μ .+ 0.5*(x/sqrt(2)).^2).*ψ .+ kinetic(ψ)) #4 field to 2 field
    #ψ = Tx*(Tx'*(wx.*ψ)) #project
    #ψ = wx.*ψ #project
end

#in place
function Lgp2!(t,ψ,dψ)
    kinetic!(ψ,dψ)
    #dψ[:] = -im*(1-im*γ)*((g*abs.(ψ).^2  .+ 0.5.*x.^2 .- μ).*ψ .+ dψ) #2 field to 2 field (non-linear works...?)
    #maybe we just need a function for nonlinear term specific to grids (aswell as potential)
    #result: at least nonlin will be a quadrature
    #get all quaratures right:
    #4 field n for P.E., K.E.
    #same trans, different weight for nlin (?doesn't work if make different):
    dψ .= -im*(1-im*γ)*((0.5.*x.^2 .- μ).*ψ .+ g*nlin(ψ) .+ dψ)
end

Lgp2!(0.,ψ0,ψ0)

tspan = (t[1],t[end])
#prob = ODEProblem(Lgp!,ψ0,tspan,callback=CallbackSet(),mass_matrix=I)
prob = ODEProblem(Lgp2!,ψ0,tspan)
alg = Tsit5()
#alg = DP5()
#abstol = 1e-12
#reltol = 1e-12
println("Started ...")
@time sol = solve(prob,alg,dt=dt,saveat=t,alg_hints=[:stiff],save_everystep=false,dense=false);
println("... Finished.")

## Transform to cartesian grid for plotting

R = sqrt(2μ)
xMax = 2*R
Nx = 200
xp = collect(linspace(-xMax,xMax,Nx))
Txp = eigmat("Hermite",M,xp)
Ntot=zeros(length(t))
for i=1:length(t)
    Ntot[i]=sum(abs.(Tx'*(wx.*sol[i])).^2)
end
figure(figsize=(6,2))
plot(t,Ntot./Ntot[1]);#ylim(0,1.2)

#Plot
f=figure(figsize=(9,2))
@manipulate for i=1:length(t) withfig(f,clear=true) do
        ψ = Txp*Tx'*(wx.*sol[i]);
        θ = unwrap(angle(ψ));
        subplot(1,2,1);plot(xp/R,g*abs.(ψ).^2,xp/R,ones(xp)*μ,"g:");
        xlim(-xMax/R,xMax/R);#ylim(-1,1.5*μ);
        xlabel(L"x/R");ylabel(L"g|ψ|^2")
        subplot(1,2,2);plot(xp/R,θ);
        xlim(-xMax/R,xMax/R);#ylim(-pi*1.1,pi*1.1);
        xlabel(L"x/R");ylabel(L"angle(ψ)");
    end
end


text:
##
# sparse sampling solvers

function estimate_v(opt::AbstractVSSOptions, m::Union{MDP,POMDP}, belief::AbstractVector, depth::Int, rng::AbstractRNG=Random.GLOBAL_RNG)
    if depth >= maxdepth(opt)
        return 0.0
    end

    a_list = Vector{Float64}[]
    q_list = Float64[]
    a_space = actions(m)

    for i in 1:max(floor(action_width(opt) * opt.action_width_decay^depth), 1)
        a = voo_sample(opt.voo, a_list, q_list, a_space)
        push!(a_list, a)
        push!(q_list, estimate_q(opt, m, belief, a, depth))
        if depth + 1 >= maxdepth(opt) && opt.last_action_null
            break
        end
    end
    return maximum(q_list)
end

# VOSS solver - Not tested
function estimate_q(opt::VOSSOptions, m::MDP, belief::AbstractVector, a, depth::Int, rng::AbstractRNG=Random.GLOBAL_RNG)
    qsum = 0.0
    children = Dict{obstype(m), Vector{statetype(m)}}()

    for i in 1:width(opt)
        s = belief[i]
        if !isterminal(m, s)
            sp, o, r = gen(m, s, a, rng)
            vp = estimate_v(opt, m, sp, depth+1)
            qsum += r + discount(m) * vp
        end
    end
    return qsum / width(opt)
end

function valuepairs(p::VSSPlanner{M}, b) where M <: POMDP
    belief = collect(rand(p.rng, b) for i in 1:p.opt.width)
    a_list = Vector{Float64}[]
    q_list = Float64[]
    a_space = actions(p.m)

    for i in 1:action_width(p.opt)
        a = voo_sample(p.opt.voo, a_list, q_list, a_space)
        push!(a_list, a)
        push!(q_list, estimate_q(p.opt, p.m, belief, a, 0))
    end
    return (a_list[i]=>q_list[i] for i in 1:action_width(p.opt))
end

function POMDPs.action(p::VSSPlanner{M}, b) where M <:POMDP
    avps = collect(valuepairs(p, b))
    best = avps[1]
    for av in avps[2:end]
        if last(av) > last(best)
            best = av
        end
    end
    return first(best)
end

# VOWSS solver
function estimate_q(opt::VOWSSOptions, m::POMDP, belief::AbstractVector, a, depth::Int, rng::AbstractRNG=Random.GLOBAL_RNG)
    q = 0.0

    predictions = Vector{statetype(m)}(undef, opt.width)
    observations = Vector{obstype(m)}(undef, opt.width)

    allterminal = true
    wsum = 0.0
    for i in 1:width(opt)
        s, w = weighted_state(belief, i)
        if !isterminal(m, s)
            allterminal = false
            sp, o, r = gen(m, s, a, rng)
            predictions[i] = sp
            observations[i] = o
            q += w * r
        end
        wsum += w
    end

    if allterminal
        return 0.0
    elseif depth + 1 >= maxdepth(opt)
        return q/wsum
    end

    nextbelief = Vector{Pair{statetype(m), Float64}}(undef, opt.width)

    for i in 1:width(opt) # needs to be a separate for loop because it needs all predictions
        s, ow = weighted_state(belief, i)
        if !isterminal(m, s)
            o = observations[i]
            for j in 1:width(opt)
                s, w = weighted_state(belief, j)
                if isterminal(m, s)
                    nextbelief[j] = s=>0.0
                else
                    sp = predictions[j]
                    nextbelief[j] = sp=>w * obs_weight(m, s, a, sp, o)
                end
            end
            vp = estimate_v(opt, m, nextbelief, depth+1)
            q += ow * discount(m) * vp
        end
    end
    return q/wsum
end

weighted_state(b::AbstractVector, i) = b[i]=>1/length(b)
weighted_state(b::AbstractVector{Pair{S,Float64}}, i) where {S} = b[i]


text:
using DataFrames, CSV

dir = @__DIR__

home_dir = joinpath(split(dir, '/')[1:end-1]...)
enzyme_list = CSV.read("/$home_dir/data/re_enzyme_list.csv", DataFrame)
enzyme_list[!,:enzyme] = convert.(String, enzyme_list[!,:enzyme] )


text:
using Distributions, CausalInference, LightGraphs
using Test, Random
# Example from Pearl, Judea. Causality. Cambridge University Press, 2009.

Random.seed!(100)
N = 2000
d = Normal()
p = 0.01

x = randn(N)
v = x + randn(N)*0.5
w = x + randn(N)*0.5
z = v + w + randn(N)*0.5
s = z + randn(N)*0.5
X = [x v w z s]
C = cor(X)
g = pcalg(size(X,2), gausscitest, (C, N), quantile(Normal(), 1-p/2))

@test collect(LightGraphs.edges(g)) == map(x->Edge(x...), [1 => 2
 1 => 3
 2 => 1
 2 => 4
 3 => 1
 3 => 4
 4 => 5])


text:
# include("../../Potentials.jl/src/Potentials.jl")
# using .Potentials
import PotentialUQ
using LinearAlgebra: I, diagm
using Statistics: mean, std
using Random
println("Starting test of sampler")
Random.seed!(1234);

n = 10
β = zeros(10)
A = 1e-1*randn(n, 10)
A[1:10, 1:10] += I(10)
b = randn(n)

β = A\b
σ = 1e-3 .* ones(n)
logσ = log.(σ)
snap = PotentialUQ.Potentials.SNAP(β, 0.9, 6, 1)
dsnap = PotentialUQ.SNAPQDistribution(snap, A, b, logσ)
# dsnap = PotentialUQ.SNAPDistribution(snap, A, b, diagm(σ))
@time samp = PotentialUQ.Sample(dsnap; num_adapts = 100)
println("Samples ", size(samp))
println("Samples ", samp[1, :, 1])
m = mean(samp[:, :, 1], dims = 1)
s = std(samp[:, :, 1], dims = 1)
println("Mean ", m, "\t")
println("Std ", s)
# show(stdout, "text/plain", dsnap.samples[1])
# println(" ")
println("Solution ", A\b)
println("MAP Solution, Optimization β =  ")
show(stdout, "text/plain", m[1:10])
println("MAP Solution, Optimization σ =  ")
show(stdout, "text/plain", exp.(m[11:end]))


text:
module OptionsPricing

export bs, crr, mc, Option, EuropeanOption, AmericanOption, AsianOption

abstract type Option end

struct EuropeanOption<:Option
    s::Float64
    k::Float64
    r::Float64
    σ::Float64
    t::Float64
    call_or_put::String
end

struct AmericanOption<:Option
    s::Float64
    k::Float64
    r::Float64
    σ::Float64
    t::Float64
    call_or_put::String
end

struct AsianOption<:Option
    s::Float64
    k::Float64
    r::Float64
    σ::Float64
    t::Float64
    call_or_put::String
    dates::Vector{DateTime}
end

include("bs.jl")
include("crr.jl")
include("mc.jl")

end # module


text:
using GaussianProcesses
using GPR
using Optim
using ConstrainedDynamicsVis
using ConstrainedDynamics
using LineSearches
using Statistics

function _experiment_cp_min(config, id; usesin = false, meandynamics = false)
    traindfs, testdfs = config["datasets"]  # Each thread operates on its own dataset -> no races
    traindf = traindfs.df[id][shuffle(1:nrow(traindfs.df[id]))[1:config["nsamples"]], :]
    testdf = testdfs.df[id][shuffle(1:nrow(testdfs.df[id]))[1:config["testsamples"]], :]
    mechanism = cartpole(1, Δt=0.01, threadlock = config["mechanismlock"])[2]  # Reset Δt to 0.01 in mechanism
    l = mechanism.bodies[2].shape.rh[2]

    xtest_future_true = [CState(x) for x in testdf.sfuture]
    # Add noise to the dataset
    for df in [traindf, testdf]
        applynoise!(df, config["Σ"], "CP", config["Δtsim"], l)
    end
    # Create train and testsets
    xtrain_old = [max2mincoordinates(CState(x), mechanism) for x in traindf.sold]
    if usesin
        xtrain_old = reduce(hcat, [[s[1], s[2], sin(s[3]), cos(s[3]), s[4]] for s in xtrain_old])
    else
        xtrain_old = reduce(hcat, xtrain_old)
    end
    xtrain_curr = [max2mincoordinates(CState(x), mechanism) for x in traindf.scurr]
    vωindices = [9, 24]
    ytrain = [[s[id] for s in xtrain_curr] for id in [2,4]]  # v12, ω21
    xtest_old = [max2mincoordinates(CState(x),mechanism) for x in testdf.sold]

    predictedstates = Vector{CState{Float64,2}}()
    params = config["params"]
    if usesin
        params = [params[1], params[2], params[3], params[4], params[4], params[5]]
    end
    gps = Vector{GPE}()

    function xtransform(x, _)
        if usesin
            x, v, sθ, cθ, ω = x
            θ = atan(sθ, cθ)
        else
            x, v, θ, ω = x
        end
        q = UnitQuaternion(RotX(θ))
        x2curr = [0, x, 0] + [0, .5sin(θ)l, -.5cos(θ)l]
        θnext = ω*0.01 + θ
        x2next = [0, x, 0] + [0, v, 0].*0.01 + [0, .5sin(θnext)l, -.5cos(θnext)l]
        v2 = (x2next - x2curr)/0.01
        return [0, x, 0, 1, 0, 0, 0, 0, v, 0, 0, 0, 0,
                x2curr..., q2vec(q)..., v2..., ω, 0, 0]
    end

    cache = MDCache()
    for (id, yi) in enumerate(ytrain)
        kernel = SEArd(log.(params[2:end]), log(params[1]))
        mean = meandynamics ? MeanDynamics(mechanism, getμ(vωindices), id, cache, xtransform=xtransform) : MeanZero()
        gp = GP(xtrain_old, yi, mean, kernel)
        GaussianProcesses.optimize!(gp, LBFGS(linesearch = BackTracking(order=2)), Optim.Options(time_limit=10.))
        push!(gps, gp)
    end

    for i in 1:length(xtest_old)
        predictedstate = predictdynamicsmin(mechanism, "CP", gps, xtest_old[i], config["simsteps"]; usesin=usesin)
        push!(predictedstates, predictedstate)
    end
    return predictedstates, xtest_future_true
end

function experiment_cp_mz_min(config, id)
    return _experiment_cp_min(config, id, usesin=false, meandynamics=false)
end

function experiment_cp_mz_min_sin(config, id)
    return _experiment_cp_min(config, id, usesin=true, meandynamics=false)
end

function experiment_cp_md_min(config, id)
    return _experiment_cp_min(config, id, usesin=false, meandynamics=true)
end

function experiment_cp_md_min_sin(config, id)
    return _experiment_cp_min(config, id, usesin=true, meandynamics=true)
end


text:
using ReactionNetworkImporters, OrdinaryDiffEq, DiffEqBiological,
      Sundials, Plots, DiffEqDevTools, ODEInterface, ODEInterfaceDiffEq,
      LSODA, TimerOutputs

gr()
prnbng = loadrxnetwork(BNGNetwork(), "BNGRepressilator",
                       joinpath(dirname(pathof(ReactionNetworkImporters)),"..","data","bcr","bcr.net"))

rn = deepcopy(prnbng.rn)
addodes!(rn; build_jac=false, build_symfuncs=false, build_paramjac=false)
tf = 100000.0
oprob = ODEProblem(rn, prnbng.u₀, (0.,tf), prnbng.p);

densejac_rn = deepcopy(prnbng.rn)
# zeroout_jac=true is needed to keep the generated expressions from being too big for the compiler
addodes!(densejac_rn; build_jac=true, zeroout_jac = true, sparse_jac = false, build_symfuncs=false, build_paramjac=false)
densejacprob = ODEProblem(densejac_rn, prnbng.u₀, (0.,tf), prnbng.p);

sparsejac_rn = deepcopy(prnbng.rn)
addodes!(sparsejac_rn; build_jac=true, sparse_jac = true, build_symfuncs=false, build_paramjac=false)
sparsejacprob = ODEProblem(sparsejac_rn, prnbng.u₀, (0.,tf), prnbng.p);

@show numspecies(rn) # Number of ODEs
@show numreactions(rn) # Apprx. number of terms in the ODE
@show numparams(rn) # Number of Parameters

const to = TimerOutput()
u₀ = prnbng.u₀
u = copy(u₀);
du = similar(u);
p = prnbng.p
@timeit to "ODERHS Eval1" rn.f(du,u,p,0.)
@timeit to "ODERHS Eval2" rn.f(du,u,p,0.)
sparsejac_rn.f(du,u,p,0.)

J = zeros(length(u),length(u))
@timeit to "DenseJac Eval1" densejac_rn.jac(J,u,p,0.)
@timeit to "DenseJac Eval2" densejac_rn.jac(J,u,p,0.)

Js = similar(sparsejac_rn.odefun.jac_prototype)
@timeit to "SparseJac Eval1" sparsejac_rn.jac(Js,u,p,0.)
@timeit to "SparseJac Eval2" sparsejac_rn.jac(Js,u,p,0.)
show(to)

sol = solve(oprob, CVODE_BDF(), saveat=tf/1000., reltol=1e-5, abstol=1e-5)
plot(sol,legend=false, fmt=:png)

@time sol = solve(oprob,CVODE_BDF(),abstol=1/10^12,reltol=1/10^12)
test_sol = TestSolution(sol)

abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (5:8);
setups = [
          #Dict(:alg=>Rosenbrock23(autodiff=false)),
          Dict(:alg=>TRBDF2(autodiff=false)),
          Dict(:alg=>CVODE_BDF()),
          #Dict(:alg=>rodas()),
          #Dict(:alg=>radau()),
          #Dict(:alg=>Rodas4(autodiff=false)),
          #Dict(:alg=>Rodas5(autodiff=false)),
          Dict(:alg=>KenCarp4(autodiff=false)),
          #Dict(:alg=>RadauIIA5(autodiff=false)),
          #Dict(:alg=>lsoda()),
          ]

wp = WorkPrecisionSet(oprob,abstols,reltols,setups;error_estimate=:l2,
                      saveat=tf/10000.,appxsol=test_sol,maxiters=Int(1e5),numruns=1)
plot(wp)

wp = WorkPrecisionSet(densejacprob,abstols,reltols,setups;error_estimate=:l2,
                      saveat=tf/10000.,appxsol=test_sol,maxiters=Int(1e5),numruns=1)
plot(wp)

setups = [
          #Dict(:alg=>Rosenbrock23(autodiff=false)),
          Dict(:alg=>TRBDF2(autodiff=false)),
          #Dict(:alg=>CVODE_BDF()),
          #Dict(:alg=>rodas()),
          #Dict(:alg=>radau()),
          #Dict(:alg=>Rodas4(autodiff=false)),
          #Dict(:alg=>Rodas5(autodiff=false)),
          Dict(:alg=>KenCarp4(autodiff=false)),
          #Dict(:alg=>RadauIIA5(autodiff=false)),
          #Dict(:alg=>lsoda()),
          ]
wp = WorkPrecisionSet(sparsejacprob,abstols,reltols,setups;error_estimate=:l2,
                      saveat=tf/10000.,appxsol=test_sol,maxiters=Int(1e5),numruns=1)
plot(wp)


text:
mime_types = Dict(
    [
        "abs" => "audio/x-mpeg"
        "ai" => "application/postscript"
        "aif" => "audio/x-aiff"
        "aifc" => "audio/x-aiff"
        "aiff" => "audio/x-aiff"
        "aim" => "application/x-aim"
        "art" => "image/x-jg"
        "asf" => "video/x-ms-asf"
        "asx" => "video/x-ms-asf"
        "au" => "audio/basic"
        "avi" => "video/x-msvideo"
        "avx" => "video/x-rad-screenplay"
        "bcpio" => "application/x-bcpio"
        "bin" => "application/octet-stream"
        "bmp" => "image/bmp"
        "body" => "text/html"
        "cdf" => "application/x-cdf"
        "cer" => "application/x-x509-ca-cert"
        "class" => "application/java"
        "cpio" => "application/x-cpio"
        "csh" => "application/x-csh"
        "css" => "text/css"
        "dib" => "image/bmp"
        "doc" => "application/msword"
        "dtd" => "text/plain"
        "dv" => "video/x-dv"
        "dvi" => "application/x-dvi"
        "eps" => "application/postscript"
        "etx" => "text/x-setext"
        "exe" => "application/octet-stream"
        "gif" => "image/gif"
        "gtar" => "application/x-gtar"
        "gz" => "application/x-gzip"
        "hdf" => "application/x-hdf"
        "hqx" => "application/mac-binhex40"
        "htc" => "text/x-component"
        "htm" => "text/html"
        "html" => "text/html"
        "ief" => "image/ief"
        "jad" => "text/vnd.sun.j2me.app-descriptor"
        "jar" => "application/octet-stream"
        "java" => "text/plain"
        "jnlp" => "application/x-java-jnlp-file"
        "jpe" => "image/jpeg"
        "jpeg" => "image/jpeg"
        "jpg" => "image/jpeg"
        "js" => "text/javascript"
        "kar" => "audio/x-midi"
        "latex" => "application/x-latex"
        "m3u" => "audio/x-mpegurl"
        "mac" => "image/x-macpaint"
        "man" => "application/x-troff-man"
        "me" => "application/x-troff-me"
        "mid" => "audio/x-midi"
        "midi" => "audio/x-midi"
        "mif" => "application/x-mif"
        "mov" => "video/quicktime"
        "movie" => "video/x-sgi-movie"
        "mp1" => "audio/x-mpeg"
        "mp2" => "audio/x-mpeg"
        "mp3" => "audio/x-mpeg"
        "mpa" => "audio/x-mpeg"
        "mpe" => "video/mpeg"
        "mpeg" => "video/mpeg"
        "mpega" => "audio/x-mpeg"
        "mpg" => "video/mpeg"
        "mpv2" => "video/mpeg2"
        "ms" => "application/x-wais-source"
        "nc" => "application/x-netcdf"
        "oda" => "application/oda"
        "pbm" => "image/x-portable-bitmap"
        "pct" => "image/pict"
        "pdf" => "application/pdf"
        "pgm" => "image/x-portable-graymap"
        "pic" => "image/pict"
        "pict" => "image/pict"
        "pls" => "audio/x-scpls"
        "png" => "image/png"
        "pnm" => "image/x-portable-anymap"
        "pnt" => "image/x-macpaint"
        "ppm" => "image/x-portable-pixmap"
        "ps" => "application/postscript"
        "psd" => "image/x-photoshop"
        "qt" => "video/quicktime"
        "qti" => "image/x-quicktime"
        "qtif" => "image/x-quicktime"
        "ras" => "image/x-cmu-raster"
        "rgb" => "image/x-rgb"
        "rm" => "application/vnd.rn-realmedia"
        "roff" => "application/x-troff"
        "rtf" => "application/rtf"
        "rtx" => "text/richtext"
        "sh" => "application/x-sh"
        "shar" => "application/x-shar"
        "smf" => "audio/x-midi"
        "snd" => "audio/basic"
        "src" => "application/x-wais-source"
        "sv4cpio" => "application/x-sv4cpio"
        "sv4crc" => "application/x-sv4crc"
        "swf" => "application/x-shockwave-flash"
        "t" => "application/x-troff"
        "tar" => "application/x-tar"
        "tcl" => "application/x-tcl"
        "tex" => "application/x-tex"
        "texi" => "application/x-texinfo"
        "texinfo" => "application/x-texinfo"
        "tif" => "image/tiff"
        "tiff" => "image/tiff"
        "tr" => "application/x-troff"
        "tsv" => "text/tab-separated-values"
        "txt" => "text/plain"
        "ulw" => "audio/basic"
        "ustar" => "application/x-ustar"
        "xbm" => "image/x-xbitmap"
        "xpm" => "image/x-xpixmap"
        "xwd" => "image/x-xwindowdump"
        "wav" => "audio/x-wav"
        "wbmp" => "image/vnd.wap.wbmp"
        "wml" => "text/vnd.wap.wml"
        "wmlc" => "application/vnd.wap.wmlc"
        "wmls" => "text/vnd.wap.wmlscript"
        "wmlscriptc" => "application/vnd.wap.wmlscriptc"
        "wrl" => "x-world/x-vrml"
        "Z" => "application/x-compress"
        "z" => "application/x-compress"
        "zip" => "application/zip"
    ]
)