text:
# Preprocessing the parse data
# Designed to read from conllu-format file,

# regular expression handles which column to take
function loadcorpus(file,v::Vocab)
    corpus = Any[]
    s = Sentence(v)
    for line in eachline(file)
        if line == "\n"
            push!(corpus, s)
            s = Sentence(v)
        elseif (m = match(r"^\d+\t(.+?)\t.+?\t(.+?)\t.+?\t.+?\t(.+?)\t(.+?)(:.+)?\t", line)) != nothing
            #                id   word   lem  upos   xpos feat head   deprel
            word = m.captures[1]
            push!(s.word, word)

            postag = get(v.postags, m.captures[2], 0)
            if postag==0
                Base.warn_once("Unknown postags")
            end
            push!(s.postag, postag)

            head = tryparse(Position, m.captures[3])
            head = isnull(head) ? -1 : head.value
            if head==-1
                Base.warn_once("Unknown heads")
            end
            push!(s.head, head)

            deprel = get(v.deprels, m.captures[4], 0)
            if deprel==0
                Base.warn_once("Unknown deprels")
            end
            push!(s.deprel, deprel)
        end
    end
    return corpus
end

function loadcorpus_v6(file,v::Vocab)
    corpus = Any[]
    s = Sentence(v)
    for line in eachline(file)
        if line == ""
            ######### Use that to experiment with the same sentence length #########
            #(length(s) == 15) && (push!(corpus, s))
            #########
            push!(corpus, s)
            s = Sentence(v)
        elseif (m = match(r"^\d+\t(.+?)\t.+?\t(.+?)\t(.+?)\t.+?\t(.+?)\t(.+?)(:.+)?\t", line)) != nothing
            #                id   word   lem  upos   xpos feat head   deprel
            word = m.captures[1]
            push!(s.word, word)

            postag = get(v.postags, m.captures[2], 0)
            if postag==0
                Base.warn_once("Unknown postags")
            end
            push!(s.postag, postag)

	    #index 3 for xpostag
	    xpostag = get(v.xpostags, m.captures[3],0)
	    if xpostag==0
		xpostag=get!(v.xpostags, m.captures[3], length(v.xpostags)+1)
	    end
            push!(s.xpostag, xpostag)

            head = tryparse(Position, m.captures[4])
            head = isnull(head) ? -1 : head.value
            if head==-1
                Base.warn_once("Unknown heads")
            end
            push!(s.head, head)

            deprel = get(v.deprels, m.captures[5], 0)
            if deprel==0
                Base.warn_once("Unknown deprels")
            end
            push!(s.deprel, deprel)
        end
    end
    ######### To assert desired length taken #########
    #@assert mapreduce(length, +, 0, corpus) / 15 == length(corpus)
    #####################
    return corpus
end

# To create vocabulary from pre-trained lstm model
function create_vocab(d)
    Vocab(d["char_vocab"],
          Dict{String, Int}(),
          d["word_vocab"],
          d["sosword"],
          d["eosword"],
          d["unkword"],
          d["sowchar"],
          d["eowchar"],
          d["unkchar"],
          get(d, "postags", UPOSTAG),
          get(d, "deprels", UDEPREL),
          get(d, "xpostags", XPOSTAG)
          )
end

function maptoint(sentences, v::Vocab)
    MAXWORD = 32
    wdict = empty!(v.idict) # it is already empty ?
    cdict = v.cdict
    unkcid = cdict[v.unkchar]
    words = Vector{Int}[]
    sents = Vector{Int}[]

    maxwordlen = 0; maxsentlen = 0;
    for w in (v.sosword, v.eosword)
        wid = get!(wdict, w, 1+length(wdict))
        word = Array(Int, length(w))
        wordi = 0 # to check 2 byte characters
        for c in w
            word[wordi+=1] = get(cdict, c, unkcid)
        end
        (wordi != length(w)) && error("Missing in single word process")
        (wordi > maxwordlen) && (maxwordlen = wordi)
        push!(words, word)
    end

    for s in sentences
        sent = Array(Int, length(s.word))
        senti = 0
        for w in s.word
            ndict = length(wdict)
            wid = get!(wdict, w, 1+ndict)
            sent[senti+=1] = wid
            if wid == 1+ndict
                word = Array(Int, length(w))
                wordi = 0
                for c in w
                    word[wordi+=1] = get(cdict, c, unkcid)
                end
                (wordi != length(w)) && error("Missing in single word process")
                if wordi > MAXWORD; wordi=MAXWORD; word = word[1:wordi]; end;
                (wordi > maxwordlen) && (maxwordlen = wordi)
                push!(words, word)
            end
        end
        @assert(senti == length(s.word))
        (senti > maxsentlen) && (maxsentlen = senti)
        push!(sents, sent)
    end
    @assert(length(wdict) == length(words))
    return words, sents, maxwordlen, maxsentlen
end

# To make ready for character based lstm, data[i] corresponds to ith time step input, to charlstm
function tokenbatch(words, maxlen, sos, eos, pad=eos)
    B = length(words) # batchsize
    T = maxlen + 2
    data = [ Array(Int, B) for t in 1:T ]
    mask = [ Array(Float32, B) for t in 1:T ]
    @inbounds for t in 1:T
        for b in 1:B
            N = length(words[b]) # wordlen
            n = t - T + N + 1 # cursor
            if n < 0
                mask[t][b] = 0
                data[t][b] = pad
            else
                mask[t][b] = 1
                if n == 0
                    data[t][b] = sos
                elseif n <= N
                    data[t][b] = words[b][n]
                elseif n == N+1
                    data[t][b] = eos
                else
                    error()
                end
            end
        end
    end
    return data, mask
end

# Gold-batch for final layer of the bilstm implementation
function goldbatch(sentences, maxlen, wdict, unkwid, pad=unkwid)
    B = length(sentences)
    T = maxlen
    data = [ Array(Int, B) for t in 1:T ]
    mask = [ Array(Float32, B) for t in 1:T ]
    for t in 1:T
        for b in 1:B
            N = length(sentences[b])
            n = t - T + N
            if n <= 0
                mask[t][b] = 0
                data[t][b] = pad
            else
                mask[t][b] = 1
                data[t][b] = get(wdict, sentences[b].word[n], unkwid)
            end
        end
    end
    return data, mask
end

# minibatching for lstm decision module implementation
function lstm_batch(corpus, batchsize; maxlen=64)
    # for now ignore sentences longer than 64 sentences
    metadata = Dict{Int, Array{Any,1}}()
    for item in corpus
        if haskey(metadata, length(item))
            push!(metadata[length(item)], item)
        elseif length(item) <= maxlen
            metadata[length(item)] = [ item ]
        else
            continue
        end
    end
    data = []
    for (k,v) in metadata
        while !isempty(v)
            endv  =  (length(v) >= batchsize ? batchsize : length(v))
            push!(data, v[1:endv])
            deleteat!(v, 1:endv)
        end
        delete!(metadata, k)
    end
    return shuffle!(data)
end


text:
using OrdinaryDiffEq, Test, DiffEqBase, DiffEqCallbacks, StaticArrays

# In-place modification of du (as array)
model(du,u,t,p) = du .= -0.1.*u
u0 = [1.0, 10.0]
tspan = (0.0, 1000.0)
prob = ODEProblem(model, u0, tspan)
sim = solve(prob, Tsit5(), callback=TerminateSteadyState())

@test all(sim(sim.t[end],Val{1}) .< 1e-8)
@test sim.t[end] < tspan[2]

# Out-of-place modification of du (as array)
smodel(u,t,p) = -0.1.*u
su0 = @SVector [1.0, 10.0]
tspan = (0.0, 1000.0)
prob = ODEProblem(smodel, su0, tspan)
sim = solve(prob, Tsit5(), callback=TerminateSteadyState())

@test all(sim(sim.t[end],Val{1}) .< 1e-8)
@test sim.t[end] < tspan[2]

ext: jl


