text:
"Take a board and turn each row and column into a set of numbers"
function compose_sets(board)
    sets = []
    for r in eachrow(board)
        push!(sets, Set(r))
    end
    for c in eachcol(board)
        push!(sets, Set(c))
    end
    return sets
end

is_winning(set) = any(s -> isempty(s), set)

function play(bingo; stop_after_first_winner=true)
    sets = compose_sets.(bingo.boards)
    number = nothing
    winner_index = nothing
    winning_number = nothing
    winning_set = nothing
    for number in bingo.numbers
        # Remove number from all sets of all boards
        for board_set in sets
            for s in board_set
                delete!(s, number)
            end
        end
        # Try to find a winner
        winner_indices = findall(is_winning, sets)
        if !isempty(winner_indices)
            winning_number = number
            winning_set = sets[winner_indices[1]]
            deleteat!(sets, winner_indices)
            stop_after_first_winner && break
        end
    end
    return winning_number, reduce(union, winning_set)
end

function day04a(bingo)
    number, remaining = play(bingo)
    return number * sum(remaining)
end

function day04b(bingo)
    number, remaining = play(bingo, stop_after_first_winner = false)
    return number * sum(remaining)
end

function day04(fpath)
    bingo = parse_bingo(fpath)
    return (day04a(bingo), day04b(bingo))
end


text:
using Test

@testset "ITensorInfiniteMPS.jl" begin
  @testset "$file" for file in readdir(@__DIR__)
    if startswith(file, "test_") && endswith(file, ".jl")
      println("Running $file")
      include(file)
    end
  end
end

ext: jl


text:
using JLD2, Interpolations
@load joinpath(@__DIR__,"src/exactcore.jld2") ψi
@load joinpath(@__DIR__,"src/ansatzcore.jld2") ψa

psii_coefs = ψi.coefs
psii_knots = ψi.knots
psia_coefs = ψa.coefs
psia_knots = ψa.knots
@save joinpath(@__DIR__,"src/data.jld2") psii_coefs psii_knots psia_coefs psia_knots

## save in new format
@load joinpath(@__DIR__,"src/data.jld2") psii_coefs psii_knots psia_coefs psia_knots
ψi = interpolate(psii_knots, psii_coefs, Gridded(Linear()))
ψa = interpolate(psia_knots, psia_coefs, Gridded(Linear()))
@save joinpath(@__DIR__,"src/cores.jld2") ψi ψa

## load interp and check type
@load joinpath(@__DIR__,"src/cores.jld2") ψi ψa


text:
"""
    DeployConfig

Abstract type which new deployment configs should be subtypes of.
"""
abstract type DeployConfig end

"""
    Documenter.documenter_key(cfg::DeployConfig)

Return the Base64-encoded SSH private key for the repository.
Defaults to reading the `DOCUMENTER_KEY` environment variable.

This method must be supported by configs that push with SSH, see
[`Documenter.authentication_method`](@ref).
"""
function documenter_key(::DeployConfig)
    return ENV["DOCUMENTER_KEY"]
end

"""
    Documenter.git_tag(cfg::DeployConfig)

Return the git tag of the build. If the build is not for a tag, return `nothing`.

This function determines the subfolder where the built docs are deployed.
Either a folder named as the tag (e.g. `vX.Y.Z`) or `devurl`
(see [`deploydocs`](@ref)) for non-tag builds.
"""
git_tag(::DeployConfig) = nothing

"""
    Documenter.should_deploy(cfg::DeployConfig; repo=repo, devbranch=devbranch)

Return `true` if the current build should deploy, and `false` otherwise.
This function is called with the `repo` and `devbranch` arguments from
[`deploydocs`](@ref).
"""
should_deploy(::DeployConfig; kwargs...) = false
should_deploy(::Nothing; kwargs...) = false # when auto-detection fails

@enum AuthenticationMethod SSH HTTPS

"""
    Documenter.authentication_method(::DeployConfig)

Return enum instance `SSH` or `HTTPS` depending on push method to be used.

Configs returning `SSH` should support [`Documenter.documenter_key`](@ref).
Configs returning `HTTPS` should support [`Documenter.authenticated_repo_url`](@ref).
"""
authentication_method(::DeployConfig) = SSH

"""
    Documenter.authenticated_repo_url(cfg::DeployConfig)

Return an authenticated URL to the upstream repository.

This method must be supported by configs that push with HTTPS, see
[`Documenter.authentication_method`](@ref).
"""
function authenticated_repo_url end

#############
# Travis CI #
#############

"""
    Travis <: DeployConfig

Default implementation of `DeployConfig`.

The following environment variables influences the build
when using the `Travis` configuration:

 - `DOCUMENTER_KEY`: must contain the Base64-encoded SSH private key for the repository.
   This variable should be set in the Travis settings for the repository. Make sure this
   variable is marked **NOT** to be displayed in the build log.

 - `TRAVIS_PULL_REQUEST`: must be set to `false`.
   This avoids deployment on pull request builds.

 - `TRAVIS_REPO_SLUG`: must match the value of the `repo` keyword to [`deploydocs`](@ref).

 - `TRAVIS_EVENT_TYPE`: may not be set to `cron`. This avoids re-deployment of existing
   docs on builds that were triggered by a Travis cron job.

 - `TRAVIS_BRANCH`: unless `TRAVIS_TAG` is non-empty, this must have the same value as
   the `devbranch` keyword to [`deploydocs`](@ref). This makes sure that only the
   development branch (commonly, the `master` branch) will deploy the "dev" documentation
   (deployed into a directory specified by the `devurl` keyword to [`deploydocs`](@ref)).

 - `TRAVIS_TAG`: if set, a tagged version deployment is performed instead; the value
   must be a valid version number (i.e. match `Base.VERSION_REGEX`). The documentation for
   a package version tag gets deployed to a directory named after the version number in
   `TRAVIS_TAG` instead.

The `TRAVIS_*` variables are set automatically on Travis. More information on how Travis
sets the `TRAVIS_*` variables can be found in the
[Travis documentation](https://docs.travis-ci.com/user/environment-variables/#default-environment-variables).
"""
struct Travis <: DeployConfig
    travis_branch::String
    travis_pull_request::String
    travis_repo_slug::String
    travis_tag::String
    travis_event_type::String
end
function Travis()
    travis_branch       = get(ENV, "TRAVIS_BRANCH",        "")
    travis_pull_request = get(ENV, "TRAVIS_PULL_REQUEST",  "")
    travis_repo_slug    = get(ENV, "TRAVIS_REPO_SLUG",     "")
    travis_tag          = get(ENV, "TRAVIS_TAG",           "")
    travis_event_type   = get(ENV, "TRAVIS_EVENT_TYPE",    "")
    return Travis(travis_branch, travis_pull_request,
        travis_repo_slug, travis_tag, travis_event_type)
end

# Check criteria for deployment
function should_deploy(cfg::Travis; repo, devbranch, kwargs...)
    ## The deploydocs' repo should match TRAVIS_REPO_SLUG
    repo_ok = occursin(cfg.travis_repo_slug, repo)
    ## Do not deploy for PRs
    pr_ok = cfg.travis_pull_request == "false"
    ## If a tag exist it should be a valid VersionNumber
    tag_ok = isempty(cfg.travis_tag) || occursin(Base.VERSION_REGEX, cfg.travis_tag)
    ## If no tag exists deploydocs' devbranch should match TRAVIS_BRANCH
    branch_ok = !isempty(cfg.travis_tag) || cfg.travis_branch == devbranch
    ## DOCUMENTER_KEY should exist (just check here and extract the value later)
    key_ok = haskey(ENV, "DOCUMENTER_KEY")
    ## Cron jobs should not deploy
    type_ok = cfg.travis_event_type != "cron"
    all_ok = repo_ok && pr_ok && tag_ok && branch_ok && key_ok && type_ok
    marker(x) = x ? "✔" : "✘"
    @info """Deployment criteria for deploying with Travis:
    - $(marker(repo_ok)) ENV["TRAVIS_REPO_SLUG"]="$(cfg.travis_repo_slug)" occurs in repo="$(repo)"
    - $(marker(pr_ok)) ENV["TRAVIS_PULL_REQUEST"]="$(cfg.travis_pull_request)" is "false"
    - $(marker(tag_ok)) ENV["TRAVIS_TAG"]="$(cfg.travis_tag)" is (i) empty or (ii) a valid VersionNumber
    - $(marker(branch_ok)) ENV["TRAVIS_BRANCH"]="$(cfg.travis_branch)" matches devbranch="$(devbranch)" (if tag is empty)
    - $(marker(key_ok)) ENV["DOCUMENTER_KEY"] exists
    - $(marker(type_ok)) ENV["TRAVIS_EVENT_TYPE"]="$(cfg.travis_event_type)" is not "cron"
    Deploying: $(marker(all_ok))
    """
    return all_ok
end

# Obtain git tag for the build
function git_tag(cfg::Travis)
    isempty(cfg.travis_tag) ? nothing : cfg.travis_tag
end

##################
# GitHub Actions #
##################

"""
    GitHubActions <: DeployConfig

Implementation of `DeployConfig` for deploying from GitHub Actions.

The following environment variables influences the build
when using the `GitHubActions` configuration:

 - `DOCUMENTER_KEY`: must contain the Base64-encoded SSH private key for the repository.
   This variable should be set in the GitHub Actions configuration file using a repository
   secret, see the documentation for
   [secret environment variables](https://help.github.com/en/articles/virtual-environments-for-github-actions#creating-and-using-secrets-encrypted-variables).

 - `GITHUB_EVENT_NAME`: must be set to `push`.
   This avoids deployment on pull request builds.

 - `GITHUB_REPOSITORY`: must match the value of the `repo` keyword to [`deploydocs`](@ref).

 - `GITHUB_REF`: must match the `devbranch` keyword to [`deploydocs`](@ref), alternatively
   correspond to a git tag.

The `GITHUB_*` variables are set automatically on GitHub Actions, see the
[documentation](https://help.github.com/en/articles/virtual-environments-for-github-actions#default-environment-variables).
"""
struct GitHubActions <: DeployConfig
    github_repository::String
    github_event_name::String
    github_ref::String
end
function GitHubActions()
    github_repository = get(ENV, "GITHUB_REPOSITORY", "") # "JuliaDocs/Documenter.jl"
    github_event_name = get(ENV, "GITHUB_EVENT_NAME", "") # "push", "pull_request" or "cron" (?)
    github_ref        = get(ENV, "GITHUB_REF",        "") # "refs/heads/$(branchname)" for branch, "refs/tags/$(tagname)" for tags
    return GitHubActions(github_repository, github_event_name, github_ref)
end

# Check criteria for deployment
function should_deploy(cfg::GitHubActions; repo, devbranch, kwargs...)
    ## The deploydocs' repo should match GITHUB_REPOSITORY
    repo_ok = occursin(cfg.github_repository, repo)
    ## Do not deploy for PRs
    pr_ok = cfg.github_event_name == "push"
    ## If a tag exist it should be a valid VersionNumber
    m = match(r"^refs/tags/(.*)$", cfg.github_ref)
    tag_ok = m === nothing ? false : occursin(Base.VERSION_REGEX, String(m.captures[1]))
    ## If no tag exists deploydocs' devbranch should match the current branch
    m = match(r"^refs/heads/(.*)$", cfg.github_ref)
    branch_ok = m === nothing ? false : String(m.captures[1]) == devbranch
    ## GITHUB_ACTOR should exist (just check here and extract the value later)
    actor_ok = haskey(ENV, "GITHUB_ACTOR")
    ## GITHUB_TOKEN should exist (just check here and extract the value later)
    token_ok = haskey(ENV, "GITHUB_TOKEN")
    all_ok = repo_ok && pr_ok && (tag_ok || branch_ok) && actor_ok && token_ok
    marker(x) = x ? "✔" : "✘"
    @info """Deployment criteria for deploying with GitHub Actions:
    - $(marker(repo_ok)) ENV["GITHUB_REPOSITORY"]="$(cfg.github_repository)" occurs in repo="$(repo)"
    - $(marker(pr_ok)) ENV["GITHUB_EVENT_NAME"]="$(cfg.github_event_name)" is "push"
    - $(marker(tag_ok || branch_ok)) ENV["GITHUB_REF"]="$(cfg.github_ref)" corresponds to a tag or matches devbranch="$(devbranch)"
    - $(marker(actor_ok)) ENV["GITHUB_ACTOR"] exists
    - $(marker(token_ok)) ENV["GITHUB_TOKEN"] exists
    Deploying: $(marker(all_ok))
    """
    return all_ok
end

authentication_method(::GitHubActions) = HTTPS
function authenticated_repo_url(cfg::GitHubActions)
    return "https://$(ENV["GITHUB_ACTOR"]):$(ENV["GITHUB_TOKEN"])@github.com/$(cfg.github_repository).git"
end

# Obtain git tag for the build
function git_tag(cfg::GitHubActions)
    m = match(r"^refs/tags/(.*)$", cfg.github_ref)
    return m === nothing ? nothing : String(m.captures[1])
end

##################
# Auto-detection #
##################
function auto_detect_deploy_system()
    if haskey(ENV, "TRAVIS_REPO_SLUG")
        return Travis()
    elseif haskey(ENV, "GITHUB_REPOSITORY")
        return GitHubActions()
    else
        return nothing
    end
end


text:
module StaticArraysLite

using Base: @propagate_inbounds

using Freeze
using LinearAlgebra

export StaticArrayLite, SArrayLite, MArrayLite

