text:
@with_kw struct Tolerances{Txtol, Tftol, Tgrtol, Tkktol}
    xtol::Txtol = 0.0
    ftol::Tftol = 1e-6
    grtol::Tgrtol = 0.0
    kkttol::Tkktol = 1e-6
end
function Base.:*(t::Tolerances, m::Real)
    return Tolerances(t.xtol * m, t.ftol * m, t.grtol * m, t.kkttol * m)
end
function (tol::Tolerances{<:Function, <:Function, <:Function, <:Function})(i)
    return Tolerances(tol.xtol(i), tol.ftol(i), tol.grtol(i), tol.kkttol(i))
end

struct DualData{TSol, Tv}
    λ::TSol
    l::Tv
    u::Tv
end
function DualData(model::Model{T, TV}) where {T, TV}
    n_i = length(constraints(model))
    λ = DualSolution(model)
    # Initial data and bounds for Optim to solve dual problem
    l = Vector(zerosof(TV, n_i))
    u = Vector(infsof(TV, n_i))
    return DualData(λ, l, u)
end

@with_kw mutable struct Options{T, Ttol <: Tolerances, TSubOptions <: Optim.Options}
    maxiter::Int = 1000
    outer_maxiter::Int = 10^8
    tol::Ttol = Tolerances()
    s_init::T = 0.5
    s_incr::T = 1.2
    s_decr::T = 0.7
    dual_caps::Tuple{T, T} = MMA.default_dual_caps(Float64)
    store_trace::Bool = false
    show_trace::Bool = false
    extended_trace::Bool = false
    subopt_options::TSubOptions = Optim.Options(allow_f_increases = false)
end
@inline function Base.getproperty(o::Options, f::Symbol)
    f === :xtol && return o.tol.xtol
    f === :ftol && return o.tol.ftol
    f === :grtol && return o.tol.grtol
    f === :kkttol && return o.tol.kkttol
    return getfield(o, f)
end

mutable struct Workspace{T, Tx, Tc, TO, TSO, TTrace <: OptimizationTrace{T}, TModel <: Model{T}, TPD <: PrimalData{T}, TXUpdater <: XUpdater{T, Tx, TPD}, Tdual_data <: DualData, TOptions, TConvCriteria, TState}
    model::TModel
    optimizer::TO
    suboptimizer::TSO
	primal_data::TPD
	asymptotes_updater::AsymptotesUpdater{T, Tx, TModel}
	variable_bounds_updater::VariableBoundsUpdater{T, Tx, TPD, TModel}
	cvx_grad_updater::ConvexApproxGradUpdater{T, Tx, TPD, TModel}
	lift_updater::LiftUpdater{T, Tc, TPD}
	lift_resetter::LiftResetter{T, Tc}
	x_updater::TXUpdater
    dual_data::Tdual_data
    dual_obj::DualObjVal{TPD, TXUpdater}
	dual_obj_grad::DualObjGrad{TPD, TXUpdater}
    tracing::Bool
    tr::TTrace
    outer_iter::Int
    iter::Int
	f_calls::Int
	g_calls::Int
    options::TOptions
    convcriteria::TConvCriteria
    convstate::TState
end
const Workspace87{T, TV1, TV2, TM, TSO, TModel, TPD} = Workspace{T, TV1, TV2, TM, MMA87, TSO, TModel, TPD}
const Workspace02{T, TV1, TV2, TM, TSO, TModel, TPD} = Workspace{T, TV1, TV2, TM, MMA02, TSO, TModel, TPD}

Base.show(io::IO, w::Workspace87) = print(io, "Workspace for the method of moving asymptotes of 1987.")
Base.show(io::IO, w::Workspace02) = print(io, "Workspace for the method of moving asymptotes of 2002.")

function PrimalData(model::Model{T, TV}, optimizer, x0) where {T, TV}
    n_i = length(constraints(model))
    n_j = dim(model)

    x, x1, x2 = copy(x0), copy(x0), copy(x0)
    α, β, σ = zerosof(TV, n_j), zerosof(TV, n_j), zerosof(TV, n_j)

    # Buffers for p0, pji, q0, qji
    p0, q0 = zerosof(TV, n_j), zerosof(TV, n_j)
    TM = MatrixOf(TV)
    p, q = zerosof(TM, n_j, n_i), zerosof(TM, n_j, n_i)
    if optimizer isa MMA87
        ρ = Vector(zerosof(TV, n_i))
    else
        ρ = Vector(onesof(TV, n_i))
    end
    r = Vector(zerosof(TV, n_i))
    ∇f_x = nansof(TV, length(x))
    f_x = T(NaN)
    f_x_previous = T(NaN)
    g = Vector(zerosof(TV, n_i))
    ∇g = zerosof(TM, n_j, n_i)

    return PrimalData(  σ, α, β, p0, q0, p, q, ρ, r,
                        zero(T), x0, x, x1, x2, f_x,
                        f_x_previous, g, ∇f_x, ∇g
                    )
end

function AsymptotesUpdater(primal_data, model, options)
    @unpack σ, x, x1, x2 = primal_data
    @unpack s_init, s_incr, s_decr = options
    return AsymptotesUpdater(model, σ, x, x1, x2, s_init, s_incr, s_decr)
end

function LiftUpdater(primal_data, model::Model{T, TV}) where {T, TV}
    @unpack ρ, g = primal_data
    n_i = length(constraints(model))
    n_j = dim(model)
    ng_approx = Vector(zerosof(TV, n_i))
    return LiftUpdater(primal_data, ρ, g, ng_approx, n_j)
end

function DualSolution(model::Model{T, TV}) where {T, TV}
    dev = whichdevice(model)
    n_i = length(constraints(model))
    if dev isa CPU
        return DualSolution{:CPU}(Vector(onesof(TV, n_i)))
    else
        return DualSolution{:GPU}(Vector(onesof(TV, n_i)))
    end
end

function Workspace( model::Model{T, TV},
                    x0::TV,
                    optimizer::TO = MMA02(),
                    suboptimizer = Optim.ConjugateGradient();
                    options = Options(),
                    convcriteria = KKTCriteria()
                ) where {T, TV, TO}

    primal_data = PrimalData(model, optimizer, x0)
    asymptotes_updater = AsymptotesUpdater(primal_data, model, options)
    variable_bounds_updater = VariableBoundsUpdater(primal_data, model, T(μ))
    cvx_grad_updater = ConvexApproxGradUpdater(primal_data, model)
    lift_updater = LiftUpdater(primal_data, model)
    x_updater = XUpdater(primal_data)
    dual_data = DualData(model)
    dual_obj = DualObjVal(primal_data, x_updater, dual_data.λ)
    dual_obj_grad = DualObjGrad(primal_data, x_updater, dual_data.λ)
    lift_resetter = LiftResetter(primal_data.ρ, T(ρmin))

    f_calls, g_calls = 0, 0
    convstate = ConvergenceState()

    # Evaluate the constraints and their gradients
    tr = OptimizationTrace{T, TO}()
    tracing = (options.store_trace || options.extended_trace || options.show_trace)

    # Iteraton counter
    outer_iter = 0
    iter = 0

    workspace = Workspace(
        model, optimizer, suboptimizer, primal_data, asymptotes_updater,
        variable_bounds_updater, cvx_grad_updater, lift_updater,
        lift_resetter, x_updater, dual_data, dual_obj, dual_obj_grad,
        tracing, tr, outer_iter, iter, f_calls, g_calls, options, convcriteria,
        convstate
    )
    update_values!(workspace)
    workspace.convstate = assess_convergence(workspace)

    return workspace
end

function update_constraints!(g, ∇g, model, x)
    for i in 1:length(g)
        g[i] = eval_constraint(model, i, x, @view(∇g[:,i]))
    end
    return g
end

@with_kw mutable struct ConvergenceState{T}
    x_converged::Bool = false
    f_converged::Bool = false
    gr_converged::Bool = false
    kkt_converged::Bool = false
    x_residual::T = Inf
    f_residual::T = Inf
    gr_residual::T = Inf
    kkt_residual::T = Inf
    f_increased::Bool = false
    converged::Bool = false
end
function ConvergenceState(::Type{T}) where T
    ConvergenceState(false, false, false, false, T(Inf), T(Inf), T(Inf), T(Inf), false, false)
end

function assess_convergence(workspace::Workspace)
    @unpack options, primal_data, dual_data, model, convcriteria = workspace
    @unpack x, x1, x2, f_x, f_x_previous, ∇f_x, ∇g, g = primal_data
    @unpack λ = dual_data
    @unpack box_max, box_min = model
    @unpack xtol, ftol, grtol, kkttol = options.tol

    if x isa CuArray
        x_residual = mapreduce((x1, x2) -> abs(x1 - x2), max, x, x1, init=zero(T))
    else
        x_residual = maxdiff(x, x1)
    end
    f_residual = abs(f_x - f_x_previous)
    gr_residual = maximum(abs, ∇f_x)
    kkt_residual = get_kkt_residual(∇f_x, g, ∇g, λ.cpu, x, box_min, box_max)

    x_converged = x_residual < xtol
    f_converged = f_residual / (abs(f_x) + ftol) < ftol
    gr_converged = gr_residual < grtol
    kkt_converged = kkt_residual < kkttol
    f_increased = f_x > f_x_previous

    if convcriteria isa DefaultCriteria
        converged = (x_converged || f_converged) && all(x -> x <= 0, g)
    else
        converged = kkt_converged# && all(x -> x <= 0, g)
    end

    return ConvergenceState(    x_converged,
                                f_converged,
                                gr_converged,
                                kkt_converged,
                                x_residual,
                                f_residual,
                                gr_residual,
                                kkt_residual,
                                f_increased,
                                converged
                            )
end


text:
module CVRPLIB

using DataStructures
using Match
using Downloads
using Random
using DelimitedFiles
import TSPLIB

    mutable struct CVRP
        name            :: String
        dimension       :: Int
        weight_type     :: String
        weights         :: Matrix{Int}
        capacity        :: Int
        distance        :: Float64
        service_time    :: Float64
        coordinates     :: Matrix{Float64}
        demand          :: Vector{Int}
        depot           :: Int
        dummy           :: Int
        customers       :: Vector{Int}
    end

    const cvrp_keys = [
        "NAME",
        "TYPE",
        "COMMENT",
        "DIMENSION",
        "EDGE_WEIGHT_TYPE",
        "EDGE_WEIGHT_FORMAT",
        "EDGE_DATA_FORMAT",
        "CAPACITY",
        "DISTANCE",
        "SERVICE_TIME",
        "NODE_COORD_TYPE",
        "DISPLAY_DATA_TYPE",
        "NODE_COORD_SECTION",
        "DEPOT_SECTION",
        "DEMAND_SECTION",
        "EDGE_DATA_SECTION",
        "FIXED_EDGES_SECTION",
        "DISPLAY_DATA_SECTION",
        "TOUR_SECTION",
        "EDGE_WEIGHT_SECTION",
        "EOF"
    ]

    include("reader.jl")
    include("writer.jl")
    include("download.jl")

    export readCVRP, readCVRPLIB, CVRP, write_cvrp, delete_cvrp

end

ext: jl


text:
module __JuliaCLIInit

project, = ARGS

@debug "Loading packages..."

import JSONRPC
import Logging
using Pkg: Pkg, TOML
using Sockets: listen

include("../rpcutils.jl")
include("../evalutils.jl")

if lowercase(get(ENV, "JULIA_CLI_REVISE", "false")) == "true"
    try
        @eval begin
            using Revise: Revise, revise
            revise_success() = isempty(Revise.queue_errors)
        end
    catch err
        @error "using Revise" exception = (err, catch_backtrace())
    end
end
if !@isdefined revise
    revise() = nothing
    revise_success() = true
end

@debug "Loading packages...DONE"

function apicall(f, commonopts)
    cwd = commonopts.cwd
    ioconfig = commonopts.ioconfig
    ignorereturn = commonopts.ignorereturn
    (result, backtrace), out, err = with_stdio(ioconfig) do
        with_newlogger() do

            revise()
            if !revise_success()
                exit(113)  # TODO: don't
                return JSONRPC.Error(ReviseError, "Revise failed"), nothing
            end

            try
                cd(f, cwd), nothing
            catch err
                err, sprint(showerror, err, catch_backtrace())
            end
        end
    end
    if backtrace !== nothing || result isa Exception
        if result isa LoadError
            result = result.error
        end
        data = (stdout = out, stderr = err, backtrace = backtrace)
        return JSONRPC.usererror(sprint(showerror, result), data)
    end
    return (result = ignorereturn ? nothing : result, stdout = out, stderr = err)
end

function getcommonoptions(input)
    cwd = get(input, "cwd", pwd())
    cwd isa AbstractString || return ArgumentError("Keyword argument `cwd` must be of type `String`.")

    _ioconfig = @return_exception ioconfig(input)

    ignorereturn = get(input, "ignorereturn", false)
    ignorereturn isa Bool || return ArgumentError("Keyword argument `ignorereturn` must be a `nothing` or `Bool`.")

    commonopts = (cwd = cwd, ioconfig = _ioconfig, ignorereturn = ignorereturn)
    return commonopts