text:
reciprocal_distretization(x, L) = reciprocal_distretization(lattice(x), L)
function reciprocal_distretization(lattice::AbstractLattice, L)
    r1, r2 = reciprocal_vectors(lattice, L)
    reciprocal_distretization(r1, r2, L)
end
function reciprocal_distretization(v1, v2, L)
    pgram = map(((i, j) for i in 0:L-1 for j in 0:L-1)) do (i, j)
        out = [0., 0.]
        temp = [0., 0.]
        min = Inf
        for n in -1:1, m in -1:1
            temp .= v1 * (n + i/L) .+ v2 * (m + j/L)
            if norm(temp) < min
                min = norm(temp)
                out .= temp
            end
        end
        out
    end
    return pgram
end

structure_factor(args...) = fourier_transform(args...)
function fourier_transform(mc::DQMC, key::Symbol, args...)
    fourier_transform(mc, mc[key], args...)
end
function fourier_transform(mc::DQMC, m::DQMCMeasurement, L::Integer, args...)
    qs = reciprocal_distretization(mc, L)
    return qs, fourier_transform(mc, qs, m, args...)
end
function fourier_transform(mc::DQMC, qs::Vector, m::DQMCMeasurement, args...)
    dirs = directions(mc)
    values = mean(m)
    fourier_transform(qs, dirs, values, args...)
end

function fourier_transform(qs::Vector, dirs::Vector, values::Vector)
    @boundscheck length(dirs) == length(qs) == length(values)
    map(qs) do q
        sum(cis(dot(q, v)) * o for (v, o) in zip(dirs, values))
    end
end

function fourier_transform(qs::Vector, dirs::Vector, values::Matrix, weights::Vector)
    @boundscheck begin
        length(weights) ≤ size(values, 2) &&
        length(dirs) == length(qs) &&
        length(dirs) == size(values, 1)
    end
    map(qs) do q
        out = zero(ComplexF64)
        for i in eachindex(dirs)
            temp = zero(ComplexF64)
            for j in eachindex(weights)
                temp += weights[j] * values[i, j]
            end
            out += cis(dot(q, dirs[i])) * temp
        end
        out
    end
end

function fourier_transform(qs::Vector, dirs::Vector, values::Array{T, 3}, weights::Vector) where {T}
    @boundscheck begin
        length(weights) ≤ size(values, 2) &&
        length(dirs) == length(qs) &&
        length(dirs) == size(values, 1)
    end
    map(qs) do q
        out = zero(ComplexF64)
        @turbo for i in eachindex(dirs)
            temp = zero(ComplexF64)
            for j in eachindex(weights), k in eachindex(weights)
                temp += weights[j] * weights[k] * values[i, j, k]
            end
            out += cis(dot(q, dirs[i])) * temp
        end
        out
    end
end

uniform_fourier(mc::DQMC, key::Symbol) = uniform_fourier(mc[key])
uniform_fourier(m::DQMCMeasurement) = sum(mean(m))

apply_symmetry(mc::DQMC, key::Symbol, weights=(1)) = apply_symmetry(mc[key], weights)
apply_symmetry(m::DQMCMeasurement, weights=(1)) = apply_symmetry(mean(m), weights)
function apply_symmetry(data::AbstractArray{T, 3}, weights=(1)) where {T}
    out = zeros(T, size(data, 1))
    @turbo for i in eachindex(out)
        for j in eachindex(weights), k in eachindex(weights)
            out[i] += weights[j] * weights[k] * data[i,j,k]
        end
    end
    out
end
function apply_symmetry(data::AbstractArray{T, 1}, weights=(1)) where {T}
    out = zero(T)
    @turbo for i in eachindex(weights)
        out += weights[i] * data[i]
    end
    out
end
function apply_symmetry(data::AbstractArray{T, 2}, weights=(1)) where {T}
    out = zeros(T, size(data, 1))
    @turbo for i in eachindex(out)
        for j in eachindex(weights)
            out[i] += weights[j] * data[i,j]
        end
    end
    out
end

# cc2superfluid density
superfluid_density(mc, key::Symbol, L = lattice(mc).L) = superfluid_density(mc, mc[key], L)
function superfluid_density(mc, m::DQMCMeasurement, L=lattice(mc).L)
    dirs = directions(mc)
    qx, qy = reciprocal_vectors(lattice(mc), L)
    superfluid_density(mean(m), dirs, qx, qy)
end
function superfluid_density(data::Array{T, 2}, dirs, qx, qy, skip_zero_distance=true) where {T}
    output = ComplexF64(0)
    for i in axes(data, 1)
        for j in 1+skip_zero_distance:size(data, 2)
            output += (cis(dot(qy, dirs[j])) - cis(dot(qx, dirs[j]))) * data[i, j]
        end
    end
    output
end


text:
###
### grow_beg!
###
# can't call grow_beg! on this
# function grow_beg!(x::DynamicAxis{T}, n::Unsigned) where {T} end

function grow_beg!(x, n::Integer)
    n < 0 && throw(ArgumentError("new length must be ≥ 0"))
    return unsafe_grow_beg!(x, n)
end
unsafe_grow_beg!(x::Vector, n) = Base._growbeg!(x, n)
unsafe_grow_beg!(x::MutableRange, n) = setfield!(x, :parent, unsafe_grow_beg(parent(x), n))

###
### grow_end!
###
""" grow_end! """
function grow_end!(x, n::Integer)
    n < 0 && throw(ArgumentError("new length must be ≥ 0"))
    unsafe_grow_end!(x, n)
    return x
end
unsafe_grow_end!(x::Vector, n) = Base._growend!(x, n)
unsafe_grow_end!(x::MutableRange, n) = setfield!(x, :parent, unsafe_grow_end(parent(x), n))

""" grow_end """
function grow_end(x, n::Integer)
    n < 0 && throw(ArgumentError("n must be positive; got $n"))
    return unsafe_grow_end(x, n)
end
function unsafe_grow_end(x::AbstractRange, n::Integer)
    if known_step(x) === 0
        s = step(x)
        return first(x):s:(last(x) + (s * n))
    else
        return first(x):(last(x) + n)
    end
end

###
### grow_to!
###
function grow_to!(x, n::Integer)
    len = length(x)
    if len < n
        unsafe_grow_to!(x, n)
        return x
    elseif len == n
        return x
    else
        throw(ArgumentError("new length must be ≥ than length of collection, got length $(length(x))."))
    end
end

unsafe_grow_to!(x::AbstractRange, n::Integer) = unsafe_grow_end!(x, n - length(x))

function grow_to(x, n::Integer)
    len = length(x)
    if len <= n
        return unsafe_grow_to(x, n)
    else
        throw(ArgumentError("new length must be ≥ than length of collection, got length $(length(x))."))
    end
end
unsafe_grow_to(x::AbstractRange, n::Integer) = unsafe_grow_end(x, n - length(x))

###
### shrink_beg!
###
function shrink_beg!(x, n::Integer)
    n > length(x) && throw(ArgumentError("new length cannot be < 0"))
    return unsafe_shrink_beg!(x, n)
end
unsafe_shrink_beg!(x::Vector, n) = Base._deletebeg!(x, n)
unsafe_shrink_beg!(x::MutableRange, n) = setfield!(x, :parent, unsafe_shrink_beg(parent(x), n))

###
### shrink_end!
###
function shrink_end!(x, n::Integer)
    n < 0 && throw(ArgumentError("new length must be ≥ 0"))
    return unsafe_shrink_end!(x, n)
end
unsafe_shrink_end!(x::Vector, n) = Base._deleteend!(x, n)
unsafe_shrink_end!(x::MutableRange, n) = setfield!(x, :parent, unsafe_shrink_end(parent(x), n))

###
### shrink_to!(x::AbstractRange, n::Integer
###
function shrink_to!(x, n::Integer)
    len = length(x)
    if len > n
        return unsafe_shrink_to!(x, n)
    elseif len == n
        return x
    else
        throw(ArgumentError("new length must be ≤ than length of collection, got length $(length(x))."))
    end
end

