text:
module TestRDA
using Test
using DataFrames
using CategoricalArrays
using RData

@testset "Loading RData files (version=$ver)" for ver in (2, 3)
    rdata_path = joinpath(dirname(@__FILE__), "data_v$ver")

    @testset "Reading minimal RData" begin
        df = DataFrame(num = [1.1, 2.2])
        min_rda = load(joinpath(rdata_path, "minimal.rda"), convert=false)
        rdf = min_rda["df"]
        @test rdf isa RData.RList
        @testset "class() and inherits()" begin
            # not SEXP
            @test_throws MethodError RData.class(5)
            @test_throws MethodError RData.inherits(5, "number")
            @test_throws MethodError RData.inherits(5, "number")
            @test_throws MethodError RData.inherits(5, ["number"])

            rnotobj = RData.RBuiltin("test") # not a ROBJ
            @inferred RData.class(rnotobj)
            @inferred RData.inherits(rnotobj, "dummy")
            @test RData.class(rnotobj) === RData.emptystrvec
            @test !RData.inherits(rnotobj, "dummy")

            @inferred RData.class(rdf)
            @inferred RData.inherits(rdf, "data.frame")
            @inferred RData.inherits(rdf, ["data.frame"])
            @test RData.class(rdf) == ["data.frame"]
            @test RData.inherits(rdf, "data.frame")
            @test RData.inherits(rdf, ["data.frame"])

            rnumvec = rdf.data[1]
            @test rnumvec isa RData.RNumericVector
            @test RData.class(rnumvec) != ["data.frame"]
            @test !RData.inherits(rnumvec, "data.frame")
            @test !RData.inherits(rnumvec, ["data.frame"])
        end
        @test sexp2julia(min_rda["df"]) == df
        @test load(joinpath(rdata_path, "minimal.rda"), convert=true)["df"] == df
        @test load(joinpath(rdata_path, "minimal_ascii.rda"))["df"] == df
    end

    @testset "Conversion to Julia types" begin
        df = DataFrame(num = [1.1, 2.2],
                       int = Int32[1, 2],
                       logi = [true, false],
                       chr = ["ab", "c"],
                       factor = categorical(["ab", "c"], compress=true),
                       cplx = [1.1+0.5im, 1.0im])
        rdf = sexp2julia(load(joinpath(rdata_path, "types.rda"), convert=false)["df"])
        @test eltype.(eachcol(rdf)) == eltype.(eachcol(df))
        @test rdf == df
        rdf_ascii = sexp2julia(load(joinpath(rdata_path, "types_ascii.rda"), convert=false)["df"])
        @test eltype.(eachcol(rdf_ascii)) == eltype.(eachcol(df))
        @test rdf_ascii == df
    end

    @testset "NAs conversion" begin
        df = DataFrame(num = Union{Float64, Missing}[1.1, 2.2],
                       int = Union{Int32, Missing}[1, 2],
                       logi = Union{Bool, Missing}[true, false],
                       chr = Union{String, Missing}["ab", "c"],
                       factor = categorical(Union{String, Missing}["ab", "c"], compress=true),
                       cplx = Union{ComplexF64, Missing}[1.1+0.5im, 1.0im])

        df[2, :] .= Ref(missing)
        push!(df, df[2, :]) # add another row
        df[3, :num] = NaN
        df[!, :cplx] = [missing, ComplexF64(1, NaN), NaN]
        @test isequal(sexp2julia(load(joinpath(rdata_path, "NAs.rda"), convert=false)["df"]), df)
        @test isequal(sexp2julia(load(joinpath(rdata_path, "NAs_ascii.rda"), convert=false)["df"]), df)
    end

    @testset "Column names conversion" begin
        rda_names = names(sexp2julia(load(joinpath(rdata_path, "names.rda"), convert=false)["df"]))
        expected_names = ["_end", "x!", "x1", "_B_C_", "x", "x_1"]
        @test rda_names == expected_names
        rda_names = names(sexp2julia(load(joinpath(rdata_path, "names_ascii.rda"), convert=false)["df"]))
        @test rda_names == expected_names
    end

    @testset "Reading RDA with complex types (environments, closures etc)" begin
        rda_envs = load(joinpath(rdata_path, "envs.rda"), convert=false)
        rda_pairlists = load(joinpath(rdata_path, "pairlists.rda"), convert=false)
        rda_closures = load(joinpath(rdata_path, "closures.rda"), convert=false)
        rda_cmpfuns = load(joinpath(rdata_path, "cmpfun.rda"), convert=false)
    end

    @testset "Proper handling of factor and ordered" begin
        f = load(joinpath(rdata_path, "ord.rda"))
        @test !isordered(f["x"])
        @test levels(f["x"]) == ["a", "b", "c"]
        @test isordered(f["y"])
        @test levels(f["y"]) == ["b", "a", "c"]
        @test f["x"] == f["y"] == ["a", "b", "c"]
    end

    @testset "List of vectors (#82)" begin
        f = load(joinpath(rdata_path, "list_of_vec.rda"))
        listofvec = f["listofvec"]
        @test listofvec isa Vector{Vector{Union{Float64, Missing}}}
        @test length(listofvec) == 3
        @test isequal(listofvec, [[1., 2., missing], [3., 4.], [5., 6., missing]])
        @test listofvec[1] isa Vector{Union{Float64, Missing}}
        @test listofvec[2] isa Vector{Union{Float64, Missing}}

        listofvec2 = f["listofvec2"]
        @test listofvec2 isa Vector{Any}
        @test length(listofvec2) == 3
        @test isequal(listofvec2, [[1, 2, missing], [3., 4.], [5., 6., missing]])
        @test listofvec2[1] isa Vector{Union{Int32, Missing}}
        @test listofvec2[2] isa Vector{Float64}

        listofvec3 = f["listofvec3"]
        @test listofvec3 isa Vector{Any}
        @test length(listofvec3) == 2
        @test isequal(listofvec3, [[1, 2], [3., 4.]])
        @test listofvec3[1] isa Vector{Int32}
        @test listofvec3[2] isa Vector{Float64}

        listofvec4 = f["listofvec4"]
        @test listofvec4 isa Vector{Vector{Float64}}
        @test length(listofvec4) == 2
        @test isequal(listofvec4, [[1., 2.], [3., 4., 5.]])
        @test listofvec4[1] isa Vector{Float64}
        @test listofvec4[2] isa Vector{Float64}

        namedlistofvec = f["namedlistofvec"]
        @test namedlistofvec isa DictoVec
        @test length(namedlistofvec) == 3
        @test namedlistofvec.name2index == Dict("A"=>1, "B"=>3)
        @test isequal(values(namedlistofvec), [[1., 2., missing], [3., 4.], [5., 6., missing]])

        testdf = f["testdf"]
        @test testdf isa DataFrame
        @test nrow(testdf) == 3
        @test eltype(testdf[!, "listascol"]) === Vector{Union{Float64, Missing}}
        @test isequal(testdf[!, "listascol"], [[1., 2., missing], [3., 4.], [5., 6., missing, 7.]])
        @test testdf[!, "listascol2"] isa Vector{Any}
        @test isequal(testdf[!, "listascol2"], [[1., 2.], [3, 4], [5., 6., 7.]])
    end # list of vectors
end # for ver in ...

@testset "Loading AltRep-containing RData files (version=3)" begin
    altrep_rda = load(joinpath("data_v3", "altrep.rda"), convert=false)
    load(joinpath("data_v3", "altrep_ascii.rda"), convert=false)
    @test length(altrep_rda) == 4
    # test AltRep objects are recognized
    @test isa(altrep_rda["longseq"], RData.RAltRep)
    # TODO test that longseq is converted into Julia UnitRange
    @test isa(altrep_rda["wrapvec"], RData.RAltRep)
    @test sexp2julia(altrep_rda["wrapvec"]) == [1.0, 2.5, 3.0]
    @test isa(altrep_rda["nonnilpairlist"], RData.RAltRep)
    @test isa(altrep_rda["factoraltrep"], RData.RAltRep)
    @test isequal(sexp2julia(altrep_rda["factoraltrep"]),
                  compress(categorical(repeat(["A", "B", missing, "C"], inner=5000))))

    # test that AltRep-based column names are converted correctly
    altrep_names_rda = load(joinpath("data_v3", "altrep_names.rda"), convert=false)
    @test isa(altrep_names_rda["altrepnames_df"].attr["names"], RData.RAltRep)
    @test sexp2julia(altrep_names_rda["altrepnames_df"]) isa DataFrame
    wide_df = sexp2julia(altrep_names_rda["altrepnames_df"])
    @test names(wide_df) == ["a", "b", "c"]

    # test automatic conversion
    altrep_conv_rda = load(joinpath("data_v3", "altrep.rda"), convert=true)
    @test altrep_conv_rda["wrapvec"] == [1.0, 2.5, 3.0]
    @test isequal(altrep_conv_rda["factoraltrep"],
                  compress(categorical(repeat(["A", "B", missing, "C"], inner=5000))))
    @test isa(altrep_conv_rda["nonnilpairlist"], Matrix{Int32})
    @test size(altrep_conv_rda["nonnilpairlist"]) == (0, 10)

    @testset "compact_xxxseq" begin
        compactseq_rda = load(joinpath("data_v3", "compact_seq.rda"), convert=true)
        @test compactseq_rda["intseq1"] === 5:200
        @test compactseq_rda["intseq2"] isa Vector{<:Integer}
        @test compactseq_rda["intseq2"] == -2:2:1000
        @test compactseq_rda["realseq1"] === 110:300
        @test compactseq_rda["realseq2"] isa Vector{Float64}
        @test compactseq_rda["realseq2"] == 110.5:300.5
        @test compactseq_rda["realseq3"] isa Vector{Float64}
        @test compactseq_rda["realseq3"] == 110.0:0.5:300.0
    end
end

@testset "Duplicate levels in factor (version=3)" begin
    dup_cat = sexp2julia(load(joinpath("data_v3", "dup_levels.rda"), convert=false)["dup_levels"])
    @test dup_cat[1] == "Paced"
    @test dup_cat[2] == "Inferior"
    @test dup_cat[end] == "Anterior"
    @test levels(dup_cat) ==
        ["Inferior", "Anterior", "LBBB", "Missing", "NoSTUp", "OtherSTUp", "Paced"]
end

end # module TestRDA


text:
@safetestset "Query params" begin

@safetestset "No query params" begin
  using Genie
  using HTTP

  port = nothing
  port = rand(8500:8900)

  route("/") do
    isempty(query()) && return ""
    isempty(params(:GET)) && return ""

    "error"
  end

  server = up(port)

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test isempty(String(response.body)) == true

  down()
  sleep(1)
  server = nothing
  port = nothing
end

@safetestset "No defaults errors out" begin
  using Genie
  using HTTP

  port = nothing
  port = rand(8500:8900)

  route("/") do
    query(:a)
  end

  server = up(port)

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 500

  down()
  sleep(1)
  server = nothing
  port = nothing
end

@safetestset "Defaults when no query params" begin
  using Genie
  using HTTP

  port = nothing
  port = rand(8500:8900)

  route("/") do
    query(:x, "10") * query(:y, "20")
  end

  server = up(port)

  # ====

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test String(response.body) == "1020"

  # ====

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port/?", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test String(response.body) == "1020"

  # ====

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port/?x", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test String(response.body) == "20"

  # ====

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port/?x&a=3", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test String(response.body) == "20"

  # ====

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port/?x&y", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test isempty(String(response.body)) == true

  down()
  sleep(1)
  server = nothing
  port = nothing
end

@safetestset "Query params processing" begin
  using Genie
  using HTTP

  port = nothing
  port = rand(8500:8900)

  route("/") do
    query(:x)
  end

  server = up(port)

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port?x=1", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test String(response.body) == "1"

  # ====

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port?x=1&x=2", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test String(response.body) == "2"

  # ====

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port?x=1&x=2&x=3", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test String(response.body) == "3"

  # ====

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port?x=1&x=2&x=3&y=0", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test String(response.body) == "3"

  # ====

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port?x=0&x[]=1&x[]=2", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test String(response.body) == "0"

  down()
  sleep(1)
  server = nothing
  port = nothing
end

@safetestset "Array query params" begin
  using Genie
  using HTTP

  port = nothing
  port = rand(8500:8900)

  route("/") do
    query(:x, "10") * join(query(Symbol("x[]"), "100"))
  end

  server = up(port)

  # ====

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test String(response.body) == "10100"

  # ====

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port/?x&x[]=1000", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test String(response.body) == "1000"

  # ====

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port/?x&x[]=1000&x[]=2000", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test String(response.body) == "10002000"

  # ====

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port/?x=9&x[]=1000&x[]=2000", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test String(response.body) == "910002000"

  # ====

  response = try
    HTTP.request("GET", "http://127.0.0.1:$port/?x=9&x[]=1000&x[]=2000&y[]=8", ["Content-Type" => "text/html"])
  catch ex
    ex.response
  end

  @test response.status == 200
  @test String(response.body) == "910002000"

  down()
  sleep(1)
  server = nothing
  port = nothing
end

end

ext: jl


text:
using Documenter, Literate, AIBECS

# Ensure that the DataDeps download work remotely
ENV["DATADEPS_ALWAYS_ACCEPT"] = true

# generate tutorials and how-to guides using Literate
src = joinpath(@__DIR__, "src")
lit = joinpath(@__DIR__, "lit")
notebooks = joinpath(src, "notebooks")

execute = true # Set to true for executing notebooks and documenter!
nb = true      # Set to true to generate the notebooks
for (root, _, files) in walkdir(lit), file in files
    splitext(file)[2] == ".jl" || continue
    ipath = joinpath(root, file)
    opath = splitdir(replace(ipath, lit=>src))[1]
    Literate.markdown(ipath, opath, documenter = execute)
    nb && Literate.notebook(ipath, notebooks, execute = execute)
end

# Documentation structure
ismd(f) = splitext(f)[2] == ".md"
pages(folder) = [joinpath(folder, f) for f in readdir(joinpath(src, folder)) if ismd(f)]

makedocs(
    sitename="AIBECS.jl",
    doctest = false, # TODO guessing I should remove that when actually deploying?
    # options
    modules = [AIBECS],
    # organisation
    pages = Any[
        "Home" => "index.md",
        "Tutorials" => pages("tutorials"),
        "How-to guides" => pages("howtos"),
        "Explanation" => pages("explanation"),
        "Reference" => pages("reference")
        ]
)

