text:
function geneigsolve(f, x₀, howmany::Int, which::Selector, alg::GolubYe)
    krylovdim = alg.krylovdim
    maxiter = alg.maxiter
    howmany > krylovdim && error("krylov dimension $(krylovdim) too small to compute $howmany eigenvalues")

    ## FIRST ITERATION: setting up
    numiter = 1
    ax₀, bx₀ = apply(geneigfun(f), x₀)
    numops = 1
    β₀ = norm(x₀)
    iszero(β₀) && throw(ArgumentError("initial vector should not have norm zero"))
    xax = dot(x₀, ax₀)/β₀^2
    xbx = dot(x₀, bx₀)/β₀^2
    T = promote_type(typeof(xax), typeof(xbx))
    invβ₀ = one(T)/β₀
    v = invβ₀ * x₀ # v = mul!(similar(x₀, T), x₀, invβ₀)
    av = mul!(similar(v), ax₀, invβ₀)
    bv = mul!(similar(v), bx₀, invβ₀)
    ρ = checkhermitian(xax) / checkposdef(xbx)
    r = axpy!(-ρ, bv, av)
    tol::typeof(ρ) = alg.tol

    # allocate storage
    HHA = fill(zero(T), krylovdim+1, krylovdim+1)
    HHB = fill(zero(T), krylovdim+1, krylovdim+1)

    # Start Lanczos iteration with A - ρ ⋅ B
    V = OrthonormalBasis([v])
    BV = [bv]
    sizehint!(V, krylovdim+1)
    sizehint!(BV, krylovdim+1)

    r, α = orthogonalize!(r, v, alg.orth) # α should be zero, otherwise ρ was miscalculated
    β = norm(r)
    m = 1
    HHA[m,m] = real(α)

    while m < krylovdim
        v = rmul!(r, 1/β)
        push!(V, r)
        HHA[m+1,m] = β
        HHA[m, m+1] = β
        βold = β
        r, α, β, bv = golubyerecurrence(f, ρ, V, βold, alg.orth)
        numops += 1
        m += 1
        n = hypot(α, β, βold)
        HHA[m,m] = checkhermitian(α, n)
        push!(BV, bv)

        if alg.verbosity > 2
            @info "Golub Ye iteration step $m: normres = $β"
        end
        β < tol && m >= howmany && break
    end

    # Process
    HA = view(HHA, 1:m, 1:m)
    HB = view(HHB, 1:m, 1:m)
    buildHB!(HB, V, BV)
    HA .+= ρ .* HB

    D, Z = geneigh!(HA, HB)
    by, rev = eigsort(which)
    p = sortperm(D, by = by, rev = rev)
    xold = V[1]

    z = view(Z, :, p[1])
    x = mul!(similar(xold), V, z)
    ax, bx = apply(geneigfun(f), x)
    numops += 1
    ρ = checkhermitian(dot(x, ax)) / checkposdef(dot(x, bx))
    r = axpy!(-ρ, bx, ax)
    normr = norm(r)

    converged = 0
    values = Vector{typeof(ρ)}(undef, 0)
    vectors = Vector{typeof(x)}(undef, 0)
    residuals = Vector{typeof(r)}(undef, 0)
    normresiduals = Vector{typeof(normr)}(undef, 0)
    while normr <= tol && converged < howmany
        push!(values, ρ)
        push!(vectors, x)
        push!(residuals, r)
        push!(normresiduals, normr)
        converged += 1

        z = view(Z, :, p[converged+1])
        x = mul!(similar(xold), V, z)
        ax, bx = apply(geneigfun(f), x)
        numops += 1
        ρ = checkhermitian(dot(x, ax)) / checkposdef(dot(x, bx))
        r = axpy!(-ρ, bx, ax)
        normr = norm(r)
    end

    if alg.verbosity > 1
        msg = "Golub-Ye generalized eigsolve in iter $numiter: "
        msg *= "$converged values converged: $ρ, normres = ("
        for i = 1:converged
            msg *= @sprintf("%.2e", normresiduals[i])
            msg *= ", "
        end
        if converged < howmany
            msg *= @sprintf("%.2e", normr)
        end
        msg *= ")"
        @info msg
    end

    ## OTHER ITERATIONS: recycle
    while numiter < maxiter && converged < howmany
        numiter += 1
        fill!(HHA, zero(T))
        fill!(HHB, zero(T))
        resize!(V, 0)
        resize!(BV, 0)

        # Start Lanczos iteration with A - ρ ⋅ B
        invβ = 1/norm(x)
        v = rmul!(x, invβ)
        bv = rmul!(bx, invβ)
        r = rmul!(r, invβ)
        r, α = orthogonalize!(r, v, alg.orth) # α should be zero, otherwise ρ was miscalculated
        β = norm(r)
        m = 1
        push!(V, v)
        HHA[m,m] = real(α)
        push!(BV, bv)

        while m < krylovdim - converged
            v = rmul!(r, 1/β)
            push!(V, r)
            HHA[m+1, m] = β
            HHA[m, m+1] = β
            βold = β
            r, α, β, bv = golubyerecurrence(f, ρ, V, βold, alg.orth)
            numops += 1
            m += 1
            n = hypot(α, β, βold)
            HHA[m,m] = checkhermitian(α, n)
            push!(BV, bv)

            if alg.verbosity > 2
                @info "Golub Ye iteration step $m: normres = $β"
            end
            β < tol && m >= howmany && break
        end

        # add xold
        v, = orthonormalize!(xold, V, alg.orth)
        av, bv = apply(geneigfun(f), v)
        numops += 1
        av = axpy!(-ρ, bv, av)
        for i = 1:m
            HHA[i, m+1] = dot(V[i], av)
            HHA[m+1, i] = conj(HHA[i, m+1])
        end
        m += 1
        HHA[m, m] = checkhermitian(dot(v, av))
        push!(V, v)
        push!(BV, bv)

        # add converged vectors
        @inbounds for i = 1:converged
            v, = orthonormalize(vectors[i], V, alg.orth)
            av, bv = apply(geneigfun(f), v)
            numops += 1
            av = axpy!(-ρ, bv, av)
            for i = 1:m
                HHA[i, m+1] = dot(V[i], av)
                HHA[m+1, i] = conj(HHA[i, m+1])
            end
            m += 1
            HHA[m, m] = checkhermitian(dot(v, av))
            push!(V, v)
            push!(BV, bv)
        end

        # Process
        HA = view(HHA, 1:m, 1:m)
        HB = view(HHB, 1:m, 1:m)
        buildHB!(HB, V, BV)
        HA .+= ρ .* HB

        D, Z = geneigh!(HA, HB)
        by, rev = eigsort(which)
        p = sortperm(D, by = by, rev = rev)
        xold = V[1]

        converged = 0
        resize!(values, 0)
        resize!(vectors, 0)
        resize!(residuals, 0)
        resize!(normresiduals, 0)
        while true
            z = view(Z, :, p[converged+1])
            x = mul!(similar(xold), V, z)
            ax, bx = apply(geneigfun(f), x)
            numops += 1
            ρ = checkhermitian(dot(x, ax)) / checkposdef(dot(x, bx))
            r = axpy!(-ρ, bx, ax)
            normr = norm(r)

            if normr > tol || converged >= howmany
                break
            end

            push!(values, ρ)
            push!(vectors, x)
            push!(residuals, r)
            push!(normresiduals, normr)
            converged += 1
        end

        if alg.verbosity > 1
            msg = "Golub-Ye generalized eigsolve in iter $numiter: "
            msg *= "$converged values converged: $ρ, normres = ("
            for i = 1:converged
                msg *= @sprintf("%.2e", normresiduals[i])
                msg *= ", "
            end
            if converged < howmany
                msg *= @sprintf("%.2e", normr)
            end
            msg *= ")"
            @info msg
        end
    end

    if converged > howmany
        howmany = converged
    end
    for k = converged+1:howmany
        z = view(Z, :, p[k])
        x = mul!(similar(xold), V, z)
        ax, bx = apply(geneigfun(f), x)
        numops += 1
        ρ = checkhermitian(dot(x, ax)) / checkposdef(dot(x, bx))
        r = axpy!(-ρ, bx, ax)
        normr = norm(r)

        push!(values, ρ)
        push!(vectors, x)
        push!(residuals, r)
        push!(normresiduals, normr)
    end

    if alg.verbosity > 0
        if converged < howmany
            @warn """GolubYe eigsolve finished without convergence after $numiter iterations:
             *  $converged eigenvalues converged
             *  norm of residuals = $((normresiduals...,))
             *  number of operations = $numops"""
        else
            @info """Lanczos eigsolve finished after $numiter iterations:
             *  $converged eigenvalues converged
             *  norm of residuals = $((normresiduals...,))
             *  number of operations = $numops"""
        end
    end

    return values, vectors, ConvergenceInfo(converged, residuals, normresiduals, numiter, numops)
end

function golubyerecurrence(f, ρ, V::OrthonormalBasis, β, orth::ClassicalGramSchmidt)
    v = V[end]
    av, bv = apply(geneigfun(f), v)
    w = axpy!(-ρ, bv, av)
    α = dot(v, w)

    w = axpy!(-β, V[end-1], w)
    w = axpy!(-α, v, w)
    β = norm(w)
    return w, α, β, bv
end
function golubyerecurrence(f, ρ, V::OrthonormalBasis, β, orth::ModifiedGramSchmidt)
    v = V[end]
    av, bv = apply(geneigfun(f), v)
    w = axpy!(-ρ, bv, av)
    w = axpy!(-β, V[end-1], w)

    w, α = orthogonalize!(w, v, orth)
    β = norm(w)
    return w, α, β, bv
end
function golubyerecurrence(f, ρ, V::OrthonormalBasis, β, orth::ClassicalGramSchmidt2)
    v = V[end]
    av, bv = apply(geneigfun(f), v)
    w = axpy!(-ρ, bv, av)
    α = dot(v, w)
    w = axpy!(-β, V[end-1], w)
    w = axpy!(-α, v, w)

    w, s = orthogonalize!(w, V, ClassicalGramSchmidt())
    α += s[end]
    β = norm(w)
    return w, α, β, bv
end
function golubyerecurrence(f, ρ, V::OrthonormalBasis, β, orth::ModifiedGramSchmidt2)
    v = V[end]
    av, bv = apply(geneigfun(f), v)
    w = axpy!(-ρ, bv, av)
    w = axpy!(-β, V[end-1], w)
    w, α = orthogonalize!(w, v, ModifiedGramSchmidt())

    s = α
    for q in V
        w, s = orthogonalize!(w, q, ModifiedGramSchmidt())
    end
    α += s
    β = norm(w)
    return w, α, β, bv
end
function golubyerecurrence(f, ρ, V::OrthonormalBasis, β, orth::ClassicalGramSchmidtIR)
    v = V[end]
    av, bv = apply(geneigfun(f), v)
    w = axpy!(-ρ, bv, av)
    α = dot(v, w)
    w = axpy!(-β, V[end-1], w)
    w = axpy!(-α, v, w)

    ab2 = abs2(α) + abs2(β)
    β = norm(w)
    nold = sqrt(abs2(β)+ab2)
    while eps(one(β)) < β < orth.η * nold
        nold = β
        w, s = orthogonalize!(w, V, ClassicalGramSchmidt())
        α += s[end]
        β = norm(w)
    end
    return w, α, β, bv
end
function golubyerecurrence(f, ρ, V::OrthonormalBasis, β, orth::ModifiedGramSchmidtIR)
    v = V[end]
    av, bv = apply(geneigfun(f), v)
    w = axpy!(-ρ, bv, av)
    w = axpy!(-β, V[end-1], w)

    w, α = orthogonalize!(w, v, ModifiedGramSchmidt())
    ab2 = abs2(α) + abs2(β)
    β = norm(w)
    nold = sqrt(abs2(β)+ab2)
    while eps(one(β)) < β < orth.η * nold
        nold = β
        s = zero(α)
        for q in V
            w, s = orthogonalize!(w, q, ModifiedGramSchmidt())
        end
        α += s
        β = norm(w)
    end
    return w, α, β, bv
end

function buildHB!(HB, V, BV)
    m = length(V)
    @inbounds for j = 1:m
        HB[j,j] = checkposdef(dot(V[j], BV[j]))
        for i = j+1:m
            HB[i,j] = dot(V[i], BV[j])
            HB[j,i] = conj(HB[i,j])
        end
    end
end

geneigfun((A,B)::Tuple{Any,Any}) = x->(apply(A,x), apply(B,x))
geneigfun(f) = f


text:
module Mapper

using ..Model
using SQLite, DBInterface, Strapping

const DB = Ref{SQLite.DB}()
getdb() = DB[]

function __init__()
    DB[] = SQLite.DB()
    DBInterface.execute(getdb(), """
        CREATE TABLE album (
            id INTEGER PRIMARY KEY,
            name TEXT,
            artist TEXT,
            year INTEGER,
            timespicked INTEGER DEFAULT 0
        )
    """)
    DBInterface.execute(getdb(), """
        CREATE TABLE songs (
            album_id INTEGER,
            name TEXT
        )
    """)
    DBInterface.execute(getdb(), """
        CREATE INDEX idx_album_id ON songs (album_id)
    """)
    return