unsafe_shrink_to!(x, n) = unsafe_shrink_end!(x, length(x) - n)

function grow_beg(x, n::Integer)
    n < 0 && throw(ArgumentError("n must be positive; got $n"))
    return unsafe_grow_beg(x, n)
end
function unsafe_grow_beg(x::AbstractRange, n::Integer)
    if known_step(x) === 0
        s = step(x)
        return (first(x) - (s * n)):s:last(x)
    else
        return (first(x) - n):last(x)
    end
end
function shrink_beg(x, n::Integer)
    n < 0 && throw(ArgumentError("n must be positive; got $n"))
    return unsafe_shrink_beg(x, n)
end
function unsafe_shrink_beg(x::AbstractRange, n::Integer)
    if known_step(x) === 0
        s = step(x)
        return (first(x) + (s * n)):s:last(x)
    else
        return (first(x) + n):last(x)
    end
end

function shrink_end(x, n::Integer)
    n < 0 && throw(ArgumentError("n must be positive; got $n"))
    return unsafe_shrink_end(x, n)
end

function unsafe_shrink_end(x::AbstractRange, n::Integer)
    if known_step(x) === 0
        s = step(x)
        return first(x):s:(last(x) - (s * n))
    else
        return first(x):(last(x) - n)
    end
end


text:
using Oceananigans.Simulations:
    stop, iteration_limit_exceeded, stop_time_exceeded, wall_time_limit_exceeded,
    TimeStepWizard, update_Δt!

function run_time_step_wizard_tests(arch)
    grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
    model = IncompressibleModel(architecture=arch, grid=grid)

    Δx = grid.Δx
    CFL = 0.45
    u₀ = 7
    Δt = 2.5
    model.velocities.u[1, 1, 1] = u₀

    wizard = TimeStepWizard(cfl=CFL, Δt=Δt, max_change=Inf, min_change=0)
    update_Δt!(wizard, model)
    @test wizard.Δt ≈ CFL * Δx / u₀

    wizard = TimeStepWizard(cfl=CFL, Δt=Δt, max_change=Inf, min_change=0.75)
    update_Δt!(wizard, model)
    @test wizard.Δt ≈ 0.75Δt

    wizard = TimeStepWizard(cfl=CFL, Δt=Δt, max_change=Inf, min_change=0, min_Δt=1.99)
    update_Δt!(wizard, model)
    @test wizard.Δt ≈ 1.99

    model.velocities.u[1, 1, 1] = u₀/100

    wizard = TimeStepWizard(cfl=CFL, Δt=Δt, max_change=1.1, min_change=0)
    update_Δt!(wizard, model)
    @test wizard.Δt ≈ 1.1Δt

    wizard = TimeStepWizard(cfl=CFL, Δt=Δt, max_change=Inf, min_change=0, max_Δt=3.99)
    update_Δt!(wizard, model)
    @test wizard.Δt ≈ 3.99

    grid_stretched = VerticallyStretchedRectilinearGrid(size = (1, 1, 1),
                                                        x = (0, 1),
                                                        y = (0, 1),
                                                        z_faces = z -> z,
                                                        halo = (1, 1, 1),
                                                        architecture=arch)

    model = IncompressibleModel(architecture=arch, grid=grid_stretched)

    Δx = grid_stretched.Δx
    CFL = 0.45
    u₀ = 7
    Δt = 2.5
    model.velocities.u .= u₀

    wizard = TimeStepWizard(cfl=CFL, Δt=Δt, max_change=Inf, min_change=0)
    update_Δt!(wizard, model)
    @test wizard.Δt ≈ CFL * Δx / u₀

    return nothing
end

function run_basic_simulation_tests(arch, Δt)
    grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
    model = IncompressibleModel(architecture=arch, grid=grid)

    model.clock.time = 0.0
    model.clock.iteration = 0

    simulation = Simulation(model, Δt=Δt, stop_iteration=1)

    # Just make sure we can construct a simulation without any errors.
    @test simulation isa Simulation

    @test iteration_limit_exceeded(simulation) == false
    @test stop(simulation) == false

    run!(simulation)

    # Just make sure run! executes without any errors.
    @test simulation isa Simulation

    # Some basic tests
    @test iteration_limit_exceeded(simulation) == true
    @test stop(simulation) == true

    t = Δt isa Number ? 3 : 5
    @test model.clock.time ≈ t
    @test model.clock.iteration == 1
    @test simulation.run_time > 0

    @test stop_time_exceeded(simulation) == false
    simulation.stop_time = 1e-12
    @test stop_time_exceeded(simulation) == true

    @test wall_time_limit_exceeded(simulation) == false
    simulation.wall_time_limit = 1e-12
    @test wall_time_limit_exceeded(simulation) == true

    # Test that simulation stops at `stop_iteration`.
    model = IncompressibleModel(architecture=arch, grid=grid)
    simulation = Simulation(model, Δt=Δt, stop_iteration=3, iteration_interval=88)
    run!(simulation)

    @test simulation.model.clock.iteration == 3

    # Test that simulation stops at `stop_time`.
    model = IncompressibleModel(architecture=arch, grid=grid)
    simulation = Simulation(model, Δt=Δt, stop_time=20.20, iteration_interval=123)
    run!(simulation)

    @test simulation.model.clock.time ≈ 20.20

    return nothing
end

function run_simulation_date_tests(arch, start_time, stop_time, Δt)
    grid = RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))

    clock = Clock(time=start_time)
    model = IncompressibleModel(architecture=arch, grid=grid, clock=clock)

    simulation = Simulation(model, Δt=Δt, stop_time=stop_time)

    @test model.clock.time == start_time
    @test simulation.stop_time == stop_time

    run!(simulation)

    @test model.clock.time == stop_time
    @test simulation.stop_time == stop_time

    return nothing
end

@testset "Time step wizard" begin
    for arch in archs
        @info "Testing time step wizard [$(typeof(arch))]..."
        run_time_step_wizard_tests(arch)
    end
end

@testset "Simulations" begin
    for arch in archs
        @info "Testing simulations [$(typeof(arch))]..."
        for Δt in (3, TimeStepWizard(Δt=5.0))
            run_basic_simulation_tests(arch, Δt)
        end

        @info "Testing simulations with DateTime [$(typeof(arch))]..."
        run_simulation_date_tests(arch, 0.0, 1.0, 0.3)
        run_simulation_date_tests(arch, DateTime(2020), DateTime(2021), 100days)
        run_simulation_date_tests(arch, TimeDate(2020), TimeDate(2021), 100days)
    end
end


text:
### Parsing config
using ArgParse
using YAML

s = ArgParseSettings()
@add_arg_table! s begin
    "--disable-display"
        help = "Use for UNIX server"
        action = :store_true
    "--expr_name"
        help = "Define expr name"
        required = false
    "--is-restart"
        help = "Continue training?"
        action = :store_true
end
parsed_args = parse_args(ARGS, s)

if parsed_args["expr_name"] != nothing
    expr_name = parsed_args["expr_name"]
    is_restart = parsed_args["is-restart"]
else
    runtime = YAML.load_file("./runtime.yaml")
    expr_name = runtime["expr_name"]
    is_restart = Bool(runtime["is_restart"])
end
conf = YAML.load_file("$expr_name/config.yaml")

## Prepare for working env
using OrdinaryDiffEq, Flux, Optim, Random, Plots
using DiffEqSensitivity
using Zygote
using ForwardDiff
using LinearAlgebra, Statistics
using ProgressBars, Printf
using Flux.Optimise: update!, ExpDecay
using Flux.Losses: mae
using Distributions
using StatsBase
using LatinHypercubeSampling
using BSON: @save, @load
using CSV
using DataFrames
using DelimitedFiles
using Plots.PlotMeasures
using Suppressor

cd(dirname(@__DIR__))
ENV["GKSwstype"] = "100"
fig_path = string("./results/", expr_name, "/figs")
ckpt_path = string("./results/", expr_name, "/checkpoint")
config_path = "./results/$expr_name/config.yaml"
pyplot()