# Deploy
deploydocs(
    repo = "github.com/JuliaOcean/AIBECS.jl.git",
    push_preview = true
)

#=
To edit locally, make sure execute is set to false and run

using LiveServer
servedocs(literate=joinpath("docs", "lit"), doc_env=true)

from the root of the package in development
=#


text:
function isprime(n)
    if n == 1
        return false
    elseif n == 2
        return true
    else
        for i = 2:(n ^ (0.5)) + 1
            if n % i == 0
                return false
            end
        end
        return true
    end
end

function prime_factors(n)
    count = 0
    final = [1]
    k = n
    if n == 1
        return []
    end
    while k != 1
        count += 1
        if k % count == 0 && isprime(count) == true
            k = k/count
            final = [final; [count]]
            count = 0
        end

    end
    return final[2:end]
end


text:
export Ray

## Standard interval

# angle is (false==0) and π (true==1)
# or ranges from (-1,1].  We use 1 as 1==true.
#orientation true means oriented out
immutable Ray{angle,T<:Number} <: IntervalDomain{T}
    center::T
    orientation::Bool
    Ray(c,o) = new(c,o)
    Ray(c) = new(c,true)
    Ray() = new(zero(T),true)
    Ray(r::Ray{angle,T}) = r
end

typealias RealRay{T} Union{Ray{false,T},Ray{true,T}}

@compat (::Type{Ray{a}}){a}(c,o) = Ray{a,typeof(c)}(c,o)
@compat (::Type{Ray{a}}){a}(c::Number) = Ray{a,typeof(c)}(c)
@compat (::Type{Ray{a}}){a}() = Ray{a,Float64}()

Base.angle{a}(d::Ray{a}) = a*π

# ensure the angle is always in (-1,1]
Ray(c,a,o) = Ray{a==0?false:(abs(a)==(1.0π)?true:mod(a/π-1,-2)+1),typeof(c)}(c,o)
Ray(c,a) = Ray(c,a,true)

Ray() = Ray{false}()

##deal with vector

function Ray(d::AbstractVector)
    @assert length(d)==2
    @assert abs(d[1])==Inf|| abs(d[2])==Inf

    if abs(d[2])==Inf
        Ray(d[1],angle(d[2]),true)
    else #d[1]==Inf
        Ray(d[2],angle(d[1]),false)
    end
end

isambiguous(d::Ray)=isnan(d.center)
Base.convert{a,T<:Number}(::Type{Ray{a,T}},::AnyDomain)=Ray{a,T}(NaN,true)
Base.convert{IT<:Ray}(::Type{IT},::AnyDomain)=Ray(NaN,NaN)

## Map interval

function mobiuspars(d::Ray)
    s=(d.orientation?1:-1)
    α=conj(cisangle(d))
    c=d.center
    s*α,-s*(1+α*c),α,1-α*c
end

for OP in (:mobius,:mobiusinv,:mobiusD,:mobiusinvD)
    @eval $OP(a::Ray,z) = $OP(mobiuspars(a)...,z)
end

ray_tocanonical(x)=(x==Inf)?1.:(x-1.)./(1+x)
ray_tocanonicalD(x)=(x==Inf)?0.:2*(1./(1+x)).^2
ray_fromcanonical(x)=(1+x)./(1-x)
ray_fromcanonicalD(x)=2*(1./(x-1.)).^2
ray_invfromcanonicalD(x)=(x-1.).^2/2

# atomatically vectorize over vector arg
@vectorize_1arg Number ray_tocanonical
@vectorize_1arg Number ray_tocanonicalD
@vectorize_1arg Number ray_fromcanonical
@vectorize_1arg Number ray_fromcanonicalD

for op in (:ray_tocanonical,:ray_tocanonicalD)
    @eval $op(o,x)=(o?1:-1)*$op(x)
end
ray_fromcanonical(o,x)=ray_fromcanonical((o?1:-1)*x)
ray_fromcanonicalD(o,x)=(o?1:-1)*ray_fromcanonicalD((o?1:-1)*x)
ray_invfromcanonicalD(o,x)=(o?1:-1)*ray_invfromcanonicalD((o?1:-1)*x)

cisangle{a}(::Ray{a})=cis(a*π)
cisangle(::Ray{false})=1
cisangle(::Ray{true})=-1

tocanonical(d::Ray,x) =
    ray_tocanonical(d.orientation,conj(cisangle(d)).*(x-d.center))
tocanonicalD(d::Ray,x) =
    conj(cisangle(d)).*ray_tocanonicalD(d.orientation,conj(cisangle(d)).*(x-d.center))
fromcanonical(d::Ray,v::AbstractArray) = eltype(d)[fromcanonical(d,vk) for vk in v]
fromcanonical(d::Ray,x) = cisangle(d)*ray_fromcanonical(d.orientation,x)+d.center
fromcanonicalD(d::Ray,x) = cisangle(d)*ray_fromcanonicalD(d.orientation,x)
invfromcanonicalD(d::Ray,x) = conj(cisangle(d))*ray_invfromcanonicalD(d.orientation,x)

arclength(d::Ray) = Inf

=={a}(d::Ray{a},m::Ray{a}) = d.center == m.center


text:
module Dashboard

using ..Model
using Dash
using DashHtmlComponents
using DashCoreComponents
using DashTable
using DashBootstrapComponents
using DelimitedFiles
using PlotlyJS
using ImageFiltering
using Agents: Agents
using JSON3
using DataFrames
using JLD2

state = Dict()

function readmapfile(filepath)
    content = read(filepath, String)
    if occursin(",", content)
        readdlm(filepath, ',')
    else
        readdlm(filepath)
    end
end

heatmapkwargs = (zauto=true, transpose=true, showscale=false)

function figure(p; kwargs...)
    figure = JSON3.read(json(p), NamedTuple{(:data, :layout, :frames)})
    return dcc_graph(; figure=figure, kwargs...)
end

function parameters_view()
    return dbc_card([#
        dbc_cardheader("Parameters"),
        dbc_cardbody([#
            dbc_formgroup([#
                dbc_label("Map:"),
                dbc_select(; id="crop", options=[], value=[]),
                html_a("Refresh"; id="refresh-map", href="javascript:void(0)"),
                html_br(),
                dbc_label("Init #BPH"),
                dbc_input(; type="number", id="nb_bph_init", value=200),
                dbc_label("Init position"),
                dbc_select(;
                    id="init_position",
                    options=[#
                        (value="corner", label="Corner"),
                        (value="random_c1", label="Random (2 patches scenario)"),
                        (value="random_c2", label="Random (9 patches scenario)"),
                        (value="border", label="Border"),
                    ],
                    value="corner",
                ),
                dbc_label("Kill probability"),
                dbc_input(;
                    type="number", id="pr_killed0", value=0.075, step=0.000001, max=1, min=0
                ),
                dbc_label("Replications"),
                dbc_input(; type="number", id="replication", value=50, max=5000, min=1),
            ]),
            dbc_cardlink([
                html_a("Run simulation"; href="javascript:void(0)", id="run-btn")
            ]),
        ]),
    ])
end

function map_view()
    return dbc_card([#
        dbc_cardheader("Map preview"),
        dcc_loading(dbc_cardbody([]; id="map-pv")),
    ])
end

function video_paramter()
    return dbc_card(
        [#
            dbc_cardheader("Video")
            dbc_cardbody(
                [
                    dbc_label("Number of steps: ")
                    dbc_input(; type="number", id="video-frames", value=2880)
                    dbc_label("Video seed: ")
                    dbc_input(; type="number", id="video-seed", value=1)
                    html_br()
                    dbc_cardlink([#
                        html_a("Generate"; id="video-btn", href="javascript:void(0)"),
                    ])
                ],
            )
        ],
    )
end

"""
    result_view()

View for exploring result (JLD2)
"""
function result_view()
    base_directory = joinpath(@__DIR__, "..")
    dir_names = readdir(base_directory)
    dir_paths = joinpath.(base_directory, dir_names)
    options = [
        (label=dir_name, value=dir_path) for
        (dir_name, dir_path) in zip(dir_names, dir_paths) if isdir(dir_path)
    ]
    return dbc_card(
        [
            dbc_cardheader("-")
            dbc_cardbody(
                [
                    dbc_label("Result folder")
                    dbc_select(; id="result-dir", options=options)
                    dbc_label("Select result")
                    dcc_loading(dbc_select(; id="result-to-view", options=[]))
                    html_p(
                        [
                            html_span("Static plot")
                            dbc_checkbox(; checked=true, id="result-static-plot")
                        ],
                    )
                    html_br()
                    dcc_loading(html_div(; id="result-plot"))
                ],
            )
        ],
    )
end

# @CALLBACKS

callbacks = Dict{Symbol,Function}()

callbacks[:result_view] = function (app, state)
    # Select folder
    callback!(
        _ -> nothing,
        app,
        Output("result-to-view", "value"),
        Input("result-dir", "value"),
    )

    # Select folder
    callback!(#
        app,
        Output("result-to-view", "options"),
        Input("result-dir", "value"),
    ) do result_path
        if isnothing(result_path)
            return []
        else
            names = readdir(result_path)
            paths = joinpath.(result_path, names)
            [
                (label=name, value=path) for (name, path) in zip(names, paths) if
                isfile(path) && endswith(name, ".jld2")
            ]
        end
    end

    # chose file + plot
    return callback!(
        app, #
        Output("result-plot", "children"),
        Input("result-to-view", "value"),
        State("result-static-plot", "checked"),
    ) do file, staticPlot
        if isnothing(file)
            return ""
        end

        data = jldopen(file) do f
            map(seed -> (seed => f[string(seed)]), 1:1000)
        end
        bph_traces = [
            scatter(;#
                x=df.step,
                y=df.count_bph,
                name="Seed $(seed)",
            ) for (seed, df) in data
        ]
        rice_traces = [
            scatter(;#
                x=df.step,
                y=df.food,
                name="Seed $(seed)",
            ) for (seed, df) in data
        ]
        bph_layout = Layout(; title="BPH")
        rice_layout = Layout(; title="Rice ≥ 50%")
        plt = hcat(plot(bph_traces, bph_layout), plot(rice_traces, rice_layout))

        # Plot layout
        metadata = jldopen(f -> f["metadata"], file)
        total_rice = data[1][2].food[1]
        passes = [seed for (seed, df) in data if df.food[end] < total_rice ÷ 2 - 100]
        npasses = length(passes)
        mapname = split(metadata.envmap, r"[\\/]")[end]
        str_passed = if npasses === 0
            "(No passed cases)"
        else
            aux = """
            Passed cases: $(npasses)/1000
            """
            aux = aux * "(" * join(Iterators.take(passes, 10), ", ") * ")"
        end
        title = """
        Map: $(mapname), #BPH: $(metadata.nb_bph_init),
        pos: $(metadata.init_position), pr_killed: $(metadata.pr_killed),
        $str_passed

        """
        relayout!(plt.plot, Layout(; showlegend=false, ymin=0, title=title))

        config = (staticPlot=staticPlot,)
        return figure(plt; config=config)
    end
end

callbacks[:refreshmap] = function (app, state)
    return callback!(
        app,#
        Output("crop", "options"),
        Output("crop", "value"),
        Input("refresh-map", "n_clicks"),
    ) do _
        map_options, map_value = let mapdir = joinpath(@__DIR__, "maps/")
            labels = readdir(mapdir)
            values = joinpath.(mapdir, labels)
            [(value=v, label=l) for (v, l) in zip(values, labels)], first(values)
        end
    end
end

callbacks[:drawmap] = function (app, state)
    callback!(
        app,
        Output("map-pv", "children"),
        Input("crop", "value"),
        Input("pr_killed0", "value"),
    ) do mappath, pr_killed0
        crop = readmapfile(mappath)
        trace = heatmap(;#
            z=crop,
            heatmapkwargs...,
            colorscale=[(i, "rgba(0.0, 1.0, 0.0, $i)") for i in 0:0.1:1],
        )
        trace2 = heatmap(;#
            z=isnan.(crop) * 1.0,
            colorscale=[(i, "rgba(1.0, 1.0, 0.0, 0.5)") for i in 0:0.1:1],
            heatmapkwargs...,
        )
        pr_killed = imfilter(isnan.(crop), Kernel.gaussian(3))
        pr_killed = pr_killed / maximum(pr_killed) * pr_killed0
        trace3 = heatmap(;
            z=pr_killed,
            colorscale=[(i, "rgba(1.0, 0.0, 0.0, $i)") for i in 0:0.1:1],
            heatmapkwargs...,
        )
        surface_flower = count(isnan, crop)
        surface_all = length(crop[:])
        surface_flower_pc = round(surface_flower / surface_all * 100; digits=2)
        layout = Layout(;
            title="Crop size: $(size(crop)). Flower/Total: $(surface_flower)/$(surface_all) ($surface_flower_pc%)",
        )
        fig = hcat(plot([trace2, trace]), plot([trace3]))
        relayout!(fig.plot, layout)
        figure(fig)
    end
end