end

function create!(album::Album)
    stmt = DBInterface.@prepare(getdb, """
        INSERT INTO album (name, artist, year) VALUES(?, ?, ?)
    """)
    cursor = DBInterface.execute(stmt, (album.name, album.artist, album.year))
    id = DBInterface.lastrowid(cursor)
    album.id = id
    DBInterface.executemany(DBInterface.@prepare(getdb, """
        INSERT INTO songs (album_id, name) VALUES (?, ?)
    """), ([id for _ = 1:length(album.songs)], album.songs))
    return
end

function update(album)
    DBInterface.execute(DBInterface.@prepare(getdb, """
        UPDATE album
        SET name = ?,
            artist = ?,
            year = ?,
            timespicked = ?
        WHERE id = ?
    """), (album.name, album.artist, album.year, album.timespicked, album.id))
    DBInterface.execute(DBInterface.@prepare(getdb, """
        DELETE FROM songs WHERE album_id = ?
    """), (album.id,))
    DBInterface.executemany(DBInterface.@prepare(getdb, """
        INSERT INTO songs (album_id, name) VALUES (?, ?)
    """), ([album.id for _ = 1:length(album.songs)], album.songs))
    return
end

function get(id)
    Strapping.construct(Album, DBInterface.execute(DBInterface.@prepare(getdb, """
        SELECT A.id, A.name, A.artist, A.year, A.timespicked, B.name as songs FROM album A
        LEFT JOIN songs B ON A.id = B.album_id
        WHERE id = ?
    """), (id,)))
end

function delete(id)
    DBInterface.execute(DBInterface.@prepare(getdb, """
        DELETE FROM album WHERE id = ?
    """), (id,))
    DBInterface.execute(DBInterface.@prepare(getdb, """
        DELETE FROM songs WHERE album_id = ?
    """), (id,))
    return
end

function getAllAlbums()
    Strapping.construct(Vector{Album}, DBInterface.execute(DBInterface.@prepare(getdb, """
        SELECT A.id, A.name, A.artist, A.year, A.timespicked, B.name as songs FROM album A
        LEFT JOIN songs B ON A.id = B.album_id
    """)))
end

end # module

ext: jl

max_stars_repo_name: vanhp/MicroServiceJulia

max_issues_repo_name: vanhp/MicroServiceJulia

max_forks_repo_name: vanhp/MicroServiceJulia


text:
using Distributions

const DISTRIBUTIONS = [Degenerate, Gaussian, Bernoulli]

# The role of these tests is to ensure that the
# the general functionality always works and
# doesn't result in NaN or Inf weights.
# Not that the RBM performance is still good.
if :integration in TEST_GROUPS
    @testset "RBM Integration" begin
        n_vis = 100
        n_hid = 10

        for T in [Float32, Float64]
            rbm = RBM(T, Gaussian, Bernoulli, n_vis, n_hid)
            Boltzmann.test(rbm; debug=true)
        end
    end
end

if :benchmark in TEST_GROUPS
    @testset "RBM Benchmark" begin
        n_vis = 100
        n_hid = 10

        suite = BenchmarkGroup()

        for T in [Float32, Float64]
            rbm = RBM(T, Gaussian, Bernoulli, n_vis, n_hid)
            suite = benchmark!(rbm, suite; debug=true)
        end

        tune!(suite)
        results = run(suite, verbose=true, seconds=10)
    end
end


text:
x = CertainValue(2.0)

test_constraints = [
    NoConstraint(),
    TruncateLowerQuantile(0.2),
    TruncateUpperQuantile(0.2),
    TruncateQuantiles(0.2, 0.8),
    TruncateMaximum(50),
    TruncateMinimum(-50),
    TruncateRange(-50, 50),
    TruncateStd(1)
]

for constraint in test_constraints
    @test constrain(x, constraint) isa CertainValue
end


text:
######################
# 4: Properties
######################

@doc Markdown.doc"""
    isnormal( v::AbstractNormalToricVariety )

Checks if the normal toric variety `v` is normal. (This function is somewhat tautological at this point.)

# Examples
```jldoctest
julia> isnormal(projective_space( 2 ))
true
```
"""
function isnormal( v::AbstractNormalToricVariety )
    return true
end
export isnormal

@doc Markdown.doc"""
    isaffine( v::AbstractNormalToricVariety )

Checks if the normal toric variety `v` is affine.

# Examples
```jldoctest
julia> isaffine( projective_space( 2 ) )
false
```
"""
function isaffine( v::AbstractNormalToricVariety )
    return pm_ntv(v).AFFINE::Bool
end
export isaffine

@doc Markdown.doc"""
    isprojective( v::AbstractNormalToricVariety )

Checks if the normal toric variety `v` is projective, i.e. if the fan of `v` is the the normal fan of a polytope.

# Examples
```jldoctest
julia> isprojective( projective_space( 2 ) )
true
```
"""
function isprojective( v::AbstractNormalToricVariety )
    return pm_ntv(v).PROJECTIVE::Bool
end
export isprojective

@doc Markdown.doc"""
    issmooth( v::AbstractNormalToricVariety )

Checks if the normal toric variety `v` is smooth.

# Examples
```jldoctest
julia> issmooth( projective_space( 2 ) )
true
```
"""
function issmooth( v::AbstractNormalToricVariety )
    return pm_ntv(v).SMOOTH::Bool
end
export issmooth

@doc Markdown.doc"""
    iscomplete( v::AbstractNormalToricVariety )

Checks if the normal toric variety `v` is complete.

# Examples
```jldoctest
julia> iscomplete( projective_space( 2 ) )
true
```
"""
function iscomplete( v::AbstractNormalToricVariety )
    return pm_ntv(v).COMPLETE::Bool
end
export iscomplete

@doc Markdown.doc"""
    has_torusfactor( v::AbstractNormalToricVariety )

Checks if the normal toric variety `v` has a torus factor.

# Examples
```jldoctest
julia> has_torusfactor( projective_space( 2 ) )
false
```
"""
function has_torusfactor( v::AbstractNormalToricVariety )
    return v.polymakeNTV.FAN_DIM < v.polymakeNTV.FAN_AMBIENT_DIM
end
export has_torusfactor

@doc Markdown.doc"""
    is_orbifold( v::AbstractNormalToricVariety )

Checks if the normal toric variety `v` is an orbifold.

# Examples
```jldoctest
julia> is_orbifold( projective_space( 2 ) )
true
```
"""
function is_orbifold( v::AbstractNormalToricVariety )
    return pm_ntv(v).SIMPLICIAL::Bool
end
export is_orbifold

@doc Markdown.doc"""
    issimplicial( v::AbstractNormalToricVariety )

Checks if the normal toric variety `v` is simplicial. Hence, this function works just as `is_orbifold`. It is implemented for user convenience.

# Examples
```jldoctest
julia> issimplicial( projective_space( 2 ) )
true
```
"""
function issimplicial( v::AbstractNormalToricVariety )
    return pm_ntv(v).SIMPLICIAL::Bool
end
export issimplicial

@doc Markdown.doc"""
    is_gorenstein( v::AbstractNormalToricVariety )

Checks if the normal toric variety `v` is Gorenstein.

# Examples
```jldoctest
julia> is_gorenstein( projective_space( 2 ) )
true
```
"""
function is_gorenstein( v::AbstractNormalToricVariety )
    return pm_ntv(v).GORENSTEIN::Bool
end
export is_gorenstein

@doc Markdown.doc"""
    is_q_gorenstein( v::AbstractNormalToricVariety )

Checks if the normal toric variety `v` is Q-Gorenstein.

# Examples
```jldoctest
julia> is_q_gorenstein( projective_space( 2 ) )
true
```
"""
function is_q_gorenstein( v::AbstractNormalToricVariety )
    return pm_ntv(v).Q_GORENSTEIN::Bool
end
export is_q_gorenstein

@doc Markdown.doc"""
    is_fano( v::AbstractNormalToricVariety )

Checks if the normal toric variety `v` is fano.

# Examples
```jldoctest
julia> is_fano( projective_space( 2 ) )
true
```
"""
function is_fano( v::AbstractNormalToricVariety )
    return pm_ntv(v).FANO::Bool
end
export is_fano


text:
mutable struct AxelrodAgent <: Agents.AbstractAgent
    id::Int
    pos::NTuple{2, Int}
    stubborn::Bool
    culture::AbstractArray
    changed_culture::Bool
end

function run_random(;
    config_name::Union{String, Nothing}=nothing,
    batch_name::String="default",
    grid_dims::Tuple=(10, 10),
    stubborn_count::Int=4,
    levels::Int=2,
    steps::Int=100,
    replicates::Int=1,
    when::Int=10,
    rnd_seed::Int=0,
    write::Bool=false
)
    if !("data" in readdir())
        mkdir("data")
    end
    if !(batch_name in readdir("data"))
        mkdir(joinpath("data", batch_name))
    end
    Random.seed!(rnd_seed)
    agent_df_list = DataFrames.DataFrame[]
    config_list = Any[]
    p = ProgressMeter.Progress(replicates, "Running Simulations ...")
    for i in 1:replicates
        config = [(rand(1:grid_dims[1]), rand(1:grid_dims[2])) for j in 1:stubborn_count]
        push!(config_list, deepcopy(config))
        model = initialize_model(grid_dims, config, levels)
        rep_df, _ = Agents.run!(
            model, agent_step!, steps, adata=[:pos, :culture, :changed_culture],
            replicates=1, when=0:when:steps, parallel=true, obtainer=deepcopy
        )
        rep_df[!, :replicate] .= i
        push!(agent_df_list, deepcopy(rep_df))
        ProgressMeter.next!(p)
    end
    agent_df = reduce(vcat, agent_df_list)
    prepare_data!(agent_df, batch_name)
    if write
        if config_name === nothing
            Feather.write(joinpath("data", batch_name, batch_name * ".feather"), agent_df)
        else
            Feather.write(joinpath("data", batch_name, config_name * ".feather"), agent_df)
        end

    end
    config_table = create_config_table(
        config_list=config_list,
        grid_dims=grid_dims,
        batch_name=batch_name,
        write=write
    )
    return agent_df, config_table
end

function run_batch(;
    config_list::Union{AbstractArray, Dict},
    batch_name::String="default",
    grid_dims::Tuple=(10, 10),
    levels::Int=2,
    steps::Int=100,
    replicates::Int=1,
    when::Int=10,
    write::Bool=false
)
    if !("data" in readdir())
        mkdir("data")
    end
    if !(batch_name in readdir("data"))
        mkdir(joinpath("data", batch_name))
    end
    results = DataFrames.DataFrame[]
    p = ProgressMeter.Progress(length(config_list), "Running Simulations ...")
    for (i, config) in enumerate(config_list)  # REWRITE FOR DICT!!!!! <- TO DO
        tmp = run_config(
            config_name=batch_name * "_" * string(i),
            batch_name=batch_name,
            grid_dims=grid_dims,
            stubborn_positions=config,
            levels=levels,
            steps=steps,
            replicates=replicates,
            when=when,
            write=false
        )
        if write
            Feather.write(
                joinpath("data", batch_name, batch_name * "_" * string(i) * ".feather"),
                tmp
            )
        end
        push!(results, deepcopy(tmp))
        ProgressMeter.next!(p)
    end
    config_table = create_config_table(
        config_list=config_list,
        grid_dims=grid_dims,
        batch_name=batch_name,
        write=write
    )
    return results, config_table
end

function run_config(;
    config_name::String,
    batch_name::String,
    grid_dims::Tuple,
    stubborn_positions::Union{AbstractArray, Dict},
    levels::Int,
    steps::Int,
    replicates::Int,
    when::Int,
    write::Bool
)
    if !("data" in readdir())
        mkdir("data")
    end
    if !(batch_name in readdir("data"))
        mkdir(joinpath("data", batch_name))
    end
    model = initialize_model(grid_dims, stubborn_positions, levels)
    agent_df, _ = Agents.run!(
        model, agent_step!, steps, adata=[:pos, :culture, :changed_culture],
        replicates=replicates, when=0:when:steps, parallel=true, obtainer=deepcopy
    )
    prepare_data!(agent_df, config_name)
    if write
        Feather.write(joinpath("data", batch_name, config_name * ".feather"), agent_df)
    end
    return agent_df
end

function initialize_model(dims::NTuple{2, Int}, stubborn_positions::Union{AbstractArray, Dict}, levels::Int)
    space = Agents.GridSpace(dims, periodic=false, moore=false)
    model = Agents.AgentBasedModel(AxelrodAgent, space, scheduler=random_activation)
    populate!(model, dims, levels)
    to_stubborn!(stubborn_positions, model)
    return model
end

function populate!(model::Agents.AgentBasedModel, dims::NTuple{2, Int}, levels::Int)
    positions = [(i, j) for i in 1:dims[1] for j in 1:dims[2]]
    for (id, pos) in enumerate(positions)
        Agents.add_agent_pos!(
            AxelrodAgent(id, pos, false, rand(0:levels - 1, 5), 0),
            model
        )
    end
    return model
end