if is_restart
    println("Continue to run $expr_name ...\n")
else
    println("Runing $expr_name ...\n")
end

fig_path = string(expr_name, "/figs")
ckpt_path = string(expr_name, "/checkpoint")

if !is_restart
    if ispath(fig_path)
        rm(fig_path, recursive=true)
    end
    if ispath(ckpt_path)
        rm(ckpt_path, recursive=true)
    end
end

if ispath(fig_path) == false
    mkdir(fig_path)
    mkdir(string(fig_path, "/conditions"))
end

if ispath(ckpt_path) == false
    mkdir(ckpt_path)
end

if haskey(conf, "grad_max")
    grad_max = conf["grad_max"]
else
    grad_max = Inf
end

ns = Int64(conf["ns"]); # number of nodes / species
tfinal = Float64(conf["tfinal"]);
ntotal = Int64(conf["ntotal"]);  # number of samples for each perturbation
nplot = Int64(conf["nplot"]);
batch_size = Int64(conf["batch_size"]);  # STEER

n_exp_train = Int64(conf["n_exp_train"]);
n_exp_val = Int64(conf["n_exp_val"]);
n_exp_test = Int64(conf["n_exp_test"]);

n_exp = n_exp_train + n_exp_val + n_exp_test;
noise = Float64(conf["noise"])
opt = ADAMW(Float64(conf["lr"]), (0.9, 0.999), Float64(conf["weight_decay"]));

n_iter_max = Int64(conf["n_iter_max"])
n_plot = Int64(conf["n_plot"]);
n_iter_buffer = Int64(conf["n_iter_buffer"])
n_iter_burnin = Int64(conf["n_iter_burnin"])
n_iter_tol = Int64(conf["n_iter_tol"])
convergence_tol = Float64(conf["convergence_tol"])


text:
using Pkg
Pkg.activate(pwd())

# # Dictionaries

d = Dict("a" => [1, 2, 3], "b" => 1)
d = Dict(:a => [1, 2, 3], :b => 1)

#+

d[:a]
d[:c]

#+

haskey(d, :c)
get(d, :c, 42)
get!(d, :c, 42)
get!(d, :d, ["hello", "world"])

d

#+

delete!(d, :d)
haskey(d, :d)

#+

pop!(d, :c)
haskey(d, :c)

#+

haskey(d, :c)
pop!(d, :c, 444)


text:
using Documenter, Blink

Blink.AtomShell.install()

makedocs(
    modules = [Blink],
    format = Documenter.HTML(),
    sitename = "Blink",
    pages = [
        "index.md",
        "guide.md",
        "Communication" => "communication.md",
        "api.md",
    ],
)