end

function eval_api(input)
    @require_keyword_arguments input
    @debug "eval_api" input

    code = get(input, "code", nothing)
    code isa AbstractString || return ArgumentError("Keyword argument `code` of type `String` is required.")

    usemain = get(input, "usemain", false)
    usemain isa Bool || return ArgumentError("Keyword argument `usemain` must be of type `Bool`.")

    args = @return_exception getstrings(input, "args")
    commonopts = @return_exception getcommonoptions(input)

    namespace = usemain ? Main : Module()
    @debug "Evaluating" code = Text(code) args namespace
    return apicall(commonopts) do
        append!(empty!(Base.ARGS), args)
        Base.include_string(namespace, code)
    end
end

function getcallable(input)
    pkgname = get(input, "pkgname", nothing)
    pkgname isa AbstractString || return ArgumentError("Keyword argument `pkgname` of type `String` is required.")

    pkguuid = get(input, "pkguuid", nothing)
    pkguuid isa AbstractString || return ArgumentError("Keyword argument `pkguuid` of type `String` is required.")
    pkguuid = @try Base.UUID(pkguuid)

    main = get(input, "main", nothing)
    main isa AbstractString || return ArgumentError("Keyword argument `main` of type `String` is required.")

    pkgid = Base.PkgId(pkguuid, pkgname)
    pkg = @try Base.require(pkgid)
    return @try mapfoldl(Symbol, getproperty, split(main, "."); init = pkg)
end

function callmain_api(input)
    @require_keyword_arguments input
    @debug "callmain_api" input

    commonopts = @return_exception getcommonoptions(input)
    f = @return_exception getcallable(input)
    args = @return_exception getstrings(input, "args")
    return apicall(commonopts) do
        Base.invokelatest(f, args)
    end
end

function callany_api(input)
    @require_keyword_arguments input
    @debug "callany_api" input

    commonopts = @return_exception getcommonoptions(input)

    args = get(input, "args", [])
    args isa AbstractVector || return ArgumentError("Keyword argument `args` must be of type `Vector`")
    args = refine_values(args)

    kwargs = get(input, "kwargs", Dict())
    kwargs isa AbstractDict || return ArgumentError("Keyword argument `kwargs` must be of type `Dict`")
    kwargs = Dict(Symbol(k) => v for (k, v) in refine_values(kwargs))

    f = @return_exception getcallable(input)
    return apicall(commonopts) do
        Base.invokelatest(f, args...; kwargs...)
    end
end

function help_api(input)
    @require_keyword_arguments input
    @debug "help_api" input
    mime = get(input, "mime", "text/plain")
    commonopts = @return_exception getcommonoptions(input)
    f = @return_exception getcallable(input)

    # Print to the terminal directly if appropriate:
    if mime === "text/plain" && commonopts.ioconfig.stdout !== nothing
        return apicall(commonopts) do
            show(stdout, mime, Docs.doc(f))
            println(stdout)
            return (result = nothing, output = nothing)
        end
    end

    docstr = sprint(show, mime, Docs.doc(f)) * "\n"
    return (result = nothing, output = docstr)
end

dispatcher = Dict(
    # JSONRPC method name => implementation
    "eval" => eval_api,
    "callmain" => callmain_api,
    "callany" => callany_api,
    "help" => help_api,
)

function project_toml_from_manifest_toml(
    projectpath::AbstractString,
    manifestpath::AbstractString,
)
    open(projectpath, write = true) do io
        project_toml_from_manifest_toml(io, manifestpath)
    end
end

function project_toml_from_manifest_toml(io::IO, manifestpath::AbstractString)
    manifest = TOML.parsefile(manifestpath)
    deps = map(keys(manifest), values(manifest)) do k, v
        k => v[1]["uuid"]
    end

    println(io, "[deps]")
    for (pkg, uuid) in deps
        println(io, pkg, " = ", repr(uuid))
    end
end

function as_proper_project(f, project)
    if startswith(project, "@") ||
       isdir(project) || basename(project) in ("Project.toml", "JuliaProject.toml")
        return f(Base.load_path_expand(project))
    end

    mktempdir() do dir
        cp(project, joinpath(dir, "Manifest.toml"))
        tomlpath = joinpath(dir, "Project.toml")
        project_toml_from_manifest_toml(tomlpath, project)
        @debug("Project.toml created", Project_toml = Text(read(tomlpath, String)))

        f(tomlpath)
    end
end

let periodstr = get(ENV, "JULIA_CLI_PERIODIC_GC", "no")
    if all(isdigit, periodstr)
        period = parse(Int, periodstr)
        @debug "Calling `GC.gc()` every $period seconds..."
        @async while true
            sleep(period)
            GC.gc()
        end
    end
end

mktempdir(prefix = "juliacli_") do dir
    path = joinpath(dir, "socket")
    @debug "Launching worker server at: $path"
    server = listen(path)
    println(path)
    close(stdout)
    redirect_stdout(stderr)  # required for `Pkg.activate`?
    try
        as_proper_project(project) do project
            Pkg.activate(project)
            JSONRPC.serve(dispatcher, server, JSONRPC.NDJSON)
        end
    catch err
        @error "JSONRPC.serve" exception = (err, catch_backtrace())
    finally
        @debug "Stopped worker server at: $path"
    end
end

end


text:
module PPSIM

using YAML
using DataFrames
using XPP

include("src/specifications.jl")
include("src/streams.jl")
include("src/models.jl")
include("src/initialise.jl")
include("src/query.jl")

end

max_stars_repo_name: el-uhu/PPSIM

max_issues_repo_name: el-uhu/PPSIM

max_forks_repo_name: el-uhu/PPSIM