callbacks[:run] = function (app, state)
    callback!(
        app,
        Output("simulation-output", "children"),
        Input("run-btn", "n_clicks_timestamp"),
        State("crop", "value"),
        State("nb_bph_init", "value"),
        State("init_position", "value"),
        State("pr_killed0", "value"),
        State("replication", "value"),
    ) do ts, map_path, nb_bph_init, init_position, pr_killed0, replication
        if isnothing(ts)
            return ""
        end
        #state["simulation-output-progress"] = 0
        adata = let bph(x) = true
            [(bph, count)]
        end
        mdata = let rice(model) = count(≥(0.5), model.food)
            [rice]
        end

        function run_simulation(seed)
            crop = readmapfile(map_path)
            init_position = Symbol(init_position)
            model = Model.init_model(
                crop, nb_bph_init, init_position, pr_killed0; seed=seed
            )
            adf, mdf = Agents.run!(
                model,
                Model.agent_step!,
                Model.model_step!,
                2880;
                adata=adata,
                mdata=mdata,
            )
            df = innerjoin(adf, mdf; on=:step)
            #state["simulation-output-progress"] = round(seed / replication * 100)
            return (seed, df)
        end
        # Run the model

        total_rice = count(!isnan, readmapfile(map_path))
        data = [run_simulation(sd) for sd in 1:replication]
        passes = [seed for (seed, df) in data if df.rice[end] < total_rice ÷ 2 - 100]
        npasses = length(passes)

        # Data traces
        bph_traces = [
            scatter(;#
                x=df.step,
                y=df.count_bph,
                name="Seed $(seed)",
            ) for (seed, df) in data
        ]
        rice_traces = [
            scatter(;#
                x=df.step,
                y=df.rice,
                name="Seed $(seed)",
            ) for (seed, df) in data
        ]
        bph_layout = Layout(; title="BPH")
        rice_layout = Layout(; title="Rice ≥ 50%")
        plt = hcat(plot(bph_traces, bph_layout), plot(rice_traces, rice_layout))

        # Plot layout
        mapname = split(map_path, r"[\\/]")[end]
        str_passed = if npasses === 0
            "(No passed cases)"
        else
            aux = """
            Passed cases: $(npasses)/$(replication)
            """
            aux = aux * "(" * join(Iterators.take(passes, 10), ", ") * ")"
        end
        title = """
        Map: $(mapname), #BPH: $(nb_bph_init),
        pos: $(init_position), pr_killed: $(pr_killed0),
        $str_passed

        """
        relayout!(plt.plot, Layout(; showlegend=false, ymin=0, title=title))

        # disable counter
        #state["simulation-output-progress"] = -1
        return figure(plt)
    end
end

#callbacks[:simulation_output_progress]
_ = function (app, state)
    return callback!(
        app,
        Output("simulation-output-progress", "value"),
        Output("simulation-output-progress", "children"),
        Output("simulation-output-progress", "style"),
        Output("simulation-output-watcher", "disabled"),
        Input("simulation-output-watcher", "n_intervals"),
        Input("run-btn", "n_clicks"),
    ) do nintervals, run_btn_clicks
        run_btn_clicks_prev = get(state, "run-btn-clicks", 0)
        run_btn_clicks = something(run_btn_clicks, 0)
        value = get(state, "simulation-output-progress", -1)
        value_str = "$(value)%"
        isloading = 100 > value ≥ 0 || run_btn_clicks > run_btn_clicks_prev
        state["run-btn-clicks"] = run_btn_clicks
        style = if isloading
            Dict()
        else
            Dict("display" => "none")
        end
        return (value, value_str, style, !isloading)
    end
end

callbacks[:run_video] = function (app, state)
    return callback!(
        app,
        Output("video-output", "children"),
        Input("video-btn", "n_clicks_timestamp"),
        State("video-frames", "value"),
        State("crop", "value"),
        State("nb_bph_init", "value"),
        State("init_position", "value"),
        State("pr_killed0", "value"),
        State("video-seed", "value"),
    ) do ts, frames, map_path, nb_bph_init, init_position, pr_killed0, seed
        if isnothing(ts)
            return ""
        end
        mapname = split(map_path, r"[/\\]")[end]
        vname = "BPH-$(mapname)-$(nb_bph_init)-$(init_position)-$(pr_killed0)-$(seed).mp4"
        vpath = normpath(joinpath(@__DIR__, "..", vname))

        crop = readmapfile(map_path)
        Model.video(
            vpath,
            crop,
            nb_bph_init,
            Symbol(init_position),
            pr_killed0;
            frames=frames,
            seed=seed,
        )

        if isfile(vpath)
            html_div("Video saved in:\n$(vpath)")
        else
            "Some thing is wrong with video: please check $vpath"
        end
    end
end

function simulation_output()
    return dbc_card(
        [
            dbc_cardheader("Output")
            dbc_cardbody([#
                #dcc_interval(; id="simulation-output-watcher"),
                #dbc_progress(;
                #    id="simulation-output-progress", value=25, striped=true, animated=true
                #),
                dcc_loading(html_div(""; id="simulation-output")),
            ],)
        ],
    )
end

function start(; host="127.0.0.1", port=8000, debug=true)
    @info "Debug = $debug"
    app = dash(; external_stylesheets=[dbc_themes.BOOTSTRAP])

    app.layout = dbc_container(
        [#
            html_h1(Model.name)
            html_h4("Input")
            dbc_row([# row
                dbc_col(parameters_view(); width=2)
                dbc_col(map_view();)
            ])
            html_br()
            html_h4("Plot output")
            simulation_output()
            html_h4("Video output")
            dbc_row(
                [
                    dbc_col(video_paramter(); width=3)
                    dbc_col(
                        dbc_card(
                            [
                                dbc_cardheader("Video")
                                dcc_loading(dbc_cardbody(; id="video-output"))
                            ],
                        );
                        width=9,
                    )
                ],
            )
            html_br()
            html_h4("Result view")
            result_view()
        ],
    )

    for (k, f!) in callbacks
        f!(app, state)
    end

    return run_server(app, host, port; debug=debug)
end

end

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
export
    dnnAlgorithmPoolingMax,
    dnnAlgorithmPoolingMin,
    dnnAlgorithmPoolingAvg

function pooling{T,N}(x::Array{T}, op, ksize::Tuple, padding::NTuple{N,Int}, strides::Tuple;
    border=dnnBorderZeros)
    attr = dnnPrimitiveAttributesCreate()
    layout = dnnLayoutCreate(x)
    cksize = Csize_t[ksize...]
    cstrides = Csize_t[strides...]
    coffsets = Cint[padding...]

    p = Ptr{Cvoid}[0]
    dnnPoolingCreateForward_F32(p, attr, op, layout, cksize, cstrides, coffsets, border)
    p_pool = p[1]

    p = Ptr{Cvoid}[0]
    dnnLayoutCreateFromPrimitive_F32(p, p_pool, dnnResourceWorkspace)
    p_workspace = p[1]
    p = Ptr{Cvoid}[0]
    dnnAllocateBuffer_F32(p, p_workspace)
    p_buf = p[1]

    dims = ntuple(i -> (size(x,i)+2padding[i]-ksize[i]) ÷ strides[i] + 1, N)
    y = Array{T}(dims..., size(x,3), size(x,4))
    resources = Ptr{Cvoid}[pointer(x),pointer(y),0,0,0,0,0,p_buf]
    dnnExecute(p_pool, resources)

    dnnPrimitiveAttributesDestroy(attr)
    dnnLayoutDelete_F32(layout)
    dnnDelete(p_pool)
    y
end


text:
"""
    PeriodicPolyacetylene(N; ϕ=π/4, l=DEFAULT_LENGTH, handedness=1, δz=nothing, θ=nothing, r0=nothing)

Molecule constructor for (twisted) Polyacetylene with periodic boundary conditions and twist angle, `ϕ`.

"""
@export struct PeriodicPolyacetylene <: AbstractMoleculeConstructor
    θ
    l
    handedness
    δz
    ϕ
    r0

    function PeriodicPolyacetylene(
        ;
        θ = π/2,
        l = DEFAULT_LENGTH,
        handedness = 1,
        δz = nothing,
        ϕ = nothing,
        r0 = nothing
    )
        return new(θ, l, handedness, δz, ϕ, r0)
    end
end

function (m::PeriodicPolyacetylene)(; kwargs...)
    θ = get(kwargs, :θ, m.θ)
    l = get(kwargs, :l, m.l)
    handedness = get(kwargs, :handedness, m.handedness)
    δz = get(kwargs, :δz, m.δz)
    ϕ = get(kwargs, :ϕ, m.ϕ)
    r0 = get(kwargs, :r0, m.r0)

    ρ = atan(2tan(θ/2))

    if δz !== nothing && θ !== nothing
        r0 = [l, 0, 0]
    end

    if δz === nothing
        δz = √3*l*sin(ρ)/2
    end
    if ϕ === nothing
        ϕ = π - 2atan(√3*cos(ρ))*handedness
    end

    if r0 === nothing
        r0 = [√(l^2 - δz^2) / 2sin(ϕ/2), 0, 0]
    end
    r0::Vector

    mols = [Molecule(), Molecule()]

    for i in 1:2
        rotation_matrix = LinearAlgebra.Givens(1, 2, cos((i-1)*ϕ), -sin((i-1)*ϕ))
        shift_vector = [0, 0, (i-1)*δz]
        axes = makeaxes(π/2 + (i-1)*ϕ, ρ*handedness)
        r1 = shift_vector + rotation_matrix*r0
        r2 = shift_vector + rotation_matrix*(r0 + [abs(l), 0, 0])
        A = Carbon(vec(r1...), orbital"2s", orbital"2p_x", orbital"2p_y", orbital"2p_z", axes=axes)
        B = Hydrogen(vec(r2...), orbital"1s", axes=axes)
        push!(mols[i], A, B)
        push!(mols[i], Bond(A, B))
    end

    return Molecules.UnitCell(mols..., Set([Bond(mols[1].atoms[1], mols[2].atoms[1])])), nothing, nothing
end


text:
# plot soundings from Chile

fp = "/media/wiebke/Samsung_T3/WiebkesDaten_fromSamsung_analysiert_Bolivia_CLOUD/Bolivia-atmo/DatafromOtherGroups/Soundings/"
antofagastafile = "$(fp)CIM00085442-data.txt"
islapasqua = open("$(fp)CIM00085469-data.txt")

i = 0
headerlines = []
open(antofagastafile) do file
    for ln in eachline(file)
	global i = i+1
        if length(ln) > 70 && occursin("#CIM00085442 2018 05", ln)
		append!(headerlines, i)
	end
    end
end
antofagasta = readlines(antofagastafile)
antofagasta[headerlines[1]+1:headerlines[2]-1]

DateTime(antofagasta[headerlines[1]][14:31], "yyyy mm dd HH HHMM")

press = Float64[]
GPH = Float64[]
temp = Float64[]
rh = Float64[]
dpdp = Float64[]
wdir = Float64[]
wspeed = Float64[]

for j in headerlines[1]+1:headerlines[2]-1
	append!(press, parse(Float64, antofagasta[j][10:15]))
	append!(GPH, parse(Float64, antofagasta[j][17:21]))
	append!(temp, parse(Float64, antofagasta[j][23:27])/10)
	append!(rh, parse(Float64, antofagasta[j][29:33])/10)
	append!(dpdp, parse(Float64, antofagasta[j][35:39])/10)
	append!(wdir, parse(Float64, antofagasta[j][41:45]))
	append!(wspeed, parse(Float64, antofagasta[j][47:51])/10)
end

figure()
plot(temp, press)


text:
using CoordinateTransformations, Rotations, StaticArrays
using RoboLib.Geom: Pose2D
export AckerParams, AckerData

struct AckerParams{SCALAR<:AbstractFloat} <: MotionParams
    car_length::SCALAR
end

mutable struct AckerData{T} <: MotionData
    pose::Pose2D{T}
    ctrl::SVector{2, T}
    # TODO(cxs): special case
    @inline function AckerData{T}(pose, ctrl) where T
        pose = Pose2D{T}(pose)
        ctrl = SVector{2, T}(ctrl)
        new{T}(pose, ctrl)
    end
    @inline function AckerData{T}() where T
        pose = Pose2D{T}()
        ctrl = zeros(SVector{2, T})
        new{T}(pose, ctrl)
    end
end
@inline AckerData(pose, ctrl) = AckerData{Float64}(pose, ctrl)
@inline AckerData() = AckerData{Float64}()

@inline function _step(x, y, theta, dt, v, delta, car_length)
    beta = atan(tan(delta) / 2)
    theta_dot = (2 * v * sin(beta)) / car_length
    theta_n = theta + theta_dot * dt

    x_dot = (car_length / (2 * sin(beta))) * (sin(theta_n + beta) - sin(theta + beta))
    y_dot = (car_length / (2 * sin(beta))) * (-cos(theta_n + beta) + cos(theta + beta))

    x_n = x + x_dot
    y_n = y + y_dot

    x_n, y_n, theta_n
end

@inline function _step_nodelta(x, y, theta, dt, v, delta, car_length)
    beta = atan(tan(delta) / 2)
    s_theta_beta, c_theta_beta = sincos(theta + beta)

    x_dot = v * c_theta_beta
    y_dot = v * s_theta_beta
    theta_dot = (2 * v * sin(beta)) / car_length

    x_n = x + x_dot * dt
    y_n = y + y_dot * dt
    theta_n = theta + theta_dot * dt

    x_n, y_n, theta_n
end

function step!(data::AckerData{T}, model::AckerParams, dt) where T
    v, delta = data.ctrl
    xc, yc, thetac = data.pose.statev
    if abs(delta) < 0.001
        x, y, theta = _step_nodelta(xc, yc, thetac, dt, v, delta, model.car_length)
    else
        x, y, theta = _step(xc, yc, thetac, dt, v, delta, model.car_length)
    end

    data.pose = Pose2D{T}(x, y, theta)
    data
end


text:
const OP_DUP = 0x76
const OP_HASH160 = 0xa9
const OP_EQUALVERIFY = 0x88
const OP_CHECKSIG = 0xac


text:
function loadConfig()
  cfg = Dict{Symbol,Any}()
  data =   YAML.load(open(joinpath(dirname(@__FILE__),"cam_cal.yml")))
  bRc = eval(Meta.parse("["*data["extrinsics"]["bRc"][1]*"]"))
  # convert to faster symbol lookup
  cfg[:extrinsics] = Dict{Symbol,Any}()
  cfg[:extrinsics][:bRc] = bRc
  cfg[:intrinsics] = Dict{Symbol,Any}()
  cfg[:intrinsics][:height] = data["left"]["intrinsics"]["height"]
  cfg[:intrinsics][:width] = data["left"]["intrinsics"]["width"]
  haskey(data["left"]["intrinsics"], "camera_matrix") ? (cfg[:intrinsics][:cam_matrix] = data["left"]["intrinsics"]["camera_matrix"]) : nothing
  cfg[:intrinsics][:cx] = data["left"]["intrinsics"]["cx"]
  cfg[:intrinsics][:cy] = data["left"]["intrinsics"]["cy"]
  cfg[:intrinsics][:fx] = data["left"]["intrinsics"]["fx"]
  cfg[:intrinsics][:fy] = data["left"]["intrinsics"]["fy"]
  cfg[:intrinsics][:k1] = data["left"]["intrinsics"]["k1"]
  cfg[:intrinsics][:k2] = data["left"]["intrinsics"]["k2"]
  cfg