function to_stubborn!(positions::Array{NTuple{2, Int}}, model::Agents.AgentBasedModel, value::Int64=0)
    for pos in positions
        stubborn_agent = Agents.get_node_agents(Agents.coord2vertex(pos, model), model)[1]
        stubborn_agent.stubborn = true
        stubborn_agent.culture = fill(value, 5)
    end
    return model
end

function to_stubborn!(positions::Dict, model::Agents.AgentBasedModel)
    for key in keys(positions)
        to_stubborn!(positions[key], model, key)
    end
    return model
end

function agent_step!(agent::Agents.AbstractAgent, model::Agents.AgentBasedModel)
    neighbors = Agents.node_neighbors(agent, model)
    interaction_partner_pos = StatsBase.sample(neighbors)
    interaction_partner = Agents.get_node_agents(
        Agents.coord2vertex(interaction_partner_pos, model), model
    )[1]
    similarity = StatsBase.mean(agent.culture .== interaction_partner.culture)
    if !(similarity == 1.0) & !agent.stubborn & (rand() <= similarity)
        assimilate!(agent, interaction_partner)
        agent.changed_culture = 1
    else
        agent.changed_culture = 0
    end
    return agent
end

function assimilate!(agent::Agents.AbstractAgent, interaction_partner::Agents.AbstractAgent)
    random_attr = rand(1:length(agent.culture))
    if !(agent.culture[random_attr] == interaction_partner.culture[random_attr])
        agent.culture[random_attr] = interaction_partner.culture[random_attr]
    else
        assimilate!(agent, interaction_partner)
    end
    return agent
end

function prepare_data!(dataframe::DataFrames.DataFrame, config_name::String)
    dataframe[!, "culture"] = [join(c) for c in dataframe[!, "culture"]]
    dataframe[!, :x] = [i[1] for i in dataframe[!, :pos]]
    dataframe[!, :y] = [i[2] for i in dataframe[!, :pos]]
    DataFrames.select!(dataframe, DataFrames.Not(:pos))
    dataframe[!, :config] .= config_name
    return dataframe
end

function create_config_table(;
    config_list::AbstractArray,
    grid_dims::Tuple=(10, 10),
    batch_name::String="default",
    write::Bool=false
)
    if !("data" in readdir())
        mkdir("data")
    end
    if !(batch_name in readdir("data"))
        mkdir(joinpath("data", batch_name))
    end
    configs = DataFrames.DataFrame([
        (i, j)
        for i in 1:grid_dims[1]
        for j in 1:grid_dims[2]
    ])
    DataFrames.names!(configs, [:x, :y])
    for (i, config) in enumerate(config_list)
        col_name = Symbol(batch_name * "_" * string(i))
        configs[!, col_name] .= 0
        for conf in config
            for row in eachrow(configs)
                if row[:x] == conf[1] && row[:y] == conf[2]
                    row[col_name] = 1
                end
            end
        end
    end
    if write
        Feather.write(joinpath("data", batch_name, "configs.feather"), configs)
    end
    return configs
end


text:
# Convenience aliasing + function maps
const RWTypes = Union{Float64,Int16,Int32,UInt16,UInt32}
const Frequency = Union{Float64,Integer}

const read_analog_fun = IdDict(
    Float64 => DAQmx.ReadAnalogF64,
    Int16   => DAQmx.ReadBinaryI16,
    Int32   => DAQmx.ReadBinaryI32,
    UInt16  => DAQmx.ReadBinaryU16,
    UInt32  => DAQmx.ReadBinaryU32)

const write_analog_fun = IdDict(
    Float64 => DAQmx.WriteAnalogF64,
    Int16   => DAQmx.WriteBinaryI16,
    Int32   => DAQmx.WriteBinaryI32,
    UInt16  => DAQmx.WriteBinaryU16,
    UInt32  => DAQmx.WriteBinaryU32)

const AI_SIGNAL_EVENTS = (DAQmx.SampleClock,
                          DAQmx.SampleCompleteEvent,
                          DAQmx.ChangeDetectionEvent,
                          DAQmx.CounterOutputEvent)

function readalloc(task::DAQTask{AI},
                   samples::Int = 1024,
                   ::Type{T} = Float64) where T <: RWTypes
    return Array{T}(undef, length(task.channels), samples)
end

function Base.read!(handle::TaskHandle, samples_perchan::Integer, dst::AbstractArray{T}) where T <: RWTypes
    samples_read = Ref{Int32}()
    buffer_samples = length(dst)
    read_analog_fun[T](handle, samples_perchan, 1.0,
                       DAQmx.GroupByScanNumber, dst,
                       buffer_samples, samples_read) |> catch_error
    return nothing
end

Base.read!(dst::AbstractArray{<:RWTypes}, task::DAQTask{AI}) = read!(task.handle, div(length(dst),length(task.channels)), dst)

function Base.read(task::DAQTask{AI}, samples::Int = 1024,
                   precision::Type{T} = Float64) where T <: RWTypes
    data = readalloc(task, samples, precision)
    read!(data, task)
    return data
end

struct DAQEventCB
    uvhandle::Ptr{Cvoid}
    task_handle::TaskHandle
    event_type::Cint
    num_samples::Cuint
end

struct DAQDoneCB
    uvhandle::Ptr{Cvoid}
    task_handle::TaskHandle
    status::Cint
end

function event_notify(task_handle::TaskHandle, event_type::Cint, num_samples::Cuint,
                      data::Ptr{Nothing})::Cint
    ptr = convert(Ptr{DAQEventCB}, data)
    uvhandle = unsafe_load(ptr).uvhandle
    val = DAQEventCB(uvhandle, task_handle, event_type, num_samples)
    unsafe_store!(ptr, val)
    ccall(:uv_async_send, Nothing, (Ptr{Nothing},), uvhandle)
    return 0
end

function done_notify(task_handle::TaskHandle, status::Cint, data::Ptr{Nothing})::Cint
    if status !== Cint(0)
        ptr = convert(Ptr{DAQDoneCB}, data)
        uvhandle = unsafe_load(ptr).uvhandle
        val = DAQDoneCB(uvhandle, task_handle, status)
        unsafe_store!(ptr, val)
        ccall(:uv_async_send, Nothing, (Ptr{Nothing},), uvhandle)
    end
    return 0
end

function clearevents!(task::DAQTask{AI})
    DAQmx.RegisterEveryNSamplesEvent(task.handle, DAQmx.Acquired_Into_Buffer,
                                     0, C_NULL, C_NULL) |> catch_error
    DAQmx.RegisterDoneEvent(task.handle, C_NULL, C_NULL) |> catch_error
    for event in AI_SIGNAL_EVENTS
        DAQmx.RegisterSignalEvent(task.handle, event, C_NULL, C_NULL) |> catch_error
    end
    return nothing
end

function start(task::DAQTask{AI}, callback::Function, samples_perchan::Integer, args...)

    num_channels = length(task.channels)

    nsamplescb = Base.AsyncCondition()
    donecb = Base.AsyncCondition()
    push!(task.conditions, nsamplescb)
    push!(task.conditions, donecb)

    data_ref = Ref(DAQEventCB(Base.unsafe_convert(Ptr{Cvoid}, nsamplescb), C_NULL, 0, 0))
    status_ref = Ref(DAQDoneCB(Base.unsafe_convert(Ptr{Cvoid}, donecb), C_NULL, 0))

    filename = Dates.format(Dates.now(), "yyyy_mm_dd_HHMMSS") * "data.h5"
    logfile = h5open(filename, "w")
    A = create_dataset(logfile, "A", Float64, ((samples_perchan*num_channels,),(-1,)), chunk=(samples_perchan * num_channels,))

    GC.@preserve nsamplescb donecb task data_ref status_ref logfile A begin

    clearevents!(task) # remove any stale callbacks
    DAQmx.RegisterEveryNSamplesEvent(task.handle, DAQmx.Acquired_Into_Buffer,
                                     samples_perchan, event_notify_c[], data_ref) |> catch_error
    DAQmx.RegisterDoneEvent(task.handle, done_notify_c[], status_ref) |> catch_error

    @async begin
        try
            start(task)
            @async while isrunning(task)
                Base.isopen(donecb) && Base.wait(donecb)
                catch_error(status_ref[].status)
                GC.safepoint()
            end
            while isrunning(task)
                Base.isopen(nsamplescb) && Base.wait(nsamplescb)
                if isrunning(task) # maybe excessively paranoid...
                    data = data_ref[]
                    callback(data.task_handle, data.num_samples, A, args...)
                end
                GC.safepoint()
            end
        catch
            rethrow()
        finally
            flush(logfile)
            close(logfile)
        end
    end
    end # GC preserve
    return nothing
end

function read_to_feed(handle::TaskHandle, num_samples::Integer, A, dst::AbstractArray, feed::Union{RemoteChannel,Channel})
    read!(handle, num_samples, dst)
    dims = HDF5.get_extent_dims(A)[1]
    new_dims = (dims[1]+length(dst),)
    HDF5.set_extent_dims(A, new_dims)
    A[(dims[1]+1):new_dims[1]] = dst
    put!(feed, dst)
    return
end

function record!(task::DAQTask{AI},
                 sampling::Frequency=20000,
                 refresh::Frequency=50;
                 feed::Union{RemoteChannel,Channel})
    num_samples_perchan::Int64 = cld(sampling, refresh)
    buffer = readalloc(task, num_samples_perchan)
    num_channels = length(task.channels)
    DAQmx.CfgSampClkTiming(task.handle, "", sampling, DAQmx.Rising,
                          DAQmx.ContSamps, num_samples_perchan) |> catch_error
    start(task, read_to_feed, num_samples_perchan, buffer, feed)
end

# Write functions

function Base.write(handle::TaskHandle, samples_perchan::Integer, wave::AbstractArray{T}) where {T <: RWTypes}

    samples_written = Ref{Int32}()
    write_analog_fun[T](handle, samples_perchan, false,
                        1.0, DAQmx.GroupByScanNumber, wave,
                        samples_written) |> catch_error

    return samples_written[]
end


text:
"""
    Network(...)
"""

mutable struct Network <: Field

end


text:
"""
    solve_op_problem!(op_problem::OperationsProblem; kwargs...)

This solves the operational model for a single instance and
outputs results of type OperationsProblemResult: objective value, time log,
a dictionary of variables and their dataframe of results, and a time stamp.

# Arguments
-`op_problem::OperationModel = op_problem`: operation model

# Examples
```julia
results = solve_op_problem!(OpModel)
```
# Accepted Key Words
-`save_path::String`: If a file path is provided the results
automatically get written to feather files
-`optimizer::OptimizerFactory`: The optimizer that is used to solve the model
"""
function solve_op_problem!(op_problem::OperationsProblem; kwargs...)

    timed_log = Dict{Symbol, Any}()

    save_path = get(kwargs, :save_path, nothing)

    if op_problem.canonical.JuMPmodel.moi_backend.state == MOIU.NO_OPTIMIZER

        if !(:optimizer in keys(kwargs))
            error("No Optimizer has been defined, can't solve the operational problem")
        end

        _, timed_log[:timed_solve_time],
        timed_log[:solve_bytes_alloc],
        timed_log[:sec_in_gc] = @timed JuMP.optimize!(op_problem.canonical.JuMPmodel,
                                                        kwargs[:optimizer])

    else

        _, timed_log[:timed_solve_time],
        timed_log[:solve_bytes_alloc],
        timed_log[:sec_in_gc] = @timed JuMP.optimize!(op_problem.canonical.JuMPmodel)

    end

    vars_result = get_model_result(op_problem)
    optimizer_log = get_optimizer_log(op_problem)
    time_stamp = get_time_stamps(op_problem)
    time_stamp = shorten_time_stamp(time_stamp)
    obj_value = Dict(:OBJECTIVE_FUNCTION => JuMP.objective_value(op_problem.canonical.JuMPmodel))
    merge!(optimizer_log, timed_log)

    results = make_results(vars_result, obj_value, optimizer_log, time_stamp)

    !isnothing(save_path) && write_results(results, save_path)

     return results

end

function _run_stage(stage::_Stage, start_time::Dates.DateTime, results_path::String; kwargs...)

    if stage.canonical.JuMPmodel.moi_backend.state == MOIU.NO_OPTIMIZER
        error("No Optimizer has been defined, can't solve the operational problem stage with key $(stage.key)")
    end
    timed_log = Dict{Symbol, Any}()
    _, timed_log[:timed_solve_time],
    timed_log[:solve_bytes_alloc],
    timed_log[:sec_in_gc] = @timed JuMP.optimize!(stage.canonical.JuMPmodel)
    model_status = JuMP.primal_status(stage.canonical.JuMPmodel)
    if model_status != MOI.FEASIBLE_POINT::MOI.ResultStatusCode
        error("Stage $(stage.key) status is $(model_status)")
    end
    retrieve_duals = get(kwargs, :duals, nothing)
    if !isnothing(retrieve_duals) && !isnothing(get_constraints(stage.canonical))
        _export_model_result(stage, start_time, results_path, retrieve_duals)
    else
        _export_model_result(stage, start_time, results_path)
    end
    _export_optimizer_log(timed_log, stage.canonical, results_path)
    stage.execution_count += 1
    return
end

