text:
const NUM = r"^[0-9]+$"
const NAME = r"^\pL\p{Xan}*$"
const FILE_PATH = [joinpath(PKGDIR, "html")]

function resolve(cmd::JusCmd, vars, str)
    if str == "?"
        ID(cmd)
    elseif (m = match(REF, str)) !== nothing
        vars[m[1]].id
    elseif match(r"^[0-9]+$", str) !== nothing
        parse(Int, str)
    elseif (m = match(PATH_COMPONENT, str)) !== nothing
        ID(m[1] == "@" ? cmd.namespace : m[1], parse(UInt, m[2]))
    else
        Symbol(str)
    end
end

const PATH_METADATA = r"^([^:]*)(?::(.*))?$"

function findvar(cmd::JusCmd, create, vars, path, metadata::Union{Nothing, Dict{Symbol}} = nothing, value = nothing)
    m = match(PATH_METADATA, path)
    if m[2] !== nothing
        path = m[1]
        metadata = parsemetadata(m[2], metadata)
    end
    path = split(path)
    @debug("PATH: $(path)")
    components = map(c-> resolve(cmd, vars, c), path)
    @debug("COMPONENTS: $(components)")
    last = components[end]
    if length(components) > 1 && components[1] isa ID
        parent = components[1]
        components = components[2:end - 1]
    else
        parent = EMPTYID
        components = components[1:end - 1]
    end
    @debug("PARENT COMPONENTS: $(repr(components))")
    @debug("LAST: $(repr(last))")
    for (i, v) in enumerate(components) # path should only be names and numbers
        isemptyid(v) && throw("'?' in the middle of $(path)")
        v isa ID && throw("ID in the middle of $(path)")
        parent == EMPTYID && throw("No parent variable for path $(path)")
        parent = cmd.config[parent][v].id
    end
    metadata === nothing && (metadata = Dict{Symbol, AbstractString}())
    @debug("PARENT: $(repr(parent))")
    @debug("CHECKING LAST...")
    if isemptyid(last)
        !create && throw("'?' without -c")
        parent != EMPTYID && throw("'?' at the end of path: $(path)")
        addvar(cmd, parent, Symbol(""), last, value, metadata), true
    elseif last isa ID
        parent != EMPTYID && throw("ID at the end of path: $(path)")
        cmd.config[last], false
    elseif parent == EMPTYID
        throw("Attempt to get path with no parent, path: $(path)")
    elseif create && last isa Union{Symbol, Number} && !haskey(cmd.config[parent], last)
        addvar(cmd, parent, last, ID(cmd), value, metadata), true
    elseif create
        throw("ERROR, attempt to create a variable that already exists: $(path)")
    else
        @debug("PATH -> $(parent).$(repr(last))")
        @debug("PARENT: $(repr(cmd.config[parent]))")
        cmd.config[parent][last], false
    end
end

function command(cmd::JusCmd{:setmeta})
    var = cmd.config[resolve(cmd, [], cmd.args[1])]
    prop = Symbol(cmd.args[2])
    value = cmd.args[3]
    parentvalue = var.parent == EMPTYID ? nothing : cmd.config[var.parent].value
    set_metadata(VarCommand(cmd, :setmeta, (), var), prop, value)
    route(parentvalue, VarCommand(cmd, :metadata, (prop,), var))
    output(cmd; result = [])
end

function command(cmd::JusCmd{:set})
    local new = []
    local vars = []
    local creating
    local metadata
    local observing
    function newset()
        creating = false
        observing = false
        metadata = Dict{Symbol, AbstractString}()
    end

    @debug("@ SET (FRED 6): $(cmd.args)")
    pos = 1
    newset()
    while pos <= length(cmd.args)
        arg = cmd.args[pos]
        if arg == "-c"
            creating = true
        elseif arg == "-m"
            metadata[cmd.args[pos + 1]] = metadata[pos + 2]
            pos += 2
        else
            value = cmd.args[pos + 1]
            var, creating = findvar(cmd, creating, vars, cmd.args[pos], metadata, value)
            @debug("FOUND VARIABLE: $(var)")
            creating && push!(new, var)
            parentvalue = var.parent == EMPTYID ? nothing : cmd.config[var.parent].value
            try
                pval = var.parent != EMPTYID ? parent_value(cmd.config, var) : nothing
                println("%%%\n%%% SETTING VARIABLE ($(pval)).$(var) FROM $(var.internal_value) TO $(value)")
                route(parentvalue, VarCommand(cmd, :set, (), var; arg = creating ? var.value : value, creating))
            catch err
                if err isa CmdException && err.type == :path
                    @debug "Error refreshing variable $(var)" exception=(err, catch_backtrace())
                    cmd.cancel = true
                else
                    rethrow(err)
                end
            end
            !cmd.cancel && push!(vars, var)
            newset()
            pos += 1
        end
        pos += 1
    end
    output(cmd; result = map(v-> json(cmd, v.id), new))
end

function command(cmd::JusCmd{:get})
    vars = []
    @debug("@ GET ARGS: $(cmd.args)")
    for path in cmd.args
        var, _ = findvar(cmd, false, vars, path)
        push!(vars, var)
    end
    # return a list of [id1, value1, id2, value2, ...]
    [e for v in vars for p in (json(cmd, v.id), json(cmd, v.value)) for e in p]
end

function command(cmd::JusCmd{:observe})
    vars = []
    for path in cmd.args
        var, _ = findvar(cmd, false, vars, path)
        push!(vars, var.id)
    end
    @debug("@ OBSERVE ARGS: $(cmd.args)")
    union!(connection(cmd).observing, vars)
    @debug("OBSERVED VARS: $(repr(map(v-> v.id, allvars(cmd.config, connection(cmd).observing...))))")
    observed = (; ((Symbol(json(id)), cmd.config[id].value) for id in connection(cmd).observing)...)
    @debug("OBSERVED: $(repr(observed))")
    for id in vars
        var = cmd.config[id]
        route(var.value, VarCommand(cmd, :observe, (), var))
    end
    output(cmd,
           result = [],
           update = Dict(json(cmd, vid) => (
               set = json(cmd, cmd.config[vid].value),
               metadata = cmd.config[vid].metadata
           ) for vid in connection(cmd).observing))
end

function output(cmd::JusCmd{NAME}; data...) where NAME
    con = connection(cmd)
    d = haskey(data, :result) ? (; data..., command = NAME) : (; data...)
    con.pending_result = merge(con.pending_result, d)
end

function finish_command(cmd::JusCmd)
    con = connection(cmd)
    println("PENDING RESULT: $(con.pending_result)")
    if haskey(con.pending_result, :result)
        refresh(cmd)
        observe(cmd)
        con.refresh_queued = false
    elseif !con.refresh_queued
        con.refresh_queued = true
        @async begin
            try
                sleep(0.2)
                #this could have been preempted by a command
                if con.refresh_queued
                    refresh(cmd)
                    observe(cmd)
                    con.refresh_queued = false
                end
            catch err
                exit(1)
            end
        end
    end
end

function observe(cmd::JusCmd)
    config = cmd.config
    isempty(config.changes) && return send_output(connection(cmd))
    for (_, connection) in config.connections
        @debug("CHECKING CONNECTION OBSERVING: $(repr(connection.observing))")
        changes = filter(e-> within(config, e[1], connection.observing), config.changes)
        if !isempty(changes)
            fmt = Dict()
            for (id, c) in changes
                var = config[id]
                if haskey(c, :set)
                    c[:set] = json(config, connection, var.value)
                end
                if haskey(c, :metadata)
                    c[:metadata] = Dict(m => var.metadata[m] for m in c[:metadata])
                end
                fmt[json(connection, id)] = Dict{Symbol, Any}(c...)
            end
            output(cmd, update = fmt)
        end
        send_output(connection)
    end
    config.changes = Dict()
