text:
using Kinetic, Plots, LinearAlgebra
using ProgressMeter: @showprogress

cd(@__DIR__)
D = Dict{Symbol,Any}()
begin
    D[:matter] = "gas"
    D[:case] = "cavity"
    D[:space] = "2d2f2v"
    D[:flux] = "kfvs"
    D[:collision] = "bgk"
    D[:nSpecies] = 1
    D[:interpOrder] = 2
    D[:limiter] = "vanleer"
    D[:boundary] = "maxwell"
    D[:cfl] = 0.3
    D[:maxTime] = 5.0

    D[:x0] = 0.0
    D[:x1] = 1.0
    D[:nx] = 41
    D[:y0] = 0.0
    D[:y1] = 1.0
    D[:ny] = 41
    D[:pMeshType] = "uniform"
    D[:nxg] = 1
    D[:nyg] = 1

    D[:umin] = -5.0
    D[:umax] = 5.0
    D[:nu] = 8
    D[:vmin] = -5.0
    D[:vmax] = 5.0
    D[:nv] = 8
    D[:vMeshType] = "rectangle"
    D[:nug] = 0
    D[:nvg] = 0

    D[:knudsen] = 0.005
    D[:mach] = 0.0
    D[:prandtl] = 1.0
    D[:inK] = 1.0
    D[:omega] = 0.81
    D[:alphaRef] = 1.0
    D[:omegaRef] = 0.5

    D[:uLid] = 0.15
    D[:vLid] = 0.0
    D[:tLid] = 1.0
end

ks = SolverSet(D)
ctr, a1face, a2face = init_fvm(ks, ks.ps, :dynamic_array; structarray = true)

t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(ks.set.maxTime ÷ dt) + 1
res = zero(ks.ib.wL)

@showprogress for iter = 1:200#nt
    #reconstruct!(ks, ctr)
    #evolve!(ks, ctr, a1face, a2face, dt; mode = Symbol(ks.set.flux), bc = Symbol(ks.set.boundary))

    # horizontal flux
    @inbounds Threads.@threads for j = 1:ks.pSpace.ny
        for i = 2:ks.pSpace.nx
            flux_gks!(
                a1face[i, j].fw,
                ctr[i-1, j].w .+ ctr[i-1, j].sw[:, 1] .* ks.ps.dx[i-1, j]/2,
                ctr[i, j].w .- ctr[i, j].sw[:, 1] .* ks.ps.dx[i, j]/2,
                ks.gas.K,
                ks.gas.γ,
                ks.gas.μᵣ,
                ks.gas.ω,
                dt,
                ks.ps.dx[i-1, j]/2,
                ks.ps.dx[i, j]/2,
                a1face[i, j].len,
                ctr[i-1, j].sw[:, 1],
                ctr[i, j].sw[:, 1],
            )
        end
    end

    # vertical flux
    @inbounds Threads.@threads for j = 2:ks.pSpace.ny
        for i = 1:ks.pSpace.nx
            wL = KitBase.local_frame(ctr[i, j-1].w, 0., 1.)
            wR = KitBase.local_frame(ctr[i, j].w, 0., 1.)
            swL = KitBase.local_frame(ctr[i, j-1].sw[:, 2], 0., 1.)
            swR = KitBase.local_frame(ctr[i, j].sw[:, 2], 0., 1.)

            flux_gks!(
                a2face[i, j].fw,
                wL .+ swL .* ks.ps.dy[i, j-1]/2,
                wR .- swR .* ks.ps.dy[i, j]/2,
                ks.gas.K,
                ks.gas.γ,
                ks.gas.μᵣ,
                ks.gas.ω,
                dt,
                ks.ps.dy[i, j-1]/2,
                ks.ps.dy[i, j]/2,
                a2face[i, j].len,
                swL,
                swR,
            )

            a2face[i, j].fw .= KitBase.global_frame(a2face[i, j].fw, 0., 1.)
        end
    end

    # boundary flux
    @inbounds Threads.@threads for j = 1:ks.pSpace.ny
        KitBase.flux_boundary_maxwell!(
            a1face[1, j].fw,
            ks.ib.bcL,
            ctr[1, j].w,
            ks.gas.K,
            ks.gas.γ,
            dt,
            ctr[1, j].dy,
            1.,
        )

        KitBase.flux_boundary_maxwell!(
            a1face[ks.pSpace.nx+1, j].fw,
            ks.ib.bcR,
            ctr[ks.pSpace.nx, j].w,
            ks.gas.K,
            ks.gas.γ,
            dt,
            ctr[ks.pSpace.nx, j].dy,
            -1.,
        )
    end

    @inbounds Threads.@threads for i = 1:ks.pSpace.nx
        w = KitBase.local_frame(ctr[i, 1].w, 0., 1.)
        KitBase.flux_boundary_maxwell!(
            a2face[i, 1].fw,
            ks.ib.bcD,
            w,
            ks.gas.K,
            ks.gas.γ,
            dt,
            ctr[i, 1].dx,
            1,
        )
        a2face[i, 1].fw .= KitBase.global_frame(a2face[i, 1].fw, 0., 1.)

        w1 = KitBase.local_frame(ctr[i, ks.pSpace.ny].w, 0., 1.)
        KitBase.flux_boundary_maxwell!(
            a2face[i, ks.pSpace.ny+1].fw,
            [1., 0.0, -0.05, 1.0],
            w1,
            ks.gas.K,
            ks.gas.γ,
            dt,
            ctr[i, ks.pSpace.ny].dy,
            -1,
        )
        a2face[i, ks.pSpace.ny+1].fw .= KitBase.global_frame(
            a2face[i, ks.pSpace.ny+1].fw,
            0.,
            1.,
        )
    end

    #update!(ks, ctr, a1face, a2face, dt, res; coll = Symbol(ks.set.collision), bc = Symbol(ks.set.boundary))

    @inbounds Threads.@threads for j = 1:ks.pSpace.ny
        for i = 1:ks.pSpace.nx
            @. ctr[i, j].w += (a1face[i, j].fw - a1face[i+1, j].fw + a2face[i, j].fw - a2face[i, j+1].fw) / ctr[i, j].dx / ctr[i, j].dy
            ctr[i, j].prim .= conserve_prim(ctr[i, j].w, ks.gas.γ)
        end
    end

    t += dt
end

plot_contour(ks, ctr)

###
# Sod along y direction
###

begin
    D[:matter] = "gas"
    D[:case] = "cavity"
    D[:space] = "2d2f2v"
    D[:flux] = "kfvs"
    D[:collision] = "bgk"
    D[:nSpecies] = 1
    D[:interpOrder] = 2
    D[:limiter] = "vanleer"
    D[:boundary] = "maxwell"
    D[:cfl] = 0.3
    D[:maxTime] = 5.0

    D[:x0] = 0.0
    D[:x1] = 1.0
    D[:nx] = 1
    D[:y0] = 0.0
    D[:y1] = 1.0
    D[:ny] = 100
    D[:pMeshType] = "uniform"
    D[:nxg] = 0
    D[:nyg] = 0

    D[:umin] = -5.0
    D[:umax] = 5.0
    D[:nu] = 8
    D[:vmin] = -5.0
    D[:vmax] = 5.0
    D[:nv] = 8
    D[:vMeshType] = "rectangle"
    D[:nug] = 0
    D[:nvg] = 0

    D[:knudsen] = 0.001
    D[:mach] = 0.0
    D[:prandtl] = 1.0
    D[:inK] = 1.0
    D[:omega] = 0.81
    D[:alphaRef] = 1.0
    D[:omegaRef] = 0.5

    D[:uLid] = 0.15
    D[:vLid] = 0.0
    D[:tLid] = 1.0
end

ks = SolverSet(D)
KS = ks
ctr = Array{ControlVolume2D}(
    undef,
    1,
    ks.ps.ny,
)
a1face = Array{Interface2D}(undef, KS.pSpace.nx + 1, KS.pSpace.ny)
a2face = Array{Interface2D}(undef, KS.pSpace.nx, KS.pSpace.ny + 1)

for j in axes(ctr, 2), i in axes(ctr, 1)
    if j <= KS.pSpace.ny ÷ 2
        prim = [1., 0., 0., 0.5]
        w = conserve_prim(prim, ks.gas.γ)

        ctr[i, j] = ControlVolume2D(
            KS.pSpace.x[i, j],
            KS.pSpace.y[i, j],
            KS.pSpace.dx[i, j],
            KS.pSpace.dy[i, j],
            w,
            prim,
        )
    else
        prim = [0.125, 0., 0., 0.625]
        w = conserve_prim(prim, ks.gas.γ)

        ctr[i, j] = ControlVolume2D(
            KS.pSpace.x[i, j],
            KS.pSpace.y[i, j],
            KS.pSpace.dx[i, j],
            KS.pSpace.dy[i, j],
            w,
            prim,
        )
    end
end

for i = 1:KS.pSpace.nx
    for j = 1:KS.pSpace.ny
        a2face[i, j] =
            Interface2D(KS.pSpace.dx[i, j], 0.0, 1.0, KS.ib.wL)
    end
    a2face[i, KS.pSpace.ny+1] =
        Interface2D(KS.pSpace.dx[i, KS.pSpace.ny], 0.0, 1.0, KS.ib.wL)
end

t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(0.2 ÷ dt) + 1
res = zero(ks.ib.wL)

@showprogress for iter = 1:nt
    #reconstruct!(ks, ctr)
    #evolve!(ks, ctr, a1face, a2face, dt; mode = Symbol(ks.set.flux), bc = Symbol(ks.set.boundary))

    # horizontal flux
    #=@inbounds Threads.@threads for j = 1:ks.pSpace.ny
        for i = 2:ks.pSpace.nx
            flux_gks!(
                a1face[i, j].fw,
                ctr[i-1, j].w .+ ctr[i-1, j].sw[:, 1] .* ks.ps.dx[i-1, j]/2,
                ctr[i, j].w .- ctr[i, j].sw[:, 1] .* ks.ps.dx[i, j]/2,
                ks.gas.K,
                ks.gas.γ,
                ks.gas.μᵣ,
                ks.gas.ω,
                dt,
                ks.ps.dx[i-1, j]/2,
                ks.ps.dx[i, j]/2,
                a1face[i, j].len,
                ctr[i-1, j].sw[:, 1],
                ctr[i, j].sw[:, 1],
            )
        end
    end=#

    # vertical flux
    @inbounds Threads.@threads for j = 2:ks.pSpace.ny
        for i = 1:ks.pSpace.nx
            wL = local_frame(ctr[i, j-1].w, 0., 1.)
            wR = local_frame(ctr[i, j].w, 0., 1.)
            swL = local_frame(ctr[i, j-1].sw[:, 2], 0., 1.)
            swR = local_frame(ctr[i, j].sw[:, 2], 0., 1.)