text:
function equals(obj::FieldGradientField, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_field(::Type{FieldGradientField}, arg0::Field, arg1::jint)
    return jcall(FieldGradientField, "getField", FieldGradientField, (Field, jint), arg0, arg1)
end

function get_one(obj::FieldGradientField)
    return jcall(obj, "getOne", FieldGradient, ())
end

function get_runtime_class(obj::FieldGradientField)
    return jcall(obj, "getRuntimeClass", Class, ())
end

function get_zero(obj::FieldGradientField)
    return jcall(obj, "getZero", FieldGradient, ())
end

function hash_code(obj::FieldGradientField)
    return jcall(obj, "hashCode", jint, ())
end


text:
module TestUserApi

using ..MyPetStore
using Swagger
using Test
using Random
using JSON
using URIs

const TEST_USER = "jlswag"
const TEST_USER1 = "jlswag1"
const TEST_USER2 = "jlswag2"

function test_404(uri)
    @info("Error handling")
    client = Swagger.Client(uri*"_invalid")
    api = UserApi(client)

    try
        loginUser(api, TEST_USER, "testpassword")
        @error("Swagger.ApiException not thrown")
    catch ex
        @test isa(ex, Swagger.ApiException)
        @test ex.status == 404
    end

    client = Swagger.Client("http://_invalid/")
    api = UserApi(client)

    try
        loginUser(api, TEST_USER, "testpassword")
        @error("Swagger.ApiException not thrown")
    catch ex
        @test isa(ex, Swagger.ApiException)
        @test startswith(ex.reason, "Could not resolve host")
    end
end

function test_set_methods()
    @info("Error handling")
    client = Swagger.Client("http://_invalid/")

    @test client.timeout[] == Swagger.DEFAULT_TIMEOUT_SECS

    Swagger.with_timeout(client, Swagger.DEFAULT_TIMEOUT_SECS + 10) do client
        @test client.timeout[] == Swagger.DEFAULT_TIMEOUT_SECS + 10
    end
    @test client.timeout[] == Swagger.DEFAULT_TIMEOUT_SECS

    api = UserApi(client)
    Swagger.with_timeout(api, Swagger.DEFAULT_TIMEOUT_SECS + 10) do api
        @test api.client.timeout[] == Swagger.DEFAULT_TIMEOUT_SECS + 10
    end
    @test client.timeout[] == Swagger.DEFAULT_TIMEOUT_SECS

    Swagger.set_timeout(client, Swagger.DEFAULT_TIMEOUT_SECS + 10)
    @test client.timeout[] == Swagger.DEFAULT_TIMEOUT_SECS + 10

    @test isempty(client.headers)
    Swagger.set_user_agent(client, "007")
    Swagger.set_cookie(client, "crumbly")
    @test client.headers["User-Agent"] == "007"
    @test client.headers["Cookie"] == "crumbly"
end

function test_login_user_hook(ctx::Swagger.Ctx)
    ctx.header["actual_password"] = "testpassword"
    ctx
end

function test_login_user_hook(resource_path::AbstractString, body::Any, headers::Dict{String,String})
    uri = URIs.parse_uri(resource_path)
    qparams = URIs.queryparams(uri)
    qparams["password"] = headers["actual_password"]
    delete!(headers, "actual_password")
    resource_path = string(URIs.URI(uri; query=escapeuri(qparams)))

    (resource_path, body, headers)
end

function test_userhook(uri)
    @info("User hook")
    client = Swagger.Client(uri; pre_request_hook=test_login_user_hook)
    api = UserApi(client)

    login_result = loginUser(api, TEST_USER, "wrongpassword")
    @test !isempty(login_result)
    result = JSON.parse(login_result)
    @test startswith(result["message"], "logged in user session")
    @test result["code"] == 200
end

function test_parallel(uri)
    @info("Parallel usage")
    client = Swagger.Client(uri)
    api = UserApi(client)

    for gcidx in 1:100
        @sync begin
            for idx in 1:10^3
                @async begin
                    @debug("[$idx] UserApi Parallel begin")
                    login_result = loginUser(api, TEST_USER, "testpassword")
                    @test !isempty(login_result)
                    result = JSON.parse(login_result)
                    @test startswith(result["message"], "logged in user session")
                    @test result["code"] == 200

                    @test_throws Swagger.ApiException getUserByName(api, randstring())
                    @test_throws Swagger.ApiException getUserByName(api, TEST_USER)

                    logout_result = logoutUser(api)
                    @test logout_result === nothing
                    @debug("[$idx] UserApi Parallel end")
                end
            end
        end
        GC.gc()
        @info("outer loop $gcidx")
    end
    nothing
end

function test(uri)
    @info("UserApi")
    client = Swagger.Client(uri)
    api = UserApi(client)

    @info("UserApi - loginUser")
    login_result = loginUser(api, TEST_USER, "testpassword")
    @test !isempty(login_result)

    #@info("UserApi - createUser")
    #user1 = User(; id=100, username=TEST_USER1, firstName="test1", lastName="user1", email="jlswag1@example.com", password="testpass1", phone="1000000001", userStatus=0)
    #@test createUser(api, user1) === nothing

    #@info("UserApi - createUsersWithArrayInput")
    #user2 = User(; id=200, username=TEST_USER2, firstName="test2", lastName="user2", email="jlswag2@example.com", password="testpass2", phone="1000000002", userStatus=0)
    #@test createUsersWithArrayInput(api, [user1, user2]) === nothing

    #@info("UserApi - createUsersWithListInput")
    #@test createUsersWithListInput(api, [user1, user2]) === nothing

    @info("UserApi - getUserByName")
    @test_throws Swagger.ApiException getUserByName(api, randstring())
    @test_throws Swagger.ApiException getUserByName(api, TEST_USER)
    #getuser_result = getUserByName(api, TEST_USER)
    #@test isa(getuser_result, User)

    #@info("UserApi - updateUser")
    #@test updateUser(api, TEST_USER2, getuser_result) === nothing
    #@info("UserApi - deleteUser")
    #@test deleteUser(api, TEST_USER2) === nothing

    @info("UserApi - logoutUser")
    logout_result = logoutUser(api)
    @test logout_result === nothing

    nothing
end

end # module TestUserApi


text:
const re_format = Dict([('a', "[A-Z]"), ('c', "[A-Z0-9]"), ('n', "[0-9]")])

const iban_global_re = r"^([A-Z]{2})([0-9]{2})([A-Z0-9]{10,34})$"

const iban_country_length = Dict([
                            ("AL", 28),  # Albania
                            ("AD", 24),  # Andorra
                            ("AE", 23),  # United Arab Emirates
                            ("AT", 20),  # Austria
                            ("AZ", 28),  # Azerbaijan
                            ("BA", 20),  # Bosnia and Herzegovina
                            ("BE", 16),  # Belgium
                            ("BG", 22),  # Bulgaria
                            ("BH", 22),  # Bahrain
                            ("BR", 29),  # Brazil
                            ("CH", 21),  # Switzerland
                            ("CR", 21),  # Costa Rica
                            ("CY", 28),  # Cyprus
                            ("CZ", 24),  # Czech Republic
                            ("DE", 22),  # Germany
                            ("DK", 18),  # Denmark
                            ("DO", 28),  # Dominican Republic
                            ("EE", 20),  # Estonia
                            ("ES", 24),  # Spain
                            ("FI", 18),  # Finland
                            ("FO", 18),  # Faroe Islands
                            ("FR", 27),  # France
                            ("GB", 22),  # United Kingdom
                            ("GE", 22),  # Georgia
                            ("GI", 23),  # Gibraltar
                            ("GL", 18),  # Greenland
                            ("GR", 27),  # Greece
                            ("GT", 28),  # Guatemala
                            ("HR", 21),  # Croatia
                            ("HU", 28),  # Hungary
                            ("IE", 22),  # Ireland
                            ("IL", 23),  # Israel
                            ("IS", 26),  # Iceland
                            ("IT", 27),  # Italy
                            ("JO", 30),  # Jordan
                            ("KZ", 20),  # Kazakhstan
                            ("KW", 30),  # Kuwait
                            ("LB", 28),  # Lebanon
                            ("LI", 21),  # Liechtenstein
                            ("LT", 20),  # Lithuania
                            ("LU", 20),  # Luxembourg
                            ("LV", 21),  # Latvia
                            ("MC", 27),  # Monaco
                            ("MD", 24),  # Moldova
                            ("ME", 22),  # Montenegro
                            ("MK", 19),  # Macedonia
                            ("MT", 31),  # Malta
                            ("MR", 27),  # Mauritania
                            ("MU", 30),  # Mauritius
                            ("NL", 18),  # Netherlands
                            ("NO", 15),  # Norway
                            ("PS", 29),  # Palestine
                            ("PK", 24),  # Pakistan
                            ("PL", 28),  # Poland
                            ("PT", 25),  # Portugal
                            ("QA", 29),  # Qatar
                            ("RO", 24),  # Romania
                            ("RS", 22),  # Serbia
                            ("SA", 24),  # Saudi Arabia
                            ("SE", 24),  # Sweden
                            ("SI", 19),  # Slovenia
                            ("SK", 24),  # Slovakia
                            ("SM", 27),  # San Marino
                            ("TN", 24),  # Tunisia
                            ("TR", 26),  # Turkey
                            ("VG", 24)  # British Virgin Islands
                             ])

const bban_format = Dict([
                          ("GB", ((4, 'a'), (14, 'n'))),
                          ("DE", ((18, 'n'), ))
                          ])


text:
using JulzDummy
using Base.Test
using TestSetExtensions

@testset ExtendedTestSet "All Tests" begin
  if ( endswith(pwd(), "/test") ) ; cd("..") ; end

  seed_int = abs(rand(Int16))
  srand(seed_int)
  println("\n Seed: $seed_int \n")

  test_dir = "$(pwd())/test"

  is_focused = Julz.check_for_focus(test_dir)

  is_sorted = Julz.tests_are_sorted

  Julz.include_all_files( test_dir,
    is_testing=true, is_focused=is_focused, is_sorted=is_sorted,
    reload_function=JulzDummy.load_input
  )

  @test true
end

return


text:
# Compare the times of simulation
using NbodyGradient, DelimitedFiles, BenchmarkTools

elements = readdlm("elements_rebound.txt", ',')
elements[:,end] .= 0.0

try
  rm("nbg_times.txt")
catch
end

for n = 2:11
  t0 = 0.0
  h = 0.05
  # Rebound integration is for 5000 code units; inner planet has
  # a period of 2π, while in our case 1-day, so differs by 5000/2π
  tmax = 2500/π

  ic = ElementsIC(t0, n, elements[1:n,:])
  s = State(ic)
  intr = Integrator(h, 0.0, tmax)

  # Time without derivatives
  t = @belapsed intr(s, grad=false) setup=(s = State($ic))

  # Time with derivatives
  grad_t = @belapsed intr(s) setup=(s = State($ic))

  #println("Nplanet: ", n-1, " No gradient: ", t, " gradient: ", grad_t, " ratio: ", grad_t/t)
  println(n-1, ",", t,",",grad_t,",", grad_t/t)
  open("nbg_times.txt", "a") do io
    writedlm(io, [n-1 t grad_t grad_t/t], ',')
  end
end

max_forks_repo_name: langfzac/nbg-papers


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/MortarContact2DAD.jl/blob/master/LICENSE

using Documenter, MortarContact2DAD

format = length(ARGS) > 0 ? Symbol(ARGS[1]) : :html

makedocs(modules=[MortarContact2DAD],
         format = format,
         checkdocs = :all,
         sitename = "MortarContact2DAD.jl",
         pages = ["index.md"]
        )


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

using Markdown
using InteractiveUtils

# ╔═╡ 79eb61d4-9635-11ec-19d4-7b61646ffe23
# hideall
begin
    using Pkg
    Pkg.activate(; temp=true)
    Pkg.add([
        "GLM"
    ])
    Pkg.develop(; path=dirname(dirname(dirname(@__DIR__))))
end

# ╔═╡ f2855237-a2dd-4c11-ad8e-2fa58a0eae9d
# hideall
try using Revise catch end

# ╔═╡ 84652bbd-f82a-4af3-ab1d-89b9a79da966
using GLM

# ╔═╡ 04507b23-fb11-43bb-aea0-02f9eddd75c5
using GLM: AbstractGLM

# ╔═╡ 4d83954a-0d01-46d2-9a80-09212bbd8c33
md"""
# GLM

This page demonstrates an example implementation of the `JModels.jl` interface for `GLM.jl` and applies `GLM.jl` via the interface.
"""

# ╔═╡ 75c323c8-9ff5-48e6-923b-0eb02f7ef9df
md"""
## Implementation of JModels

This is an example implementation which could live inside `GLM.jl`.
"""

# ╔═╡ 1e6c881f-ae0a-4393-ac42-cd212428f20b
import JModels

# ╔═╡ 78d45aff-63ff-468d-8ca5-b70669dce442
JModels.ismodel(AbstractGLM) = true

# ╔═╡ a3dff183-54ae-4215-9f0c-0456bbe9f30b
function JModels.fit(
        mt::Type{GeneralizedLinearModel},
        data;
        settings = (;
			formula=FormulaTerm(term(:y), term(:x)),
        	family=Binomial(),
	        link=LogitLink()
		)
    )
    return glm(settings.formula, data, settings.family, settings.link)
end

# ╔═╡ 369f5d90-f796-44c7-9e6c-f7396078053a
function JModels.predict(
        model::Union{AbstractGLM,LinearModel,RegressionModel},
        data;
        settings=NamedTuple()
    )
    return GLM.predict(model, data; settings...)
end

# ╔═╡ 2868f696-a892-47e9-b6c9-9bdf53aeed39
md"""
## Example usage
"""

# ╔═╡ 318832b4-a04c-4709-9783-442da4f4fc9f
data = (; x = 1:3, y = rand(3));

# ╔═╡ 1e2e759f-8a1f-43a4-98c7-81dee72f5160
model = JModels.fit(GeneralizedLinearModel, data)

# ╔═╡ 14ae9587-e62b-479b-a32b-5686272df916
JModels.predict(model, data)

# ╔═╡ Cell order:
# ╠═4d83954a-0d01-46d2-9a80-09212bbd8c33
# ╠═79eb61d4-9635-11ec-19d4-7b61646ffe23
# ╠═f2855237-a2dd-4c11-ad8e-2fa58a0eae9d
# ╠═84652bbd-f82a-4af3-ab1d-89b9a79da966
# ╠═04507b23-fb11-43bb-aea0-02f9eddd75c5
# ╠═75c323c8-9ff5-48e6-923b-0eb02f7ef9df
# ╠═1e6c881f-ae0a-4393-ac42-cd212428f20b
# ╠═78d45aff-63ff-468d-8ca5-b70669dce442
# ╠═a3dff183-54ae-4215-9f0c-0456bbe9f30b
# ╠═369f5d90-f796-44c7-9e6c-f7396078053a
# ╠═2868f696-a892-47e9-b6c9-9bdf53aeed39
# ╠═318832b4-a04c-4709-9783-442da4f4fc9f
# ╠═1e2e759f-8a1f-43a4-98c7-81dee72f5160
# ╠═14ae9587-e62b-479b-a32b-5686272df916


text:
## `GroupElement` interface

### Obligatory methods for `GroupElement`

Base.parent(g::GroupElement) =
    throw(InterfaceNotImplemented(:Group, "Base.parent(::$(typeof(g)))"))

@doc Markdown.doc"""
    parent_type(element_type)
Given the type of an element return the type of its parent.
"""
AbstractAlgebra.parent_type(GEl::Type{<:GroupElement}) = throw(
    InterfaceNotImplemented(:Group, "GroupsCore.parent_type(::Type{$GEl})"),
)

@doc Markdown.doc"""
    istrulyequal(g::GEl, h::GEl) where {GEl<:GroupElement}
Return the mathematical equality of group elements.

This function may not return due to e.g. unsolvable word problem in groups.
"""
istrulyequal(g::GEl, h::GEl) where {GEl<:GroupElement} = throw(
    InterfaceNotImplemented(:Group, "GroupsCore.istrulyequal(::$GEl, ::$GEl)"),
)

@doc Markdown.doc"""
    hasorder(g::GroupElement)
Return `true` if `g` has finite order (without computing it).
"""
hasorder(g::GroupElement) = throw(
    InterfaceNotImplemented(:Group, "GroupsCore.hasorder(::$(typeof(g)))"),
)

@doc Markdown.doc"""
    deepcopy_internal(g::GroupElement, ::IdDict)
Return a completely intependent copy of group element `g` **without copying its parent**.

That is `parent(g) === parent(deepcopy(g))` must be satisfied. There is no need to implement this method if `g` is `isbits`.
"""
Base.deepcopy_internal(g::GroupElement, stackdict::IdDict) = throw(
    InterfaceNotImplemented(
        :Group,
        "Base.deepcopy_internal(::$(typeof(g)), ::IdDict)",
    ),
)
# TODO: Technically, it is not necessary to implement `deepcopy_internal` method if `parent(g)` can be reconstructed exactly from `g` (i.e. either it's cached, or a singleton). However by defining this fallback we force everybody to implement it, except isbits group elements.
@doc Markdown.doc"""
    inv(g::GroupElement)
Return the group inverse of `g`
"""
Base.inv(g::GroupElement) =
    throw(InterfaceNotImplemented(:Group, "Base.inv(::$(typeof(g)))"))
@doc Markdown.doc"""
    (*)(g::GEl, h::GEl) where GEl<:GroupElement
Return the result of group binary operation on `g` and `h`.
"""
Base.:(*)(g::GEl, h::GEl) where {GEl<:GroupElement} = throw(
    InterfaceNotImplemented(
        :Group,
        "Base.:(*)(::$(typeof(g)), ::$(typeof(g)))",
    ),
)

### Default implementations for `GroupElement`

#### Modification not recommended
@doc Markdown.doc"""
    one(g::GroupElement)
Return the identity element of the parent group of `g`.
"""
Base.one(g::GroupElement) = one(parent(g))

@doc Markdown.doc"""
    order([BigInt, ]g::GroupElement)
    order(I::Type{<:Integer}, g::GroupElement)
Return the order of `g` as an instance of `I`.

Only arbitrary size integers are required to return mathematically correct answer.
"""
AbstractAlgebra.order(g::GroupElement) = order(BigInt, g)

@doc Markdown.doc"""
    conj(g::GEl, h::GEl) where {GEl<:GroupElement}
Return conjugation of `g` by `h`, i.e. `h^-1*g*h`.
"""
Base.conj(g::GEl, h::GEl) where {GEl<:GroupElement} = conj!(similar(g), g, h)

@doc Markdown.doc"""
    ^(g::GEl, h::GEl) where {GEl<:GroupElement}
Conjugation action of `h` on `g`. See `conj`.
"""
Base.:(^)(g::GEl, h::GEl) where {GEl<:GroupElement} = conj(g, h)

@doc Markdown.doc"""
    comm(g::GEl, h::GEl[, Vararg{GEl}...) where GEl<:GroupElement
Return the commutator `inv(g)*inv(h)*g*h` of `g` and `h`.

The `Vararg` version returns the repeated (`foldl`) commutator, i.e.
`comm(g, h, k) == comm(comm(g, h), k)`.
"""
function comm(g::GEl, h::GEl, tail::GEl...) where {GEl<:GroupElement}
    res = comm!(similar(g), g, h)
    for k in tail
        res = comm!(res, res, k)
    end
    return res
end

Base.literal_pow(::typeof(^), g::GroupElement, ::Val{-1}) = inv(g)

Base.:(/)(g::GEl, h::GEl) where {GEl<:GroupElement} =
    div_right!(similar(g), g, h)

#### Modification possible for performance reasons

@doc Markdown.doc"""
    similar(g::GroupElement)
Return an arbitrary (and possibly uninitialized) group element sharing the parent with `g`.
"""
Base.similar(g::GroupElement) = one(g)
# optimization: determine triviality of `g` without constructing `one(g)`.
Base.isone(g::GroupElement) = g == one(g)

# TODO: semantic clash: isequal is weaker in julia than `==`, we need it the other way round here → istrulyequal
@doc Markdown.doc"""
    ==(g::GEl, h::GEl) where {GEl<:GroupElement}
The "best effort" equality for group elements.

Depending on the group this might, or might not be the correct equality, which can be obtained using `Base.isequal`. Nonetheless the implication `g == h` → `istrulyequal(g, h)` must be always satisfied, i.e. "best effort" equality might return `false` even when group elements are equal.

For example in a finitely presented group, `==` may return the equality of words.
"""
Base.:(==)(g::GEl, h::GEl) where {GEl<:GroupElement} = istrulyequal(g, h)

function Base.:^(g::GroupElement, n::Integer)
    n == 0 && return one(g)
    n < 0 && return inv(g)^-n
    return Base.power_by_squaring(g, n)
end

#### Modification RECOMMENDED for performance reasons

function AbstractAlgebra.order(::Type{I}, g::GroupElement) where {I<:Integer}
    hasorder(g) || throw("$g does not seem to have finite order")
    isone(g) && return I(1)
    o = I(1)
    gg = deepcopy(g)
    out = similar(g)
    while !isone(gg)
        o += I(1)
        gg = mul!(out, gg, g)
    end
    return o
end

Base.hash(g::GroupElement, h::UInt) = hash(typeof(g), h)

##### Mutable API
@doc Markdown.doc"""
    one!(g::GroupElement)
Return `one(g)`, possibly modifying `g`.
"""
one!(g::GroupElement) = one(parent(g))

@doc Markdown.doc"""
    inv!(out::GEl, g::GEl) where GEl<:GroupElement
Return `inv(g)`, possibly modifying `out`.

Aliasing of `g` with `out` is allowed.
"""
AbstractAlgebra.inv!(out::GEl, g::GEl) where {GEl<:GroupElement} = inv(g)

@doc Markdown.doc"""
    mul!(out::GEl, g::GEl, h::GEl) where GEl<:GroupElement
Return `g*h`, possibly modifying `out`.

Aliasing of `g` or `h` with `out` is allowed.
"""
AbstractAlgebra.mul!(out::GEl, g::GEl, h::GEl) where {GEl<:GroupElement} = g * h

@doc Markdown.doc"""
    div_right!(out::GEl, g::GEl, h::GEl) where GEl<:GroupElement
Return `g*h^-1`, possibly modifying `out`.

Aliasing of `g` or `h` with `out` is allowed;
"""
div_right!(out::GEl, g::GEl, h::GEl) where {GEl<:GroupElement} =
    mul!(out, g, inv(h))

@doc Markdown.doc"""
    div_left!(out::GEl, g::GEl, h::GEl) where GEl<:GroupElement
Return `h^-1*g`, possibly modifying `out`.

Aliasing of `g` or `h` with `out` is allowed;
"""
function div_left!(out::GEl, g::GEl, h::GEl) where {GEl<:GroupElement}
    out = (out === g || out === h) ? inv(h) : inv!(out, h)
    return mul!(out, out, g)
end
@doc Markdown.doc"""
    conj!(out::GEl, g::GEl, h::GEl) where GEl<:GroupElement
Return `h^-1*g*h, `possibly modifying `out`.

Aliasing of `g` or `h` with `out` is allowed.
"""
function conj!(out::GEl, g::GEl, h::GEl) where {GEl<:GroupElement}
    out = (out === g || out === h) ? inv(h) : inv!(out, h)
    out = mul!(out, out, g)
    return mul!(out, out, h)
end

@doc Markdown.doc"""
    comm!(out::GEl, g::GEl, h::GEl) where GEl<:GroupElement
Return `g^-1*h^-1*g*h, `possibly modifying `out`.

Aliasing of `g` or `h` with `out` is allowed.
"""
# TODO: can we make comm! with 3 arguments without allocation??
function comm!(out::GEl, g::GEl, h::GEl) where {GEl<:GroupElement}
    out = conj!(out, g, h)
    return div_left!(out, out, g)
end


text:
module CompTime
export @ct_enable, runtime, comptime, generate_code, debug

"""
Notes:

Currently does not work with kwargs. These are tricky to handle correctly, so as a proof of concept, we do not handle them.
"""

using MLStyle
using MacroTools: splitdef, combinedef, @capture, postwalk, prewalk

function runtime end
function comptime end
function generate_code end

function debug(fn, args...)
  generate_code(fn, map(typeof, args)...)
end

macro ct_enable(def::Expr)
  @assert def.head == :function
  parts = FunctionParts(def)
  esc(quote
    $(make_apply_def(parts))

    $(make_generate_code_def(parts))

    $(make_comptime_def(parts))

    $(make_runtime_def(parts))
  end)
end

struct FunctionParts
  name::Any
  args::Vector{Expr}
  whereparams::Any
  body::Any
  function FunctionParts(def::Expr)
    parts = splitdef(def)
    @assert parts[:kwargs] == []
    new(parts[:name], normalize_args(parts[:args]), parts[:whereparams], parts[:body])
  end
  function FunctionParts(name, args, whereparams, body)
    new(name, args, whereparams, body)
  end
end

function arg_names(parts::FunctionParts)
  map(parts.args) do arg
    @match arg begin
      Expr(:(::), name, _) => name
    end
  end
end

function arg_types(parts::FunctionParts)
  map(parts.args) do arg
    @match arg begin
      Expr(:(::), _, type) => type
    end
  end
end

function make_function(parts::FunctionParts)
  combinedef(Dict(
    :name=>parts.name,
    :args=>parts.args,
    :whereparams=>parts.whereparams,
    :body=>parts.body,
    :kwargs=>[]
  ))
end

function normalize_args(args::Vector)
  map(args) do arg
    @match arg begin
      Expr(:(::), name, type) => Expr(:(::), name, type)
      Expr(:(::), type) => Expr(:(::), gensym(), type)
      name::Symbol => Expr(:(::), name, Any)
      _ => error("unsupported argument format $arg")
    end
  end
end

function make_runtime_def(parts::FunctionParts)
  make_function(FunctionParts(
    GlobalRef(CompTime, :runtime),
    [Expr(:(::), Expr(:call, typeof, parts.name)), parts.args...],
    parts.whereparams,
    postwalk(strip_ct_macros, parts.body)
  ))
end

function make_comptime_def(parts::FunctionParts)
  make_function(FunctionParts(
    GlobalRef(CompTime, :comptime),
    [Expr(:(::), Expr(:call, typeof, parts.name)), parts.args...],
    parts.whereparams,
    quote
      if $(Expr(:generated))
        $(generate_code)($(parts.name), $(arg_types(parts)...))
      else
        $(Expr(:meta, :generated_only))
        return
      end
    end
  ))
end

function make_generate_code_def(parts::FunctionParts)
  type_args = map(arg_types(parts)) do type
    Expr(:(::), Expr(:curly, Type, type))
  end

  make_function(FunctionParts(
    GlobalRef(CompTime, :generate_code),
    [Expr(:(::), Expr(:call, typeof, parts.name)), type_args...],
    parts.whereparams,
    comptime_expr(parts.body)
  ))
end

function make_apply_def(parts::FunctionParts)
  make_function(FunctionParts(
    parts.name,
    parts.args,
    parts.whereparams,
    Expr(:call, comptime, parts.name, arg_names(parts)...)
  ))
end

function comptime_loop(mkloop, body)
  lines = gensym("lines")
  mkbody = comptime_expr(body)
  quote
    $lines = []
    $(mkloop(:(push!($lines, $mkbody))))
    Expr(:block, $lines...)
  end
end

function comptime_if(parts)
  branches = normalize_if(parts...)
  (cond, thenpart) = branches[1]
  rest = branches[2:end]
  Expr(:if, cond, comptime_expr(thenpart), comptime_else(rest))
end

function comptime_else(branches)
  if length(branches) > 0
    (cond, thenpart) = branches[1]
    Expr(:elseif, cond, comptime_expr(thenpart), comptime_else(branches[2:end]))
  else
    nothing
  end
end

function normalize_if(cond, thenpart)
  [(cond, thenpart)]
end

function normalize_if(cond, thenpart, elsepart)
  [(cond, thenpart), normalize_else(elsepart)...]
end

function normalize_else(elsepart)
  @match elsepart begin
    Expr(:elseif, cond, thenpart) =>
      [(cond, thenpart)]
    Expr(:elseif, cond, thenpart, moreelse) =>
      [(cond, thenpart), normalize_else(moreelse)...]
    _ => [(true, elsepart)]
  end
end

q(s::Symbol) = Expr(:quote, s)
q(e::Expr) = Expr(:quote, e)
q(x) = x

function ct_control(expr)
  generator = @match expr begin
    Expr(:for, head, body) =>
      comptime_loop(gen -> Expr(:for, head, gen), body)
    Expr(:if, parts...) =>
      comptime_if(parts)
    Expr(:while, cond, body) =>
      comptime_loop(gen -> Expr(:while, cond, gen), body)
    _ => error("unsupported control structure $expr")
  end
  Expr(:$, generator)
end

function expand_ct_macros(expr)
  @match expr begin
    Expr(:macrocall, mname, _, e) =>
      if mname == Symbol("@ct")
        Expr(:$, Expr(:call, q, e))
      elseif mname == Symbol("@ct_ctrl")
        ct_control(e)
      else
        expr
      end
    _ => expr
  end
end

function strip_ct_macros(expr)
  @match expr begin
    Expr(:macrocall, mname, _, e) =>
      if mname == Symbol("@ct")
        e
      elseif mname == Symbol("@ct_ctrl")
        e
      else
        expr
      end
    _ => expr
  end
end

function comptime_expr(expr)
  Expr(:quote, postwalk(expand_ct_macros, expr))
end

end


text:
using FourierCore, FourierCore.Grid, FourierCore.Domain
using FFTW, LinearAlgebra, BenchmarkTools, Random, JLD2
using GLMakie
rng = MersenneTwister(1234)
# Random.seed!(123456789)
Random.seed!(12)
# jld_name = "high_order_timestep_spatial_tracer_"
jld_name = "blocky"
include("transform.jl")
include("random_phase_kernel.jl")
# using GLMakie
using CUDA
arraytype = CuArray
Ω = S¹(4π)^2
N = 2^7 # number of gridpoints
grid = FourierGrid(N, Ω, arraytype=arraytype)
nodes, wavenumbers = grid.nodes, grid.wavenumbers

x = nodes[1]
y = nodes[2]
kˣ = wavenumbers[1]
kʸ = wavenumbers[2]
# construct filter
kxmax = maximum(kˣ)
kymax = maximum(kˣ)
filter = @. (kˣ)^2 + (kʸ)^2 ≤ ((kxmax / 2)^2 + (kymax / 2)^2)

# now define the random field
wavemax = 3
𝓀 = arraytype(collect(-wavemax:0.5:wavemax))
𝓀ˣ = reshape(𝓀, (length(𝓀), 1))
𝓀ʸ = reshape(𝓀, (1, length(𝓀)))
A = @. 0.1 * (𝓀ˣ * 𝓀ˣ + 𝓀ʸ * 𝓀ʸ)^(-11 / 12)
A[A.==Inf] .= 0.0
φ = arraytype(2π * rand(size(A)...))
field = arraytype(zeros(N, N))

##
# Fields
# velocity
ψ = arraytype(zeros(ComplexF64, N, N))
u = similar(ψ)
v = similar(ψ)

# theta
θ = similar(ψ)
∂ˣθ = similar(ψ)
∂ʸθ = similar(ψ)
κΔθ = similar(ψ)
θ̇ = similar(ψ)
s = similar(ψ)
θ̅ = similar(ψ)
k₁ = similar(ψ)
k₂ = similar(ψ)
k₃ = similar(ψ)
k₄ = similar(ψ)
θ̃ = similar(ψ)

# source
s = similar(ψ)
@. s = cos(kˣ[5] * x)

# phase
φ̇ = similar(A)

# operators
∂x = im * kˣ
∂y = im * kʸ
Δ = @. ∂x^2 + ∂y^2

# plan ffts
P = plan_fft!(ψ)
P⁻¹ = plan_ifft!(ψ)

# number of gridpoints in transition is about λ * N / 2
bump(x; λ=10 / N, width = π/2) = 0.5 * (tanh((x + width / 2) / λ) - tanh((x - width / 2) / λ))
bumps(x; λ=20 / N, width = 1.0) = 0.25 * (bump(x, λ = λ, width = width) + bump(x, λ = λ, width = 2.0*width) + bump(x, λ = λ, width = 3.0*width) + bump(x, λ = λ, width = 4.0*width))

##
Δx = x[2] - x[1]
κ = Δx^2
# κ = 2 / 2^8 # fixed diffusivity
# κ = 2e-4
Δt = (Δx) / (4π) * 5

# take the initial condition as negative of the source
tic = Base.time()

# save some snapshots
ψ_save = typeof(real.(Array(ψ)))[]
θ_save = typeof(real.(Array(ψ)))[]

r_A = Array(@. sqrt((x - 2π)^2 + (y - 2π)^2))
θ_A = [bumps(r_A[i, j]) for i in 1:N, j in 1:N]
θ .= CuArray(θ_A)
# @. θ = bump(sqrt(x^2 + y^2)) # scaling so that source is order 1
θclims = extrema(Array(real.(θ))[:])
P * θ # in place fft
@. κΔθ = κ * Δ * θ
P⁻¹ * κΔθ # in place fft
s .= -κΔθ * 0.0
P⁻¹ * θ # in place fft
θ̅ .= 0.0

t = [0.0]
tend = 5

iend = ceil(Int, tend / Δt)

params = (; ψ, A, 𝓀ˣ, 𝓀ʸ, x, y, φ, u, v, ∂ˣθ, ∂ʸθ, s, P, P⁻¹, filter)

size_of_A = size(A)

realizations = 100
for j in 1:realizations
    # new realization of flow
    rand!(rng, φ) # between 0, 1
    φ .*= 2π # to make it a random phase
    event = stream_function!(ψ, A, 𝓀ˣ, 𝓀ʸ, x, y, φ)
    wait(event)
    t[1] = 0

    θ .= CuArray(θ_A)
    for i = 1:iend
        # fourth order runge kutta on deterministic part
        # keep ψ frozen is the correct way to do it here

        # the below assumes that φ is just a function of time
        θ_rhs_new!(k₁, θ, params)
        @. θ̃ = θ + Δt * k₁ * 0.5

        φ_rhs_normal!(φ̇, φ, rng)
        @. φ += sqrt(Δt / 2) * φ̇

        θ_rhs_new!(k₂, θ̃, params)
        @. θ̃ = θ + Δt * k₂ * 0.5
        θ_rhs_new!(k₃, θ̃, params)
        @. θ̃ = θ + Δt * k₃

        φ_rhs_normal!(φ̇, φ, rng)
        @. φ += sqrt(Δt / 2) * φ̇

        θ_rhs_new!(k₄, θ̃, params)
        @. θ += Δt / 6 * (k₁ + 2 * k₂ + 2 * k₃ + k₄)

        # update stochastic part
        # φ_rhs_normal!(φ̇, φ, rng)
        # @. φ += sqrt(Δt) * φ̇

        t[1] += Δt
        # save output

        #=
        if i % div(iend, 10) == 0
            println("Saving at i=", i)
            push!(ψ_save, Array(real.(ψ)))
            push!(θ_save, Array(real.(θ)))
            println("extrema are ", extrema(θ_save[end]))
            println("time is t = ", t[1])
        end

        if t[1] >= tstart
            θ̅ .+= Δt * θ
        end

        if i % div(iend, 100) == 0
            println("time is t = ", t[1])
            local toc = Base.time()
            println("the time for the simulation is ", toc - tic, " seconds")
            println("extrema are ", extrema(real.(θ)))
            println("on wavenumber index ", index_choice)
        end
        =#

    end
    println("finished realization ", j)
    @. θ̅ += θ / realizations
end

toc = Base.time()
println("the time for the simulation was ", toc - tic, " seconds")

x_A = Array(x)[:] .- 2π
θ_F = Array(real.(θ))
θ̅_F = Array(real.(θ̅))

fig = Figure(resolution=(2048, 512))
ax1 = Axis(fig[1, 1], title="t = 0")
ax2 = Axis(fig[1, 2], title="instantaneous t = " * string(tend))
ax3 = Axis(fig[1, 4], title="ensemble average t = " * string(tend))
println("the extrema of the end field is ", extrema(θ_F))
println("the extrema of the ensemble average is ", extrema(θ̅_F))
colormap = :bone_1
colormap = :nipy_spectral
heatmap!(ax1, x_A, x_A, θ_A, colormap=colormap, colorrange=(0.0, 1.0), interpolate=true)
hm = heatmap!(ax2, x_A, x_A, θ_F, colormap=colormap, colorrange=(0.0, 1.0), interpolate=true)
hm_e = heatmap!(ax3, x_A, x_A, θ̅_F, colormap=colormap, colorrange=(0.0, 0.2), interpolate=true)
Colorbar(fig[1, 3], hm, height=Relative(3 / 4), width=25, ticklabelsize=30, labelsize=30, ticksize=25, tickalign=1,)
Colorbar(fig[1, 5], hm_e, height=Relative(3 / 4), width=25, ticklabelsize=30, labelsize=30, ticksize=25, tickalign=1,)
display(fig)

# indices = abs.(θ_F[:]) .> 1e-3;
# hist(θ_F[indices])
#=
println("saving ", jld_name * string(index_choice) * ".jld2")
θ̅a = Array(real.(θ̅))
xnodes = Array(x)[:]
ynodes = Array(y)[:]
kˣ_wavenumbers = Array(kˣ)[:]
kʸ_wavenumbers = Array(kˣ)[:]
source = Array(s)
jldsave(jld_name * string(index_choice) * ".jld2"; ψ_save, θ_save, θ̅a, κ, xnodes, ynodes, kˣ_wavenumbers, kʸ_wavenumbers, source)
=#

max_stars_repo_name: sandreza/FourierCore

max_issues_repo_name: sandreza/FourierCore

max_forks_repo_name: sandreza/FourierCore


text:
"""
    load(fpath::String;
         gamma::Union{Nothing,Float64}=nothing, expand_paletted::Bool=false, background=false)
    load(s::IO;
         gamma::Union{Nothing,Float64}=nothing, expand_paletted::Bool=false, background=false)

Read a PNG image as a julia `Array`.

# Arguments

- `fpath`: the path to the png file to be loaded.
- `s`: an IO stream containing the png file.

# Keywords

- `gamma`: the end-to-end coefficient for gamma correction, can be used to override the
    automatic gamma correction, a value of `1.0` means no gamma correction. By default,
    gamma correction is applied according to a `gAMA` or `sRGB` chunks, if present;
    if neither chunk is present a value of `0.45455` for the file gamma is assumed.
    Screen gamma is currently always assumed to be `2.2`.
- `background`: can be used to set a background color for transparent images.
    Accepted values are `false` for no background, `true` for loading the image
    with a solid background if `bKGD` chunk is set, as well as user-provided backgrounds
    of type `UInt8`, `Gray` or `RGB`, which represent a palette index, gray and true color
    backgrounds respectively.
- `expand_paletted`: when reading in simple paletted images, i.e. having a `PLTE` chunk and
   an 8 bit depth, the image will be represented as an `IndirectArray` with `OffsetArray`
   `values` field. To always get back a plain `Matrix` of colorants, use `expand_paletted=true`.

# Returns

- `Matrix{Gray{N0f8}}`: for grayscale images (without transparency) and bit depth lower or equal to 8.
- `Matrix{Gray{N0f16}}`: for grayscale images (without transparency) and bit depth equal to 16.
- `Matrix{GrayA{N0f8}}`: for grayscale images (with transparency) and bit depth lower or equal to 8.
- `Matrix{GrayA{N0f16}}`: for grayscale images (with transparency) and bit depth equal to 16.
- `Matrix{RGB{N0f8}}`: for true color images (without transparency) and bit depth equal to 8.
    Also for palleted images without transparency when `expand_paletted=true`.
- `Matrix{RGB{N0f16}}`: for true color images (without transparency) and bit depth equal to 16.
- `Matrix{RGBA{N0f8}}`: for true color images (with transparency) and bit depth equal to 8.
    Also for palleted images without transparency when `expand_paletted=true`.
- `Matrix{RGBA{N0f16}}`: for true color images (with transparency) and bit depth equal to 16.
- `IndirectArray{RGB{N0f8}, 2, UInt8, Matrix{UInt8}, OffsetVector{RGB{N0f8}, Vector{RGB{N0f8}}}}`:
    for palleted images with transparency when `expand_paletted=false`.
- `IndirectArray{RGBA{N0f8}, 2, UInt8, Matrix{UInt8}, OffsetVector{RGBA{N0f8}, Vector{RGBA{N0f8}}}}`:
    for palleted images with transparency when `expand_paletted=false`.
"""
function load(fpath::String; gamma::Union{Nothing,Float64}=nothing, expand_paletted::Bool=false, background=false)
    fp = open_png(fpath)
    png_ptr = create_read_struct()
    @debug "Load PNG File:" fpath png_ptr
    info_ptr = create_info_struct(png_ptr)
    png_init_io(png_ptr, fp)
    png_set_sig_bytes(png_ptr, PNG_BYTES_TO_CHECK)
    out = _load(png_ptr, info_ptr, gamma=gamma, expand_paletted=expand_paletted, background=background)
    close_png(fp)
    return out
end

maybe_lock(f, io::IO) = lock(f, io)
# IOStream doesn't support locking...
maybe_lock(f, io::IOStream) = f()
maybe_lock(f, io::IOBuffer) = f()
maybe_lock(f, io::Base64EncodePipe) = f()
maybe_lock(f, io::IOContext) = maybe_lock(f, io.io)

if VERSION < v"1.6.0-DEV.1652"
    _iswritable(x) = Base.iswritable(x)
    _iswritable(pipe::Base64EncodePipe) = Base.iswritable(pipe.io)
else
    _iswritable(x) = iswritable(x)
end

function load(s::IO; gamma::Union{Nothing,Float64}=nothing, expand_paletted::Bool=false, background=false)
    isreadable(s) || throw(ArgumentError("read failed, IOStream is not readable"))
    Base.eof(s) && throw(EOFError())

    png_ptr = create_read_struct()
    info_ptr = create_info_struct(png_ptr)

    maybe_lock(s) do
        if s isa IOBuffer
            png_set_read_fn(png_ptr, pointer_from_objref(s), readcallback_iobuffer_c[])
        else
            png_set_read_fn(png_ptr, s.handle, readcallback_c[])
        end
        # https://stackoverflow.com/questions/22564718/libpng-error-png-unsigned-integer-out-of-range
        png_set_sig_bytes(png_ptr, 0)
        return _load(png_ptr, info_ptr, gamma=gamma, expand_paletted=expand_paletted, background=background)
    end
end

function _readcallback(png_ptr::png_structp, data::png_bytep, length::png_size_t)::Cvoid
    a = png_get_io_ptr(png_ptr)
    ccall(:ios_readall, Csize_t, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), a, data, length)
    return
