text:
function testonicwb(m; home=joinpath("d:\\", "icwb2-data"))
        for source in ["pku", "msr", "cityu", "as"]
                test_file = joinpath(home, "testing", "$(source)_test.utf8")
                truth_file = joinpath(home, "gold", source == "as" ? "as_testing_gold.utf8" : "$(source)_test_gold.utf8")
                dict_file = joinpath(home, "gold", "$(source)_training_words.utf8")
                x = [line for line in eachline(test_file) if length(line) > 0]
                z = [split(line) for line in eachline(truth_file) if length(line) > 0]
                y = m(x)
                #=
                debug_file = joinpath(pathof(KongYiji), "..", "..", "test", "wrong.cases.$(source).2.txt")
                open(debug_file, "w") do io
                        for i in 1:length(truth)
                                o, s = output[i], truth[i]
                                ms = eachmatch(o, s)
                                for i = 1:length(ms) - 1
                                        lo, ls = ms[i] .+ 1
                                        ro, rs = ms[i + 1] .- 1
                                        if lo <= ro
                                                print(io, "S : "); for j in ls:rs print(io, s[j], '\t') end; println(io)
                                                print(io, "O : "); for j in lo:ro print(io, o[j], '\t') end; println(io)
                                                println(io, "-------------------------")
                                        end
                                end
                                println(io, "=================================")
                        end
                end
                =#
                dict = Set([word for word in eachline(dict_file) if length(word) > 0])
                println(KongYiji.IcwbScoreTable(dict, z, y, source))
                println("==================================")

        end
end


text:
# Type aliases
const SFloat{T<:AbstractFloat} = Union{T, Complex{T}}
const SMatrix{T<:SFloat} = Union{SparseMatrixCSC{T,Int}, Diagonal{T}, Symmetric{T}, Hermitian{T}, Matrix{T}}
const SVector{T<:SFloat} = Union{SparseVector{T,Int}, Vector{T}}

const AbstractVecOrTuple{T} = Union{AbstractVector{<:T},NTuple{N,T} where N}

# Base QuObject abstract types
abstract type QuObject end
abstract type QuVector <: QuObject end
abstract type QuMatrix <: QuObject end

"""
    Ket(x, dims=(length(x),))

Construct a ket state vector from the vector `x`. A vector of length `N` will by default be assumed to be an element of a single Hilbert space of dimension `N`. If the vector is an element of a tensor product of Hilbert spaces, the dimensions can be defined manually by passing a tuple of subspace dimensions `dims`. In that case, `prod(dims)` must equal `length(x)`. By default, the vector is stored in sparse format.

The `Ket` type has two fields, `data` and `dims`, which store the vector data and the subspace dimensions. A `Ket`, like a [`Bra`](@ref) or an [`Operator`](@ref) is parameterized by the number of subspaces it lives in. Two different kets must have the same system dimensions in order to be added together.

It is possible to normalize the ket vector after construction with the `normalize!` function.

# Example
```jldoctest
julia> ψ = normalize!(Ket([1,1]))
2-d Ket{Array{Float64,1},1} with dimensions 2
0.71∠0°|0⟩ + 0.71∠0°|1⟩
```
"""
struct Ket{T<:SVector,D} <: QuVector
    data::T
    dims::Dims{D}
    function Ket{T,D}(x, dims::Dims{D}) where {T<:SVector,D}
        checkdims(dims,size(x))
        return new{T,D}(convert(T,x),dims)
    end
end
Ket(x::T, dims::Dims{D}=size(x)) where {T<:SVector,D} = Ket{T,D}(x,dims)
Ket(x::AbstractVector,dims::Dims=size(x)) = Ket(float.(x),dims)

"""
    Bra(x, dims=(length(x),))

Bra vector type. The dual vector to the `Ket`.

The `Bra` type has two fields, `data` and `dims`, which store the vector data and the subspace dimensions. A `Bra`, like a [`Ket`](@ref) or an [`Operator`](@ref) is parameterized by the number of subspaces it lives in. Two different kets must have the same system dimensions in order to be added together.

It is possible to normalize the bra vector after construction with the `normalize!` function.
"""
struct Bra{T<:SVector,D} <: QuVector
    data::T
    dims::Dims{D}
    function Bra{T,D}(x,dims::Dims{D}) where {T<:SVector,D}
        checkdims(dims,size(x))
        return new{T,D}(convert(T,x), dims)
    end
end
Bra(x::T, dims::Dims{D}=size(x)) where {T<:SVector,D} = Bra{T,D}(x,dims)
Bra(x::AbstractVector,dims::Dims=size(x)) = Bra(float.(x),dims)

"""
    Operator(B, dims=(size(B,1),))

Construct a linear operator from the matrix `B`. An `N`×`N` matrix will by default be assumed to describe an operator that acts on a single Hilbert space of dimension `N`. If the matrix represents a linear operator on a tensor product of Hilbert spaces, the dimensions can be defined manually by passing a tuple of subspace dimensions `dims`. In that case, `prod(dims)` must equal `size(B,1)`.

The `Operator` type has two fields, `data` and `dims`, which store the matrix data and the subspace dimensions. An `Operator`, like a [`Ket`](@ref) or a [`Bra`](@ref), is parameterized by the number of subspaces it lives in. Two different density matrices must have the same system dimensions in order to be added together. An `Operator` may or may not be Hermitian.

# Example
```jldoctest
julia> σ = Operator([0 -im ; im 0])
2×2 Operator{Array{Complex{Float64},2},1} with dimensions 2
 0.0+0.0im  0.0-1.0im
 0.0+1.0im  0.0+0.0im
```
"""
struct Operator{T<:SMatrix,D} <: QuMatrix
    data::T
    dims::Dims{D}
    function Operator{T,D}(B,dims::Dims{D}) where {T<:SMatrix,D}
        checkdims(dims,size(B))
        return new{T,D}(convert(T,B), dims)
    end
end
Operator(B::T,dims::Dims{D}=(size(B,1),)) where {T<:SMatrix,D} = Operator{T,D}(B,dims)
Operator(B::AbstractMatrix,dims::Dims=(size(B,1),)) = Operator(float.(B),dims)