"""
    execute!(sim::Simulation; verbose::Bool = false, kwargs...)

Solves the simulation model for sequential Simulations
and populates a nested folder structure created in Simulation()
with a dated folder of featherfiles that contain the results for
each stage and step.

# Arguments
- `sim::Simulation=sim`: simulation object created by Simulation()

# Example
```julia
sim = Simulation("test", 7, stages, "/Users/lhanig/Downloads/";
verbose = true, system_to_file = false)
execute!!(sim::Simulation; verbose::Bool = false, kwargs...)
```

# Accepted Key Words
`dual_constraints::Vector{Symbol}`: if dual variables are desired in the
results, include a vector of the variable names to be included
"""

function execute!(sim::Simulation; verbose::Bool = false, kwargs...)
    _prepare_workspace!(sim.ref, sim.base_name, sim.simulation_folder)
    if sim.ref.reset
        sim.ref.reset = false
    elseif sim.ref.reset == false
        error("Reset the simulation")
    end
    variable_names = Dict()
    steps = get_steps(sim)
    for s in 1:steps
        verbose && println("Step $(s)")
        for (ix, stage) in enumerate(sim.stages)
            verbose && println("Stage $(ix)")
            interval = stage.interval
            for run in 1:stage.executions
                sim.ref.current_time = sim.ref.date_ref[ix]
                verbose && println("Simulation TimeStamp: $(sim.ref.current_time)")
                raw_results_path = joinpath(sim.ref.raw,"step-$(s)-stage-$(ix)",replace_chars("$(sim.ref.current_time)",":","-"))
                mkpath(raw_results_path)
                update_stage!(stage, s, sim)
                dual_constraints = get(kwargs, :dual_constraints, nothing)
                _run_stage(stage, sim.ref.current_time, raw_results_path; duals = dual_constraints)
                sim.ref.run_count[s][ix] += 1
                sim.ref.date_ref[ix] = sim.ref.date_ref[ix] + interval
            end
            @assert stage.executions == stage.execution_count
            stage.execution_count = 0 # reset stage execution_count
        end

    end
    date_run = convert(String,last(split(dirname(sim.ref.raw),"/")))
    ref = make_references(sim, date_run)
    return ref
end


text:
module PMSudoku

# package code goes here

end # module


text:
A = sparse([2,4,3],[2,2,4],[0.2,0.4,0.6])
Af = Matrix(A)

let k
k = 0
for I in eachindex(stored(A))
    @test A[I] == A.nzval[k+=1]
end

k = 0
for j = inds(A, 2)
    for I in eachindex(stored(A, :, j))
        @test A[I] == A.nzval[k+=1]
    end
end

k = 0
for v in each(stored(A))
    @test v == A.nzval[k+=1]
end

k = 0
for j = inds(A, 2)
    for v in each(stored(A, :, j))
        @test v == A.nzval[k+=1]
    end
end

k = 0
for I in each(index(A))
    @test A[I] == Af[k+=1]
end

k = 0
for j = inds(A, 2)
    for I in each(index(A, :, j))
        @test A[I] == Af[k+=1]
    end
end

k = 0
for v in each(A)
    @test v == Af[k+=1]
end

k = 0
for j = inds(A, 2)
    for v in each(A, :, j)
        @test v == Af[k+=1]
    end
end
end

# Sparse matrix-vector multiplication
function matvecmul_ind!(b::AbstractVector, A::AbstractMatrix, x::AbstractVector)
    fill!(b, 0)
    inds(A, 2) == inds(x, 1) || throw(DimensionMismatch("inds(A, 2) = $(inds(A, 2)) does not agree with inds(x, 1) = $(inds(x, 1))"))
    for j in inds(A, 2)
        xj = x[j]
        for (ib, iA) in sync(Follower(index(b)), index(stored(A, :, j)))
            b[ib] += A[iA]*xj
        end
    end
    b
end
function matvecmul_val!(b::AbstractVector, A::AbstractMatrix, x::AbstractVector)
    fill!(b, 0)
    inds(A, 2) == inds(x, 1) || throw(DimensionMismatch("inds(A, 2) = $(inds(A, 2)) does not agree with inds(x, 1) = $(inds(x, 1))"))
    for j in inds(A, 2)
        xj = x[j]
        for (ib, a) in sync(Follower(index(b)), stored(A, :, j))
            b[ib] += a*xj
        end
    end
    b
end

x = [1,-5,7,-13]
btrue = A*x
b = similar(btrue)
matvecmul_ind!(b, A, x)
@test b ≈ btrue
matvecmul_val!(b, A, x)
@test b ≈ btrue
@test_throws DimensionMismatch matvecmul_ind!(b, A, [1])
@test_throws DimensionMismatch matvecmul_ind!([0.1,0.2], A, x)
@test_throws DimensionMismatch matvecmul_val!(b, A, [1])
@test_throws DimensionMismatch matvecmul_val!([0.1,0.2], A, x)


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ e7925da0-08b0-11eb-10d8-e37b1e342978
begin
	using Revise
	using VolterraGP
	using Plots
	using Distributions
	using LinearAlgebra
end

# ╔═╡ 8b6f1218-08b3-11eb-2db7-13b175c26316
import Pkg; Pkg.add(Pkg.PackageSpec(path="/Users/magnus/Documents/phd/code/repos/VolterraGP"))

# ╔═╡ 20e87940-08b1-11eb-3618-a92d2228cb3a
md"
Volterra Kernels
================

An iteractive notebook demonstrating the properties of Volterra kernels
"

# ╔═╡ ab3de578-08b4-11eb-1fe9-2f18a33760a0
function k(τ, τp, C, Gpars, upars; mout=false)
	pars = DiffableParameters([0.1], Gpars, upars)
	if mout
		gp = GaussianProcess(threeEQs, 2, C, 1, Data([], [[]]), pars)
		return VolterraGP.kernel(τp, τ, 1, 2, gp)
	else
		gp = GaussianProcess(threeEQs, 1, C, 1, Data([], [[]]), pars)
		return VolterraGP.kernel(τp, τ, 1, 1, gp)
	end
end

# ╔═╡ 71a180c4-08b7-11eb-09cd-adf1fdfa7f14
md"
C=2 case, so 3 smoothing kerenels and one base kerenel, all of EQ form.
"

# ╔═╡ adcce676-08b5-11eb-393e-f3e2ed46246a
@bind g1 html"<input type='range' step='0.01' min='0' max='3'>"

# ╔═╡ 997768e4-08b6-11eb-1131-9b5ab9222922
@bind g2 html"<input type='range' step='0.01' min='0' max='3'>"

# ╔═╡ 9977aba6-08b6-11eb-16bd-fd79436e020a
@bind g3 html"<input type='range' step='0.01' min='0' max='3'>"

# ╔═╡ 9977f3e0-08b6-11eb-08f1-f1c92afa2158
@bind u html"<input type='range' step='0.01' min='0' max='3'>"

# ╔═╡ 40b9b7fe-08b4-11eb-37bc-8ff804c302da
begin
	τ = collect(-10:0.1:10)
	C = 2
	GparsC2 = reshape([g1, g2, g3], (1, 3, 1))
	uparsC2 = [u]
	yC2 = [k(t, 0., C, GparsC2, uparsC2) for t in τ]
	plot(τ, yC2)
end

# ╔═╡ 301084ee-0d48-11eb-0fa6-058f6c71c2fe
begin
	K = reduce(hcat, [[k(t, tp, C, GparsC2, uparsC2) for t in τ] for tp in τ])
	dist = MvNormal(K + 1e-5*I)
	samp = rand(dist)
	plot(samp)
end

# ╔═╡ 6be3f44e-08b7-11eb-3827-fbbf2fff84b8
md"
Adjustable C with with homogeneous smoothing kernels, for 1 output
"

# ╔═╡ d20b4a9e-08b7-11eb-3093-9572d4fdbeca
@bind Cv html"<input type='range' step='1' min='1' max='5'>"

# ╔═╡ d34dcca4-08b7-11eb-38b7-21de4471ae2b
@bind g html"<input type='range' step='0.01' min='0' max='3'>"

# ╔═╡ d4d813c4-08b7-11eb-2fd9-ad141b799b1a
@bind uCv html"<input type='range' step='0.01' min='0' max='3'>"

# ╔═╡ e65c08f8-08b7-11eb-1680-b535c915caf0
begin
	GparsCv = fill(g, (1, sum(1:Cv), 1))
	uparsCv = [uCv]
	yCv = [k(t, 0., Cv, GparsCv, uparsCv) for t in τ]
	plot(τ, yCv)
end

# ╔═╡ 04e644ca-08bb-11eb-1c44-79d43cef287f
md"
Adjustable C with with homogeneous smoothing kernels, for different ouputs
"

# ╔═╡ 8e408b6a-08b9-11eb-046a-cfe34a818fb5
@bind Cvd html"<input type='range' step='1' min='1' max='5'>"

# ╔═╡ 9d98bb00-08b9-11eb-26de-159d482901cd
@bind gd1 html"<input type='range' step='0.01' min='0' max='3'>"

# ╔═╡ a5dd7256-08b9-11eb-257a-7b599ac29688
@bind gd2 html"<input type='range' step='0.01' min='0' max='3'>"

# ╔═╡ add40150-08b9-11eb-3b17-17b9a10668e3
@bind uCvd html"<input type='range' step='0.01' min='0' max='3'>"

# ╔═╡ c348ddb2-08b9-11eb-0d70-57bb6e917392
begin
	GparsCvd = [fill(gd1, (1, sum(1:Cvd), 1)) ; fill(gd2, (1, sum(1:Cvd), 1))]
	uparsCvd= [uCvd]
	yCvd = [k(t, 0., Cvd, GparsCvd, uparsCvd; mout=true) for t in τ]
	plot(τ, yCvd)
end

# ╔═╡ Cell order:
# ╟─8b6f1218-08b3-11eb-2db7-13b175c26316
# ╟─e7925da0-08b0-11eb-10d8-e37b1e342978
# ╟─20e87940-08b1-11eb-3618-a92d2228cb3a
# ╟─ab3de578-08b4-11eb-1fe9-2f18a33760a0
# ╟─71a180c4-08b7-11eb-09cd-adf1fdfa7f14
# ╠═adcce676-08b5-11eb-393e-f3e2ed46246a
# ╠═997768e4-08b6-11eb-1131-9b5ab9222922
# ╠═9977aba6-08b6-11eb-16bd-fd79436e020a
# ╠═9977f3e0-08b6-11eb-08f1-f1c92afa2158
# ╟─40b9b7fe-08b4-11eb-37bc-8ff804c302da
# ╟─301084ee-0d48-11eb-0fa6-058f6c71c2fe
# ╠═6be3f44e-08b7-11eb-3827-fbbf2fff84b8
# ╠═d20b4a9e-08b7-11eb-3093-9572d4fdbeca
# ╠═d34dcca4-08b7-11eb-38b7-21de4471ae2b
# ╠═d4d813c4-08b7-11eb-2fd9-ad141b799b1a
# ╠═e65c08f8-08b7-11eb-1680-b535c915caf0
# ╟─04e644ca-08bb-11eb-1c44-79d43cef287f
# ╠═8e408b6a-08b9-11eb-046a-cfe34a818fb5
# ╠═9d98bb00-08b9-11eb-26de-159d482901cd
# ╠═a5dd7256-08b9-11eb-257a-7b599ac29688
# ╠═add40150-08b9-11eb-3b17-17b9a10668e3
# ╠═c348ddb2-08b9-11eb-0d70-57bb6e917392


text: e=:(println("e=:($(string(e)));eval(e)"));eval(e)

max_stars_repo_name: MakeNowJust/quine

max_issues_repo_name: makenowjust/quine

max_forks_repo_name: makenowjust/quine


text:
using TDD
using Test

@testset "TDD.jl" begin
    # Write your tests here.
    include("interval.jl")
end


text:
using ClimateMachine, MPI
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Topologies
using ClimateMachine.MPIStateArrays

import ClimateMachine.Mesh.Grids: DiscontinuousSpectralElementGrid
import ClimateMachine.Mesh.Topologies: StackedBrickTopology
export Interval, Periodic

import LinearAlgebra: ×
import Base: *, getindex

abstract type AbstractDomain end

struct IntervalDomain{S, B} <: AbstractDomain
    min::S
    max::S
    periodic::B

    function IntervalDomain(a, b, periodic)
        a, b = promote(a, b)

        return new{typeof(a), typeof(periodic)}(a, b, periodic)
    end
end

function Interval(a, b)
    return IntervalDomain(a, b, false)
end

function Periodic(a, b)
    return IntervalDomain(a, b, true)
end

function Base.show(io::IO, Ω::IntervalDomain)
    min = Ω.min
    max = Ω.max
    periodic = Ω.periodic

    printstyled(io, "[", color = 226)
    printstyled(io, "$min, $max", color = 7)
    if periodic
        printstyled(io, ")", color = 226)
    else
        printstyled(io, "]", color = 226)
    end

    return nothing
end

struct ProductDomain{T} <: AbstractDomain
    domains::T
end

function ×(a::IntervalDomain, b::IntervalDomain)
    return ProductDomain((a, b))
end

function ×(a::ProductDomain, b::IntervalDomain)
    return ProductDomain((a.domains..., b))
end