end

function _readcallback_iobuffer(png_ptr::png_structp, data::png_bytep, length::png_size_t)::Cvoid
    a = png_get_io_ptr(png_ptr)
    io = unsafe_pointer_to_objref(a)
    unsafe_read(io, Ptr{UInt8}(data), length)
    return
end

function _load(png_ptr, info_ptr; gamma::Union{Nothing,Float64}=nothing, expand_paletted::Bool=false, background=false)
    png_read_info(png_ptr, info_ptr)
    flags = png_get_valid(
        png_ptr, info_ptr,
        PNG_INFO_sRGB | PNG_INFO_gAMA | PNG_INFO_tRNS | PNG_INFO_cHRM | PNG_INFO_bKGD | PNG_INFO_PLTE
    )
    valid_sRGB = (flags & PNG_INFO_sRGB) != 0
    valid_gAMA = (flags & PNG_INFO_gAMA) != 0
    valid_tRNS = (flags & PNG_INFO_tRNS) != 0
    valid_cHRM = (flags & PNG_INFO_cHRM) != 0
    valid_bKGD = (flags & PNG_INFO_bKGD) != 0
    valid_PLTE = (flags & PNG_INFO_PLTE) != 0

    width = png_get_image_width(png_ptr, info_ptr)
    height = png_get_image_height(png_ptr, info_ptr)
    color_type = color_type_orig = png_get_color_type(png_ptr, info_ptr)
    bit_depth = bit_depth_orig = png_get_bit_depth(png_ptr, info_ptr)
    num_channels = png_get_channels(png_ptr, info_ptr)
    interlace_type = png_get_interlace_type(png_ptr, info_ptr)
    background_color = nothing
    is_transparent = valid_tRNS | (color_type & PNG_COLOR_MASK_ALPHA) != 0
    is_gray = (color_type == PNG_COLOR_TYPE_GRAY_ALPHA) | (color_type == PNG_COLOR_TYPE_GRAY)

    read_as_paletted = !expand_paletted && color_type == PNG_COLOR_TYPE_PALETTE && bit_depth == 8 && valid_PLTE
    read_with_background = _check_background_load(background, is_gray, valid_PLTE, read_as_paletted)

    screen_gamma = PNG_DEFAULT_sRGB
    image_gamma = Ref{Cdouble}(-1.0)
    intent = Ref{Cint}(-1)
    if isnothing(gamma)
        if valid_sRGB
            if png_get_sRGB(png_ptr, info_ptr, intent) != 0
                png_set_gamma(png_ptr, screen_gamma, PNG_DEFAULT_sRGB)
            else
                if png_get_gAMA(png_ptr, info_ptr, image_gamma) == 0
                    image_gamma[] = 0.45455
                end
                png_set_gamma(png_ptr, screen_gamma, image_gamma[])
            end
        elseif valid_gAMA
            if png_get_gAMA(png_ptr, info_ptr, image_gamma) == 0
                image_gamma[] = 0.45455
            end
            png_set_gamma(png_ptr, screen_gamma, image_gamma[])
        end
    elseif gamma != 1
        image_gamma[] = 2.2inv(gamma)
        png_set_gamma(png_ptr, screen_gamma, image_gamma[])
    end

    if !read_as_paletted
        if color_type == PNG_COLOR_TYPE_PALETTE
            png_set_palette_to_rgb(png_ptr)
            color_type = PNG_COLOR_TYPE_RGB
        end

        if color_type == PNG_COLOR_TYPE_GRAY && bit_depth < 8
            png_set_expand_gray_1_2_4_to_8(png_ptr)
            png_set_packing(png_ptr)
            bit_depth = 8
        end

        if valid_tRNS
            png_set_tRNS_to_alpha(png_ptr)
            if color_type == PNG_COLOR_TYPE_GRAY || color_type == PNG_COLOR_TYPE_RGB
                color_type |= PNG_COLOR_MASK_ALPHA
            end
        end

        if read_with_background & valid_bKGD & is_transparent
            png_set_filler(png_ptr, 0xff, PNG_FILLER_AFTER)
        end
        buffer_eltype = _buffer_color_type(color_type, bit_depth)

        bit_depth == 16 && png_set_swap(png_ptr)
    else
        buffer_eltype = UInt8
    end

    if read_with_background
        background_color = process_background(png_ptr, info_ptr, _adjust_background_bitdepth(background, bit_depth))
        is_transparent || @warn("Background color for non-transparent image: $(background_color)")
    end

    n_passes = png_set_interlace_handling(png_ptr)
    png_read_update_info(png_ptr, info_ptr)

    # Gamma correction is applied to a palette after `png_read_update_info` is called
    if read_as_paletted
        # TODO: Figure out the length of palette before calling png_get_PLTEs
        #    `png_get_palette_max(png_ptr, info_ptr)` seems like a good option
        #    I can't make it work (I only ever get it to return zero). Note that it
        #    has to be called after png_read_image / png_read_png.
        palette_length = Ref{Cint}()
        palette_buffer = Vector{RGB{N0f8}}(undef, PNG_MAX_PALETTE_LENGTH)
        GC.@preserve palette_buffer begin
            png_get_PLTE(png_ptr, info_ptr, pointer_from_objref(palette_buffer), palette_length)
            palette = palette_buffer[1:palette_length[]]
        end
        if valid_tRNS
            alpha_buffer = Vector{N0f8}(undef, palette_length[])
            alphas_cnt = Ref{Cint}()
            GC.@preserve alpha_buffer begin
                png_get_tRNS(png_ptr, info_ptr, pointer_from_objref(alpha_buffer), alphas_cnt, C_NULL)
            end
            alpha_buffer[alphas_cnt[]+1:palette_length[]] .= one(N0f8)  # palette entries are opaque by default
            palette = map(RGBA, palette, alpha_buffer)
        end
        buffer_eltype = UInt8
    end

    # We transpose to work around libpng expecting row-major arrays
    buffer = Array{buffer_eltype}(undef, width, height)

    @debug(
        "Read PNG info:",
        PNG_HEADER_VERSION_STRING,
        png_ptr,
        height,
        width,
        color_type_orig,
        color_type,
        bit_depth_orig,
        bit_depth,
        num_channels,
        interlace_type,
        gamma,
        image_gamma[],
        screen_gamma,
        background,
        intent[],
        read_as_paletted,
        n_passes,
        buffer_eltype,
        valid_sRGB,
        valid_gAMA,
        valid_tRNS,
        valid_cHRM,
        valid_bKGD,
        valid_PLTE,
    )

    GC.@preserve buffer begin
        buffer = Base.invokelatest(_load!, buffer, png_ptr, info_ptr)
    end

    if expand_paletted || color_type != PNG_COLOR_TYPE_PALETTE
        return buffer
    else
        # We got 0-based indices back from libpng and converting to 1-based could overflow UInt8.
        # Using UInt16 for index would cost us large part of the savings provided by IndirectArray.
        return IndirectArray(buffer, OffsetArray(palette, -1))
    end