abstract type StaticArrayLite{T,N} <: AbstractArray{T, N} end

include("mapreduce.jl")
include("linalg.jl")

struct SArrayLite{S<:Tuple,T,N,L} <: StaticArrayLite{T, N}
    data::NTuple{L,T}
end

Base.length(::SArrayLite{<:Any, <:Any, <:Any, L}) where {L} = L
Base.size(::SArrayLite{Tuple{}}) = ()
Base.size(::SArrayLite{Tuple{S1}}) where S1 = (S1,)
Base.size(::SArrayLite{Tuple{S1, S2}}) where {S1, S2} = (S1, S2)
Base.axes(::SArrayLite{Tuple{}}) = ()
Base.axes(::SArrayLite{Tuple{S1}}) where S1 = (Base.OneTo(S1),)
Base.axes(::SArrayLite{Tuple{S1, S2}}) where {S1, S2} = (Base.OneTo(S1), Base.OneTo(S2))
Base.IndexStyle(::SArrayLite) = Base.IndexLinear()
Freeze.issettable(::SArrayLite) = false

@propagate_inbounds function Base.getindex(a::SArrayLite, i::Int)
    a.data[i]
end

mutable struct MArrayLite{S<:Tuple,T,N,L} <: StaticArrayLite{T, N}
    data::NTuple{L,T}

    function MArrayLite{S,T,N,L}(x::NTuple{L,T}) where {S,T,N,L}
        new{S,T,N,L}(x)
    end

    function MArrayLite{S,T,N,L}(::UndefInitializer) where {S,T,N,L}
        new{S,T,N,L}()
    end
end

Base.length(::MArrayLite{<:Any, <:Any, <:Any, L}) where {L} = L
Base.size(::MArrayLite{Tuple{}}) = ()
Base.size(::MArrayLite{Tuple{S1}}) where S1 = (S1,)
Base.size(::MArrayLite{Tuple{S1, S2}}) where {S1, S2} = (S1, S2)
Base.axes(::MArrayLite{Tuple{}}) = ()
Base.axes(::MArrayLite{Tuple{S1}}) where S1 = (Base.OneTo(S1),)
Base.axes(::MArrayLite{Tuple{S1, S2}}) where {S1, S2} = (Base.OneTo(S1), Base.OneTo(S2))
Base.IndexStyle(::MArrayLite) = Base.IndexLinear()
Freeze.issettable(::MArrayLite) = true

@propagate_inbounds function Base.getindex(a::MArrayLite, i::Int)
    @boundscheck checkbounds(a,i)
    T = eltype(a)

    if isbitstype(T)
        return GC.@preserve a unsafe_load(Base.unsafe_convert(Ptr{T}, pointer_from_objref(a)), i)
    end
    a.data[i]
end

@propagate_inbounds function Base.setindex!(a::MArrayLite, val, i::Int)
    @boundscheck checkbounds(a,i)
    T = eltype(a)

    if isbitstype(T)
        GC.@preserve a unsafe_store!(Base.unsafe_convert(Ptr{T}, pointer_from_objref(a)), convert(T, val), i)
    else
        # This one is unsafe (#27)
        # unsafe_store!(Base.unsafe_convert(Ptr{Ptr{Nothing}}, pointer_from_objref(v.data)), pointer_from_objref(val), i)
        error("setindex!() with non-isbitstype eltype is not supported by StaticArrays. Consider using SizedArray.")
    end

    return val
end

Freeze.thaw(a::SArrayLite{S,T,N,L}) where {S,T,N,L} = MArrayLite{S,T,N,L}(a.data)
Freeze.freeze(a::MArrayLite{S,T,N,L}) where {S,T,N,L} = SArrayLite{S,T,N,L}(a.data)
Base.similar(a::StaticArrayLite, ::Type{T}, sz::Tuple{}) where {T} = MArrayLite{Tuple{},T,0,1}(undef)
@inline Base.similar(a::StaticArrayLite, ::Type{T}, sz::Tuple{Int}) where {T} = MArrayLite{Tuple{sz[1]},T,1,sz[1]}(undef)
@inline Base.similar(a::StaticArrayLite, ::Type{T}, sz::Tuple{Int,Int}) where {T} = MArrayLite{Tuple{sz[1], sz[2]},T,2,sz[1]*sz[2]}(undef)

end # module


text:
module boustrophedon

export Node, nodepush!, readboustrophedon

"""
A node in a tree structure. Contains a value, a left node, and a right node.
"""
mutable struct Node{T}
    value::T
    left::Union{Node{T}, Nothing}
    right::Union{Node{T}, Nothing}
end

"A more terse form of the default constructor, no need to specify left/right."
function Node(value::T) where T <: Any
    Node(value, nothing, nothing)
end

"""
    nodepush!(node, value)

Push a new value onto a tree, starting with `node`. Will attempt to insert the new
value as a node into the first empty 'slot' on the tree, searching top to bottom,
left to right.
"""
function nodepush!(node::Node{T}, value::T) where T <: Any
    current_row = [node]
    while true
        next_row = []
        for node in current_row
            if isnothing(node.left)
                node.left = Node(value)
                return
            else
                push!(next_row, node.left)
            end
            if isnothing(node.right)
                node.right = Node(value)
                return
            else
                push!(next_row, node.right)
            end
        end # for
        current_row = next_row
    end # while
end

"""
    readboustrophedon(head)

Given a link to a tree node `head`, read all the descendant values in that tree in
'boustrophedon' order, that is, alternating left-to-right and right-to-left down the
levels of the tree.
"""
function readboustrophedon(head::Node{T})::Vector{T} where T <: Any
    output = []
    current_row = [head]
    ltr = true

    while !isempty(current_row)
        next_row = []
        buffer = []
        for node in current_row
            # Builds a buffer in the direction needed for the output, if outputting
            # this row left-to-right, append values to the right side of the buffer.
            # If outputting this row right-to-left, push values to the left side of
            # the buffer.
            if ltr; push!(buffer, node.value); else; pushfirst!(buffer, node.value); end
            isnothing(node.left) || push!(next_row, node.left)
            isnothing(node.right) || push!(next_row, node.right)
        end
        append!(output, buffer)
        ltr = !ltr
        current_row = next_row
    end
    return output
end

end # module

module test

using Test;
using ..boustrophedon;

@testset "Sample Test" begin
    input = Node('a')
    nodepush!(input, 'b')
    nodepush!(input, 'c')
    nodepush!(input, 'd')
    nodepush!(input, 'e')
    nodepush!(input, 'f')
    nodepush!(input, 'g')
    @test readboustrophedon(input) == collect("acbdefg")

    input = Node('p')
    foreach(c -> nodepush!(input, c), collect("orgnosnoitacit"))
    @test readboustrophedon(input) == collect("prognostication")
end

end # module

max_stars_repo_name: JenniferToops/CodingDojo

max_issues_repo_name: JenniferToops/CodingDojo

max_forks_repo_name: JenniferToops/CodingDojo


text:
#####
##### Dipole kernel
#####

"""
    dipole_kernel(sz, vsz; kwargs...) =
        dipole_kernel(Float64, sz, vsz; kwargs...)

    dipole_kernel(
        ::Type{T<:AbstractFloat},
        sz::NTuple{3, Integer},
        vsz::NTuple{3, Real};
        bdir::NTuple{3, Real} = (0, 0, 1),
        method::Symbol = :kspace,
        dsz::NTuple{3, Integer} = sz,
        transform::Union{Nothing, Symbol} = nothing,
        shift::Bool = false
    ) -> Array{T, 3}

Dipole kernel.

By default the dipole kernel is constructed in k-space and centered at index `(1,1,1)`.

### Arguments
- `sz::NTuple{3, Integer}`: array size
- `vsz::NTuple{3, Real}`: voxel size

### Keywords
- `bdir::NTuple{3, Real} = (0, 0, 1)`: unit vector of B field direction
- `method::Symbol = :kspace`: create in image space `(:i, :ispace)` or
    kspace `(:k, :kspace)`
- `dsz::NTuple{3, Integer} = sz`:
    - `method ∈ (:k, :kspace)`: unused
    - `method ∈ (:i, :ispace)`: dipole kernel size in image space
- `transform::Union{Nothing, Symbol} = nothing`:
    - `method ∈ (:k, :kspace)`: create `:rfft` or `:fft` kspace dipole kernel
    - `method ∈ (:i, :ispace)`: transform image space dipole kernel `(:rfft, :fft)`
- `shift::Bool = false`:
    - `method ∈ (:k, :kspace)`:
        kernel centered at `(1,1,1)` (`false`) or `N÷2+1` (`true`)
    - `method ∈ (:i, :ispace) and transform = nothing`:
        kernel centered at `N÷2+1` (`false`) or `(1,1,1)` (`true`)
    - `method ∈ (:i, :ispace) and transform ∈ (:rfft, :fft)`:
        kernel centered at `(1,1,1)` (`false`) or `N÷2+1` (`true`)

### Returns
- `Array{T, 3}`: dipole kernel
"""
dipole_kernel(sz, vsz; kwargs...) =
    dipole_kernel(Float64, sz, vsz; kwargs...)

function dipole_kernel(
    ::Type{T},
    sz::NTuple{3, Integer},
    vsz::NTuple{3, Real};
    bdir::NTuple{3, Real} = (0, 0, 1),
    method::Symbol = :kspace,
    dsz::NTuple{3, Integer} = sz,
    transform::Union{Nothing, Symbol} = nothing,
    shift::Bool = false
) where {T<:AbstractFloat}
    !any(<=(0), sz)     || throw(ArgumentError("invalid array size"))
    !any(<=(0), vsz)    || throw(DomainError(vsz, "vsz must be > 0"))
    !iszero(norm(bdir)) || throw(DomainError(bdir, "norm of bdir must not be zero"))
    !any(<(0), dsz)     || throw(DomainError(dsz, "dsz must be ≥ 0"))

    method ∈ (:k, :kspace, :i, :ispace) ||
        throw(ArgumentError("method must be one of :k, :kspace, :i, :ispace, got :$(method)"))

    transform === nothing || transform ∈ (:fft, :rfft) ||
        throw(ArgumentError("transform must be one of nothing, :rfft, :fft, got :$(transform)"))

    if method == :k || method == :kspace
        D = Array{T, 3}(undef, transform == :rfft ? (sz[1]>>1 + 1, sz[2], sz[3]) : sz)
        return _dipole_kernel!(D, sz, vsz, bdir, :k; shift=shift)

    elseif method == :i || method == :ispace
        if transform === nothing
            d = Array{T, 3}(undef, sz)
            return _dipole_kernel!(d, dsz, vsz, bdir, :i; shift=shift)

        elseif transform == :fft
            FFTW.set_num_threads(num_cores())
            D = Array{T, 3}(undef, sz)
            D̂ = Array{complex(T), 3}(undef, sz)
            P = plan_fft!(D̂)

            D = _dipole_kernel!(D, D̂, dsz, vsz, bdir, P)
            return shift ? fftshift(D) : D

        elseif transform == :rfft
            FFTW.set_num_threads(num_cores())
            D = Array{T, 3}(undef, (sz[1]>>1 + 1, sz[2], sz[3]))
            D̂ = Array{complex(T), 3}(undef, (sz[1]>>1 + 1, sz[2], sz[3]))
            d = Array{T, 3}(undef, sz)
            P = plan_rfft(d)

            D = _dipole_kernel!(D, D̂, d, dsz, vsz, bdir, P)
            return shift ? fftshift(D) : D
        end
    end

    return nothing
end

function _dipole_kernel!(
    D::AbstractArray{T, 3},
    D̂::AbstractArray{Complex{T}, 3},
    d::AbstractArray{T, 3},
    sz::NTuple{3, Integer},
    vsz::NTuple{3, Real},
    bdir::NTuple{3, Real},
    P::Union{FFTW.rFFTWPlan{T, -1}, FFTW.cFFTWPlan{Complex{T}, -1}},
    method::Symbol,
    transform::Symbol,
) where {T<:AbstractFloat}
    if method == :k
        D = _dipole_kernel!(D, size(d), vsz, bdir, method)
    elseif method == :i && transform == :fft
        D = _dipole_kernel!(D, D̂, sz, vsz, bdir, P)
    elseif method == :i && transform == :rfft
        D = _dipole_kernel!(D, D̂, d, sz, vsz, bdir, P)
    end
    return D
end

function _dipole_kernel!(
    D::AbstractArray{T, 3},
    D̂::AbstractArray{Complex{T}, 3},
    d::AbstractArray{T, 3},
    sz::NTuple{3, Integer},
    vsz::NTuple{3, Real},
    bdir::NTuple{3, Real},
    P::Union{FFTW.rFFTWPlan{T, -1}, FFTW.cFFTWPlan{Complex{T}, -1, false}}
) where {T<:AbstractFloat}
    d = _dipole_kernel!(d, sz, vsz, bdir, :i, shift=true)

    D̂ = mul!(D̂, P, d)
    D = _tcopyto!(real, D, D̂)

    return D
end

function _dipole_kernel!(
    D::AbstractArray{T, 3},
    D̂::AbstractArray{Complex{T}, 3},
    sz::NTuple{3, Integer},
    vsz::NTuple{3, Real},
    bdir::NTuple{3, Real},
    P::FFTW.cFFTWPlan{Complex{T}, -1, true}
) where {T<:AbstractFloat}
    D̂ = _dipole_kernel!(D̂, sz, vsz, bdir, :i, shift=true)

    D̂ = P*D̂
    D = _tcopyto!(real, D, D̂)

    return D
end