end

function send_output(con::Connection)
    if con.pending_result != (;)
        output(con.ws, con.pending_result)
        con.pending_result = (;)
    end
end

function refresh(cmd::JusCmd)
    for (_, v) in cmd.config.vars
        try
            v.parent == EMPTYID && refresh(cmd, v)
        catch err
            @error "Error refreshing" exception=(err, catch_backtrace())
        end
    end
end

function refresh(cmd::JusCmd, var::Var)
    if has_path(var) && var.readable
        parent = parent_value(cmd.config, var)
        if parent !== nothing
            old = var.json_value
            var.refresh_exception = nothing
            try
                route(parent, VarCommand(:get, (); var, config = cmd.config, connection = connection(cmd)))
                old !== var.json_value && changed(cmd.config, var)
            catch err
                var.refresh_exception = err
                var.error_count += 1
                refresh_error(var, err, catch_backtrace())
            end
        end
    end
    for (_, v) in var.namedchildren
        refresh(cmd, v)
    end
    for v in var.indexedchildren
        refresh(cmd, v)
    end
end

"""
    refresh_error

Handle errors during refresh. Log path errors at the debug level.
"""
refresh_error(var::Var, ex::CmdException{:path}, bt) = @debug "Error refreshing variable $(var)" exception=(ex, bt)
refresh_error(var::Var, ex, bt) = @error "Error refreshing variable $(var)" exception=(ex, bt)

function serve(config::Config, ws)
    (; namespace, secret) = input(ws)
    if haskey(config.namespaces, namespace)
        if config.namespaces[namespace].secret !== secret
            @debug("Bad attempt to connect for $(namespace)")
            output(ws, error = "Wrong secret for $(namespace)")
            return
        end
    else
        config.namespaces[namespace] = Namespace(; name=namespace, secret)
    end
    con = Connection(; ws, namespace)
    config.connections[ws] = con
    config.init_connection(con)
    @debug("Connection for $(namespace)")
    while !eof(ws)
        cmd = nothing
        try
            !isopen(ws) && break
            cmd = [input(ws)...]
            jcmd = JusCmd(config, ws, namespace, cmd)
            command(jcmd)
            finish_command(jcmd)
        catch err
            if !(err isa Base.IOError || err isa HTTP.WebSockets.WebSocketError)
                err isa ArgumentError && println(err)
                @error "Error handling comand $(cmd) $(err)" exception=(err, catch_backtrace())
            end
            break
        end
    end
    close(config, con)
    @debug("CLIENT CLOSED: $(repr(namespace))")
end

function close(config::Config, con::Connection)
    Base.close(con.ws)
    for v in con.vars
        delete!(config.vars, v.id)
        delete!(config.changes, v.id)
    end
    delete!(config.connections, con.ws)
end

const FILE_EXT_PAT = r"^(.*)\.([^.]*)$"

const MIME_TYPES_FOR_EXTENSIONS = Dict(
    "js" => "text/javascript",
    "mjs" => "text/javascript",
    "css" => "text/css",
    "html" => "text/html",
    "png" => "image/png",
    "jpg" => "image/jpg",
    "gif" => "image/gif",
)

function mime_type(filename::AbstractString)
    m = match(FILE_EXT_PAT, filename)
    m === nothing && throw("Unknown MIME type for filename $(filename)")
    _, ext = m
    return get(MIME_TYPES_FOR_EXTENSIONS, ext) do
        throw("Unknown MIME type for filename $(filename) with extension $(ext)")
    end
end

add_file_dir(dirname) = push!(FILE_PATH, realpath(dirname))

function find_file(path)
    for dir in FILE_PATH
        dpath = joinpath(dir, path)
        println("TRYING PATH $dpath")
        if isdir(dpath)
            return :directory, dpath
        elseif isfile(dpath)
            return :file, dpath
        end
    end
    return :missing, path
end

"""
    serve_file

taken from [yig's file server](https://gist.github.com/yig/f65e86b7730019d4060449f24342fcb4)
enhanced to handle mime types
"""
function serve_file(req::HTTP.Request)
    ## If it's not a path inside the jail (current working directory), return a 403.
    jail = realpath(FILE_PATH[1])
    # println( "Jail path: ", jail )

    ## Convert the request path to a real path (remove symlinks, remove . and ..)
    ## UPDATE: realpath() assumes that the input points to something that exists.
    ##         We can't use it, since we don't know if the requested path exists.
    ##         We'll use normpath() instead.
    ## Requests start with "/" referring to the server as the root.
    target = replace(req.target, r"^([^?]*)\?.*"=> s"\1")
    @assert length(target) > 0 && target[1] == '/'
    ## Drop the leading "/".
    request_path = normpath(joinpath(jail, target[2:end]))
    ## Convert it to a relative path.
    relative_path = relpath(request_path, jail)
    println("Requested path: ", request_path)
    println("Relative path: ", relative_path)
    ## Return forbidden if the request is above the current working directory.
    length(relative_path) > 0 && splitpath(relative_path)[1] == ".." &&
        @show return HTTP.Response(403)
    type, filepath = find_file(relative_path)
    ## Return not implemented if the request is for a directory.
    if type == :directory
        ## Is there an index.html?
        index = joinpath(filepath, "index.html")
        if isfile(index)
            HTTP.Response(200, ["Content-Type" => mime_type(target)], body = read(index))
        else
            @show HTTP.Response(501)
        end
        ## Return the contents for a file.
    elseif type == :file
        HTTP.Response(200, ["Content-Type" => mime_type(target)], body = read(filepath))
    else
        ## If it's not a file, return a 404.
        @show HTTP.Response(404)
    end
end

function server(config::Config, socket = nothing)
    router = HTTP.Router()
    function handle_request(http)
        if http.message.target == "/ws"
            HTTP.WebSockets.upgrade(http) do ws
                config.serverfunc(config, ws)
            end
        else
            HTTP.handle(router, http)
        end
    end

    println("SERVER ON $(config.host):$(config.port)")
    HTTP.@register(router, "GET", "/ws", r-> serve_websocket(config, r))
    HTTP.@register(router, "GET", "/", serve_file)
    if socket === nothing
        HTTP.listen(handle_request, config.host, config.port)
    else
        HTTP.listen(handle_request; server=socket)
    end
    @debug("HTTP SERVER FINISHED")
end


text:
module TestAqua

import Aqua
import BenchPerf

test() = Aqua.test_all(
    BenchPerf;
    # Compat.jl is conditionally loaded:
    stale_deps = if BenchPerf.Internal.USE_COMPAT
        true
    else
        (; ignore = [:Compat])
    end,
)

end  # module

ext: jl


text:
#=
                            Largest prime factor

    The prime factors of 13195 are 5, 7, 13 and 29.

    What is the largest prime factor of the number 600851475143 ?

                            Answer: 6857
=#

include("../project_euler.jl")
using project_euler

print(pop!(prime_factors(600851475143)), "\n")

max_stars_repo_name: Fazendaaa/project-euler

max_issues_repo_name: Fazendaaa/project-euler

max_forks_repo_name: Fazendaaa/project-euler