end

function _load!(buffer::Matrix{T}, png_ptr, info_ptr) where T    # separate to support precompilation of permutedims
    png_read_image(png_ptr, map(pointer, eachcol(buffer)))
    png_read_end(png_ptr, info_ptr)
    png_destroy_read_struct(Ref{Ptr{Cvoid}}(png_ptr), Ref{Ptr{Cvoid}}(info_ptr), C_NULL)
    return permutedims(buffer, (2, 1))
end

function _buffer_color_type(color_type, bit_depth)
    bit_depth = Int(bit_depth)
    if color_type == PNG_COLOR_TYPE_GRAY
        colors_type = Gray{bit_depth > 8 ? Normed{UInt16,bit_depth} : Normed{UInt8,bit_depth}}
    elseif color_type == PNG_COLOR_TYPE_PALETTE
        colors_type = RGB{bit_depth == 16 ? N0f16 : N0f8}
    elseif color_type == PNG_COLOR_TYPE_RGB
        colors_type = RGB{bit_depth == 16 ? N0f16 : N0f8}
    elseif color_type == PNG_COLOR_TYPE_RGB_ALPHA
        colors_type = RGBA{bit_depth == 16 ? N0f16 : N0f8}
    elseif color_type == PNG_COLOR_TYPE_GRAY_ALPHA
        colors_type = GrayA{bit_depth > 8 ? Normed{UInt16,bit_depth} : Normed{UInt8,bit_depth}}
    else
        throw(error("Unknown color type: $color_type"))
    end
    return colors_type