function ×(a::IntervalDomain, b::ProductDomain)
    return ProductDomain((a, b.domains...))
end

function *(a::AbstractDomain, b::AbstractDomain)
    return a × b
end

function Base.show(io::IO, Ω::ProductDomain)
    for (i, domain) in enumerate(Ω.domains)
        print(domain)

        if i != length(Ω.domains)
            printstyled(io, " × ", color = 118)
        end
    end

    return nothing
end

getindex(Ω::ProductDomain, i::Int) = Ω.domains[i]

ndims(Ω::IntervalDomain) = 1
ndims(Ω::ProductDomain) = +(ndims.(Ω.domains)...)

function periodicityof(Ω::ProductDomain)
    periodicity = ones(Bool, ndims(Ω))

    for i in 1:ndims(Ω)
        periodicity[i] = Ω[i].periodic
    end

    return Tuple(periodicity)
end

```
CubedSphere for GCM
```

struct DeepSphericalShellDomain{S} <: AbstractDomain
    radius::S
    height::S

    function DeepSphericalShellDomain(; radius = nothing, height = nothing)
        radius, height = promote(radius, height)

        return new{typeof(radius)}(radius, height)
    end
end

struct OceanDomain{S} <: AbstractDomain
    radius::S
    depth::S

    function OceanDomain(; radius = nothing, depth = nothing)
        radius, height = promote(radius, depth)

        return new{typeof(radius)}(radius, depth)
    end
end

```
Sphere helper functions for GCM
```
rad(x, y, z) = sqrt(x^2 + y^2 + z^2)
lat(x, y, z) = asin(z / rad(x, y, z)) # ϕ ∈ [-π/2, π/2]
lon(x, y, z) = atan(y, x) # λ ∈ [-π, π)

r̂ⁿᵒʳᵐ(x, y, z) = norm([x, y, z]) ≈ 0 ? 1 : norm([x, y, z])^(-1)
ϕ̂ⁿᵒʳᵐ(x, y, z) = norm([x, y, 0]) ≈ 0 ? 1 : (norm([x, y, z]) * norm([x, y, 0]))^(-1)
λ̂ⁿᵒʳᵐ(x, y, z) = norm([x, y, 0]) ≈ 0 ? 1 : norm([x, y, 0])^(-1)

r̂(x, y, z) = r̂ⁿᵒʳᵐ(x, y, z) * @SVector([x, y, z])
ϕ̂(x, y, z) = ϕ̂ⁿᵒʳᵐ(x, y, z) * @SVector [x * z, y * z, -(x^2 + y^2)]
λ̂(x, y, z) = λ̂ⁿᵒʳᵐ(x, y, z) * @SVector [-y, x, 0]

"""
function DiscontinuousSpectralElementGrid(Ω::ProductDomain; elements = nothing, polynomialorder = nothing)
# Description
Computes a DiscontinuousSpectralElementGrid as specified by a product domain
# Arguments
-`Ω`: A product domain object
# Keyword Arguments TODO: Add brickrange and topology as keyword arguments
-`elements`: A tuple of integers ordered by (Nx, Ny, Nz) for number of elements
-`polynomialorder`: A tupe of integers ordered by (npx, npy, npz) for polynomial order
-`FT`: floattype, assumed Float64 unless otherwise specified
-`mpicomm`: default = MPI.COMM_WORLD
-`array`: default = Array, but should generally be ArrayType
# Return
A DiscontinuousSpectralElementGrid object
"""
function DiscontinuousSpectralElementGrid(
    Ω::ProductDomain{FT},
    elements,
    polynomialorder,
    mpicomm = MPI.COMM_WORLD,
    boundary = (1, 2),
    array = Array,
) where {FT}
    if elements == nothing
        error_message = "Please specify the number of elements as a tuple whose size is commensurate with the domain,"
        error_message = "e.g., a 3 dimensional domain would need a specification like elements = (10,10,10)."
        @error(error_message)
        return nothing
    end

    if polynomialorder == nothing
        error_message = "Please specify the polynomial order as a tuple whose size is commensurate with the domain,"
        error_message = "e.g., a 3 dimensional domain would need a specification like polynomialorder = (3,3,3)."
        @error(error_message)
        return nothing
    end

    dimension = ndims(Ω)

    if (dimension < 2) || (dimension > 3)
        error_message = "SpectralElementGrid only works with dimensions 2 or 3. "
        error_message *= "The current dimension is " * string(ndims(Ω))
        println("The domain is ", Ω)
        @error(error_message)
        return nothing
    end

    if ndims(Ω) != length(elements)
        @error("Specified too many elements for the dimension of the domain")
        return nothing
    end

    if ndims(Ω) != length(polynomialorder)
        @error("Specified too many polynomialorders for the dimension of the domain")
        return nothing
    end

    periodicity = periodicityof(Ω)
    tuple_ranges = []

    for i in 1:dimension
        push!(tuple_ranges, range(FT(Ω[i].min); length = elements[i] + 1, stop = FT(Ω[i].max)))
    end

    brickrange = Tuple(tuple_ranges)
    if boundary == nothing
        boundary = (ntuple(j -> (1, 2), dimension - 1)..., (3, 4))
    end

    topology = StackedBrickTopology(mpicomm, brickrange; periodicity = periodicity, boundary = boundary)

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = array,
        polynomialorder = polynomialorder,
    )
    return grid
end

```
for CubedSphere
```
# using CLIMAParameters
# using CLIMAParameters.Planet: MSLP, R_d, day, grav, Omega, planet_radius
# struct EarthParameterSet <: AbstractEarthParameterSet end
# const param_set = EarthParameterSet()
# _a::Float64 = planet_radius(param_set)
# atmos_height::Float64 = 30e3

# ```
# grid = DiscontinuousSpectralElementGrid(Ω=, elements=, polynomialorder=)
# ```

function DiscontinuousSpectralElementGrid(
    Ω::DeepSphericalShellDomain{FT},
    elements,
    polynomialorder,
    mpicomm = MPI.COMM_WORLD,
    boundary = (1, 2),
    array = Array,
) where {FT}
    Rrange = grid1d(Ω.radius, Ω.radius + Ω.height, nelem = elements.vertical)

    topl = StackedCubedSphereTopology(mpicomm, elements.horizontal, Rrange, boundary = boundary)

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = array,
        polynomialorder = (polynomialorder.horizontal, polynomialorder.vertical),
        meshwarp = ClimateMachine.Mesh.Topologies.equiangular_cubed_sphere_warp,
    )
    return grid
end
function DiscontinuousSpectralElementGrid(
    Ω::OceanDomain{FT},
    elements,
    polynomialorder,
    mpicomm = MPI.COMM_WORLD,
    boundary = (1, 2),
    array = Array,
) where {FT}
    Rrange = grid1d(Ω.radius + Ω.depth, Ω.radius, nelem = elements.vertical)

    topl = StackedCubedSphereTopology(mpicomm, elements.horizontal, Rrange, boundary = boundary)

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = array,
        polynomialorder = (polynomialorder.horizontal, polynomialorder.vertical),
        meshwarp = ClimateMachine.Mesh.Topologies.equiangular_cubed_sphere_warp,
    )
    return grid
end

function is_surface(xc, yc, zc)
    # Sphere case - could dispatch on domain type, maybe?
    height_from_surface = (xc^2 + yc^2 + zc^2)^0.5 - planet_radius(param_set)
    return height_from_surface ≈ 0
end
# end of module

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
module CastepReader

    CastepFortranFile(io) = FortranFile(io, convert="big-endian")

    const IntF = Int32

    """
    Module for electronic Enum types
    """
    module Orbitals
        @enum Spin spinup=1 spindown=2
        @enum Orbital begin
            s=0
            px=10
            py=11
            pz=12
            dxy=20
            dyz=21
            dz2=22
            dxz=23
            dx2=24
            f_xxx=30
            f_yyy=31
            f_zzz=32
            f_xyz=33
            f_z_xx_yy=34
            f_y_zz_xx=35
            f_x_yy_zz=36
        end
        function l(x::Orbital)
            div(Integer(x), 10)
        end
        function m(x::Orbital)
            Integer(x) % 10
        end

        const _castep_orb_names = Dict{Orbital, String}(
            s => "S",
            px => "Px",
            py => "Py",
            pz => "Pz",
            dxy => "Dxy",
            dyz => "Dzy",
            dz2 => "Dzz",
            dxz => "Dzx",
            dx2 => "Dxx-yy",
            f_xxx => "Fxxx",
            f_yyy => "Fyyy",
            f_zzz => "Fzzz",
            f_xyz => "Fxyz",
            f_z_xx_yy => "Fz(xx-yy)",
            f_y_zz_xx => "Fy(zz-xx)",
            f_x_yy_zz => "Fx(yy-zz)",
        )
        function castep_name(orb::Orbital)
            _castep_orb_names[orb]
        end
        const castep_orbital_order = [
            [Orbitals.s],
            [Orbitals.px, Orbitals.py, Orbitals.pz],
            [Orbitals.dz2, Orbitals.dyz, Orbitals.dxz, Orbitals.dx2, Orbitals.dxy],
            [Orbitals.f_xxx, Orbitals.f_yyy, Orbitals.f_zzz, Orbitals.f_xyz, Orbitals.f_z_xx_yy, Orbitals.f_y_zz_xx, Orbitals.f_x_yy_zz]
        ]
        export Spin, Orbital, l, m, castep_name
    end

    using .Orbitals
    export Orbitals

    ## Constants
    using Unitful
    using UnitfulAtomic

    """Factor to convert Hatree to eV"""
    const Ha2eV = ustrip(uconvert(u"eV", 1. * u"Eh_au"))
    """Factor to convert Bohr to Å"""
    const Bohr2Ang = ustrip(uconvert(u"Å", 1. * u"a0_au"))

    include("castep_bin.jl")
    include("wave.jl")
    include("optical_matrix.jl")
    include("pdos_bin.jl")
    include("bands.jl")
    include("cell.jl")
    include("writers/chgcar.jl")
    export read_castep_check, read_cst_ome, read_ome_bin, read_pdos_bin
end


text:
using MASW, JLD2, Plots
Plots.pyplot()

#===
# example to run masw(). Example and SampleData are based on matlab script written by John Schuh.

Matlab script: Copyright (c) 2018, John Schuh All rights reserved.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