text:
"""
    DirichletCategorical(α,n)

A *Dirichlet-categorical distribution* is the compound distribution of the [`Categorical`](@ref) distribution where the probability of each category is distributed according to the [`Dirichlet`](@ref). It has two parameters: the concentration parameters `α` and the optional observation counts `n`.

```math
p(x = i \\vert n, \\alpha) = \\int\\limits_p p(x \\vert p) p(p \\vert n, \\alpha) dp = \\frac{n_i + \\alpha_1}{\\sum\\limits_k n_k + α_k}
```

```julia
DirichletCategorical(α, n)      # DirichcletCategorical distribution with concentration parameters α and observation counts n.

params(d)       # Get the parameters, i.e. (α, n)
update!(d, o)   # Update the counts d.n with observation(s) o.
```

External links:

* [Categorical distribution conjugate prior on Wikipedia](https://en.wikipedia.org/wiki/Categorical_distribution#Bayesian_inference_using_conjugate_prior)
"""
struct DirichletCategorical{T <: Real, U <: Integer} <: DiscreteUnivariateDistribution
    α::Vector{T}
    α0::T
    n::Vector{U}

    function DirichletCategorical(α::Vector{T}, n::Vector{U}) where {T, U}
        α0 = sum(abs, α)
        N = sum(abs, n)
        all(x -> x > 0, α) || throw(ArgumentError("all elemnts in α must be greater than zero"))
        all(x -> x >= 0, n) || throw(ArgumentError("all elements in n must be greater or equal than zero"))
        length(n) == length(α) || throw(ArgumentError("n and α must be of same length"))
        new{T,U}(α, α0, n)
    end
end

DirichletCategorical(α::Vector{<: Integer}, n::Vector{<: Integer})= DirichletCategorical(float(α), n)
DirichletCategorical(k::Integer) = DirichletCategorical(ones(k), zeros(Int, k))
DirichletCategorical(α::Vector{<: Real}) = DirichletCategorical(α, zeros(Int, length(α)))

Base.show(io::IO, d::DirichletCategorical) = show(io, d, (:α, :n))

ncategories(d::DirichletCategorical) = length(d.α)
length(d::DirichletCategorical) = ncategories(d)
params(d::DirichletCategorical) = (d.α, d.n)
@inline partype(d::DirichletCategorical{T}) where {T} = T

function rand(rng::AbstractRNG, d::DirichletCategorical)
    quantile(d, rand(rng))
end

function pdf(d::DirichletCategorical, x::Int)
    if insupport(d, x)
        @inbounds (d.α[x] + d.n[x]) / (d.α0 + sum(d.n))
    else
        zero(x)
    end
end

function logpdf(d::DirichletCategorical, x::Int)
    if insupport(d, x)
        @inbounds log(d.α[x] + d.n[x]) - log(d.α0 + sum(d.n))
    else
        -Inf
    end
end

function cdf(d::DirichletCategorical, x::Int)
    if x < minimum(d)
        zero(x)
    elseif x > maximum(d)
        one(x)
    else
        @inbounds sum((d.α[1:x] + d.n[1:x]) / (d.α0 + sum(d.n)))
    end
end

function quantile(d::DirichletCategorical, q::Real)
    0 <= q <= 1 || throw(DomainError(q, "q ∈ [0,1]"))
    p = probs(d)
    cp = p[1]
    i = 1
    while cp < q
        i += 1
        @inbounds cp += p[i]
    end
    i
end

probs(d::DirichletCategorical) = (d.n + d.α) .* inv(d.α0 + sum(d.n))

function minimum(d::DirichletCategorical)
    return 1
end

function maximum(d::DirichletCategorical)
    return length(d.α)
end

function insupport(d::DirichletCategorical, x::Real)
    return x >= minimum(d) && x <= maximum(d)
end

function mode(d::DirichletCategorical)
    m, i = findmax(dirichlet_mode(d.α + d.n, d.α0 + sum(d.n)))
    i
end

function update!(d::DirichletCategorical, obs::Int)
    insupport(d, obs) || throw(ArgumentError("out of bounds observation"))

    @inbounds d.n[obs] += 1
end

function update!(d::DirichletCategorical, obs::Vector{Int})
    all(map(x -> insupport(d, x), obs)) || throw(ArgumentError("out of bounds observation"))

    for o in obs
        @inbounds d.n[o] += 1
    end
end


text:
# @show LOAD_PATH
# push!(LOAD_PATH, expanduser("~/doc/SearchBall/"))
# @show LOAD_PATH
import SearchBall; const sb = SearchBall

# config = "five_vs_five1"
config = "startposition01"
# config = "one_vs_many1"
# config = "symmetric1"
# config = "test"
# config = "2dlinNN"
# config = "2dlinNN01"
# config = "two_vs_three"

#TODO: make runnable from everywhere
configFile = joinpath(pwd(), "configs", string(config, ".config"))
# configFile = joinpath("/home/tobias/doc/SearchBall/SearchBall.jl/", "configs", string(config, ".config"))
configFile = normpath(configFile)
sb.main(configFile)


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 06b847e0-fcb2-11ea-1880-0543a2db1b1c
begin
	using Dierckx
	using Plots
end

# ╔═╡ 24e001e0-fcb2-11ea-0559-1b3da141b5a3
rAti = [0, 10, 20, 50, 100, 200, 500, 1e3, 2e3, 5e3, 1e4]

# ╔═╡ 9670e220-fcb2-11ea-37be-4f9bfa6437b1
zAti = [sin.(r/1e3) for r ∈ rAti]

# ╔═╡ a228e9f0-fcb2-11ea-1cd0-a7ef10379149
plot(rAti, zAti)

# ╔═╡ 2d598b10-fcb3-11ea-1579-4709585b2f96
# itpAti = LinearInterpolation(rAti, zAti)
splAti = Spline1D(rAti, zAti)

# ╔═╡ de70a04e-fcb3-11ea-1e68-3be191449dfb
# zAtiFcn(r) = itpAti(r)
zAtiFcn(r) = splAti(r)

# ╔═╡ 42fec5be-fcb3-11ea-2d2b-7df1958c48cf
begin
	r = range(0, rAti[end], length = 101)
	scatter(r, zAtiFcn)
	plot!(rAti, zAti)
end

# ╔═╡ fae441b0-fcb3-11ea-2baf-09505779f9a3
begin
	function InterpolatingFunction(rng, val)
		Itp = LinearInterpolation(rng, val)
		return ItpFcn(r) = Itp(r)
	end

	function InterpolatingFunction(rng, dpt, val)
		Itp = LinearInterpolation((dpt, rng), val)
		return ItpFcn(r, z) = Itp(z, r)
	end
end

# ╔═╡ 4ac0e5d0-fcb4-11ea-1477-0be39085827b
zAtiFcn2 = InterpolatingFunction(rAti, zAti)

# ╔═╡ 5cf89810-fcb4-11ea-194b-235bb415165c
scatter(r, zAtiFcn2)

# ╔═╡ 66914250-fcb4-11ea-20e3-07319a79a510
begin
	rc = rAti
	zc = [0, 10, 20, 100, 500, 1e3]
# 	c = 1550 .+ sin.(rc'/1e3) .+ sin.(zc/1e3)
	c = [1550 + sin(r/1e3) + sin(z/1e3) for z ∈ zc, r ∈ rc]
end

# ╔═╡ b45b8130-fcb4-11ea-0801-bd0afb438aca
heatmap(rc, zc, c)

# ╔═╡ cd252d60-fcb4-11ea-3112-d7f36a0f4c04
begin
	cItp = LinearInterpolation((zc, rc), c)
	cFcn(r, z) = cItp(z, r)
end

# ╔═╡ f3661d90-fcb4-11ea-0352-a5284b9dc242
begin
	z = range(0, zc[end], length = 101)
	heatmap(r, z, cFcn)