end

### Write ##########################################################################################

const SupportedPaletteColor = Union{
    AbstractRGB{<:Union{N0f8,AbstractFloat}},
    TransparentRGB{T,<:Union{N0f8,AbstractFloat}} where T,
}
"""
    save(fpath::String, image::AbstractArray;
         compression_level::Int=0, compression_strategy::Int=3, filters::Int=4)
    save(s::IO, image::AbstractArray;
         compression_level::Int=0, compression_strategy::Int=3, filters::Int=4)

Write out a julia `Array` as a PNG image.

# Arguments
- `fpath`: the filesystem path where to store the resulting PNG image.
- `s`: the IO stream to which the resulting PNG image should be written to.
- `image`: the julia `Array` which should be encoded as a PNG image. The type, eltype, and
    dimensionality jointly determine the way the PNG file is encoded. The following table
    roughly summarises the encoding rules:
| **dims**    | **type**                                       | **eltype** (`T`)                                     | **PNG color type**          | **PNG bit depth**                                |
|:------------|:-----------------------------------------------|:-----------------------------------------------------|:----------------------------|-------------------------------------------------:|
| `[h, w]`    | `AbstracArray{T,2}`                            | `Union{AbstractFloat, Unsigned, Normed, Gray}`       | `PNG_COLOR_TYPE_GRAY`       | 8 by default, 16 when `N0f16` or `N4f12` is used |
| `[h, w, 1]` | `AbstracArray{T,3}`                            | `Union{AbstractFloat, Unsigned, Normed, Gray}`       | `PNG_COLOR_TYPE_GRAY`       | 8 by default, 16 when `N0f16` or `N4f12` is used |
| `[h, w, 2]` | `AbstracArray{T,3}`                            | `Union{AbstractFloat, Unsigned, Normed, GrayA}`      | `PNG_COLOR_TYPE_GRAY_ALPHA` | 8 by default, 16 when `N0f16` or `N4f12` is used |
| `[h, w, 3]` | `AbstracArray{T,3}`                            | `Union{AbstractFloat, Unsigned, Normed, RGB, BGR}`   | `PNG_COLOR_TYPE_RGB`        | 8 by default, 16 when `N0f16` or `N4f12` is used |
| `[h, w, 4]` | `AbstracArray{T,3}`                            | `Union{AbstractFloat, Unsigned, Normed, RGBA, ARGB}` | `PNG_COLOR_TYPE_RGB_ALPHA`  | 8 by default, 16 when `N0f16` or `N4f12` is used |
| `[h, w]`    | `IndirectArray{T, 2, S, Matrix{S}, Vector{T}}` | `Union{RGB, BGR, RGBA, ARGB}`                        | `PNG_COLOR_TYPE_PALETTE`    |                                                8 |

``Note``: Only `UInt8` and `UInt16` are supported for `Unsigned`, `Float32` and `Float64` for `AbstractFloat`;
    `N0f8`, `N0f16` and `N4f12` for `Normed`. These types are also the only valid types for colorants.

``Note``: `Bool`s and `BitArray` are also supported, but the bit depth of the image will still be 8.

``Note``: When `image` is an `IndirectArray` with up to 256 unique `RGB` colors, the result is encoded as a paletted image.
    Palletes with 16 bit depths are not supported. The palette (`values` field of the `IndirectArray`) could also be represented with an `OffsetArray`.

# Keywords

- `compression_level`: `0` (`Z_NO_COMPRESSION`), `1` (`Z_BEST_SPEED`), ..., `9` (`Z_BEST_COMPRESSION`)
- `compression_strategy`: 0 (`Z_DEFAULT_STRATEGY`), 1 (`Z_FILTERED`), 2 (`Z_HUFFMAN_ONLY`),
    3 (`Z_RLE`), 4 (`Z_FIXED`)
- `filters`: specify a type of preprocessing applied to each row which can increase its compressability.
    Valid values are 0 (`None`), 1 (`Sub`), 2 (`Up`), 3 (`Average`), 4 (`Paeth`).
- `file_gamma`: the value governing the gamma encoding of the image. When `nothing`,
    the image stored as `sRGB`, otherwise the gamma value provided will populate a `gAMA`
    chunk of the image.
- `background`: optional background color to be stored in the `bKGD` chunk. Only meaningful for transparent images.
    Valid values are `nothing` for no background, `UInt8` as a palette index for palleted images,
    `Gray` for grayscale images and `RGB` for true color images.

# Returns
- `nothing`
"""
function save(
    fpath::String,
    image::S;
    compression_level::Integer = Z_BEST_SPEED,
    compression_strategy::Integer = Z_RLE,
    filters::Integer = Int(PNG_FILTER_PAETH),
    file_gamma::Union{Nothing,Float64} = nothing,
    background::Union{Nothing,UInt8,AbstractGray,AbstractRGB} = nothing,
) where {
    T,
    S<:Union{AbstractMatrix{T},AbstractArray{T,3}}
}
    @assert Z_DEFAULT_STRATEGY <= compression_strategy <= Z_FIXED
    @assert Z_NO_COMPRESSION <= compression_level <= Z_BEST_COMPRESSION
    @assert 2 <= ndims(image) <= 3
    @assert size(image, 3) <= 4

    fp = ccall(:fopen, Ptr{Cvoid}, (Cstring, Cstring), fpath, "wb")
    fp == C_NULL && error("Could not open $(fpath) for writing")

    png_ptr = create_write_struct()
    @debug "Save PNG File:" fpath png_ptr
    info_ptr = create_info_struct(png_ptr)

    png_init_io(png_ptr, fp)

    _save(png_ptr, info_ptr, image,
        compression_level=compression_level,
        compression_strategy=compression_strategy,
        filters=filters,
        file_gamma=file_gamma,
        background=background
    )

    close_png(fp)
    return
end
function save(
    s::IO,
    image::S;
    compression_level::Integer = Z_BEST_SPEED,
    compression_strategy::Integer = Z_RLE,
    filters::Integer = Int(PNG_FILTER_PAETH),
    file_gamma::Union{Nothing,Float64} = nothing,
    background::Union{Nothing,UInt8,AbstractGray,AbstractRGB} = nothing,
) where {
    T,
    S<:Union{AbstractMatrix,AbstractArray{T,3}}
}
    @assert Z_DEFAULT_STRATEGY <= compression_strategy <= Z_FIXED
    @assert Z_NO_COMPRESSION <= compression_level <= Z_BEST_COMPRESSION
    @assert 2 <= ndims(image) <= 3
    @assert size(image, 3) <= 4
    _iswritable(s) || throw(ArgumentError("write failed, IOStream is not writeable"))

    png_ptr = create_write_struct()
    info_ptr = create_info_struct(png_ptr)
    maybe_lock(s) do
        r = Ref{Any}(s)
        GC.@preserve r begin
            png_set_write_fn(png_ptr, r, writecallback_c[], C_NULL)

            _save(png_ptr, info_ptr, image,
                compression_level=compression_level,
                compression_strategy=compression_strategy,
                filters=filters,
                file_gamma=file_gamma,
                background=background,
            )
        end
    end
    return
end