function _dipole_kernel!(
    D::AbstractArray{T, 3},
    sz::NTuple{3, Integer},
    vsz::NTuple{3, Real},
    bdir::NTuple{3, Real},
    method::Symbol;
    shift::Bool = false,
) where {T<:Union{AbstractFloat, Complex{<:AbstractFloat}}}
    nx, ny, nz = size(D)

    dx = convert(NTuple{3, Float64}, vsz)
    B = SVector{3, Float64}(bdir)

    if !(norm(B) ≈ 1)
        B = B ./ norm(B)
    end

    if method == :k
        X, Y, Z = _grid(sz, dx, fft=true, shift=shift)
        a = 1 / 3

        @inbounds @batch for k in 1:nz
            for j in 1:ny
                for i in 1:nx
                    k̂ = SVector{3, Float64}(X[i], Y[j], Z[k])
                    kz = k̂⋅B
                    k2 = k̂⋅k̂
                    D[i,j,k] = a - (kz*kz)/k2
                end
            end
        end
        # set DC term (k̂⋅k̂ := 0)
        D[1] = 0

    elseif method == :i
        dsz = map(n -> iseven(n) ? n-1 : n, sz)
        Xd, Yd, Zd = _grid(dsz, dx)
        X1, X2 = extrema(Xd)
        Y1, Y2 = extrema(Yd)
        Z1, Z2 = extrema(Zd)

        X, Y, Z = _grid(size(D), dx, shift=shift)
        a = prod(dx) * inv(4*pi)

        @inbounds @batch for k in 1:nz
            z = Z[k]
            zout = z < Z1 || z > Z2
            @fastpow for j in 1:ny
                y = Y[j]
                zyout = zout || y < Y1 || y > Y2
                for i in 1:nx
                    x = X[i]
                    zyxout = zyout || x < X1 || x > X2

                    if zyxout
                        D[i,j,k] = 0
                    else
                        r = SVector{3, Float64}(x, y, z)
                        rz = r⋅B
                        r2 = r⋅r
                        D[i,j,k] = a * (3*rz^2-r2)/sqrt(r2^5)
                    end
                end
            end
        end
        # Lorentz sphere correction
        D[1] = 0
    end

    return D
end

#####
##### Spherical Mean Value kernel
#####

"""
    smv_kernel(sz, vsz, r; kwargs...) =
        smv_kernel(Float64, sz, vsz, r; kwargs...)

    smv_kernel(
        ::Type{T<:AbstractFloat},
        sz::NTuple{3, Integer},
        vsz::NTuple{3, Real}
        r::Real;
        transform::Union{Nothing, Symbol} = nothing,
        shift::Bool = false
    ) -> Array{T, 3}

Spherical mean value kernel (SMV).

### Arguments
- `sz::NTuple{3, Integer}`: array size
- `vsz::NTuple{3, Real}`: voxel size
- `r::Real`: radius of sphere in units of `vsz`

### Keywords
- `transform::Union{Nothing, Symbol} = nothing`:
    transform SMV kernel `(:rfft, :fft)`
- `shift::Bool = false`:
    - `transform = nothing`:
        sphere centered at `N÷2+1` (`false`) or `(1,1,1)` (`true`)
    - `transform ∈ (:rfft, :fft)`:
        sphere centered at `(1,1,1)` (`false`) or `N÷2+1` (`true`)

### Returns
- `Array{T, 3}`: SMV kernel
"""
smv_kernel(sz, vsz, r; kwargs...) =
    smv_kernel(Float64, sz, vsz, r; kwargs...)

function smv_kernel(
    ::Type{T},
    sz::NTuple{3, Integer},
    vsz::NTuple{3, Real},
    r::Real;
    transform::Union{Nothing, Symbol} = nothing,
    shift::Bool = false
) where {T<:AbstractFloat}
    !any(<=(0), sz)  || throw(ArgumentError("invalid array size"))
    !any(<=(0), vsz) || throw(DomainError(vsz, "vsz must be > 0"))
    r >= 0           || throw(DomainError(r, "r must be ≥ 0"))

    transform === nothing || transform ∈ (:fft, :rfft) ||
        throw(ArgumentError("transform must be one of nothing, :rfft, :fft, got :$(transform)"))

    if transform === nothing
        s = Array{T, 3}(undef, sz)
        return _smv_kernel!(s, vsz, r, shift=shift)

    elseif transform == :fft
        FFTW.set_num_threads(num_cores())
        S = Array{T, 3}(undef, sz)
        Ŝ = Array{complex(T), 3}(undef, sz)
        P = plan_fft!(Ŝ)

        S  = _smv_kernel!(S, Ŝ, vsz, r, P)
        return shift ? fftshift(S) : S

    elseif transform == :rfft
        FFTW.set_num_threads(num_cores())
        S = Array{T, 3}(undef, (sz[1]>>1 + 1, sz[2], sz[3]))
        Ŝ = Array{complex(T), 3}(undef, (sz[1]>>1 + 1, sz[2], sz[3]))
        s = Array{T, 3}(undef, sz)
        P = plan_rfft(s)

        S = _smv_kernel!(S, Ŝ, s, vsz, r, P)
        return shift ? fftshift(S) : S
    end

    return nothing
end

function _smv_kernel!(
    s::AbstractArray{T, 3},
    vsz::NTuple{3, Real},
    r::Real;
    shift::Bool = false
) where {T<:AbstractFloat}
    s = _sphere!(s, vsz, r, shift=shift)

    # normalizing
    a = inv(sum(s))
    s = _tcopyto!(x -> a*x, s, s)

    return s
end

function _smv_kernel!(
    S::AbstractArray{T, 3},
    Ŝ::AbstractArray{Complex{T}, 3},
    s::AbstractArray{T, 3},
    vsz::NTuple{3, Real},
    r::Real,
    P::Union{FFTW.rFFTWPlan{T, -1}, FFTW.cFFTWPlan{Complex{T}, -1, false}}
) where {T<:AbstractFloat}
    s = _sphere!(s, vsz, r, shift=true)

    # normalization factor
    a = inv(sum(s))

    # fft, discard imaginary (even function -> imag = 0), and normalize
    Ŝ = mul!(Ŝ, P, s)
    S = _tcopyto!(x -> a*real(x), S, Ŝ)

    return S
end

function _smv_kernel!(
    S::AbstractArray{T, 3},
    Ŝ::AbstractArray{Complex{T}, 3},
    vsz::NTuple{3, Real},
    r::Real,
    P::FFTW.cFFTWPlan{Complex{T}, -1, true}
) where {T<:AbstractFloat}
    Ŝ = _sphere!(Ŝ, vsz, r, shift=true)

    # normalization factor
    a = inv(sum(Ŝ))

    # fft, discard imaginary (even function -> imag = 0), and normalize
    Ŝ = P*Ŝ
    S = _tcopyto!(x -> a*real(x), S, Ŝ)

    return S
end

function _sphere!(
    s::AbstractArray{T, 3},
    vsz::NTuple{3, Real},
    r::Real;
    shift::Bool = false
) where {T<:Union{AbstractFloat, Complex{<:AbstractFloat}}}
    sz = size(s)
    nx, ny, nz = sz

    dx = convert(NTuple{3, Float64}, vsz)
    r2 = convert(Float64, r*r)

    X, Y, Z = _grid(sz, dx, shift=shift)

    _zero = zero(T)
    _one = one(T)

    @inbounds for k in 1:nz
        for j in 1:ny
            for i in 1:nx
                R = SVector{3, Float64}(X[i], Y[j], Z[k])
                s[i,j,k] = ifelse(R⋅R <= r2, _one, _zero)
            end
        end
    end

    return s
end

#####
##### Laplace kernel
#####

"""
    laplace_kernel(vsz) =
        laplace_kernel(Float64, vsz)

    laplace_kernel(sz, vsz; kwargs...) =
        laplace_kernel(Float64, sz, vsz; kwargs...)

    laplace_kernel(::Type{T<:AbstractFloat}, vsz::NTuple{3, Real}) =
        laplace_kernel(T, (3, 3, 3), vsz, transform=nothing, shift=false)

    laplace_kernel(
        ::Type{T<:AbstractFloat},
        sz::NTuple{3, Integer},
        vsz::NTuple{3, Real};
        negative::Bool = false,
        transform::Union{Nothing, Symbol} = nothing,
        shift::Bool = false
    ) -> Array{T, 3}

Discrete 7-point stencil Laplacian kernel.

### Arguments
- `sz::NTuple{3, Integer}`: array size
- `vsz::NTuple{3, Real}`: voxel size

### Keywords
- `negative::Bool = false`: construct negative Laplacian (`true`)
- `transform::Union{Nothing, Symbol} = nothing`:
    transform Laplacian kernel `(:rfft, :fft)`
- `shift::Bool = false`:
    - `transform = nothing`:
        sphere centered at `N÷2+1` (`false`) or `(1,1,1)` (`true`)
    - `transform ∈ (:rfft, :fft)`:
        sphere centered at `(1,1,1)` (`false`) or `N÷2+1` (`true`)

### Returns
- `Array{T, 3}`: Laplacian kernel
"""
laplace_kernel(vsz) =
    laplace_kernel(Float64, vsz)

laplace_kernel(sz, vsz; kwargs...) =
    laplace_kernel(Float64, sz, vsz; kwargs...)

laplace_kernel(::Type{T}, vsz::NTuple{3, Real}) where {T<:AbstractFloat} =
    laplace_kernel(T, (3, 3, 3), vsz, transform=nothing, shift=false)

function laplace_kernel(
    ::Type{T},
    sz::NTuple{3, Integer},
    vsz::NTuple{3, Real};
    negative::Bool = false,
    transform::Union{Nothing, Symbol} = nothing,
    shift::Bool = false
) where {T<:AbstractFloat}
    !any(<=(0), sz)  || throw(ArgumentError("invalid array size"))
    !any(<=(0), vsz) || throw(DomainError(vsz, "vsz must be > 0"))

    transform === nothing || transform ∈ (:fft, :rfft) ||
        throw(ArgumentError("transform must be one of nothing, :rfft, :fft, got :$(transform)"))

    if transform === nothing
        Δ = Array{T, 3}(undef, sz)
        return _laplace_kernel!(Δ, vsz, negative=negative, shift=shift)

    elseif transform == :fft
        FFTW.set_num_threads(num_cores())
        L = Array{T, 3}(undef, sz)
        L̂ = Array{complex(T), 3}(undef, sz)
        P = plan_fft!(L̂)

        L = _laplace_kernel!(L, L̂, vsz, P, negative=negative)
        return shift ? fftshift(L) : L

    elseif transform == :rfft
        FFTW.set_num_threads(num_cores())
        L = Array{T, 3}(undef, (sz[1]>>1 + 1, sz[2], sz[3]))
        L̂ = Array{complex(T), 3}(undef, (sz[1]>>1 + 1, sz[2], sz[3]))
        Δ = Array{T, 3}(undef, sz)
        P = plan_rfft(Δ)

        L = _laplace_kernel!(L, L̂, Δ, vsz, P, negative=negative)
        return shift ? fftshift(L) : L
    end

    return nothing
end

function _laplace_kernel!(
    L::AbstractArray{T, 3},
    L̂::AbstractArray{Complex{T}, 3},
    Δ::AbstractArray{T, 3},
    vsz::NTuple{3, Real},
    P::Union{FFTW.rFFTWPlan{T, -1}, FFTW.cFFTWPlan{Complex{T}, -1, false}};
    negative::Bool = false,
) where {T<:AbstractFloat}
    Δ = _laplace_kernel!(Δ, vsz, negative=negative, shift=true)

    L̂ = mul!(L̂, P, Δ)
    L = _tcopyto!(real, L, L̂)

    return L
end

function _laplace_kernel!(
    L::AbstractArray{T, 3},
    L̂::AbstractArray{Complex{T}, 3},
    vsz::NTuple{3, Real},
    P::FFTW.cFFTWPlan{Complex{T}, -1, true};
    negative::Bool = false
) where {T<:AbstractFloat}
    L̂ = _laplace_kernel!(L̂, vsz, negative=negative, shift=true)

    L̂ = P*L̂
    L = _tcopyto!(real, L, L̂)

    return L
end

function _laplace_kernel!(
    Δ::AbstractArray{T, 3},
    vsz::NTuple{3, Real};
    negative::Bool = false,
    shift::Bool = false,
) where {T<:Union{AbstractFloat, Complex{<:AbstractFloat}}}
    hx, hy, hz = convert(NTuple{3, Float64}, inv.(vsz.*vsz))
    D = -2*(hx + hy + hz)

    if negative
        D = -D
        hx = -hx
        hy = -hy
        hz = -hz
    end

    Δ = tfill!(Δ, zero(T))

    if shift
        Δ[1,1,1]   = D
        Δ[2,1,1]   = hx
        Δ[end,1,1] = hx
        Δ[1,2,1]   = hy
        Δ[1,end,1] = hy
        Δ[1,1,2]   = hz
        Δ[1,1,end] = hz
    else
        i, j, k = size(Δ).>>1 .+ 1
        Δ[i,j,k-1] = hz
        Δ[i,j-1,k] = hy
        Δ[i-1,j,k] = hx
        Δ[i,j,k]   = D
        Δ[i+1,j,k] = hx
        Δ[i,j+1,k] = hy
        Δ[i,j,k+1] = hz
    end

    return Δ
end

#####
##### Utility
#####