end

# ╔═╡ 28f73070-fcb5-11ea-0148-91f62d5ef3d7
begin
	zInd = 6
	zVal = zc[zInd]
	scatter(rc, c[zInd, :])
	plot!(r, r -> cFcn(r, zVal))
end

# ╔═╡ 5472ef00-fcb5-11ea-0640-672ede9c99b3
cFcn2 = InterpolatingFunction(rc, zc, c)

# ╔═╡ bd6c4be0-fcb6-11ea-1343-c1584a8f72d9
heatmap(r, z, cFcn2)

# ╔═╡ 984269a2-fcb9-11ea-04f1-0927fe41dfa3
cMat = vcat(hcat(0, transpose(rc)), hcat(zc, c))

# ╔═╡ 118e5a80-fcba-11ea-3a27-671c712ec008
begin
	rVec = [r for r ∈ cMat[1, 2:end]]
	zVec = [z for z ∈ cMat[2:end, 1]]
	cVec = cMat[2:end, 2:end]
end

# ╔═╡ b9c1e4ae-fcba-11ea-2d72-4778c0a74263
size(rVec), size(zVec), size(cVec)

# ╔═╡ 6adf09e0-fcba-11ea-2d8b-e5bbbece2d99
cFcn3 = InterpolatingFunction(rVec, zVec, cVec)

# ╔═╡ 4648ee10-fcbb-11ea-2247-735e9cb6e01d
heatmap(r, z, cFcn3)

# ╔═╡ Cell order:
# ╠═06b847e0-fcb2-11ea-1880-0543a2db1b1c
# ╠═24e001e0-fcb2-11ea-0559-1b3da141b5a3
# ╠═9670e220-fcb2-11ea-37be-4f9bfa6437b1
# ╠═a228e9f0-fcb2-11ea-1cd0-a7ef10379149
# ╠═2d598b10-fcb3-11ea-1579-4709585b2f96
# ╠═de70a04e-fcb3-11ea-1e68-3be191449dfb
# ╠═42fec5be-fcb3-11ea-2d2b-7df1958c48cf
# ╠═fae441b0-fcb3-11ea-2baf-09505779f9a3
# ╠═4ac0e5d0-fcb4-11ea-1477-0be39085827b
# ╠═5cf89810-fcb4-11ea-194b-235bb415165c
# ╠═66914250-fcb4-11ea-20e3-07319a79a510
# ╠═b45b8130-fcb4-11ea-0801-bd0afb438aca
# ╠═cd252d60-fcb4-11ea-3112-d7f36a0f4c04
# ╠═f3661d90-fcb4-11ea-0352-a5284b9dc242
# ╠═28f73070-fcb5-11ea-0148-91f62d5ef3d7
# ╠═5472ef00-fcb5-11ea-0640-672ede9c99b3
# ╠═bd6c4be0-fcb6-11ea-1343-c1584a8f72d9
# ╠═984269a2-fcb9-11ea-04f1-0927fe41dfa3
# ╠═118e5a80-fcba-11ea-3a27-671c712ec008
# ╠═b9c1e4ae-fcba-11ea-2d72-4778c0a74263
# ╠═6adf09e0-fcba-11ea-2d8b-e5bbbece2d99
# ╠═4648ee10-fcbb-11ea-2247-735e9cb6e01d


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

"""
    SaveSolutionCallback(; interval=0,
                           save_initial_solution=true,
                           save_final_solution=true,
                           output_directory="out",
                           solution_variables=cons2prim)

Save the current numerical solution every `interval` time steps. `solution_variables` can be any
callable that converts the conservative variables at a single point to a set of solution variables.
The first parameter passed to `solution_variables` will be the set of conservative variables and the
second parameter is the equation struct.
"""
mutable struct SaveSolutionCallback{SolutionVariables}
  interval::Int
  save_initial_solution::Bool
  save_final_solution::Bool
  output_directory::String
  solution_variables::SolutionVariables
end

function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:SaveSolutionCallback})
  @nospecialize cb # reduce precompilation time

  save_solution_callback = cb.affect!
  print(io, "SaveSolutionCallback(interval=", save_solution_callback.interval, ")")
end

function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{<:Any, <:SaveSolutionCallback})
  @nospecialize cb # reduce precompilation time

  if get(io, :compact, false)
    show(io, cb)
  else
    save_solution_callback = cb.affect!

    setup = [
             "interval" => save_solution_callback.interval,
             "solution variables" => save_solution_callback.solution_variables,
             "save initial solution" => save_solution_callback.save_initial_solution ? "yes" : "no",
             "save final solution" => save_solution_callback.save_final_solution ? "yes" : "no",
             "output directory" => abspath(normpath(save_solution_callback.output_directory)),
            ]
    summary_box(io, "SaveSolutionCallback", setup)
  end
end

function SaveSolutionCallback(; interval=0,
                                save_initial_solution=true,
                                save_final_solution=true,
                                output_directory="out",
                                solution_variables=cons2prim)

  # FIXME: Deprecations introduced in v0.3
  if solution_variables isa Symbol
    Base.depwarn("Providing the keyword argument `solution_variables` as a `Symbol` is deprecated." *
                 "Use functions such as `cons2cons` or `cons2prim` instead.", :SaveSolutionCallback)
    if solution_variables == :conservative
      solution_variables = cons2cons
    elseif solution_variables == :primitive
      solution_variables = cons2prim
    else
      error("Unknown `solution_variables` $solution_variables.")
    end
  end

  solution_callback = SaveSolutionCallback(interval, save_initial_solution, save_final_solution,
                                           output_directory, solution_variables)

  DiscreteCallback(solution_callback, solution_callback, # the first one is the condition, the second the affect!
                   save_positions=(false,false),
                   initialize=initialize!)
end

function initialize!(cb::DiscreteCallback{Condition,Affect!}, u, t, integrator) where {Condition, Affect!<:SaveSolutionCallback}
  solution_callback = cb.affect!

  mpi_isroot() && mkpath(solution_callback.output_directory)

  semi = integrator.p
  mesh, _, _, _ = mesh_equations_solver_cache(semi)
  @trixi_timeit timer() "I/O" begin
    if mesh.unsaved_changes
      mesh.current_filename = save_mesh_file(mesh, solution_callback.output_directory)
      mesh.unsaved_changes = false
    end
  end

  if solution_callback.save_initial_solution
    solution_callback(integrator)
  end

  return nothing
end

# this method is called to determine whether the callback should be activated
function (solution_callback::SaveSolutionCallback)(u, t, integrator)
  @unpack interval, save_final_solution = solution_callback

  # With error-based step size control, some steps can be rejected. Thus,
  #   `integrator.iter >= integrator.destats.naccept`
  #    (total #steps)       (#accepted steps)
  # We need to check the number of accepted steps since callbacks are not
  # activated after a rejected step.
  return interval > 0 && (
    ((integrator.destats.naccept % interval == 0) && !(integrator.destats.naccept == 0 && integrator.iter > 0)) ||
    (save_final_solution && isfinished(integrator)))
end