end

cfg = loadConfig()

cw = cfg[:intrinsics][:cx] + parsed_args["cxoffset"]
ch = cfg[:intrinsics][:cy] + parsed_args["cyoffset"]
fx = fy = parsed_args["focalscale"]*cfg[:intrinsics][:fx]
camK = [fx 0 cw; 0 fy ch; 0 0 1]
tagsize = parsed_args["tagscale"]*0.172
# k1,k2 = cfg[:intrinsics][:k1], cfg[:intrinsics][:k2] # makes it worse
k1, k2 = 0.0, 0.0

# tag extrinsic rotation
Rx = RotX(-pi/2)
Rz = RotZ(-pi/2)
bTc= LinearMap(Rz) ∘ LinearMap(Rx)

datadir = joinpath(ENV["HOME"],"data","racecar")

imgfolder = "images"
datafolder = ENV["HOME"]*"/data/racecar/$(parsed_args["folder_name"])/";

# TODO: automate construction of camidxs from data folder itself
global camidxs = 0:5:5
if parsed_args["folder_name"] == "labrun1"
  camidxs =  100:5:1060
elseif parsed_args["folder_name"] == "labrun2"
  camidxs =  0:5:1625
elseif parsed_args["folder_name"] == "labrun3"
  camidxs =  205:5:1165
elseif parsed_args["folder_name"] == "labrun4"
  camidxs =  0:5:1715
elseif parsed_args["folder_name"] == "labrun5"
  camidxs =  0:5:1020
elseif parsed_args["folder_name"] == "labrun6"
  camidxs =  0:5:1795
elseif parsed_args["folder_name"] == "labrun7"
  camidxs =  0:5:1750 # 0:5:2135
elseif parsed_args["folder_name"] == "labrun8"
    camidxs =  0:5:1765
elseif parsed_args["folder_name"] == "straightrun3"
  camidxs =  175:5:370
else
  @error "unknown Racecar data folder, $(parsed_args["folder_name"])"
end

# datafolder = ENV["HOME"]*"/data/racecar/straightrun3/"  # 175:5:370
# datafolder = ENV["HOME"]*"/data/racecar/labrun2/; camidxs =  0:5:1625
# datafolder = ENV["HOME"]*"/data/racecar/labrun3/"; # camidxs =
# datafolder = ENV["HOME"]*"/data/racecar/$(parsed_args["folder_name"])/"; camidxs =  0:5:1020
# datafolder = ENV["HOME"]*"/data/racecar/labrun6/"; camidxs =  0:5:1795
# datafolder = ENV["HOME"]*"/data/racecar/labrun7/"; camidxs =  0:5:2135
# datafolder = ENV["HOME"]*"/data/racecar/labrun8/"; camidxs =  0:5:


text:
#!/usr/bin/env julia
################################################################################
#
#  Copyright (c) 2015 Wojciech Migda
#  All rights reserved
#  Distributed under the terms of the MIT license
#
################################################################################
#
#  Filename: data_to_h5.jl
#
#  Decription:
#       Convert BMP data and put it along feature labels into H5 storages.
#
#  Authors:
#       Wojciech Migda
#
################################################################################
#
#  History:
#  --------
#  Date         Who  Ticket     Description
#  ----------   ---  ---------  ------------------------------------------------
#  2015-11-22   wm              Initial version
#  2015-12-09   wm              Redesign. +CLI, +image scaling, +parallel
#
################################################################################

addprocs(div(CPU_CORES, 2))

@everywhere using Images
@everywhere using ColorTypes
@everywhere using HDF5
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        #"--opt1"
        #    help = "an option with an argument"
        "--data-dir", "-D"
            help = "directory with input CSV files, BMP 'train' and 'test' subfolders, and where H5 will be stored"
            arg_type = String
            required = true
        "--do-x-labels"
            help = "store X features file labels"
            arg_type = Bool
            default = false
        "--do-y"
            help = "store 'y' train vector"
            arg_type = Bool
            default = false
        "--image-size", "-s"
            help = "size to which BMP files will be resized (in pixels). 0 or negative value to skip generating feature vectors"
            arg_type = Int
            default = 0

        #"arg1"
        #    help = "a positional argument"
        #    required = true
    end

    return parse_args(s)
end

@everywhere function image_producer(DATA_PATH, FILE_NAMES)

    function core()
        for img_name in FILE_NAMES
            produce(imread("$(DATA_PATH)/$(img_name).Bmp"))
        end
    end

    Task(() -> core())

end

@everywhere function image_resizer(IMAGES, SHAPE)

    function core()
        for img in IMAGES
            produce(Images.imresize(img, SHAPE))
        end
    end

    Task(() -> core())

end

@everywhere function image_converter(IMAGES, VLEN)

    function core()
        for img in IMAGES
            img |>
                img -> convert(Array{Gray}, float32(img)) |>
                img -> convert(Array{Float32}, img) |>
                img -> reshape(img, 1, VLEN) |>
                produce
        end
    end

    Task(() -> core())

end

function transform(
    DATA_PATH::String,
    do_x_labels::Bool,
    do_y::Bool,
    IMSZ::Int
    )

    const TRAIN_FEAT_LABELS = readdlm(DATA_PATH * "/trainLabels.csv", ',', ASCIIString, header=true)[1] # [2] is the header
    const TEST_FEAT_LABELS = readdlm(DATA_PATH * "/sampleSubmission.csv", ',', ASCIIString, header=true)[1]
    const train_y = map(x -> int(x[1]), TRAIN_FEAT_LABELS[:, 2])

    if do_y
        h5write("$(DATA_PATH)/trainY.h5", "gsv/train/Y/data", train_y)
    end

    if do_x_labels
        h5write("$(DATA_PATH)/testXlabels.h5", "gsv/test/X/labels", TEST_FEAT_LABELS[:, 1])
    end

    if IMSZ > 0
        const work_descr = [
            {
                "subfolder" => "train",
                "files_labels" => TRAIN_FEAT_LABELS[:, 1],
                "h5prefix" => "train",
                "dataset" => "gsv/train/X/data"
            }
            {
                "subfolder" => "test",
                "files_labels" => TEST_FEAT_LABELS[:, 1],
                "h5prefix" => "test",
                "dataset" => "gsv/test/X/data"}
        ]

        @sync @parallel for work = work_descr
            const X = image_producer(DATA_PATH * "/" * work["subfolder"], work["files_labels"]) |>
                i -> image_resizer(i, (IMSZ, IMSZ)) |>
                i -> image_converter(i, IMSZ ^ 2) |>
                i -> vcat(i...)

            h5write(DATA_PATH * "/" * work["h5prefix"] * "X_$(IMSZ).h5", work["dataset"], X)
        end
    end

end

function main()
    parsed_args = parse_commandline()
    println("Parsed args:")
    for (arg,val) in parsed_args
        println("  $arg  =>  $val")
    end

    transform(
        parsed_args["data-dir"],
        parsed_args["do-x-labels"],
        parsed_args["do-y"],
        parsed_args["image-size"]
    )
end

if ~isinteractive()
    main()
end


text:
@testset "Points on the sphere" begin
    lmax = 10

    N = lmax + 1
    Θ, Φ = sph_points(N)
    @test length(Θ) == N
    M = length(Φ)
    @test all(0 < θ < π for θ in Θ)
    @test all(0 ≤ ϕ < 2π for ϕ in Φ)
end


text:
"""
Methods for position space Chebyshev methods.
"""
module Chebyshev

export chep_pts, mat_X, mat_D1, mat_D2, mat_fd_D1, mat_fd_D2, to_cheb, to_real

include("CustomTypes.jl")

using .CustomTypes
using SparseArrays

"""
    cheb_pts(nx::myI)::Vector{myF}

Computes Generates nx Chebyshev points in [-1,1]
"""
function cheb_pts(nx::myI)::Vector{myF}
   return [cos(pi*i/(nx-tomyF(1))) for i=0:(nx-1)]
end

"""
    cheb_pts(xmin::myF, xmax::myF, nx::myI)::Vector{myF}

Computes Chebyshev points on interval [xmin,xmax]
"""
function cheb_pts(xmin::myF, xmax::myF, nx::myI)::Vector{myF}
   pts = cheb_pts(nx)
   m = (xmax - xmin)/tomyF(2)
   b = (xmax + xmin)/tomyF(2)
   return [m*pts[i] + b for i=1:nx]
end

"""
    mat_X(
      xmin::myF,
      xmax::myF,
      nx::myI
    )::SparseMatrixCSC{myF, myI}

Computes matrix for multiplication of x in real space.
"""
function mat_X(
      xmin::myF,
      xmax::myF,
      nx::myI
   )::SparseMatrixCSC{myF, myI}
   @assert nx>4

   X = Vector{myI}(undef,0)
   Y = Vector{myI}(undef,0)
   V = Vector{myF}(undef,0)

   pts = cheb_pts(xmin, xmax, nx)

   for i=1:nx
      push!(X,i)
      push!(Y,i)
      push!(V,pts[i])
   end

   return sparse(X,Y,V)
end

"""
    mat_D1(
      xmin::myF,
      xmax::myF,
      nx::myI
    )::Matrix{myF}

Computes derivative matrix D1 in real space.
"""
function mat_D1(
      xmin::myF,
      xmax::myF,
      nx::myI
   )::Matrix{myF}
   @assert nx>4

   M = Matrix{myF}(undef,nx,nx)
   n = nx-1
   pts = cheb_pts(nx)

   M[1,1]   = (tomyF(2)*(n^2) + tomyF(1)) / tomyF(6)
   M[nx,nx] = -M[1,1]

   M[1,nx] = tomyF(0.5)*((-1.0)^n)
   M[nx,1] = -M[1,nx]

   for i=2:(nx-1)
      M[1 ,i] = tomyF(2.0) *((-1)^(i+1 )/(tomyF(1) - pts[i]))
      M[nx,i] = tomyF(-2.0)*((-1)^(i+nx)/(tomyF(1) + pts[i]))
      M[i, 1] = tomyF(-0.5)*((-1)^(i+1 )/(tomyF(1) - pts[i]))
      M[i,nx] = tomyF(0.5) *((-1)^(i+nx)/(tomyF(1) + pts[i]))

      for j=2:(nx-1)
         if i!=j
            M[i,j] = ((tomyF(-1))^(i+j))/(pts[i] - pts[j])
         end
      end
   end

   for i=1:(nx-1)
      M[i,i] = 0.0
      for j=1:nx
         if i!=j
            M[i,i] -= M[i,j]
         end
      end
   end

   M .*= tomyF(2)/(xmax-xmin)

   return M
end

"""
    mat_D2(
      xmin::myF,
      xmax::myF,
      nx::myI
    )::Matrix{myF}

Computes derivative matrix D1 in real space.
"""
function mat_D2(
      xmin::myF,
      xmax::myF,
      nx::myI
   )::Matrix{myF}
   @assert nx>4

   M = Matrix{myF}(undef,nx,nx)
   n = nx-1
   pts = cheb_pts(nx)

   M[1,1]   = (tomyF(2)*(n^2) + tomyF(1)) / tomyF(6)
   M[nx,nx] = -M[1,1]

   M[1,nx] = tomyF(0.5)*((-1.0)^n)
   M[nx,1] = -M[1,nx]

   for i=2:(nx-1)
      M[1 ,i] = tomyF(2.0) *((-1)^(i+1 )/(tomyF(1) - pts[i]))
      M[nx,i] = tomyF(-2.0)*((-1)^(i+nx)/(tomyF(1) + pts[i]))
      M[i, 1] = tomyF(-0.5)*((-1)^(i+1 )/(tomyF(1) - pts[i]))
      M[i,nx] = tomyF(0.5) *((-1)^(i+nx)/(tomyF(1) + pts[i]))

      for j=2:(nx-1)
         if i!=j
            M[i,j] = ((tomyF(-1))^(i+j))/(pts[i] - pts[j])
         end
      end
   end

   for i=1:(nx-1)
      M[i,i] = 0.0
      for j=1:nx
         if i!=j
            M[i,i] -= M[i,j]
         end
      end
   end

   M .*= tomyF(2)/(xmax-xmin)

   return M
end

"""
    mat_fd_D1(
      xmin::myF,
      xmax::myF,
      nx::myI
    )::Matrix{myF}

Compute 2nd order finite difference in Chebyshev points.
"""
function mat_fd_D1(
      xmin::myF,
      xmax::myF,
      nx::myI
   )::Matrix{myF}
   @assert nx>4

   P = zeros(myF,nx,nx)
   pts = cheb_pts(xmin,xmax,nx)

   h(i) = pts[i+1] - pts[i]

   X = Vector{myF}()
   Y = Vector{myF}()
   V = Vector{myF}()

   append!(X,1)
   append!(Y,1)
   append!(V,
      (2*pts[1]-pts[2]-pts[3])
      /
      ((pts[1]-pts[2])*(pts[1]-pts[3]))
     )
   append!(X,1)
   append!(Y,2)
   append!(V,
      (1.0/(pts[2]-pts[1])) + (1.0/(pts[3]-pts[2]))
     )
   append!(X,1)
   append!(Y,3)
   append!(V,
      (pts[1]-pts[2])
      /
      ((pts[3]-pts[1])*(pts[3]-pts[2]))
     )

   for i=2:(nx-1)
      append!(X,i)
      append!(Y,i-1)
      append!(V,-(h(i)/h(i-1))/(h(i-1)+h(i)))
      append!(X,i)
      append!(Y,i)
      append!(V,(1.0/h(i-1)) - (1.0/h(i)))
      append!(X,i)
      append!(Y,i+1)
      append!(V,(h(i-1)/h(i))/(h(i-1)+h(i)))
   end

   append!(X,nx)
   append!(Y,nx)
   append!(V,
      (1.0/(pts[nx]-pts[nx-2])) + (1.0/(pts[nx]-pts[nx-1]))
     )
   append!(X,nx)
   append!(Y,nx-1)
   append!(V,
      (pts[nx-2]-pts[nx])
      /
      ((pts[nx-2]-pts[nx-1])*(pts[nx-1]-pts[nx]))
     )
   append!(X,nx)
   append!(Y,nx-2)
   append!(V,
      (pts[nx]-pts[nx-1])
      /
      ((pts[nx-2]-pts[nx-1])*(pts[nx-2]-pts[nx]))
     )

   return sparse(X,Y,V)