function _save(png_ptr, info_ptr, image::S;
    compression_level::Integer = Z_BEST_SPEED,
    compression_strategy::Integer = Z_RLE,
    filters::Integer = Int(PNG_FILTER_PAETH),
    file_gamma::Union{Nothing,Float64} = nothing,
    background::Union{Nothing,UInt8,AbstractGray,AbstractRGB} = nothing,
) where {
    T,
    S<:Union{AbstractMatrix{T},AbstractArray{T,3}}
}
    image = _enforce_dense_cpu_array(image)
    height, width = size(image)[1:2]
    bit_depth = _get_bit_depth(image)
    color_type = _get_color_type(image)
    is_transparent = (color_type & PNG_COLOR_MASK_ALPHA) != 0
    approx_bytes = round(Int, (height + 1) * width * bit_depth / 8 * (((color_type | PNG_COLOR_MASK_COLOR > 0) ? 3 : 1) + (color_type | PNG_COLOR_MASK_ALPHA > 0)))

    png_set_filter(png_ptr, PNG_FILTER_TYPE_BASE, UInt32(filters))
    png_set_compression_level(png_ptr, compression_level)
    png_set_compression_strategy(png_ptr, compression_strategy)
    png_set_compression_window_bits(png_ptr, min(15, max(8, _nextpow2exp(approx_bytes))))

    if color_type == PNG_COLOR_TYPE_PALETTE
        # TODO: 1, 2, 4 bit-depth indices for palleted
        _png_check_paletted(image)
        palette = image.values
        is_transparent = eltype(palette) <: TransparentRGB
        color_count = length(palette)
        if is_transparent
            alphas = _palette_alpha(palette)
            alpha_count = color_count
            while (alpha_count > 0) && (alphas[alpha_count] == 1)
                alpha_count -= 1
            end
            png_set_PLTE(png_ptr, info_ptr, _standardize_palette(color.(palette)), color_count)
            png_set_tRNS(png_ptr, info_ptr, alphas, alpha_count, C_NULL)
        else
            png_set_PLTE(png_ptr, info_ptr, _standardize_palette(palette), color_count)
        end
    else
        image_eltype = eltype(image)
        if (image_eltype <: BGR || image_eltype <: BGRA || image_eltype <: ABGR || image_eltype <: ARGB32)
            png_set_bgr(png_ptr)
        end

        if (image_eltype <: ABGR || image_eltype <: ARGB)
            png_set_swap_alpha(png_ptr)
        end

        if color_type == PNG_COLOR_TYPE_GRAY && bit_depth < 8
            png_set_packing(png_ptr)
            bit_depth = 8  # TODO: support 1, 2, 4 bit-depth gray images
        end
    end
    if !(background === nothing)
        _check_background_save(background, color_type)
        _png_set_bKGD(png_ptr, info_ptr, _adjust_background_bitdepth(background, bit_depth))
    end

    if file_gamma === nothing
        # gAMA and cHRM chunks should be always present for compatibility with older systems
        png_set_sRGB_gAMA_and_cHRM(png_ptr, info_ptr, PNG_sRGB_INTENT_PERCEPTUAL)
    else
        png_set_gAMA(png_ptr, info_ptr, file_gamma)
    end

    @debug(
        "Write PNG info:",
        png_ptr,
        height,
        width,
        bit_depth,
        color_type,
        filters,
        compression_level,
        compression_strategy,
        palette,
        typeof(image)
    )

    # TODO: on error this throws an abort signal because we currently don't handle `longjmp`
    png_set_IHDR(
        png_ptr,
        info_ptr,
        width,
        height,
        bit_depth,
        color_type,
        PNG_INTERLACE_NONE,
        PNG_COMPRESSION_TYPE_BASE,
        PNG_FILTER_TYPE_BASE,
    )

    png_write_info(png_ptr, info_ptr)
    # Handles endianness for 16 bit, must be set after png_write_info
    bit_depth == 16 && png_set_swap(png_ptr)

    # We transpose to work around libpng expecting row-major arrays
    _write_image(permutedims(_prepare_buffer(image), (2, 1)), png_ptr, info_ptr)

    png_destroy_write_struct(Ref(png_ptr), Ref(info_ptr))
end

function _writecallback(png_ptr::png_structp, data::png_bytep, length::png_size_t)::Csize_t
    a = png_get_io_ptr(png_ptr)
    io = unsafe_load(Ptr{Any}(a))
    unsafe_write(io, data, length)
end

function _write_image(buf::AbstractArray{T,2}, png_ptr::Ptr{Cvoid}, info_ptr::Ptr{Cvoid}) where {T}
    GC.@preserve buf begin
        ccall(
            (:png_write_image, libpng),
            Cvoid,
            (Ptr{Cvoid}, Ptr{Ptr{T}}),
            png_ptr,
            map(pointer, eachcol(buf)),
        )
    end
    png_write_end(png_ptr, info_ptr)
end

function _png_check_paletted(image)
    palette = image.values
    color_count = length(palette)
    color_type = eltype(palette)

    ndims(image) != 2 && throw(ArgumentError("Only 2D `IndirectArrays` are supported"))
    color_count > 256 && throw(ArgumentError("Maximum size of `image.velues` is 256 colors"))
    if !(color_type <: SupportedPaletteColor)
        throw(ArgumentError(
            "Only 8-bit (transparent) RGB colors are supported for paletted images"
        ))
    end
end

_prepare_buffer(x::IndirectArray) = UInt8.(x.index .- first(axes(x.values, 1)))
_prepare_buffer(x::BitArray) = _prepare_buffer(collect(x))
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:Colorant{<:Normed}} = x
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:UInt8} =  reinterpret(Gray{N0f8}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:UInt16} = reinterpret(Gray{N0f16}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:Normed} = reinterpret(Gray{T}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:Gray{Bool}} =  Gray{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:Gray{<:AbstractFloat}} =  Gray{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:GrayA{<:AbstractFloat}} = GrayA{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:RGB{<:AbstractFloat}} =   RGB{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:RGBA{<:AbstractFloat}} =  RGBA{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:BGR{<:AbstractFloat}} =   BGR{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:BGRA{<:AbstractFloat}} =  BGRA{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:ARGB{<:AbstractFloat}} =  ARGB{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:ABGR{<:AbstractFloat}} =  ABGR{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:Union{AbstractFloat,Bool}} = reinterpret(Gray{N0f8}, N0f8.(x))
_prepare_buffer(x::AbstractArray{T,3}) where {T<:Union{AbstractFloat,Bool}} = __prepare_buffer(N0f8.(x))
_prepare_buffer(x::AbstractArray{T,3}) where {T<:Union{UInt8,Int8}} = __prepare_buffer(reinterpret(N0f8, x))
_prepare_buffer(x::AbstractArray{T,3}) where {T<:Union{UInt16,Int16}} = __prepare_buffer(reinterpret(N0f16, x))
_prepare_buffer(x::AbstractArray{T,3}) where {T<:Normed} = __prepare_buffer(x)

function __prepare_buffer(x::AbstractArray{T,3}) where {T}
    nchannels = size(x, 3)
    if nchannels == 1
        ifelse(ndims(x) == 3, _prepare_buffer(dropdims(x, dims=3)), x)
    elseif nchannels == 2
        GrayA.(colorview(GrayA, view(x, :, :, 1), view(x, :, :, 2)))
    elseif nchannels == 3
        RGB.(colorview(RGB, view(x, :, :, 1), view(x, :, :, 2), view(x, :, :, 3)))
    elseif nchannels == 4
        RGBA.(colorview(RGBA, view(x, :, :, 1), view(x, :, :, 2), view(x, :, :, 3), view(x, :, :, 4)))
    else
        error("Unsupported number of channels $(nchannels) in input. Only <= 4 is expected.")
    end
end

# On performance: if the array type has efficient convert method to Array then this is
# almost a no-op
function _enforce_dense_cpu_array(img::AbstractArray)
    if Base.has_offset_axes(img)
        convert(Array, OffsetArrays.no_offset_view(img))
    else
        convert(Array, img)
    end
end
_enforce_dense_cpu_array(img::DenseArray) = img
_enforce_dense_cpu_array(img::OffsetArray) = _enforce_dense_cpu_array(parent(img))
_enforce_dense_cpu_array(img::IndirectArray) = img # PNGFiles has built-in support for this type

_get_color_type(x::AbstractArray{<:Gray}) = PNG_COLOR_TYPE_GRAY
_get_color_type(x::AbstractArray{<:GrayA}) = PNG_COLOR_TYPE_GRAY_ALPHA
_get_color_type(x::AbstractArray{<:AbstractRGB}) = PNG_COLOR_TYPE_RGB
_get_color_type(x::AbstractArray{<:AbstractARGB}) = PNG_COLOR_TYPE_RGBA
_get_color_type(x::AbstractArray{<:AbstractRGBA}) = PNG_COLOR_TYPE_RGBA
_get_color_type(x::IndirectArray) = PNG_COLOR_TYPE_PALETTE
function _get_color_type(
        x::AbstractArray{T, N}
    ) where {
        T<:Union{Normed,Unsigned,Bool,AbstractFloat},
        N
    }
    if N == 2
        return PNG_COLOR_TYPE_GRAY
    elseif N == 3
        d = size(x, 3)
        d == 1 && (return PNG_COLOR_TYPE_GRAY)
        d == 2 && (return PNG_COLOR_TYPE_GRAY_ALPHA)
        d == 3 && (return PNG_COLOR_TYPE_RGB)
        d == 4 && (return PNG_COLOR_TYPE_RGBA)
    end
    error("Number of dimensions $(N) in image not supported (only 2D or 3D Arrays are expected).")
end

_standardize_palette(p::AbstractVector) = _enforce_dense_cpu_array(__standardize_palette(p))
__standardize_palette(p::AbstractVector{<:RGB}) = p
__standardize_palette(p::AbstractVector{<:AbstractRGB}) = RGB.(p)
__standardize_palette(p::AbstractVector{<:RGB{<:AbstractFloat}}) = RGB{N0f8}.(p)
__standardize_palette(p::AbstractVector{<:AbstractRGB{<:AbstractFloat}}) = RGB{N0f8}.(p)
__standardize_palette(p::OffsetArray) = __standardize_palette(parent(p))

_palette_alpha(p::OffsetArray) where {T} = _palette_alpha(collect(p))
_palette_alpha(p::AbstractVector{<:TransparentRGB{T,N0f8}}) where {T} = alpha.(p)
_palette_alpha(p::AbstractVector{<:TransparentRGB{T,<:AbstractFloat}}) where {T} = N0f8.(alpha.(p))


text:
import Pkg
Pkg.activate(".")

using MCDTS
using Test
using Random
using DynamicalSystemsBase
using DelayEmbeddings
import MCDTS.L

# Check Lorenz System
Random.seed!(1234)
ds = Systems.lorenz()
data = trajectory(ds,200)
data = data[2001:2:3750,:]

w1 = DelayEmbeddings.estimate_delay(data[:,1],"mi_min")
w2 = DelayEmbeddings.estimate_delay(data[:,2],"mi_min")
w3 = DelayEmbeddings.estimate_delay(data[:,3],"mi_min")

w = maximum(hcat(w1,w2,w3))

@testset "MCDTS embedding tests" begin

    include("MCDTS_attractor_reconstruction_tests.jl")
    include("MCDTS_convergent_cross_mapping_tests.jl")
    include("MCDTS_nearest_neighbor_prediction_test.jl")

end


text:
using Zygote

struct GaussNum{T<:Real} <: Real
    μ::T
    σ::T
end
mu(x::GaussNum) = x.μ
sig(x::GaussNum) = x.σ
GaussNum(x,y) = GaussNum(promote(x,y)...)
±(x,y) = GaussNum(x,y)
Base.convert(::Type{T}, x::T) where T<:GaussNum = x
Base.convert(::Type{GaussNum{T}}, x::Number) where T = GaussNum(x,zero(T))
Base.promote_rule(::Type{GaussNum{T}}, ::Type{S}) where {T,S} = GaussNum{T}
Base.promote_rule(::Type{GaussNum{T}}, ::Type{GaussNum{T}}) where T = GaussNum{T}

# convert(GaussNum{Float64}, 1.0) |> display
# promote(GaussNum(1.0,1.0), 2.0) |> display
# error()

#+(x::GaussNum{T},a::T) where T =GaussNum(x.μ+a,x.σ)
#+(a::T,x::GaussNum{T}) where T =GaussNum(x.μ+a,x.σ)
#-(x::GaussNum{T},a::T) where T =GaussNum(x.μ-a,x.σ)
#-(a::T,x::GaussNum{T}) where T =GaussNum(x.μ-a,x.σ)
#*(x::GaussNum{T},a::T) where T =GaussNum(x.μ*a,a*x.σ)
#*(a::T,x::GaussNum{T}) where T =GaussNum(x.μ*a,a*x.σ)

# function Base.:*(x1::GaussNum, x2::GaussNum)
#     f(x1,x2) = x1 * x2
#     s1 = Zygote.gradient(μ -> f(μ,x2.μ), x1.μ)[1]^2 * x1.σ^2
#     s2 = Zygote.gradient(μ -> f(x1.μ,μ), x2.μ)[1]^2 * x2.σ^2
#     GaussNum(f(x1.μ,x2.μ), sqrt(s1+s2))
# end

function _uncertain(f, args::GaussNum...)
    μs  = [x.μ for x in args]
    dfs = Zygote.gradient(f,μs...)
    σ   = map(zip(dfs,args)) do (df,x)
        df^2 * x.σ^2
    end |> sum |> sqrt
    GaussNum(f(μs...), σ)
end

function _uncertain(expr::Expr)
    if expr.head == :call
        :(_uncertain($(expr.args[1]), $(expr.args[2:end]...)))
    else
        error("Expression has to be a :call")
    end
end

macro uncertain(expr)
    _uncertain(expr)
end

getmodule(f) = first(methods(f)).module

function _register(func::Symbol)
    mod = getmodule(eval(func))
    :($(mod).$(func)(args::GaussNum...) = _uncertain($func, args...))