# this method is called when the callback is activated
function (solution_callback::SaveSolutionCallback)(integrator)
  u_ode = integrator.u
  @unpack t, dt = integrator
  iter = integrator.destats.naccept
  semi = integrator.p
  mesh, _, _, _ = mesh_equations_solver_cache(semi)

  @trixi_timeit timer() "I/O" begin
    @trixi_timeit timer() "save mesh" if mesh.unsaved_changes
      mesh.current_filename = save_mesh_file(mesh, solution_callback.output_directory, iter)
      mesh.unsaved_changes = false
    end

    element_variables = Dict{Symbol, Any}()
    @trixi_timeit timer() "get element variables" begin
      get_element_variables!(element_variables, u_ode, semi)
      callbacks = integrator.opts.callback
      if callbacks isa CallbackSet
        for cb in callbacks.continuous_callbacks
          get_element_variables!(element_variables, u_ode, semi, cb; t=integrator.t, iter=integrator.destats.naccept)
        end
        for cb in callbacks.discrete_callbacks
          get_element_variables!(element_variables, u_ode, semi, cb; t=integrator.t, iter=integrator.destats.naccept)
        end
      end
    end

    @trixi_timeit timer() "save solution" save_solution_file(u_ode, t, dt, iter, semi, solution_callback, element_variables)
  end

  # avoid re-evaluating possible FSAL stages
  u_modified!(integrator, false)
  return nothing
end

@inline function save_solution_file(u_ode, t, dt, iter,
                                    semi::AbstractSemidiscretization, solution_callback,
                                    element_variables=Dict{Symbol,Any}())
  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
  u = wrap_array_native(u_ode, mesh, equations, solver, cache)
  save_solution_file(u, t, dt, iter, mesh, equations, solver, cache, solution_callback, element_variables)
end

# TODO: Taal refactor, move save_mesh_file?
# function save_mesh_file(mesh::TreeMesh, output_directory, timestep=-1) in io/io.jl

include("save_solution_dg.jl")

end # @muladd


text:
# ccall wrapper for calling functions in libraries that might not be available

export @runtime_ccall, decode_ccall_function, @checked, @workspace, @argout

"""
    @runtime_ccall((function_name, library), returntype, (argtype1, ...), argvalue1, ...)

Extension of `ccall` that performs the lookup of `function_name` in `library` at run time.
This is useful in the case that `library` might not be available, in which case a function
that performs a `ccall` to that library would fail to compile.

After a slower first call to load the library and look up the function, no additional
overhead is expected compared to regular `ccall`.
"""
macro runtime_ccall(target, args...)
    # decode ccall function/library target
    Meta.isexpr(target, :tuple) || error("Expected (function_name, library) tuple")
    function_name, library = target.args

    # global const ref to hold the function pointer
    @gensym fptr_cache
    @eval __module__ begin
        # uses atomics (release store, acquire load) for thread safety.
        # see https://github.com/JuliaGPU/CUDAapi.jl/issues/106 for details
        const $fptr_cache = Threads.Atomic{UInt}(0)
    end

    return quote
        # use a closure to hold the lookup and avoid code bloat in the caller
        @noinline function cache_fptr!()
            library = Libdl.dlopen($(esc(library)))
            $(esc(fptr_cache))[] = Libdl.dlsym(library, $(esc(function_name)))

            $(esc(fptr_cache))[]
        end

        fptr = $(esc(fptr_cache))[]
        if fptr == 0        # folded into the null check performed by ccall
            fptr = cache_fptr!()
        end

        ccall(reinterpret(Ptr{Cvoid}, fptr), $(map(esc, args)...))
    end

    return
end

# decode `ccall` or `@runtime_ccall` invocations and extract the function that is called
function decode_ccall_function(ex)
    # check is used in front of `ccall` or `@runtime_ccall`s that work on a tuple (fun, lib)
    if Meta.isexpr(ex, :call)
        @assert ex.args[1] == :ccall
        @assert Meta.isexpr(ex.args[2], :tuple)
        fun = String(ex.args[2].args[1].value)
    elseif Meta.isexpr(ex, :macrocall)
        @assert ex.args[1] == Symbol("@runtime_ccall")
        @assert Meta.isexpr(ex.args[3], :tuple)
        fun = String(ex.args[3].args[1].value)
    else
        error("@check should prefix ccall or @runtime_ccall")
    end

    # strip any version tag (e.g. cuEventDestroy_v2 -> cuEventDestroy)
    m = match(r"_v\d+$", fun)
    if m !== nothing
        fun = fun[1:end-length(m.match)]
    end

    return fun
end

"""
    @checked function foo(...)
        rv = ...
        return rv
    end

Macro for wrapping a function definition returning a status code. Two versions of the
function will be generated: `foo`, with the function body wrapped by an invocation of the
`@check` macro (to be implemented by the caller of this macro), and `unsafe_foo` where no
such invocation is present and the status code is returned to the caller.
"""
macro checked(ex)
    # parse the function definition
    @assert Meta.isexpr(ex, :function)
    sig = ex.args[1]
    @assert Meta.isexpr(sig, :call)
    body = ex.args[2]
    @assert Meta.isexpr(body, :block)

    # generate a "safe" version that performs a check
    safe_body = quote
        @check $body
    end
    safe_sig = Expr(:call, sig.args[1], sig.args[2:end]...)
    safe_def = Expr(:function, safe_sig, safe_body)

    # generate a "unsafe" version that returns the error code instead
    unsafe_sig = Expr(:call, Symbol("unsafe_", sig.args[1]), sig.args[2:end]...)
    unsafe_def = Expr(:function, unsafe_sig, body)

    return esc(:($safe_def, $unsafe_def))
end

"""
    @argout call(..., out(output_arg), ...)

Change the behavior of a function call, returning the value of an argument instead.
A common use case is to pass a newly-created Ref and immediately dereference that output:

    @argout(some_getter(Ref{Int}()))[]

If no output argument is specified, `nothing` will be returned.
Multiple output arguments return a tuple.

"""
macro argout(ex)
    Meta.isexpr(ex, :call) || throw(ArgumentError("@argout macro should be applied to a function call"))

    block = quote end

    # look for output arguments (`out(...)`)
    output_vars = []
    args = ex.args[2:end]
    for (i,arg) in enumerate(args)
        if Meta.isexpr(arg, :call) && arg.args[1] == :out
            # allocate a variable
            @gensym output_val
            push!(block.args, :($output_val = $(ex.args[i+1].args[2]))) # strip `output(...)`
            push!(output_vars, output_val)

            # replace the argument
            ex.args[i+1] = output_val
        end
    end

    # generate a return
    push!(block.args, ex)
    if isempty(output_vars)
        push!(block.args, :(nothing))
    elseif length(output_vars) == 1
        push!(block.args, :($(output_vars[1])))
    else
        push!(block.args, :(tuple($(output_vars...))))
    end

    esc(block)
end


text:
function export_shape_functions(ap::AbstractSeg;inte::Int=100)
    # coordinates = Dict{String,Vector{Float64}}()
    # push!(coordinates,"x"=>zeros(inte+1))
    # push!(coordinates,"y"=>zeros(inte+1))
    # push!(coordinates,"z"=>zeros(inte+1))
    # datas = Dict{Int,Vector{Float64}}()
    x = zeros(inte+1)
    max_id = maximum(ap.id)
    # push!(datas,i=>zeros(inte+1) for i in 1:max_id)
    datas = zeros(inte+1,max_id)
    for j in 1:inte+1
        ξⱼ = -1+2/inte*(j-1)
        xⱼ = get_coordinates(ap,ξⱼ)
        # coordinates["x"][j] = xᵢ[1]
        # coordinates["y"][j] = xᵢ[2]
        # coordinates["z"][j] = xᵢ[3]
        x[j] = xⱼ[1]
        N = get_shape_functions(ap,ξⱼ,:∂1)
        for i in 1:get_number_of_indices
            I = get_global_indice(ap,i)
            # datas[I][j] = N[i]
            datas[j,I] = N[i]
        end
    end
    # return coordinates, datas
    return x, datas