function _grid(
    sz::NTuple{N, Integer},
    dx::NTuple{N, Real};
    fft::Bool = false,
    shift::Bool = false
) where {N}
    return ntuple(Val(N)) do i
        if fft
            # AbstractFFTs.fftfreq(n, fs) = Frequencies((n+1) >> 1, n, fs/n)
            X = Frequencies((sz[i]+1) >> 1, sz[i], inv(sz[i]*dx[i]))
            collect(shift ? fftshift(X) : X)
        else
            X = Frequencies((sz[i]+1) >> 1, sz[i], dx[i])
            collect(shift ? X : fftshift(X))
        end
    end
end


text:
function table_lookup(datin,datout,in)
    out=datout[1]
    if in<=datin[1]
        done=1
    else
        done=0
    end
    k=1
    while done==0
        k=k+1
        if in<=datin[k]
            out=datout[k-1]+(datout[k]-datout[k-1])/(datin[k]-datin[k-1])*(in-datin[k-1])
            done=1
        elseif k>=length(datin)
            out=datout[k]
            done=1
        end
    end
    return out
end

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text:
# mapobs

struct MappedData{F,D} <: AbstractDataContainer
    f::F
    data::D
end

Base.show(io::IO, data::MappedData) = print(io, "mapobs($(data.f), $(summary(data.data)))")
Base.show(io::IO, data::MappedData{F,<:AbstractArray}) where {F} =
    print(io, "mapobs($(data.f), $(ShowLimit(data.data, limit=80)))")
Base.length(data::MappedData) = numobs(data.data)
Base.getindex(data::MappedData, idx::Int) = data.f(getobs(data.data, idx))
Base.getindex(data::MappedData, idxs::AbstractVector) = data.f.(getobs(data.data, idxs))

"""
    mapobs(f, data)

Lazily map `f` over the observations in a data container `data`.
```julia
data = 1:10
getobs(data, 8) == 8
mdata = mapobs(-, data)
getobs(mdata, 8) == -8
```
"""
mapobs(f, data) = MappedData(f, data)
mapobs(f::typeof(identity), data) = data

"""
    mapobs(fs, data)

Lazily map each function in tuple `fs` over the observations in data container `data`.
Returns a tuple of transformed data containers.
"""
mapobs(fs::Tuple, data) = Tuple(mapobs(f, data) for f in fs)

struct NamedTupleData{TData,F} <: AbstractDataContainer
    data::TData
    namedfs::NamedTuple{F}
end

Base.length(data::NamedTupleData) = numobs(getfield(data, :data))

function Base.getindex(data::NamedTupleData{TData,F}, idx::Int) where {TData,F}
    obs = getobs(getfield(data, :data), idx)
    namedfs = getfield(data, :namedfs)
    return NamedTuple{F}(f(obs) for f in namedfs)
end

Base.getproperty(data::NamedTupleData, field::Symbol) =
    mapobs(getproperty(getfield(data, :namedfs), field), getfield(data, :data))

Base.show(io::IO, data::NamedTupleData) =
    print(io, "mapobs($(getfield(data, :namedfs)), $(getfield(data, :data)))")

"""
    mapobs(namedfs::NamedTuple, data)

Map a `NamedTuple` of functions over `data`, turning it into a data container
of `NamedTuple`s. Field syntax can be used to select a column of the resulting
data container.

```julia
data = 1:10
nameddata = mapobs((x = sqrt, y = log), data)
getobs(nameddata, 10) == (x = sqrt(10), y = log(10))
getobs(nameddata.x, 10) == sqrt(10)
```
"""
function mapobs(namedfs::NamedTuple, data)
    return NamedTupleData(data, namedfs)
end

# filterobs

"""
    filterobs(f, data)

Return a subset of data container `data` including all indices `i` for
which `f(getobs(data, i)) === true`.

```julia
data = 1:10
numobs(data) == 10
fdata = filterobs(>(5), data)
numobs(fdata) == 5
```
"""
function filterobs(f, data; iterfn = _iterobs)
    return obsview(data, [i for (i, obs) in enumerate(iterfn(data)) if f(obs)])
end

_iterobs(data) = [getobs(data, i) for i = 1:numobs(data)]

# groupobs

"""
    groupobs(f, data)

Split data container data `data` into different data containers, grouping
observations by `f(obs)`.

```julia
data = -10:10
datas = groupobs(>(0), data)
length(datas) == 2
```
"""
function groupobs(f, data)
    groups = Dict{Any,Vector{Int}}()
    for i = 1:numobs(data)
        group = f(getobs(data, i))
        if !haskey(groups, group)
            groups[group] = [i]
        else
            push!(groups[group], i)
        end
    end
    return Dict(group => obsview(data, idxs) for (group, idxs) in groups)
end

# joinumobs

struct JoinedData{T,N} <: AbstractDataContainer
    datas::NTuple{N,T}
    ns::NTuple{N,Int}
end

JoinedData(datas) = JoinedData(datas, numobs.(datas))

Base.length(data::JoinedData) = sum(data.ns)

function Base.getindex(data::JoinedData, idx)
    for (i, n) in enumerate(data.ns)
        if idx <= n
            return getobs(data.datas[i], idx)
        else
            idx -= n
        end
    end
end

"""
    joinobs(datas...)

Concatenate data containers `datas`.

```julia
data1, data2 = 1:10, 11:20
jdata = joinumobs(data1, data2)
getobs(jdata, 15) == 15
```
"""
joinobs(datas...) = JoinedData(datas)

"""
    shuffleobs([rng], data)

Return a "subset" of `data` that spans all observations, but
has the order of the observations shuffled.

The values of `data` itself are not copied. Instead only the
indices are shuffled. This function calls [`obsview`](@ref) to
accomplish that, which means that the return value is likely of a
different type than `data`.

```julia
# For Arrays the subset will be of type SubArray
@assert typeof(shuffleobs(rand(4,10))) <: SubArray

# Iterate through all observations in random order
for x in eachobs(shuffleobs(X))
    ...
end
```

The optional parameter `rng` allows one to specify the
random number generator used for shuffling. This is useful when
reproducible results are desired. By default, uses the global RNG.
See `Random` in Julia's standard library for more info.

For this function to work, the type of `data` must implement
[`numobs`](@ref) and [`getobs`](@ref). See [`ObsView`](@ref)
for more information.
"""
shuffleobs(data) = shuffleobs(Random.GLOBAL_RNG, data)

function shuffleobs(rng::AbstractRNG, data)
    obsview(data, randperm(rng, numobs(data)))
end

max_stars_repo_name: JuliaML/MLBase

max_issues_repo_name: JuliaML/MLBase


text:
"""
    get_ripsers()

Download latest versions of Ripser and Cubical Ripser from GitHub and build them.
"""
function get_ripsers()
    exedir = joinpath(@__DIR__, "../exes")
    if isfile("$exedir/ripser-mod2") && isfile("$exedir/ripser-coef")
        @info "ripser exists"
    else
        rm("ripser-mod2"; force=true)
        rm("ripser-coef"; force=true)
        run(`git clone https://github.com/Ripser/ripser`)
        cd("ripser")
        run(`make all`)
        mv("ripser", "$exedir/ripser-mod2")
        mv("ripser-coeff", "$exedir/ripser-coef")
        cd("..")
        rm("ripser", force=true, recursive=true)
    end
    if isfile("$exedir/CR2") && isfile("$exedir/CR3")
        @info "cubical ripser exists"
    else
        for d in (2, 3)
            rm("CR$d"; force=true)
            run(`git clone https://github.com/CubicalRipser/CubicalRipser_$(d)dim`)
            cd("CubicalRipser_$(d)dim")
            run(`make`)
            mv("CR$d", "$exedir/CR$d")
            cd("..")
            rm("CubicalRipser_$(d)dim"; force=true, recursive=true)
        end
    end
end

"""
    ripsers_cmd(filename, string=false; dim_max, modulus, threshold)

Generate a Ripser or Cubical Ripser command for benchmarking, taking file extension and
keyword args into account. If `string=false`, return a `Cmd`, else return a vector of
`String`s.
"""
function ripsers_cmd(filename, string=false; dim_max=1, modulus=2, threshold=nothing)
    ext = splitext(filename)[2]
    exedir = joinpath(@__DIR__, "../exes")
    if ext == ".dipha"
        if dim_max == 1
            cmd_strings = ["$exedir/CR2", filename]
        elseif dim_max == 2
            cmd_strings = ["$exedir/CR3", filename]
        else
            error("invalid dim_max for cubical ripser")
        end
    else
        thresh = isnothing(threshold) ? String[] : ["--threshold", "$threshold"]
        if ext == ".dist"
            format = "distance"
        elseif ext == ".spdist"
            format = "sparse"
        elseif ext == ".pts"
            format = "point-cloud"
        else
            error("invalid extension `$ext`")
        end
        ripser = modulus == 2 ? "$exedir/ripser-mod2" : "$exedir/ripser-coef"
        mod = modulus == 2 ? String[] : ["--modulus", "modulus"]
        cmd = "$ripser --dim $dim_max $mod $thresh $format $(filename)"
        cmd_strings = [
            [ripser, "--dim", "$dim_max"]; mod; thresh; ["--format", format, filename]
        ]
    end
    if string
        return cmd_strings
    else
        return Cmd(cmd_strings)
    end
end

function ripser_benchmark(b::Benchmark)
    return @benchmarkable read(
        $(ripsers_cmd(b.filename; b.kwargs..., b.threshold))
    ) gcsample=true samples=5 seconds=1000
end


text:
using FortranStrings
using Test

@testset "FortranString{T} creating         " begin

    @testset "FortranString{UInt8} " begin

        s = FortranStrings.FortranString{UInt8}("abc")

        @test length(s) == 3
        @test s         == F8"abc"
        @test s         == F"abc"
        @test s         == "abc"
        @test F"abc"    == s
        @test F8"abc"   == s
        @test "abc"     == s
        @test s         == FortranStrings.FortranString{UInt8}(UInt8[0x61, 0x62, 0x63])
        @test string(s) == "abc"
        @test String(s) == "abc"
        @test s[1:1]    == 'a'
        @test s[2:2]    == 'b'
        @test s[3:3]    == 'c'
        @test s.data    == UInt8[0x61, 0x62, 0x63]

        io = IOBuffer()
        show(io, s)
        @test String(take!(io)) == "F8\"abc\""

        show(io, typeof(s))
        @test String(take!(io)) == "FortranString{UInt8}"
    end

    @testset "FortranString{Char}  " begin

        s = FortranStrings.FortranString{Char}("abc")

        @test length(s) == 3
        @test s         == F8"abc"
        @test s         == F"abc"
        @test s         == "abc"
        @test F"abc"    == s
        @test F8"abc"   == s
        @test "abc"     == s
        @test s         == FortranStrings.FortranString{Char}(['a', 'b', 'c'])
        @test string(s) == "abc"
        @test String(s) == "abc"
        @test s[1:1]    == 'a'
        @test s[2:2]    == 'b'
        @test s[3:3]    == 'c'
        @test s.data    == ['a', 'b', 'c']

        io = IOBuffer()
        show(io, s)
        @test String(take!(io)) == "F\"abc\""

        show(io, typeof(s))
        @test String(take!(io)) == "FortranString{Char}"
    end

    @testset "FortranString{UInt64}" begin

        s = FortranStrings.FortranString{UInt64}("abc")

        @test length(s) == 3
        @test s         == F8"abc"
        @test s         == F"abc"
        @test s         == "abc"
        @test F"abc"    == s
        @test F8"abc"   == s
        @test "abc"     == s
        @test s         == FortranStrings.FortranString{UInt64}(UInt64[0x0061, 0x0062, 0x0063])
        @test string(s) == "abc"
        @test_throws MethodError String(s) == "abc" # String(<:AbstractVector{UInt64}) not exist
        @test s[1:1]    == 'a'
        @test s[2:2]    == 'b'
        @test s[3:3]    == 'c'
        @test s.data    == UInt64[0x0061, 0x0062, 0x0063]

        io = IOBuffer()
        show(io, s)
        @test String(take!(io)) == "FUInt64\"abc\""

        show(io, typeof(s))
        @test String(take!(io)) == "FortranString{UInt64}"
    end

    @testset "Quoting              " begin
        for T in (Char, UInt8, UInt64)
            @eval S = FortranString{$T}
            @eval begin
                @test       $S("''abc''")             ==             $S([0x27, 0x27, 0x61, 0x62, 0x63, 0x27, 0x27])
                @test     $S("\"\"abc\"\"")           ==             $S([0x22, 0x22, 0x61, 0x62, 0x63, 0x22, 0x22])
                @test       $S("''abc''", "q")        ==                   $S([0x27, 0x61, 0x62, 0x63, 0x27])
                @test     $S("''''abc''''", "q")      ==             $S([0x27, 0x27, 0x61, 0x62, 0x63, 0x27, 0x27])
                @test       $S("''abc''", "d")        ==             $S([0x27, 0x27, 0x61, 0x62, 0x63, 0x27, 0x27])
                @test     $S("''''abc''''", "d")      == $S([0x27, 0x27, 0x27, 0x27, 0x61, 0x62, 0x63, 0x27, 0x27, 0x27, 0x27])
                @test       $S("''abc''", "qq")       ==             $S([0x27, 0x27, 0x61, 0x62, 0x63, 0x27, 0x27])
                @test     $S("''''abc''''", "qq")     == $S([0x27, 0x27, 0x27, 0x27, 0x61, 0x62, 0x63, 0x27, 0x27, 0x27, 0x27])
                @test     $S("\"\"abc\"\"", "q")      ==             $S([0x22, 0x22, 0x61, 0x62, 0x63, 0x22, 0x22])
                @test $S("\"\"\"\"abc\"\"\"\"", "q")  == $S([0x22, 0x22, 0x22, 0x22, 0x61, 0x62, 0x63, 0x22, 0x22, 0x22, 0x22])
                @test     $S("\"\"abc\"\"", "d")      ==                   $S([0x22, 0x61, 0x62, 0x63, 0x22])
                @test $S("\"\"\"\"abc\"\"\"\"", "d")  ==             $S([0x22, 0x22, 0x61, 0x62, 0x63, 0x22, 0x22])
                @test     $S("\"\"abc\"\"", "qq")     ==                   $S([0x22, 0x61, 0x62, 0x63, 0x22])
                @test $S("\"\"\"\"abc\"\"\"\"", "qq") ==             $S([0x22, 0x22, 0x61, 0x62, 0x63, 0x22, 0x22])

                @test_throws ArgumentError $S("abc", "qqq")
                @test_throws ArgumentError $S("abc", "dd")
                @test_throws ArgumentError $S("abc", "z")
            end
        end
    end