#=
            flux_gks!(
                a2face[i, j].fw,
                wL .+ swL .* ks.ps.dy[i, j-1]/2,
                wR .- swR .* ks.ps.dy[i, j]/2,
                ks.gas.K,
                ks.gas.γ,
                ks.gas.μᵣ,
                ks.gas.ω,
                dt,
                ks.ps.dy[i, j-1]/2,
                ks.ps.dy[i, j]/2,
                a2face[i, j].len,
                swL,
                swR,
            )=#

            flux_roe!(a2face[i, j].fw, wL, wR, ks.gas.γ, dt)

            a2face[i, j].fw .= KitBase.global_frame(a2face[i, j].fw, 0., 1.) .* ctr[i, j].dx
        end
    end

    #update!(ks, ctr, a1face, a2face, dt, res; coll = Symbol(ks.set.collision), bc = Symbol(ks.set.boundary))

    @inbounds Threads.@threads for j = 2:ks.pSpace.ny-1
        for i = 1:ks.pSpace.nx
            #@. ctr[i, j].w += (a1face[i, j].fw - a1face[i+1, j].fw + a2face[i, j].fw - a2face[i, j+1].fw) / ctr[i, j].dx / ctr[i, j].dy
            @. ctr[i, j].w += (a2face[i, j].fw - a2face[i, j+1].fw) / ctr[i, j].dx / ctr[i, j].dy
            ctr[i, j].prim .= conserve_prim(ctr[i, j].w, ks.gas.γ)
        end
    end

    t += dt
end

begin
    sol = zeros(ks.ps.ny, 4)
    for i in axes(sol, 1)
        sol[i, :] .= ctr[1, i].prim
        sol[i, 4] = 1 / sol[i, 4]
    end

    plot(ks.ps.y[:], sol[:, :])
end

###
# Sod along x direction
###

begin
    D[:matter] = "gas"
    D[:case] = "cavity"
    D[:space] = "2d2f2v"
    D[:flux] = "kfvs"
    D[:collision] = "bgk"
    D[:nSpecies] = 1
    D[:interpOrder] = 2
    D[:limiter] = "vanleer"
    D[:boundary] = "maxwell"
    D[:cfl] = 0.3
    D[:maxTime] = 5.0

    D[:x0] = 0.0
    D[:x1] = 1.0
    D[:nx] = 100
    D[:y0] = 0.0
    D[:y1] = 1.0
    D[:ny] = 1
    D[:pMeshType] = "uniform"
    D[:nxg] = 0
    D[:nyg] = 0

    D[:umin] = -5.0
    D[:umax] = 5.0
    D[:nu] = 8
    D[:vmin] = -5.0
    D[:vmax] = 5.0
    D[:nv] = 8
    D[:vMeshType] = "rectangle"
    D[:nug] = 0
    D[:nvg] = 0

    D[:knudsen] = 0.001
    D[:mach] = 0.0
    D[:prandtl] = 1.0
    D[:inK] = 1.0
    D[:omega] = 0.81
    D[:alphaRef] = 1.0
    D[:omegaRef] = 0.5

    D[:uLid] = 0.15
    D[:vLid] = 0.0
    D[:tLid] = 1.0
end

ks = SolverSet(D)
KS = ks
ctr = Array{ControlVolume2D}(
    undef,
    ks.ps.nx,
    1,
)
a1face = Array{Interface2D}(undef, KS.pSpace.nx + 1, KS.pSpace.ny)

for j in axes(ctr, 2), i in axes(ctr, 1)
    if i <= KS.pSpace.nx ÷ 2
        prim = [1., 0., 0., 0.5]
        w = conserve_prim(prim, ks.gas.γ)

        ctr[i, j] = ControlVolume2D(
            KS.pSpace.x[i, j],
            KS.pSpace.y[i, j],
            KS.pSpace.dx[i, j],
            KS.pSpace.dy[i, j],
            w,
            prim,
        )
    else
        prim = [0.125, 0., 0., 0.625]
        w = conserve_prim(prim, ks.gas.γ)

        ctr[i, j] = ControlVolume2D(
            KS.pSpace.x[i, j],
            KS.pSpace.y[i, j],
            KS.pSpace.dx[i, j],
            KS.pSpace.dy[i, j],
            w,
            prim,
        )
    end
end

for i = 1:KS.pSpace.nx
    for j = 1:KS.pSpace.ny
        a1face[i, j] =
            Interface2D(KS.pSpace.dx[i, j], 1.0, 0.0, KS.ib.wL)
    end
    a1face[ks.ps.nx+1, 1] =
        Interface2D(KS.pSpace.dx[i, 1], 1.0, 0.0, KS.ib.wL)
end

t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(0.2 ÷ dt) + 1
res = zero(ks.ib.wL)

@showprogress for iter = 1:nt
    #reconstruct!(ks, ctr)
    #evolve!(ks, ctr, a1face, a2face, dt; mode = Symbol(ks.set.flux), bc = Symbol(ks.set.boundary))

    # horizontal flux
    @inbounds Threads.@threads for j = 1:ks.pSpace.ny
        for i = 2:ks.pSpace.nx
            #=flux_gks!(
                a1face[i, j].fw,
                ctr[i-1, j].w .+ ctr[i-1, j].sw[:, 1] .* ks.ps.dx[i-1, j]/2,
                ctr[i, j].w .- ctr[i, j].sw[:, 1] .* ks.ps.dx[i, j]/2,
                ks.gas.K,
                ks.gas.γ,
                ks.gas.μᵣ,
                ks.gas.ω,
                dt,
                ks.ps.dx[i-1, j]/2,
                ks.ps.dx[i, j]/2,
                a1face[i, j].len,
                ctr[i-1, j].sw[:, 1],
                ctr[i, j].sw[:, 1],
            )=#

            flux_roe!(a1face[i, j].fw, ctr[i-1, j].w, ctr[i, j].w, ks.gas.γ, dt)
            a1face[i, j].fw .*= ctr[i, j].dy
        end
    end

    #update!(ks, ctr, a1face, a2face, dt, res; coll = Symbol(ks.set.collision), bc = Symbol(ks.set.boundary))

    @inbounds Threads.@threads for j = 1:ks.pSpace.ny
        for i = 2:ks.pSpace.nx-1
            #@. ctr[i, j].w += (a1face[i, j].fw - a1face[i+1, j].fw + a2face[i, j].fw - a2face[i, j+1].fw) / ctr[i, j].dx / ctr[i, j].dy
            @. ctr[i, j].w += (a1face[i, j].fw - a1face[i+1, j].fw) / ctr[i, j].dx / ctr[i, j].dy
            ctr[i, j].prim .= conserve_prim(ctr[i, j].w, ks.gas.γ)
        end
    end

    t += dt
end

begin
sol = zeros(ks.ps.nx, 4)
for i in axes(sol, 1)
    sol[i, :] .= ctr[i, 1].prim
    sol[i, 4] = 1 / sol[i, 4]
end

plot(ks.ps.x[:], sol[:, :])
end

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
"""
    parse_author(node)::NamedTuple

This parses the email, name, and ID of the author node.
"""
parse_author(node) = (email = escape_string(node.email),
                      name = escape_string(node.name),
                      id = isnothing(node.user) ? missing : escape_string(node.user.id))
"""
    parse_commit(branch, node)::NamedTuple

This parses a commit node and adds the branch it queried.
"""
function parse_commit(branch, node)
    if isnothing(node)
        println(branch)
        throw(ErrorException("Weird thing going on"))
    end
    authors = parse_author.(getproperty.(node.authors.edges, :node))
    emails = getproperty.(authors, :email)
    names = getproperty.(authors, :name)
    users = getproperty.(authors, :id)
    (branch = branch,
     id = node.id,
     sha1 = node.oid,
     committed_ts = replace(node.committedDate, "Z" => ""),
     emails = emails,
     names = names,
     users = users,
     additions = node.additions,
     deletions = node.deletions)
end
"""
    query_commits(branch::AbstractString)::Nothing
"""
function query_commits(branch::AbstractString)::Nothing
    @unpack conn, schema = GHOST.PARALLELENABLER
    since = execute(conn, "SELECT MIN(committedat) AS since FROM $(schema).commits WHERE branch = '$branch';") |>
        (obj -> only(getproperty.(obj, :since)))
    since = coalesce(since, GHOST.GH_FIRST_REPO_TS)
    output = DataFrame(vcat(fill(String, 4), fill(Vector{Union{Missing,String}}, 3), fill(Int, 2)),
                       [:branch, :id, :sha1, :committed_ts, :emails, :names, :users, :additions, :deletions],
                       0)
    query = String(read(joinpath(pkgdir(GHOST), "src", "assets", "graphql", "04_commits_single.graphql"))) |>
        (obj -> replace(obj, r"\s+" => " ")) |>
        (obj -> replace(obj, r"\s+(\{|\}|\:)\s*" => s"\1")) |>
        (obj -> replace(obj, r"(:|,|\.{3})\s*" => s"\1")) |>
        strip |>
        string
    vars = Dict("since" => string(since, "Z"),
                "until" => "2020-01-01T00:00:00Z",
                "node" => branch,
                "first" => 32,
                )
    success = false
    json = try
        while !success
            result = graphql(query, vars = vars, max_retries = 1)
            json = JSON3.read(result.Data)
            if haskey(json, :errors)
                if first(json.errors).type == "NOT_FOUND"
                    execute(conn, "UPDATE $schema.repos SET status = 'NOT_FOUND' WHERE branch = '$branch';")
                    return
                end
            end
            try
                json = json.data.node.target.history
                success = !isempty(json.edges)
            catch err
                vars["first"] == 1 && throw(ErrorException("$branch is not playing nice."))
                vars["first"] ÷= 2
                sleep(0.25)
            end
        end
        json
    catch err
        throw(ErrorException("$branch is not playing nice ($first)."))
    end
    for edge in json.edges
        push!(output, GHOST.parse_commit(branch, edge.node))
    end
    execute(conn, "BEGIN;")
    load!(output,
          conn,
          string("INSERT INTO $(schema).commits VALUES (",
                 join(("\$$i" for i in 1:size(output, 2)), ','),
                 ") ON CONFLICT ON CONSTRAINT commits_pkey DO NOTHING;"))
    execute(conn, "COMMIT;")
    while json.pageInfo.hasNextPage
        sleep(0.25)
        vars["until"] = string(DateTime(output[end,:committed_ts]), "Z")
        vars["first"] = 64
        success = false
        json = try
            while !success
                result = graphql(query, vars = vars, max_retries = 1)
                json = JSON3.read(result.Data)
                if haskey(json, :errors)
                    if first(json.errors).type == "NOT_FOUND"
                        execute(conn, "UPDATE $schema.repos SET status = 'NOT_FOUND' WHERE branch = '$branch';")
                        return
                    end
                end
                try
                    json = json.data.node.target.history
                    success = !isempty(json.edges)
                catch err
                    vars["first"] == 1 && throw(ErrorException("$branch is not playing nice."))
                    vars["first"] ÷= 2
                    sleep(0.25)
                end
            end
            json
        catch err
            throw(ErrorException("$branch is not playing nice ($first)."))
        end
        output = DataFrame(vcat(fill(String, 4), fill(Vector{Union{Missing,String}}, 3), fill(Int, 2)),
                       [:branch, :id, :sha1, :committed_ts, :emails, :names, :users, :additions, :deletions],
                       0)
        for edge in json.edges
            push!(output, GHOST.parse_commit(branch, edge.node))
        end
        execute(conn, "BEGIN;")
        load!(output,
              conn,
              string("INSERT INTO $schema.commits VALUES (",
                     join(("\$$i" for i in 1:size(output, 2)), ','),
                     ") ON CONFLICT ON CONSTRAINT commits_pkey DO NOTHING;"))
        execute(conn, "COMMIT;")
    end
    execute(conn,
            """
            UPDATE $(schema).repos
            SET status = 'Done'
            WHERE branch = '$branch'
            ;
            """)
    println("$branch done at $(now())")
    nothing