end

export_shape_functions(aps::Vector{AbstractSeg},gs::Symbol...;inte::Int=100)= export_shape_functions(aps,inte,Val.(gs)...)

function export_shape_functions(aps::Vector{AbstractSeg},inte::Int,::Val{:∂1})
    nₑ = length(aps)
    x = zeros(nₑ*(inte+1))
    max_id = 0
    for ap in aps
        max_temp = maximum(ap.id)
        max_id = max_id > max_temp ? max_id : max_temp
    end
    N = zeros(nₑ*(inte+1),max_id)
    for n in 1:nₑ
        ap = aps[n]
        for j in 1:inte+1
            ξⱼ = -1+2/inte*(j-1)
            xⱼ = get_coordinates(ap,ξⱼ)
            x[(inte+1)*(n-1)+j] = xⱼ[1]
            N_ = get_shape_functions(ap,ξⱼ,:∂1)
            for i in 1:get_number_of_indices(ap)
                I = get_global_indice(ap,i)
                N[(inte+1)*(n-1)+j,I] = N_[i]
            end
        end
    end
    return x, N
end

function export_shape_functions(aps::Vector{AbstractSeg},inte::Int,::Val{:∂1},::Val{:∂x})
    nₑ = length(aps)
    x = zeros(nₑ*(inte+1))
    max_id = 0
    for ap in aps
        max_temp = maximum(ap.id)
        max_id = max_id > max_temp ? max_id : max_temp
    end
    N = zeros(nₑ*(inte+1),max_id)
    ∂N∂x = zeros(nₑ*(inte+1),max_id)
    for n in 1:nₑ
        ap = aps[n]
        for j in 1:inte+1
            ξⱼ = -1+2/inte*(j-1)
            xⱼ = get_coordinates(ap,ξⱼ)
            x[(inte+1)*(n-1)+j] = xⱼ[1]
            N_, ∂N∂x_ = get_shape_functions(ap,ξⱼ,:∂1,:∂x)
            for i in 1:get_number_of_indices(ap)
                I = get_global_indice(ap,i)
                N[(inte+1)*(n-1)+j,I] = N_[i]
                ∂N∂x[(inte+1)*(n-1)+j,I] = ∂N∂x_[i]
            end
        end
    end
    return x, N, ∂N∂x
end

function export_shape_functions(aps::Vector{AbstractSeg},inte::Int,::Val{:∂1},::Val{:∂x},::Val{:∂y},::Val{:∂z})
    nₑ = length(aps)
    x = zeros(nₑ*(inte+1))
    max_id = 0
    for ap in aps
        max_temp = maximum(ap.id)
        max_id = max_id > max_temp ? max_id : max_temp
    end
    N = zeros(nₑ*(inte+1),max_id)
    ∂N∂x = zeros(nₑ*(inte+1),max_id)
    ∂N∂y = zeros(nₑ*(inte+1),max_id)
    ∂N∂z = zeros(nₑ*(inte+1),max_id)
    for n in 1:nₑ
        ap = aps[n]
        for j in 1:inte+1
            ξⱼ = -1+2/inte*(j-1)
            xⱼ = get_coordinates(ap,ξⱼ)
            x[(inte+1)*(n-1)+j] = xⱼ[1]
            N_, ∂N∂x_, ∂N∂y_, ∂N∂z_ = get_shape_functions(ap,ξⱼ,:∂1,:∂x,:∂y,:∂z)
            for i in 1:get_number_of_indices(ap)
                I = get_global_indice(ap,i)
                N[(inte+1)*(n-1)+j,I] = N_[i]
                ∂N∂x[(inte+1)*(n-1)+j,I] = ∂N∂x_[i]
                ∂N∂y[(inte+1)*(n-1)+j,I] = ∂N∂y_[i]
                ∂N∂z[(inte+1)*(n-1)+j,I] = ∂N∂z_[i]
            end
        end
    end
    return x, N, ∂N∂x, ∂N∂y, ∂N∂z
end

## VTK
const cell_type = Dict{Any,Int}(AbstractPoi=>1,AbstractSeg=>3,AbstractTri=>5,AbstractQuad=>9)

mutable struct VTKExport
    filename::String
    topic::String
    celltype::String
    pointdata::Union{Dict{String,Symbol},Nothing}
    celldata::Union{Dict{String,Symbol},Nothing}
    parameters::Union{Dict{String,Float64},Nothing}
end
VTKExport(;filename::String="default.vtk",topic::String="vtk file",celltype::String="UNSTRUCTURED_GRID",pointdata::Union{Dict{String,Symbol},Nothing}=nothing,celldata::Union{Dict{String,Symbol},Nothing}=nothing,parameters::Union{Dict{String,Float64},Nothing}=nothing) = VTKExport(filename,topic,celltype,pointdata,celldata,parameters)

function (op::VTKExport)(aps::Vector{Approximator},d::AbstractVector{Float64})
    fid = open(op.filename,"w")
    nₑ = length(aps)
    nₚ = length(aps[1].nodes)
    write(fid,"# vtk DataFile Version 2.0\n")
    write(fid,op.topic*"\n")
    write(fid,"ASCII\n")
    write(fid,"DATASET "*op.celltype*"\n")
    # POINTS
    write(fid,"POINTS $nₚ float\n")
    for node in aps[1].nodes
        x = node.x
        y = node.y
        z = node.z
        write(fid,"$x $y $z\n")
    end

    # CELLS
    op(fid,aps,Val(Meta.parse(op.celltype)))

    # POINT_DATA
    if op.pointdata ≠ nothing
        write(fid,"POINT_DATA $nₚ\n")
        for (name,datatype) in op.pointdata
            op(fid,aps,d,name,Val(datatype))
        end
    end
    # CELL_DATA
    if op.celldata ≠ nothing
        write(fid,"CELL_DATA $nₑ\n")
        for (name,datatype) in op.celldata
            op(fid,aps,d,name,Val(datatype))
        end
    end
    close(fid)
end

function (op::VTKExport)(fid::IO,aps::Vector{Approximator},::Val{:UNSTRUCTURED_GRID})
    nₑ = length(aps)
    n = nₑ
    for ap in aps
        n += get_number_of_indices(ap)
    end
    write(fid,"CELLS $nₑ $n\n")
    for ap in aps
        nᵢ = get_number_of_indices(ap)
        str = "$nᵢ"
        for I in ap.id .- 1
            str *= " $I"
        end
        str *= "\n"
        write(fid,str)
    end
    write(fid,"CELL_TYPES $nₑ\n")
    for ap in aps
        ntag = cell_type[supertype(typeof(ap))]
        write(fid,"$ntag\n")
    end
end

function (::VTKExport)(fid::IO,aps::Vector{Approximator},d::AbstractVector{Float64},name::String,::Val{:Scale})
    write(fid,"SCALARS "*name*" float 1\n")
    write(fid,"LOOKUP_TABLE "*name*"\n")
    for u in d
        write(fid,"$u\n")
    end
end

function (::VTKExport)(fid::IO,aps::Vector{Approximator},d::AbstractVector{Float64},name::String,::Val{:u_2D})
    nₚ = length(aps[1].nodes)
    write(fid,"VECTORS "*name*" float\n")
    for i in 1:nₚ
        u = d[2*i-1]
        v = d[2*i]
        w = 0.0
        write(fid,"$u $v $w\n")
    end
end