end

@testset "FortranString{T} broadcasting     " begin

    for T in (Char, UInt8, UInt64)
        @eval S = FortranString{$T}

        @eval begin
            #
            @test (incr(x)=x+0x01 ; r = incr.($S("abc"))               ; (r, typeof(r))) == ($S("bcd"), $S)
            @test (incr(x)=x+0x01 ; r = incr.(incr.($S("abc")))        ; (r, typeof(r))) == ($S("cde"), $S)
            @test (incr(x)=x+0x01 ; r = incr.(incr.(incr.($S("abc")))) ; (r, typeof(r))) == ($S("def"), $S)
            @test_throws DimensionMismatch $S("abc") .== $S("ab")
            @test_throws DimensionMismatch $S("abc") .== $S("abcd")
            @test (r = $S("abc").==$S("abc")          ;  r) == BitVector([1,1,1])
            @test (r = $S("abc").==$S("abd")          ;  r) == BitVector([1,1,0])
            @test (r = $S("abc").!=$S("abc")          ;  r) == BitVector([0,0,0])
            @test (r = $S("abc").!=$S("abd")          ;  r) == BitVector([0,0,1])
            @test (r = $S("abc").<=$S("abc")          ;  r) == BitVector([1,1,1])
            @test (r = $S("abc").<=$S("abd")          ;  r) == BitVector([1,1,1])
            @test (r = $S("abc").<$S("abc")           ;  r) == BitVector([0,0,0])
            @test (r = $S("abc").<$S("abd")           ;  r) == BitVector([0,0,1])
            @test (r = $S("abc").>=$S("abc")          ;  r) == BitVector([1,1,1])
            @test (r = $S("abc").>=$S("abd")          ;  r) == BitVector([1,1,0])
            @test (r = $S("abc").>$S("abc")           ;  r) == BitVector([0,0,0])
            @test (r = $S("abc").>$S("abd")           ;  r) == BitVector([0,0,0])
            @test (r = $S("abc").==$S("abd").==$S("cbd").==$S("cba");  r) == BitVector([0,1,0])
            @test (r = $S("abc").==$S("abd").==$S("cbd").=="cba".==$S("cba");  r) == BitVector([0,1,0])
            $T != Char && @test (r = $S("abc").+$S("ABC").-$S("abc"); (r, typeof(r))) == ($S("ABC"), $S)
            @test (r = $S("abc").+FortranString{UInt64}("ABC").-FortranString{UInt8}("abc");  r) == $S("ABC")
            $T != Char && @test (r = $S("abc").+   'A'   .-$S("abc");  r) == $S("AAA")
            @test (r = $S("abc").+   0x41  .-$S("abc");  r) == $S("AAA")
            @test (r = $S("abc").+   65    .-$S("abc");  r) == $S("AAA")

            if $T != Char
                s1 = $S(collect(1+1:6+1))
                s2 = $S(collect(1+2:6+2))
                s3 = $S(collect(1+3:6+3))
                s4 = $S(collect(1+4:6+4))
                s5 = $S(collect(1+5:6+5))
                s6 = $S(collect(1+6:6+6))
                s7 = $S(collect(1+7:6+7))
                s8 = $S(collect(1+8:6+8))
                @test (s1.+s2)                         == (s1+s2)
                @test (s1.+s2.+s3)                     == (s1+s2+s3)
                @test (s1.+s2.+s3.+s4)                 == (s1+s2+s3+s4)
                @test (s1.+s2.+s3.+s4.+s5)             == (s1+s2+s3+s4+s5)
                @test (s1.+s2.+s3.+s4.+s5.+s6)         == (s1+s2+s3+s4+s5+s6)
                @test (s1.+s2.+s3.+s4.+s5.+s6.+s7)     == (s1+s2+s3+s4+s5+s6+s7)
                @test (s1.+s2.+s3.+s4.+s5.+s6.+s7.+s8) == (s1+s2+s3+s4+s5+s6+s7+s8)
            end
        end

        for STR in (String, FortranString{Char}, FortranString{UInt8}, FortranString{UInt64})
            @eval begin
                @test (s = $S("abc"); s      .= $STR("")               ; (s, typeof(s))) == ($S("   "), $S)
                @test (s = $S("abc"); s[1:3] .= $STR("")               ; (s, typeof(s))) == ($S("   "), $S)
                @test (s = $S("abc"); s[1:2] .= $STR("")               ; (s, typeof(s))) == ($S("  c"), $S)
                @test (s = $S("abc"); s[2:2] .= $STR("")               ; (s, typeof(s))) == ($S("a c"), $S)
                @test (s = $S("abc"); s      .= $STR("1")              ; (s, typeof(s))) == ($S("1  "), $S)
                @test (s = $S("abc"); s[2:3] .= $STR("1")              ; (s, typeof(s))) == ($S("a1 "), $S)
                @test (s = $S("abc"); s[1:3] .= $STR("1")              ; (s, typeof(s))) == ($S("1  "), $S)
                @test (s = $S("abc"); s      .= $STR("123")            ; (s, typeof(s))) == ($S("123"), $S)
                @test (s = $S("abc"); s      .= $STR("123456789")      ; (s, typeof(s))) == ($S("123"), $S)
                @test (s = $S("abc"); s[1:3] .= $STR("123")            ; (s, typeof(s))) == ($S("123"), $S)
                @test (s = $S("abc"); s[1:2] .= $STR("123")            ; (s, typeof(s))) == ($S("12c"), $S)
                @test (s = $S("abc"); s[2:2] .= $STR("123")            ; (s, typeof(s))) == ($S("a1c"), $S)
                @test (s = $S("abc"); s[2:3] .= $STR("123")            ; (s, typeof(s))) == ($S("a12"), $S)
                @test (s = $S("abc"); s      .= $STR("1")[1:1]         ; (s, typeof(s))) == ($S("1  "), $S)
                @test (s = $S("abc"); s[2:3] .= $STR("1")              ; (s, typeof(s))) == ($S("a1 "), $S)
                @test (s = $S("abc"); s[1:3] .= $STR("1")[1:1]         ; (s, typeof(s))) == ($S("1  "), $S)
                @test (s = $S("abc"); s      .= $STR("123")[1:1]       ; (s, typeof(s))) == ($S("1  "), $S)
                @test (s = $S("abc"); s      .= $STR("123456789")[5:7] ; (s, typeof(s))) == ($S("567"), $S)
                @test (s = $S("abc"); s[2:3] .= $STR("123456789")[5:7] ; (s, typeof(s))) == ($S("a56"), $S)
                @test (s = $S("abc"); s[1:3] .= $STR("123")[2:3]       ; (s, typeof(s))) == ($S("23 "), $S)
                @test (s = $S("abc"); s[1:2] .= $STR("123")[2:3]       ; (s, typeof(s))) == ($S("23c"), $S)
                @test (s = $S("abc"); s[2:2] .= $STR("123")[2:3]       ; (s, typeof(s))) == ($S("a2c"), $S)
                @test (s = $S("abc"); s[2:3] .= $STR("123")[2:3]       ; (s, typeof(s))) == ($S("a23"), $S)
            end
        end
    end

end

@testset "Some functions on FortranString{T}" begin

    @testset "concatenating    " begin

        s8 = F8"abc"
        s  = F"123"

        @test "abc" * s8 * "abc"             == F8"abcabcabc"
        @test '1' * s8 * '3'                 == F8"1abc3"
        @test typeof("abc" * s8 * "abc")     == FortranString{UInt8}
        @test typeof("abc" * s * "abc")      == FortranString{Char}
        @test typeof("abc" * s * s8 * "abc") == FortranString{Char}
        @test typeof("abc" * s8 * s * "abc") == FortranString{UInt8}

    end

    @testset "assigning        " begin

        s = F8"abc"

        @test (s[2] =  '2'; s) == F8"a2c"
        @test (s[3] = 0x20; s) == F8"a2 "
        @test (s[3] = 0   ; s) == FortranStrings.FortranString{UInt8}(UInt8[0x61, 0x32, 0x00])

    end

    @testset "vector functions " begin

        s = F8"abc"

        for T in (Char, UInt8, UInt64)
            @eval S = FortranString{$T}
            @eval begin
                s = $S("abc")
                @test map(Char, s) == ['a', 'b', 'c']
                @test collect(s)   == $T['a', 'b', 'c']
            end
        end
    end

    @testset "Fortran functions" begin

        for S in (String, FortranString{Char}, FortranString{UInt8}, FortranString{UInt64})
            @eval begin
                s = $S("     ")
                @test LEN(s)               == 5
                @test LEN_TRIM(s)          == 0
                @test LNBLNK(s)            == 0
                @test TRIM(s)              == $S("")
                @test isa(TRIM(s), $S==String ? AbstractString : $S)
                @test strip(s)             == $S("")
                @test isa(strip(s), $S==String ? AbstractString : $S)

                s = $S(" abc ")
                @test LEN(s)               == 5
                @test LEN_TRIM(s)          == 4
                @test LNBLNK(s)            == 4
                @test TRIM(s)              == $S(" abc")
                @test isa(TRIM(s), $S==String ? AbstractString : $S)
                @test strip(s)             == $S("abc")
                @test isa(strip(s), $S==String ? AbstractString : $S)

                @test REPEAT(s, 2)         == $S(" abc  abc ")
                @test isa(REPEAT(s, 2), $S)

                @test (r = uppercase($S("abc")); (r, typeof(r))) == ($S("ABC"), $S)
                @test (r = lowercase($S("ABC")); (r, typeof(r))) == ($S("abc"), $S)

                s = $S(" abc  abc ")
                $S==String || @test split(s, $S("c"))    == [$S(" ab"), $S("  ab"), $S(" ")]
                $S==String || @test split(s, r"c")       == [$S(" ab"), $S("  ab"), $S(" ")]
                @test occursin($S("c"), s) == true
                @test occursin($S("d"), s) == false
                @test (m = match(r"(c)", s); length(m.captures)) == 1
                @test (m = match(r"(d)", s); m===nothing)        == true

                @test INDEX($S("12341234"), $("23")           ) == 2
                @test INDEX($S("12341234"), $("23"), true     ) == 6
                @test INDEX($S("12341234"), $("23"), back=true) == 6
                @test INDEX($S("12341234"), $("00")           ) == 0
                @test INDEX($S("12341234"), $("00"), true     ) == 0

                @test SCAN($S("12341234"), $("23")            ) == 2
                @test SCAN($S("12341234"), $("23"), true      ) == 7
                @test SCAN($S("12341234"), $("23"), back=true ) == 7
                @test SCAN($S("12341234"), $("ab")            ) == 0
                @test SCAN($S("12341234"), $("ab"), true      ) == 0

            end
        end

        # check trim spaces but not '\0'
        s = F8"abc"
        @test (s[3] = 0x20; s) == F8"ab "
        @test LEN(s)           == 3
        @test LEN_TRIM(s)      == 2
        @test (s[3] = 0   ; s) == FortranStrings.FortranString{UInt8}(UInt8[0x61, 0x62, 0x00])
        @test LEN(s)           == 3
        @test LEN_TRIM(s)      == 3

        @test LSAME(F8"TRI", F8"tre") == true
        @test LSAME(F8"TRI", F"tre")  == true
        @test LSAME(F"TRI", F8"tre")  == true
        @test LSAME("TRI", F8"tre")   == true
        @test LSAME(F8"TRI", "tre")   == true
        @test LSAME(F8"TRI", 't')     == true
        @test LSAME(F8"TRI", F"upp")  == false
        @test LSAME(F"TRI", F8"upp")  == false
        @test LSAME(F8"TRI", 'u')     == false

    end

end


text:
"""
    abstract type AbstractEntity

Entity of geometrical nature. Identifiable throught its `(dim,tag)` key.
"""
abstract type AbstractEntity end

"""
    struct ElementaryEntity <: AbstractEntity

The most basic representation of an `AbstractEntity`.

# Fields:
- `dim::UInt8`: the geometrical dimension of the entity (e.g. line has `dim=1`, surface has `dim=2`, etc)
- `tag::Int64`: an integer tag associated to the entity
- `boundary::Vector{ElementaryEntity}`: the entities of dimension `dim-1` forming the entity's boundary
"""
struct ElementaryEntity <: AbstractEntity
    dim::UInt8
    tag::Int64
    boundary::Vector{ElementaryEntity}
    function ElementaryEntity(d::Integer, t::Integer, boundary::Vector{ElementaryEntity})
        msg = "an elementary entities in the boundary has the wrong dimension"
        for b in boundary
            @assert geometric_dimension(b) == d-1 msg
        end
        ent = new(d, t, boundary)
        # every entity gets added to a global variable ENTITIES so that we can
        # ensure the (d,t) pair is a UUID for an entity, and to easily retrieve
        # different entities.
        _global_add_entity!(ent)
        return ent
    end