end

"""
    mat_fd_D2(
      xmin::myF,
      xmax::myF,
      nx::myI
    )::Matrix{myF}

Compute 2nd order finite difference in Chebyshev points.
"""
function mat_fd_D2(
      xmin::myF,
      xmax::myF,
      nx::myI
   )::Matrix{myF}
   @assert nx>4

   P = zeros(myF,nx,nx)
   pts = cheb_pts(xmin,xmax,nx)

   h(i) = pts[i+1] - pts[i]

   X = Vector{myF}()
   Y = Vector{myF}()
   V = Vector{myF}()

   append!(X,1)
   append!(Y,1)
   append!(V,
      2*(3*pts[1]-pts[2]-pts[3]-pts[4])
      /
      ((pts[1]-pts[2])*(pts[1]-pts[3])*(pts[1]-pts[4]))
     )
   append!(X,1)
   append!(Y,2)
   append!(V,
      2*(-2*pts[1]+pts[3]+pts[4])
      /
      ((pts[1]-pts[2])*(pts[2]-pts[3])*(pts[2]-pts[4]))
     )
   append!(X,1)
   append!(Y,3)
   append!(V,
      2*(-2*pts[1]+pts[2]+pts[4])
      /
      ((pts[1]-pts[3])*(-pts[2]+pts[3])*(pts[3]-pts[4]))
     )
   append!(X,1)
   append!(Y,4)
   append!(V,
      2*(-2*pts[1]+pts[2]+pts[3])
      /
      ((pts[2]-pts[4])*(-pts[1]+pts[4])*(-pts[3]+pts[4]))
     )

   for i=2:(nx-1)
      append!(X,i)
      append!(Y,i-1)
      append!(V,(2.0/h(i-1))/(h(i-1)+h(i)))
      append!(X,i)
      append!(Y,i)
      append!(V,- 2.0/(h(i-1)*h(i)))
      append!(X,i)
      append!(Y,i+1)
      append!(V,(2.0/h(i))/(h(i-1)+h(i)))
   end
   append!(X,nx)
   append!(Y,nx)
   append!(V,
      2*(3*pts[nx]-pts[nx-1]-pts[nx-2]-pts[nx-3])
      /
      ((pts[nx]-pts[nx-1])*(pts[nx]-pts[nx-2])*(pts[nx]-pts[nx-3]))
     )
   append!(X,nx)
   append!(Y,nx-1)
   append!(V,
      2*(-2*pts[nx]+pts[nx-2]+pts[nx-3])
      /
      ((pts[nx]-pts[nx-1])*(pts[nx-1]-pts[nx-2])*(pts[nx-1]-pts[nx-3]))
     )
   append!(X,nx)
   append!(Y,nx-2)
   append!(V,
      2*(-2*pts[nx]+pts[nx-1]+pts[nx-3])
      /
      ((pts[nx]-pts[nx-2])*(-pts[nx-1]+pts[nx-2])*(pts[nx-2]-pts[nx-3]))
     )
   append!(X,nx)
   append!(Y,nx-3)
   append!(V,
      2*(-2*pts[nx]+pts[nx-1]+pts[nx-2])
      /
      ((pts[nx-1]-pts[nx-3])*(-pts[nx]+pts[nx-3])*(-pts[nx-2]+pts[nx-3]))
     )

   return sparse(X,Y,V)
end

"""
    to_cheb(f::Vector{T})::Vector{T} where T<:Number

Convert to Chebyshev space.
We assume we are working with Chebyshev-Gauss-Lobatto points.
"""
function to_cheb(f::Vector{T})::Vector{T} where T<:Number

   N = length(f) - 1
   c = zeros(T,N+1)

   for i=1:(N+1)
      n = i-1

      c[i] += f[1]/tomyF(N)
      c[i] += ((-1)^n)*f[end]/tomyF(N)

      for k=2:N
         c[i] += (2.0/N)*f[k]*cos(n*(k-1)*pi/tomyF(N))
      end
   end
   c[1] /= tomyF(2) ## from normalization of inner product

   return c
end

"""
    to_real(c::Vector{T})::Vector{T} where T<:Number

Convert to Real space.
We assume we are working with Chebyshev-Gauss-Lobatto points.
"""
function to_real(c::Vector{T})::Vector{T} where T<:Number

   N = length(c) - 1
   f = zeros(T,N+1)

   for i=1:(N+1)
      n = i-1
      for j=1:(N+1)
         f[j] += c[i]*(cos(n*(j-1)*pi/tomyF(N)))
      end
   end
   return f
end

end


text:
(function (checkfermat a b c n)
  (if (&& (== (+ (^ a n) (^ b n)) (^ c n)) (> n 2))
      (println "Holy smokes, Fermat was wrong!")
      (println "No, that doesn't work.")))

(checkfermat 1 2 3 4)
(checkfermat 1 2 3 2)


text:
using Revise # reduce recompilation time
using Plots
# using Documenter
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using UnPack

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D

using SetupDG

push!(LOAD_PATH, "./examples/EntropyStableEuler.jl/src")
using EntropyStableEuler
using EntropyStableEuler.Fluxes1D

function wavespeed_1D(rho,rhou,E)
    p = pfun_nd(rho,(rhou,),E)
    cvel = @. sqrt(γ*p/rho)
    return @. abs(rhou/rho) + cvel
end
unorm(U) = sum(map((x->x.^2),U))
function pfun_nd(rho, rhoU, E)
    rhoUnorm2 = unorm(rhoU)./rho
    return @. (γ-1)*(E - .5*rhoUnorm2)
end

function primitive_to_conservative_hardcode(rho,U,p)
    rhoU = rho.*U
    Unorm = unorm(U)
    E = @. p/(γ-1) + .5*rho*Unorm
    return (rho,rhoU,E)
end