function (::VTKExport)(fid::IO,aps::Vector{Approximator},d::AbstractVector{Float64},name::String,::Val{:u})
    nₚ = length(aps[1].nodes)
    write(fid,"VECTORS "*name*" float\n")
    for i in 1:nₚ
        u = d[3*i-2]
        v = d[3*i-1]
        w = d[3*i]
        write(fid,"$u $v $w\n")
    end
end

function (op::VTKExport)(fid::IO,aps::Vector{Approximator},d::AbstractVector{Float64},name::String,::Val{:σ_PlaneStress})
    write(fid,"TENSORS "*name*" float\n")
    E = op.parameters["E"]
    ν = op.parameters["ν"]
    Cᵢᵢᵢᵢ = E/(1-ν^2)
    Cᵢᵢⱼⱼ = E*ν/(1-ν^2)
    Cᵢⱼᵢⱼ = E/2/(1+ν)
    for ap in aps
        B₁,B₂ = get_shape_functions(ap,SVector{3,Float64}(0.,0.,0.),Val(:∂x),Val(:∂y))
        ε₁₁ = 0.
        ε₂₂ = 0.
        ε₁₂ = 0.
        for i in 1:get_number_of_indices(ap)
            I = get_global_indice(ap,i)
            ε₁₁ += B₁[i]*d[2*I-1]
            ε₂₂ += B₂[i]*d[2*I]
            ε₁₂ += B₁[i]*d[2*I] + B₂[i]*d[2*I-1]
        end
        σ₁₁ = Cᵢᵢᵢᵢ*ε₁₁ + Cᵢᵢⱼⱼ*ε₂₂
        σ₂₂ = Cᵢᵢᵢᵢ*ε₂₂ + Cᵢᵢⱼⱼ*ε₁₁
        σ₁₂ = Cᵢⱼᵢⱼ*ε₁₂
        write(fid,"$σ₁₁ $σ₁₂ 0.0\n")
        write(fid,"$σ₁₂ $σ₂₂ 0.0\n")
        write(fid,"0.0 0.0 0.0\n")
        write(fid,"\n")
    end
end

function (op::VTKExport)(fid::IO,aps::Vector{Approximator},d::AbstractVector{Float64},name::String,::Val{:σ_PlaneStrain})
    write(fid,"TENSORS "*name*" float\n")
    E = op.parameters["E"]
    ν = op.parameters["ν"]
    Cᵢᵢᵢᵢ = E*(1-ν)/(1+ν)/(1-ν^2)
    Cᵢᵢⱼⱼ = E*ν/(1+ν)/(1-ν^2)
    Cᵢⱼᵢⱼ = E/2/(1+ν)
    for ap in aps
        B₁,B₂ = get_shape_functions(ap,SVector{3,Float64}(0.,0.,0.),Val(:∂x),Val(:∂y))
        ε₁₁ = 0.
        ε₂₂ = 0.
        ε₁₂ = 0.
        for i in 1:get_number_of_indices(ap)
            I = get_global_indice(ap,i)
            ε₁₁ += B₁[i]*d[2*I-1]
            ε₂₂ += B₂[i]*d[2*I]
            ε₁₂ += B₁[i]*d[2*I] + B₂[i]*d[2*I-1]
        end
        σ₁₁ = Cᵢᵢᵢᵢ*ε₁₁ + Cᵢᵢⱼⱼ*ε₂₂
        σ₂₂ = Cᵢᵢᵢᵢ*ε₂₂ + Cᵢᵢⱼⱼ*ε₁₁
        σ₁₂ = Cᵢⱼᵢⱼ*ε₁₂
        σ₃₃ = Cᵢᵢⱼⱼ*(ε₁₁+ε₂₂)
        write(fid,"$σ₁₁ $σ₁₂ 0.0\n")
        write(fid,"$σ₁₂ $σ₂₂ 0.0\n")
        write(fid,"0.0 0.0 $σ₃₃\n")
        write(fid,"\n")
    end
end

function (op::VTKExport)(fid::IO,aps::Vector{Approximator},d::AbstractVector{Float64},name::String,::Val{:σ})
    write(fid,"TENSORS "*name*" float\n")
    E = op.parameters["E"]
    ν = op.parameters["ν"]
    Cᵢᵢᵢᵢ = E*(1-ν)/(1+ν)/(1-ν^2)
    Cᵢᵢⱼⱼ = E*ν/(1+ν)/(1-ν^2)
    Cᵢⱼᵢⱼ = E/2/(1+ν)
    for ap in aps
        B₁,B₂,B₃ = get_shape_functions(ap,SVector{3,Float64}(0.,0.,0.),Val(:∂x),Val(:∂y),Val(:∂z))
        ε₁₁ = 0.
        ε₂₂ = 0.
        ε₃₃ = 0.
        ε₁₂ = 0.
        ε₁₃ = 0.
        ε₂₃ = 0.
        for i in 1:get_number_of_indices(ap)
            I = get_global_indice(ap,i)
            ε₁₁ += B₁[i]*d[3*I-2]
            ε₂₂ += B₂[i]*d[3*I-1]
            ε₃₃ += B₃[i]*d[3*I]
            ε₁₂ += B₁[i]*d[3*I-1] + B₂[i]*d[3*I-2]
            ε₁₃ += B₁[i]*d[3*I] + B₃[i]*d[3*I-2]
            ε₂₃ += B₂[i]*d[3*I] + B₃[i]*d[3*I-1]
        end
        σ₁₁ = Cᵢᵢᵢᵢ*ε₁₁ + Cᵢᵢⱼⱼ*ε₂₂ + Cᵢᵢⱼⱼ*ε₃₃
        σ₂₂ = Cᵢᵢⱼⱼ*ε₁₁ + Cᵢᵢᵢᵢ*ε₂₂ + Cᵢᵢⱼⱼ*ε₃₃
        σ₃₃ = Cᵢᵢⱼⱼ*ε₁₁ + Cᵢᵢⱼⱼ*ε₂₂ + Cᵢᵢᵢᵢ*ε₃₃
        σ₁₂ = Cᵢⱼᵢⱼ*ε₁₂
        σ₁₃ = Cᵢⱼᵢⱼ*ε₁₃
        σ₂₃ = Cᵢⱼᵢⱼ*ε₂₃
        write(fid,"$σ₁₁ $σ₁₂ $σ₁₃\n")
        write(fid,"$σ₁₂ $σ₂₂ $σ₂₃\n")
        write(fid,"$σ₁₃ $σ₂₃ $σ₃₃\n")
        write(fid,"\n")
    end
end
## Phase field export
function (op::VTKExport)(aps::Vector{Approximator},d::AbstractVector{Float64},dᵥ::AbstractVector{Float64},::Val{:PFM_PlaneStress})
    fid = open(op.filename,"w")
    nₑ = length(aps)
    nₚ = length(aps[1].nodes)
    write(fid,"# vtk DataFile Version 2.0\n")
    write(fid,op.topic*"\n")
    write(fid,"ASCII\n")
    write(fid,"DATASET "*op.celltype*"\n")
    # POINTS
    write(fid,"POINTS $nₚ float\n")
    for node in aps[1].nodes
        x = node.x[1]
        y = node.x[2]
        z = node.x[3]
        write(fid,"$x $y $z\n")
    end

    # CELLS
    op(fid,aps,Val(Meta.parse(op.celltype)))

    # POINT_DATA
    write(fid,"POINT_DATA $nₚ\n")
    op(fid,aps,d,"disp",Val(:u_2D))
    op(fid,aps,dᵥ,"damage",Val(:Scale))

    # CELL_DATA
    write(fid,"CELL_DATA $nₑ\n")
    op(fid,aps,d,"stress",Val(:σ_PlaneStress))
    close(fid)