# Conversion between different QuObjects
Bra(x::Ket) = Bra(conj(x.data),x.dims)
Ket(x::Bra) = Ket(conj(x.data),x.dims)
Operator(x::Ket) = Operator(x.data*x.data',x.dims)
Operator(x::Bra) = Operator(conj(x.data)*tranpose(x.data),x.dims)

# function to check if the dims given are consistent with the object size
function checkdims(dims::Dims,size::NTuple{1,Integer})
    if prod(dims) != (m=size[1])
        throw(DimensionMismatch("subspace dimensions $dims are not consistent with a vector of length $m"))
    end
end
function checkdims(dims::Dims,size::NTuple{2,Integer})
    if size[2] != (m=size[1])
        throw(DimensionMismatch("matrix is not square, size is $m×$(size[2])"))
    end
    if prod(dims) != m
        throw(DimensionMismatch("subspace dimensions $dims are not consistent with a matrix of size $m×$m"))
    end
end

# Density type remnants
struct Density{T<:SMatrix,D} <: QuMatrix
    data::T
    dims::Dims{D}
    function Density{T,D}(A::T,dims::Dims{D}) where {T<:SMatrix,D}
        checkdims(dims,size(A))
        isapproxhermitian(A) || throw(ArgumentError("a density matrix must be Hermitian"))
        return new{T,D}(A, dims)
    end
end

ext: jl


text:
using OrderedCollections: OrderedDict
using PropertyDicts
using Test

@testset "PropertyDicts" begin
    d = Dict("foo"=>1, :bar=>2)
    _keys = collect(keys(d))
    pd = PropertyDict(d)

    str_props = PropertyDict("foo" => 1, "bar" => 2)
    sym_props = PropertyDict(:foo => 1, :bar => 2)

    nt = (d =1, )
    pnt = pairs(nt)
    ntpd = PropertyDict(pnt)
    @test keys(PropertyDict(ntpd)) === keys(nt)
    @test values(PropertyDict(ntpd)) === values(pnt)
    @test propertynames(PropertyDict(ntpd)) === propertynames(nt)
    @test empty!(PropertyDict(Dict("foo"=>1, :bar=>2))) isa PropertyDict

    if isdefined(Base, :hasproperty)
        @test hasproperty(sym_props, "bar")
        @test hasproperty(str_props, :bar)
        @test hasproperty(pd, :foo)
        @test hasproperty(pd, "bar")
    end

    @testset "convert" begin
        expected = OrderedDict
        result = convert(expected, pd)

        @test result isa expected
    end

    @testset "get" begin
        @testset "default value" begin
            default = "baz"

            @test get(pd, "DNE", default) == default
            @test get(() -> 3, str_props, :foo) == 1
            @test get(() -> 3, str_props, :baz) == 3
            @test get(str_props, :baz, 3) == 3
            @test get(sym_props, "baz", 3) == 3
            @test get!(str_props, :baz, 3) == 3
            @test get!(sym_props, "baz", 3) == 3
            @test get!(() -> 4, str_props, :baz) == 3
            @test get!(() -> 4, sym_props, "baz") == 3
            @test get!(() -> 4, sym_props, "buz") == 4
        end

        @testset "$(typeof(key))" for key in _keys
            @test get(pd, key, "DNE") == d[key]
        end
    end

    @testset "getindex - $key" for key in _keys
        @test getindex(pd, key) == getindex(d, key)
    end

    @testset "getproperty" begin
        @test pd.foo == 1
        @test pd.bar == 2
        sym_props."spam" = 4
        @test sym_props."spam" == 4

        str_props.spam = 4
        @test str_props.spam == 4
    end

    @testset "iterate" begin
        @test iterate(pd) == iterate(d)
        @test iterate(pd, 1) == iterate(d, 1)
        @test iterate(pd, 2) == iterate(d, 2)
    end

    @testset "iteratorsize" begin
        @test Base.IteratorSize(pd) == Base.IteratorSize(d)
    end

    @testset "iteratoreltype" begin
        @test Base.IteratorEltype(pd) == Base.IteratorEltype(d)
    end

    @testset "length" begin
        @test length(pd) == length(d)
    end

    @testset "string" begin
        @test string(pd) == string(d)
    end

    push!(pd, :buz => 10)
    @test pop!(pd, :buz, 20) == 10
    @test pop!(pd, :buz, 20) == 20
    @test sizehint!(pd, 5) === pd
    @test get(pd, delete!(pd, "foo"), 10) == 10
end


text:
PositiveFactorizations.floattype(::Type{Union{T,Missing}}) where {T} = float(T)

function newton_raphson(fgh!, x::AbstractArray{T}; ρ=one(T)/2, c=1e-4, tol=1e-4, max_iter=1000) where T
    n = length(x)
    R = Base.promote_op(/, T, T)
    grad = zeros(R, n)
    hes = zeros(R, n, n)
    y = zero(R)
    for i = 1:max_iter
        y = fgh!(x, grad, hes, true)
        search_dir = -(cholesky(Positive, hes) \ grad)
        norm(search_dir) > tol || return x, y, grad, hes
        step_size = one(T)
        while fgh!(x + search_dir * step_size, grad, hes, false) > y + c * step_size * grad'search_dir
            step_size *= ρ
            step_size > 1e-10 || error("Line search failed! Problematic Hessian or Gradient?")
        end
        x .+= search_dir * step_size
    end
    throw(ConvergenceException(max_iter))
end


text:
module Chain

export @chain

is_aside(x) = false
is_aside(x::Expr) = x.head == :macrocall && x.args[1] == Symbol("@aside")

insert_first_arg(symbol::Symbol, firstarg) = Expr(:call, symbol, firstarg)
insert_first_arg(any, firstarg) = error("Can't insert an argument to $any. Needs to be a Symbol or a call expression")

function insert_first_arg(e::Expr, firstarg)

    # f(a, b) --> f(firstarg, a, b)
    if e.head == :call && length(e.args) > 1
        Expr(e.head, e.args[1], firstarg, e.args[2:end]...)

    # @. somesymbol --> somesymbol.(firstarg)
    elseif e.head == :macrocall && length(e.args) == 3 && e.args[1] == Symbol("@__dot__") &&
            e.args[2] isa LineNumberNode && e.args[3] isa Symbol
        Expr(:., e.args[3], Expr(:tuple, firstarg))

    # @macro(a, b) --> @macro(firstarg, a, b)
    elseif e.head == :macrocall && e.args[1] isa Symbol && e.args[2] isa LineNumberNode
        Expr(e.head, e.args[1], e.args[2], firstarg, e.args[3:end]...)

    else
        error("Can't prepend first arg to expression $e that isn't a call.")
    end
end

function rewrite(expr, replacement)
    aside = is_aside(expr)
    if aside
        length(expr.args) != 3 && error("Malformed @aside macro")
        expr = expr.args[3] # 1 is macro symbol, 2 is LineNumberNode
    end

    had_underscore, new_expr = replace_underscores(expr, replacement)

    if !aside
        if !had_underscore
            new_expr = insert_first_arg(new_expr, replacement)
        end
        replacement = gensym()
        new_expr = Expr(Symbol("="), replacement, new_expr)
    end

    (new_expr, replacement)
end

rewrite(l::LineNumberNode, replacement) = (l, replacement)

function rewrite_chain_block(firstpart, block)
    if !(block isa Expr && block.head == :block)
        error("Second argument of @chain must be a begin / end block")
    end

    block_expressions = block.args
    isempty(block_expressions) && error("No expressions found in chain block.")

    rewritten_exprs = []
    replacement = firstpart

    for expr in block_expressions
        rewritten, replacement = rewrite(expr, replacement)
        push!(rewritten_exprs, rewritten)
    end

    result = Expr(:let, Expr(:block), Expr(:block, rewritten_exprs...))

    :($(esc(result)))
end

"""
    @chain(initial_value, block::Expr)

Rewrites a block expression to feed the result of each line into the next one.
The initial value is given by the first argument.

In all lines, underscores are replaced by the previous line's result.
If there are no underscores and the expression is a symbol, the symbol is rewritten
to a function call with the previous result as the only argument.
If there are no underscores and the expression is a function call or a macrocall,
the call has the previous result prepended as the first argument.

Example:

```
x = @chain [1, 2, 3] begin
    filter(!=(2), _)
    sqrt.(_)
    sum
end
x == sum(sqrt.(filter(!=(2), [1, 2, 3])))
```
"""
macro chain(initial_value, block::Expr)
    rewrite_chain_block(initial_value, block)
end

function rewrite_chain_block(block)
    if !(block isa Expr && block.head == :block)
        error("Only argument of single-argument @chain must be a begin / end block")
    end

    block_expressions = block.args
    isempty(block_expressions) && error("No expressions found in chain block.")

    # assign first line to first gensym variable
    firstvar = gensym()
    rewritten_exprs = []
    replacement = firstvar

    did_first = false
    for expr in block_expressions
        # could be and expression first or a LineNumberNode, so a bit convoluted
        # we just do the firstvar transformation for the first non LineNumberNode
        # we encounter
        if !(did_first || expr isa LineNumberNode)
            expr = Expr(Symbol("="), firstvar, expr)
            did_first = true
            push!(rewritten_exprs, expr)
            continue
        end

        rewritten, replacement = rewrite(expr, replacement)
        push!(rewritten_exprs, rewritten)
    end

    result = Expr(:let, Expr(:block), Expr(:block, rewritten_exprs...))

    :($(esc(result)))
end

"""
    @chain(block::Expr)

Rewrites a block expression to feed the result of each line into the next one.
The first line serves as the initial value and is not rewritten.

In all other lines, underscores are replaced by the previous line's result.
If there are no underscores and the expression is a symbol, the symbol is rewritten
to a function call with the previous result as the only argument.
If there are no underscores and the expression is a function call or a macrocall,
the call has the previous result prepended as the first argument.

Example:

```
x = @chain begin
    [1, 2, 3]
    filter(!=(2), _)
    sqrt.(_)
    sum
end
x == sum(sqrt.(filter(!=(2), [1, 2, 3])))
```
"""
macro chain(block::Expr)
    rewrite_chain_block(block)
end

function replace_underscores(expr::Expr, replacement)
    found_underscore = false

    # if a @chain macrocall is found, only its first arg can be replaced if it's an
    # underscore, otherwise the macro insides are left untouched
    if expr.head == :macrocall && expr.args[1] == Symbol("@chain")
        length(expr.args) != 4 && error("Malformed nested @chain macro")
        expr.args[2] isa LineNumberNode || error("Malformed nested @chain macro")
        arg3 = if expr.args[3] == Symbol("_")
            found_underscore = true
            replacement
        else
            expr.args[3]
        end
        newexpr = Expr(:macrocall, Symbol("@chain"), expr.args[2], arg3, expr.args[4])
    # for all other expressions, their arguments are checked for underscores recursively
    # and replaced if any are found
    else
        newargs = map(x -> replace_underscores(x, replacement), expr.args)
        found_underscore = any(first.(newargs))
        newexpr = Expr(expr.head, last.(newargs)...)
    end
    return found_underscore, newexpr
end

function replace_underscores(x, replacement)
    if x == Symbol("_")
        true, replacement
    else
        false, x
    end
end

end


text:
using FDFDSolver
using Documenter

DocMeta.setdocmeta!(FDFDSolver, :DocTestSetup, :(using FDFDSolver); recursive=true)

makedocs(;
    modules=[FDFDSolver],
    authors="Mohamed Kamal AbdElrahman",
    repo="https://github.com/MKAbdElrahman/FDFDSolver.jl/blob/{commit}{path}#{line}",
    sitename="FDFDSolver.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://MKAbdElrahman.github.io/FDFDSolver.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/MKAbdElrahman/FDFDSolver.jl",
    devbranch="main",
)


text:
using DataFrames
using CSV

cd("/Users/kiristern/Documents/GitHub/SDM/")

predators = CSV.read("data/preds/originals/Canis_latrans.csv")
# # create new column "new_sp" for subspecies naming
# predators.new_sp = copy(predators.species)
# for i in 1:length(df)
#     if predators.taxonRank[i] == "SUBSPECIES"
#         predators.new_sp[i] = string(predators.species[i], "_", df.infraspecificEpithet[i])
#     end
# end

predators = select(predators, [:species, :decimalLatitude, :decimalLongitude, :year])

# select only ursus arctos from imported all predators csv file
#df = predators[predators.species .== "Ursus arctos", :]

p1 = CSV.read("data/prey_sp/Thamnophis_atratus.csv", delim="\t")
p1 = select(p1,[:species, :decimalLatitude, :decimalLongitude, :year])

p2 = CSV.read("data/prey_sp/Taricha_torosa.csv", delim="\t")
p2 = select(p2,[:species, :decimalLatitude, :decimalLongitude, :year])

p3 = CSV.read("data/prey_sp/Rana_draytonii.csv", delim="\t")
p3 = select(p3,[:species, :decimalLatitude, :decimalLongitude, :year])

p4 = CSV.read("data/prey_sp/Spermophilus_beecheyi.csv", delim="\t")
p4 = select(p4,[:species, :decimalLatitude, :decimalLongitude, :year])

p5 = CSV.read("data/prey_sp/Anaxyrus_boreas.csv", delim="\t")
p5 = select(p5,[:species, :decimalLatitude, :decimalLongitude, :year])

p6 = CSV.read("data/prey_sp/Fulica_americana.csv", delim="\t")
p6 = select(p6,[:species, :decimalLatitude, :decimalLongitude, :year])

p7 = CSV.read("data/prey_sp/Otospermophilus_beecheyi.csv", delim="\t")
p7 = select(p7,[:species, :decimalLatitude, :decimalLongitude, :year])

p8 = CSV.read("data/prey_sp/Odocoileus_virginianus.csv", delim="\t")
p8 = select(p8,[:species, :decimalLatitude, :decimalLongitude, :year])

p9 = CSV.read("data/prey_sp/Thomomys_bottae.csv", delim="\t")
p9 = select(p9,[:species, :decimalLatitude, :decimalLongitude, :year])

# p10 = CSV.read("data/prey_sp/Boiga_irregularis.csv", delim="\t")
# p10 = select(p10,[:species, :decimalLatitude, :decimalLongitude, :year])
#
# p11 = CSV.read("data/prey_sp/Trichosurus_vulpecula.csv", delim="\t")
# p11 = select(p11,[:species, :decimalLatitude, :decimalLongitude, :year])
#
# p12 = CSV.read("data/prey_sp/Tachyglossus_aculeatus.csv", delim="\t")
# p12 = select(p12,[:species, :decimalLatitude, :decimalLongitude, :year])

#create one dataframe with predators/prey dataframes
pred_prey = vcat(predators, p1, p2, p3, p4, p5,p6,p7,p8,p9)

pred_prey = dropmissing!(pred_prey, [:decimalLatitude, :decimalLongitude, :year])

# rename!(pred_prey, :decimalLatitude => :latitude)
# rename!(pred_prey, :decimalLongitude => :longitude)

CSV.write("data/pred_prey/Canis_latrans.csv", pred_prey)

#select for one species from dataframe
species = [pred_prey[pred_prey.species .== u,:] for u in unique(pred_prey.species)]
occ = species[1]

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
"""
    timeout(f::Function, seconds::Real) -> Union{Some,Nothing}

Executes the provided function `f` and returns the result wrapped in a `Some`. If the given
function takes longer than `seconds` the function is terminated and `nothing` is returned.

Note that the timeout will only work correctly if the called function `f` yields which will
occur automatically if performing an I/O operation.
"""
function timeout(f::Function, seconds::Real)
    result = nothing
    c = Condition()

    try
        # Execute the given function as a task so we can interrupt it if necessary
        t = @async begin
            result = Some(f())
            notify(c)
        end

        # Create a timeout task which aborts the wait early if we hit the timeout.
        # Note: We use an async task here to ensure when `f` finishes before the timeout it
        # is not blocked by the timer.
        @async begin
            # Note: The `pollint` will may cause abort events to occur after the specified
            # timeout. e.g. if `seconds = 1.5` then the abort would occur at 2 seconds.
            timedwait(float(seconds); pollint=1.0) do
                istaskdone(t)
            end
            notify(c)
        end

        wait(c)

        # Kill the function task if it is still executing
        istaskdone(t) || @async Base.throwto(t, InterruptException())
        wait(t)
    catch e
        # Unwrap the TaskFailedException
        if VERSION >= v"1.3.0-alpha.110" && e isa TaskFailedException
            e = e.task.exception
        end

        # Ignore the kill exception
        e isa InterruptException || rethrow(e)
    end

    return result
end


text:
module RedyFlavoured
using MLStyle.AbstractPatterns
using MLStyle.Err: PatternCompilationError

Config = NamedTuple{(:type, :ln)}
Scope = ChainDict{Symbol, Symbol}
ViewCache = ChainDict{Pair{TypeObject, Any}, Tuple{Symbol, Bool}}
Terminal = Dict{Int, Any}
function update_parent!(view_cache::ViewCache)
    parent = view_cache.init[]
    for (typed_viewer, (sym, _)) in view_cache.cur
        parent[typed_viewer] = (sym, false)
    end
end

struct CompileEnv
    # Dict(user_defined_name => actual_name). mangling for scope safety
    scope::Scope
    # Dict(view => (viewed_cache_symbol => guarantee_of_defined?))
    view_cache::ViewCache
    terminal::Dict{Int, Any}
    hygienic::Bool
    ret::Symbol
    final::Symbol
end

function CompileEnv(terminal::Terminal, hygienic::Bool, ret::Symbol, final::Symbol)
    CompileEnv(Scope(), ViewCache(), terminal, hygienic, ret, final)
end

function (env::CompileEnv)(;
    scope::Union{Nothing, Scope} = nothing,
    view_cache::Union{Nothing, ViewCache} = nothing,
)
    scope === nothing && (scope = env.scope)
    view_cache === nothing && (view_cache = env.view_cache)
    CompileEnv(scope, view_cache, env.terminal, env.hygienic, env.ret, env.final)
end

abstract type Cond end
struct AndCond <: Cond
    left::Cond
    right::Cond
end

struct OrCond <: Cond
    left::Cond
    right::Cond
end

struct TrueCond <: Cond
    stmt::Any
end

TrueCond() = TrueCond(true)

struct CheckCond <: Cond
    expr::Any
end

"""
build to disjunctive forms
"""
function build_readable_expression!(
    exprs::Vector{Any},
    following_stmts::Vector{Any},
    cond::CheckCond,
)
    expr = cond.expr
    if !isempty(following_stmts)
        expr = Expr(:block, following_stmts..., expr)
        empty!(following_stmts)
    end
    push!(exprs, expr)
end

function build_readable_expression!(
    exprs::Vector{Any},
    following_stmts::Vector{Any},
    cond::TrueCond,
)
    cond.stmt isa Union{Bool, Int, Float64, Nothing} && return #= shall contain more literal typs =#
    push!(following_stmts, cond.stmt)
    nothing
end

function build_readable_expression!(
    exprs::Vector{Any},
    following_stmts::Vector{Any},
    cond::AndCond,
)
    build_readable_expression!(exprs, following_stmts, cond.left)
    build_readable_expression!(exprs, following_stmts, cond.right)
end

function build_readable_expression!(
    exprs::Vector{Any},
    following_stmts::Vector{Any},
    cond::OrCond,
)
    exprs′ = []
    following_stmts′ = []
    build_readable_expression!(exprs′, following_stmts′, cond.left)
    left = to_expression(exprs′, following_stmts′)

    empty!(exprs′)
    empty!(following_stmts′)

    build_readable_expression!(exprs′, following_stmts′, cond.right)
    right = to_expression(exprs′, following_stmts′)

    empty!(exprs′)
    empty!(following_stmts′)

    bool_or = Expr(:||, left, right)
    if !isempty(following_stmts)
        bool_or = Expr(:block, following_stmts..., bool_or)
        empty!(following_stmts)
    end

    push!(exprs, bool_or)
end

function to_expression(cond::Cond)
    exprs = []
    following_stmts = []
    build_readable_expression!(exprs, following_stmts, cond)
    to_expression(exprs, following_stmts)
end

function to_expression(exprs::Vector{Any}, following_stmts::Vector)
    bool_and(a, b) = Expr(:&&, a, b)
    if isempty(following_stmts)
        isempty(exprs) && return true
        foldr(bool_and, exprs)
    else
        init = Expr(:block, following_stmts..., true)
        foldr(bool_and, exprs, init = init)
    end
end

allsame(xs::Vector) = any(e -> e == xs[1], xs)

const CACHE_NO_CACHE = 0
const CACHE_MAY_CACHE = 1
const CACHE_CACHED = 2

function static_memo(
    f::Function,
    view_cache::ViewCache,
    op::APP;
    ty::TypeObject,
    depend::Union{Nothing, APP} = nothing,
)
    if op isa NoPre
        nothing
    elseif op isa NoncachablePre
        cached_sym = nothing
        f(cached_sym, CACHE_NO_CACHE)
    else
        cache_key = depend === nothing ? op : (depend => op)
        cache_key = Pair{TypeObject, Any}(ty, cache_key)
        cached = get(view_cache, cache_key) do
            nothing
        end::Union{Tuple{Symbol, Bool}, Nothing}
        if cached === nothing
            cached_sym = gensym("cache")
            computed_guarantee = false
        else
            (cached_sym, computed_guarantee) = cached
        end
        if !computed_guarantee
            f(cached_sym, CACHE_MAY_CACHE)
            view_cache[cache_key] = (cached_sym, true)
            cached_sym
        else
            f(cached_sym, CACHE_CACHED)
        end
    end
end

function init_cache(view_cache::ViewCache)
    block = Expr(:block)
    cache_syms = block.args
    # TODO: OPT
    for_chaindict_dup(view_cache) do _, (view_cache_sym, _)
        push!(cache_syms, :($view_cache_sym = nothing))
    end
    if isempty(cache_syms)
        true
    else
        block
    end
end

function myimpl()
    function cache(f)
        function apply(env::CompileEnv, target::Target{true})::Cond
            target′ = target.with_repr(gensym(), Val(false))
            AndCond(TrueCond(:($(target′.repr) = $(target.repr))), f(env, target′))
        end
        function apply(env::CompileEnv, target::Target{false})::Cond
            f(env, target)
        end
        apply
    end

    wildcard(::Config) = (::CompileEnv, target::Target) -> TrueCond()

    literal(v, config::Config) = function ap_literal(::CompileEnv, target::Target)::Cond
        ty = typeof(v)
        if v isa Symbol
            v = QuoteNode(v)
        end
        (isprimitivetype(ty) || ty.size == 0 && !ty.mutable) ?
        CheckCond(:($(target.repr) === $v)) : CheckCond(:($(target.repr) == $v))
    end

    function and(ps::Vector{<:Function}, config::Config)
        @assert !isempty(ps)
        function ap_and_head(env::CompileEnv, target::Target{false})::Cond
            hd = ps[1]::Function
            tl = view(ps, 2:length(ps))
            init = hd(env, target)

            # the first conjuct must be executed, so the computation can get cached:
            # e.g.,
            #   match val with
            #   | View1 && View2 ->
            # and we know `View1` must be cached.
            (computed_guarantee′, env′, ret) = foldl(
                tl,
                init = (true, env, init),
            ) do (computed_guarantee, env, last), p
                # `TrueCond` means the return expression must be evaluated to `true`
                computed_guarantee′ = computed_guarantee && last isa TrueCond
                if !computed_guarantee′ && computed_guarantee
                    view_cache = env.view_cache
                    view_cache′ = child(view_cache)
                    view_cache_change = view_cache′.cur
                    env = env(view_cache = view_cache′)
                end
                computed_guarantee′, env, AndCond(last, p(env, target))
            end

            if !computed_guarantee′
                update_parent!(env′.view_cache)
            end
            ret
        end |> cache
    end

    function or(ps::Vector{<:Function}, config::Config)
        @assert !isempty(ps)
        function ap_or(env::CompileEnv, target::Target{false})::Cond
            or_checks = Cond[]
            scope = env.scope
            view_cache = env.view_cache
            scopes = Dict{Symbol, Symbol}[]
            n_ps = length(ps)
            for p in ps
                scope′ = child(scope)
                env′ = env(; scope = scope′, view_cache = view_cache)
                push!(or_checks, p(env′, target.clone))
                push!(scopes, scope′.cur)
            end

            # check the change of scope discrepancies for all branches
            intersected_new_names = reduce(
                intersect!,
                (keys(scope) for scope in scopes[2:n_ps]),
                init = Set(keys(scopes[1])),
            )

            for key in intersected_new_names
                refresh = gensym(key)
                for i in eachindex(or_checks)
                    check = or_checks[i]
                    old_name = scopes[i][key]
                    or_checks[i] =
                        AndCond(or_checks[i], TrueCond(:($refresh = $old_name)))
                end
                scope[key] = refresh
            end
            foldr(OrCond, or_checks)
        end |> cache
    end

    # C(p1, p2, .., pn)
    # pattern = (target: code, remainder: code) -> code
    function decons(comp::PComp, extract::Function, ps::Vector, config::Config)
        ty = config.type
        ln = config.ln

        function ap_decons(env, target::Target{false})::Cond
            # type check
            chk = if target.type <: ty
                TrueCond()
            else
                target.type_narrow!(ty)
                CheckCond(:($(target.repr) isa $ty))
            end

            scope = env.scope
            # compute pattern viewing if no guarantee of being computed
            view_cache = env.view_cache

            target_sym::Symbol = target.repr
            viewed_sym::Any = target_sym
            static_memo(view_cache, comp.guard1; ty = ty) do cached_sym, cache_flag
                if cache_flag === CACHE_CACHED
                    chk = AndCond(chk, CheckCond(:($cached_sym.value)))
                    return
                end

                if cache_flag === CACHE_NO_CACHE
                    guard_expr = comp.guard1(target_sym)
                    chk = AndCond(chk, CheckCond(guard_expr))
                elseif cache_flag === CACHE_MAY_CACHE
                    guard_expr = comp.guard1(target_sym)
                    do_cache = Expr(
                        :if,
                        :($cached_sym === nothing),
                        :($cached_sym = Some($guard_expr)),
                    )
                    chk = AndCond(
                        chk,
                        AndCond(TrueCond(do_cache), CheckCond(:($cached_sym.value))),
                    )
                else
                    error("impossible: invalid flag")
                end
                nothing
            end

            static_memo(view_cache, comp.view; ty = ty) do cached_sym, cache_flag
                if cache_flag === CACHE_NO_CACHE
                    viewed_sym = gensym()
                    viewed_expr = comp.view(target_sym)
                    chk = AndCond(TrueCond(:($viewed_sym = $viewed_expr)))
                elseif cache_flag === CACHE_CACHED
                    viewed_sym = :($cached_sym.value)
                elseif cache_flag === CACHE_MAY_CACHE
                    viewed_expr = comp.view(target_sym)
                    do_cache = Expr(
                        :if,
                        :($cached_sym === nothing),
                        :($cached_sym = Some($viewed_expr)),
                    )
                    chk = AndCond(chk, TrueCond(do_cache))
                    viewed_sym = :($cached_sym.value)
                else
                    error("impossible: invalid flag")
                end
                nothing
            end

            static_memo(
                view_cache,
                comp.guard2;
                ty = ty,
                depend = comp.view,
            ) do cached_sym, cache_flag
                if cache_flag === CACHE_CACHED
                    chk = AndCond(chk, CheckCond(:($cached_sym.value)))
                    return
                end

                if cache_flag === CACHE_NO_CACHE
                    guard_expr = comp.guard2(viewed_sym)
                    chk = AndCond(chk, CheckCond(guard_expr))
                elseif cache_flag === CACHE_MAY_CACHE
                    guard_expr = comp.guard2(viewed_sym)
                    do_cache = Expr(
                        :if,
                        :($cached_sym === nothing),
                        :($cached_sym = Some($guard_expr)),
                    )
                    chk = AndCond(
                        chk,
                        AndCond(TrueCond(do_cache), CheckCond(:($cached_sym.value))),
                    )
                else
                    error("impossible: invalid flag")
                end
                nothing
            end

            for i in eachindex(ps)
                p = ps[i]::Function
                field_target = Target{true}(
                    extract(viewed_sym, i, scope, ln),
                    Ref{TypeObject}(Any),
                )
                view_cache′ = ViewCache()
                env′ = env(; view_cache = view_cache′)
                elt_chk = p(env′, field_target)
                chk = AndCond(chk, AndCond(TrueCond(init_cache(view_cache′)), elt_chk))
            end
            chk
        end |> cache
    end

    function guard(pred::Function, config::Config)
        function ap_guard(env, target::Target{false})::Cond
            expr = pred(target.repr, env.scope, config.ln)
            expr === true ? TrueCond() : CheckCond(expr)
        end |> cache
    end

    function effect(perf::Function, config::Config)
        function ap_effect(env, target::Target{false})::Cond
            TrueCond(perf(target.repr, env.scope, config.ln))
        end |> cache
    end

    (
        and = and,
        or = or,
        literal = literal,
        wildcard = wildcard,
        decons = decons,
        guard = guard,
        effect = effect,
    )
end

const redy_impl = myimpl()

function compile_spec!(
    env::CompileEnv,
    suite::Vector{Any},
    x::Shaped,
    target::Target{IsComplex},
) where {IsComplex}
    if IsComplex && !(x.case isa Leaf)
        sym = gensym()
        push!(suite, :($sym = $(target.repr)))
        target = target.with_repr(sym, Val(false))
    end
    mkcond = re_tagless(x.pattern, x.ln)(redy_impl)
    ln = x.ln
    push!(suite, ln)
    cond = mkcond(env, target)
    conditional_expr = to_expression(cond)
    true_clause = Expr(:block)
    compile_spec!(env, true_clause.args, x.case, target)
    push!(
        suite,
        conditional_expr === true ? true_clause :
        Expr(:if, conditional_expr, true_clause),
    )
end

function compile_spec!(env::CompileEnv, suite::Vector{Any}, x::Leaf, target::Target)
    body = env.terminal[x.cont]
    ret = env.ret
    if env.hygienic
        bound = Expr(:block)
        let_expr = Expr(:let, bound, body)
        for_chaindict(env.scope) do k, v
            push!(bound.args, Expr(:(=), k, v))
        end
        push!(suite, Expr(:(=), ret, let_expr))
    else
        for_chaindict(env.scope) do k, v
            push!(suite, :($k = $v))
        end
        push!(suite, :($ret = $body))
    end
    push!(suite, CFGJump(env.final))
end

function compile_spec!(
    env::CompileEnv,
    suite::Vector{Any},
    x::SwitchCase,
    target::Target{IsComplex},
) where {IsComplex}
    if IsComplex
        sym = gensym()
        push!(suite, :($sym = $(target.repr)))
        target = target.with_repr(sym, Val(false))
    else
        sym = target.repr
    end

    for (ty, case) in x.cases
        true_clause = Expr(:block)
        # create new `view_cache` as only one case will be executed
        view_cache′ = child(env.view_cache)
        env′ = env(; scope = child(env.scope), view_cache = view_cache′)
        compile_spec!(env′, true_clause.args, case, target.with_type(ty))
        update_parent!(view_cache′)
        push!(suite, Expr(:if, :($sym isa $ty), true_clause))
    end
end

function compile_spec!(
    env::CompileEnv,
    suite::Vector{Any},
    x::EnumCase,
    target::Target{IsComplex},
) where {IsComplex}
    if IsComplex
        sym = gensym()
        push!(suite, :($sym = $(target.repr)))
        target = target.with_repr(sym, Val(false))
    end
    for case in x.cases
        # use old view_cache:
        # cases are tried in order,
        # hence `view_cache` can inherit from the previous case
        env′ = env(; scope = child(env.scope))
        compile_spec!(env′, suite, case, target.clone)
    end
end

function compile_spec(
    env::CompileEnv,
    target::Any,
    case::AbstractCase,
    ln::LineNumberNode,
)
    target = Target{true}(target, Ref{TypeObject}(Any))
    ret = Expr(:block)
    suite = ret.args
    if env.hygienic
        push!(suite, :($(env.ret) = nothing))
    end
    view_cache = env.view_cache

    compile_spec!(env, suite, case, target)
    pushfirst!(suite, init_cache(view_cache))

    msg = "matching non-exhaustive, at $ln"
    push!(suite, Expr(:call, error, msg))
    push!(suite, CFGLabel(env.final))
    push!(suite, env.ret)
    if env.hygienic
        # during match process can have functions that write variables.
        # hence,
        # 1. if the match expression is at global scope, write operations not allowed
        # 2. those write operations might affect outside scope.
        # to address this:
        ret = Expr(:let, Expr(:block), ret)
        # this applies to @match, but not to @switch
        # when using @switch, take care about unexpected write operations and
        # global variable issues.
    end
    CFGSpec(ret)
end

"""compile a series of `Term => Symbol`(match clauses/branches) to a Julia expression
"""
function backend(
    expr_to_match::Any,
    clauses::Vector{Pair{Function, Tuple{LineNumberNode, Int}}},
    terminal::Terminal,
    ln::LineNumberNode;
    hygienic::Bool = true,
)
    spec = spec_gen(clauses)
    env = CompileEnv(terminal, hygienic, gensym(:return), gensym(:final))
    compile_spec(env, expr_to_match, spec, ln)
end
end


text:
using ImageCore, Colors, ColorVectorSpace
using Test, Statistics

# Different access patterns (getindex)
function mysum_elt_boundscheck(A)
    s = zero(eltype(A))
    for a in A
        s += a
    end
    s
end
function mysum_index_boundscheck(A)
    s = zero(eltype(A))
    for I in eachindex(A)
        s += A[I]
    end
    s
end
function mysum_elt_inbounds(A)
    s = zero(eltype(A))
    @inbounds for a in A
        s += a
    end
    s
end
function mysum_index_inbounds_simd(A)
    s = zero(eltype(A))
    @inbounds @simd for I in eachindex(A)
        s += A[I]
    end
    s
end
# setindex!
function myfill1!(A, val)
    f = convert(eltype(A), val)
    for I in eachindex(A)
        A[I] = f
    end
    A
end
function myfill2!(A, val)
    f = convert(eltype(A), val)
    @inbounds @simd for I in eachindex(A)
        A[I] = f
    end
    A
end

# Rather than using BenchmarkTools (and thus run one test repeatedly,
# accumulating timings), we run the same test interleaving the two
# array types. This is designed to reduce the risk of spurious
# failure, particularly on shared machines like Travis where they may
# get "distracted" by other tasks
function test_getindex(f, ar, cv, n)
    t_ar = Array{Float64}(undef, n)
    t_cv = Array{Float64}(undef, n)
    # Store the results to prevent the compiler from eliding the call
    f_ar = Ref(f(ar))
    f_cv = Ref(f(cv))
    @test f_ar[] ≈ f_cv[]  # but this also gives us a chance to test correctness
    for i = 1:n
        t_ar[i] = (tstart = time(); f_ar[] = f(ar); time()-tstart)
        t_cv[i] = (tstart = time(); f_cv[] = f(cv); time()-tstart)
    end
    median(t_ar), median(t_cv)
end
function test_setindex(f, ar, cv, n)
    t_ar = Array{Float64}(undef, n)
    t_cv = Array{Float64}(undef, n)
    for i = 1:n
        t_ar[i] = @elapsed f(ar, zero(eltype(ar)))
        t_cv[i] = @elapsed f(cv, zero(eltype(cv)))
    end
    median(t_ar), median(t_cv)
end

cc_getindex_funcs = (mysum_elt_boundscheck,
                     mysum_index_boundscheck,
                     mysum_elt_inbounds,
                     mysum_index_inbounds_simd)
cc_setindex_funcs = (myfill1!,
                     myfill2!)

# Performance tolerances
isfast = VERSION >= v"1.6.0-DEV.1083"
chanvtol = Dict{Any,Int}(mysum_index_inbounds_simd => isfast ? 3 : 20,
                         mysum_elt_boundscheck => isfast ? 3 : 20,
                         myfill1! => 20,
                         myfill2! => isfast ? 3 : 20)
chanvdefault = isfast ? 3 : 10
colvtol = Dict{Any,Int}(mysum_elt_boundscheck=>isfast ? 3 : 5,
                        mysum_index_boundscheck=>isfast ? 3 : 5)
colvdefault = 3

ssz = (1000,300)

@info "Benchmark tests are warnings for now"
# @testset "benchmarks" begin
for T in (Float32, Float64)
    c = rand(RGB{T}, ssz...)
    a = copy(reinterpretc(T, c))
    vchan = channelview(c)
    vcol = colorview(RGB, a)

    # view versions
    rview = 2:ssz[1]-1
    csub = view(c, rview, :)
    asub = view(a, :, rview, :)
    vchansub = channelview(csub)
    vcolsub = colorview(RGB, asub)

    for (suite, testf) in ((cc_getindex_funcs, test_getindex),
                           (cc_setindex_funcs, test_setindex))
        for f in suite
            # channelview
            t_ar, t_cv = testf(f, a, vchan, 30)
            tol = haskey(chanvtol, f) ? chanvtol[f] : chanvdefault
            if t_cv >= tol*t_ar
                @warn "channelview1: failed on $f with eltype $T, time ratio $(t_cv/t_ar), tol $tol"
            end

            t_ar, t_cv = testf(f, asub, vchansub, 30)
            tol = haskey(chanvtol, f) ? chanvtol[f] : chanvdefault
            if t_cv >= tol*t_ar
                @warn "channelview2: failed on $f with eltype $T, time ratio $(t_cv/t_ar), tol $tol"
            end

            # colorview
            t_ar, t_cv = testf(f, c, vcol, 30)
            tol = haskey(colvtol, f) ? colvtol[f] : colvdefault
            if t_cv >= tol*t_ar
                @warn "colorview1: failed on $f with eltype $T, time ratio $(t_cv/t_ar), tol $tol"
            end

            t_ar, t_cv = testf(f, csub, vcolsub, 30)
            tol = haskey(colvtol, f) ? colvtol[f] : colvdefault
            if t_cv >= tol*t_ar
                @warn "colorview2: failed on $f with eltype $T, time ratio $(t_cv/t_ar), tol $tol"
            end
        end
    end
end
# end


text:
"""
    PointVortex(x,y,q)

Construct a point vortex of charge `q` at `x`,`y`.

See also: [`vortex_array`](@ref), [`ScalarVortex`](@ref), [`vortex!`](@ref), [`rand_pointvortex`](@ref)
"""
function PointVortex(v::Array{Float64,2})
    ~isempty(v) ? (return PointVortex.(v[:,1],v[:,2],v[:,3])) : (return Array{PointVortex}(undef,0))
end
function PointVortex(vort::Array{ScalarVortex{T},1}) where T
    pv = PointVortex[]
    for j in eachindex(vort)
        push!(pv,vort[j].vort)
    end
    return pv
end

"""
    varray = vortex_array(v::PointVortex)

Return `varray<:Array{Float64,2}` of vortex data.

See also: [`PointVortex`](@ref)
"""
vortex_array(v::PointVortex) = [v.xv v.yv v.qv]

function vortex_array(v::Array{PointVortex,1})
    ~isempty(v) > 0 ? (return reduce(vcat,vortex_array.(v))) : (return Array{Float64}(undef, 0, 0))
end

function uniform(a,b)
    @assert a<b
    return a + (b-a)*rand()
end

function uniform(a,b,n)
    return uniform.(a*ones(n),b)
end

uniform() = uniform(-.5,.5)
uniform(n) = uniform(-.5,.5,n)

rand_charge() = rand([-1 1],)
rand_charge(n) = rand([-1 1],n)

"""
    rand_pointvortex()
    rand_pointvortex(n::Int)
    rand_pointvortex(ψ::Field)
    rand_pointvortex(n::Int,ψ::Field)

Sample `n::Int` random point vortices on spatial domain specified by `ψ::Field`.

See also: [`PointVortex`](@ref), [`randVortex`](@ref)
"""
rand_pointvortex() = PointVortex(uniform(),uniform(),rand_charge())
rand_pointvortex(n) = PointVortex.(uniform(n), uniform(n), rand_charge(n))

function rand_pointvortex(psi::Field)
    @unpack x,y = psi; a,b = first(x),last(x); c,d = first(y), last(y)
    return PointVortex(uniform(a,b,),uniform(c,d,),rand_charge())
end

function rand_pointvortex(n,psi::Field)
    @unpack x,y = psi;
    dx = x[2]-x[1]; dy = y[2] - y[1]
    a,b = first(x)+2dx,last(x)-2dx; c,d = first(y)+2dy, last(y)-2dy
    return PointVortex.(uniform(a,b,n),uniform(c,d,n),rand_charge(n))
end

charge(v::PointVortex) = v.qv
xpos(v::PointVortex) = v.xv
ypos(v::PointVortex) = v.yv
vpos(v::PointVortex) = [v.xv v.yv]

xpos(v::Array{PointVortex,1}) = xpos.(v)
ypos(v::Array{PointVortex,1}) = ypos.(v)
charge(v::Array{PointVortex,1}) = charge.(v)
pos(v::PointVortex)= [xpos(v) ypos(v)]
pos(v::Array{PointVortex,1}) = [xpos.(v) ypos.(v)]

pos(v::Dipole) = [xpos(v.vp) ypos(v.vp)]


text:
"""
```
    NonmaximaSuppression <: AbstractEdgeThinningAlgorithm
    NonmaximaSuppression(; threshold::Union{Number, Percentile} = Percentile(20))

    f = NonmaximaSuppression()
    f(out::AbstractArray, mag::AbstractArray, g₁::AbstractArray, g₂::AbstractArray, f::NonmaximaSuppression)
```

Isolates local maxima of gradient magnitude `mag` along the local gradient
direction and stores the result in `out`.  The arguments `g₁` and `g₂` represent
the  gradient in the first spatial dimension (y), and the second spatial
dimension (x), respectively.

# Details

TODO

# Example

```julia

using TestImages, FileIO, ImageView, ImageEdgeDetection, ImageFiltering

img =  testimage("mandril_gray")

# Gradient in the first and second spatial dimension
g₁, g₂ = imgradients(img, KernelFactors.scharr)

# Gradient magnitude
mag = hypot.(g₁, g₂)

nms = zeros(eltype(mag), axes(mag))
# Instantiate the NonmaximaSuppression functor.
f = NonmaximaSuppression()

# Suppress the non-maximal gradient magnitudes and store the result in `nms`.
f(nms, mag, g₁, g₂)

imshow(img)
imshow(mag)
imshow(nms)
```

# References
J. Canny, "A Computational Approach to Edge Detection," in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. PAMI-8, no. 6, pp. 679-698, Nov. 1986, doi: 10.1109/TPAMI.1986.4767851.
"""
@with_kw struct NonmaximaSuppression{ T <: Union{Real,AbstractGray, Percentile}} <: AbstractEdgeThinningAlgorithm
    threshold::T = Percentile(20)
end

function (f::NonmaximaSuppression)(out::AbstractArray, mag::AbstractArray, g₁::AbstractArray, g₂::AbstractArray)
    @unpack threshold = f

    low_threshold =  typeof(threshold) <: Percentile ? StatsBase.percentile(vec(mag), threshold.p) : threshold

    # Isolate local maxima of gradient magnitude by “non-maximum suppression”
    # along the local gradient direction.
    suppress_non_maxima!(out, mag, g₁, g₂, low_threshold)

    return out
end

"""
```
    suppress_non_maxima!(out::AbstractArray, mag::AbstractArray, gx::AbstractArray, gy::AbstractArray, threshold::Number)
```

Isolates local maxima of gradient magnitude by “non-maximum suppression” along the local gradient direction.

"""
function suppress_non_maxima!(out::AbstractArray, mag::AbstractArray, g₁::AbstractArray, g₂::AbstractArray, threshold::Number)
    itp = interpolate(mag, BSpline(Linear()))
    etp = extrapolate(itp, Flat())
    m₁ = zero(eltype(mag))
    m₂ = zero(eltype(mag))
    @inbounds for i in CartesianIndices(mag)
        r, c = i.I
        d₁ = gray(g₁[i])
        d₂ = gray(g₂[i])
        mc = mag[i]
        if mc < threshold || mc == 0 || isnan(mc)
            out[r,c] = zero(eltype(mag))
        else
            # Ensure the vector 𝐝 = [d₁, d₂] has unit norm.
            d₁ = d₁ / mc
            d₂ = d₂ / mc
            m₁ = etp(r + d₁, c + d₂)
            m₂ = etp(r - d₁, c - d₂)
            out[r,c] = (m₁ <= mc) && (mc >= m₂) ? mag[r,c] : zero(eltype(mag))
        end
    end
    return nothing
end