deploydocs(
    repo = "github.com/JuliaGizmos/Blink.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
function compute!(I::IntegralHelper, entry::String)

    if entry == "S"
        compute_S!(I)
    elseif entry == "T"
        compute_T!(I)
    elseif entry == "V"
        compute_V!(I)
    elseif entry == "ERI"
        compute_ERI!(I)
    else
        throw(FermiException("Invalid key for IntegralHelper: $(entry)."))
    end
end

function compute_S!(I::IntegralHelper{T, E, AtomicOrbitals}) where {T<:AbstractFloat, E<:AbstractERI}
    bs = I.orbitals.basisset
    I.cache["S"] = FermiMDArray(overlap(bs, T))
end

function compute_T!(I::IntegralHelper{T, E, AtomicOrbitals}) where {T<:AbstractFloat, E<:AbstractERI}
    bs = I.orbitals.basisset
    I.cache["T"] = FermiMDArray(kinetic(bs, T))
end

function compute_V!(I::IntegralHelper{T, E, AtomicOrbitals}) where {T<:AbstractFloat, E<:AbstractERI}
    bs = I.orbitals.basisset
    I.cache["V"] = FermiMDArray(nuclear(bs, T))
end

function compute_ERI!(I::IntegralHelper{T, E, AtomicOrbitals}) where {T<:AbstractFloat, E<:AbstractDFERI}

    bs = I.orbitals.basisset
    auxbs = I.eri_type.basisset
    J = FermiMDArray(ERI_2e2c(auxbs, T))
    Pqp = FermiMDArray(ERI_2e3c(bs, auxbs, T))
    Jh = similar(J)
    Jh .= Array(real(J^(-1/2)))
    @tensor b[Q,p,q] := Pqp[p,q,P]*Jh[P,Q]
    I.cache["ERI"] = b
end

function compute_ERI!(I::IntegralHelper{T, Chonky, AtomicOrbitals}) where T<:AbstractFloat
    bs = I.orbitals.basisset
    I.cache["ERI"] = FermiMDArray(ERI_2e4c(bs, T))
end

function compute_ERI!(I::IntegralHelper{T, SparseERI, AtomicOrbitals}) where T<:AbstractFloat
    bs = I.orbitals.basisset
    I.cache["ERI"] = Fermi.FermiSparse(sparseERI_2e4c(bs, T)...)
end


text:
abstract type AbstractParam end

"""
    Param

Holds information about one potentially calibrated parameter (array).
Default value must always be set. Determines size of inputs.
Everything else can be either empty or must have the same size.
"""
mutable struct Param{T1 <: Any} <: AbstractParam
    name :: Symbol
    description :: String
    symbol :: String
    value :: T1
    defaultValue :: T1
    "Value bounds"
    lb :: T1
    ub :: T1
    isCalibrated :: Bool
end

## ----------  ParamVector

"""
    ParamVector

Vector containing all of a model's potentially calibrated parameters.
Parameters contain values, not just default values
They are kept in sync with values in object

Intended workflow:
    See `SampleModel`
    Create a model object with parameters as fields
        Otherwise the code gets too cumbersome
        Constructor initializes ParamVector with defaults (or user inputs)
    During calibration
        Each object generates a Dict of calibrated parameters
        Make this into a vector of Floats that can be passed to the optimizer.
        Optimization algorithm changes the floats
        Make floats back into Dict
        Copy back into model objects

Going from a vector of Dicts to a vector of Floats and back:
    `make_guess`
    `set_params_from_guess!`
    These are called on the top level model object

# ToDo: Make the process of going from model -> vector and vice versa more robust.
    Currently, the user has to ensure that the ordering of ParamVectors and model
    objects never changes.
"""
Base.@kwdef mutable struct ParamVector
    "ObjectId of the ModelObject. To ensure that no mismatches occur."
    objId :: ObjectId
    # A Dict would be natural, but it helps to preserve the order of the params
    # pv :: Vector{Param} = Vector{Param{Any}}()
    pv :: OrderedDict{Symbol, AbstractParam} = OrderedDict{Symbol, AbstractParam}()
    "Governs scaling of parameters into guess vectors for optimization"
    pTransform :: ParamTransformation = LinearTransformation{ValueType}()
end

"""
	$(SIGNATURES)

Collection of `ParamVector` for several `ModelObject`s.

Supports iteration like a `Dict`:

```julia
for (objId, pvec) in pvCollection
    @assert objId isa ObjectId
    @assert pvec isa ParamVector
end
```
"""
mutable struct PVectorCollection
    d :: OrderedDict{ObjectId, ParamVector}
end

## ----------  Parameters

"""
	IncreasingVector

Encodes an increasing vector of fixed length. Its values are calibrated.
"""
mutable struct IncreasingVector{T1} <: ModelObject
	objId :: ObjectId
	pvec :: ParamVector
	x0 :: T1
	dxV :: Vector{T1}
end

"""
	$(SIGNATURES)

Retrieve values of an `IncreasingVector`.
"""
values(iv :: IncreasingVector{T1}) where T1 =
	iv.x0 .+ cumsum(vcat(zero(T1), iv.dxV));

values(iv :: IncreasingVector{T1}, idx) where T1 =
    values(iv)[idx];

Base.length(iv :: IncreasingVector) = Base.length(iv.dxV) + 1;

# Displays parameters in levels, not as intercept and increments.
function param_table(iv :: IncreasingVector{T1}, isCalibrated :: Bool) where T1
    if isCalibrated
        # This is where we get the description and symbol from
        p = iv.pvec[2];
        pt = ParamTable(1);
        set_row!(pt, 1, string(p.name), p.symbol, p.description,
            formatted_value(values(iv)));
    else
        pt = nothing;
    end
    return pt
end

"""
	$(SIGNATURES)

Increasing or decreasing vector with bounds.
The special case where the vector is of length 1 is supported.

A `BoundedVector` is typically constructed with an empty `ParamVector`. Then [`set_pvector!`](@ref) is used to initialize the `ParamVector`.
The `ParamVector` contains a single entry which must be named `:dxV`. It sets the values for the eponymous `BoundedVector` field. The `dxV` are typically in [0, 1]. They represent the increments in the vector.

# Example:
```julia
bv = BoundedVector(objId, pvec, true, 5.0, 10.0, [0.5, 0.4]);
values(bv, 1) == 7.5;
values(bv, 2) == 7.5 + 0.4 * (10.0 - 7.5);
```
"""
mutable struct BoundedVector{T1} <: ModelObject
    objId :: ObjectId
    pvec :: ParamVector
    increasing :: Symbol
    # Values are in these bounds
    lb :: T1
    ub :: T1
    # Increments, typically in [0, 1]
    dxV :: Vector{T1}
end

mutable struct BVector{T1} <: AbstractParam
    name :: Symbol
    description :: String
    symbol :: String
    # Increments relative to bounds. In (0, 1)
    dxV :: Vector{T1}
    defaultDxV :: Vector{T1}
    lb :: T1
    ub :: T1
    isCalibrated :: Bool
    increasing :: Symbol
end

# """
# 	$(SIGNATURES)

# Array with some fixed and some calibrated parameters.
# """
# mutable struct CalibratedArraySwitches{T1, N} <: ModelSwitches
#     pvec :: ParamVector
#     defaultValueM :: Array{T1, N}
#     lbM :: Array{T1, N}
#     ubM :: Array{T1, N}
#     isCalM :: Array{Bool, N}
# end

# mutable struct CalibratedArray{T1, N} <: ModelObject
#     objId :: ObjectId
#     switches :: CalibratedArraySwitches{T1, N}
#     calValueV :: Vector{T1}
#     # This gets allocated once and then updated with calibrated
#     # values from calValueV
#     valueM :: Array{T1, N}
# end

mutable struct CalArray{T1, N} <: AbstractParam
    name :: Symbol
    description :: String
    symbol :: String
    value :: Array{T1, N}
    defaultValue :: Array{T1, N}
    lb :: Array{T1, N}
    ub :: Array{T1, N}
    isCalM :: Array{Bool, N}
end

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


text:
"""
    @option struct Blueprint

The general type `Blueprint` contains a list of other blueprints.

# Fields

- `name::String`: the name of this blueprint
- `config_file::Maybe{String}`: the corresponding config file name of this blueprint.
- `meta:Dict{String, Any}`: meta info from config file or other places.
- `blueprints::Vector{Any}`: list of blueprint objects to compile.
"""
@option struct Blueprint
    name::String # name of this blueprint
    config_file::Maybe{String} # can be nothing if no config file
    meta::Dict{String, Any} = Dict{String, Any}()
    blueprints::Vector{Any} = []

    function Blueprint(
        name::String, config_file::Maybe{String},
        meta::Dict{String, Any}, blueprints::Vector{Any}
    )
        new(name, config_file, meta, sort(blueprints, by=priority))
    end
end

Base.show(io::IO, ::MIME"text/plain", x::Blueprint) = GarishPrint.pprint_struct(io, x)

"""
    priority(blueprint)

Return the priority of `blueprint`. Default is `1000`.
"""
priority(blueprint) = 1000 # defaul priority

"""
    @option struct Context

Type for storing the compile context.

# Fields

- `args::Tuple`: the input arguments.
- `kwargs::NamedTuple`: the input keyword arguments.
- `stack::NTuple{N, Blueprint}`: [`Blueprint`](@ref) stack, topmost is the current blueprint under compile.
- `state::Dict{Symbol, Any}`: compile states that each blueprint can pass in.
"""
@option struct Context{N, Args <: Tuple, Kwargs <: NamedTuple}
    # global inputs that is
    # not from config meta
    args::Args = ()
    kwargs::Kwargs = NamedTuple()

    # blueprint stack
    stack::NTuple{N, Blueprint} = ()
    # state
    state::Dict{Symbol, Any} = Dict{Symbol, Any}()
end

function Base.getindex(ctx::Blueprint, name::String)
    haskey(ctx, name) || error("expect $name in $(ctx.config_file)")
    return ctx.meta[name]
end

compile(bp::Blueprint, args...; kwargs...) =
    compile(bp, Context(args, NamedTuple(kwargs), (bp, ), Dict{Symbol, Any}()))

function compile(bp::Blueprint, ctx::Context)
    new_ctx = Context(ctx.args, ctx.kwargs, (bp, ctx.stack...), ctx.state)
    for bp in bp.blueprints
        compile(bp, new_ctx)
    end
    return
end

function has_blueprint(ctx::Context, ::Type{T}) where T
    for bp in ctx.stack, each in bp.blueprints
        each isa T && return true
    end
    return false
end

function get_blueprint(ctx::Context, ::Type{T}) where T
    for bp in ctx.stack, each in bp.blueprints
        each isa T && return each
    end
    return
end

abstract type MustacheBlueprint end
function blueprint_view(::MustacheBlueprint, ctx::Context)
    return Dict{String, Any}()
end

function src_file(bp::MustacheBlueprint, ::Context)
    return joinpath(bp.template.root, bp.template.path, bp.template.file)
end

function dst_file(bp::MustacheBlueprint, ctx::Context)
    return joinpath(ctx.kwargs.path, bp.template.path, bp.template.file)
end

function compile(bp::MustacheBlueprint, ctx::Context)
    file = src_file(bp, ctx)::String
    text = render_from_file(file, blueprint_view(bp, ctx))
    dst = dst_file(bp, ctx)::String
    mkpath(dirname(dst))
    write(dst, text)
    return
end


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

using Markdown
using InteractiveUtils

# ╔═╡ bb427b30-47ae-11eb-28e6-19ba8afc227d
begin
using Flux
images = Flux.Data.MNIST.images();
labels = Flux.Data.MNIST.labels();
end

# ╔═╡ 96433c12-47f1-11eb-1589-bb257b66361c
begin
	using Random
	x=rand(1:60000)
	md"""
	$(images[x])
	# $(labels[x])"""
end

# ╔═╡ 03e5e240-4800-11eb-0dac-254f8ca6e0f8
begin
	using Flux: onehotbatch
	y=onehotbatch(labels, 0:9)
end

# ╔═╡ 13575c8a-49d2-11eb-35d5-678446f5c47e
using Flux:@epochs

# ╔═╡ a8485fde-47fb-11eb-0c26-05f1e3f300e3
foo=images[800]

# ╔═╡ ee91bc14-47fc-11eb-216c-61f1357586eb
reshape(foo,:)

# ╔═╡ c4977514-47fe-11eb-21f1-6da498eb7362
float(reshape(foo,:))

# ╔═╡ a7e7e35a-47ff-11eb-18df-13c13cf2690f
X=hcat(float.(reshape.(images,:))...)

# ╔═╡ a0362f42-49d1-11eb-0b18-cb34ea364cd9
m = Chain(
  Dense(28*28,40, relu),
  Dense(40, 10),
  softmax)

# ╔═╡ ab272c08-49d1-11eb-218e-fbd72b5323fd
begin
	using Flux:crossentropy,onecold,throttle
	loss(X, y) = crossentropy(m(X), y)
	opt = ADAM()
end

# ╔═╡ dae68f32-49d4-11eb-0f5f-b9caf19a35a2
begin
	using Statistics:mean
	accuracy(X, y) = mean(onecold(m(X)) .== onecold(y))
	md"""
	# Training accuracy: $(accuracy(X, y))
	# Test accuracy: $(accuracy(hcat(float.(reshape.(Flux.Data.MNIST.images(:test), :))...), onehotbatch(Flux.Data.MNIST.labels(:test), 0:9)))"""
end

# ╔═╡ d288829a-49d1-11eb-330c-b3ad51d54203
progress = () -> @show(loss(X, y))

# ╔═╡ d9026228-49d1-11eb-31ed-370d7411fa45
@epochs 100 Flux.train!(loss, params(m),[(X,y)], opt, cb = throttle(progress, 10))

# ╔═╡ 86716a10-49d4-11eb-0c5f-c15aa095dd26
begin
	foobar=rand(1:10000)
	md"""
	$(Flux.Data.MNIST.images(:test)[foobar])
	# Prediction: $(onecold(m(hcat(float.(reshape.(Flux.Data.MNIST.images(:test), :))...)[:,foobar])) - 1)
	"""
end

# ╔═╡ Cell order:
# ╠═bb427b30-47ae-11eb-28e6-19ba8afc227d
# ╠═96433c12-47f1-11eb-1589-bb257b66361c
# ╠═a8485fde-47fb-11eb-0c26-05f1e3f300e3
# ╠═ee91bc14-47fc-11eb-216c-61f1357586eb
# ╠═c4977514-47fe-11eb-21f1-6da498eb7362
# ╠═a7e7e35a-47ff-11eb-18df-13c13cf2690f
# ╠═03e5e240-4800-11eb-0dac-254f8ca6e0f8
# ╠═a0362f42-49d1-11eb-0b18-cb34ea364cd9
# ╠═ab272c08-49d1-11eb-218e-fbd72b5323fd
# ╠═d288829a-49d1-11eb-330c-b3ad51d54203
# ╠═13575c8a-49d2-11eb-35d5-678446f5c47e
# ╠═d9026228-49d1-11eb-31ed-370d7411fa45
# ╠═86716a10-49d4-11eb-0c5f-c15aa095dd26
# ╠═dae68f32-49d4-11eb-0f5f-b9caf19a35a2

max_stars_repo_name: punkmath/Misc

max_issues_repo_name: punkmath/Misc

max_forks_repo_name: punkmath/Misc


text:
# push!(LOAD_PATH, "/home/renee/phd/src/ebit-evolution.project/ebit-ode-server/")
include("/home/renee/phd/src/ebit-evolution.project/ebit-ode-server/EbitODEMessages.jl")
include("/home/renee/phd/src/ebit-evolution.project/ebit-ode-server/EbitCloudSpatialExtend.jl")
include("/home/renee/phd/src/ebit-evolution.project/ebit-ode-server/EbitSolver.jl")
include("/home/renee/phd/src/ebit-evolution.project/ebit-ode-server/EbitODEServer.jl")

module EbitScratching

# using Revise
using Revise
using Main.EbitSolver
using Main.EbitODEServer
using Main.EbitODEMessages
using ProtoBuf

function test()
    msg = Any
    open("testdata/simple_test.kairos") do file
        msg = EbitSolver.solve_ode(ProtoBuf.readproto(file, EbitODEMessages.Message()).ode_problem, ::Any -> false)
    end

    open("testdata/simpe_test_result.kairos", "w+") do file
        ProtoBuf.writeproto(file, msg)
    end
end

end


text:
module CellSegmenter

using Images, Distances, LinearAlgebra, Statistics, DataStructures

include("slic.jl")
include("normalizePixels.jl")
include("gaussianKernel.jl")
include("dbscan.jl")
include("spectralClustering.jl")
include("initQuadTree.jl")
include("buildTree.jl")
include("branchTree.jl")

mutable struct Cluster
  l::Float32
  a::Float32
  b::Float32
  y::Int
  x::Int
end

struct Pixel
  x::Int
  y::Int
  Pixel(x,y) = new(x,y)
end

struct Rectangle
  x::Int
  y::Int
  width::Int
  height::Int

  Rectangle(x,y,w,h) = new(x,y,w,h)
end

mutable struct QuadTree
  boundary::Rectangle
  capacity::Float64
  branch::Int
  leaf::Bool
  pixels::Array{Pixel}

  QuadTree(boundary, capacity, branch) = new(boundary, capacity,
    branch, false, [])
end

export slic, normalizePixels, gaussianKernel, dbscan, spectralClustering, buildTree, branchTree, initQuadTree

end

max_stars_repo_name: neuro-myoung/CellSegmenter

max_issues_repo_name: neuro-myoung/CellSegmenter

max_forks_repo_name: neuro-myoung/CellSegmenter


text:
try
    ts = time_ns()
    let
        test = ByteStreamTestHarness("construction", 15)
        execute(test, InputEnded(false))
        execute(test, BufferEmpty(true))
        execute(test, Eof(false))
        execute(test, BytesRead(0))
        execute(test, BytesWritten(0))
        execute(test, RemainingCapacity(15))
        execute(test, BufferSize(0))
    end

    let
        test = ByteStreamTestHarness("construction-end", 15)
        execute(test, EndInput())
        execute(test, InputEnded(true))
        execute(test, BufferEmpty(true))
        execute(test, Eof(true))
        execute(test, BytesRead(0))
        execute(test, BytesWritten(0))
        execute(test, RemainingCapacity(15))
        execute(test, BufferSize(0))
    end
    tc = (Int(time_ns() - ts)) / 1e9
    @info "byte_stream_construction passed......... Time: $tc"
catch e
    println(stderr, e)
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
@export function E_Gd_L1L2(Gd_spins, L1, L2, J1, J2)
    return   Gd_spins[1].v ⋅ ( J2 * L1 + J1 * L2) +
             Gd_spins[2].v ⋅ ( J1 * L1 + J2 * L2) +
             Gd_spins[3].v ⋅ (-J1 * L1 + J2 * L2) +
             Gd_spins[4].v ⋅ (-J2 * L1 + J1 * L2) +
             Gd_spins[5].v ⋅ (-J2 * L1 - J1 * L2) +
             Gd_spins[6].v ⋅ (-J1 * L1 - J2 * L2) +
             Gd_spins[7].v ⋅ ( J1 * L1 - J2 * L2) +
             Gd_spins[8].v ⋅ ( J2 * L1 - J1 * L2)
end

#Gd under magnetic field
@export E_Gd_H(Gd_spins, H, χ) = χ*mapreduce(x ->  x.v⋅H, -, Gd_spins, init=0.0)

#L1 & L2 under magnetic field
@export E_L_H(L, H, χ) = 0.5χ * (L⋅H)^2

#L1 L2 interaction
@export E_L1L2(L1, L2, γ) = 0.5γ * (L1⋅L2)^2

#L1 & L2 interaction with their easy axis
@export E_L_easy(L, easy, K_L) = - 0.5K_L * (L⋅easy)^2

#Gd interaction with their easy axis
@export function E_Gd_easy(Gd_spins, easy, K_Gd )
    out = 0.0
    for i in eachindex(Gd_spins)
        out -= (Gd_spins[i] ⋅ easy[i])^2
    end
    return K_Gd * out
end

@export Landau_F(ϕ, H; α=10, β=1, Hc=1) = α*(H-Hc) * cos(ϕ-π/4)^2 + sin(2ϕ)^2 + β*(H-Hc)*sin(ϕ)^2


text:
using Test;
include("../src/pointInPolyhedron.jl");

@testset "Other Geometric Stuff" begin
    @testset "areCoplanar" begin
        x = [1.0; 0.0; 0.0];
        y = [0.0; 1.0; 0.0];
        z = [0.0; 0.0; 1.0];

        @test areCoplanar(x, x, y, z) == true;
        @test areCoplanar(x, y, y, z) == true;
        @test areCoplanar(x, y, z, z) == true;
        @test areCoplanar(x, 2*x, y, 2*y) == true;
        @test areCoplanar(x, 2*x, 3*x, y) == true;
        @test areCoplanar(x, 2*x, y, x+y) == true;
        @test areCoplanar(x, x+y, y, 2*x+3*y) == true;
        @test areCoplanar(x, z, x+z, 2*x-z) == true;

        @test areCoplanar(x, y, z, x+2*y) == false;
        @test areCoplanar(x, y, z, x-x) == false;
        @test areCoplanar(x, y, x+y, z) == false;
        @test areCoplanar(x, y, x+y, z*10e-7) == false;
        @test areCoplanar(x, y, x-y, z*10e-8) == false;
        @test areCoplanar(x, y, z, (1-10e-8)*z) == false;
        @test areCoplanar(x, y, 2*x, x-z) == false;
        @test areCoplanar(x, 2*x, 2*y, z) == false;
    end
end

@testset "Determine Simplex Solid Angle" begin
    p = [0.0; 0.0; 0.0];
    x = [1.0; 0.0; 0.0];
    y = [0.0; 1.0; 0.0];
    z = [0.0; 0.0; 1.0];

    @testset "Quadrants" begin
        @test isapprox(simplexSolidAngle(p, [+x -x +z])/π, 1.0);
        @test isapprox(simplexSolidAngle(p, [+x -x+y*1e-10 +z])/π, 1.0);
    end

    @testset "Octants" begin

        @test isapprox(simplexSolidAngle(p, [+x +y +z])/π, 0.5);
        @test isapprox(simplexSolidAngle(p, [+x +y -z])/π, 0.5);
        @test isapprox(simplexSolidAngle(p, [+x -y +z])/π, 0.5);
        @test isapprox(simplexSolidAngle(p, [+x -y -z])/π, 0.5);
        @test isapprox(simplexSolidAngle(p, [-x +y +z])/π, 0.5);
        @test isapprox(simplexSolidAngle(p, [-x +y -z])/π, 0.5);
        @test isapprox(simplexSolidAngle(p, [-x -y +z])/π, 0.5);
        @test isapprox(simplexSolidAngle(p, [-x -y -z])/π, 0.5);

    end

    @testset "Hexant" begin
        @testset "xy-Hexant" begin

            @test isapprox(simplexSolidAngle(p, [+x +x+y +z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+x+y +y +z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+y +y-x +z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+y-x -x +z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-x -x-y +z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-x-y -y +z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-y -y+x +z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-y+x +x +z])/π, 0.25);

            @test isapprox(simplexSolidAngle(p, [+x +x+y -z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+x+y +y -z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+y +y-x -z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+y-x -x -z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-x -x-y -z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-x-y -y -z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-y -y+x -z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-y+x +x -z])/π, 0.25);
        end

        @testset "xz-Hexant" begin
            @test isapprox(simplexSolidAngle(p, [+x +y +x+z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+x+z +y +z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+z +y +z-x])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+z-x +y -x])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-x +y -x-z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-x-z +y -z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-z +y -z+x])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-z+x +y +x])/π, 0.25);

            @test isapprox(simplexSolidAngle(p, [+x -y +x+z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+x+z -y +z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+z -y +z-x])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+z-x -y -x])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-x -y -x-z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-x-z -y -z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-z -y -z+x])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [-z+x -y +x])/π, 0.25);
        end

        @testset "yz-Hexant" begin
            @test isapprox(simplexSolidAngle(p, [+x +y +y+z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+x +y+z +z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+x +z +z-y])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+x +z-y -y])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+x -y -y-z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+x -y-z -z])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+x -z -z+y])/π, 0.25);
            @test isapprox(simplexSolidAngle(p, [+x -z+y +y])/π, 0.25);
        end
    end

    @testset "Other Angles" begin
        @test isapprox(simplexSolidAngle(p, [+x +y-x +z])/π, 0.75)
        @test isapprox(simplexSolidAngle(p, [+x +y-x -z])/π, 0.75)
        @test isapprox(simplexSolidAngle(p, [+x -y-x +z])/π, 0.75)
        @test isapprox(simplexSolidAngle(p, [+x -y-x -z])/π, 0.75)
        @test isapprox(simplexSolidAngle(p, [-x +y+x +z])/π, 0.75)
        @test isapprox(simplexSolidAngle(p, [-x +y+x -z])/π, 0.75)
        @test isapprox(simplexSolidAngle(p, [-x -y+x +z])/π, 0.75)
        @test isapprox(simplexSolidAngle(p, [-x -y+x -z])/π, 0.75)
    end