end

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
module test_OnlineAI

# reload("OnlineAI")
using OnlineAI
# include("../src/utils.jl")

function testxor(maxiter::Int)

  # all xor inputs and results
  inputs = float([0 0; 0 1; 1 0; 1 1])
  targets = Float64[(sum(row(inputs,i))==1.0)*0.8+0.1 for i in 1:size(inputs,1)]

  # all sets are the same
  data = buildSolverData(inputs, targets)
  datasets = DataSets(data, data, data)

  net = buildNeuralNet([2,2,1]; η=0.5, μ=0.1)

  params = buildSolverParams(maxiter=maxiter, minerror=1e-6)
  solve!(net, params, datasets)

  for d in data
    output = feedforward!(net, d.input)
    println("Result: input=$(d.input) target=$(d.target) output=$output")
  end

  net
end

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

using Qwt

sinfunction(x) = sin(x) / 2 + sin(x*4) / 4
# sinfunction(x) = 0.9 * x^2 / maxx^2
# sinfunction(x) = x^2

function updatesinplot(net::NeuralNet, params::SolverParams, datasets::DataSets, stats::SolverStats)

  setdata(plt.lines[1], xxx, Float64[predict(net, [f])[1] for f in xxx])
  setdata(plt.lines[2], xxx, map(sinfunction,xxx))
  title(plt, @sprintf("Update: %4d  Error: %10.4f", stats.numiter, stats.validationError))

  push!(plt2, 1, stats.numiter, stats.validationError)

  refresh(plt)
  refresh(plt2)

  shoulddogif && saveframe(anim)
end

function buildRegNN1()
  buildNeuralNet([1,10,10,1]; η=0.005, μ=0.1, activation=TanhMapping())
end

function buildRegNN2()
  η=0.02
  μ=0.2
  NeuralNet(Layer[
    buildLayer(1, 10),
    buildLayer(10, 10),
    buildLayer(10, 1, IdentityMapping())
    ], η, μ)
end

function testsin(maxiter::Int, dogif = false)

  global maxx = 10.0
  global sp = subplot(zeros(0,2), ncols=1)
  global plt = oplot(sp.plots[1], zeros(1,2), labels=["estimate","target"])
  global plt2 = oplot(sp.plots[2], zeros(1), xlabel="Iteration", ylabel="Network Error"); empty!(plt2)
  global xxx = collect(-maxx:.1:maxx)
  global anim = animation(sp, "/tmp/png")
  global shoulddogif = dogif

  inputs = collect(linspace(-maxx,maxx,1000))
  inputs = reshape(inputs, length(inputs), 1)
  targets = map(sinfunction, inputs)

  # all sets are the same
  data = buildSolverData(inputs, targets)
  datasets = DataSets(data, sample(data,30), data)

  # net = buildNeuralNet([1,10,10,1]; η=0.005, μ=0.1, activation=TanhMapping())
  net = buildRegressionNet(1, 1, [30,30]; hiddenMapping = ReLUMapping(), solver=NNetSolver(η = 0.05, dropout=DropoutStrategy(on=false)))

  params = SolverParams(maxiter=maxiter, onbreak=updatesinplot, displayiter=10000, minerror=1e-3)
  solve!(net, params, datasets)

  shoulddogif && makegif(anim)
  net
end

end


text:
module MuttsTest

using MuttsInterface
using Test
using MacroTools

#mutable struct Foo <: Mutt
#    x :: Int
#    y :: Int
#    __mutt_mutable :: Bool
#end
#
#Foo(x :: Int, y :: Int) = Foo(x, y, true)
#Foo(f :: Foo) = Foo(f.x, f.y)
#
#f = Foo(3,5)
#@assert is_mutts_mutable(f)
#f.x = 4
#
#Base.copy(v :: Foo) = Foo(v.x, v.y)
#g = branch!(f)
#@assert !is_mutts_mutable(f)
#@assert is_mutts_mutable(g)
#mark_immutable!(g)
#
#@assert !is_mutts_mutable(g)

#mutable struct Bar <: Mutt
#    f :: Foo
#    z :: Float64
#    __mutt_mutable :: Bool
#end
#
#Bar(f :: Foo, z :: Float64) = Bar(f, z, true)
#Bar(b :: Bar) = Bar(f, z)
#
#

# Mutts macro expansion

@testset "macro parsing" begin
    @test @macroexpand(@mutt struct S end) isa Expr
    @test @macroexpand(@mutt struct S x end) isa Expr

    # Bad parses
    @test_throws Exception @macroexpand(@mutt 1)
    @test_throws Exception @macroexpand(@mutt begin
        struct A
        end
        struct B
        end
    end)
end
@testset "Simple Macro Usage" begin
    @eval begin
        @mutt struct S
            x
        end
        s = S(1)
        @test is_mutts_mutable(s)
        mark_immutable!(s)
        @test !is_mutts_mutable(s)
        MuttsInterface.make_mutable_copy(s::S) = S(s.x)

        s = S(1)
        @test is_mutts_mutable(branch!(s))
        @test !is_mutts_mutable(s)

        # Assignment from mark_immutable! (PR #4)
        s1 = mark_immutable!(s)
        @test s1 == s
    end
end

@mutt struct M
    x
end
@testset "Can't modify immutable Mutts" begin
    m = M(1)
    m.x = 2
    @test m.x == 2
    mark_immutable!(m)
    @test_throws Exception m.x = 3
end

@testset "Inner constructors" begin
    # (eval required to create structs inside a Testset)
    @eval begin
        # Test Inner Constructors
        @mutt struct SimpleFields
            x
            y
        end
        @mutt struct NoCustomInner
            x :: Int
            y
        end
        @mutt struct WithCustomInners
            x :: Int
            y
            WithCustomInners(x, y=2) = new(x,y)
        end
        function make end
        @mutt struct WithInnerFunctions
            x :: Int
            y
            WithInnerFunctions(x, y=2) = new(x,y)
            # keyword arguments constructor
            WithInnerFunctions(; x=1, y=2) = new(x,y)
            # Custom inner function
            @__MODULE__().make() = new(1,2)
        end
        @mutt struct DefaultsWithTypeParams{A,B}
            x :: A
            y :: B
        end
        @mutt struct CustomWithTypeParams{A,B}
            x :: A
            y :: B
            CustomWithTypeParams() = new{Int,Int}(1,2)
            CustomWithTypeParams{A,B}(x,y) where {A,B} = new{A,B}(x,y)
        end
    end
    @eval begin
        @test is_mutts_mutable(SimpleFields(1,2))
        @test is_mutts_mutable(NoCustomInner(1,2))
        @test is_mutts_mutable(WithCustomInners(1,2))
        @test is_mutts_mutable(WithCustomInners(1))

        @test is_mutts_mutable(WithInnerFunctions(1))
        @test is_mutts_mutable(WithInnerFunctions(y=1))
        @test is_mutts_mutable(make())

        @test is_mutts_mutable(DefaultsWithTypeParams{Int,String}(1,""))
        @test is_mutts_mutable(CustomWithTypeParams())
        @test is_mutts_mutable(CustomWithTypeParams{Int,String}(1,""))
    end
end

@testset "structs with supertypes" begin
    @eval begin
        @mutt struct MuttsFloat <: AbstractFloat
            v :: Float64
        end
        @test MuttsFloat <: AbstractFloat
        v = MuttsFloat(2.0)
        @test v isa AbstractFloat

        v.v = 3.0
        @test v.v == 3.0
        mark_immutable!(v)
        @test_throws Exception v.v = 4.0
    end
end

end