end
"""
    query_commits(branches::AbstractVector{<:AbstractString}, batch_size::Integer)::Nothing
"""
function query_commits(branches::AbstractVector{<:AbstractString}, batch_size::Integer)::Nothing
    @unpack conn, schema = PARALLELENABLER
    output = DataFrame(vcat(fill(String, 4), fill(Vector{Union{Missing,String}}, 3), fill(Int, 2)),
                       [:branch, :id, :sha1, :committed_ts, :emails, :names, :users, :additions, :deletions],
                       0)
    query = String(read(joinpath(pkgdir(GHOST), "src", "assets", "graphql", "03_commits.graphql"))) |>
        (obj -> replace(obj, r"\s+" => " ")) |>
        (obj -> replace(obj, r"\s+(\{|\}|\:)\s*" => s"\1")) |>
        (obj -> replace(obj, r"(:|,|\.{3})\s*" => s"\1")) |>
        strip |>
        string
    vars = Dict("until" => "$(floor(now(), Year))Z",
                "nodes" => branches,
                "first" => batch_size)
    result = graphql(query, vars = vars)
    json = try
        json = JSON3.read(result.Data)
        json.data
    catch err
        if length(branches) == 1
            execute(conn, "UPDATE $(schema).repos SET status = 'FOR_LATER' WHERE branch = '$(only(branches))';")
        else
            query_commits(view(branches, 1:length(branches) ÷ 2), batch_size)
            query_commits(view(branches, length(branches) ÷ 2 + 1:lastindex(branches)), batch_size)
        end
        return
    end
    for (branch, nodes) in zip(branches, values(json.nodes))
        if isnothing(nodes)
            execute(conn, "UPDATE $(schema).repos SET status = 'NOT_FOUND' WHERE branch = '$branch';")
        else
            for edge in nodes.target.history.edges
                nodes = values(edge)
                if any(isnothing, nodes)
                    execute(conn, "UPDATE $(schema).repos SET status = 'SERVICE_UNAVAILABLE' WHERE branch = '$branch';")
                else
                    for node in nodes
                        push!(output, parse_commit(branch, node))
                    end
                end
            end
        end
    end
    try
        execute(conn, "BEGIN;")
        load!(output,
              conn,
              string("INSERT INTO $(schema).commits VALUES (",
                     join(("\$$i" for i in 1:size(output, 2)), ','),
                     ") ON CONFLICT ON CONSTRAINT commits_pkey DO NOTHING;"))
        execute(conn, "COMMIT;")
    catch err
        println(branches)
        throw(err)
    end
    execute(conn,
            """
            UPDATE $(schema).repos
            SET status = 'Done'
            WHERE branch = ANY('{$(join(unique(output.branch), ','))}'::text[])
            ;
            """)
    nothing
end

ext: jl


text:
# -------------- Implementation of Product Basis

struct Product1pBasis{NB, TB <: Tuple, NSYM, SYMS, T} <: OneParticleBasis{T}
   bases::TB
   spec::Vector{NamedTuple{SYMS, NTuple{NSYM, Int}}}
   indices::Vector{NTuple{NB, Int}}
   _typeT::Type{T}
end

function Product1pBasis(bases;
                        SYMS = _symbols_prod(bases),
                        T = promote_type(fltype.(bases)...) )
   # TODO: discuss whether to construct an optimal ordering, e.g.
   #       should the discrete bases come first once we implement the
   #       "strongzero" method?
   NSYM = length(SYMS)
   NB = length(bases)
   Product1pBasis( tuple(bases...),
                   NamedTuple{SYMS, NTuple{NSYM, Int}}[],
                   NTuple{NB, Int}[],
                   T )
end

import Base.*
*(B1::OneParticleBasis, B2::OneParticleBasis) =
      Product1pBasis((B1, B2))
*(B1::Product1pBasis, B2::OneParticleBasis) =
      Product1pBasis((B1.bases..., B2))
*(B1::OneParticleBasis, B2::Product1pBasis) =
      Product1pBasis((B1, B2.bases...))
*(B1::Product1pBasis, B2::Product1pBasis) =
      Product1pBasis((B1.bases..., B2.bases...))

_numb(b::Product1pBasis{NB}) where {NB} = NB

Base.length(basis::Product1pBasis) = length(basis.spec)

fltype(basis::Product1pBasis) = promote_type(fltype.(basis.bases)...)

alloc_temp(basis::Product1pBasis, args...) =
      (
         B = alloc_B.(basis.bases),
         tmp = alloc_temp.(basis.bases)
      )

alloc_temp_d(basis::Product1pBasis, args...) =
      (
         B = alloc_B.(basis.bases),
         tmp = alloc_temp.(basis.bases),
         dB = alloc_dB.(basis.bases),
         tmpd = alloc_temp_d.(basis.bases)
      )

function gradtype(basis::Product1pBasis)
   # get the types of the sub-bases ...
   B = alloc_B.(basis.bases)
   dB = alloc_dB.(basis.bases)
   # ... do some artificial arithmetic  ...
   x = dB[1][1] * prod( B[a][1] for a = 2:length(B) )
   for a = 2:length(B)
      y = dB[a][1]
      for b = 1:length(B)
         b == a && continue
         y *= B[b][1]
      end
      x += y
   end
   # ... to find out the end-result of the gradient calculation
   return typeof(x)
end

@generated function add_into_A!(A, tmp, basis::Product1pBasis{NB}, X) where {NB}
   quote
      Base.Cartesian.@nexprs $NB i -> evaluate!(tmp.B[i], tmp.tmp[i], basis.bases[i], X)
      for (iA, ϕ) in enumerate(basis.indices)
         t = one(eltype(A))
         Base.Cartesian.@nexprs $NB i -> (t *= tmp.B[i][ϕ[i]])
         A[iA] += t
      end
      return nothing
   end
end

@generated function add_into_A_dA!(A, dA, tmpd, basis::Product1pBasis{NB}, X
                                   ) where {NB}
   quote
      Base.Cartesian.@nexprs($NB, i -> begin   # for i = 1:NB
         evaluate_ed!(tmpd.B[i], tmpd.dB[i], tmpd.tmpd[i], basis.bases[i], X)
      end)
      for (iA, ϕ) in enumerate(basis.indices)
         # evaluate A
         t = one(eltype(A))
         Base.Cartesian.@nexprs($NB, i -> begin   # for i = 1:NB
            t *= tmpd.B[i][ϕ[i]]
         end)
         A[iA] += t

         # evaluate dA
         dA[iA] = zero(eltype(dA))
         Base.Cartesian.@nexprs($NB, a -> begin  # for a = 1:NB
            dt = tmpd.dB[a][ϕ[a]]
            Base.Cartesian.@nexprs($NB, b -> begin  # for b = 1:NB
               if b != a
                  dt *= tmpd.B[b][ϕ[b]]
               end
            end)
            dA[iA] += dt
         end)
      end
      return nothing
   end
end

_symbols_prod(bases) = tuple(union( symbols.(bases)... )...)

symbols(basis::Product1pBasis{NB, TB, NSYM, SYMS}
            ) where {NB, TB, NSYM, SYMS} = SYMS

function indexrange(basis::Product1pBasis)
   allsyms = tuple(symbols(basis)...)
   rg = Dict{Symbol, Vector{Int}}([ sym => Int[] for sym in allsyms]...)
   for b in basis.bases
      rgb = indexrange(b)
      for sym in allsyms
         if haskey(rgb, sym)
            rg[sym] = union(rg[sym], rgb[sym])
         end
      end
   end
   # HACK: fix the m range based on the maximal l-range
   #       this needs to be suitably generalised if we have multiple
   #       (l, m) pairs, e.g. (l1, m1), (l2, m2)
   if haskey(rg, :m)
      maxl = maximum(rg[:l])
      rg[:m] = collect(-maxl:maxl)
   end

   # convert the range into a named tuple so that we remember the order!!
   return NamedTuple{allsyms}(ntuple(i -> rg[allsyms[i]], length(allsyms)))
end

isadmissible(b, basis::Product1pBasis) = all(isadmissible.(Ref(b), basis.bases))

function set_spec!(basis::Product1pBasis{NB}, spec) where {NB}
   empty!(basis.spec)
   append!(basis.spec, spec)
   empty!(basis.indices)
   for b in basis.spec
      inds = ntuple(i -> get_index(basis.bases[i], b), NB)
      push!(basis.indices, inds)
   end
   return basis
end

get_spec(basis::Product1pBasis) = basis.spec

get_spec(basis::Product1pBasis, i::Integer) = basis.spec[i]

degree(b, basis::Product1pBasis) = sum( degree(b, B) for B in basis.bases )

# TODO: this looks like a horrible hack ...
function rand_radial(basis::Product1pBasis)
   for B in basis.bases
      if B isa ScalarACEBasis
         return rand_radial(B)
      end
   end
   return nothing
end


text:
using DataFrames
using DelimitedFiles
using Glob
using CSV

function joinDFS(dfs)
  tsv = dfs[1];
  for i=2:length(dfs)
    tsv = innerjoin(tsv,dfs[i][:,1:end .!=2],on=:barcode,makeunique=true);
  end
  return tsv;
end

const re = Regex("CRE.{1}");