end

@testset "Point of a Cube" begin
    V = [
        0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0;
        0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0;
        0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0
    ]
    EV = [
        [1, 2], [2, 4], [4, 3], [3, 1],
        [5, 6], [6, 8], [8, 7], [7, 5],
        [1, 5], [2, 6], [4, 8], [3, 5]
    ];
    FV = [
        [1, 3, 4, 2], [5, 6, 8, 7],
        [1, 5, 7, 3], [2, 4, 8, 6],
        [1, 2, 6, 5], [3, 7, 8, 4],
    ];
    triFV = [
        [1, 3, 4], [4, 2, 1], [5, 6, 8], [8, 7, 5],
        [1, 5, 7], [7, 3, 1], [2, 4, 8], [8, 6, 2],
        [1, 2, 6], [6, 5, 1], [3, 7, 8], [8, 4, 3],
    ];

    @testset "Internal Points" begin
        @test pointInPolyhedron3D([0.5; 0.5; 0.5], V, EV, triFV) == 1;
        @test pointInPolyhedron3D([0.01; 0.01; 0.01], V, EV, triFV) == 1;
        @test pointInPolyhedron3D([0.01; 0.01; 0.00000001], V, EV, triFV) == 1;
        @test pointInPolyhedron3D([0.9999999; 0.5; 10^-5], V, EV, triFV) == 1;
    end

    @testset "Outer Points" begin
        @test pointInPolyhedron3D([0.01; 0.01; -0.0000001], V, EV, triFV) == -1;
        @test pointInPolyhedron3D([0.5; 0.5; 1.5], V, EV, triFV) == -1;
        @test pointInPolyhedron3D([0.5; 0.5; 1.0000000001], V, EV, triFV) == -1;
        @test pointInPolyhedron3D([0.0; 0.0; -0.00000001], V, EV, triFV) == -1;
    end

    @testset "Edge Points" begin
        @test pointInPolyhedron3D([0.5; 0.5; 0.0], V, EV, triFV) == 0;
        @test pointInPolyhedron3D([0.5; 0.0; 0.0], V, EV, triFV) == 0;
        @test pointInPolyhedron3D([0.0; 0.0; 0.0], V, EV, triFV) == 0;
        @test pointInPolyhedron3D([1.0; 0.5; 0.5], V, EV, triFV) == 0;
    end