function euler_fluxes_1D(rhoL,uL,betaL,rhologL,betalogL,
    rhoR,uR,betaR,rhologR,betalogR)

    rholog = logmean.(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean.(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = (@. .5*(rhoL+rhoR))
    uavg   = (@. .5*(uL+uR))

    unorm = (@. uL*uR)
    pa    = (@. rhoavg/(betaL+betaR))
    f4aux = (@. rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm)

    FxS1 = (@. rholog*uavg)
    FxS2 = (@. FxS1*uavg + pa)
    FxS3 = (@. f4aux*uavg)

    return (FxS1,FxS2,FxS3)
end

function euler_fluxes(rhoL,uL,betaL,rhoR,uR,betaR)
    rhologL,betalogL,rhologR,betalogR = map(x->log.(x),(rhoL,betaL,rhoR,betaR))
    return euler_fluxes_1D(rhoL,uL,betaL,rhologL,betalogL,
                           rhoR,uR,betaR,rhologR,betalogR)
end

function dopri45_coeffs()
    rk4a = [0.0             0.0             0.0             0.0             0.0             0.0         0.0
            0.2             0.0             0.0             0.0             0.0             0.0         0.0
            3.0/40.0        9.0/40.0        0.0             0.0             0.0             0.0         0.0
            44.0/45.0      -56.0/15.0       32.0/9.0        0.0             0.0             0.0         0.0
            19372.0/6561.0 -25360.0/2187.0  64448.0/6561.0  -212.0/729.0    0.0             0.0         0.0
            9017.0/3168.0  -355.0/33.0      46732.0/5247.0  49.0/176.0      -5103.0/18656.0 0.0         0.0
            35.0/384.0      0.0             500.0/1113.0    125.0/192.0     -2187.0/6784.0  11.0/84.0   0.0 ]

    rk4c = vec([0.0 0.2 0.3 0.8 8.0/9.0 1.0 1.0 ])

    # coefficients to evolve error estimator = b1-b2
    rk4E = vec([71.0/57600.0  0.0 -71.0/16695.0 71.0/1920.0 -17253.0/339200.0 22.0/525.0 -1.0/40.0 ])

    return rk4a,rk4E,rk4c
end

const TOL = 1e-16
"Approximation parameters"
N = 4 # The order of approximation
K = 2^7
T = 0.1

# Becker viscous shocktube
const γ = 1.4
const M_0 = 3.0
const mu = 0.1
const lambda = 2/3*mu
const Pr = 3/4
const cp = γ/(γ-1)
const cv = 1/(γ-1)
const kappa = mu*cp/Pr

const v_inf = 0.2
const rho_0 = 1.0
const v_0 = 1.0
const m_0 = rho_0*v_0
const v_1 = (γ-1+2/M_0^2)/(γ+1)
const v_01 = sqrt(v_0*v_1)

# const uL = v_0+v_inf
# const uR = v_1+v_inf
# const rhoL = m_0/v_0
# const rhoR = m_0/v_1
# const eL = 1/(2*γ)*((γ+1)/(γ-1)*v_01^2-v_0^2)
# const eR = 1/(2*γ)*((γ+1)/(γ-1)*v_01^2-v_1^2)
# const pL = (γ-1)*rhoL*eL
# const pR = (γ-1)*rhoR*eR
# const EL = pL/(γ-1)+0.5*rhoL*uL^2
# const ER = pR/(γ-1)+0.5*rhoR*uR^2

const XL = -2.0#-1.0
const XR = 2.0
gr(size=(300,300),ylims=(0,5.0),legend=false,markerstrokewidth=1,markersize=2)
plot()

"Mesh related variables"
VX = LinRange(XL,XR,K+1)
EToV = transpose(reshape(sort([1:K; 2:K+1]),2,K))
rd = init_reference_interval(N)

# High order mesh
@unpack r,M,Dr,Vq,Pq,Vf,LIFT,wq,VDM = rd
V1 = vandermonde_1D(1,r)/vandermonde_1D(1,[-1;1])
x = V1*VX[transpose(EToV)]
xf = Vf*x
mapM = reshape(1:2*K,2,K)
mapP = copy(mapM)
mapP[1,2:end] .= mapM[2,1:end-1]
mapP[2,1:end-1] .= mapM[1,2:end]

"""Geometric factors"""
J = repeat(transpose(diff(VX)/2),N+1,1)
nrJ = [-1;1]
wf  = [1;1]
nxJ = repeat([-1;1],1,K)
rxJ = 1.0

"""Geometric factors"""
J = repeat(transpose(diff(VX)/2),N+1,1)
nxJ = repeat([-1;1],1,K)
rxJ = 1.0

# construct hybridized SBP operators
Qr = Pq'*M*Dr*Pq
Ef = Vf*Pq
Br = diagm(wf.*nrJ)
Qrh = .5*[Qr-Qr' Ef'*Br;
         -Br*Ef  Br]

Vh = [Vq;Vf]
Ph = M\transpose(Vh)
VhP = Vh*Pq

# make sparse skew symmetric versions of the operators"
# precompute union of sparse ids for Qr, Qs
Qrh_skew = .5*(Qrh-transpose(Qrh))

"""Initial condition"""
function bisection_solve_velocity(x,max_iter,tol)
    v_L = v_1
    v_R = v_0
    num_iter = 0

    L_k = kappa/m_0/cv
    f(v) = -x + 2*L_k/(γ+1)*log((v_0-v)^((v_0/(v_0-v_1)))*(v-v_1)^(-(v_1)/(v_0-v_1)))

    v_new = (v_L+v_R)/2
    while num_iter < max_iter
        v_new = (v_L+v_R)/2

        if abs(f(v_new)) < tol
            return v_new
        elseif sign(f(v_L)) == sign(f(v_new))
            v_L = v_new
        else
            v_R = v_new
        end
        num_iter += 1
    end

    return v_new
end

const max_iter = 10000
const tol = 1e-16

function exact_sol_viscous_shocktube(x,t)
    u   = bisection_solve_velocity(x-v_inf*t,max_iter,tol)
    rho = m_0/u
    e   = 1/(2*γ)*((γ+1)/(γ-1)*v_01^2-u^2)
    return rho, rho*(v_inf+u), rho*(e+1/2*(v_inf+u)^2)
end

U = exact_sol_viscous_shocktube.(x,0.0)
U = ([x[1] for x in U], [x[2] for x in U], [x[3] for x in U])

function flux_differencing!(UF,Uh,Qrh_skew,rxJ,Nh,Nq,Nc,K)
    for k = 1:K
        for j = 1:Nh # col idx
            for i = j:Nh # row idx
                if i <= Nq || j <= Nq # Skip lower right block
                    Fx = euler_fluxes(Uh[1][i,k],Uh[2][i,k],Uh[3][i,k],Uh[1][j,k],Uh[2][j,k],Uh[3][j,k])
                    for c = 1:Nc
                        update_val = 2*rxJ*Qrh_skew[i,j]*Fx[c]
                        UF[c][i,k] += update_val
                        UF[c][j,k] -= update_val
                    end
                end
            end
        end
    end
end

function rhs_inviscid(U,K,VhP,Vq,LIFT,Ph,Qrh_skew,rxJ,mapP,nxJ,rhoL,uL,betaL,rhoR,uR,betaR)
    J = (XR-XL)/K/2 # assume uniform interval
    Nc = 3
    Nh,Nq = size(VhP)

    # Preallocate flux differencing array
    UF = [zeros(Float64,Nh,K) for i = 1:Nc]

    # entropy var projection
    VU = v_ufun((x->Vq*x).(U)...)
    VU = (x->VhP*x).(VU)
    Uh = u_vfun(VU...)

    # convert to rho,u,v,beta vars
    (rho,rhou,E) = Uh
    beta = betafun(rho,rhou,E)
    Ubh = (rho, rhou./rho, beta) # redefine Q = (rho,u,v,β)

    # compute face values
    UM = (x->x[Nq+1:Nh,:]).(Ubh)
    UP = (x->x[mapP]).(UM)
    UP[1][1] = rhoL
    UP[2][1] = uL
    UP[3][1] = betaL
    UP[1][end] = rhoR
    UP[2][end] = uR
    UP[3][end] = betaR

    # simple lax friedrichs dissipation
    Uf =  (x->x[Nq+1:Nh,:]).(Uh)
    (rhoM,rhouM,EM) = Uf
    rhoUM_n = @. (rhouM*nxJ)
    lam = abs.(wavespeed_1D(rhoM,rhoUM_n,EM))
    LFc = .25*max.(lam,lam[mapP])

    fSx = euler_fluxes(UM[1],UM[2],UM[3],UP[1],UP[2],UP[3])
    normal_flux(fx,u) = fx.*nxJ - LFc.*(u[mapP]-u)
    flux = normal_flux.(fSx,Uf)
    rhsU = (x->LIFT*x).(flux)

    flux_differencing!(UF,Ubh,Qrh_skew,rxJ,Nh,Nq,Nc,K)
    rhsU = (x->Ph*x).(UF) .+ rhsU
    rhsU = (x->-x./J).(rhsU)

    return rhsU
end

function rhs_viscous(U,K,N,Vq,Vf,Pq,Dr,LIFT,mapP,nxJ,v1L,v2L,v3L,v1R,v2R,v3R)
    J = (XR-XL)/K/2 # assume uniform interval
    Nc = 3

    # entropy var projection
    VU = v_ufun((x->Vq*x).(U)...)
    VU = (x->Pq*x).(VU)

    # compute and interpolate to quadrature
    VUf = (x->Vf*x).(VU)
    VUP = (x->x[mapP]).(VUf)
    VUP[1][1] = v1L
    VUP[2][1] = v2L
    VUP[3][1] = v3L
    VUP[1][end] = v1R
    VUP[2][end] = v2R
    VUP[3][end] = v3R

    surfx(uP,uf) = LIFT*(@. .5*(uP-uf)*nxJ)
    # Strong form, dv/dx
    VUx = (x->Dr*x).(VU)
    VUx = VUx .+ surfx.(VUP,VUf)
    VUx = VUx./J

    VUx = (x->Vq*x).(VUx)
    VUq = (x->Vq*x).(VU)

    # σ = K dv/dx
    sigma = zero.(VUx)
    for k = 1:K
        for i = 1:size(Vq,1)
            Kx = zeros(3,3)
            v1 = VUq[1][i,k]
            v2 = VUq[2][i,k]
            v4 = VUq[3][i,k]
            Kx[2,2] = -(2*mu-lambda)/v4
            Kx[2,3] = (2*mu-lambda)*v2/v4^2
            Kx[3,2] = Kx[2,3]
            Kx[3,3] = -(2*mu-lambda)*v2^2/v4^3+kappa/cv/v4^2

            sigma[2][i,k] += Kx[2,2]*VUx[2][i,k] + Kx[2,3]*VUx[3][i,k]
            sigma[3][i,k] += Kx[3,2]*VUx[2][i,k] + Kx[3,3]*VUx[3][i,k]
        end
    end
    sigma = (x->Pq*x).(sigma)

    sxf = (x->Vf*x).(sigma)
    sxP = (x->x[mapP]).(sxf)
    sxP[1][1] = sxf[1][1]
    sxP[2][1] = sxf[2][1]
    sxP[3][1] = sxf[3][1]
    sxP[1][end] = sxf[1][end]
    sxP[2][end] = sxf[2][end]
    sxP[3][end] = sxf[3][end]

    # strong form, dσ/dx
    penalization(uP,uf) = LIFT*(@. -.5*(uP-uf))
    sigmax = (x->Dr*x).(sigma)
    sigmax = sigmax .+ surfx.(sxP,sxf)
    sigmax = sigmax./J

    return sigmax
end

function rhs_ESDG(U,K,VhP,Vq,Vf,Pq,Dr,LIFT,Ph,Qrh_skew,rxJ,mapP,nxJ,t)

    rhoL,rhouL,EL = exact_sol_viscous_shocktube(XL,t)
    rhoR,rhouR,ER = exact_sol_viscous_shocktube(XR,t)
    uL = rhouL/rhoL
    uR = rhouR/rhoR
    pL = (γ-1)*(EL-.5*rhoL*uL^2)
    pR = (γ-1)*(ER-.5*rhoR*uR^2)
    betaL = rhoL/(2*pL)
    betaR = rhoR/(2*pR)
    v1L,v2L,v3L = v_ufun(rhoL,rhoL*uL,EL)
    v1R,v2R,v3R = v_ufun(rhoR,rhoR*uR,ER)

    rhsI = rhs_inviscid(U,K,VhP,Vq,LIFT,Ph,Qrh_skew,rxJ,mapP,nxJ,rhoL,uL,betaL,rhoR,uR,betaR)
    rhsV = rhs_viscous(U,K,N,Vq,Vf,Pq,Dr,LIFT,mapP,nxJ,v1L,v2L,v3L,v1R,v2R,v3R)
    return rhsI .+ rhsV
end

# Time stepping
"Time integration"
t = 0.0
U = collect(U)
resU = [zeros(size(x)),zeros(size(x)),zeros(size(x))]
resW = [zeros(size(x)),zeros(size(x)),zeros(size(x))]
resZ = [zeros(size(x)),zeros(size(x)),zeros(size(x))]

Vp = vandermonde_1D(N,LinRange(-1,1,10))/VDM

anim = Animation()
const ptL = XL+(XR-XL)/K/(N+1)/2
const ptR = XR-(XR-XL)/K/(N+1)/2
const hplot = (XR-XL)/K/(N+1)
i = 1

while t < T

    # SSPRK(3,3)
    # dt = min(1e-4,T-t)
    CN = (N+1)*(N+2)/2  # estimated trace constant
    dt = min(2/CN/K^2,T-t)
    rhsU = rhs_ESDG(U,K,VhP,Vq,Vf,Pq,Dr,LIFT,Ph,Qrh_skew,rxJ,mapP,nxJ,t)
    @. resW = U + dt*rhsU
    rhsU = rhs_ESDG(resW,K,VhP,Vq,Vf,Pq,Dr,LIFT,Ph,Qrh_skew,rxJ,mapP,nxJ,t)
    @. resZ = resW+dt*rhsU
    @. resW = 3/4*U+1/4*resZ
    rhsU = rhs_ESDG(resW,K,VhP,Vq,Vf,Pq,Dr,LIFT,Ph,Qrh_skew,rxJ,mapP,nxJ,t)
    @. resZ = resW+dt*rhsU
    @. U = 1/3*U+2/3*resZ

    global t = t + dt
    global i = i + 1
    println("Current time $t with time step size $dt, and final time $T")
    if mod(i,1000) == 1
        plot(x,U[1])
        # for k = 1:K
        #     plot!(ptL+(k-1)*hplot*(N+1):hplot:ptL+k*hplot*(N+1), 1 .-L_plot[:,k],st=:bar,alpha=0.2)
        # end
        frame(anim)
    end

    # if i == 10
    #     break
    # end
end

# rka,rkE,rkc = dopri45_coeffs()
# bcopy!(x,y) = x .= y

# # DOPRI storage
# Utmp = similar.(U)
# rhsUrk = ntuple(x->zero.(U),length(rkE))

# errEst = 0.0
# prevErrEst = 0.0

# t = 0.0
# CN = (N+1)*(N+2)/2  # estimated trace constant
# dt = min(2/CN/K^2,T-t)
# dt0 = dt
# i = 0
# interval = 5

# dthist = Float64[dt]
# thist = Float64[0.0]
# errhist = Float64[0.0]
# wsJ = diagm(wf)

# rhsU = rhs_ESDG(U,K,VhP,Vq,Vf,Pq,Dr,LIFT,Ph,Qrh_skew,rxJ,mapP,nxJ,t)
# bcopy!.(rhsUrk[1],rhsU) # initialize DOPRI rhs (FSAL property)

# while t < T
#     # DOPRI step and
#     rhstest = 0.0
#     for INTRK = 2:7
#         k = zero.(Utmp)
#         for s = 1:INTRK-1
#             bcopy!.(k, @. k + rka[INTRK,s]*rhsUrk[s])
#         end
#         bcopy!.(Utmp, @. U + dt*k)
#         rhsU = rhs_ESDG(Utmp,K,VhP,Vq,Vf,Pq,Dr,LIFT,Ph,Qrh_skew,rxJ,mapP,nxJ,t)#rhsRK(Qtmp,rd,md,ops,euler_fluxes)
#         bcopy!.(rhsUrk[INTRK],rhsU)
#     end
#     errEstVec = zero.(Utmp)
#     for s = 1:7
#         bcopy!.(errEstVec, @. errEstVec + rkE[s]*rhsUrk[s])
#     end

#     errTol = 1e-5
#     errEst = 0.0
#     for field = 1:length(Utmp)
#         errEstScale = @. abs(errEstVec[field]) / (errTol*(1+abs(U[field])))
#         errEst += sum(errEstScale.^2) # hairer seminorm
#     end
#     errEst = sqrt(errEst/(length(U[1])*4))
#     if errEst < 1.0 # if err small, accept step and update
#             bcopy!.(U, Utmp)
#             global t += dt
#             bcopy!.(rhsUrk[1], rhsUrk[7]) # use FSAL property
#     end
#     order = 5
#     dtnew = .8*dt*(.9/errEst)^(.4/(order+1)) # P controller
#     if i > 0 # use PI controller if prevErrEst available
#             dtnew *= (prevErrEst/max(1e-14,errEst))^(.3/(order+1))
#     end
#     global dt = max(min(10*dt0,dtnew),1e-9) # max/min dt
#     global prevErrEst = errEst

#     push!(dthist,dt)
#     push!(thist,t)

#     global i = i + 1  # number of total steps attempted
#     if i%interval==0
#         println("i = $i, t = $t, dt = $dtnew, errEst = $errEst, rhstest = $rhstest")
#     end
# end

# plot(x,U[1])
# gif(anim,"~/Desktop/tmp.gif",fps=15)

xq = Vq*x
exact_U = @. exact_sol_viscous_shocktube.(xq,T)
exact_rho = [x[1] for x in exact_U]
exact_rhou = [x[2] for x in exact_U]
exact_u = exact_rhou./exact_rho
exact_E = [x[3] for x in exact_U]

rho = Vq*U[1]
u = Vq*(U[2]./U[1])
rhou = Vq*U[2]
E = Vq*U[3]
p = Vq*pfun_nd.(U[1],U[2],U[3])
J = (XR-XL)/K/2

# Linferr = maximum(abs.(exact_rho-rho))/maximum(abs.(exact_rho)) +
#           maximum(abs.(exact_u-u))/maximum(abs.(exact_u)) +
#           maximum(abs.(exact_E-E))/maximum(abs.(exact_E))

# L1err = sum(J*wq.*abs.(exact_rho-rho))/sum(J*wq.*abs.(rho)) +
#         sum(J*wq.*abs.(exact_u-u))/sum(J*wq.*abs.(u)) +
#         sum(J*wq.*abs.(exact_E-E))/sum(J*wq.*abs.(E))

Linferr = maximum(abs.(exact_rho-rho))/maximum(abs.(exact_rho)) +
          maximum(abs.(exact_rhou-rhou))/maximum(abs.(exact_rhou)) +
          maximum(abs.(exact_E-E))/maximum(abs.(exact_E))

L1err = sum(J*wq.*abs.(exact_rho-rho))/sum(J*wq.*abs.(rho)) +
        sum(J*wq.*abs.(exact_rhou-rhou))/sum(J*wq.*abs.(rhou)) +
        sum(J*wq.*abs.(exact_E-E))/sum(J*wq.*abs.(E))

L2err = sqrt(sum(J*wq.*abs.(exact_rho-rho).^2))/sqrt(sum(J*wq.*abs.(rho).^2)) +
        sqrt(sum(J*wq.*abs.(exact_rhou-rhou).^2))/sqrt(sum(J*wq.*abs.(rhou).^2)) +
        sqrt(sum(J*wq.*abs.(exact_E-E).^2))/sqrt(sum(J*wq.*abs.(E).^2))

println("N = $N, K = $K")
println("L1 error is $L1err")
println("L2 error is $L2err")
println("Linf error is $Linferr")

max_stars_repo_name: yiminllin/ESDG-CNS

max_issues_repo_name: yiminllin/ESDG-CNS

max_forks_repo_name: yiminllin/ESDG-CNS


text:
using HTTP
using HTTP.Sockets
using Workspace
s = Workspace.Session()
router = Workspace.make_router(s)
HTTP.serve(router, Sockets.localhost, 8081)


text:
module TestRoot
using Test

include("../src/Tree/root.jl")

@testset "root:  tree10                  " begin
    pi = Int[0, 0, 1, 2, 3, 0, 7, 8, 3, 8] .+ 1
    @test find_root(pi) == 1
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
"""
```
B = smf(A; ω, slope, dhₘ, dh₀, cellsize)
```
Applies the simple morphological filter by *Pingel et al. (2013)* [^pingel2013] to `A`.

# Output
- `B::Array{Float64,2}` A filtered version of A

# Arguments
- `A::Array{T,2}` Input Array
- `ω::Float64=18.` Maximum window size [m]
- `slope::Float64=0.01` Terrain slope [m/m]
- `cellsize::Float64=1.` Cellsize in [m]

[^pingel2013]: Pingel, Thomas J., Keith C. Clarke, and William A. McBride. 2013. ‘An Improved Simple Morphological Filter for the Terrain Classification of Airborne LIDAR Data’. ISPRS Journal of Photogrammetry and Remote Sensing 77 (March): 21–30. <https://doi.org/10.1016/j.isprsjprs.2012.12.002>.
"""
function smf(A::AbstractMatrix{T};
    ω::Real = 17.0,
    slope::Real = 0.01,
    cellsize::Real = 1.0) where {T<:Real}

    lastsurface = -copy(A)
    is_low = falses(size(A))
    radii = 3:2:round(Int, ω / cellsize)

    for radius ∈ radii
        threshold = slope * radius * cellsize
        thissurface = opening_circ(lastsurface, radius)
        is_low .|= ((lastsurface - thissurface) .> threshold)
        lastsurface = thissurface
    end

    lastsurface = copy(A)
    is_obj = falses(size(A))
    for radius ∈ radii
        threshold = slope * radius * cellsize
        thissurface = opening_circ(lastsurface, radius)
        is_obj .|= ((lastsurface - thissurface) .> threshold)
        lastsurface = thissurface
    end
    out = Array{Union{Missing,T}}(A)
    out[is_low.|is_obj] .= missing
    return out
end


text:
struct NotFoundError <: Exception
    msg
end

Base.showerror(io::Core.IO, e::NotFoundError) = print(io, "NotFoundError: $(e.msg)")

function convert(field::NTuple{256,UInt8}, ::Type{String})
    return String(filter(x->x != 0, UInt8[field...]))
end

function String(field::NTuple{256,UInt8})
    return String(filter(x->x != 0, UInt8[field...]))
end

toversion(version::Cuint) = VersionNumber(VK_VERSION_MAJOR(version),
										  VK_VERSION_MINOR(version),
										  VK_VERSION_PATCH(version))

function Base.show(io::IO, pdp::VkPhysicalDeviceProperties)
	println(io, "Physical Device Properties: ")
	println(io, "    API Version: ", toversion(pdp.apiVersion))
	println(io, "    Driver Version: ", toversion(pdp.driverVersion))

	println(io, "    Vendor ID ", pdp.vendorID)
	println(io, "    Device ID: ", pdp.deviceID)
	println(io, "    Device Type: ", pdp.deviceType)
	println(io, "    Device Name: ", String(pdp.deviceName))
	println(io, "    Pipeline Cache UUID: ", String(collect(pdp.pipelineCacheUUID)))
	println(io, "    Limits: ", pdp.limits)
	println(io, "    Sparse Properties: \n    ", pdp.sparseProperties)
end

function create_appinfo(; name = "Vulkan Instance", exts = C_NULL, version = VK_MAKE_VERSION(0, 0, 1), engine_name = "NoEngine", engine_version = VK_MAKE_VERSION(0, 0, 1), api_version = VK_API_VERSION_1_0())
    capp_info = Container(VkApplicationInfo(VK_STRUCTURE_TYPE_APPLICATION_INFO, exts, pointer(name), version, pointer(engine_name), engine_version, api_version), Any[Ref(name), Ref(engine_name)])
    return capp_info
end

function available_extensions()
    ext_names = GLFW.GetRequiredInstanceExtensions()
    avail_exts_count = Ref{UInt32}()
    @vkcheck vkEnumerateInstanceExtensionProperties(C_NULL, avail_exts_count, C_NULL)
    avail_exts = Array{VkExtensionProperties}(undef, avail_exts_count[])
    @vkcheck vkEnumerateInstanceExtensionProperties(C_NULL, avail_exts_count, avail_exts)
    return avail_exts
end

function available_layers()
    layer_count = Ref{UInt32}(0)
    @vkcheck vkEnumerateInstanceLayerProperties(layer_count, C_NULL)
    avail_layers = Array{VkLayerProperties}(undef, layer_count[])
    @vkcheck vkEnumerateInstanceLayerProperties(layer_count, avail_layers)
    return avail_layers
end

function check_extensions(ext_names::AbstractArray{T} where T <: AbstractString, avail_exts)
    avail_ext_names = map(x->convert(x.extensionName, String), avail_exts)
    @assert issubset(ext_names, avail_ext_names) "The following extensions are not available: $(symdiff(intersect(avail_layer_names, layer_names), layer_names))"
end

function check_layers(layer_names::AbstractArray{T} where T <: AbstractString, avail_layers)
    avail_layer_names = map(x->convert(x.layerName, String), avail_layers)
    @assert issubset(layer_names, avail_layer_names) "The following layers are not available: $(symdiff(intersect(avail_layer_names, layer_names), layer_names))"
end

@inline function required_extensions(glfw, debug)
    avail_exts = available_extensions()
    ext_names = []
    if glfw
        ext_names = [ext_names..., GLFW.GetRequiredInstanceExtensions()...]
    end
    if debug
        push!(ext_names, "VK_EXT_debug_utils")
    end
    if ext_names != []
        check_extensions(ext_names, available_extensions())
        @info "Enabled extensions:$(map(x->"\n          $x", ext_names)...)"
        exts = Base.unsafe_convert(Ptr{Cstring}, pointer.(ext_names))
        return exts, length(ext_names)
    else
        return C_NULL, 0
    end
end

@inline function required_layers(debug)
    layer_names = []
    if debug
        push!(layer_names, "VK_LAYER_KHRONOS_validation")
    end
    if layer_names != []
        check_layers(Vector{String}(layer_names), available_layers())
        @info "Enabled layers:$(map(x->"\n          $x", layer_names)...)"
        layers = Base.unsafe_convert(Ptr{Cstring}, pointer.(layer_names))
        layers_count = length(layers)
    else
        layers = C_NULL
        layers_count = 0
    end
    return layers, layers_count
end

function create_instance(; glfw = true, debug = true)
    capp_info = create_appinfo()
    exts, exts_count = required_extensions(glfw, debug)
    layers, layers_count = required_layers(debug)
    app_info_ref = Ref(capp_info[])
    cinfo = @safe_contain VkInstanceCreateInfo(VK_STRUCTURE_TYPE_INSTANCE_CREATE_INFO, C_NULL, 0, unsafe_pointer(app_info_ref), layers_count, layers, exts_count, exts)
    append!(cinfo, capp_info)
    inst_ref = Ref{VkInstance}()
    @vkcheck vkCreateInstance(Ref(cinfo[]), C_NULL, inst_ref)
    cinst = Container(inst_ref[], [cinfo])
    finalizer(x->vkDestroyInstance(x[], C_NULL), cinst)
    return cinst
end

function available_physical_devices(inst)
    device_count = Ref{UInt32}(0)
    @vkcheck vkEnumeratePhysicalDevices(inst, device_count, C_NULL)
    devices = Array{VkPhysicalDevice}(undef, device_count[])
    @vkcheck vkEnumeratePhysicalDevices(inst, device_count, devices)
    return devices
end

function properties(devices::AbstractArray{VkPhysicalDevice})
    return map(properties, devices)
end

function properties(device::VkPhysicalDevice)
    device_props = Ref{VkPhysicalDeviceProperties}()
    vkGetPhysicalDeviceProperties(device, device_props)
    return device_props[]
end

function select_physical_device(f, pdps)
    for pdp in pdps
        if Bool(f(pdp))
            return pdp
        end
    end
end

function initialize(; glfw = true, debug = true)
    cinst = create_instance(;glfw = glfw, debug = debug)
    # # not working yet
    if debug
        cmessenger = debug_utils_messenger(cinst[])
        append!(cinst, cmessenger)
        @info "Activated debug callback"
    end
    pdevices = available_physical_devices(cinst[])
    pdps = properties(pdevices)
    @debug "Physical devices found:$(map(x->"\n          $(String(x.deviceName))", pdps)...)"
    pdp = select_physical_device(pdp->pdp.deviceType == VK_PHYSICAL_DEVICE_TYPE_DISCRETE_GPU, pdps)
    @info "Selected device \"$(String(pdp.deviceName))\""
    pdevice = pdevices[findfirst(x->x == pdp, pdps)]

    cinst, pdevice
end


text:
import Telegrambot, Telegrambot.InlineQueryResultArticle
import UUIDs
import Scryfall
botApi = "bot552445015:AAFm1fTBaFO9jK42fuMdQe_OMtMj7pbNr9o"

function welcomeMsg(name::AbstractString)
    msg = """
    Welcome to using Scryfall bot, powered by https://github.com/Moelf/Telegrambot.jl
    Basic usage:
    - /searchtext <fuzzy_card_name>
    - /searchimg <fuzzy_card_name>
    """
    return msg
end

function getImgVer(fuzzyName::AbstractString)
    nameList = split(fuzzyName, " ")
    if occursin("ver=", nameList[end])
        vercode = string(fuzzyName[end-2:end])
        name = join(nameList[1:end-1])
        print("name: $name, vercode: $vercode")
    else
        name = fuzzyName
        vercode=" "
    end
    #= print(name, verCode) =#
    return Scryfall.getImgurl(name; setCode=vercode)
end

function echo(incoming::AbstractString)
    return incoming
end

txtCmds = Dict()
txtCmds["searchtext"] = Scryfall.getOracle
txtCmds["searchimg"] = getImgVer
txtCmds["start"] = welcomeMsg
txtCmds["repeat_msg"] = echo #this will respond to '/repeat_msg <any thing>'
inlineOpts = Dict() #Title, result pair
inlineOpts["Make Uppercase"] = uppercase #this will generate an pop-up named Make Uppercase and upon tapping return uppercase(<user_input>)
Telegrambot.startBot(botApi; textHandle = txtCmds, inlineQueryHandle=inlineOpts)


text:
include("split_explicit.jl")

function SplitConfig(
    name,
    resolution,
    dimensions,
    coupling,
    rotation = Fixed();
    boundary_conditions = (
        OceanBC(Impenetrable(FreeSlip()), Insulating()),
        OceanBC(Penetrable(FreeSlip()), Insulating()),
    ),
    solver = SplitExplicitSolver,
    dt_slow = 90 * 60,
)
    mpicomm = MPI.COMM_WORLD
    ArrayType = ClimateMachine.array_type()

    N, Nˣ, Nʸ, Nᶻ = resolution
    Lˣ, Lʸ, H = dimensions

    xrange = range(FT(0); length = Nˣ + 1, stop = Lˣ)
    yrange = range(FT(0); length = Nʸ + 1, stop = Lʸ)
    zrange = range(FT(-H); length = Nᶻ + 1, stop = 0)

    brickrange_2D = (xrange, yrange)
    topl_2D = BrickTopology(
        mpicomm,
        brickrange_2D,
        periodicity = (true, true),
        boundary = ((0, 0), (0, 0)),
    )
    grid_2D = DiscontinuousSpectralElementGrid(
        topl_2D,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )

    brickrange_3D = (xrange, yrange, zrange)
    topl_3D = StackedBrickTopology(
        mpicomm,
        brickrange_3D;
        periodicity = (true, true, false),
        boundary = ((0, 0), (0, 0), (1, 2)),
    )
    grid_3D = DiscontinuousSpectralElementGrid(
        topl_3D,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )

    problem = SimpleBox{FT}(
        dimensions...;
        BC = boundary_conditions,
        rotation = rotation,
    )

    dg_3D, dg_2D = setup_models(
        solver,
        problem,
        grid_3D,
        grid_2D,
        param_set,
        coupling,
        dt_slow,
    )

    return SplitConfig(name, dg_3D, dg_2D, solver, mpicomm, ArrayType)
end

function setup_models(
    ::Type{SplitExplicitSolver},
    problem,
    grid_3D,
    grid_2D,
    param_set,
    coupling,
    _,
)

    model_3D = HydrostaticBoussinesqModel{FT}(
        param_set,
        problem;
        coupling = coupling,
        cʰ = FT(1),
        αᵀ = FT(0),
        κʰ = FT(0),
        κᶻ = FT(0),
    )

    model_2D = ShallowWaterModel{FT}(
        param_set,
        problem,
        ShallowWater.ConstantViscosity{FT}(model_3D.νʰ),
        nothing;
        coupling = coupling,
        c = FT(1),
    )

    N = polynomialorders(grid_3D)
    Nvert = N[end]
    vert_filter = CutoffFilter(grid_3D, Nvert - 1)
    exp_filter = ExponentialFilter(grid_3D, 1, 8)

    integral_model = DGModel(
        VerticalIntegralModel(model_3D),
        grid_3D,
        CentralNumericalFluxFirstOrder(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    dg_2D = DGModel(
        model_2D,
        grid_2D,
        CentralNumericalFluxFirstOrder(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    Q_2D = init_ode_state(dg_2D, FT(0); init_on_cpu = true)

    modeldata = (
        dg_2D = dg_2D,
        Q_2D = Q_2D,
        vert_filter = vert_filter,
        exp_filter = exp_filter,
        integral_model = integral_model,
    )

    dg_3D = DGModel(
        model_3D,
        grid_3D,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient();
        modeldata = modeldata,
    )

    return dg_3D, dg_2D

end

function setup_models(
    ::Type{SplitExplicitLSRK2nSolver},
    problem,
    grid_3D,
    grid_2D,
    param_set,
    _,
    dt_slow,
)
    add_fast_substeps = 2
    numImplSteps = 5
    numImplSteps > 0 ? ivdc_dt = dt_slow / FT(numImplSteps) : ivdc_dt = dt_slow
    model_3D = OceanModel{FT}(
        param_set,
        problem,
        cʰ = FT(1),
        αᵀ = FT(0),
        κʰ = FT(0),
        κᶻ = FT(0),
        add_fast_substeps = add_fast_substeps,
        numImplSteps = numImplSteps,
        ivdc_dt = ivdc_dt,
    )

    model_2D = BarotropicModel(model_3D)

    dg_2D = DGModel(
        model_2D,
        grid_2D,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    Q_2D = init_ode_state(dg_2D, FT(0); init_on_cpu = true)

    dg_3D = OceanDGModel(
        model_3D,
        grid_3D,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient();
        modeldata = (dg_2D, Q_2D),
    )

    return dg_3D, dg_2D

end


text:
module MessagePassingTest

using Test
using ForneyLab
using ForneyLab: generateId, addNode!, associate!, summaryPropagationSchedule, <<, matches, flatten, Cluster

@testset "Message" begin
    msg = Message(Univariate, GaussianMeanVariance, m=0.0, v=1.0)
    @test isa(msg, Message{GaussianMeanVariance})
    @test isa(msg, Message{GaussianMeanVariance, Univariate})
    @test !isa(msg, Message{PointMass})
    @test !isa(msg, Message{GaussianMeanPrecision})
    @test !isa(msg, Message{GaussianMeanVariance, Multivariate})
    @test !isa(msg, Message{GaussianMeanVariance, MatrixVariate})
    @test msg.dist == ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0)

    @test Message(Univariate, PointMass, m=0.0) == Message(Univariate, PointMass, m=0.0)
    @test Message(Univariate, PointMass, m=0.0) != Message(Univariate, GaussianMeanVariance, m=0.0, v=1.0)
    @test Message(Univariate, PointMass, m=0.0) != Message(Univariate, PointMass, m=1.0)
end

@testset "<<" begin
    @test <<(Message{Gaussian}, Message)
    @test <<(Message{Gaussian, Univariate}, Message)
    @test <<(Message{Gaussian, Univariate}, Message{Gaussian, Univariate})
    @test <<(Message{GaussianMeanVariance, Univariate}, Message{Gaussian, Univariate})
    @test !<<(Message{GaussianMeanVariance, Univariate}, Message{Gaussian, Multivariate})
    @test !<<(Message{GaussianMeanVariance, Univariate}, Message{PointMass, Univariate})
    @test <<(Message{GaussianMeanVariance, Univariate}, Message{Gaussian})
    @test <<(Message{Gaussian}, Message{Gaussian})
    @test <<(Message{GaussianMeanVariance}, Message{Gaussian})
    @test !<<(Nothing, Message{Gaussian})
    @test !<<(Nothing, Message{GaussianMeanVariance})
    @test <<(Message{Gamma, Univariate}, Message{Union{Gamma, Wishart}, Univariate})
    @test <<(Message{Gamma}, Message{Union{Gamma, Wishart}})
    @test !<<(Message{GaussianMeanVariance, Univariate}, ProbabilityDistribution{Univariate, GaussianMeanVariance})
    @test !<<(ProbabilityDistribution{Univariate, GaussianMeanVariance}, Message{GaussianMeanVariance, Univariate})
end

@testset "matches" begin
    @test matches(Message{Gaussian}, Message{FactorFunction})
    @test matches(Message{FactorFunction}, Message{Gaussian})
    @test !matches(Nothing, Message{Gaussian})
end

# Integration helper
mutable struct MockNode <: FactorNode
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Int,Interface}

    function MockNode(vars::Vector{Variable}; id=generateId(MockNode))
        n_interfaces = length(vars)
        self = new(id, Array{Interface}(undef, n_interfaces), Dict{Int,Interface}())
        addNode!(currentGraph(), self)

        for idx = 1:n_interfaces
            self.i[idx] = self.interfaces[idx] = associate!(Interface(self), vars[idx])
        end

        return self
    end
end

@testset "ScheduleEntry" begin
    FactorGraph()
    nd = MockNode([Variable(), Variable()])

    # Schedule should be constructable by hand
    schedule = [ScheduleEntry(nd.i[1], Nothing), ScheduleEntry(nd.i[2], Nothing)]
    @test isa(schedule, Schedule)
end

@testset "summaryPropagationSchedule" begin
    # Construct
    #
    #    a    b    c
    # [1]--[2]--[3]--[4]
    #            |d
    #           [5]

    FactorGraph()

    a = Variable()
    b = Variable()
    c = Variable()
    d = Variable()

    n1 = MockNode([a])
    n2 = MockNode([a, b])
    n3 = MockNode([b, c, d])
    n4 = MockNode([c])
    n5 = MockNode([d])

    cl = Cluster(n1, [a.edges; b.edges])

    schedule = summaryPropagationSchedule(Variable[d], Cluster[cl])

    @test length(schedule) == 7
    @test ScheduleEntry(n1.i[1], Nothing) in schedule
    @test ScheduleEntry(n2.i[2], Nothing) in schedule
    @test ScheduleEntry(n4.i[1], Nothing) in schedule
    @test ScheduleEntry(n3.i[3], Nothing) in schedule
    @test ScheduleEntry(n5.i[1], Nothing) in schedule
    @test ScheduleEntry(n3.i[1], Nothing) in schedule
    @test ScheduleEntry(n2.i[1], Nothing) in schedule
end

@testset "flatten" begin
    FactorGraph()

    a = Variable()
    b = Variable()
    c = Variable()

    n1 = MockNode([a])
    n2 = MockNode([a, b])
    n3 = MockNode([b, c])
    n4 = MockNode([Variable()])

    schedule = [ScheduleEntry(n1.i[1], Nothing),
                ScheduleEntry(n2.i[2], Nothing),
                ScheduleEntry(n3.i[2], Nothing)]
    schedule[2].internal_schedule = [   ScheduleEntry(n1.i[1], Nothing),
                                        ScheduleEntry(n4.i[1], Nothing)]

    flat_schedule = flatten(schedule)

    @test length(flat_schedule) == 3
    @test flat_schedule[1] == schedule[1]
    @test flat_schedule[2] == schedule[2].internal_schedule[2]
    @test flat_schedule[3] == schedule[3]
end

end # module


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function test_ellipse(fname)
    pagewidth, pageheight = 1200, 1400
    Drawing(pagewidth, pageheight, fname)
    origin() # move 0/0 to center
    background("ivory")
    setopacity(0.5)
    setline(0.3)

    @test_throws ErrorException Partition(0, pageheight, 200, 200)
    @test_throws ErrorException Partition(pagewidth, 0, 200, 200)

    pagetiles = Partition(pagewidth, pageheight, 200, 200)
    for (pos, n) in pagetiles
        randomhue()
        ellipse(pos, rand(50:pagetiles.tilewidth), rand(50:pagetiles.tileheight), :strokepreserve)
        clip()
        for i in 1:10
            ellipse(pos, rand(50:pagetiles.tilewidth), rand(50:pagetiles.tileheight), :fill)
            randomhue()
            ellipse(pos, Point(pos.x + rand(50:pagetiles.tilewidth), pos.y), rand(50:pagetiles.tileheight), :stroke)
        end
        setline(1)
        sethue("black")
        fillstroke()
        ellipse(pos, 5, 5, :fill)
        clipreset()
        # test three point version
        circle.((pos - (50, 0), pos + (50, 0), pos + (0, 2n)), 4, :fill)
        ellipse(pos - (50, 0), pos + (50, 0), pos + (0, 2n), :stroke)
    end

    # test ellipseinquad
    pg = ngon(O, 250, 6, π/6, vertices=true)

    top = vcat(O, pg[[3, 4, 5]])
    left = vcat(O, pg[[1, 2, 3]])
    right = vcat(O, pg[[5, 6, 1]])

    for (n, side) in enumerate((top, left, right))
        sethue([Luxor.julia_green, Luxor.julia_red, Luxor.julia_purple, Luxor.julia_blue][n])
        poly(side, :path, reversepath=true, close=true)
        results = ellipseinquad(polyscale!(side, 0.5, center = polycentroid(side)), :fill)
        @test results[1] != O
        @test !iszero(results[2])
        @test !iszero(results[3])
        @test !iszero(results[4])
    end

    @test finish() == true
end

fname = "ellipse-test1.png"
test_ellipse(fname)
println("...finished test: output in $(fname)")


text:
abstract type AbstractSimpleBoxProblem <: AbstractOceanProblem end

"""
    ocean_init_aux!(::HBModel, ::AbstractSimpleBoxProblem)

save y coordinate for computing coriolis, wind stress, and sea surface temperature

# Arguments
- `m`: model object to dispatch on and get viscosities and diffusivities
- `p`: problem object to dispatch on and get additional parameters
- `A`: auxiliary state vector
- `geom`: geometry stuff
"""
function ocean_init_aux!(::HBModel, ::AbstractSimpleBoxProblem, A, geom)
    FT = eltype(A)
    @inbounds A.y = geom.coord[2]

    # needed for proper CFL condition calculation
    A.w = -0
    A.pkin = -0
    A.wz0 = -0

    A.uᵈ = @SVector [-0, -0]
    A.ΔGᵘ = @SVector [-0, -0]

    return nothing
end

function ocean_init_aux!(::OceanModel, ::AbstractSimpleBoxProblem, A, geom)
    FT = eltype(A)
    @inbounds A.y = geom.coord[2]

    # needed for proper CFL condition calculation
    A.w = -0
    A.pkin = -0
    A.wz0 = -0

    A.u_d = @SVector [-0, -0]
    A.ΔGu = @SVector [-0, -0]

    return nothing
end

function ocean_init_aux!(::SWModel, ::AbstractSimpleBoxProblem, A, geom)
    @inbounds A.y = geom.coord[2]

    A.Gᵁ = @SVector [-0, -0]
    A.Δu = @SVector [-0, -0]

    return nothing
end

function ocean_init_aux!(::BarotropicModel, ::AbstractSimpleBoxProblem, A, geom)
    @inbounds A.y = geom.coord[2]

    A.Gᵁ = @SVector [-0, -0]
    A.U_c = @SVector [-0, -0]
    A.η_c = -0
    A.U_s = @SVector [-0, -0]
    A.η_s = -0
    A.Δu = @SVector [-0, -0]
    A.η_diag = -0
    A.Δη = -0

    return nothing
end

"""
    coriolis_parameter

northern hemisphere coriolis

# Arguments
- `m`: model object to dispatch on and get coriolis parameters
- `y`: y-coordinate in the box
"""
@inline coriolis_parameter(
    m::Union{HBModel, OceanModel},
    ::AbstractSimpleBoxProblem,
    y,
) = m.fₒ + m.β * y
@inline coriolis_parameter(
    m::Union{SWModel, BarotropicModel},
    ::AbstractSimpleBoxProblem,
    y,
) = m.fₒ + m.β * y

############################
# Basic box problem        #
# Set up dimensions of box #
############################

abstract type AbstractRotation end
struct Rotating <: AbstractRotation end
struct Fixed <: AbstractRotation end

"""
    SimpleBoxProblem <: AbstractSimpleBoxProblem

Stub structure with the dimensions of the box.
Lˣ = zonal (east-west) length
Lʸ = meridional (north-south) length
H  = height of the ocean
"""
struct SimpleBox{R, T, BC} <: AbstractSimpleBoxProblem
    rotation::R
    Lˣ::T
    Lʸ::T
    H::T
    boundary_conditions::BC
    function SimpleBox{FT}(
        Lˣ, # m
        Lʸ, # m
        H;  # m
        rotation = Fixed(),
        BC = (
            OceanBC(Impenetrable(FreeSlip()), Insulating()),
            OceanBC(Penetrable(FreeSlip()), Insulating()),
        ),
    ) where {FT <: AbstractFloat}
        return new{typeof(rotation), FT, typeof(BC)}(rotation, Lˣ, Lʸ, H, BC)
    end
end

@inline coriolis_parameter(
    m::Union{HBModel, OceanModel},
    ::SimpleBox{R},
    y,
) where {R <: Fixed} = -0
@inline coriolis_parameter(
    m::Union{SWModel, BarotropicModel},
    ::SimpleBox{R},
    y,
) where {R <: Fixed} = -0

@inline coriolis_parameter(
    m::Union{HBModel, OceanModel},
    ::SimpleBox{R},
    y,
) where {R <: Rotating} = m.fₒ
@inline coriolis_parameter(
    m::Union{SWModel, BarotropicModel},
    ::SimpleBox{R},
    y,
) where {R <: Rotating} = m.fₒ

function ocean_init_state!(
    m::Union{SWModel, BarotropicModel},
    p::SimpleBox,
    Q,
    A,
    local_geometry,
    t,
)
    coords = local_geometry.coord

    k = (2π / p.Lˣ, 2π / p.Lʸ, 2π / p.H)
    ν = viscosity(m)

    gH = gravity_speed(m)
    @inbounds f = coriolis_parameter(m, p, coords[2])

    U, V, η = barotropic_state!(p.rotation, (coords..., t), ν, k, (gH, f))

    Q.U = @SVector [U, V]
    Q.η = η

    return nothing
end

viscosity(m::SWModel) = (m.turbulence.ν, m.turbulence.ν, -0)
viscosity(m::BarotropicModel) = (m.baroclinic.νʰ, m.baroclinic.νʰ, -0)

gravity_speed(m::SWModel) = grav(parameter_set(m)) * m.problem.H
gravity_speed(m::BarotropicModel) =
    grav(parameter_set(m)) * m.baroclinic.problem.H

function ocean_init_state!(
    m::Union{HBModel, OceanModel},
    p::SimpleBox,
    Q,
    A,
    local_geometry,
    t,
)

    coords = local_geometry.coord

    k = (2π / p.Lˣ, 2π / p.Lʸ, 2π / p.H)
    ν = (m.νʰ, m.νʰ, m.νᶻ)

    gH = grav(parameter_set(m)) * p.H
    @inbounds f = coriolis_parameter(m, p, coords[2])

    U, V, η = barotropic_state!(p.rotation, (coords..., t), ν, k, (gH, f))
    u°, v° = baroclinic_deviation(p.rotation, (coords..., t), ν, k, f)

    u = u° + U / p.H
    v = v° + V / p.H

    Q.u = @SVector [u, v]
    Q.η = η
    Q.θ = -0

    return nothing
end

function barotropic_state!(
    ::Fixed,
    (x, y, z, t),
    (νˣ, νʸ, νᶻ),
    (kˣ, kʸ, kᶻ),
    params,
)
    gH, _ = params

    M = @SMatrix [-νˣ*kˣ^2 gH*kˣ; -kˣ 0]
    A = exp(M * t) * @SVector [1, 1]

    U = A[1] * sin(kˣ * x)
    V = -0
    η = A[2] * cos(kˣ * x)

    return (U = U, V = V, η = η)
end

function baroclinic_deviation(
    ::Fixed,
    (x, y, z, t),
    (νˣ, νʸ, νᶻ),
    (kˣ, kʸ, kᶻ),
    f,
)
    λ = νˣ * kˣ^2 + νᶻ * kᶻ^2

    u° = exp(-λ * t) * cos(kᶻ * z) * sin(kˣ * x)
    v° = -0

    return (u° = u°, v° = v°)
end

function barotropic_state!(
    ::Rotating,
    (x, y, z, t),
    (νˣ, νʸ, νᶻ),
    (kˣ, kʸ, kᶻ),
    params,
)
    gH, f = params

    M = @SMatrix [-νˣ*kˣ^2 f gH*kˣ; -f -νˣ*kˣ^2 0; -kˣ 0 0]
    A = exp(M * t) * @SVector [1, 1, 1]

    U = A[1] * sin(kˣ * x)
    V = A[2] * sin(kˣ * x)
    η = A[3] * cos(kˣ * x)

    return (U = U, V = V, η = η)
end

function baroclinic_deviation(
    ::Rotating,
    (x, y, z, t),
    (νˣ, νʸ, νᶻ),
    (kˣ, kʸ, kᶻ),
    f,
)
    λ = νˣ * kˣ^2 + νᶻ * kᶻ^2

    M = @SMatrix[-λ f; -f -λ]
    A = exp(M * t) * @SVector[1, 1]

    u° = A[1] * cos(kᶻ * z) * sin(kˣ * x)
    v° = A[2] * cos(kᶻ * z) * sin(kˣ * x)

    return (u° = u°, v° = v°)
end

@inline kinematic_stress(p::SimpleBox, y) = @SVector [-0, -0]