files = glob("*.counts"); #get all counts files
dfs = [DataFrame() for _ in 1:length(files)]; #might be faster just to push

if length(files) > 1
  #tables = [Array() for _ in 1:length(files)];

  for i=1:length(files)  #turn files into matricies/tables
    tsv = readdlm(files[i]);
    dfs[i] = DataFrame(barcode=tsv[:,1],name=tsv[:,2],count=tsv[:,3]);
    m = match(re,files[i]);
    mString = m.match;
    rename!(dfs[i],:count => mString);
  end

  #tsv = dfs[1];
  #for i=2:length(dfs)
  #  tsv = innerjoin(tsv,dfs[i][:,1:end .!=2],on=:barcode,makeunique=true);
  #end
  tsv = joinDFS(dfs);
  CSV.write("combined.tsv",tsv,delim="\t");

#for i=1:length(files)
#  tsv = readdlm(files[i]);
#  dfs[i] = DataFrame(barcode=tsv[:,1],name=tsv[:,2],count=tsv[:,3]);
#end

end

max_stars_repo_name: HenryCWong/genomic-barcode-counter

max_issues_repo_name: HenryCWong/genomic-barcode-counter

max_forks_repo_name: HenryCWong/genomic-barcode-counter


text:
using Random
using LinearAlgebra
using Statistics
using Optim
using DataFrames
using DataFramesMeta
using CSV
using HTTP
using GLM

# 1. read in and reshape
cd("C:/Users/jakek/Documents/econ_phd/structural_modeling_course/forked_class/ProblemSets/PS5-ddc/")
include("create_grids.jl")

# transform to long.
df = CSV.read("busdata.csv")
df = @transform(df, bus_id = 1:size(df,1))
dfy = @select(df, :bus_id,:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20,:RouteUsage,:Branded)
dfy_long = DataFrames.stack(dfy, Not([:bus_id,:RouteUsage,:Branded]))
rename!(dfy_long, :value => :Y)
dfy_long = @transform(dfy_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfy_long, Not(:variable))