end

@testset "Determine Planar Angles" begin
    @testset "2D Single Angle" begin
        @test isapprox(rot2Dangle([+1.0;  0.0]), 0.00 * π);
        @test isapprox(rot2Dangle([+1.0; +1.0]), 0.25 * π);
        @test isapprox(rot2Dangle([ 0.0; +1.0]), 0.50 * π);
        @test isapprox(rot2Dangle([-1.0; +1.0]), 0.75 * π);
        @test isapprox(rot2Dangle([-1.0;  0.0]), 1.00 * π);
        @test isapprox(rot2Dangle([-1.0; -1.0]), 1.25 * π);
        @test isapprox(rot2Dangle([ 0.0; -1.0]), 1.50 * π);
        @test isapprox(rot2Dangle([+1.0; -1.0]), 1.75 * π);
    end

    @testset "2D Angle" begin
        @test isapprox(rot2Dangle([1.;0.], [+1.; 0.]), 0.00 * π);
        @test isapprox(rot2Dangle([1.;0.], [+1.;+1.]), 0.25 * π);
        @test isapprox(rot2Dangle([1.;0.], [ 0.;+1.]), 0.50 * π);
        @test isapprox(rot2Dangle([1.;0.], [-1.;+1.]), 0.75 * π);
        @test isapprox(rot2Dangle([1.;0.], [-1.;+0.]), 1.00 * π);
        @test isapprox(rot2Dangle([1.;0.], [-1.;-1.]), 1.25 * π);
        @test isapprox(rot2Dangle([1.;0.], [ 0.;-1.]), 1.50 * π);
        @test isapprox(rot2Dangle([1.;0.], [+1.;-1.]), 1.75 * π);

        @test isapprox(rot2Dangle([1.;1.], [+1.; 1.]), 0.00 * π);
        @test isapprox(rot2Dangle([1.;1.], [ 0.;+1.]), 0.25 * π);
        @test isapprox(rot2Dangle([1.;1.], [-1.;+1.]), 0.50 * π);
        @test isapprox(rot2Dangle([1.;1.], [-1.; 0.]), 0.75 * π);
        @test isapprox(rot2Dangle([1.;1.], [-1.;-1.]), 1.00 * π);
        @test isapprox(rot2Dangle([1.;1.], [ 0.;-1.]), 1.25 * π);
        @test isapprox(rot2Dangle([1.;1.], [+1.;-1.]), 1.50 * π);
        @test isapprox(rot2Dangle([1.;1.], [+1.; 0.]), 1.75 * π);

        @test isapprox(rot2Dangle([+1.;-1.], [+1.;-1.]), 0.00 * π);
        @test isapprox(rot2Dangle([+1.;-1.], [+1.; 0.]), 0.25 * π);
        @test isapprox(rot2Dangle([+1.;-1.], [+1.;+1.]), 0.50 * π);
        @test isapprox(rot2Dangle([+1.;-1.], [ 0.;+1.]), 0.75 * π);
        @test isapprox(rot2Dangle([+1.;-1.], [-1.;+1.]), 1.00 * π);
        @test isapprox(rot2Dangle([+1.;-1.], [-1.; 0.]), 1.25 * π);
        @test isapprox(rot2Dangle([+1.;-1.], [-1.;-1.]), 1.50 * π);
        @test isapprox(rot2Dangle([+1.;-1.], [ 0.;-1.]), 1.75 * π);

        @test isapprox(rot2Dangle([-1.;+0.], [-1.; 0.]), 0.00 * π);
        @test isapprox(rot2Dangle([-1.;+0.], [-1.;-1.]), 0.25 * π);
        @test isapprox(rot2Dangle([-1.;+0.], [ 0.;-1.]), 0.50 * π);
        @test isapprox(rot2Dangle([-1.;+0.], [+1.;-1.]), 0.75 * π);
        @test isapprox(rot2Dangle([-1.;+0.], [+1.; 0.]), 1.00 * π);
        @test isapprox(rot2Dangle([-1.;+0.], [+1.;+1.]), 1.25 * π);
        @test isapprox(rot2Dangle([-1.;+0.], [ 0.;+1.]), 1.50 * π);
        @test isapprox(rot2Dangle([-1.;+0.], [-1.;+1.]), 1.75 * π);
    end