end

function _register(funcs::Expr)
    Expr(:block, map(_register, funcs.args)...)
end

macro register(funcs)
    _register(funcs)
end

@register - + *

f(x,y) = x+y*x

# @register *
# @register +
# @register -
# @register f

asdf(x1::GaussNum{T},x2::GaussNum{T}) where T =GaussNum(x1.μ*x2.μ, sqrt((x2.μ*x1.σ).^2 + (x1.μ * x2.σ).^2))
gggg(x1::GaussNum{T},x2::GaussNum{T}) where T =GaussNum(x1.μ+x2.μ, sqrt(x1.σ.^2 + x2.σ.^2))

x1 = GaussNum(rand(),rand())
x2 = GaussNum(rand(),rand())

display(x1*x2)
display(asdf(x1,x2))
display(_uncertain(*,x1,x2))
display(@uncertain x1*x2)

display(x1-x2)
display(x1+x2)
display(f(x1,x2))
#error()

using Plots
using JLD2

abstract type AbstractODEProblem end

struct ODEProblem{F,T,U,P} <: AbstractODEProblem
    f::F
    tspan::T
    u0::U
    θ::P
end

abstract type ODESolver end
struct Euler{T} <: ODESolver
    dt::T
end
struct RK2{T} <: ODESolver
    dt::T
end

function f(x,θ)
    α, β, γ, δ = θ
    x₁, x₂ = x

    dx₁ = α*x₁ - β*x₁*x₂
    dx₂ = δ*x₁*x₂ - γ*x₂

    [dx₁, dx₂]
end

function solve(prob::AbstractODEProblem, solver::ODESolver)
    t = prob.tspan[1]; u = prob.u0
    us = [u]; ts = [t]
    while t < prob.tspan[2]
        (u,t) = solver(prob, u, t)
        push!(us,u)
        push!(ts,t)
    end
    ts, reduce(hcat,us)
end

function (solver::Euler)(prob::ODEProblem, u, t)
    f, θ, dt  = prob.f, prob.θ, solver.dt
    (u + dt*f(u,θ), t+dt)
end

function (solver::RK2)(prob::ODEProblem, u, t)
    f, θ, dt  = prob.f, prob.θ, solver.dt
    uh = u + f(u,θ)*dt
    u + dt/2*(f(u,θ) + f(uh,θ)), t+dt
end

@recipe function plot(ts::AbstractVector, xs::AbstractVector{<:GaussNum})
    # you can set a default value for an attribute with `-->`
    # and force an argument with `:=`
    μs = [x.μ for x in xs]
    σs = [x.σ for x in xs]
    @series begin
        :seriestype := :path
        # ignore series in legend and color cycling
        primary := false
        linecolor := nothing
        fillcolor := :lightgray
        fillalpha := 0.5
        fillrange := μs .- σs
        # ensure no markers are shown for the error band
        markershape := :none
        # return series data
        ts, μs .+ σs
    end
    ts, μs
end

θ = [0.1,0.2,0.3,0.2]
u0 = [GaussNum(1.0,0.1),GaussNum(1.0,0.1)]
tspan = (0.,100.)
dt = 0.1
prob = ODEProblem(f,tspan,u0,θ)

t,X=solve(prob, RK2(0.2))
p1 = plot(t, X[1,:], label="x", lw=3)
plot!(p1, t, X[2,:], label="y", lw=3)

display(p1)


text:
using DynamicPolynomials
# include("D:\\Programs\\blockpop\\TSSOS\\src\\TSSOS.jl")
using TSSOS

# unconstrained optimization using the TSSOS hierarchy
@polyvar x[1:3]
f = 1+x[1]^4+x[2]^4+x[3]^4+x[1]*x[2]*x[3]+x[2]
opt,sol,data = tssos_first(f, x, newton=true, reducebasis=true, TS="block", solution=true, QUIET=true)
# optimum = 0.475274778453039

opt,sol,data = tssos_higher!(data, TS="MD", solution=true, QUIET=true)
# optimum = 0.4752747815682939

# constrained optimization using the TSSOS hierarchy
@polyvar x[1:3]
f = 1+x[1]^4+x[2]^4+x[3]^4+x[1]*x[2]*x[3]+x[2]
g = 1-x[1]^2-2*x[2]^2
h = x[2]^2+x[3]^2-1
pop = [f, g, h]
opt,sol,data = tssos_first(pop, x, 2, numeq=1, TS="block", solution=true, QUIET=true)
# optimum = 0.6861925732538557

opt,sol,data = tssos_higher!(data, TS="MD", solution=true, QUIET=true)
# optimum = 0.6861925611033903

# constrained optimization using the CS-TSSOS hierarchy
n = 6
@polyvar x[1:n]
f = 1+sum(x.^4)+x[1]*x[2]*x[3]+x[3]*x[4]*x[5]+x[3]*x[4]*x[6]+x[3]*x[5]*x[6]+x[4]*x[5]*x[6]
pop = [f, 1-sum(x[1:3].^2), 1-sum(x[3:6].^2)]
d = 2
opt,sol,data = cs_tssos_first(pop, x, d, numeq=1, TS="block", solution=true, QUIET=true)
# optimum = 0.7174253329322346

opt,sol,data = cs_tssos_higher!(data, TS="MD", solution=true, QUIET=true)
# optimum = 0.7174253409780793

println("Run successfully!")


text:
# This should also demonstrate traits, which is
# harder to understand than simple dispatch, but allows for more re-use of code.

# TO DO: Don't change the type when State changes.
# Just use two arguments. The first can be a type.
abstract type AbstractThing end
abstract type Kind end
abstract type State end
Base.@kwdef mutable struct Thing{K<:Kind, S<:State} <: AbstractThing
    x::Length = 0.0m   # Could potentially be a function of 'time' between animation panels.
    y::Length = 0.0m
    w::Length = (1 + √5)m / 2
    h::Length = 1.0m
    name
end
Thing{K}(;kwargs...) where {K<:Kind} = Thing{K}{Transit}(;kwargs...)
struct Section <: Kind end # We may add more types later, imperative style...
struct Mounted <: State end
struct Transit <: State end

height(th::Thing{K, S}) where {K, S} = th.h

Base.@kwdef mutable struct Container{K<:Kind, V <: AbstractVector} <: AbstractThing
    children::Vector = Vector()::V # Intending AbstractThings here
    x::Length = 0.0∙m
    y::Length = 0.0∙m
    w::Length = (1 + √5)m      # Intending masking width
    h::Length = 2.0m           # Intending masking height
    name=""
end
Container{K}(;kwargs...) where {K<:Kind} = Container{K}{Vector}(;kwargs...)

# Interfaces: Use Container.children directly....
Base.length(c::Container{K, V}) where {K, V} = length(c.children)
Base.getindex(c::Container{K, V}, I) where {K, V} = c.children[I]
Base.size(c::Container{K, V}) where {K, V} = length(c.children)
Base.IndexStyle(::Type{<:Container}) = IndexLinear()
Base.iterate(c::Container{K, V}) where {K, V} = isempty(c.children) ? nothing : (c.children[1], 2)
function Base.iterate(c::Container{K, V}, s) where {K, V}
    if (s > length(c))
        return
    end
    return c[s], s + 1
end
Base.lastindex(c::Container{K, V}) where {K, V} = lastindex(c.children)

putin!(c::Container{K, V}, t) where {K, V} = push!(c.children, t)
pop!(c::Container{K, V}) where {K, V} = pop!(c.children)

moveto_x!(th, x) = th.x = x
moveto_y!(th, y) = th.y = y
moveto!(th, x, y) = moveto_x!(th, x), moveto_y!(th, y)
moveto!(th, pt) = moveto!(th, pt.x * scale_pt_to_unit(m), -pt.y * scale_pt_to_unit(m))

move_x!(th, Δx) = th.x += Δx
move_y!(th, Δy) = th.y += Δy
move!(th, Δx, Δy) = move_x!(th, Δx), move_y!(th, Δy)
move!(th, Δpt) = move!(th, Δpt.x * scale_pt_to_unit(m), -Δpt.y * scale_pt_to_unit(m))

set_w!(th, w) = th.w = w
set_h!(th, h) = th.h = h
set_width_height!(th, w, h) = set_w!(th, w), set_h!(th, h)
set_width_height!(th, pt) = set_width_height!(th, pt.x * scale_pt_to_unit(m), -pt.y * scale_pt_to_unit(m))

struct Panel <: Kind end
height(c::Container{Panel, V}) where V = c.h
function drawit(c::Container{Panel, V}) where V
    name = c.name
    x, y, w, h = promote(c.x, c.y, c.w, c.h)
    Δx = w / 2
    @layer begin
        translate(Point(x, y))
        # text box, pixels
        tw = 1.2 * pixelwidth(name)
        # pixel vertical offset, +y is down
        ph = -scale_to_pt(h)
        #
        p1 = Point(-Δx, 0h)
        p2 = Point(-Δx, h)
        p3 = Point(Δx, h)
        tpxr = p3.x - EM
        tpx = tpxr - tw / 2
        tpxl = max(p2.x + EM, tpx - tw / 2)
        p2_1 = Point(tpxl, ph)
        p2_2 = Point(tpx, ph)
        p2_3 = Point(tpxr, ph)
        p4 = Point(Δx, 0h)
        @layer begin
            sethue(PALETTE[10])
            brush(p1, p2)
            brush(p2, p2_1)
            brush(p2_3, p3)
            brush(p3, p4)
            brush(p4, p1)
            text(name, p2_2; valign = :center, halign = :center)
        end
        # add mask
         box(p2 + (5, 5), p4 + (-5, -5), :clip; vertices=false)
        # draw children, with this layer's origo
        for thing in c
            drawit(thing)
        end
    end
end

struct Ground <: Kind end
height(c::Container{Ground, V}) where V = 0.0m
function drawit(c::Container{Ground, V}) where V
    name = c.name
    x, y, w, h = promote(c.x, c.y, c.w, c.h)
    Δx = w / 2
    pw = scale_to_pt(w)
    @layer begin
        translate(Point(x, y))
        ptul = Point(-Δx, 0h)
        groundbox = BoundingBox(rect(ptul, w, h))
        @layer begin
            grc = PALETTE[6]
            grcb = color_with_lumin(PALETTE[6], 10)
            gradient = mesh(
                box(groundbox, vertices=true),
                [grcb, grc, grc, grcb])
            setmesh(gradient)
            paint()
            sethue(grc)
            hatchx = range(0, pw, step = 20)
            ptx = [ptul+ (x, 0.0) for x in hatchx]
            rule.(ptx, π / 4; boundingbox = groundbox)
        end
        # draw children, with this layer's origo
        for thing in c
            drawit(thing)
        end
    end
end

function drawname(pt, name, maxwidth)
    snam = string(name)
    width = pixelwidth(snam)
    if width < scale_to_pt(maxwidth)
        height = row_height()
        draw_background(pt - (width / 2, height / 2), width, height)
        sethue(PALETTE[6])
        text(snam, pt, halign=:center, valign = :middle)
    end
end
struct Stack <: Kind end
height(c::Container{Stack, V}) where V = sum(height, c; init = 0.0m)
"""
putin!(c::Container{Stack}, th::Thing{K}) where K <: Kind

Much like 'push!', but swithes the original object with a 'Mounted' type.
"""
function putin!(c::Container{Stack, V}, th::Thing{K}) where {V, K <: Kind}
    t = Thing{K}{Mounted}(;name = th.name,
                           x = th.x,
                           y = height(c),
                           w = th.w,
                           h = th.h)
    push!(c.children, t)
end

function drawit(c::Container{Stack, V}) where V
    name = c.name
    x, y = promote(c.x, c.y)
    @layer begin
        translate(Point(x, y))
        drawname(Point(0, EM), name, 50EM)
        for thing in c
            drawit(thing)
        end
    end
end

function drawit(th::Thing{Section, S}) where S<:State
    name = th.name
    x, y, w, h = promote(th.x, th.y, th.w, th.h)
    @layer begin
        pt = Point(x, y)  + (x, y)
        center = Point(x, y + h / 2)
        pt0 = pt - ( w / 2, 0.0w )
        pt1 = pt + ( w / 2, 0.0w )
        ble = cylinderblend(pt0, pt1, S)
        setblend(ble)
        squircle(center, w / 2, h / 2; rt = 0.15, action = :fill)
        drawname(center, name, w)
    end
end

function cylinderblend(pt0, pt1, ::Type{T}) where T
    ble = blend(pt0, pt1)
    # Blend colours
    vc = (parse(RGBA, s) for s in ("gold4", "gold3", "gold2", "gold1"))
    c1, c2, c3, c4 = if T <: Mounted
        vc
    else
        map(c-> color_with_saturation.(c, 0.01), vc)
    end
    bcs = [c1, c2, c3, c4, c3, c2, c1]
    # Corresponding path length ∈(0..1)
    bls  = 0.5 .+ 0.5 .* [-1, -0.71, 0.5, 0, 0.5, 0.71, 1]
    for (l, c) in zip(bls, bcs)
        addstop(ble, l, c)
    end
    ble
end