•    John Schuh (2020). MASW Dispersion Curve (https://www.mathworks.com/matlabcentral/fileexchange/65138-masw-dispersion-curve), MATLAB Central File Exchange.
      Retrieved January 14, 2020.
===#

curdir = pwd()
fi = jldopen(curdir*"/../example/SampleData.jld", "r")
U = fi["suma_sg"]
close(fi)

fs=1000; # sampling frequency Hz
min_x = 1; # min offset [m] of reciever spread
max_x = 100; # max offset [m] of reciever spread
d_x = 2; # distance between receivers [m]

x = min_x:d_x:max_x; # spatial sampling vector (receiver position)
t=1/fs:1/fs:length(U[:,1])/fs;  # time vector in seconds

figdir = "./"
figname = "dispersionimage.png"

ct=300:1:1200;
freqlimits=(15.0, 80.0);

(f, ct, Udisp, DispersionCurve)=masw(U, x, t, ct, plotfreqlimits=freqlimits,figdir="./",figname="functiontest.png")


text:
# duplication hex example
# add more julia processes
using Distributed
nprocs() < 3 ? addprocs(4-nprocs()) : nothing

using Caesar
using GraffSDK
using UUIDs
## Inter-operating visualization packages for Caesar/RoME/IncrementalInference exist
# using RoMEPlotting

# Graff.....
# 1. Import the initialization code.
cd(joinpath(dirname(pathof(GraffSDK)), "..", "examples"))
# 1a. Create a Configuration
config = loadGraffConfig("synchronyConfig1.json")
#Create a hexagonal sessions
config.sessionId = "Hexagonal_"*replace(string(uuid4()), "-" => "")
println(getGraffConfig())
# 1b. Check the credentials and the service status
printStatus()

# Make sure that the session and robot exist.
if !isRobotExisting()
    @info "Robot $(config.robotId) doesn't exist, creating.."
    newRobot = RobotRequest(config.robotId, "My New Bot", "Description of my neat robot", "Active");
    robot = addRobot(newRobot);
else
    @info "Robot $(config.robotId) already exists, carry on!"
end
# ----- Done!

## Fetch existing marginals from previous sessions
# TODO: new GraffSDK call to fetch desired landmark marginals
existingSessions = getSessions()
# TODO: Get LANDMARK MAP_Ests
# use these LANDMARK estimates as priors in the local version when the same tags seen again
# getNodes()
# Decide which ones we want...

# Make a new session
session = nothing
if !isSessionExisting()
    println(" -- Session '$(config.sessionId)' doesn't exist for robot '$(config.robotId)', creating it...")
    newSessionRequest = SessionDetailsRequest(config.sessionId, "Hexagonal duplication example.", "Pose2", false) #Don't initialize it.
    session = addSession(newSessionRequest)
else
    session = getSession()
end

## LOCAL MEMORY
# start with an empty factor graph object
fg = initfg()
fg.solverParams.isfixedlag = true
fg.qfl = 10

# Local: Add the first pose :x0
addVariable!(fg, :x0, Pose2)
# Graff: Add x0 to Graff
addVariable(:x0, Pose2, String[])

# Add at a fixed location PriorPose2 to pin :x0 to a starting location (0, 0, 0)
prior = IIF.Prior( MvNormal([0; 0; 0], Matrix(Diagonal([0.1;0.1;0.05].^2)) ))
# Local: Add a prior to x0
addFactor!(fg, [:x0], prior)
# Graff: Add the prior to x0
addFactor([:x0], prior)

# Drive around in a hexagon
for i in 0:5
  # define previous and next pose symbols
  psym = Symbol("x$i")
  nsym = Symbol("x$(i+1)")
  deltaMeasurement = [10.0;0;pi/3]
  pOdo = Float64[0.1 0 0; 0 0.1 0; 0 0 0.1]
  pp = Pose2Pose2(MvNormal(deltaMeasurement, pOdo.^2))

  ## LOCAL MEMORY
  # Local: add next pose variable +odometry to local memory version
  addVariable!(fg, nsym, Pose2)
  addFactor!(fg, [psym;nsym], pp )
  # Graff: add odometry
  addVariable(nsym, Pose2, String[])
  addFactor([psym;nsym], pp)
end

# Graff: Let the cloud instance know that the graph is ready
putReady(true)

# perform inference on local version, and remember first runs are slower owing to Julia's just-in-time compiling
batchSolve!(fg)

# No need to call batchSolve on cloud - it picks up on changes and
# 8. Let's check on the solver updates.
sessionLatest = getSession()
while session.lastSolvedTimestamp != sessionLatest.lastSolvedTimestamp
  println("Comparing latest session solver timestamp $(sessionLatest.lastSolvedTimestamp) with original $(session.lastSolvedTimestamp) - still the same so sleeping for 2 seconds")
  sleep(2)
  sessionLatest = getSession()
end
# get cloud nodes
getNodes()

# For Juno/Jupyter style use
pl = drawPoses(fg)
# For scripting use-cases you can export the image
Gadfly.draw(Gadfly.PDF("/tmp/test1.pdf", 20cm, 10cm),pl)  # or PNG(...)

# Add landmarks with Bearing range measurements
addVariable!(fg, :l1, Point2, labels=["LANDMARK"])
addVariable(:l1, Point2, ["LANDMARK"])
p2br = Pose2Point2BearingRange(Normal(0,0.1), Normal(20.0,1.0))
addFactor!(fg, [:x0; :l1], p2br)
addFactor([:x0, :l1], p2br)

# Initialize :l1 numerical values but do not rerun solver
ensureAllInitialized!(fg)
pl = drawPosesLandms(fg)
Gadfly.draw(Gadfly.PDF("/tmp/test2.pdf", 20cm, 10cm),pl)  # or PNG(...)

# Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.1), Normal(20.0,1.0))
addFactor!(fg, [:x6; :l1], p2br2)
addFactor([:x6, :l1], p2br2)

# Graff: PutReady to tell Graff to solve
putReady(true)
# Local: Solve
batchSolve!(fg)

# If you like you can request a full graph solve
requestSessionSolve()

# redraw
pl = drawPosesLandms(fg)
Gadfly.draw(Gadfly.PDF("/tmp/test3.pdf", 20cm, 10cm),pl)  # or PNG(...)

#


text:
module NiLangCore
using MatchCore
using TupleTools

include("lens.jl")

include("utils.jl")
include("symboltable.jl")
include("stack.jl")
include("Core.jl")
include("vars.jl")

include("instr.jl")
include("dualcode.jl")
include("preprocess.jl")
include("variable_analysis.jl")
include("compiler.jl")
include("checks.jl")

if Base.VERSION >= v"1.4.2"
    include("precompile.jl")
    _precompile_()
end

end # module


text:
function train_vKITTI_Conv(model,predict,loss,w; lr=1e-5,epochs=100,batchsize=24,testInterval=1,
    trainVariations=["30-deg-left","30-deg-right","clone","fog","rain"],
    testVariations=["30-deg-left","30-deg-right","clone","fog","rain"],
    trainMode="odd",
    testMode="even"
    )
    @time begin
        w = map(KnetArray,w)
        m = map(zeros,w)
        v = map(zeros,w)
        beta1 = 0.9; beta2 = 0.999; adamT = 1; eps = 1e-8;
        printlnFlush("lr=$(lr)\tepochs=$(epochs)\tbatchsize=$(batchsize)\tmodel=Conv_$(model)\ttestInterval=$(testInterval)")
        printlnFlush(string("Number of parameters: ",sum(map(length,w))," (",Int(floor(sum(map(length,w))/1e6))," million)"))
        printlnFlush("Loading data...")
        @time xtrain,ytrain,colors = loadData_vKITTI(worlds,trainVariations,model,trainMode)
        @time xtest,ytest,colors = loadData_vKITTI(worlds,testVariations,model,testMode)
        printlnFlush("Training...")
        grdnt = gradloss(loss)
        trainbatchCount = Int64(ceil(size(ytrain,4)/batchsize))
        testbatchCount = Int64(ceil(size(ytest,4)/batchsize))
        maxacc = 0
        for epoch=1:epochs
            trainloss = testloss = trainncorrect = testncorrect = trainncount = testncount = 0
            for batch=1:trainbatchCount
                printUpdate(@sprintf("Batch: %5d/%-5d",batch,trainbatchCount))
                #Calculate current indices
                startIndex = 1+(batch-1)*batchsize
                endIndex = min(batch*batchsize,size(ytrain,4))
                #Minibatch
                batchx = KnetArray(xtrain[:,:,:,startIndex:endIndex])
                batchy = ytrain[:,:,:,startIndex:endIndex]
                g,l = grdnt(w,batchx,batchy)
                trainloss += l
                #Adam
                for i in 1:length(w)
                    m[i] = beta1*m[i] + (1-beta1)*g[i]
                    v[i] = beta2*v[i] + (1-beta2)*(g[i].*g[i])
                    mhat = m[i] ./ (1-beta1^adamT)
                    vhat = v[i] ./ (1-beta2^adamT)
                    w[i] = w[i] - (lr ./ (sqrt(vhat)+eps)).*mhat
                end
                adamT += 1
            end
            #Testing accuracy
            if epoch%testInterval == 0 || epoch == 1
                for batch=1:testbatchCount
                    printUpdate(@sprintf("Batch: %5d/%-5d",batch,testbatchCount))
                    #Calculate current indices
                    startIndex = 1+(batch-1)*batchsize
                    endIndex = min(batch*batchsize,size(ytest,4))
                    #Minibatch
                    batchx = KnetArray(xtest[:,:,:,startIndex:endIndex])
                    batchy = ytest[:,:,:,startIndex:endIndex]
                    #Accuracy
                    ypred = Array(predict(w,batchx))
                    testloss += pixelwiseSoftloss(ypred,batchy)
                    testAcc = pixelwiseAccuracy(ypred,batchy)
                    testncorrect += testAcc[1]
                    testncount += testAcc[2]
                end
            end
            testaccuracy = testncorrect/testncount*100
            if isnan(testaccuracy)
                printlnFlush(@sprintf("Epoch: %-6d\tTrain Loss: %0.10f",epoch,trainloss/trainbatchCount))
            else
                printlnFlush(@sprintf("Epoch: %-6d\tTrain Loss: %0.10f\tTest Loss: %0.10f\tTest Accuracy: %0.4f%%",
                    epoch,trainloss/trainbatchCount,testloss/testbatchCount,testaccuracy))
            end
            if testaccuracy > maxacc
                maxacc = testaccuracy
                saveNetwork(w,"vKITTI_Conv_$(model)")
            end
            if epoch == epochs
                printlnFlush(@sprintf("Maximum Test Accuracy: %0.4f%%",maxacc))
            end
        end
    end
end

function train_vKITTI_RNN(model,predict,loss,w; lr=1e-5,epochs=100,testInterval=1,
    trainVariations=["30-deg-left","30-deg-right","clone","fog","rain"],
    testVariations=["30-deg-left","30-deg-right","clone","fog","rain"],
    trainMode="odd",
    testMode="even"
    )
    @time begin
        D = 1
        w = map(KnetArray,w)
        m = map(zeros,w)
        v = map(zeros,w)
        h = KnetArray(zeros(Float32,28,28,512,1))
        beta1 = 0.9; beta2 = 0.999; adamT = 1; eps = 1e-8;
        printlnFlush("lr=$(lr)\tepochs=$(epochs)\tmodel=RNN_$(model)\ttestInterval=$(testInterval)")
        printlnFlush(string("Number of parameters: ",sum(map(length,w))," (",Int(floor(sum(map(length,w))/1e6))," million)"))
        printlnFlush("Loading data...")
        @time xtrain,ytrain,colors = loadData_vKITTI(worlds,trainVariations,model,trainMode)
        @time xtest,ytest,colors = loadData_vKITTI(worlds,testVariations,model,testMode)
        printlnFlush("Training...")
        grdnt = gradloss(loss)
        trainbatchCount = size(ytrain,4)
        testbatchCount = size(ytest,4)
        maxacc = 0
        for epoch=1:epochs
            trainloss = testloss = trainncorrect = testncorrect = trainncount = testncount = 0
            #Train
            h = zeros(h)
            for t=1:trainbatchCount
                printUpdate(@sprintf("Time: %5d/%-5d",t,trainbatchCount))
                g,l = grdnt(w,h,xtrain,ytrain,t,D)
                trainloss += l
                #Adam
                for i in 1:length(w)
                    m[i] = beta1*m[i] + (1-beta1)*g[i]
                    v[i] = beta2*v[i] + (1-beta2)*(g[i].*g[i])
                    mhat = m[i] ./ (1-beta1^adamT)
                    vhat = v[i] ./ (1-beta2^adamT)
                    w[i] = w[i] - (lr ./ (sqrt(vhat)+eps)).*mhat
                end
                adamT += 1
            end
            #Testing accuracy
            if epoch%testInterval == 0 || epoch == 1
                h = zeros(h)
                for t=1:testbatchCount
                    printUpdate(@sprintf("Time: %5d/%-5d",t,testbatchCount))
                    ypred = Array(predict(w,h,KnetArray(xtest[:,:,:,t:t])))
                    testloss += pixelwiseSoftloss(ypred,ytest[:,:,:,t:t])
                    testAcc = pixelwiseAccuracy(ypred,ytest[:,:,:,t:t])
                    testncorrect += testAcc[1]
                    testncount += testAcc[2]
                end
            end
            testaccuracy = testncorrect/testncount*100
            if isnan(testaccuracy)
                printlnFlush(@sprintf("Epoch: %-6d\tTrain Loss: %0.10f",epoch,trainloss/trainbatchCount))
            else
                printlnFlush(@sprintf("Epoch: %-6d\tTrain Loss: %0.10f\tTest Loss: %0.10f\tTest Accuracy: %0.4f%%",
                    epoch,trainloss/trainbatchCount,testloss/testbatchCount,testaccuracy))
            end
            if testaccuracy > maxacc
                maxacc = testaccuracy
                saveNetwork(w,"vKITTI_RNN_$(model)")
            end
            if epoch == epochs
                printlnFlush(@sprintf("Maximum Test Accuracy: %0.4f%%",maxacc))
            end
        end
    end
end

function train_vKITTI_LSTM(model,predict,loss,w; lr=1e-5,epochs=100,testInterval=1,
    trainVariations=["30-deg-left","30-deg-right","clone","fog","rain"],
    testVariations=["30-deg-left","30-deg-right","clone","fog","rain"],
    trainMode="odd",
    testMode="even"
    )
    @time begin
        D = 1
        w = map(KnetArray,w)
        m = map(zeros,w)
        v = map(zeros,w)
        h = KnetArray(zeros(Float32,28,28,512,1))
        c = KnetArray(zeros(Float32,28,28,512,1))
        beta1 = 0.9; beta2 = 0.999; adamT = 1; eps = 1e-8;
        printlnFlush("lr=$(lr)\tepochs=$(epochs)\tmodel=LSTM_$(model)\ttestInterval=$(testInterval)")
        printlnFlush(string("Number of parameters: ",sum(map(length,w))," (",Int(floor(sum(map(length,w))/1e6))," million)"))
        printlnFlush("Loading data...")
        @time xtrain,ytrain,colors = loadData_vKITTI(worlds,trainVariations,model,trainMode)
        @time xtest,ytest,colors = loadData_vKITTI(worlds,testVariations,model,testMode)
        printlnFlush("Training...")
        grdnt = gradloss(loss)
        trainbatchCount = size(ytrain,4)
        testbatchCount = size(ytest,4)
        maxacc = 0
        for epoch=1:epochs
            trainloss = testloss = trainncorrect = testncorrect = trainncount = testncount = 0
            #Train
            h = zeros(h)
            c = zeros(c)
            for t=1:trainbatchCount
                printUpdate(@sprintf("Time: %5d/%-5d",t,trainbatchCount))
                g,l = grdnt(w,h,c,xtrain,ytrain,t,D)
                trainloss += l
                #Adam
                for i in 1:length(w)
                    m[i] = beta1*m[i] + (1-beta1)*g[i]
                    v[i] = beta2*v[i] + (1-beta2)*(g[i].*g[i])
                    mhat = m[i] ./ (1-beta1^adamT)
                    vhat = v[i] ./ (1-beta2^adamT)
                    w[i] = w[i] - (lr ./ (sqrt(vhat)+eps)).*mhat
                end
                adamT += 1
            end
            #Testing accuracy
            if epoch%testInterval == 0 || epoch == 1
                h = zeros(h)
                c = zeros(c)
                for t=1:testbatchCount
                    printUpdate(@sprintf("Time: %5d/%-5d",t,testbatchCount))
                    ypred = Array(predict(w,h,c,KnetArray(xtest[:,:,:,t:t])))
                    testloss += pixelwiseSoftloss(ypred,ytest[:,:,:,t:t])
                    testAcc = pixelwiseAccuracy(ypred,ytest[:,:,:,t:t])
                    testncorrect += testAcc[1]
                    testncount += testAcc[2]
                end
            end
            testaccuracy = testncorrect/testncount*100
            if isnan(testaccuracy)
                printlnFlush(@sprintf("Epoch: %-6d\tTrain Loss: %0.10f",epoch,trainloss/trainbatchCount))
            else
                printlnFlush(@sprintf("Epoch: %-6d\tTrain Loss: %0.10f\tTest Loss: %0.10f\tTest Accuracy: %0.4f%%",
                    epoch,trainloss/trainbatchCount,testloss/testbatchCount,testaccuracy))
            end
            if testaccuracy > maxacc
                maxacc = testaccuracy
                saveNetwork(w,"vKITTI_LSTM_$(model)")
            end
            if epoch == epochs
                printlnFlush(@sprintf("Maximum Test Accuracy: %0.4f%%",maxacc))
            end
        end
    end
end

function train_vKITTI_GRU(model,predict,loss,w; lr=1e-5,epochs=100,testInterval=1,
    trainVariations=["30-deg-left","30-deg-right","clone","fog","rain"],
    testVariations=["30-deg-left","30-deg-right","clone","fog","rain"],
    trainMode="odd",
    testMode="even"
    )
    @time begin
        D = 1
        w = map(KnetArray,w)
        m = map(zeros,w)
        v = map(zeros,w)
        h = KnetArray(zeros(Float32,28,28,512,1))
        beta1 = 0.9; beta2 = 0.999; adamT = 1; eps = 1e-8;
        printlnFlush("lr=$(lr)\tepochs=$(epochs)\tmodel=GRU_$(model)\ttestInterval=$(testInterval)")
        printlnFlush(string("Number of parameters: ",sum(map(length,w))," (",Int(floor(sum(map(length,w))/1e6))," million)"))
        printlnFlush("Loading data...")
        @time xtrain,ytrain,colors = loadData_vKITTI(worlds,trainVariations,model,trainMode)
        @time xtest,ytest,colors = loadData_vKITTI(worlds,testVariations,model,testMode)
        printlnFlush("Training...")
        grdnt = gradloss(loss)
        trainbatchCount = size(ytrain,4)
        testbatchCount = size(ytest,4)
        maxacc = 0
        for epoch=1:epochs
            trainloss = testloss = trainncorrect = testncorrect = trainncount = testncount = 0
            #Train
            h = zeros(h)
            for t=1:trainbatchCount
                printUpdate(@sprintf("Time: %5d/%-5d",t,trainbatchCount))
                g,l = grdnt(w,h,xtrain,ytrain,t,D)
                trainloss += l
                #Adam
                for i in 1:length(w)
                    m[i] = beta1*m[i] + (1-beta1)*g[i]
                    v[i] = beta2*v[i] + (1-beta2)*(g[i].*g[i])
                    mhat = m[i] ./ (1-beta1^adamT)
                    vhat = v[i] ./ (1-beta2^adamT)
                    w[i] = w[i] - (lr ./ (sqrt(vhat)+eps)).*mhat
                end
                adamT += 1
            end
            #Testing accuracy
            if epoch%testInterval == 0 || epoch == 1
                h = zeros(h)
                for t=1:testbatchCount
                    printUpdate(@sprintf("Time: %5d/%-5d",t,testbatchCount))
                    ypred = Array(predict(w,h,KnetArray(xtest[:,:,:,t:t])))
                    testloss += pixelwiseSoftloss(ypred,ytest[:,:,:,t:t])
                    testAcc = pixelwiseAccuracy(ypred,ytest[:,:,:,t:t])
                    testncorrect += testAcc[1]
                    testncount += testAcc[2]
                end
            end
            testaccuracy = testncorrect/testncount*100
            if isnan(testaccuracy)
                printlnFlush(@sprintf("Epoch: %-6d\tTrain Loss: %0.10f",epoch,trainloss/trainbatchCount))
            else
                printlnFlush(@sprintf("Epoch: %-6d\tTrain Loss: %0.10f\tTest Loss: %0.10f\tTest Accuracy: %0.4f%%",
                    epoch,trainloss/trainbatchCount,testloss/testbatchCount,testaccuracy))
            end
            if testaccuracy > maxacc
                maxacc = testaccuracy
                saveNetwork(w,"vKITTI_GRU_$(model)")
            end
            if epoch == epochs
                printlnFlush(@sprintf("Maximum Test Accuracy: %0.4f%%",maxacc))
            end
        end
    end
end

#E1 : Degrees
#E2 : Half Train/Test

#   Conv
#E1
function run_vKITTI_Conv_D_E1_train()
    train_vKITTI_Conv("D_VGG",vKITTI_Conv_D,vKITTI_Conv_D_loss,vKITTI_Conv_D_w(); lr=1e-5,epochs=100,batchsize=24,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","clone"],
        testVariations=["15-deg-left","15-deg-right"],
        trainMode="full",
        testMode="full"
    )
end
function run_vKITTI_Conv_RGB_E1_train()
    train_vKITTI_Conv("RGB_VGG",vKITTI_Conv_RGB,vKITTI_Conv_RGB_loss,vKITTI_Conv_RGB_w(); lr=1e-5,epochs=100,batchsize=24,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","clone"],
        testVariations=["15-deg-left","15-deg-right"],
        trainMode="full",
        testMode="full"
    )
end
function run_vKITTI_Conv_RGBD_E1_train()
    train_vKITTI_Conv("RGBD_VGG",vKITTI_Conv_RGBD,vKITTI_Conv_RGBD_loss,vKITTI_Conv_RGBD_w(); lr=1e-5,epochs=50,batchsize=16,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","clone"],
        testVariations=["15-deg-left","15-deg-right"],
        trainMode="full",
        testMode="full"
    )
end
#E2
function run_vKITTI_Conv_D_E2_train()
    train_vKITTI_Conv("D_VGG",vKITTI_Conv_D,vKITTI_Conv_D_loss,vKITTI_Conv_D_w(); lr=1e-5,epochs=100,batchsize=24,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","15-deg-left","15-deg-right","clone","fog","morning","overcast","rain","sunset"],
        testVariations=["30-deg-left","30-deg-right","15-deg-left","15-deg-right","clone","fog","morning","overcast","rain","sunset"],
        trainMode="firstHalf",
        testMode="secondHalf"
    )
end
function run_vKITTI_Conv_RGB_E2_train()
    train_vKITTI_Conv("RGB_VGG",vKITTI_Conv_RGB,vKITTI_Conv_RGB_loss,vKITTI_Conv_RGB_w(); lr=1e-5,epochs=100,batchsize=24,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","15-deg-left","15-deg-right","clone","fog","morning","overcast","rain","sunset"],
        testVariations=["30-deg-left","30-deg-right","15-deg-left","15-deg-right","clone","fog","morning","overcast","rain","sunset"],
        trainMode="firstHalf",
        testMode="secondHalf"
    )
end
function run_vKITTI_Conv_RGBD_E2_train()
    train_vKITTI_Conv("RGBD_VGG",vKITTI_Conv_RGBD,vKITTI_Conv_RGBD_loss,vKITTI_Conv_RGBD_w(); lr=1e-5,epochs=50,batchsize=16,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","15-deg-left","15-deg-right","clone","fog","morning","overcast","rain","sunset"],
        testVariations=["30-deg-left","30-deg-right","15-deg-left","15-deg-right","clone","fog","morning","overcast","rain","sunset"],
        trainMode="firstHalf",
        testMode="secondHalf"
    )
end

#   RNN
function run_vKITTI_RNN_D_E1_train()
    train_vKITTI_RNN("D_VGG",vKITTI_RNN_D,vKITTI_RNN_D_loss,vKITTI_RNN_D_w(); lr=1e-5,epochs=25,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","clone"],
        testVariations=["15-deg-left","15-deg-right"],
        trainMode="full",
        testMode="full"
    )
end
function run_vKITTI_RNN_RGB_E1_train()
    train_vKITTI_RNN("RGB_VGG",vKITTI_RNN_RGB,vKITTI_RNN_RGB_loss,vKITTI_RNN_RGB_w(); lr=1e-5,epochs=25,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","clone"],
        testVariations=["15-deg-left","15-deg-right"],
        trainMode="full",
        testMode="full"
    )
end
function run_vKITTI_RNN_RGBD_E1_train()
    train_vKITTI_RNN("RGBD_VGG",vKITTI_RNN_RGBD,vKITTI_RNN_RGBD_loss,vKITTI_RNN_RGBD_w(); lr=1e-5,epochs=25,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","clone"],
        testVariations=["15-deg-left","15-deg-right"],
        trainMode="full",
        testMode="full"
    )
end

#   LSTM
#E1
function run_vKITTI_LSTM_D_E1_train()
    train_vKITTI_LSTM("D_VGG",vKITTI_LSTM_D,vKITTI_LSTM_D_loss,vKITTI_LSTM_D_w(); lr=1e-5,epochs=25,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","clone"],
        testVariations=["15-deg-left","15-deg-right"],
        trainMode="full",
        testMode="full"
    )
end
function run_vKITTI_LSTM_RGB_E1_train()
    train_vKITTI_LSTM("RGB_VGG",vKITTI_LSTM_RGB,vKITTI_LSTM_RGB_loss,vKITTI_LSTM_RGB_w(); lr=1e-5,epochs=25,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","clone"],
        testVariations=["15-deg-left","15-deg-right"],
        trainMode="full",
        testMode="full"
    )
end
function run_vKITTI_LSTM_RGBD_E1_train()
    train_vKITTI_LSTM("RGBD_VGG",vKITTI_LSTM_RGBD,vKITTI_LSTM_RGBD_loss,vKITTI_LSTM_RGBD_w(); lr=1e-5,epochs=25,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","clone"],
        testVariations=["15-deg-left","15-deg-right"],
        trainMode="full",
        testMode="full"
    )
end
#E2
function run_vKITTI_LSTM_D_E2_train()
    train_vKITTI_LSTM("D_VGG",vKITTI_LSTM_D,vKITTI_LSTM_D_loss,vKITTI_LSTM_D_w(); lr=1e-5,epochs=25,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","15-deg-left","15-deg-right","clone","fog","morning","overcast","rain","sunset"],
        testVariations=["30-deg-left","30-deg-right","15-deg-left","15-deg-right","clone","fog","morning","overcast","rain","sunset"],
        trainMode="firstHalf",
        testMode="secondHalf"
    )
end
function run_vKITTI_LSTM_RGB_E2_train()
    train_vKITTI_LSTM("RGB_VGG",vKITTI_LSTM_RGB,vKITTI_LSTM_RGB_loss,vKITTI_LSTM_RGB_w(); lr=1e-5,epochs=25,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","15-deg-left","15-deg-right","clone","fog","morning","overcast","rain","sunset"],
        testVariations=["30-deg-left","30-deg-right","15-deg-left","15-deg-right","clone","fog","morning","overcast","rain","sunset"],
        trainMode="firstHalf",
        testMode="secondHalf"
    )
end
function run_vKITTI_LSTM_RGBD_E2_train()
    train_vKITTI_LSTM("RGBD_VGG",vKITTI_LSTM_RGBD,vKITTI_LSTM_RGBD_loss,vKITTI_LSTM_RGBD_w(); lr=1e-5,epochs=25,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","15-deg-left","15-deg-right","clone","fog","morning","overcast","rain","sunset"],
        testVariations=["30-deg-left","30-deg-right","15-deg-left","15-deg-right","clone","fog","morning","overcast","rain","sunset"],
        trainMode="firstHalf",
        testMode="secondHalf"
    )
end

#   GRU
function run_vKITTI_GRU_D_E1_train()
    train_vKITTI_GRU("D_VGG",vKITTI_GRU_D,vKITTI_GRU_D_loss,vKITTI_GRU_D_w(); lr=1e-5,epochs=25,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","clone"],
        testVariations=["15-deg-left","15-deg-right"],
        trainMode="full",
        testMode="full"
    )
end
function run_vKITTI_GRU_RGB_E1_train()
    train_vKITTI_GRU("RGB_VGG",vKITTI_GRU_RGB,vKITTI_GRU_RGB_loss,vKITTI_GRU_RGB_w(); lr=1e-5,epochs=25,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","clone"],
        testVariations=["15-deg-left","15-deg-right"],
        trainMode="full",
        testMode="full"
    )
end
function run_vKITTI_GRU_RGBD_E1_train()
    train_vKITTI_GRU("RGBD_VGG",vKITTI_GRU_RGBD,vKITTI_GRU_RGBD_loss,vKITTI_GRU_RGBD_w(); lr=1e-5,epochs=25,testInterval=1,
        trainVariations=["30-deg-left","30-deg-right","clone"],
        testVariations=["15-deg-left","15-deg-right"],
        trainMode="full",
        testMode="full"
    )
end

#   Conv
#E1
export run_vKITTI_Conv_D_E1_train
export run_vKITTI_Conv_RGB_E1_train
export run_vKITTI_Conv_RGBD_E1_train
#E2
export run_vKITTI_Conv_D_E2_train
export run_vKITTI_Conv_RGB_E2_train
export run_vKITTI_Conv_RGBD_E2_train
#   RNN
export run_vKITTI_RNN_D_E1_train
export run_vKITTI_RNN_RGB_E1_train
export run_vKITTI_RNN_RGBD_E1_train
#   LSTM
#E1
export run_vKITTI_LSTM_D_E1_train
export run_vKITTI_LSTM_RGB_E1_train
export run_vKITTI_LSTM_RGBD_E1_train
#E2
export run_vKITTI_LSTM_D_E2_train
export run_vKITTI_LSTM_RGB_E2_train
export run_vKITTI_LSTM_RGBD_E2_train
#   GRU
export run_vKITTI_GRU_D_E1_train
export run_vKITTI_GRU_RGB_E1_train
export run_vKITTI_GRU_RGBD_E1_train

max_stars_repo_name: eeyrdkl/MScThesis

max_issues_repo_name: eeyrdkl/MScThesis

max_forks_repo_name: eeyrdkl/MScThesis


text:
module DotMaps

export DotMap

"""
   DotMaps.DotMap(::AbstractDict)

Constructs a DotMap from a Dict.  This provides the same functionaliity as dictionaries, but allows indexing with `.` instead of (or in addition to) `[""]`.
"""
struct DotMap
    __dict__ ::Dict{Symbol,Any}
    DotMap() = new(Dict{Symbol,Any}())
end

function DotMap(d::AbstractDict)
   dm = DotMap()
   for (k, v) in d
      dm.__dict__[Symbol(k)] = DotMap(v)
   end
   return dm
end

DotMap(d::Any) = d

"""
   DotMaps.todict(::DotMap; keys_as_strings=false)

Constructs a Dict from a DotMap. If `keys_as_strings`, the keys will be `String` instead of `Symbol`.
"""
function todict(obj::DotMap; keys_as_strings::Bool = false)
   dict = Dict()
   for (k, v) in obj
      nk = keys_as_strings ? string(k) : k
      dict[nk] = todict(v, keys_as_strings = keys_as_strings)
   end

   return dict
end

# return at leaves
todict(obj::Any; keys_as_strings::Bool = false) = obj

# make dots work
function Base.getproperty(obj::DotMap, name::Symbol)
   if name == :__dict__
      return getfield(obj, name)
   else
      return obj.__dict__[name]
   end
end

# make dictionary indexing work
Base.getindex(obj::DotMap, name::Symbol) = Base.getindex(obj.__dict__, name)
Base.getindex(obj::DotMap, name::Any) = Base.getindex(obj, Symbol(name))

Base.setindex!(obj::DotMap, x, name::Symbol) = Base.setindex!(obj.__dict__, x, name)
Base.setindex!(obj::DotMap, x::Dict, name::Symbol) = Base.setindex!(obj.__dict__, DotMap(x), name)
Base.setindex!(obj::DotMap, x, name) = Base.setindex!(obj, x, Symbol(name))

# assignment with dots
Base.setproperty!(obj::DotMap, name::Symbol, x) = Base.setindex!(obj, x, name)

# iteration
Base.iterate(obj::DotMap) = Base.iterate(obj.__dict__)
Base.iterate(obj::DotMap, s::Any) = Base.iterate(obj.__dict__, s)
Base.length(obj::DotMap) = Base.length(obj.__dict__)
Base.firstindex(obj::DotMap) = Base.firstindex(obj.__dict__)
Base.lastindex(obj::DotMap) = Base.lastindex(obj.__dict__)

# dictionary methods
Base.keys(obj::DotMap) = Base.keys(obj.__dict__)
Base.values(obj::DotMap) = Base.values(obj.__dict__)
Base.collect(obj::DotMap) = Base.collect(obj.__dict__)

# retrieval/modification
Base.get(obj::DotMap, key::Symbol, default) = Base.get(obj.__dict__, key, default)
Base.get(obj::DotMap, key::Any, default) = Base.get(obj.__dict__, Symbol(key), default)
Base.get!(obj::DotMap, key::Symbol, default) = Base.get!(obj.__dict__, key, default)
Base.get!(obj::DotMap, key::Any, default) = Base.get!(obj.__dict__, Symbol(key), default)
Base.getkey(obj::DotMap, key::Symbol, default) = Base.getkey(obj.__dict__, key, default)
Base.getkey(obj::DotMap, key::Any, default) = Base.getkey(obj.__dict__, Symbol(key), default)
Base.pop!(obj::DotMap, key::Symbol) = Base.pop!(obj.__dict__, key)
Base.pop!(obj::DotMap, key::Any) = Base.pop!(obj.__dict__, Symbol(key))
Base.delete!(obj::DotMap, key::Symbol) = Base.pop!(obj.__dict__, key)
Base.delete!(obj::DotMap, key::Any) = Base.pop!(obj.__dict__, Symbol(key))
Base.filter!(pred, obj::DotMap) = Base.filter!(pred, obj.__dict__)
Base.filter(pred, obj::DotMap) = DotMap(Base.filter(pred, obj.__dict__))

# about this dict
Base.propertynames(obj::DotMap) = Base.keys(obj.__dict__)
Base.isempty(obj::DotMap) = Base.isempty(obj.__dict__)
Base.show(obj::DotMap) = Base.show(obj.__dict__)

# containment
Base.in(key::Symbol, obj::DotMap) = Base.in(key, obj.__dict__)
Base.in(key::Any, obj::DotMap) = Base.in(Symbol(key), obj.__dict__)
Base.haskey(obj::DotMap, key::Symbol) = Base.haskey(obj.__dict__, key)
Base.haskey(obj::DotMap, key::Any) = Base.haskey(obj.__dict__, Symbol(key))

end


text:
baremodule Maybe

export @something
# export @?

macro something end

function length end
function eltype end
function first end
function last end
function get end
function getindex end
function getproperty end

function _break end

"""
    Maybe.Extras

A namespace for extra API; this is for preserving `Maybe.*` namespace
for (mainly) `Base`-compatible API.

Since there is no name clash with `Base` API, `using Maybe.Extras`
imports the API defined in `Maybe.Extras`.
"""
baremodule Extras
export getnested, ifnothing, maybe, defaultto
function maybe end
function ifnothing end
function defaultto end
function getnested end
end

const X = Extras

module Implementations
using ..Maybe: Maybe
using Base: @propagate_inbounds
using Base.Meta: isexpr
using ExprTools: combinedef, splitdef

import ..Maybe: @something

if !@isdefined hasproperty
    using Compat: hasproperty
end

if !@isdefined isnothing
    using Compat: isnothing
end

if !isdefined(Base, Symbol("@var_str"))
    macro var_str(x)
        return Symbol(x)
    end
end

include("utils.jl")
include("lift.jl")
include("something.jl")
include("core.jl")
include("extras.jl")
finalize_implementations()
end

Implementations.finalize_package()

end


text:
using InstantiateFromURL

# activate the QuantEcon environment
activate_github("QuantEcon/QuantEconLecturePackages", tag = "v0.9.0")

using LinearAlgebra, Statistics, Compat # load common packages

max_stars_repo_name: leopoldomig/lecture-source-jl

max_issues_repo_name: leopoldomig/lecture-source-jl

max_forks_repo_name: leopoldomig/lecture-source-jl


text:
using CSV, DataFrames, XLSX, CairoMakie, Jedi, ColorSchemes, Statistics

CairoMakie.activate!()
# Set Style
Jedi.styles.default_makie!()

sequential_cmaps = [
    "Blues",
    "Reds",
    "Greens",
    "Oranges",
    "Purples",
]

##
# Find home directory
dir = @__DIR__
home_dir = joinpath(split(dir, '/')[1:end-4]...)

# Find date
workdir = split(dir, '/')[end]
DATE = parse(Int64, split(workdir, '_')[1])
RUN_NO = parse(Int64, split(workdir, '_')[2][end])
##

# Load the data.
file = "/$(home_dir)/data/plate_reader/$workdir.txt"

xf = XLSX.readxlsx(dir * "/$(DATE)_plate_layout.xlsx")

# see all sheet names
layout = XLSX.sheetnames(xf)

# Read layout information
layout_info = DataFrames.DataFrame()
# Loop through layout info
for table in layout
    flattened = vec(hcat(XLSX.getdata(xf[table])...))
    DataFrames.insertcols!(layout_info, table => flattened)
end

# Define rows and columns read on plate
DataFrames.insertcols!(layout_info, "row" => [x[1] for x in layout_info[!, "well"]])
DataFrames.insertcols!(layout_info, "col" => [x[2:end] for x in layout_info[!, "well"]])

##
# Initialize plot

fig = Figure(resolution=(900, 300 * ((length(layout) + 1) ÷ 2)))

# Loop through features
i=0
for feature in layout
    cmap = rand(sequential_cmaps)

    # If well, don't add information
    if (feature == "well")
        continue
    end
    i += 1

    println(feature)
    ax = Axis(fig[(i-1) ÷ 2 + 1, ((i-1) % 2) + 1 ])
    # Extract strain data and pivot dataframe to obtain proper dimensions
    data = layout_info[!, ["row", "col", feature]]
    # Add code for categorical data
    cat_dict = Dict(unique(data[!, feature]) .=> collect(1:length(unique(data[!, feature]))))
    data[!, "feature_code"] = [cat_dict[x] for x in data[!, feature]]
    row_dict = Dict(unique(data[!, "row"]) .=> collect(1:length(unique(data[!, "row"]))))
    data[!, "row_code"] = [row_dict[x] for x in data[!, "row"]]
    col_dict = Dict(unique(data[!, "col"]) .=> collect(1:length(unique(data[!, "col"]))))
    data[!, "col_code"] = [col_dict[x] for x in data[!, "col"]]

    # Pivot dataframe for both code and label to generate heatmap
    cmap = colorschemes[Symbol("$(sequential_cmaps[i])_$(max(3, length(unique(data.feature_code))))")]

    # Generate colormap
    heatmap!(ax, data.col_code, data.row_code, data.feature_code, strokewidth = 1, strokecolor="black", colormap=cmap)
    text!(ax, data[!, feature], position = Point.(data.col_code, data.row_code), align = (:center, :center),
    offset = (0, 0), color = :black, textsize=3)
    # Set plot title
    ax.title = feature
    ax.xgridcolor="white"
    ax.ygridcolor="white"
    ax.xgridwidth=5
    ax.xticks=(collect(1:12), string.(collect(1:12)))
    ax.yticks=(collect(1:8), ["A", "B", "C", "D", "E", "F", "G", "H"])
    ax.yreversed = true
end

# Save it to the output file
if ~ispath("$dir/output")  # Check if directyr exists
    mkdir("$dir/output")  # Generate directory if required
end
fig
save("$dir/output/plate_layout.pdf", fig)


text:
# ============================================================
# Theories
# ============================================================

export prim_conserve,
       conserve_prim,
       mixture_prim_conserve,
       mixture_conserve_prim,
       em_coefficients,
       advection_flux,
       burgers_flux,
       euler_flux,
       euler_jacobi
export gauss_moments,
       mixture_gauss_moments,
       moments_conserve,
       diatomic_moments_conserve,
       mixture_moments_conserve,
       pdf_slope,
       mixture_pdf_slope,
       moments_conserve_slope,
       mixture_moments_conserve_slope,
       discrete_moments,
       pressure,
       stress,
       heat_flux,
       maxwellian,
       maxwellian!,
       mixture_maxwellian,
       mixture_maxwellian!,
       shakhov,
       shakhov!,
       rykov!,
       reduce_distribution,
       full_distribution,
       ref_vhs_vis,
       vhs_collision_time,
       rykov_zr,
       aap_hs_collision_time,
       aap_hs_prim,
       aap_hs_diffeq!,
       shift_pdf!,
       hs_boltz_kn,
       kernel_mode,
       boltzmann_fft,
       boltzmann_fft!,
       boltzmann_ode!,
       boltzmann_nuode!,
       bgk_ode!,
       chapman_enskog
export heat_capacity_ratio,
       sound_speed
export sample_maxwell,
       next_collision_time,
       boundary_time

include("theory_continuum.jl")
include("theory_plasma.jl")
include("theory_maxwellian.jl")
include("theory_shakhov.jl")
include("theory_diatomic.jl")
include("theory_atom.jl")
include("theory_pdf.jl")
include("theory_fsm.jl")
include("theory_moments_pure.jl")
include("theory_moments_mixture.jl")
include("theory_thermo.jl")
include("theory_particle.jl")