end

@testset "Point in a Square" begin
    V = [
        0.0 1.0 1.0 0.0;
        0.0 0.0 1.0 1.0
    ];
    EV = [ [1, 2], [2, 3], [3, 4], [4, 1] ];

    @test pointInPolyhedron2D([0.5;  0.5], V, EV) == +1;
    @test pointInPolyhedron2D([0.1;  0.1], V, EV) == +1;
    @test pointInPolyhedron2D([0.1;  0.9], V, EV) == +1;
    @test pointInPolyhedron2D([0.5; 1e-7], V, EV) == +1;

    @test pointInPolyhedron2D([0.5;  1.5], V, EV) == -1;
    @test pointInPolyhedron2D([1.0;  1.5], V, EV) == -1;
    @test pointInPolyhedron2D([0.5;-1e-7], V, EV) == -1;
    @test pointInPolyhedron2D([1.5;  0.5], V, EV) == -1;

    @test pointInPolyhedron2D([1.0;  0.5], V, EV) ==  0;
    @test pointInPolyhedron2D([1.0;  1.0], V, EV) ==  0;
    @test pointInPolyhedron2D([0.5; 1e-9], V, EV) ==  0;
    @test pointInPolyhedron2D([0.5;-1e-9], V, EV) ==  0;
end

@testset "Point in a Octagon" begin
    V = [
        +1.0 +0.6  0.0 -0.6 -1.0 -0.6  0.0 +0.6;
         0.0 +0.6 +1.0 +0.6  0.0 -0.6 -1.0 -0.6
    ];
    EV = [
        [1, 2], [2, 3], [3, 4], [4, 5],
        [5, 6], [6, 7], [7, 8], [8, 1]
    ];

    @test pointInPolyhedron2D([0.5;     0.5], V, EV) == +1;
    @test pointInPolyhedron2D([0.1;     0.1], V, EV) == +1;
    @test pointInPolyhedron2D([0.0;     0.9], V, EV) == +1;
    @test pointInPolyhedron2D([0.0; 1-1e-16], V, EV) == +1;

    @test pointInPolyhedron2D([0.5;     1.5], V, EV) == -1;
    @test pointInPolyhedron2D([1.0;     1.5], V, EV) == -1;
    @test pointInPolyhedron2D([1.0;  -1e-23], V, EV) == -1;
    @test pointInPolyhedron2D([1.4;    -0.6], V, EV) == -1;

    @test pointInPolyhedron2D([0.6;     0.6], V, EV) ==  0;
    @test pointInPolyhedron2D([1.0;     0.0], V, EV) ==  0;
    @test pointInPolyhedron2D([0.0; 1-1e-17], V, EV) ==  0;
    @test pointInPolyhedron2D([0.8;     0.3], V, EV) ==  0;
end


text:
module SDLGamesLib

using SimpleDirectMediaLayer
SDL2 = SimpleDirectMediaLayer

include("objects.jl")
include("display.jl")
include("timing.jl")
include("window.jl")
include("events.jl")
include("anim.jl")
#include("pixel.jl")

end # module


text:
#===============
    STRUCTS
===============#

mutable struct RTMaterial
    color::RayTracer.rt_color
    ambient::Float64
    diffuse::Float64
    specular::Float64
    shininess::Float64
end

#===============
    FUNCTIONS
===============#

function material(
        color::RayTracer.rt_color,ambient::Float64,
        diffuse::Float64,specular::Float64,shininess::Float64
    )
    return RTMaterial(color,ambient,diffuse,specular,shininess)
end

function material()
    return RTMaterial(RayTracer.color(1,1,1),0.1,0.9,0.9,200)
end

function Base.isapprox(a::RTMaterial,b::RTMaterial,atol::Real=0)
    return (≈(a.color,b.color)) && (≈(a.ambient,b.ambient,atol=atol)) &&
        (≈(a.diffuse,b.diffuse,atol=atol)) && (≈(a.specular,b.specular,atol=atol)) &&
        (≈(a.shininess,b.shininess,atol=atol))
end

function lighting(
        material::RTMaterial,light::RayTracer.RTLight,p::RayTracer.rt_tuple,
        eye_v::RayTracer.rt_tuple,normal_v::RayTracer.rt_tuple
    )
    # Combining light and material color
    effective_color = material.color * light.intensity

    # Direction of light source
    light_v = RayTracer.normalize(light.position - p)

    # Ambient contribution to final color
    ambient = effective_color*material.ambient

    light_normal_component = RayTracer.dot(light_v,normal_v)
    if light_normal_component<0
        diffuse = RayTracer.color(0,0,0)
        specular = RayTracer.color(0,0,0)
    else
        # Diffuse contribution
        diffuse = effective_color*material.diffuse*light_normal_component

        # Specular contribution
        light_v_reflected = RayTracer.reflect(-light_v,normal_v)
        reflection_eye_component = RayTracer.dot(light_v_reflected,eye_v)

        if reflection_eye_component<0
            specular = RayTracer.color(0,0,0)
        else
            factor = reflection_eye_component^material.shininess
            specular = light.intensity*material.specular*factor
        end
    end

    return ambient + diffuse + specular
end