# next reshape the odometer variable
dfx = @select(df, :bus_id,:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
dfx_long = DataFrames.stack(dfx, Not([:bus_id]))
rename!(dfx_long, :value => :Odometer)
dfx_long = @transform(dfx_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfx_long, Not(:variable))

# join reshaped df's back together
df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
sort!(df_long,[:bus_id,:time])

Zst = convert(Vector, df[:,:Zst])
B = convert(Vector, df[:,:Branded])
X = df[:, 43:62]

## 2. Run flexible logit
beta_glm = glm(@formula(Y ~ Odometer*Odometer*RouteUsage*RouteUsage*Branded*time*time ), df_long, Binomial(), LogitLink())
println(beta_glm)

### 3
## a. state transitions
zval,zbin,xval,xbin,xtran = create_grids()

## b.
Helper = DataFrame(Odometer = kron(ones(zbin),xval),RouteUsage = kron(ones(xbin),zval), Branded = vec(zeros(size(xtran,1),1)),time = vec(zeros(size(xtran,1),1)))
beta=0.9

function ccp_estim(glm_est, helper, zval,zbin,xval,xbin,xtran,beta,df_long,df, B, Zst,X)
	FV = repeat(zeros(size(xtran,1),2 ),1,1,21)
	FVT1 = zeros(size(df,1),20 )
	for t=2:20
		for b=0:1
			helper[:,:time] = t .*(ones(size(helper,1)))
			helper[:,:Branded] = b.* ones(size(helper,1))
			p0 = 1 .- predict(glm_est,helper)
			FV[:,b+1,t] = -beta .* log.(p0)
		end
	end

	for i=1:size(df,1)
		for t=1:20
			row0 = 1+(Zst[i]-1)*xbin
			row1 = X[i,t]+(Zst[i]-1)*xbin
			FVT1[i,t] = (xtran[row1,:] .- xtran[row0,:])' * FV[row0:row0+xbin-1,B[i]+1,t+1]
		end
	end
	FVT1 = @transform(convert(DataFrame, FVT1), bus_id = 1:size(df,1))
	FVT1 = DataFrames.stack(FVT1, Not([:bus_id]))
	FVT1 = @transform(FVT1, time = kron(collect([1:20]...),ones(size(df,1))))
	sort!(FVT1,[:bus_id,:time])
	return FVT1
end

fvt1 = ccp_estim(beta_glm, Helper, zval,zbin,xval,xbin,xtran,beta,df_long,df,B,Zst, X)

fv=convert(Vector, fvt1[:,:value])
df_long = @transform(df_long,fv=fv )

## final estimation

theta_hat_ccp_glm = glm(@formula(Y ~ Odometer + Branded), df_long, Binomial(), LogitLink(), offset=df_long.fv)


text:
#This file will contain some demo functions for DDPG algorithms:
include("src/init.jl")# Include all relevant functions
include("src/DDPG.jl") #Contains Deep deterministic Policy Gradient functions

"""
This file contains the script to execute the DDPG algorithm:
https://arxiv.org/abs/1509.02971
"""

res = zeros(params.n_iter)
prec = zeros(params.n_iter)

#1. Randomly Initialize
#-------------------------------------------------------------
#1.1 Actor
θa  = init_θ_single_norm(params.structure_net_pl)
θa_mod = copy(θa)
opt_para_a = optimizers(θa, Adam; lr = params.lr_pl)
ĝ_a = copy(θa)

#1.2 Critic
Q  = init_θ_single_norm(params.structure_net_qf)
Q_mod = copy(Q)
opt_para_q = optimizers(Q, Adam; lr = params.lr_qf)
ĝ_v = copy(Q)

#1.3 Replay Buffer and Latest Trajectories
buffer = Traj_Data()
buffer.traj = Array{Trajectory,1}(undef,0)
buffer.timesteps = 0
buffer.n_traj = 0

all_traj  = init_all_traj(params)

#main function running DDPG
function DDPG()
    weights_old = zeros(length(net_to_vec(θa)),2)
    a_t_old = all_traj.traj[1].a_t

    #1.4 Init running mean std
    rms_all = init_rms(params.s_init)

    #Iterate
    for i_iter = 1:params.n_iter
        #Init Progress bar
        progress = Progress(params.n_iter,1)

        if floor(i_iter/100) == ceil(i_iter/100)
            #make_video(all_traj,params, i_iter)
            save_data(res, prec, θa, Q, all_traj, i_iter)
        end

        #2.Collect Data using current stochastic Policy
        #-------------------------------------------------------------
        res[i_iter], prec[i_iter] = parallel_episodes(params,θa,θa,all_traj,rms_all)
        #(ĝ_a, ĝ_v, res[i_iter], prec[i_iter]) = parallel_episodes_skeleton(params,θa,θa,all_traj,rms_all)
        prec[i_iter] = loss_Q(Q,Q_mod,θa,θa_mod,all_traj.traj[1],params)

        #3.Pump data in Replay Buffer:
        #-------------------------------------------------------------
        #3.1 Generate minibatches from current data_batch
        mini_batch = build_minibatches(all_traj,params)
        #3.2 Pump the data into the buffer
        pump_buffer!(buffer,mini_batch,params)

        #4.Train all the stuff now:
        #-------------------------------------------------------------
        println(" ")
        @time parallel_DDPG_update(buffer,Q,Q_mod,θa,θa_mod,opt_para_a,opt_para_q,params)

        #5.Visualize Progress
        ProgressMeter.update!(progress,i_iter) #update progress bar
        if floor(i_iter/5) == ceil(i_iter/5)
            println("new visualization: ")
            a_t_old ,weights_old = plot_stuff(params, all_traj, a_t_old,
            weights_old, i_iter,θa, Q, res, prec/params.n_traj, ĝ_a, ĝ_v)
        end
        println(res[i_iter])
    end
end

DDPG()

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
module Dev

using Distributions
include("grad.jl")

mutable struct VDNormal{F <: AbstractFloat}
  m::F
  s::F
  log_s::F
  VDNormal(m::F, s::F) where F = new{F}(m, s, log(s))
end

rand(d::VDNormal) = randn() * d.s + d.m
grad_logq(d::VDNormal, x) = grad_logq_normal(d.m, d.log_s, x)

function update!(d::VDNormal, grads, lr)
  g_m, g_log_s = grads
  d.m += lr * g_m
  d.log_s += lr * g_log_s
  d.s = exp(d.log_s)
  return d
end

end # Dev

#= Test
vd = Dev.VDNormal(1., .5)
x = Dev.rand(vd)
grads = Dev.grad_logq(vd, x)
Dev.update!(vd, grads, .01)
=#


text:
using Test,TightBinding
using Plots
using LinearAlgebra

@testset "Real to k" begin
    @testset "1D" begin
        @time energies = TightBinding.test_1D()
        sumt = 26.21913117088203
        @test round(sum(abs.(energies))*100) == round(sumt*100)
    end
    @testset "2Dsquare" begin
        sumt = 52.438262341764066
        @test round(sum(abs.(TightBinding.test_2Dsquare()))*100) == round(sumt*100)
    end
    @testset "Graphene" begin
        @time hamk = TightBinding.test_2DGraphene()
        sumt = 5.975026051875501
        @test round(sum(abs.(hamk))*100) == round(sumt*100)
#        @time test_Graphene()
    end
    @testset "Iron pnictides" begin
        @time hamk =  TightBinding.test_pnictides()
        sumt = 2.680985878019015
        @test round(sum(abs.(hamk))*100) == round(sumt*100)

        @time hamk = TightBinding.test_pnictides_5orbitals()
        sumt = 3.0141792528265308
        @test round(sum(abs.(hamk))*100) == round(sumt*100)

        #@time test_Fe2()
        #@time test_Fe5()
    end

end

@testset "k to real" begin
    @testset "Topological insulator" begin
        @time mat_e = TightBinding.test_surface()
        sumt = 17011.010972794524
        @test round(sum(abs.(mat_e))*100) == round(sumt*100)
    end
end


text:
##############################################################
function softparzen(X, S_centre, r, γ = 1e-6)
##############################################################

  N, D = size(X)
  M = length(S_centre)

  # Usually kernel has also some normalisation constant
  # but in the following calculations it cancels out anyway
  K = zeros(M,N)
  for j in 1:M
    for n in 1:N
      @inbounds K[j,n] = exp(-0.5*norm(X[S_centre[j],:] - X[n,:])^2 / (r*r))
    end
  end

  # Soft version for coefficients
  Q = zeros(M)
  sumK = sum(K)
  for j in 1:M
    @inbounds Q[j] = sum(K[j,:]) / sumK
  end

  # responsibilities - note different to standard mixture models
  # See equation (6)
  resp = zeros(M,N)
  for j in 1:M
    sum_Kj_over_n = sum(K[j,:])
    for n in 1:N
      @inbounds resp[j,n] = K[j,n] / sum_Kj_over_n
    end
  end

  # Soft version for means
  # See equation (7)
  mu = zeros(M, D)
  for j=1:M
    for n=1:N
      for i=1:D
        @inbounds mu[j,i] += resp[j,n] * X[n,i]
      end
    end
  end

  # Soft version for covariances
  C = [zeros(D, D) for j in 1:M]

  for j in 1:M
    for i in 1:D
      for d in 1:D
        for n in 1:N
          @inbounds C[j][i,d] += resp[j,n] * ((X[n,i] - mu[j,i]) * (X[n,d] - mu[j,d]))
        end
      end
    end

    C[j] = (C[j] + C[j]')*0.5 + γ*I

  end

  return Q, mu, C

end


text:
# TODO: this probably should be removed and all the functionality moved into `sextic.jl`

module Elasticity_110

using JuLIP: AbstractAtoms, AbstractCalculator, calculator,
         stress

using JuLIPMaterials.CLE: elastic_moduli, voigt_moduli

const Tensor{T} = Array{T,4}

function four_to_two_index(i, j)

if (i == 1 && j == 1)
   dumb = 1
elseif (i == 2 && j ==2)
   dumb = 2
elseif (i == 3 && j == 3)
   dumb = 3
elseif ((i == 3 && j == 2) || (i == 2 && j == 3))
   dumb = 4
elseif ((i == 3 && j == 1) || (i == 1 && j == 3))
   dumb = 5
elseif ((i == 2 && j == 1) || (i == 1 && j == 2))
   dumb = 6
end

return dumb

end

function two_to_four_index(k)

it = 0
jit = 0

if k == 1
   it = 1;
   jit = 1;
elseif k == 2
   it = 2;
   jit = 2;
elseif k == 3
   it = 3;
   jit = 3;
elseif k == 4
   it = 3;
   jit = 2;
elseif k == 5
   it = 3;
   jit = 1;
elseif k == 6
   it = 2;
   jit = 1;
end

return it, jit

end

function little_a(D::Array{T,2},r,s) where {T}

   C = zeros(3,3,3,3)
   Chat = zeros(3,3,3,3)

   #Convert back to Tensor notation
   for i=1:3, j = 1:3, k = 1:3, l = 1:3
   	m = four_to_two_index(i,j)
        n = four_to_two_index(k,l)
        C[i,j,k,l] = D[m,n]
   end

   #Rotate the tensor to correct orientation
   Tr = [1/sqrt(2) -1/sqrt(2) 0; 0 0 1; 1/sqrt(2) 1/sqrt(2) 0]
   Q = zeros(3,3,3,3)
   for i=1:3, j=1:3, k=1:3, l=1:3
	Q[i,j,k,l] = Tr[k,i]*Tr[l,j]
   end

   for i=1:3, j=1:3, k=1:3, l=1:3, g=1:3, h=1:3, m=1:3, n=1:3
	Chat[i,j,k,l] = Chat[i,j,k,l] + Q[g,h,i,j]*C[g,h,m,n]*Q[m,n,k,l]
   end

   #M = zeros(3,1)
   #M(1,1) = Chat[r,1,s,1]
   #M(2,1) = Chat[r,1,s,2] + Chat[r,2,s,1]
   #M(3,1) = Chat[r,2,s,2]
   #return M(1,1), M(2,1), M(3,1)
   p0 = Chat[r,1,s,1]
   p1 = Chat[r,1,s,2] + Chat[r,2,s,1]
   p2 = Chat[r,2,s,2]
   return p0, p1, p2
end

function fourth_order_basis(D::Array{T,2},a;
            Tr = [1/sqrt(2) -1/sqrt(2) 0; 0 0 1; 1/sqrt(2) 1/sqrt(2) 0]
				) where {T}
   C = zeros(3,3,3,3)
   Chat = zeros(3,3,3,3)

   #Convert back to Tensor notation
   for i=1:3, j = 1:3, k = 1:3, l = 1:3
   	m = four_to_two_index(i,j)
        n = four_to_two_index(k,l)
        C[i,j,k,l] = D[m,n]
   end

   #Rotate the tensor to correct orientation
   Q = zeros(3,3,3,3)
   for i=1:3, j=1:3, k=1:3, l=1:3
	Q[i,j,k,l] = Tr[k,i]*Tr[l,j]
   end

   for i=1:3, j=1:3, k=1:3, l=1:3, g=1:3, h=1:3, m=1:3, n=1:3
	Chat[i,j,k,l] = Chat[i,j,k,l] + Q[g,h,i,j]*C[g,h,m,n]*Q[m,n,k,l]
   end

   M = zeros(6,6)
   #Convert the tensor back to 6 by 6
   for i=1:6, j=1:6
        m, n = two_to_four_index(i)
        p, q = two_to_four_index(j)
        M[i,j] = Chat[m,n,p,q]
   end

  return M
end

function A_coefficients(p::Array{Complex{Float64},1},D::Array{T,2}) where {T}

  A = Complex{Float64}[0 0 0; 0 0 0; 0 0 0]
  x = Complex{Float64}[0; 0 ; 0]
  y = Complex{Float64}[0; 0 ; 0]
  w = Complex{Float64}[0; 0 ; 0]
  z = Complex{Float64}[0; 0 ; 0]

  for i=1:3
    x[i] = D[1,4]^2 - D[4,4]*D[5,5]-(D[4,4]^2+D[2,2]*D[5,5])*( (real(p[i]))^2 - (imag(p[i]))^2) - D[4,4]*D[2,2]*( ((real(p[i]))^2 - (imag(p[i]))^2)^2 - 4*(real(p[i]))^2*(imag(p[i]))^2  )   #real B'
    y[i] = -2*(D[2,2]*D[5,5]+D[4,4]^2)*real(p[i])*imag(p[i])-4*D[2,2]*D[4,4]*real(p[i])*imag(p[i])*( (real(p[i]))^2 - (imag(p[i]))^2)    #image B'
    w[i] = 2*D[2,2]*D[1,4]*real(p[i])*( (real(p[i]))^2 - (imag(p[i]))^2)+D[1,4]*(D[4,4]-D[1,2])*real(p[i])-4*D[1,4]*D[2,2]*real(p[i])*(imag(p[i]))^2   #real B''
    z[i] = 4*D[1,4]*D[2,2]*(real(p[i]))^2*imag(p[i])+2*D[2,2]*D[1,4]*imag(p[i])*( (real(p[i]))^2 - (imag(p[i]))^2)  +D[1,4]*(D[4,4]-D[1,2])*imag(p[i])  #image B''
  end

  for i=1:3
    A[1,i] = (x[i]*w[i]+y[i]*z[i])/(w[i]^2+z[i]^2)+ ((y[i]*w[i]-x[i]*z[i])/ (w[i]^2+z[i]^2))*im
    A[2,i] = 2*(imag(p[i])*imag(A[1,i]) - real(A[1,i])*real(p[i]))-(D[5,5]+D[4,4]*( (real(p[i]))^2 - (imag(p[i]))^2)  )/D[4,4] - (2*D[4,4]*real(p[i])*imag(p[i]) + 2*D[1,4]*(real(p[i])*imag(A[1,i])-imag(p[i])*real(A[1,i]) )  )*(1/D[1,4])*im
    A[3,i] = 1+0*im
  end

  return A

end

function D_coefficients(p::Array{Complex{Float64},1},D::Array{T,2},
						      A::Array{Complex{Float64},2}, b) where {T}

  alpha = zeros(6,6)
  v = zeros(6,1)
  v[1] = b #first three components of v are burgers vector
  for i=1:3
    for j=1:6
      l = ceil(j/2)
      l = convert(Int,l)
      if mod(j,2) == 0
        alpha[i,j] = -imag(A[i,l])
      else
        alpha[i,j] = real(A[i,l])
      end
    end
  end

  for j=1:3
    k = 2*j-1
    m = 2*j
    l = ceil(k/2)
    l = convert(Int,l)
    alpha[4,k] = D[6,6]*(real(A[1,l])*real(p[l])-imag(A[1,l])*imag(p[l])+real(A[2,l]) ) + D[5,6]
    alpha[5,k] = D[1,2]*real(A[1,l])+D[2,2]*(real(A[2,l])*real(p[l]) - imag(A[2,l])*imag(p[l])  )
    alpha[6,k] = D[1,4]*real(A[1,l])+D[4,4]*real(p[l])

    alpha[4,m] = -D[6,6]*(real(A[1,l])*imag(p[l])+imag(A[1,l])*real(p[l])+imag(A[2,l])  )
    alpha[5,m] = -D[1,2]*imag(A[1,l])-D[2,2]*(real(A[2,l])*imag(p[l])+imag(A[2,l])*real(p[l]))
    alpha[6,m] = -D[1,4]*imag(A[1,l])-D[4,4]*imag(p[l])
  end
  D = \(alpha,v)
  return D
end

function sextic_roots(D::Array{T,2}) where {T}

  inter = D[1,1]*D[2,2]-2*D[1,2]*D[6,6]-D[1,2]^2
  lead = D[2,2]*D[4,4]*D[6,6]
  k_0 = D[5,5]*D[6,6]*D[1,1]/lead

  k_2 = (D[4,4]*D[1,1]*D[6,6] + D[5,5]*inter)/lead

  k_4 = (D[4,4]*inter+D[5,5]*D[2,2]*D[6,6])/lead

#Compute the roots p^2 = r of the sextic polynomial using general solution of cubic
  Q =  (3*k_2 - k_4^2)/9
  R =  (9*k_2*k_4-27*k_0-2*k_4^3)/54
  E =  Q^3 + R^2
  #print("E: ")
  #print(E)
  S = cbrt(R+sqrt(E))
  U = cbrt(R-sqrt(E))
  r_1 = -1/3*k_4 + (S + U) + 0*im
  r_2 = -1/3*k_4-1/2*(S+U)+1/2*im*sqrt(3)*(S-U)
  r_3 = -1/3*k_4-1/2*(S+U)-1/2*im*sqrt(3)*(S-U)
#Now compute the roots p = \pm sqrt(r) making sure to take those with positive imaginary part
  p = Complex{Float64}[0; 0 ; 0]

  p[1] = sqrt(r_1)
  p[2] = sqrt(r_2)
  p[3] = sqrt(r_3)

  for i=1:3
    if imag(p[i]) < 0
      p[i] = -p[i]
    end
  end

  return p
end

end


text:
# CUDA Initialization and Context Management

export context, context!, device!, device_reset!

## initialization

"""
    CUDAnative.prepare_cuda_call()

Prepare state for calling CUDA API functions.

Many CUDA APIs, like the CUDA driver API used by CUDAnative, use implicit thread-local state
to determine, e.g., which device to use. With Julia however, code is grouped in tasks.
Execution can switch between them, and tasks can be executing on (and in the future migrate
between) different threads. To synchronize these two worlds, call this function before any
CUDA API call to update thread-local state based on the current task and its context.

If you need to maintain your own thread-local state, subscribe to context and task switch
events using [`CUDAnative.atcontextswitch`](@ref) and [`CUDAnative.attaskswitch`](@ref) for
proper invalidation.
"""
@inline function prepare_cuda_call()
    tid = Threads.threadid()
    task = current_task()

    # detect when a different task is now executing on a thread
    if @inbounds thread_tasks[tid] != task
        switched_tasks(tid, task)
    end

    # initialize a CUDA context when first executing on a thread
    if @inbounds thread_contexts[tid] === nothing
        initialize_thread(tid)
    end

    check_exceptions()

    return
end

# the default device unitialized tasks will use, set when switching devices.
# this behavior differs from the CUDA Runtime, where device 0 is always used.
# this setting won't be used when switching tasks on a pre-initialized thread.
const default_device = Ref{Union{Nothing,CuDevice}}(nothing)

# CUDA uses thread-bound contexts, but calling CuCurrentContext all the time is expensive,
# so we maintain our own thread-local state keeping track of the current context.
const thread_contexts = Union{Nothing,CuContext}[]
@noinline function initialize_thread(tid::Int)
    ctx = CuCurrentContext()
    if ctx === nothing
        dev = something(default_device[], CuDevice(0))
        device!(dev)
    else
        # compatibility with externally-initialized contexts
        thread_contexts[tid] = ctx
    end

end

# Julia executes with tasks, so we need to keep track of the active task for each thread
# in order to detect task switches and update the thread-local state accordingly.
# doing so using task_local_storage is too expensive.
const thread_tasks = Union{Nothing,WeakRef}[]
@noinline function switched_tasks(tid::Int, task::Task)
    thread_tasks[tid] = WeakRef(task)
    _attaskswitch(tid, task)

    # switch contexts if task switched to was already bound to one
    ctx = get(task_local_storage(), :CuContext, nothing)
    if ctx !== nothing
        context!(ctx)
    end
    # NOTE: deactivating the context in the case ctx===nothing would be more correct,
    #       but that confuses CUDA and leads to invalid contexts later on.
end

"""
    CUDAnative.attaskswitch(f::Function)

Register a function to be called after switching tasks on a thread. The function is passed
two arguments: the thread ID, and the task switched to.

Use this hook to invalidate thread-local state that depends on the current task.
"""
attaskswitch(f::Function) = (pushfirst!(task_hooks, f); nothing)
const task_hooks = []
_attaskswitch(tid, task) = foreach(f->Base.invokelatest(f, tid, task), task_hooks)

## context-based API

"""
    context()::CuContext

Get or create a CUDA context for the current thread (as opposed to
`CUDAdrv.CuCurrentContext` which may return `nothing` if there is no context bound to the
current thread).
"""
@inline function context()
    tid = Threads.threadid()

    prepare_cuda_call()
    ctx = @inbounds thread_contexts[tid]::CuContext

    if Base.JLOptions().debug_level >= 2
        @assert ctx == CuCurrentContext()
    end

    ctx
end

"""
    context!(ctx::CuContext)

Bind the current host thread to the context `ctx`.

Note that the contexts used with this call should be previously acquired by calling
[`context`](@ref), and not arbirary contexts created by calling the `CuContext` constructor.
"""
function context!(ctx::CuContext)
    # update the thread-local state
    tid = Threads.threadid()
    thread_ctx = @inbounds thread_contexts[tid]
    if thread_ctx != ctx
        thread_contexts[tid] = ctx
        activate(ctx)
        _atcontextswitch(tid, ctx)
    end

    # update the task-local state
    task_local_storage(:CuContext, ctx)

    return
end

"""
    context!(f, ctx)

Sets the active context for the duration of `f`.
"""
function context!(f::Function, ctx::CuContext)
    old_ctx = CuCurrentContext()
    try
        context!(ctx)
        f()
    finally
        if old_ctx != nothing
            context!(old_ctx)
        end
    end
end

"""
    CUDAnative.atcontextswitch(f::Function)

Register a function to be called after switching contexts on a thread. The function is
passed two arguments: the thread ID, and the context switched to.

If the new context is `nothing`, this indicates that the context is being unbound from this
thread (typically during device reset).

Use this hook to invalidate thread-local state that depends on the current device or context.
"""
atcontextswitch(f::Function) = (pushfirst!(context_hooks, f); nothing)
const context_hooks = []
_atcontextswitch(tid, ctx) = foreach(f->Base.invokelatest(f, tid, ctx), context_hooks)

## device-based API

"""
    device!(dev::Integer)
    device!(dev::CuDevice)

Sets `dev` as the current active device for the calling host thread. Devices can be
specified by integer id, or as a `CuDevice` (slightly faster).

Although this call is fairly cheap (50-100ns), it is only intended for interactive use, or
for initial set-up of the environment. If you need to switch devices on a regular basis,
work with contexts instead and call [`context!`](@ref) directly (5-10ns).

If your library or code needs to perform an action when the active context changes,
add a hook using [`CUDAnative.atcontextswitch`](@ref).
"""
function device!(dev::CuDevice, flags=nothing)
    tid = Threads.threadid()

    # configure the primary context
    pctx = CuPrimaryContext(dev)
    if flags !== nothing
        @assert !isactive(pctx) "Cannot set flags for an active device. Do so before calling any CUDA function, or reset the device first."
        CUDAdrv.cuDevicePrimaryCtxSetFlags(dev, flags)
    end

    # bail out if switching to the current device
    if @inbounds thread_contexts[tid] !== nothing && dev == device()
        return
    end

    # have new threads use this device as well
    default_device[] = dev

    # activate a context
    ctx = CuContext(pctx)
    context!(ctx)
end
device!(dev::Integer, flags=nothing) = device!(CuDevice(dev), flags)

"""
    device!(f, dev)

Sets the active device for the duration of `f`.
"""
function device!(f::Function, dev::CuDevice)
    old_ctx = CuCurrentContext()
    try
        device!(dev)
        f()
    finally
        if old_ctx != nothing
            context!(old_ctx)
        end
    end
end
device!(f::Function, dev::Integer) = device!(f, CuDevice(dev))

"""
    device_reset!(dev::CuDevice=device())

Reset the CUDA state associated with a device. This call with release the underlying
context, at which point any objects allocated in that context will be invalidated.
"""
function device_reset!(dev::CuDevice=device())
    pctx = CuPrimaryContext(dev)
    ctx = CuContext(pctx)

    # unconditionally reset the primary context (don't just release it),
    # as there might be users outside of CUDAnative.jl
    unsafe_reset!(pctx)

    # wipe the context handles for all threads using this device
    for (tid, thread_ctx) in enumerate(thread_contexts)
        if thread_ctx == ctx
            thread_contexts[tid] = nothing
            _atcontextswitch(tid, nothing)
        end
    end

    return
end


text:
using DataStructures
using AxisArrays

import Base: getindex, keys, convert, show, IO, MIME

"""
AxisArrays provides the equivalent of an xarray `DataArray`, this wrapper
adds the functionality of an xarray `DataSet` by containing multiple
`AxisArray`s.

Object should be created by calling the relevant
`convert($(FUNCTIONNAME), dataset::DataSet)` method on a [`DataSet`](@ref DataSet).

# See also

[`Backend`](@ref Backend), [`ArrayWrapper`](@ref ArrayWrapper)
"""
struct AxisArrayWrapper <: ArrayWrapper
    "Dictionary of `DimensionName::String => DimensionLength::Int`"
    dimensions::OrderedDict
    "Named tuple of `DatasetName::Symbol => Dataset::AxisArray`"
    datasets::T where T <: NamedTuple
    "Dictionary of `AttributeName::String => AttributeValue::Any`"
    attributes::OrderedDict
    "Dictionary containing encoding information (usually `source`, `filter_by_keys`, and `ecode_cf`)"
    encoding::Dict

    #  Manually define inner constructor here so that it does not appear twice
    #  in the docs
    AxisArrayWrapper(dimensions, datasets, attributes, encoding) = new(dimensions, datasets, attributes, encoding)
end

getindex(obj::AxisArrayWrapper, key) = getfield(obj, :datasets)[key]
keys(obj::AxisArrayWrapper) = keys(getfield(obj, :datasets))
haskey(obj::AxisArrayWrapper, key) = key in keys(obj)

function Base.getproperty(obj::AxisArrayWrapper, key::Symbol)
    if key in keys(obj)
        return getindex(obj, key)
    else
        return getfield(obj, key)
    end
end

function convert(::Type{AxisArrayWrapper}, dataset::DataSet)::AxisArrayWrapper
    dimensions = dataset.dimensions
    attributes = dataset.attributes
    encoding = dataset.encoding

    multidimensional_idx =  (
        [size(v.data) for v in values(dataset.variables)]
        .|> length
        .|> x -> x > 1)
    multidimensional_keys = collect(keys(dataset.variables))[multidimensional_idx]
    multidimensional_values = [dataset.variables[k] for k in multidimensional_keys]

    shared_dimensions = [dataset.variables[k] for k in keys(dataset.dimensions)]
    shared_axis = [Axis{Symbol(k)}(dataset.variables[k].data) for k in keys(dataset.dimensions)]

    datasets = NamedTuple{Tuple(Symbol.(multidimensional_keys))}((
        AxisArray(CfGRIB.convert(Array, dataset.variables[k].data), shared_axis...)
        for k in multidimensional_keys
    ))

    return AxisArrayWrapper(dimensions, datasets, attributes, encoding)
end

convert(::Type{AxisArray}, dataset::DataSet)::AxisArrayWrapper = convert(AxisArrayWrapper, dataset)

function Base.show(io::IO, mime::MIME"text/plain", da::AxisArrayWrapper)
    dimensions_list = join(["$k: $v" for (k, v) in pairs(da.dimensions)], ", ")
    str_dimensions = " Dimensions ($(length(da.dimensions))):\n  $dimensions_list"

    dataset_list = [summary(ds) for ds in da.datasets]
    axes_list = split(dataset_list[1], "\n")[2:end - 1]
    axes_list = replace.(axes_list, "    " => "  ")
    #  When arrays are printed their length is limited by the current displaysize
    #  so adding text before/after will make the text overflow to the next line
    #  here we do a hacky fix for that
    display_width = displaysize(io)[2]
    for (i, line) in enumerate(axes_list)
        if length(line) <= display_width
            continue
        end

        new_line = split(line, ",")

        #  dot dot dot (…) index - index of pair of values that has ellipses
        #  e.g. `15.0 … 342.0`
        ddd_idx = findfirst(occursin.("…", split(line, ",")))
        while sum(length.(new_line)) + length(new_line) + 3 > display_width
            deleteat!(new_line, ddd_idx)
            str_new_ellipses = "$(new_line[ddd_idx - 1])  … $(new_line[ddd_idx])"
            deleteat!(new_line, ddd_idx - 1)
            new_line[ddd_idx - 1] = str_new_ellipses
            ddd_idx = ddd_idx - 1
        end

        axes_list[i] = join(new_line, ",")
    end

    axes_list = join(axes_list, "\n")
    str_axes = " Axes: \n$axes_list"

    dataset_data = [split(dl, "\n")[end] for dl in dataset_list]
    dataset_data = [
        "   " * replace(s, "And data, a" => string(keys(da.datasets)[i]) * ",")
        for (i, s)
        in enumerate(dataset_data)
    ]
    dataset_variables = join(dataset_data, "\n")

    str_variables = " Variables:\n$dataset_variables"

    str_show = join([str_dimensions, str_axes, str_variables], "\n")

    println(io, "AxisArrayWrapper with $(length(da.datasets)) dataset(s)")

    println(str_show)
    print(" Attributes: ")
    show(io, mime, da.attributes)
end


text:
using ApproxFun, SpectralMeasures, Test, LinearAlgebra
    import ApproxFun: ldiv_coefficients
############
### Tests
############

@testset "Spectral measures" begin
    @testset "Chebyshev U" begin
        a = Float64[]; b= Float64[]
        @time μ=spectralmeasure(a,b)
        @test μ.(-.99:.01:.99) ≈ sqrt.(1 .-(-.99:.01:.99).^2)*2/π
    end

    @testset "Chebyshev T" begin
        a = [0.]; b=[1/sqrt(2),.5]
        @time μ=spectralmeasure(a,b)
        @test μ.(-.99:.01:.99) ≈ sqrt.(1 .-(-.99:.01:.99).^2).^(-1)/π
    end

    @testset "Legendre" begin
        n=100; a=zeros(n); b=(1:n-1)./sqrt.(4*(1:n-1).^2 .-1)
        μ=spectralmeasure(a,b)
        @test μ.(-.99:.01:.99) ≈ 0.5one.(-.99:.01:.99) atol=0.001
    end

    @testset "A strange hump" begin
        n=10;a=zeros(n); b=sqrt.(1:(n-1))./(2sqrt(n))
        @time μ=spectralmeasure(a,b)
        @test μ(0.1) ≈ 2.031454229315879 atol=1E-9 # empirical
    end

    @testset "Jacobi polynomials" begin
        α = 1.2; β = 1.1
        n=107;a = @. (β^2-α^2)/((2(0:n)+α+β)*(2*(1:n+1)+ (α+β)))
        b = @. 2sqrt(((1:n)*(α+(1:n))*(β+(1:n))*(α+β+(1:n)))/
                ((2(1:n)+α+β-1)*((2(1:n)+α+β)^2)*(2*(1:n) + α+β+1)))
        μ=spectralmeasure(a,b)
        x=Fun()
        ν=(1-x)^α*(x+1)^β
        ν/=sum(ν)

        @test μ.(-.99:.01:.99) ≈ ν.(-.99:.01:.99) atol=0.001
    end

    @testset "Simplest perturbation" begin
        k = 17
        a = [k/20];b=Float64[]
        μ=spectralmeasure(a,b)

        k = -19
        a = [k/20];b=Float64[]
        ν=spectralmeasure(a,b)
    end

    @testset "Check robustness for extra zeros" begin
        for n=3:6,m=0:2
            L=freejacobioperator()
            K=SymTriOperator(-[ones(5);zeros(m)],zeros(n))

            J=L+0.5K
            @test isa(J,SymTriPertToeplitz)
            @test Matrix(J[1:3,1:6]) == [-0.5  0.5  0   0 0 0;
                                        0.5 -0.5  0.5 0 0 0;
                                        0    0.5 -0.5 0.5 0 0]

            A=J-3.0*I
            Q,L=ql(A)

            @test norm((Q*L-A)[1:30,1:30]) < 100eps()

            @test isa(L*Q,SpectralMeasures.SymTriPertToeplitz)
        end
    end
end

@testset "eig ops" begin
    L=freejacobioperator()
    K=SymTriOperator(-ones(5),zeros(4))

    J=L+0.5K

    D,Q=eigen(J)

    r = ldiv_coefficients(Q, mul_coefficients(Q,[1.0]))
    @test r ≈ [1;zeros(length(r)-1)]
    @test coefficients(Q\(Q*[1.0])) ≈ [1;zeros(length(r)-1)]

    @test Q\(Q*[1.0]) ≈ Fun([1;zeros(length(r)-1)],domainspace(Q))
    @test Q\(D*(Q*[1.0])) ≈ J*[1.0]
end

@testset "zero b" begin
    J=SymTriPertToeplitz(ones(0),Float64[],0.0,0.5) + 1.01I
    D,Q = eigen(J)

    @test Q\(D*(Q*[1.0])) ≈ J*[1.0]

    J = -freejacobioperator()
    D,Q = eigen(J)

    @test Q\(D*(Q*[1.0])) ≈ J*[1.0]

    J=SymTriPertToeplitz(ones(0),Float64[],0.0,-0.5) + 1.01I

    Q,L = ql(J)
    @test norm((Q*L-J)[1:100,1:100]) ≤ 100eps()

    @time D,Q = eigen(J)
    @test Q\(D*(Q*[1.0])) ≈ J*[1.0]
end

@testset "Zeros in QL (#19)" begin
    J = SymTriPertToeplitz([2.0], Float64[0.0000000001], -2, 0.5)
    Q,L = ql(J)
    @test (Q*L)[1:10,1:10] ≈ J[1:10,1:10]

    J = SymTriPertToeplitz([2.0], Float64[0.0], -2, 0.5)
    Q,L = ql(J)
    @test (Q*L)[1:10,1:10] ≈ J[1:10,1:10]
end

@testset "Diagonal operator functions" begin
    Δ = freejacobioperator()
    K=SymTriOperator(-2ones(3),zeros(4))
    J = Δ + K
    D,U = eigen(J)

    v = randn(10)
    @test pad((U\(exp(D)*U*v)).coefficients,100) ≈ exp(Matrix(J[1:100,1:100]))*pad(v,100)
end


text:
# String with Length
# This allows to compute length once and only once
struct StringWithLength{T<:AbstractString} <: AbstractString
    s::T
    l::Int
end
string_with_length(s::AbstractString) = StringWithLength(s, length(s))
string_with_length(s::StringWithLength) = s
Base.length(s::StringWithLength) = s.l
Base.iterate(s::StringWithLength, i::Integer = firstindex(s.s)) = iterate(s.s, i)
Base.isequal(s1::StringWithLength, s2::AbstractString) = isequal(s.s1, s2)
Base.isequal(s1::AbstractString, s2::StringWithLength) = isequal(s1, s2.s)
Base.nextind(s::StringWithLength, i::Int, n::Int = 1) = nextind(s.s, i, n)
Base.ncodeunits(s::StringWithLength) = ncodeunits(s.s)
Base.isvalid(s::StringWithLength, i::Int) = isvalid(s.s, i)
function reorder(s1::AbstractString, s2::AbstractString)
    s1 = string_with_length(s1)
    s2 = string_with_length(s2)
    if length(s1) <= length(s2)
         return s1, s2
    else
        return s2, s1
    end
end

## Find common prefixes (up to lim. -1 means Inf)
function remove_prefix(s1::AbstractString, s2::AbstractString, lim::Integer = -1)
    l = 0
    x1 = iterate(s1)
    x2 = iterate(s2)
    while (x1 !== nothing) & (x2 !== nothing) & (l < lim || lim < 0)
        ch1, state1 = x1
        ch2, state2 = x2
        ch1 != ch2 && break
        x1 = iterate(s1, state1)
        x2 = iterate(s2, state2)
        l += 1
    end
    return l, x1, x2
end

# Return start of commn substring in s1, start of common substring in s2, and length of substring
# Indexes refer to character number, not index
function longest_common_substring(s1::AbstractString, s2::AbstractString, len1::Integer, len2::Integer)
    if len1 > len2
        start2, start1, len = longest_common_substring(s2, s1, len2, len1)
    else
        start1, start2, len = 0, 0, 0
        p = zeros(Int, len2)
        i1 = 0
        for ch1 in s1
            i1 += 1
            i2 = 0
            oldp = 0
            for ch2 in s2
                i2 += 1
                newp = 0
                if ch1 == ch2
                    newp = oldp > 0 ? oldp : i2
                    currentlength = (i2 - newp + 1)
                    if currentlength > len
                        start1, start2, len = i1 - currentlength + 1, newp, currentlength
                    end
                end
                p[i2], oldp = newp, p[i2]
            end
        end
    end
    return start1, start2, len
end


text:
function render_all()
    return dbc_container([
        dbc_card([
            dbc_cardheader(
                dbc_tabs(
                    [
                        render_data_tab(),
                        render_nn_tab(),
                        render_prediction_tab()
                    ],
                    id="card-tabs",
                    card=true,
                    active_tab="tab-1",
                ),
            )
        ])
    ],
    fluid=true
    )
end
function render_info()
  return dbc_jumbotron([
          html_h1("Dash Flux App", className="dflx-display-3"),
          html_p(
              "Web app to generate, train, test neural networks",
              className="lead",
          ),
          html_hr(className="my-2"),
          html_p("Users can upload data, select paramteres, create neurl network layour,
          trace the NN training and predict using new data")
      ])
end
function render_data_tab()
  return dbc_tab([
      dbc_cardbody([
          html_div([
            render_info(),
              dbc_row([
                dbc_col(md=3,[
                  html_div(
                    id="dflx-def-upload-data",
                    [
                      html_label("Select a default dataset:")
                      dcc_dropdown(
                        id="dflx-def-upload-data-drop",
                        options =[
                            Dict("label" => "Select Data", "value" => "0"),
                            Dict("label" => "Real Estate", "value" => "1"),
                            Dict("label" => "Solar", "value" => "2"),
                            Dict("label" => "House Price", "value" => "3")
                        ],
                        value="0"
                        )
                    ],
                  )
                ]),
                dbc_col([
                  html_h3("OR")
                ], md=1),
                dbc_col(md=6, [
                  dcc_upload(
                    id="dflx-upload-data",
                    children=html_div([
                        "Drag and Drop or ",
                        html_a("Select Data Files")
                    ]),
                    style=Dict(
                        "width" => "100%",
                        "height" => "60px",
                        "lineHeight" => "60px",
                        "borderWidth" => "1px",
                        "borderStyle" => "dashed",
                        "borderRadius" => "5px",
                        "textAlign" => "center",
                        "margin" => "10px"
                    ),
                    # Allow multiple files to be uploaded
                    multiple=false
                )
                ])
              ]),
              dbc_row([
                  dbc_col([html_div(id="dflx-output-data-upload")], md=6),
                  dbc_col([html_div(id="dflx-output-plot")], md=6),
              ],
              align="center",
              )

          ])
      ]),
  ],
  id="dflx-data-tab", label="Data", tab_id="tab-1")

end
function render_nn_tab()
  dbc_tab([
      html_div([
          dbc_row([
              dbc_col([
                  dbc_card([
                      dbc_cardbody(
                          [
                              html_h4("Select Inputs", className="card-title"),
                              html_div(id="dflx-input-name-list"),
                          ]
                      )
                  ])
                  ],
                  md=2
              ),
              dbc_col([
                  dbc_card([
                      dbc_cardbody(
                          [
                              html_h4("Hidden Layers", className="card-title"),
                              html_div([
                                  dbc_row([
                                      dbc_col([html_i(id="dflx-add-hidden-layer",className="fa fa-plus-square fa-2x")],md=1,style=Dict("color" => "green")),
                                      dbc_col([html_i(id="dflx-remove-hidden-layer",className="fa fa-minus-square fa-2x")],md=1,style=Dict("color" => "red"))
                                  ], style=Dict("cursor" => "pointer"))
                              ]),
                              dcc_store(id="dflx-nn-layer-count"),
                              dbc_row(id = "dflx-hidden-layers",
                                render_nn_layer(1)),
                              html_div(id="dflx-nnlayer-count", style=Dict("display" => "none"))
                          ]
                      )
                  ])
                  ],
                  md=8
              ),
              dbc_col([
                  dbc_card([
                      dbc_cardbody(
                          [
                              html_h4("Select Outputs", className="card-title"),
                              html_div(id="dflx-output-name-list"),
                          ]
                      )
                  ])
                  ],
                  md=2
              )
          ],
          align="center",
          ),
          dbc_row([
              dbc_col([
                  daq_numericinput(id="dflx-nntrain-percent", label="Percentage of training dataset", value=70, max = 99, min = 1)
              ]),

              dbc_col([
                dcc_dropdown(
                    id="dflx-nntrain-optimizer",
                    options = [
                        (label = "Gradient Descent", value = "descent"),
                        (label = "Momentum", value = "momentum"),
                        (label = "Nesterov", value = "nesterov"),
                        (label = "RMSProp", value = "rmsprop"),
                        (label = "ADAM", value = "adam"),
                        (label = "Rectified ADAM", value = "radam"),
                        (label = "AdaMax", value = "adamax"),
                    ],
                    value = "adam",
                    placeholder="Select Optimizer",
                )
              ]),
              dbc_col([
                  dcc_input(id="dflx-nntrain-lrate", placeholder="Learning Rate", value=0.1, type="number")
              ]),
              dbc_col([
                daq_numericinput(id="dflx-nntrain-epoch", label="Epochs", value=4, min = 1, max=10000)
              ]),
              dbc_col([
                  dbc_button(
                      "Start Training", id="dflx-nnlayer-submit",
                      color="success", className="mr-1",
                      style = Dict(
                          "margin-top" => "20px",
                          "display" => "flex",
                          "justify-content" => "center"
                      )
                  ),
              ])
          ]),
      ]),
      html_div(id="dflx-training-updates")
  ],
  style = Dict("margin-top" => "20px"),
  id="dflx-nn-tab", label="NN Config", tab_id="tab-2", disabled="true")
end
function render_prediction_tab()
  return dbc_tab(
    id="dflx-pred-tab", label="Prediction",
    tab_id="tab-3", disabled="true",
    style = Dict("margin-top" => "20px"),
  )
end
function render_prediction_componnets(labels_in)
  if labels_in isa Nothing
    in_card = [
        html_div([
          dcc_input(
              id = (index = 1, type = "preddone"),
              placeholder = "nothing",
              type = "number",
              disabled=true
          )
      ])
    ]
  elseif labels_in isa Vector
      in_card = [
          html_div([
            dcc_input(
                id = (index = i, type = "preddone"),
                placeholder = val,
                type = "number"
            )
        ]) for (i,val) in enumerate(labels_in)
      ]
  elseif labels_in isa String
      in_card = [
            html_div([
              dcc_input(
                  id = (index = 1, type = "preddone"),
                  placeholder = labels_in,
                  type = "number"
              )
          ])
        ]
  end
  return dbc_row([
      dbc_col([
        dbc_card(
            dbc_cardbody([
                html_h5("Input", className="card-title"),
                html_div(
                  in_card
                )
            ]),
        )
      ],md=4),
      dbc_col([
        dbc_card(
            dbc_cardbody([
                html_h5("Neural network", className="card-title")
            ]),
        )
      ],md=4),
      dbc_col([
        dbc_card(
            dbc_cardbody([
                html_h5("Output", className="card-title"),
                html_div(
                  id="pred_out"
                )
            ]),
        )
      ],md=4)
  ], align="center")
end
function render_table(df, filename; n_rows = 1)
  df = df[!, Not("id")]
  return html_div([
    html_h3(filename),
    DashTable.dash_datatable(
      data = Dict.(pairs.(eachrow(df))),
      columns=[Dict("name" =>c, "id" => c) for c in names(df)],
      virtualization=true,
      fixed_rows=Dict("headers" => true),
      style_cell=Dict(
          "minWidth" =>  95, "maxWidth" =>  95, "width" =>  95
      ),
      style_table=Dict("height" => 300),  # default is 500
      page_action="none"
  )])
end
function render_graph_title(df)
  row_names = names(df)
  row_names = row_names[row_names .!= "id"]
  return html_div([
    html_h5("Plots"),
    dbc_row([
      dbc_col([
        render_single_graphtitle(row_names)
      ],
      md=6),
      dbc_col([
        render_mult_graphtitle(row_names)
      ],
      md=6)
    ],
    align="center"
    ),
    dcc_graph(id="dflx-plot-graph", figure =  render_graph(df, row_names[end], row_names[end]))
  ])
end
function render_graph(df, x_xolumn, y_columns)
  if !(y_columns isa Vector)
    dt = [Dict(
      "x" => df[!,Symbol(x_xolumn)],
      "y" => df[!,Symbol(y_columns)],
      "type" => "scatter",
      "mode" => "markers"
    )]
  else
        dt = [Dict(
      "x" => df[!,Symbol(x_xolumn)],
      "y" => df[!,Symbol(yc)],
      "type" => "scatter",
      "mode" => "markers"
    ) for yc in y_columns]
  end
  return Dict(
      "data" => dt,
      "layout" => Dict(
          "title" => "Dash Data Visualization"
      )
  )
end
function render_mult_graphtitle(row_names;id_val="dflx-graph-dropdown-y")
  return dcc_dropdown(
    id=id_val,
    options = [Dict("label" => rn, "value" => rn) for rn in row_names],
    value = row_names[end],
    multi=true
  )
end
function render_single_graphtitle(row_names;id_val="dflx-graph-dropdown-x")
  return dcc_dropdown(
    id=id_val,
    options = [Dict("label" => rn, "value" => rn) for rn in row_names],
    value = row_names[end],
    multi=false
  )
end

function render_nn_layer(n_click)
  return [dbc_col([
        dbc_card([
            dbc_cardbody(
                [
                    html_h4("Layer $i", className="card-title"),
                    daq_numericinput(label="Output count", value=2, max = 10000, min = 1),
                    dcc_dropdown(
                        id="dflx-nntrain-activation",
                        options = [
                            (label = "No Activation", value = "no"),
                            (label = "ReLU", value = "relu"),
                            (label = "σ", value = "sigmoid"),
                            (label = "softmax", value = "softmax"),
                            (label = "tanh", value = "tanh"),
                            (label = "RReLU", value = "rrelu"),
                            (label = "CeLU", value = "celu"),
                            (label = "eLU", value = "elu"),
                            (label = "GeLU", value = "gelu")
                        ],
                        value="no",
                        placeholder="Select Activation",
                    )
                ]
            )
        ])
        ],
        md=3
    ) for i in 1:n_click]

end
function render_training()
  global x_graph
  global y_graph
  return dbc_row([
      dbc_col([
          dcc_interval(
              id="interval-component",
              interval=1000, # in milliseconds
              n_intervals=0,
              disabled=false
          ),
          dcc_graph(id="live-update-graph",
            figure=Dict(
                "data" => [
                    Dict("x" => x_graph, "y" => y_graph),
                ]
            )
          )
      ]),
      dbc_col([
         html_div(id="dflx-testing-section")
      ]),
  ])
end
function render_testing(m, X_test,y_test)
  yₚ, maerr, mserr, crenpy = (nothing,nothing,nothing,nothing)
  try
    yₚ, maerr, mserr, crenpy = test_nn(m,X_test,y_test)
  catch e
  end
  row1 = html_tr([html_td("Mean Absolute Error"), html_td(maerr)])
  row2 = html_tr([html_td("Mean Squared Error"), html_td(mserr)])
  row3 = html_tr([html_td("Crossentropy"), html_td(crenpy)])
  table_header = [html_thead(html_tr([html_th("Measure"), html_th("Value")]))]

  table_body = [html_tbody([row1, row2, row3])]
  return html_div([
    html_h3("Test Results")
    dbc_table([table_header ; table_body], bordered=true)
  ])
end

max_stars_repo_name: DEICIDIA/fluxdash

max_issues_repo_name: DEICIDIA/fluxdash

max_forks_repo_name: DEICIDIA/fluxdash


text:
module HHModelCompoenents

function inf_tau_h_rm(V)
    i = (1+exp((V+65)/6))^(-1);
    t = 100*(7*exp((V+60)/11)+10*exp(-(V+60)/25))^(-1)+0.6;
    i, t
end

function inf_tau_m_rm(V)
    i = (1+exp(-(V+38)/7))^(-1);
    t = 10*(5*exp((V+60)/18)+36*exp(-(V+60)/25))^(-1)+0.04;
    i, t
end

function inf_tau_n_rm(V)
    a, b = alpha_beta_n(V);
    i = a / (a+b);
    t = 1 / (a+b);
    i, t
end

function alpha_beta_n(V)
    if V != 60
        a = -0.01*(V+60)/(exp(-(V+60)/10)-1);
    else
        a = 0.1;
    end
    b = 0.125*exp(-(V+70)/80);
    a, b
end

end