end

"""
    ElementaryEntity(dim,tag)

Construct an [`ElementaryEntity`](@ref) with an empty boundary .
"""
function ElementaryEntity(dim,tag)
    ElementaryEntity(dim,tag,ElementaryEntity[])
end
ElementaryEntity(dim) = ElementaryEntity(dim,_new_tag(dim))

geometric_dimension(ω::ElementaryEntity) = ω.dim

"""
    tag(ω::ElementaryEntity)

Return the unique tag (for a given dimension) of the elementary entity.
"""
key(ω::ElementaryEntity) = (geometric_dimension(ω), ω.tag)

"""
    boundary(ω::ElementaryEntity)

Return the vector of elementary entities making the boundary.
"""
boundary(ω::ElementaryEntity) = ω.boundary

"""
    ==(Ω1::ElementaryEntity,Ω2::ElementaryEntity)

Two elementary entities are considered equal if their `dim` and `tag` fields
match.

Notice that this implies `dim` and `tag` of an elementary entity should uniquely
define it, and therefore global variables like [`TAGS`](@ref) are needed to make
sure newly created `AbstractEntities` have a new `(dim,tag)` identifier.
"""
function Base.:(==)(Ω1::AbstractEntity, Ω2::AbstractEntity)
    d1,t1 = key(Ω1)
    d2,t2 = key(Ω2)
    d1 == d2  || return false
    abs(t1) == abs(t2) || return false
    # boundary(Ω1) == boundary(Ω2) || return false # this should not be needed
    return true
end

"""
    ParametricEntity{D,F}

A geometric entity given by an explicit `parametrization::F` with domain given
by a singleton type `D`.

This differs from an `ElementaryEntity` in that the underlying representation of
the geometric entity is available.
"""
struct ParametricEntity{D,F} <: AbstractEntity
    dim::UInt8
    tag::Int
    parametrization::F
    function ParametricEntity{D,F}(dim,tag,f::F) where {D,F}
        ent = new{D,F}(dim,tag,f)
        # every entity gets added to a global variable ENTITIES so that we can
        # ensure the (d,t) pair is a UUID for an entity, and to easily retrieve
        # different entities.
        _global_add_entity!(ent)
        return ent
    end
end

"""
    const ParametricCurve{F} = ParametricEntity{ReferenceLine,F}

A curve defined by the function `f::F`.
"""
const ParametricCurve{F} = ParametricEntity{ReferenceLine,F}

function Base.reverse(ent::ParametricEntity{D,F}) where {D,F}
    ParametricEntity{D,F}(ent.dim,-ent.tag,ent.parametrization)
end

function ParametricEntity{D,F}(f::F) where {D,F}
    d = geometric_dimension(D)
    t = _new_tag(d) # automatically generate a new (valid) tag
    ParametricEntity{D,F}(d,t,f)
end

function ParametricEntity(f,d::AbstractReferenceShape)
    F = typeof(f)
    D = typeof(d)
    return ParametricEntity{D,F}(f)
end
ParametricEntity{D}(f) where {D} = ParametricEntity(f,D())

domain(p::ParametricEntity{D}) where {D<:AbstractReferenceShape} = D()
parametrization(p::ParametricEntity) = p.parametrization

key(ent::ParametricEntity) = (ent.dim,ent.tag)

# FIXME: what should this be?
boundary(ent::ParametricEntity) = ElementaryEntity[]

function Base.eltype(p::ParametricEntity)
    # NOTE: this relies on the brittle promote_op
    d = domain(p)
    x = center(d)
    f = parametrization(p)
    T = Base.promote_op(f,typeof(x))
end

ambient_dimension(p::ParametricEntity)   = length(eltype(p))
geometric_dimension(p::ParametricEntity) = geometric_dimension(domain(p))

(par::ParametricEntity)(x) = par.parametrization(x)

jacobian(psurf::ParametricEntity,s::SVector)  = ForwardDiff.jacobian(psurf.parametrization,s::SVector)
jacobian(psurf::ParametricEntity,s)           = jacobian(psurf,SVector(s))

function normal(ent::AbstractEntity,u)
    dim,tag = key(ent)
    s = sign(tag)
    N    = ambient_dimension(ent)
    M    = geometric_dimension(ent)
    msg  = "computing the normal vector requires the entity to be of co-dimension one."
    @assert N-M == 1 msg
    if M == 1 # a line in 2d
        t = jacobian(ent,u)
        ν = SVector(t[2],-t[1])
        return s*ν/norm(ν)
    elseif M == 2 # a surface in 3d
        j  = jacobian(ent,u)
        t₁ = j[:,1]
        t₂ = j[:,2]
        ν  = s*cross(t₁,t₂)
        return ν/norm(ν)
    else
        notimplemented()
    end
end

"""
    line(a,b)

Create a straight line connecting points `a` and `b`. This returns an instance
of [`ParametricCurve`](@ref).
"""
function line(a::SVector,b::SVector)
    f = (u) -> a + u[1]*(b-a)
    ParametricCurve(f)
end
line(a,b) = line(SVector(a),SVector(b))

"""
    Point{N,T} <: AbstractEntity
"""
struct Point{N,T} <: AbstractEntity
    coords::SVector{N,T}
end

geometric_dimension(::Type{<:Point}) = 0
geometric_dimension(pt::Point) = geometric_dimension(typeof(pt))
ambient_dimension(::Point{N}) where N = N
ambient_dimension(::Type{<:Point{N}}) where N = N

geometric_dimension(::Type{<:SVector}) = 0
geometric_dimension(pt::SVector) = geometric_dimension(typeof(pt))
ambient_dimension(::SVector{N}) where N = N
ambient_dimension(::Type{<:SVector{N}}) where N = N

#####################################################################

# Variables and functions to globally keep track of entities

#####################################################################

"""
    const TAGS::OrderedDict{Int,Vector{Int}}

Global dictionary storing the used entity tags (the value) for a given dimension
(the key).
"""
const TAGS     = OrderedDict{Int,Vector{Int}}()

"""
    const ENTITIES

Global dictionary storing the used entity tags (the value) for a given dimension
(the key).
"""
const ENTITIES = OrderedDict{Tuple{Int,Int},AbstractEntity}()

function _global_add_entity!(ent::AbstractEntity)
    d,t = key(ent)
    _add_tag!(d,t) # add this tag to global list to make sure it is not used again
    msg = "overwriting ENTITIES: value in key ($d,$t) will be replaced"
    haskey(ENTITIES,(d,t)) && (@warn msg)
    ENTITIES[(d,t)] = ent
    return d,t
end

"""
    _new_tag(dim)

Generate a unique tag for an `AbstractEntity` of dimension `dim`.

The implementation consists of adding one to the maximum value of `TAGS[dim]`

# See also: [`TAGS`](@ref).

"""
function _new_tag(dim)
    if !haskey(TAGS,dim)
        return 1
    else
        tnew = maximum(TAGS[dim]) + 1
        return tnew
    end
end

function _add_tag!(dim,tag)
    if is_new_tag(dim,tag)
        # now add key
        if haskey(TAGS,dim)
            push!(TAGS[dim],tag)
        else
            TAGS[dim] = [tag,]
        end
    else
        # print warning but don't add duplicate tag
        msg  = "entity of dimension $dim and tag $tag already exists in TAGS.
        Creating a possibly duplicate entity."
        @warn msg
    end
    return TAGS
end

function is_new_tag(dim,tag)
    if haskey(TAGS,dim)
        existing_tags = TAGS[dim]
        if in(tag,existing_tags)
            return false
        end
    end
    return true
end

clear_tags!() = empty!(TAGS)
clear_entities!() = empty!(ENTITIES)
function clear!()
    clear_tags!()
    clear_entities!()
    nothing
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
#=
variables = Dict{Symbol, DataFrames.DataFrame}()
variables[:P__ThermalStandard] = DataFrames.DataFrame(
    :one => [1, 2, 3, 2, 1],
    :two => [3, 2, 1, 2, 3],
    :three => [1, 2, 3, 2, 1],
)
variables[:P__RenewableDispatch] = DataFrames.DataFrame(
    :one => [3, 2, 3, 2, 3],
    :two => [1, 2, 1, 2, 1],
    :three => [3, 2, 3, 2, 3],
)

parameters = Dict{Symbol, DataFrames.DataFrame}()
parameters[:parameter_P_FixedGeneration] = DataFrames.DataFrame(
    :one => [2, 2, 1, 2, 2],
    :two => [3, 4, 1, 2, 2],
    :three => [3, 2, 3, 1, 1],
)

parameters[:parameter_P_PowerLoad] = DataFrames.DataFrame(
    :one => [3, 1, 3, 2, 1],
    :two => [1, 2, 1, 1, 1],
    :three => [3, 3, 3, 2, 3],
)
optimizer_log = Dict()
objective_value = Dict()
dual_values = Dict{Symbol, Any}()
base_power = 100.0
right_now = round(Dates.now(), Dates.Hour)
timestamp =
    DataFrames.DataFrame(:Range => right_now:Dates.Hour(1):(right_now + Dates.Hour(4)))
res = PG.Results(
    base_power,
    variables,
    optimizer_log,
    objective_value,
    timestamp,
    dual_values,
    parameters,
)

generators = Dict("Coal" => [:one; :two], "Wind" => [:three])
=#

function run_test_sim(result_dir::String)
    sim_name = "results_sim"
    sim_path = joinpath(result_dir, sim_name)
    if ispath(sim_path)
        executions = tryparse.(Int, readdir(sim_path))
        sim = joinpath(sim_path, string(maximum(executions)))
        @info "Reading results from last execution" sim
    else
        mkpath(result_dir)
        GLPK_optimizer =
            optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF)
        c_sys5_hy_uc = build_system("c_sys5_hy_uc", add_reserves = true)
        c_sys5_hy_ed = build_system("c_sys5_hy_ed")

        branches = Dict()
        services = Dict(
            :ReserveDown => ServiceModel(VariableReserve{ReserveDown}, RangeReserve),
            :ReserveUp => ServiceModel(VariableReserve{ReserveUp}, RangeReserve),
        )
        devices = Dict(
            :Generators => DeviceModel(ThermalStandard, ThermalStandardUnitCommitment),
            :Ren => DeviceModel(RenewableDispatch, RenewableFullDispatch),
            :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
            :ILoads => DeviceModel(InterruptibleLoad, DispatchablePowerLoad),
            :HydroEnergyReservoir =>
                DeviceModel(HydroEnergyReservoir, HydroDispatchReservoirStorage),
            :HydroROR => DeviceModel(HydroDispatch, FixedOutput),
        )
        template_hydro_st_uc =
            OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)

        devices = Dict(
            :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
            :Ren => DeviceModel(RenewableDispatch, RenewableFullDispatch),
            :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
            :ILoads => DeviceModel(InterruptibleLoad, DispatchablePowerLoad),
            :HydroEnergyReservoir =>
                DeviceModel(HydroEnergyReservoir, HydroDispatchReservoirStorage),
            :HydroROR => DeviceModel(HydroDispatch, FixedOutput),
        )
        template_hydro_st_ed =
            OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, Dict())
        stages_definition = Dict(
            "UC" => Stage(
                GenericOpProblem,
                template_hydro_st_uc,
                c_sys5_hy_uc,
                GLPK_optimizer,
            ),
            "ED" => Stage(
                GenericOpProblem,
                template_hydro_st_ed,
                c_sys5_hy_ed,
                GLPK_optimizer,
                constraint_duals = [:CopperPlateBalance],
            ),
        )

        sequence_cache = SimulationSequence(
            step_resolution = Hour(24),
            order = Dict(1 => "UC", 2 => "ED"),
            feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
            horizons = Dict("UC" => 24, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(24), Consecutive()),
                "ED" => (Hour(1), Consecutive()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => SemiContinuousFF(
                    binary_source_stage = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
                ("ED", :devices, :HydroEnergyReservoir) => IntegralLimitFF(
                    variable_source_stage = PSI.ACTIVE_POWER,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            cache = Dict(
                ("UC",) => TimeStatusChange(PSY.ThermalStandard, PSI.ON),
                ("UC", "ED") => StoredEnergy(PSY.HydroEnergyReservoir, PSI.ENERGY),
            ),
            ini_cond_chronology = InterStageChronology(),
        )
        sim = Simulation(
            name = "results_sim",
            steps = 2,
            stages = stages_definition,
            stages_sequence = sequence_cache,
            simulation_folder = result_dir,
        )
        build!(sim)
        execute!(sim)
    end

    results = SimulationResults(sim)
    results_uc = get_stage_results(results, "UC")
    results_ed = get_stage_results(results, "ED")

    return results_uc, results_ed
end


text:
"""
This is an abstract type for each kind
of ITEAD Nextion display objects.
"""
abstract type AbstractNexObject end

"""
    PageID(pid)

Page identifier.
"""
struct PageID
    _val::UInt8
end
PageID() = PageID(0)

"""
    ComponentID(pid)

Component identifier.
"""
struct ComponentID
    _val::UInt8
end
ComponentID() = ComponentID(0)

"""
    Name(s)

Name of component.
"""
struct Name
    _s::String
end
String(name::Name) = name._s

"""
    NexID(nexSerial, pid, cid, name)

An identifier for a Nextion object.

It's composed of 4 fields:

- `_nexSerial`: Nextion Serial struct with COM port
- `_name`: an unique name
- `_pid`: Page ID
- `_cid`: Component ID
"""
struct NexID{T <: AbstractNexSerial}
    _nexSerial::T

    _name::Name
    _pid::PageID
    _cid::ComponentID
end

"""
    PageID(nid) -> PageID

Return `PageID` from `NexID` `nid`.
"""
function PageID(nid::NexID)::PageID
    nid._pid
end
String(pid::PageID) = string(pid._val)

"""
    PageID(nexSerial)

Return PageID of current page.

It send `sendme` command to Nextion and wait response.
"""
function PageID(nexSerial::NexSerial)::PageID
    send(nexSerial, "sendme")
    sp = nexSerial._serial
    timeout_ms = 1000
    r = my_readuntil(sp, v_uint8_eoc, timeout_ms)
    #@info r
    if length(r) != 5
        error("Nextion must return a 5 bytes length message")
    end
    if r[1] != Integer(Return.Code.CURRENT_PAGE_ID_HEAD)
        error("Nextion must return $(Return.Code.CURRENT_PAGE_ID_HEAD)")
    end
    PageID(r[2])
end

"""
    ComponentID(nid) -> ComponentID

Return `ComponentID` from `NexID` `nid`.
"""
function ComponentID(nid::NexID)::ComponentID
    nid._cid
end
String(cid::ComponentID) = string(cid._val)

"""
    Name(nid) -> Name

Return `Name` from `NexID` `nid`.
"""
function Name(nid::NexID)::Name
    nid._name
end

"""
    obj.visible = val

Display Nextion object `obj` when `val` is `true`.
Hide it when `val` is `false`.

    obj.text = text_value

Set text from string value contained in `text_value` to Nextion object `obj`.

    obj.value = value

Set numerical (or boolean) value from `value` to Nextion object `obj`.

    obj.font = new_font

Set font from `Font` struct contained in `new_font` to Nextion object `obj`.

"""
function _setcommonproperty!(obj::AbstractNexObject, property::Symbol, new_val)

    ab = abilities(obj)

    # IViewable
    if property == :visible && :viewable in ab
        obj.viewable.visible = new_val

    # IStringValued
    elseif property == :text && :stringvalued in ab
        obj.stringvalued.value = new_val

    # INumericalValued
    elseif property == :value && :numericalvalued in ab
        if obj.numericalvalued.rn === nothing
            obj.numericalvalued.value = new_val
        else
            new_val = Number(obj.numericalvalued.rn(new_val))
            obj.numericalvalued.value = new_val
        end

    # IBooleanValued
    elseif property == :value && :booleanvalued in ab
        obj.booleanvalued.value = new_val

    # IFontStyleable
    elseif property == :font && :fontstyleable in ab
        obj.fontstyleable.font = new_val

    # IColourable
    elseif property == :backcolor && :colourable in ab
        obj.colourable.backcolor = new_val

    elseif property == :forecolor && :colourable in ab
        obj.colourable.forecolor = new_val

    # "IPicturable"
    elseif property == :picture && :picturable in ab
        obj.picturable.picture = new_val

    # "IWidthable"
    elseif property == :width && :widthable in ab
        obj.widthable.width = new_val

    # "IHeigthable"
    elseif property == :height && :heightable in ab
        obj.heightable.height = new_val

    # setfield!
    else
        setfield!(obj, property, new_val)

    end
end

"""
    obj.visible

Return True if `obj` is visible; return False otherwise.

    obj.text

Return `text`` property value.

    obj.value

Return `value`` property (numerical valued or boolean valued).
"""
function _getcommonproperty(obj::AbstractNexObject, property::Symbol)
    ab = abilities(obj)

    # IViewable
    if property == :visible && :viewable in ab
        obj.viewable.visible

    # IStringValued
    elseif property == :text && :stringvalued in ab
        obj.stringvalued.value

    # INumericalValued
    elseif property == :value && :numericalvalued in ab
        obj.numericalvalued.value

    # IBooleanValued
    elseif property == :value && :booleanvalued in ab
        obj.booleanvalued.value

    # IFontStyleable
    elseif property == :alignment && :fontstyleable in ab
        obj.fontstyleable.alignment

    # IColourable
    elseif property == :backcolor && :colourable in ab
        obj.colourable.backcolor
    elseif property == :forecolor && :colourable in ab
        obj.colourable.forecolor

    # "IPicturable"
    elseif property == :picture && :picturable in ab
        obj.picturable.picture

    # "IWidthable"
    elseif property == :width && :widthable in ab
        obj.widthable.width

    # "IHeigthable"
    elseif property == :height && :heightable in ab
        obj.heightable.height

    # getfield
    else
        getfield(obj, property)

    end
end

"""
    NexSerial(nid) -> NexSerial

Return NexSerial from NexID `nid`
"""
function NexSerial(nid::NexID)
    nid._nexSerial
end

"""
    NexID(obj) -> NexID

Return `NexID` from Nextion object obj.
"""
function NexID(obj::AbstractNexObject)
    obj._nid
end

#function name(obj::AbstractNexObject)
#    @deprecate name(obj) get(obj, Name)
#    obj._nid._name
#end

#=
# obj need to be mutable

function set!(obj::AbstractNexObject, pid::PageID)
  obj._nid._pid = pid
end

function set!(obj::AbstractNexObject, cid::ComponentID)
  obj._nid._cid = cid
end

function rename!(obj::AbstractNexObject, name::String)
  obj._nid._name = name
end
=#


text:
# A port of IDONE to Julia
# Mostly to get an idea of the performance overhead
# of using Python instead.

using LinearAlgebra
using Distributions

struct IDONEModel
    # Dimensionality
    d :: Int64
    # Bounds
    lb :: Vector{Float64}
    ub :: Vector{Float64}
    # Regularization parameter
    reg :: Float64
    # Weights of ReLU basis functions
    # As IDONE is for integers, these weights as well as biases
    # are integers as well. Store as float anyways.
    W :: Array{Float64, 2}
    # Biases of ReLU basis functions
    b :: Array{Float64, 1}
    # Basis function weights.
    # Updated using Recursive Least Squares
    c :: Vector{Float64}
    Δc :: Vector{Float64}
    # Keep track of the last x update.
    last_x_upd :: Vector{Float64}
    # Recursive Least Squares covariance matrix
    P :: Array{Float64, 2}
end

function generateWb(d :: Int64, lb :: Vector{Float64}, ub :: Vector{Float64})
    @assert d == length(lb)
    @assert d == length(ub)

    # Model offset
    n_basis = 1
    # One variable dependent (incl correction for d = 1)
    n_basis += (2 * round(Int64, sum(ub .- lb)))
    # Subsequent variable dependent
    n_basis += (2 * round(Int64, sum(1 .+ (ub[1:end-1] .- lb[2:end]) .- (lb[2:end] .- ub[1:end-1]))))

    # Actually construct the matrix
    W = zeros(Float64, (n_basis, d))
    b = zeros(Float64, n_basis)

    # Model offset
    fni = 1
    b[fni] = 1
    fni += 1

    # One variable dependent
    for k in 1:d
        for i in lb[k]:ub[k]
            if i == lb[k]
                W[fni, k] = 1
                b[fni] = -i
                fni += 1
            elseif i == ub[k]
                W[fni, k] = -1
                b[fni] = i
                fni += 1
            else
                W[fni, k] = -1
                b[fni] = i
                W[fni+1, k] = 1
                b[fni+1] = -i
                fni += 2
            end
        end
    end

    # Two subsequent variables
    for k in 2:d
        for i in (lb[k] - ub[k - 1]):(ub[k] - lb[k - 1] + 1)
            if i == lb[k] - ub[k-1]
                W[fni, k] = 1
                W[fni, k - 1] = -1
                b[fni] = -i
                fni += 1
            elseif i == ub[k] - lb[k - 1] + 1
                W[fni, k] = -1
                W[fni, k - 1] = 1
                b[fni] = i
                fni += 1
            else
                W[fni, k] = -1
                W[fni, k - 1] = 1
                b[fni] = i
                W[fni+1, k] = 1
                W[fni+1, k - 1] = -1
                b[fni+1] = -i
                fni += 2
            end
        end
    end

    return W, b
end

function generate_model(d :: Int64, lb :: Vector{Float64}, ub :: Vector{Float64})
    W, b = generateWb(d, lb, ub)
    D = length(b)
    c = ones(D)
    Δc = zeros(D)
    # No initial model offset
    c[1] = 0
    # reg = 1e-3
    reg = 1e-8
    P = Matrix(I / reg, D, D)
    return IDONEModel(d, lb, ub, reg, W, b, c, Δc, zeros(d), P)
end

function update_model!(model :: IDONEModel, x :: Vector{Float64}, y :: Float64)
    Zx = Z(model, x)
    model.last_x_upd .= x
    PZx = model.P * Zx
    g = PZx ./ (1 + Zx' * PZx)
    model.P .= model.P .- g * PZx'
    model.Δc .= (y - Zx' * model.c) .* g
    # model.c .= model.c .+ (y .- Zx' * model.c) .* g
    model.c .+= model.Δc
end

# Model evaluation
function id_relu(x)
    return max(0, x)
end

function ∇id_relu(x)
    return (x > 0) + 0.5 * (x == 0)
end

function Z(model :: IDONEModel, x :: Vector{Float64})
    return id_relu.(model.W * x .+ model.b)
end

function ∇Z(model :: IDONEModel, x :: Vector{Float64})
    return ∇id_relu.(model.W * x .+ model.b)
end

function ∇Z!(g :: Vector{Float64}, model :: IDONEModel, x :: Vector{Int64})
    resize!(g, length(model.c))
    mul!(g, model.W, x)
    g .+= model.b
    g .= ∇id_relu.(g)
end

function predict(model :: IDONEModel, x :: Vector{Float64})
    y_pred = dot(model.c, Z(model, x))
    return y_pred
end

function predict_v_and_diff(model :: IDONEModel, x :: Vector{Float64})
    Zx = Z(model, x)
    y_pred = dot(model.c, Zx)
    y_pred_diff = dot(model.Δc, Zx)
    return y_pred, y_pred_diff
end

function ∇predict(model :: IDONEModel, x :: Vector{Float64})
    return model.c' * Diagonal(∇Z(model, x)) * model.W
end

function ∇predict!(g :: Vector{Float64}, model :: IDONEModel, x :: Vector{Float64})
    gz = ∇Z(model, x)
    gz .*= model.c
    # gz' * W = W' * gz
    # Because apparently the other way around is not recognized.
    # Maybe a bug in Julia?
    mul!(g, model.W', gz)
end

function predict_dist_m2(model :: IDONEModel, x :: Vector{Float64})
    y_pred, y_pred_diff = predict_v_and_diff(model, x)
    return Normal(y_pred, abs(y_pred_diff) / sqrt(2))
end

function predict_EI_fg!(g :: Union{Vector{Float64}, Nothing}, model :: IDONEModel, x :: Vector{Float64}, f_best :: Float64)
    # Mixed Zx, gradient.
    vs = model.W * x .+ model.b
    Zx = id_relu.(vs)
    gZx = ∇id_relu.(vs)

    # Predict y value
    y_pred = dot(model.c, Zx)
    # Predict y variance.
    y_pred_diff = dot(model.Δc, Zx)
    # Predict gradient for y value.
    if !isnothing(g)
        gZx .*= model.c
        g_y_pred = model.W' * gZx
        # And for its' corresponding variance.
        gZx .*= model.Δc
        g_y_pred_diff = model.W' * gZx
    end
    # Estimated distribution
    μ = y_pred

    if x == model.last_x_upd
        #  Special case. Define σ = 0.
        if !isnothing(g)
            g .= g_y_pred
        end
        return y_pred - f_best
    end

    σ = abs(y_pred_diff) / sqrt(2)
    # Update gradient for σ
    if !isnothing(g)
        g_y_pred_diff .*= sign(y_pred_diff)
        g_y_pred_diff .*= sqrt(2)
    end

    # Calculate Expected Improvement and its corresponding gradient.
    n01 = Normal()

    # Calculate expected improvement
    # Calculate difference between computed mu and f_best
    # update gradient g_y_pred to be the gradient of diff_f_best_μ
    diff_f_best_μ = f_best - μ

    if !isnothing(g)
        g_y_pred .*= -1
    end
    #
    diff_f_best_μ_div_sigma = diff_f_best_μ / σ

    # cdf_diff_f_best_μ = cdf(n01, diff_f_best_μ_div_sigma)
    # pdf_diff_f_best_μ = pdf(n01, diff_f_best_μ_div_sigma)
    # 1 - to flip (we want to minimize the value, not maximize.)
    # integral goes from (f_best − μ)/σ to ∞ instead of -∞ to (f_best − μ)/σ.
    cdf_diff_f_best_μ = 1 - cdf(n01, diff_f_best_μ_div_sigma)
    pdf_diff_f_best_μ = - pdf(n01, diff_f_best_μ_div_sigma)
    ei = diff_f_best_μ * cdf_diff_f_best_μ + σ * pdf_diff_f_best_μ

    if !isnothing(g)
        g .=
            # Grad of diff_f_best_μ * cdf_diff_f_best_μ
            g_y_pred .* cdf_diff_f_best_μ .+
            (diff_f_best_μ * pdf_diff_f_best_μ) .* g_y_pred ./ g_y_pred_diff .+
            # Grad of σ * pdf_diff_f_best_μ
            g_y_pred_diff .* pdf_diff_f_best_μ .+
            (σ * -diff_f_best_μ_div_sigma * pdf_diff_f_best_μ) .* g_y_pred ./ g_y_pred_diff
    end
    return ei#, g
end

##

function max_EI(μ, σ, f_best)
    if σ < 1e-8
        # Small sigma -> full certainty.
        # Which means μ will be sampled with 100% certainty:
        # eg. it is either a full on improvement (if μ - f_best > 0)
        # or no improvement at all.
        return max(μ - f_best, 0.)
    end
    u = Normal()
    f_best_min_μ = f_best - μ
    f_best_min_μ_div_σ = f_best_min_μ / σ
    return f_best_min_μ * cdf(u, f_best_min_μ_div_σ) + σ * pdf(u, f_best_min_μ_div_σ)
end

function max_EI(n :: Normal, f_best)
    return max_EI(n.μ, n.σ, f_best)
end

function min_EI(μ, σ, f_best)
    if σ < 1e-8
        # Small sigma -> full certainty.
        # Which means μ will be sampled with 100% certainty:
        # eg. it is either a full on improvement (if μ - f_best > 0)
        # or no improvement at all.
        return min(μ - f_best, 0.)
    end
    u = Normal()
    f_best_min_μ = f_best - μ
    f_best_min_μ_div_σ = f_best_min_μ / σ
    return f_best_min_μ * (1. - cdf(u, f_best_min_μ_div_σ)) - σ * pdf(u, f_best_min_μ_div_σ)
end

function min_EI(n :: Normal, f_best)
    return min_EI(n.μ, n.σ, f_best)
end

##
using NLSolversBase
using LineSearches
using Optim

struct IDONE
    model :: IDONEModel
    f :: Function
    best_x :: Vector{Float64}
    best_y :: Ref{Float64}
    y0 :: Float64
    max_iter :: Int64
    current_iter :: Ref{Int64}
end

function scale_f(y0 :: Float64, y :: Float64)
    if abs(y0) > 1e-8
        return (y - y0) / abs(y0)
    else
        return (y - y0)
    end
end

function inv_scale_f(y0 :: Float64, y :: Float64)
    if abs(y0) > 1e-8
        return y * abs(y0) + y0
    else
        return y + y0
    end
end

function optimize_surrogate(model :: IDONEModel, initial_x :: Vector{Float64}, f_best :: Float64)

    df = OnceDifferentiable(
        (x) -> predict(model, x),
        (g, x) -> ∇predict!(g, model, x),
        ones(model.d)
        )

    # function fg!(F, G, x)
    #     return predict_EI_fg!(G, model, x, f_best)
    # end

    # df = OnceDifferentiable(only_fg!(fg!), ones(model.d))

    # -- Optimize using Optim.jl
    # This seems to occasonally run into errors, likely related
    # to their implementation, as manual inspection indicates the
    # model is well behaved at the point of error.
    # Backtracking Linesearch has issues with NaN values as well.
    # Which seems to be caused by being used together with Fminbox.
    optimizer = Fminbox()
    # optimizer = Fminbox(ParticleSwarm())
    # optimizer = Fminbox(LBFGS(linesearch=LineSearches.BackTracking()))
    options = Optim.Options(iterations=50)

    surrogate_solution = nothing
    try
        # surrogate_solution = optimize(df, model.lb, model.ub, initial_x, optimizer, options)
        surrogate_solution = optimize(df, model.lb .- 0.2, model.ub .+ 0.2, initial_x, optimizer, options)
    catch e
        # println("Optimizer failed: $e")
    end
    retries = 0
    while (isnothing(surrogate_solution) || any(isnan.(surrogate_solution.minimizer)) || isnan(surrogate_solution.minimum)) && retries < 10
        try
            surrogate_solution = optimize(df, model.lb, model.ub, initial_x, optimizer, options)
        catch e
            # println("Optimizer failed: $e")
        end
        retries += 1
    end

    if retries >= 10
        return
    end
    # println("Iteration [$(idone.current_iter[])/$(idone.max_iter)] Model indicates $(inv_scale(idone.y0, surrogate_solution.minimum)) at $(surrogate_solution.minimizer).")
    next_x = round.(surrogate_solution.minimizer)
    pred_y = surrogate_solution.minimum

    return next_x, pred_y
end

function variation!(next_x :: Vector{Float64}, idone :: IDONE)
    # Explore if not the last evaluation.
    if idone.current_iter != idone.max_iter
        # Random fuzzy exploration.
        for j in 1:idone.model.d
            # 1/d probability of exploring.
            rand() > 1/idone.model.d && continue
            lbounded = next_x[j] == idone.model.lb[j]
            ubounded = next_x[j] == idone.model.lb[j]
            # Bounds only allow one value.
            lbounded && ubounded && continue
            # Exploration direction is fixed.
            lbounded && (next_x[j] += 1; continue)
            ubounded && (next_x[j] -= 1; continue)
            # Pick randomly otherwise
            next_x[j] += rand(Bool) * 2 - 1
        end
    end
end

function step!(idone :: IDONE)

    # Obtain next point according to model.
    # next_x, pred_y = optimize_surrogate(idone.model, idone.best_x)
    result = optimize_surrogate(idone.model, idone.best_x, idone.best_y[])
    if isnothing(result)
        return
    end

    next_x, pred_y = result

    # Clamp it to be certain.
    next_x .= clamp.(next_x, idone.model.lb, idone.model.ub)
    # Apply variation
    variation!(next_x, idone)

    # @assert all(idone.model.lb .<= next_x)
    # @assert all(next_x .<= idone.model.ub)
    next_x .= clamp.(next_x, idone.model.lb, idone.model.ub)

    # Evaluate function at suggested point next_x
    next_y = idone.f(next_x)
    scaled_next_y = scale_f(idone.y0[], next_y)

    # println("Iteration [$(idone.current_iter[])/$(idone.max_iter)] Scaled $(scaled_next_y) found $(pred_y) predicted")

    # println("Predicted before update: $(predict(idone.model, next_x))")
    update_model!(idone.model, next_x, scaled_next_y)
    # println("Predicted after update: $(predict(idone.model, next_x)), should be $(scaled_next_y)")

    idone.current_iter[] % 1 == 0 &&
        println("Iteration [$(idone.current_iter[])/$(idone.max_iter)] Found value $(next_y) at $(next_x)")

    # Update best
    if idone.best_y[] > next_y
        idone.best_x .= next_x
        idone.best_y[] = next_y
    end

    # Next iteration
    idone.current_iter[] += 1
end

function optimize_IDONE(f :: Function, d :: Int64, lb :: Vector{Float64}, ub :: Vector{Float64}, n_eval=20)
    model = generate_model(d, lb, ub)
    # Start in the middle
    best_x = round.(lb .+ (ub .- lb) ./ 2)
    # println("$(best_x)")
    y0 = f(best_x)
    best_y = y0
    # Note, scaled_y = (y - y0) / abs(y0)
    update_model!(model, best_x, 0.0)

    #
    id = IDONE(model, f, best_x, Ref(best_y), y0, n_eval-1, Ref(1))

    for _ in 1:n_eval-1
        step!(id)
    end

    return id.best_y[], id.best_x, id.model
end

rosen_wrong(x) = sum(100.0 .* (x[1:end] .- x[end:-1:1] .^ 2.0) .^ 2.0 .+ (1 .- x[end:-1:1]) .^ 2.0)
rosen(x) = sum(100.0 .* (x[2:end] .- x[1:end-1] .^ 2.0) .^ 2.0 .+ (1 .- x[1:end-1]) .^ 2.0)

function nqueens(x)
    n = length(x)
    if !(all(x .> 0) && all(x .<= n))
        error("Incorrect indices in $x")
    end
    queens_in_column = zeros(length(x))
    queens_in_primary_diagonal = zeros(2 * length(x) - 1)
    queens_in_perp_diagonal = zeros(2 * length(x) - 1)
    score = 0.
    for (i, p) in enumerate(round.(Int64, x))
        score += queens_in_column[p]
        queens_in_column[p] += 1
        Pdiag = p - i + n
        score += queens_in_primary_diagonal[Pdiag]
        queens_in_primary_diagonal[Pdiag] += 1
        Pdiag = i + p - 1
        score += queens_in_perp_diagonal[Pdiag]
        queens_in_perp_diagonal[Pdiag] += 1
    end
    return score
end

max_stars_repo_name: AlgTUDelft/ExpensiveOptimBenchmark

max_issues_repo_name: WholeG/ExpensiveOptimBenchmark

max_forks_repo_name: WholeG/ExpensiveOptimBenchmark


text:
# The API client calling functions hosted by the server.
# Runs client when invoked directly with "--runclnt" argument.
# Call `run_clnt` otherwise.
using JuliaWebAPI
using Logging
using ZMQ
using Base.Test
using Requests

Logging.configure(level=INFO)

const NCALLS = 100
const APIARGS = randperm(NCALLS*4)

function printresp(apiclnt, testname, resp)
    hresp = httpresponse(apiclnt.format, resp)
    println("\t$(testname): $(hresp)")
    println("\t\tdata: $(hresp.data)")
    println("\t\thdrs: $(hresp.headers)")
end

function run_clnt(fmt, tport)
    ctx = Context()
    apiclnt = APIInvoker(tport, fmt)

    println("testing httpresponse...")
    resp = apicall(apiclnt, "testfn1", 1, 2, narg1=3, narg2=4)
    printresp(apiclnt, "testfn1", resp)

    resp = apicall(apiclnt, "testfn2", 1, 2, narg1=3, narg2=4)
    printresp(apiclnt, "testfn1", resp)

    resp = apicall(apiclnt, "testbinary", 10)
    printresp(apiclnt, "testbinary", resp)

    tic()
    for idx in 1:100
        arg1,arg2,narg1,narg2 = APIARGS[(4*idx-3):(4*idx)]
        resp = apicall(apiclnt, "testfn1", arg1, arg2; narg1=narg1, narg2=narg2)
        @test fnresponse(apiclnt.format, resp)["data"] == (arg1 * narg1) + (arg2 * narg2)
    end
    t = toc();
    println("time for $NCALLS calls to testfn1: $t secs @ $(t/NCALLS) per call")

    tic()
    for idx in 1:100
        arg1,arg2,narg1,narg2 = APIARGS[(4*idx-3):(4*idx)]
        resp = apicall(apiclnt, "testfn2", arg1, arg2; narg1=narg1, narg2=narg2)
        @test fnresponse(apiclnt.format, resp) == (arg1 * narg1) + (arg2 * narg2)
    end
    t = toc();
    println("time for $NCALLS calls to testfn2: $t secs @ $(t/NCALLS) per call")

    tic()
    for idx in 1:100
        arrlen = APIARGS[idx]
        resp = apicall(apiclnt, "testbinary", arrlen)
        @test isa(fnresponse(apiclnt.format, resp), Array)
    end
    t = toc();
    println("time for $NCALLS calls to testbinary: $t secs @ $(t/NCALLS) per call")

    # Test Array invocation
    println("testing array invocation...")

    resp = apicall(apiclnt, "testArray", Float64[1.0 2.0; 3.0 4.0])
    @test fnresponse(apiclnt.format, resp) == 12

    close(ctx)
    close(tport)
end

function run_httpclnt()
    println("starting http rpc tests.")

    resp = get("http://localhost:8888/")
    @test resp.status == 404

    resp = get("http://localhost:8888/invalidapi")
    @test resp.status == 404

    resp = JSON.parse(readstring(get("http://localhost:8888/testfn1/1/2")))
    @test resp["code"] == 0
    @test resp["data"] == 5

    resp = JSON.parse(readstring(get("http://localhost:8888/testfn1/1/2"; data=Dict(:narg1=>3, :narg2=>4))))
    @test resp["code"] == 0
    @test resp["data"] == 11

    println("testing file upload...")
    filename = "a.txt"
    postdata = """------WebKitFormBoundaryIabcPsAlNKQmowCx\r\nContent-Disposition: form-data; name="filedata"; filename="a.txt"\r\nContent-Type: text/plain\r\n\r\nLorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum.\n\r\n------WebKitFormBoundaryIabcPsAlNKQmowCx\r\nContent-Disposition: form-data; name="filename"\r\n\r\na.txt\r\n------WebKitFormBoundaryIabcPsAlNKQmowCx--\r\n"""
    headers = Dict("Content-Type"=>"multipart/form-data; boundary=----WebKitFormBoundaryIabcPsAlNKQmowCx")
    resp = JSON.parse(readstring(post("http://localhost:8888/testFile"; headers=headers, data=postdata)))
    @test resp["code"] == 0
    @test resp["data"] == "5,446"

    println("testing preprocessor...")
    resp = get("http://localhost:8888/testfn1/1/2"; headers = Dict("juliawebapi"=>"404"))
    @test resp.status == 404
    println("finished http rpc tests.")
end

# run client if invoked with run flag
!isempty(ARGS) && (ARGS[1] == "--runclnt") && run_clnt(JuliaWebAPI.JSONMsgFormat(), JuliaWebAPI.ZMQTransport("127.0.0.1", 9999, ZMQ.REQ, false, ctx))
!isempty(ARGS) && (ARGS[1] == "--runhttpclnt") && run_httpclnt()

max_stars_repo_name: tanmaykm/JuliaWebAPIPlugins

max_issues_repo_name: tanmaykm/JuliaWebAPIPlugins

max_forks_repo_name: tanmaykm/JuliaWebAPIPlugins