text:
"""
Represents a state in an set of abitrary (yet low-rank) subspaces of a set of FockConfigs.
e.g.

    v[FockConfig][TuckerConfig] => Tucker Decomposed Tensor

# Data
- `clusters::Vector{Cluster}`
- `data::OrderedDict{FockConfig,OrderedDict{TuckerConfig,Tucker}}`
- `p_spaces::Vector{ClusterSubspace}`
- `q_spaces::Vector{ClusterSubspace}`
"""
struct CompressedTuckerState{T,N}
    clusters::Vector{Cluster}
    data::OrderedDict{FockConfig{N},OrderedDict{TuckerConfig{N},Tucker{T,N}}}
    p_spaces::Vector{ClusterSubspace}
    q_spaces::Vector{ClusterSubspace}
end
Base.haskey(ts::CompressedTuckerState, i) = return haskey(ts.data,i)
Base.getindex(ts::CompressedTuckerState, i) = return ts.data[i]
Base.setindex!(ts::CompressedTuckerState, i, j) = return ts.data[j] = i
Base.iterate(ts::CompressedTuckerState, state=1) = iterate(ts.data, state)
normalize!(ts::CompressedTuckerState) = scale!(ts, 1/sqrt(orth_dot(ts,ts)))

"""
    CompressedTuckerState(clusters::Vector{Cluster},
        fconfig::FockConfig{N},
        cluster_bases::Vector{ClusterBasis}) where {N}

Constructor using only a single FockConfig. This allows us to turn the CMF state into a CompressedTuckerState.
As such, it chooses the ground state of each cluster in the Fock sector specified by `FockConfig` to be the
P space, and then the Q space is defined as the orthogonal complement of this state within the available basis,
specified by `cluster_bases`.
# Arguments
- `clusters`: vector of clusters types
- `fconfig`: starting FockConfig
- `cluster_basis`: list of ClusterBasis types - needed to know the dimensions of the q-spaces
# Returns
- `CompressedTuckerState`
"""
function CompressedTuckerState(clusters::Vector{Cluster},
        fconfig::FockConfig{N},
        cluster_bases::Vector{ClusterBasis}) where {N}
    #={{{=#
    T = Float64

    #
    # start by building the P and Q spaces needed
    p_spaces = Vector{ClusterSubspace}()
    q_spaces = Vector{ClusterSubspace}()
    # define p spaces
    for ci in clusters
        tss = ClusterSubspace(ci)
        tss[fconfig[ci.idx]] = 1:1
        push!(p_spaces, tss)
    end

    # define q spaces
    for tssp in p_spaces
        tss = get_ortho_compliment(tssp, cluster_bases[tssp.cluster.idx])
        push!(q_spaces, tss)
    end

    data = OrderedDict{FockConfig{N},OrderedDict{TuckerConfig{N},Tucker{T,N}} }()
    state = CompressedTuckerState(clusters, data, p_spaces, q_spaces)
    add_fockconfig!(state, fconfig)

    factors = []
    for ci in clusters
        dim = length(p_spaces[ci.idx][fconfig[ci.idx]])
        push!(factors, 1.0Matrix(I, dim, 1))
    end
    factors = tuple(factors...)

    tconfig = TuckerConfig([p_spaces[ci.idx].data[fconfig[ci.idx]] for ci in clusters])
    tdata = Tucker(reshape([1.0], tuple(ones(Int64, N)...)), factors)

    state[fconfig][tconfig] = tdata
    return state
#=}}}=#
end

"""
    CompressedTuckerState(clusters::Vector{Cluster},
        p_spaces::Vector{FermiCG.ClusterSubspace},
        q_spaces::Vector{FermiCG.ClusterSubspace})

Constructor - specify input p and q spaces
# Arguments
- `clusters`: vector of clusters types
- `p_spaces`: list of p space ranges for each cluster
- `q_spaces`: list of q space ranges for each cluster
# Returns
- `CompressedTuckerState`
"""
function CompressedTuckerState(clusters::Vector{Cluster},
        p_spaces::Vector{FermiCG.ClusterSubspace},
        q_spaces::Vector{FermiCG.ClusterSubspace})
    #={{{=#

    N = length(clusters)
    T = Float64

    #factors = []
    #for ci in clusters
    #    dim = length(p_spaces[ci.idx][init_fspace[ci.idx]])
    #    push!(factors, 1.0Matrix(I, dim, 1))
    #end
    #factors = tuple(factors...)

    #tconfig = TuckerConfig([p_spaces[ci.idx].data[init_fspace[ci.idx]] for ci in clusters])
    #fconfig = FockConfig(init_fspace)
    #tdata = Tucker(reshape([1.0], tuple(ones(Int64, N)...)), factors)
    data = OrderedDict{FockConfig{N},OrderedDict{TuckerConfig{N},Tucker{T,N}} }()

    #data[fconfig][tconfig] = tdata
    return CompressedTuckerState(clusters, data, p_spaces, q_spaces)
#=}}}=#
end

"""
    CompressedTuckerState(ts::TuckerState; thresh=-1, max_number=nothing, verbose=0)

Create a `CompressedTuckerState` from a `TuckerState`
# Arguments
- `ts::TuckerState`
- `thresh=-1`: discard singular values smaller than `thresh`
- `max_number=nothing`: if != `nothing`, only keep up to `max_number` singular vectors per SVD
- `verbose=0`: print level
# Returns
- `CompressedTuckerState`
"""
function CompressedTuckerState(ts::TuckerState{T,N}; thresh=-1, max_number=nothing, verbose=0) where {T,N}
#={{{=#
    # make all AbstractState subtypes parametric
    nroots = nothing
    for (fock,configs) in ts
        for (config,coeffs) in configs
            if nroots == nothing
                nroots = last(size(coeffs))
            else
                nroots == last(size(coeffs)) || error(" mismatch in number of roots")
            end
        end
    end

    nroots == 1 || error(" Conversion to CompressedTuckerState can only have 1 root")

    data = OrderedDict{FockConfig{N},OrderedDict{TuckerConfig{N},Tucker{T,N} }}()
    for (fock, tconfigs) in ts.data
        for (tconfig, coeffs) in tconfigs

            #
            # Since TuckerState has extra dimension for state index, remove that
            tuck = Tucker(reshape(coeffs,size(coeffs)[1:end-1]), thresh=thresh, max_number=max_number, verbose=verbose)
            if length(tuck) > 0
                if haskey(data, fock)
                    data[fock][tconfig] = tuck
                else
                    data[fock] = OrderedDict(tconfig => tuck)
                end
            end
        end
    end
    return CompressedTuckerState(ts.clusters, data, ts.p_spaces, ts.q_spaces)
end
#=}}}=#

"""
    compress(ts::CompressedTuckerState{T,N}; thresh=-1, max_number=nothing, verbose=0) where {T,N}

Compress state via HOSVD
# Arguments
- `ts::CompressedTuckerState`
- `thresh = -1`: threshold for compression
- `max_number`: only keep certain number of vectors per TuckerConfig
- `verbose=0`: print level
# Returns
- `CompressedTuckerState`
"""
function compress(ts::CompressedTuckerState{T,N}; thresh=-1, max_number=nothing, verbose=0) where {T,N}
    d = OrderedDict{FockConfig{N}, OrderedDict{TuckerConfig{N}, Tucker{T,N}}}()
    for (fock, tconfigs) in ts.data
        for (tconfig, coeffs) in tconfigs
            tmp = compress(ts.data[fock][tconfig], thresh=thresh, max_number=max_number)
            if length(tmp) == 0
                continue
            end
            if haskey(d, fock)
                d[fock][tconfig] = tmp
            else
                d[fock] = OrderedDict(tconfig => tmp)
            end
        end
    end
    return CompressedTuckerState(ts.clusters, d, ts.p_spaces, ts.q_spaces)
end

"""
    orth_add!(ts1::CompressedTuckerState, ts2::CompressedTuckerState)

Add coeffs in `ts2` to `ts1`

Note: this assumes `t1` and `t2` have the same compression vectors
"""
function orth_add!(ts1::CompressedTuckerState, ts2::CompressedTuckerState)
#={{{=#
    for (fock,configs) in ts2
        if haskey(ts1, fock)
            for (config,coeffs) in configs
                if haskey(ts1[fock], config)
                    ts1[fock][config].core .+= ts2[fock][config].core
                else
                    ts1[fock][config] = ts2[fock][config]
                end
            end
        else
            ts1[fock] = ts2[fock]
        end
    end
#=}}}=#
end

"""
    nonorth_add!(ts1::CompressedTuckerState, ts2::CompressedTuckerState)

Add coeffs in `ts2` to `ts1`

Note: this does not assume `t1` and `t2` have the same compression vectors
"""
function nonorth_add!(ts1::CompressedTuckerState, ts2::CompressedTuckerState)
#={{{=#
    for (fock,configs) in ts2
        if haskey(ts1, fock)
            for (config,coeffs) in configs
                if haskey(ts1[fock], config)
                    ts1[fock][config] = ts1[fock][config] + ts2[fock][config] # note this is non-trivial work here
                else
                    ts1[fock][config] = ts2[fock][config]
                end
            end
        else
            ts1[fock] = ts2[fock]
        end
    end
#=}}}=#
end

"""
    add_fockconfig!(s::CompressedTuckerState, fock::FockConfig)
"""
function add_fockconfig!(s::CompressedTuckerState{T,N}, fock::FockConfig) where {T,N}
    s.data[fock] = OrderedDict{TuckerConfig, Tucker{T,N}}()
end

"""
    Base.length(s::CompressedTuckerState)
"""
function Base.length(s::CompressedTuckerState)
    l = 0
    for (fock,tconfigs) in s.data
        for (tconfig, tuck) in tconfigs
            l += length(tuck)
        end
    end
    return l
end
"""
    prune_empty_fock_spaces!(s::AbstractState)

remove fock_spaces that don't have any configurations
"""
function prune_empty_fock_spaces!(s::AbstractState)
    focklist = keys(s.data)
    for fock in focklist
        if length(s.data[fock]) == 0
            delete!(s.data, fock)
        end
    end
    focklist = keys(s.data)
    for (fock,tconfigs) in s.data
        for (tconfig,coeff) in tconfigs
        end
    end
end
"""
    prune_empty_TuckerConfigs!(s::T) where T<:Union{TuckerState, CompressedTuckerState}

remove fock_spaces that don't have any configurations
"""
function prune_empty_TuckerConfigs!(s::T) where T<:Union{TuckerState, CompressedTuckerState}
    focklist = keys(s.data)
    for fock in focklist
        tconflist = keys(s.data[fock])
        for tconf in tconflist
            if length(s.data[fock][tconf]) == 0
                delete!(s.data[fock], tconf)
            end
        end
    end
    for (fock,tconfigs) in s.data
        for (tconfig,coeff) in tconfigs
        end
    end
    prune_empty_fock_spaces!(s)
end

"""
    get_vector(s::CompressedTuckerState)

Return a vector of the variables. Note that this is the core tensors being returned
"""
function get_vector(cts::CompressedTuckerState)

    v = zeros(length(cts), 1)
    idx = 1
    for (fock, tconfigs) in cts
        for (tconfig, tuck) in tconfigs
            dims = size(tuck.core)

            dim1 = prod(dims)
            v[idx:idx+dim1-1,:] = copy(reshape(tuck.core,dim1))
            idx += dim1
        end
    end
    return v
end
"""
    set_vector!(s::CompressedTuckerState)
"""
function set_vector!(ts::CompressedTuckerState, v)

    #length(size(v)) == 1 || error(" Only takes vectors", size(v))
    nbasis = size(v)[1]

    idx = 1
    for (fock, tconfigs) in ts
        for (tconfig, tuck) in tconfigs
            dims = size(tuck)

            dim1 = prod(dims)
            ts[fock][tconfig].core .= reshape(v[idx:idx+dim1-1], size(tuck.core))
            idx += dim1
        end
    end
    nbasis == idx-1 || error("huh?", nbasis, " ", idx)
    return
end
"""
    zero!(s::CompressedTuckerState)
"""
function zero!(s::CompressedTuckerState)
    for (fock, tconfigs) in s
        for (tconfig, tcoeffs) in tconfigs
            fill!(s[fock][tconfig].core, 0.0)
        end
    end
end

"""
    Base.display(s::CompressedTuckerState; thresh=1e-3)

Pretty print
"""
function Base.display(s::CompressedTuckerState; thresh=1e-3)
#={{{=#
    println()
    @printf(" --------------------------------------------------\n")
    @printf(" ---------- # Fockspaces -------------------: %5i  \n",length(keys(s.data)))
    @printf(" ---------- # Configs    -------------------: %5i  \n",length(s))
    @printf(" --------------------------------------------------\n")
    @printf(" Printing contributions greater than: %f", thresh)
    @printf("\n")
    @printf(" %-20s%-10s%-10s%-20s\n", "Weight", "# configs", "(full)", "(α,β)...")
    @printf(" %-20s%-10s%-10s%-20s\n", "-------","---------", "---------", "----------")
    for (fock,configs) in s.data
        prob = 0
        len = 0

        lenfull = 0
        for (config, tuck) in configs
            prob += sum(tuck.core .* tuck.core)
            len += length(tuck.core)
            lenfull += prod(dims_large(tuck))
        end
        if prob > thresh
        #if lenfull > 0
            #@printf(" %-20.3f%-10i%-10i", prob,len, lenfull)
            @printf(" %-20.3f%-10s%-10s", prob,"","")
            for sector in fock
                @printf("(%2i,%-2i)", sector[1],sector[2])
            end
            println()

            #@printf("     %-16s%-20s%-20s\n", "Weight", "", "Subspaces")
            #@printf("     %-16s%-20s%-20s\n", "-------", "", "----------")
            for (config, tuck) in configs
                probi = sum(tuck.core .* tuck.core)
                @printf("     %-16.3f%-10i%-10i", probi,length(tuck.core),prod(dims_large(tuck)))
                for range in config
                    @printf("%7s", range)
                end
                println()
            end
            #println()
            @printf(" %-20s%-20s%-20s\n", "---------", "", "----------")
        end
    end
    print(" --------------------------------------------------\n")
    println()
#=}}}=#
end
"""
    print_fock_occupations(s::CompressedTuckerState; thresh=1e-3)

Pretty print
"""
function print_fock_occupations(s::CompressedTuckerState; thresh=1e-3)
#={{{=#

    println()
    @printf(" --------------------------------------------------\n")
    @printf(" ---------- # Fockspaces -------------------: %5i  \n",length(keys(s.data)))
    @printf(" ---------- # Configs    -------------------: %5i  \n",length(s))
    @printf(" --------------------------------------------------\n")
    @printf(" Printing contributions greater than: %f", thresh)
    @printf("\n")
    @printf(" %-20s%-10s%-10s%-20s\n", "Weight", "# configs", "(full)", "(α,β)...")
    @printf(" %-20s%-10s%-10s%-20s\n", "-------","---------", "---------", "----------")
    for (fock,configs) in s.data
        prob = 0
        len = 0
        lenfull = 0
        for (config, tuck) in configs
            prob += sum(tuck.core .* tuck.core)
            len += length(tuck.core)
            lenfull += prod(dims_large(tuck))
        end
        if prob > thresh
            @printf(" %-20.3f%-10i%-10i", prob,len,lenfull)
            for sector in fock
                @printf("(%2i,%-2i)", sector[1],sector[2])
            end
            println()
        end
    end
    print(" --------------------------------------------------\n")
    println()
#=}}}=#
end

"""
    dot(ts1::FermiCG.CompressedTuckerState, ts2::FermiCG.CompressedTuckerState)

Dot product between `ts2` and `ts1`

Warning: this assumes both `ts1` and `ts2` have the same tucker factors for each `TuckerConfig`
"""
function orth_dot(ts1::CompressedTuckerState, ts2::CompressedTuckerState)
#={{{=#
    overlap = 0.0
    for (fock,configs) in ts2
        haskey(ts1, fock) || continue
        for (config,coeffs) in configs
            haskey(ts1[fock], config) || continue
            overlap += sum(ts1[fock][config].core .* ts2[fock][config].core)
        end
    end
    return overlap
#=}}}=#
end

"""
    nonorth_dot(ts1::FermiCG.CompressedTuckerState, ts2::FermiCG.CompressedTuckerState; verbose=0)

Dot product between 1ts2` and `ts1` where each have their own Tucker factors
"""
function nonorth_dot(ts1::CompressedTuckerState, ts2::CompressedTuckerState; verbose=0)
#={{{=#
    overlap = 0.0
    for (fock,configs) in ts2
        haskey(ts1, fock) || continue
        verbose == 0 || display(fock)
        for (config,coeffs) in configs
            haskey(ts1[fock], config) || continue
            verbose == 0 || display(config)
            overlap += dot(ts1[fock][config] , ts2[fock][config])
            verbose == 0 || display(dot(ts1[fock][config] , ts2[fock][config]))
        end
    end
    return overlap
#=}}}=#
end

"""
    scale!(ts::FermiCG.CompressedTuckerState, a::T<:Number)

Scale `ts` by a constant
"""
function scale!(ts::CompressedTuckerState, a::T) where T<:Number
    #={{{=#
    for (fock,configs) in ts
        for (config,tuck) in configs
            ts[fock][config].core .*= a
        end
    end
    #=}}}=#
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
@testset "Image Partitioning" begin
    img = testimage("mandrill")
    # Window size divides image size without remainder.
    subdivision = partition_image(AllowContraction(), img, 32)
    size(subdivision) == (16,16)
    for i in subdivision
        @test length.(i) == (32,32)
    end

    # Window size divides image size with remainder.
    subdivision = partition_image(AllowContraction(), img, 31)
    for r = 1:18
        for c = 1:18
            t = subdivision[r,c]
            if c == 1 || c == 18
                @test length(t[2]) == 8
            end
            if  1 < c < 18
                @test length(t[2]) == 31
            end
            if r == 1 || r == 18
                @test length(t[1]) == 8
            end
            if  1 < r < 18
                @test length(t[1]) == 31
            end
        end
    end
end

ext: jl


text:
using Metatheory
using Metatheory.EGraphs
using Metatheory.Library
using Rewriters

using Test

Cat = @theory begin
    id(A).(A→A) ⋅ f.(A→B) == f.(A→B)
    f.(A→B)⋅id(B).(B→B) == f.(A→B)
    (f.(A→B)⋅g.(B→C))⋅h.(C→D) == f.(A→B)⋅(g.(B→C)⋅h.(C→D))
end

# DOES NOT FIXPOINT!
tag_matcher_t = @theory begin
    sigma(A, B) => σ(A, B).(A⊗B→B⊗A)
    1(A) => id(A).(A→A)
    :f |> :(f.(A → B))
    :g |> :(g.(B → C))
    :h |> :(h.(B → C))
    :j |> :(j.(B → B))
    :k |> :(j.(C → D))
end

tag_matcher(x) = Chain(tag_matcher_t)(x)

# FIXME TAGGING skip . operator
function tag(x)
    r = Postwalk(PassThrough(If(x -> !istree(x) || operation(x) != :., tag_matcher)))(x)
    println("tagged $x to $(r)")
    r
end

macro areequal_tag(t, exprs...)
    exprs = map(tag, exprs)

    :(areequal($t, $(exprs)...))
end

@testset "Cat" begin
    @test @areequal_tag Cat 1(B) id(B).(B→B)
    @test @areequal_tag Cat (1(A) ⋅ f) f
    @test @areequal_tag Cat (1(A)⋅f) f
    @test @areequal_tag Cat (f⋅1(B)) f
    @test @areequal_tag Cat (f⋅j)⋅j f⋅(j⋅j)
    @test @areequal_tag Cat (f⋅j)⋅h f⋅(j⋅h)
end

MonCat = Cat ∪ monoid(:(⊗), :(:munit)) ∪ @theory begin
    f.(A→B) ⊗ g.(C→D) == (f⊗g).(A⊗C→B⊗D)
    id(A⊗B).(A⊗B→A⊗B) == id(A).(A→A)⊗id(B).(B→B)
end

@testset "MonCat" begin
    @test @areequal_tag MonCat 1(A)⊗1(B) 1(A⊗B)
    @test @areequal_tag MonCat (f ⊗ g) ((f⊗g).(A⊗B→B⊗C))
    @test @areequal_tag MonCat 1(A⊗B)⋅(f⊗g) (f⊗g)
    @test @areequal_tag MonCat ((f ⊗ g)⋅1(B⊗C))  (f⊗g)
    @test @areequal_tag MonCat 1(A⊗B)⋅(f⊗g) (f⊗g).(A⊗B→B⊗C)
end

SymMonCat = MonCat ∪ @theory begin
    σ(A, B).(A⊗B→B⊗A) ⋅ σ(B, A).(B⊗A→A⊗B) == id(A⊗B).(A⊗B→A⊗B)

    σ(A, B⊗C).(A⊗(B⊗C)→(B⊗C)⊗A) == σ(A⊗B, C).((A⊗B)⊗C→C⊗(A⊗B))
    (f.(A→B)⊗g.(C→D))⋅σ(B,D).(B⊗D→D⊗B) == σ(A,C).(A⊗C→C⊗A) ⋅ (g.(C→D)⊗f.(A→B))
    σ(A,C).(A⊗C→C⊗A) ⋅ (g.(C→D)⊗f.(A→B)) == (f.(A→B)⊗g.(C→D))⋅σ(B,D).(B⊗D→D⊗B)
end

@testset "SymMonCat" begin
    @test @areequal_tag SymMonCat 1(A)⊗1(B) 1(A⊗B)
    @test @areequal_tag SymMonCat (f ⊗ g) (f ⊗ g)
    @test @areequal_tag SymMonCat 1(A⊗B) ⋅(f ⊗ g) (f ⊗ g)
    @test @areequal_tag SymMonCat (f ⊗ g) ⋅ 1(B⊗C) (f ⊗ g)
    @test @areequal_tag SymMonCat 1(A⊗B)⋅(f⊗g) (f⊗g).(A⊗B→B⊗C)

    # println("==========================================")

    @test falseormissing(@areequal_tag SymMonCat sigma(A,B)⋅sigma(A,B)   1(A⊗B))
    @test @areequal_tag SymMonCat sigma(A,B) σ(A,B).(A⊗B→B⊗A)
    @test @areequal_tag SymMonCat sigma(A,B)⋅sigma(B,A) 1(A⊗B)
    @test @areequal_tag SymMonCat sigma(B,A)⋅sigma(A,B) 1(B⊗A)
    @test @areequal_tag SymMonCat (f ⊗ k)⋅sigma(B,D) sigma(A,C)⋅(k ⊗ f)
    @test @areequal_tag SymMonCat sigma(A,A)⋅(f.(A→A)⊗g.(A→A))⋅sigma(A,A) g.(A→A)⊗f.(A→A)
    @test @areequal_tag SymMonCat sigma(B,A)⋅(f.(A→A)⊗g.(B→A))⋅sigma(A,A) g.(B→A)⊗f.(A→A)
    @test @areequal_tag SymMonCat sigma(B,A)⋅(f.(A→C)⊗g.(B→D))⋅sigma(C,D) g.(B→D)⊗f.(A→C)
end


text:
export CuForwardEulerTimeStepper, CuFilteredForwardEulerTimeStepper,
       CuRK4TimeStepper, CuFilteredRK4TimeStepper

struct CuForwardEulerTimeStepper{T,dim} <: AbstractForwardEulerTimeStepper
  dt::Float64
  N::CuArray{Complex{T},dim}    # Explicit linear and nonlinear terms
  CuForwardEulerTimeStepper(dt::T, N::CuArray{Complex{T},dim}) where {T,dim} = new{T,dim}(
    dt, CuArray(zeros(size(N))))
end

CuForwardEulerTimeStepper(dt, N::Array{Complex{T},dim}) where {T,dim} = CuForwardEulerTimeStepper(dt, CuArray(N))

struct CuFilteredForwardEulerTimeStepper{T,dim} <: AbstractFilteredForwardEulerTimeStepper
  dt::Float64
  N::CuArray{Complex{T},dim}    # Explicit linear and nonlinear terms
  filter::CuArray{T,dim}        # Filter for solution
end

function CuFilteredForwardEulerTimeStepper(dt, LC, g; filterkwargs...)
  N = CuArray(zeros(cxeltype(LC), size(LC)))
  filter = makefilter(g, typeof(dt), size(LC); filterkwargs...)
  CuFilteredForwardEulerTimeStepper{ndims(N)}(dt, N, filter)
end

struct CuRK4TimeStepper{T,dim} <: AbstractRK4TimeStepper
  dt::Float64
  sol₁::CuArray{T,dim}
  RHS₁::CuArray{T,dim}
  RHS₂::CuArray{T,dim}
  RHS₃::CuArray{T,dim}
  RHS₄::CuArray{T,dim}
end

function CuRK4TimeStepper(dt, LC)
  @createcuarrays cxeltype(LC) size(LC) sol₁ RHS₁ RHS₂ RHS₃ RHS₄
  CuRK4TimeStepper{cxeltype(LC),ndims(LC)}(dt, sol₁, RHS₁, RHS₂, RHS₃, RHS₄)
end

struct CuFilteredRK4TimeStepper{T,dim} <: AbstractFilteredRK4TimeStepper
  dt::Float64
  sol₁::CuArray{T,dim}
  RHS₁::CuArray{T,dim}
  RHS₂::CuArray{T,dim}
  RHS₃::CuArray{T,dim}
  RHS₄::CuArray{T,dim}
  filter::CuArray{T,dim}    # Filter for solution
end

function CuFilteredRK4TimeStepper(dt, LC, g; filterkwargs...)
  @createcuarrays cxeltype(LC) size(LC) sol₁ RHS₁ RHS₂ RHS₃ RHS₄
  filter = makefilter(g, typeof(dt), size(LC); filterkwargs...)
  CuFilteredRK4TimeStepper{cxeltype(LC),ndims(LC)}(dt, sol₁, RHS₁, RHS₂, RHS₃, RHS₄, filter)
end

ForwardEulerTimeStepper(dt, LC::CuArray{T,dim}) where {T,dim} = CuForwardEulerTimeStepper(dt, LC)
RK4TimeStepper(dt, LC::CuArray{T,dim}) where {T,dim} = CuRK4TimeStepper(dt, LC)

FilteredForwardEulerTimeStepper(
  dt, LC::CuArray{T,dim}, g; filterkwargs...) where {T,dim} = CuFilteredForwardEulerTimeStepper(dt, LC, g;
                                                                                filterkwargs...)
FilteredRK4TimeStepper(
  dt, LC::CuArray{T,dim}, g; filterkwargs...) where {T,dim} = CuFilteredRK4TimeStepper(dt, LC, g; filterkwargs...)


text:
# include("elements/figure.jl")

# Use github flavored markdown as our starting point
mrbl = deepcopy(Markdown.github)

include("block.jl")
include("inline.jl")
# Inserted infront of the paragraph element because that's the catch-all case
insert!(mrbl.regular, length(mrbl.regular), YAMLelement)
insert!(mrbl.regular, length(mrbl.regular), Markdown.horizontalrule)
mrbl.inner['['] = push!(mrbl.inner['['], interpreted_inline)
mrbl.inner['\$'] = [inline_math]
mrbl.inner['`'] = [inline_code]
Markdown.flavors[:mrbl] = mrbl

### DEBUG INFO. DO NOT DELETE YET ###
# Markdown.@flavor mrbl [Markdown.list, Markdown.indentcode, Markdown.blockquote, Markdown.fencedcode, Markdown.hashheader,
#                 Markdown.github_table, Markdown.github_paragraph,
#
#                 Markdown.linebreak, Markdown.escapes, Markdown.en_dash, Markdown.inline_code, Markdown.asterisk_bold,
#                 Markdown.asterisk_italic, Markdown.image, Markdown.link]
#
# const mrbl = Markdown.config([Markdown.list, Markdown.indentcode, Markdown.blockquote, Markdown.fencedcode, Markdown.hashheader,
#                 Markdown.github_table, Markdown.github_paragraph,
#
#                 Markdown.linebreak, Markdown.escapes, Markdown.en_dash, Markdown.inline_code, Markdown.asterisk_bold,
#                 Markdown.asterisk_italic, Markdown.image, Markdown.link]...)
# flavors[$(Expr(:quote, name))] = $(esc(name))
#
# Markdown.github
#
# Base.Markdown.Config(
#     Function[],
#
#     [Base.Markdown.list,Base.Markdown.indentcode,Base.Markdown.blockquote,Base.Markdown.fencedcode,Base.Markdown.hashheader,Base.Markdown.github_table,Base.Markdown.github_paragraph,Base.Markdown.linebreak,Base.Markdown.escapes,Base.Markdown.en_dash,Base.Markdown.inline_code,Base.Markdown.asterisk_bold,Base.Markdown.asterisk_italic,Base.Markdown.image,Base.Markdown.link],
#
#     Dict{Char,Array{Function,1}}()
# )
#
# Base.Markdown.Config(
#     [Base.Markdown.blockquote,Base.Markdown.fencedcode,Base.Markdown.hashheader],
#
#     [Base.Markdown.list,Base.Markdown.indentcode,Base.Markdown.github_table,Base.Markdown.github_paragraph],
#
#     Dict(
#     '['=>[Base.Markdown.link],
#     '\\'=>[Base.Markdown.linebreak,Base.Markdown.escapes],
#     '*'=>[Base.Markdown.asterisk_bold,Base.Markdown.asterisk_italic],
#     '`'=>[Base.Markdown.inline_code],
#     '!'=>[Base.Markdown.image],
#     '-'=>[Base.Markdown.en_dash])
# )


text:
using Test
using ITensors
using LinearAlgebra

@testset "Ops to MPO" begin
  ∑H = Sum{Op}()
  ∑H += 1.2, "X", 1, "X", 2
  ∑H += 2, "Z", 1
  ∑H += 2, "Z", 2

  @test ∑H isa Sum{Scaled{Prod{Op},Float64}}

  s = siteinds("Qubit", 2)
  H = MPO(∑H, s)

  Id(n) = Op(I, n)
  X(n) = Op("X", n)
  Z(n) = Op("Z", n)
  T(o) = ITensor(o, s)
  Hfull = 1.2 * T(X(1)) * T(X(2)) + 2 * T(Z(1)) * T(Id(2)) + 2 * T(Id(1)) * T(Z(2))

  @test prod(H) ≈ Hfull

  @test prod(MPO(X(1), s)) ≈ T(X(1)) * T(Id(2))
  @test prod(MPO(2X(1), s)) ≈ 2T(X(1)) * T(Id(2))
  @test prod(MPO(X(1) * Z(2), s)) ≈ T(X(1)) * T(Z(2))
  @test prod(MPO(3.5X(1) * Z(2), s)) ≈ 3.5T(X(1)) * T(Z(2))
  @test prod(MPO(X(1) + Z(2), s)) ≈ T(X(1)) * T(Id(2)) + T(Id(1)) * T(Z(2))
  @test prod(MPO(X(1) + 3.3Z(2), s)) ≈ T(X(1)) * T(Id(2)) + 3.3T(Id(1)) * T(Z(2))
  @test prod(MPO((X(1) + Z(2)) / 2, s)) ≈ 0.5T(X(1)) * T(Id(2)) + 0.5T(Id(1)) * T(Z(2))
end

function heisenberg_old(N)
  os = OpSum()
  for j in 1:(N - 1)
    os += "Sz", j, "Sz", j + 1
    os += 0.5, "S+", j, "S-", j + 1
    os += 0.5, "S-", j, "S+", j + 1
  end
  return os
end

function heisenberg(N)
  os = Sum{Op}()
  for j in 1:(N - 1)
    os += "Sz", j, "Sz", j + 1
    os += 0.5, "S+", j, "S-", j + 1
    os += 0.5, "S-", j, "S+", j + 1
  end
  return os
end

@testset "OpSum comparison" begin
  N = 4
  s = siteinds("S=1/2", N)
  os_old = heisenberg_old(N)
  os_new = heisenberg(N)
  @test os_old isa OpSum
  @test os_new isa Sum{Scaled{Prod{Op},Float64}}
  Hold = MPO(os_old, s)
  Hnew = MPO(os_new, s)
  @test prod(Hold) ≈ prod(Hnew)
end

@testset "Square Hamiltonian" begin
  N = 4
  ℋ = heisenberg(N)
  ℋ² = expand(ℋ^2)
  s = siteinds("S=1/2", N)
  H = MPO(ℋ, s)
  H² = MPO(ℋ², s)
  @test norm(replaceprime(H' * H, 2 => 1) - H²) ≈ 0 atol = 1e-14
end


text:
module Tprint

import Term: unescape_brackets, escape_brackets, has_markup
import Term: highlight as highlighter

import ..Renderables: AbstractRenderable
import ..Style: apply_style
import ..Layout: hstack

export tprint, tprintln

"""
    tprint

Similar to standard lib's `print` function but with added
**styling functionality**

!!! tip "highlighting"
    Set `highlight=true` to automatically highlight the output.
"""
function tprint end

tprint(x; highlight = true) = tprint(stdout, x; highlight = highlight)

tprint(io::IO, x; highlight = true) = tprint(io, string(x); highlight = highlight)

tprint(io::IO, ::MIME"text/html", x; highlight = true) =
    tprint(io, x; highlight = highlight)

"""
---
    tprint(x::AbstractString)

Apply style to a string and print it
"""
tprint(io::IO, x::AbstractString; highlight = true) =
    print(io, (highlight ? apply_style ∘ highlighter : apply_style)(x))

"""
---
    tprint(x::AbstractRenderable)

Print an `AbstractRenderable`.

Equivalent to `println(x)`
"""
tprint(io::IO, x::AbstractRenderable; highlight = true) =
    print(io, x; highlight = highlight)

function tprint(io::IO, args...)
    for (n, arg) in enumerate(args)
        tprint(io, arg)

        if n < length(args)
            args[n + 1] isa AbstractRenderable || print(io, " ")
        end
    end
    return nothing
end

function tprint(args...; highlight = true)
    for (n, arg) in enumerate(args)
        tprint(arg; highlight = highlight)

        if n < length(args)
            args[n + 1] isa AbstractRenderable || print(" ")
        end
    end
    return nothing
end

"""
---
    tprintln

Similar to standard lib's `println` function but with added
styling functionality.
"""

tprintln(args...; highlight = true) = tprint(args..., "\n"; highlight = highlight)

end


text:
# Autogenerated wrapper script for GoogleSparseHash_jll for any
export dense_hash_map, dense_hash_set, densehashtable_h, hashtable_common_h, libc_allocator_with_realloc_h, sparse_hash_map, sparse_hash_set, sparsehashtable_h, sparsetable, traits

JLLWrappers.@generate_wrapper_header("GoogleSparseHash")
JLLWrappers.@declare_file_product(dense_hash_map)
JLLWrappers.@declare_file_product(dense_hash_set)
JLLWrappers.@declare_file_product(densehashtable_h)
JLLWrappers.@declare_file_product(hashtable_common_h)
JLLWrappers.@declare_file_product(libc_allocator_with_realloc_h)
JLLWrappers.@declare_file_product(sparse_hash_map)
JLLWrappers.@declare_file_product(sparse_hash_set)
JLLWrappers.@declare_file_product(sparsehashtable_h)
JLLWrappers.@declare_file_product(sparsetable)
JLLWrappers.@declare_file_product(traits)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_file_product(
        dense_hash_map,
        "include/sparsehash/dense_hash_map",
    )

    JLLWrappers.@init_file_product(
        dense_hash_set,
        "include/sparsehash/dense_hash_set",
    )

    JLLWrappers.@init_file_product(
        densehashtable_h,
        "include/sparsehash/internal/densehashtable.h",
    )

    JLLWrappers.@init_file_product(
        hashtable_common_h,
        "include/sparsehash/internal/hashtable-common.h",
    )

    JLLWrappers.@init_file_product(
        libc_allocator_with_realloc_h,
        "include/sparsehash/internal/libc_allocator_with_realloc.h",
    )

    JLLWrappers.@init_file_product(
        sparse_hash_map,
        "include/sparsehash/sparse_hash_map",
    )

    JLLWrappers.@init_file_product(
        sparse_hash_set,
        "include/sparsehash/sparse_hash_set",
    )

    JLLWrappers.@init_file_product(
        sparsehashtable_h,
        "include/sparsehash/internal/sparsehashtable.h",
    )

    JLLWrappers.@init_file_product(
        sparsetable,
        "include/sparsehash/sparsetable",
    )

    JLLWrappers.@init_file_product(
        traits,
        "include/sparsehash/traits",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# This file is part of RLinearAlgebra.jl

using Test

@testset verbose=true "RLinearAlbera.jl" begin
    for file in readlines(joinpath(@__DIR__, "testgroups.txt"))
        include(file * ".jl")
    end
end


text:
module titus
using Random

mersenne = MersenneTwister(12346)
# [(char, weight)]
sparkle_chars = ['*', '*', '*', '~', '~', '~', '~', '~', '(',  ')']

# he's a dog, so food emoji right?
food_emoji = ["🍏", "🍎", "🍐", "🍊", "🍋", "🍌", "🍉", "🍇", "🍓", "🍈", "🍒", "🍑", "🥭", "🍍",
							"🥥", "🥝", "🍅", "🍆", "🥑", "🥦", "🥬", "🥒", "🌽", "🥕", "🧄", "🧅", "🥔", "🍠",
							"🥐", "🥯", "🍞", "🥖", "🥨", "🧀", "🥚", "🍳", "🧈", "🥞", "🧇", "🥓", "🥩", "🍗",
							"🍖", "🦴", "🌭", "🍔", "🍟", "🍕", "🥪", "🥙", "🧆", "🌮", "🌯", "🥗", "🥘", "🥫",
							"🍝", "🍜", "🍲", "🍛", "🍣", "🍱", "🥟", "🦪", "🍤", "🍙", "🍚", "🍘", "🍥", "🥠",
							"🥮", "🍢", "🍡", "🍧", "🍨", "🍦", "🥧", "🧁", "🍰", "🎂", "🍮", "🍭", "🍬", "🍫",
							"🍿", "🍩", "🍪", "🌰", "🥜", "🍯", "🥛", "🍼", "☕", "🍵", "🧃", "🥤", "🍶", "🍺",
							"🍻", "🥂", "🍷", "🥃", "🍸", "🍹", "🧉", "🍾", "🧊", "🥄", "🍴", "🥣", "🥡", "🥢",
							"🧂"]

function random_bool()
	n = rand(1)[1]
	n > 0.5
end

function randomize_upper(s)
	k = ""
	for char in s
		if random_bool()
			k *= uppercase(char)
		else
			k *= char
		end
	end
	return k
end

function rand_between(lower, upper)
	n = floor(rand(1)[1] * upper)
  while n > upper || n < lower
		n = floor(rand(1)[1] * upper)
	end
	return n
end

function get_random(arr)
	rand(mersenne, arr, 1)[1]
end

function random_emoji(str)
	period = r"\."
	s = str
	while !isnothing(match(period, s))
		replacement = get_random(food_emoji)
		s = replace(s, "."=>replacement, count=1)
	end
	s
end

function sparklez()
	sparkle_length = rand_between(5, 11)
	sparkle = ""
	while sparkle_length > 0
		sparkle *= get_random(sparkle_chars)
		sparkle_length -= 1
	end
	sparkle
end

println(map(x->sparklez()*random_emoji(randomize_upper(x))*sparklez(), ARGS))

end # module

max_stars_repo_name: pcapel/talk-like-titus

max_issues_repo_name: pcapel/talk-like-titus

max_forks_repo_name: pcapel/talk-like-titus


text:
using Photometry
using BenchmarkTools
using CSV
using DataFrames
using ProgressLogging

AP = [
    CircularAperture,
    CircularAnnulus,
    EllipticalAperture,
    EllipticalAnnulus,
    RectangularAperture,
    RectangularAnnulus,
]
PARAMS = [
    (3),
    (3, 5),
    (3, 3, 0),
    (3, 5, 4, 0),
    (3, 5, 0),
    (3, 5, 4, 0)
]

data = randn(512, 512) .+ 10

rows = []

@progress for N in [1, 10, 50, 100], method in [:exact, :center, (:subpixel, 5)], (A, par) in zip(AP, PARAMS)
    aper = A(256, 256, par...)
    apers = repeat([aper], N)
    func(a, d) = aperture_photometry(a, d, method=method)
    time = @belapsed $func($apers, $data)
    m = method isa Tuple ? "$(String(method[1]))-$(method[2])" : String(method)
    push!(rows, (N=N, aperture=A, method=m, time=time))
end

results = DataFrame(rows)
CSV.write(joinpath(@__DIR__, "julia_apertures.csv"), results)


text:
function eyler(n)
    s = 0.0
    a = 1.0
    for k in 1:n+1
    # границы значений индекса k по отношению к прежнему варианту теперь смещены на 1 вправо
        s += a
        a /= k
    end
    s
end # сложность O(n)

function eyler_2()
    s = 0
    a = 1
    k = 0
    while s + a != s
        k += 1
        s += a
        a /= k
    end
    return s
end
#при вычислении `s+a` должны будут складываться соответствующие 53-х битные
# мантиссы (здесь речь идет о формате типа Float64). Но сначала порядок меньшего
# слагаемого должен быть приведен к порядку большего слагаемого путем "сдвига" его
#мантиссы вправо с заменой "освобождающихся" битов нулями. Таким образом, если сдвиг
#составит 53 бита или больше, то мантисса меньшего числа окажется полностью замененной нулями.

function sin_(x)
    a = x
    s = 0
    e = 0.000000001
    k = 1
    while abs(a) > e
        s += a
        k += 1
        a *= -x*x/(2*k-1)/(2*k-2)
    end
    s
end

function sin_(x)
    xx=x^2
    a=x
    m=2
    s=typeof(x)(0)
    while s+a != s
        s+=a
        a=-a*xx/m/(m+1)
        m+=2
    end
    #УТВ: |sin(x)-s| <= 1е-16
    return s
end
#цикл завершится, как только разность двоитчных порядков величин `s` и `a` достигнет
# `53` (что примерно соответствует шестнадцати десятичный порядков).

function exp_(x, eps)
    sum = 1;
    f = x
    k = 2
    while f > eps
        sum += f
        f *= x/k
        k += 1
    end
    sum
end

function harmonic_sum()
    s=0.0
    k=1
    a=1.0
    while s+a != s
        a=1/k
        s+=a
        k+=1
    end
    return s
end

function cos_(x, n)
    mx2 = -x*x
    part = 1
    sum = 1
    for i = 2:2:n
        part *= mx2/((i-1)*i)
        sum += part
    end
    sum
end

function cos_m(x, e = 0.0000000000001)
    s = 0
    an = 1
    n = 0
    while abs(an) > e
        s += an
        n += 1
        an *= -x*x/(2*n-1)/(2*n)
    end
    s
end

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

using Plots

#
x = 0:0.1:10
plot(x, cos_.(x, 4))
#
p = plot() # получена ссылка на пустой график
for n in 2:2:8
    plot!(p, x, cos_.(x,n)) # в текущий график добавлена новая кривая
end
display(p)
#

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
using Test
using SPCSpectra

@testset "SPCSpectra" begin
    dir = joinpath(pkgdir(SPCSpectra), "test", "data")
    filename = "4d_map.spc"
    path = joinpath(dir, filename)
    spc = SPC(path)
    # The data is from a Nicolet FT-IR spectrometer.
    # That is, a Fourier infrared spectrometer.
    @test spc.param_dict["SRC"] == "IR Source"
    @test spc.param_dict["MODEL"] == "Nicolet"
end


text:
function count_questions(group_lines)
    return length(Set(replace(group_lines, '\n' => "")))
end

function run(s)
    # Your code here
    return sum(count_questions(group_lines) for group_lines in split(s, "\n\n"))
end

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

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()


text:
#Convert number in [0,1)^2 to [0,1)
function xy2d(x::Float64, y:: Float64)
    N = 2^30
    xint = floor(Int64, x*N)
    yint = floor(Int64, y*N)
    # d will be in [0,N^2)
    d = xy2d(N, xint, yint)/N^2
end

function d2xy(d:: Float64)
    N::Int64 = 2^30
    dint = floor(Int64, d*N^2)
    x::Float64, y::Float64 = d2xy(N^2, dint)
    x/N, y/N
end

#Converted from wikipedia
#convert (x,y) to d
function xy2d(n::Integer, x::Integer, y::Integer)
    rx, ry, s, d = 0, 0, 0, 0
    s = div(n,2)
    while s>0
        rx = (x & s) > 0
        ry = (y & s) > 0
        d += s * s * ((3*rx) $ ry)
        x, y = rot(s, x, y, rx, ry)
        s = div(s,2)
    end
    return d::Integer
end

#convert d to (x,y)
function d2xy(n::Integer, d::Integer)
    rx, ry, s, t = 0, 0, 0, d
    x, y = 0, 0
    s = 1
    while s<n
        rx = 1 & (div(t,2))
        ry = 1 & (t $ rx)
        x, y = rot(s, x, y, rx, ry)
        x += s*rx
        y += s*ry
        t = div(t,4)
        s *= 2
    end
    return x::Integer, y::Integer
end

#rotate/flip a quadrant appropriately
function rot(n::Integer, x::Integer, y::Integer, rx::Integer, ry::Integer)
    if ry == 0
        if rx == 1
            x = n-1 - x
            y = n-1 - y
        end
        #Swap x and y
        t = x
        x = y
        y = t
    end
    return x::Integer, y::Integer
end


text:
"""The module provides signals, a mechanism that models data retrievable from or
written to CAN-bus data. A signal models one data entity, e.g., one variable inside the
CAN-bus data.
"""
module Signals
import Base

"""
"""
abstract type AbstractSignal{T} end

"""
"""
abstract type UnnamedSignal{T} <: AbstractSignal{T} end

"""
"""
abstract type AbstractIntegerSignal{T <: Integer} <: UnnamedSignal{T} end

"""
"""
abstract type AbstractFloatSignal{T <: AbstractFloat} <: UnnamedSignal{T} end

"""
"""
struct Bit <: AbstractIntegerSignal{Bool}
    start::UInt16
end

"""
"""
function Bit(start::Integer)
    start = convert(UInt16, start)
    return Bit(start)
end

"""
"""
function Bit(; start::Integer=0)
    return Bit(start)
end

"""
"""
function start(signal::Bit)::UInt16
    return signal.start
end

"""
"""
function Base.length(signal::Bit)::UInt16
    return 1
end

"""
"""
function byte_order(signal::Bit)::Symbol
    return :little_endian
end

"""
"""
function Base.:(==)(lhs::Bit, rhs::Bit)
    return start(lhs) == start(rhs)
end

"""
"""
struct Unsigned{T} <: AbstractFloatSignal{T}
    start::UInt16
    length::UInt16
    factor::T
    offset::T
    byte_order::Symbol

    """
    """
    function Unsigned(start::UInt16,
                      length::UInt16,
                      factor::T,
                      offset::T,
                      byte_order::Symbol) where {T <: AbstractFloat}
          if byte_order != :little_endian && byte_order != :big_endian
              throw(DomainError(byte_order, "Byte order not supported"))
          end

          if length == 0
              throw(DomainError(length, "The length has to be greater or equal to 1"))
          end

          return new{T}(start, length, factor, offset, byte_order)
    end
end

"""
"""
function Unsigned(start::Integer,
                  length::Integer,
                  factor::Union{Integer, AbstractFloat},
                  offset::Union{Integer, AbstractFloat},
                  byte_order::Symbol)
    start = convert(UInt16, start)
    length = convert(UInt16, length)
    if factor isa Integer && offset isa Integer
        factor = convert(Float64, factor)
        offset = convert(Float64, offset)
    else
        factor, offset = promote(factor, offset)
    end
    return Unsigned(start, length, factor, offset, byte_order)
end

"""
"""
function Unsigned(; start::Integer,
                    length::Integer,
                    factor::Union{Integer, AbstractFloat},
                    offset::Union{Integer, AbstractFloat},
                    byte_order::Symbol=:little_endian)
    return Unsigned(start, length, factor, offset, byte_order)
end

"""
"""
function Unsigned{T}(start::Integer,
                length::Integer;
                factor::Union{Integer, AbstractFloat}=one(T),
                offset::Union{Integer, AbstractFloat}=zero(T),
                byte_order::Symbol=:little_endian) where {T}
    factor = convert(T, factor)
    offset = convert(T, offset)
    return Unsigned(start, length, factor, offset, byte_order)
end

"""
"""
function Unsigned{T}(; start::Integer,
                       length::Integer,
                       factor::Union{Integer, AbstractFloat}=one(T),
                       offset::Union{Integer, AbstractFloat}=zero(T),
                       byte_order::Symbol=:little_endian) where {T}
    factor = convert(T, factor)
    offset = convert(T, offset)
    return Unsigned(start, length, factor, offset, byte_order)
end

"""
"""
function start(signal::Unsigned{T})::UInt16 where {T}
    return signal.start
end

"""
"""
function Base.length(signal::Unsigned{T})::UInt16 where {T}
    return signal.length
end

"""
"""
function factor(signal::Unsigned{T})::T where {T}
    return signal.factor
end

"""
"""
function offset(signal::Unsigned{T})::T where {T}
    return signal.offset
end

"""
"""
function byte_order(signal::Unsigned{T})::Symbol where {T}
    return signal.byte_order
end

"""
"""
function Base.:(==)(lhs::F, rhs::F) where {T, F <: AbstractFloatSignal{T}}
    if start(lhs) != start(rhs)
        return false
    end

    if length(lhs) != length(rhs)
        return false
    end

    if factor(lhs) != factor(rhs)
        return false
    end

    if offset(lhs) != offset(rhs)
        return false
    end

    if byte_order(lhs) != byte_order(rhs)
        return false
    end

    return true
end

struct Signed{T} <: AbstractFloatSignal{T}
    start::UInt16
    length::UInt16
    factor::T
    offset::T
    byte_order::Symbol

    function Signed(start::UInt16,
                    length::UInt16,
                    factor::T,
                    offset::T,
                    byte_order::Symbol) where {T <: AbstractFloat}
        if byte_order != :little_endian && byte_order != :big_endian
            throw(DomainError(byte_order, "Byte order not supported"))
        end

        if length == 0
            throw(DomainError(length, "The length has to be greater or equal to 1"))
        end

        return new{T}(start, length, factor, offset, byte_order)
    end
end

"""
"""
function Signed(start::Integer,
                length::Integer,
                factor::Union{Integer, AbstractFloat},
                offset::Union{Integer, AbstractFloat},
                byte_order::Symbol)
    start = convert(UInt16, start)
    length = convert(UInt16, length)
    if factor isa Integer && offset isa Integer
        factor = convert(Float64, factor)
        offset = convert(Float64, offset)
    else
        factor, offset = promote(factor, offset)
    end
    return Signed(start, length, factor, offset, byte_order)
end

"""
"""
function Signed(; start::Integer,
                  length::Integer,
                  factor::Union{Integer, AbstractFloat},
                  offset::Union{Integer, AbstractFloat},
                  byte_order::Symbol=:little_endian)
    return Signed(start, length, factor, offset, byte_order)
end

"""
"""
function Signed{T}(start::Integer,
                   length::Integer;
                   factor::Union{Integer, AbstractFloat}=one(T),
                   offset::Union{Integer, AbstractFloat}=zero(T),
                   byte_order::Symbol=:little_endian) where {T}
    factor = convert(T, factor)
    offset = convert(T, offset)
    return Signed(start, length, factor, offset, byte_order)
end

"""
"""
function Signed{T}(; start::Integer,
                     length::Integer,
                     factor::Union{Integer, AbstractFloat}=one(T),
                     offset::Union{Integer, AbstractFloat}=zero(T),
                     byte_order::Symbol=:little_endian) where {T}
    factor = convert(T, factor)
    offset = convert(T, offset)
    return Signed(start, length, factor, offset, byte_order)
end

"""
"""
function start(signal::Signed{T})::UInt16 where {T}
    return signal.start
end

"""
"""
function Base.length(signal::Signed{T})::UInt16 where {T}
    return signal.length
end

"""
"""
function factor(signal::Signed{T})::T where {T}
    return signal.factor
end

"""
"""
function offset(signal::Signed{T})::T where {T}
    return signal.offset
end

"""
"""
function byte_order(signal::Signed{T})::Symbol where {T}
    return signal.byte_order
end

"""
"""
struct FloatSignal{T} <: AbstractFloatSignal{T}
    start::UInt16
    factor::T
    offset::T
    byte_order::Symbol

    function FloatSignal(start::UInt16, factor::T, offset::T,
                         byte_order::Symbol) where {T <: AbstractFloat}
        new{T}(start, factor, offset, byte_order)
    end
end

"""
"""
function FloatSignal(start::Integer,
                     factor::Union{Integer,AbstractFloat},
                     offset::Union{Integer,AbstractFloat},
                     byte_order::Symbol)
    start = convert(UInt16, start)
    if factor isa Integer && offset isa Integer
        factor = convert(Float64, factor)
        offset = convert(Float64, offset)
    else
        factor, offset = promote(factor, offset)
    end
    return FloatSignal(start, factor, offset, byte_order)
end

"""
"""
function FloatSignal(; start::Integer,
                       factor::Union{Integer,AbstractFloat},
                       offset::Union{Integer,AbstractFloat},
                       byte_order::Symbol)
    return FloatSignal(start, factor, offset, byte_order)
end

"""
"""
function FloatSignal{T}(start::Integer;
                        factor::Union{Integer,AbstractFloat}=one(T),
                        offset::Union{Integer,AbstractFloat}=zero(T),
                        byte_order::Symbol=:little_endian) where {T}
    factor = convert(T, factor)
    offset = convert(T, offset)
    return FloatSignal(start, factor, offset, byte_order)
end

function FloatSignal{T}(; start::Integer,
                          factor::Union{Integer,AbstractFloat}=one(T),
                          offset::Union{Integer,AbstractFloat}=zero(T),
                          byte_order::Symbol=:little_endian) where {T}
    factor = convert(T, factor)
    offset = convert(T, offset)
    return FloatSignal(start, factor, offset, byte_order)
end

const Float16Signal = FloatSignal{Float16}
const Float32Signal = FloatSignal{Float32}
const Float64Signal = FloatSignal{Float64}

"""
"""
function start(signal::FloatSignal{T})::UInt16 where {T}
    return signal.start
end

function Base.length(signal::FloatSignal{T})::UInt16 where {T}
    return 8sizeof(T)
end

"""
"""
function factor(signal::FloatSignal{T})::T where {T}
    return signal.factor
end

"""
"""
function offset(signal::FloatSignal{T})::T where {T}
    return signal.offset
end

"""
"""
function byte_order(signal::FloatSignal{T})::Symbol where {T}
    return signal.byte_order
end

"""
"""
struct Raw <: AbstractIntegerSignal{UInt64}
    start::UInt16
    length::UInt16
    byte_order::Symbol

    """
    """
    function Raw(start::UInt16, length::UInt16,byte_order::Symbol)
        if length == 0
            throw(DomainError(length, "The length has to be greater or equal to 1"))
        end

        if byte_order != :little_endian && byte_order != :big_endian
            throw(DomainError(byte_order, "Byte order not supported"))
        end

        return new(start, length, byte_order)
    end
end

"""
"""
function Raw(start::Integer, length::Integer, byte_order::Symbol)
    start = convert(UInt16, start)
    length = convert(UInt16, length)
    return Raw(start, length, byte_order)
end

"""
"""
function Raw(; start::Integer,
               length::Integer,
               byte_order::Symbol=:little_endian) where {T}
    return Raw(start, length, byte_order)
end

"""
"""
function start(signal::Raw)::UInt16
    return signal.start
end

"""
"""
function Base.length(signal::Raw)::UInt16
    return signal.length
end

"""
"""
function byte_order(signal::Raw)::Symbol
    return signal.byte_order
end

"""
"""
struct NamedSignal{T} <: AbstractSignal{T}
    name::String
    unit::Union{Nothing,String}
    default::Union{Nothing,T}
    signal::UnnamedSignal{T}

    function NamedSignal(name::String,
                         unit::Union{Nothing,String},
                         default::Union{Nothing,T},
                         signal::UnnamedSignal{T}) where {T}
        if name == ""
            throw(DomainError(name, "name cannot be the empty string"))
        end

        return new{T}(name, unit, default, signal)
    end
end

"""
"""
function NamedSignal(; name::String,
                     unit::Union{Nothing,String}=nothing,
                     default::Union{Nothing,T}=nothing,
                     signal::UnnamedSignal{T}) where {T}
    return NamedSignal(name, unit, default, signal)
end

"""
"""
function name(signal::NamedSignal{T})::String where {T}
    return signal.name
end

"""
"""
function unit(signal::NamedSignal{T})::Union{Nothing,String} where {T}
    return signal.unit
end

"""
"""
function default(signal::NamedSignal{T})::Union{Nothing,T} where {T}
    return signal.default
end

"""
"""
function signal(signal::NamedSignal{T})::UnnamedSignal{T} where {T}
    return signal.signal
end

const Signal = NamedSignal

"""
"""
struct Bits
    bits::Set{UInt16}
end

"""
"""
function Bits(bits::Integer...)
    Bits(Set(UInt16[bits...]))
end

"""
"""
function Bits(signal::AbstractFloatSignal{T}) where {T}
    bits = Set{UInt16}()
    start_bit = start(signal)
    if byte_order(signal) == :little_endian
        for i=0:length(signal)-1
            bit_pos = start_bit + i
            push!(bits, bit_pos)
        end
    elseif byte_order(signal) == :big_endian
        for j=0:length(signal)-1
            push!(bits, start_bit)
            if start_bit % 8 == 0
                start_byte = div(start_bit,8)
                start_bit = 8 * (start_byte + 1) + 7
            else
                start_bit -= 1
            end
        end
    end
    return Bits(bits)
end

"""
"""
function Bits(signal::AbstractIntegerSignal{T}) where {T}
    bits = Set{UInt16}()
    start_bit = start(signal)
    if byte_order(signal) == :little_endian
        for i=0:length(signal)-1
            bit_pos = start_bit + i
            push!(bits, bit_pos)
        end
    elseif byte_order(signal) == :big_endian
        for j=0:length(signal)-1
            push!(bits, start_bit)
            if start_bit % 8 == 0
                start_byte = div(start_bit,8)
                start_bit = 8 * (start_byte + 1) + 7
            else
                start_bit -= 1
            end
        end
    end
    return Bits(bits)
end

function Bits(sig::NamedSignal{T}) where {T}
    return Bits(signal(sig))
end

function Base.:(==)(lhs::Bits, rhs::Bits)::Bool
    return (lhs.bits) == (rhs.bits)
end

"""
"""
function share_bits(lhs::Bits, rhs::Bits)::Bool
    return !isdisjoint(lhs.bits, rhs.bits)
end

"""
"""
function overlap(lhs::AbstractSignal{R}, rhs::AbstractSignal{S})::Bool where {R,S}
    lhs_bits = Bits(lhs)
    rhs_bits = Bits(rhs)
    return share_bits(lhs_bits, rhs_bits)
end

function check(signal::UnnamedSignal{T}, available_bytes::UInt8)::Bool where {T}
    bits = Bits(signal)
    max_byte = max(UInt8[div(bit,8) for bit in bits.bits]...)
    if max_byte < available_bytes
        return true
    else
        return false
    end
end

"""
"""
function check(sig::NamedSignal{T}, available_bytes::UInt8)::Bool where {T}
    return check(signal(sig), available_bytes)
end

export Bit, Unsigned, Signed, Raw, Float16Signal, Float32Signal, Float64Signal
export Signal, FloatSignal
export NamedSignal
export start, factor, offset, byte_order
export name, unit, default, signal
end


text:
struct RadarDataJPL{T}
    object::String
    utcepoch::DateTime
    delay::T
    delay_sigma::T
    delay_units::String
    doppler::T
    doppler_sigma::T
    doppler_units::String
    freq::T
    rcvr::Int
    xmit::Int
    bouncepoint::String
    #Inner constructor
    function RadarDataJPL{T}(
        object::String,
        utcepoch::DateTime,
        delay::T,
        delay_sigma::T,
        delay_units::String,
        doppler::T,
        doppler_sigma::T,
        doppler_units::String,
        freq::T,
        rcvr::Int,
        xmit::Int,
        bouncepoint::String
    ) where {T<:Number}
        return new{T}(object, utcepoch, delay, delay_sigma, delay_units,
            doppler, doppler_sigma, doppler_units, freq, rcvr, xmit, bouncepoint)
    end
end

# Outer constructors
RadarDataJPL(r::RadarDataJPL{T}) where {T<:Number} = r
function RadarDataJPL(object::String,
    utcepoch::DateTime,
    delay::T,
    delay_sigma::T,
    delay_units::String,
    doppler::T,
    doppler_sigma::T,
    doppler_units::String,
    freq::T,
    rcvr::Int,
    xmit::Int,
    bouncepoint::String
) where {T<:Number}
    return RadarDataJPL{T}(object, utcepoch, delay, delay_sigma, delay_units,
        doppler, doppler_sigma, doppler_units, freq, rcvr, xmit, bouncepoint)
end

utcepoch(rdata::RadarDataJPL) = rdata.utcepoch
delay(rdata::RadarDataJPL) = rdata.delay
delay_sigma(rdata::RadarDataJPL) = rdata.delay_sigma
delay_units(rdata::RadarDataJPL) = rdata.delay_units
doppler(rdata::RadarDataJPL) = rdata.doppler
doppler_sigma(rdata::RadarDataJPL) = rdata.doppler_sigma
doppler_units(rdata::RadarDataJPL) = rdata.doppler_units
freq(rdata::RadarDataJPL) = rdata.freq
rcvr(rdata::RadarDataJPL) = rdata.rcvr
xmit(rdata::RadarDataJPL) = rdata.xmit
bouncepoint(rdata::RadarDataJPL) = rdata.bouncepoint

#monostatic mode: check that each receiver and transmitter are the same
ismonostatic(rdata::RadarDataJPL) = rdata.rcvr == rdata.xmit

# check whether a given transmission frequency (in MHz) belongs to the S band (IEEE nomenclature)
issband(f_MHz::T) where {T<:Real} = 2000.0 ≤ f_MHz ≤ 4000.0

# check whether a given transmission frequency (in MHz) belongs to the X band (IEEE nomenclature)
isxband(f_MHz::T) where {T<:Real} = 8000.0 ≤ f_MHz ≤ 12000.0

function process_radar_data_jpl(radar_data_jpl_file_path::String)
    radar_data_jpl = readdlm(radar_data_jpl_file_path, '\t')
    dataformat_jpl = "y-m-d H:M:S"
    jpl_datetime = DateTime.(radar_data_jpl[:,2], dataformat_jpl)
    del_ind = findall(x->x=="us", radar_data_jpl[:,5])
    dop_ind = findall(x->x=="Hz", radar_data_jpl[:,5])
    jpl_datetime_unrep = sort(
    union(
        DateTime.(radar_data_jpl[del_ind,2], dataformat_jpl),
        DateTime.(radar_data_jpl[dop_ind,2], dataformat_jpl)
        )
    )
    radar_data_jpl_processed = Array{Any}(undef, length(jpl_datetime_unrep), 12)
    for i in eachindex(jpl_datetime_unrep)
        inds = findall(x->x==jpl_datetime_unrep[i], jpl_datetime)
        if length(inds) == 1
            if radar_data_jpl[inds[1],5] == "us"
                x = vcat(radar_data_jpl[inds[1],1], jpl_datetime_unrep[i], radar_data_jpl[inds[1],3:5], [NaN, NaN, "X"], radar_data_jpl[inds[1], 6:end])
            elseif radar_data_jpl[inds[1],5] == "Hz"
                x = vcat(radar_data_jpl[inds[1],1], jpl_datetime_unrep[i], [NaN, NaN, "X"], radar_data_jpl[inds[1],3:5], radar_data_jpl[inds[1], 6:end])
            else
                @warn "Measurement units not us nor Hz"
            end
            radar_data_jpl_processed[i,:] .= x
        elseif length(inds) == 2
            if radar_data_jpl[inds[1],5] == "us"
                y = vcat(radar_data_jpl[inds[1],1], jpl_datetime_unrep[i], radar_data_jpl[inds[1],3:5], radar_data_jpl[inds[2],3:5], radar_data_jpl[inds[1], 6:end])
            elseif radar_data_jpl[inds[1],5] == "Hz"
                y = vcat(radar_data_jpl[inds[1],1], jpl_datetime_unrep[i], radar_data_jpl[inds[2],3:5], radar_data_jpl[inds[1],3:5], radar_data_jpl[inds[1], 6:end])
            else
                @warn "Measurement units not us nor Hz"
            end
            radar_data_jpl_processed[i,:] .= y
        else
            @warn "length(inds) != 1 && length(inds) != 2"
        end
    end
    radar_data_vec = RadarDataJPL.(
        String.(radar_data_jpl_processed[:,1]),
        DateTime.(radar_data_jpl_processed[:,2]),
        Float64.(radar_data_jpl_processed[:,3]),
        Float64.(radar_data_jpl_processed[:,4]),
        String.(radar_data_jpl_processed[:,5]),
        Float64.(radar_data_jpl_processed[:,6]),
        Float64.(radar_data_jpl_processed[:,7]),
        String.(radar_data_jpl_processed[:,8]),
        Float64.(radar_data_jpl_processed[:,9]),
        Int.(radar_data_jpl_processed[:,10]),
        Int.(radar_data_jpl_processed[:,11]),
        String.(radar_data_jpl_processed[:,12])
    )
    return radar_data_vec
end


text:
@testset "moments" begin

# Build histogram from an image correctly and compare the threshold to the answer
# from another implementation.
img = testimage("cameraman")
imgg = gray.(img)
imgg_uint8 = reinterpret.(imgg)
edges, counts = build_histogram(imgg_uint8, 0:1:255)
t1 = find_threshold(Moments(),counts[1:end], 0:1:255)
@test t1 == 111

# Build a histogram from an array and compare the threshold to that in another
# implementation.
arr = [10  8 10  9 20 21 32 30 40 41 41 40;
       12 10 11 10 19 20 30 28 38 40 40 39;
       10  9 10  8 20 21 30 29 42 40 40 39;
       11 10  9 11 19 21 31 30 40 42 38 40]
edges, counts = build_histogram(arr, 0:255)
t2 = find_threshold(Moments(), counts[1:end], edges)
@test t2 == 27

end

max_stars_repo_name: AbbeyMac/HistogramThresholding

max_issues_repo_name: AbbeyMac/HistogramThresholding

max_forks_repo_name: AbbeyMac/HistogramThresholding


text:
import Base: length, getindex
export SampleSet, MetricNN, QuasiMetricNN
export NearNeighborCache, MutableNNC, ImmutableNNC
export ControlCache, EmptyControlCache, MutableControlCache, ImmutableControlCache
export DistanceDataStructure, EmptyDistanceDS, BruteDistanceDS, TreeDistanceDS
export filter_neighborhood, addpoints
# export resetcache!
# export saveNN, loadNN!
export inball, knn, inball!, knn!, mutualknn!         # for metrics only (symmetric)
export inballF, knnF, inballF!, knnF!, mutualknnF!    # forwards
export inballB, knnB, inballB!, knnB!, mutualknnB!    # backwards

### Typedefs
abstract SampleSet
abstract NearNeighborCache{T}; begin
    ## NearNeighborCache implementations
    # TODO: consider EmptyNNC for online algorithms or streaming-batch
    immutable MutableNNC{T} <: NearNeighborCache{T}
        D::Vector{SparseVector{T,Int}}
        r::Vector{T}
    end
    MutableNNC{T}(N::Int, ::Type{T}) = MutableNNC{T}(Array(SparseVector{T,Int}, N), zeros(T, N))
    immutable ImmutableNNC{T} <: NearNeighborCache{T}
        D::SparseMatrixCSC{T,Int}
        r::Vector{T}
    end
    ImmutableNNC(M::MutableNNC) = ImmutableNNC(hcat(M.D...), M.r)
    MutableNNC{T}(I::ImmutableNNC{T}) = MutableNNC(size(I.D,1), T) # clear the cache
end

abstract ControlCache{U<:ControlInfo}; begin
    ## ControlCache implementations
    immutable EmptyControlCache{U} <: ControlCache{U} end
    EmptyControlCache{U}(::Type{U}) = EmptyControlCache{U}()
    EmptyControlCache() = EmptyControlCache(NullControl)
    immutable MutableControlCache{U} <: ControlCache{U}
        US::Vector{SparseVector{U,Int}}
    end
    MutableControlCache{U}(N::Int, ::Type{U}) = MutableControlCache{U}(Array(SparseVector{U,Int}, N))
    immutable ImmutableControlCache{U} <: ControlCache{U}
        US::SparseMatrixCSC{U,Int}
    end
    ImmutableControlCache(M::MutableControlCache) = ImmutableControlCache(hcat(M.U...))
    MutableControlCache{U}(I::ImmutableControlCache{U}) = MutableControlCache(size(I.US,1), U) # clear the cache
end

abstract DistanceDataStructure{T}; begin
    ## DistanceDataStructure implementations
    immutable EmptyDistanceDS{T} <: DistanceDataStructure{T} end
    EmptyDistanceDS{T}(::Type{T}) = EmptyDistanceDS{T}()
    immutable BruteDistanceDS{A<:AbstractMatrix,T} <: DistanceDataStructure{T}
        Dmat::A
    end
    BruteDistanceDS{T}(Dmat::AbstractMatrix{T}) = BruteDistanceDS{typeof(Dmat),T}(Dmat)
    immutable TreeDistanceDS{NNT<:NearestNeighbors.NNTree,T} <: DistanceDataStructure{T}
        tree::NNT
    end
    TreeDistanceDS{V}(tree::NearestNeighbors.NNTree{V}) = TreeDistanceDS{typeof(tree),eltype(V)}(tree)
end

### Construction
type MetricNN{M<:SymmetricDistance,S<:State,N<:NearNeighborCache,D<:DistanceDataStructure,U<:ControlCache} <: SampleSet
    V::Vector{S}
    dist::M
    init::S
    cache::N
    DS::D
    US::U
end
function MetricNN{S<:State}(V::Vector{S}, dist::SymmetricDistance, init::S)
    N = length(V)
    T = eltype(S)
    MetricNN(V, dist, init, MutableNNC(N,T), helper_data_structures(V, dist)...)
end

type QuasiMetricNN{M<:AsymmetricDistance,S<:State,N<:NearNeighborCache,
                   D<:DistanceDataStructure,U<:ControlCache} <: SampleSet
    V::Vector{S}
    dist::M
    init::S
    cacheF::N   # Everything's just pointers, but F and B may be different (e.g. CSC and CSR), trading memory for speed
    cacheB::N
    DSF::D
    USF::U
    DSB::D
    USB::U
end
function QuasiMetricNN{S<:State}(V::Vector{S}, dist::AsymmetricDistance, init::S)
    N = length(V)
    T = eltype(S)
    QuasiMetricNN(V, dist, init, MutableNNC(N,T), MutableNNC(N,T), helper_data_structures(V, dist)...)
end

## Fallbacks (should usually be extended)
helper_data_structures{S}(V::Vector{S}, d::SymmetricDistance) = (EmptyDistanceDS(eltype(S)),
                                                                 EmptyControlCache(controltype(d)))
function helper_data_structures{S}(V::Vector{S}, d::AsymmetricDistance)
    DS, US = EmptyDistanceDS(eltype(S)), EmptyControlCache(controltype(d))
    (DS, US, DS, US)
end

### Utilities

@inline function filter_neighborhood(n::AbstractSparseVector, f::BitVector)
    inds, ds = nonzeroinds(n), nonzeros(n)
    SparseVector(length(n), inds[f[inds]], ds[f[inds]])
end
addpoints(NN::MetricNN, W) = MetricNN([NN.V; W], NN.dist, NN.init)
addpoints(NN::QuasiMetricNN, W) = QuasiMetricNN([NN.V; W], NN.dist, NN.init)
length(NN::SampleSet) = length(NN.V)
getindex(NN::SampleSet, i::Int) = i > 0 ? NN.V[i] : NN.init
getindex(NN::SampleSet, I::AbstractVector) = [NN[i] for i in I]

### Serialization

### Near Neighbor Computation
# TODO: forwards/backwards should be probably moved into the datastructure; perhaps a type parameterization
@inline inball!(NN::MetricNN, v::Int, r, f::BitVector) = filter_neighborhood(inball!(NN, v, r), f)
@inline inballF!(NN::QuasiMetricNN, v::Int, r, f::BitVector) = filter_neighborhood(inballF!(NN, v, r), f)
@inline inballB!(NN::QuasiMetricNN, v::Int, r, f::BitVector) = filter_neighborhood(inballB!(NN, v, r), f)
for direction in ("", "F", "B")
    NNsym = (direction == "" ? :MetricNN : :QuasiMetricNN)
    inballD! = Symbol(:inball, direction, "!")
    cacheD = Symbol(:cache, direction)
    DSD = Symbol(:DS, direction)
    @eval $inballD!{M,S,N<:ImmutableNNC}(NN::$NNsym{M,S,N}, v::Int, r) = viewcol(NN.$cacheD.D, v)
    @eval function $inballD!{M,S,N<:MutableNNC}(NN::$NNsym{M,S,N}, v::Int, r)
        if !isdefined(NN.$cacheD.D, v)
            NN.$cacheD.D[v] = inball(NN.V, NN.dist, NN.$DSD, v, r, $(direction != "B"))
            NN.$cacheD.r[v] = r
        end
        NN.$cacheD.D[v]
    end
end

function inball(V::Vector, dist::PreMetric, DS::DistanceDataStructure, v::Int, r, forwards::Bool = true)
    allds = forwards ? colwise(dist, V[v], V) : colwise(dist, V, V[v])
    N = length(V)
    inds = Vector{Int}(0)
    ds = similar(allds, 0)
    @inbounds for i in 1:N
        if i != v && allds[i] <= r
            push!(inds, i)
            push!(ds, allds[i])
        end
    end
    SparseVector(N, inds, ds)
end

function inball{A<:Matrix}(V::Vector, dist::PreMetric, DS::BruteDistanceDS{A}, v::Int, r, forwards::Bool = true)
    N = length(V)
    inds = Vector{Int}(0)
    ds = similar(DS.Dmat, 0)
    @inbounds for i in 1:N
        if i != v && DS.Dmat[i,v] <= r
            push!(inds, i)
            push!(ds, allds[i])
        end
    end
    SparseVector(N, inds, ds)
end

function inball{A<:SparseMatrixCSC}(V::Vector, dist::PreMetric, DS::BruteDistanceDS{A}, v::Int, r, forwards::Bool = true)
    vcol = viewcol(DS.Dmat, v)
    allinds, allds = nonzeroinds(vcol), nonzeros(vcol)
    inds = Vector{Int}(0)
    ds = Vector{eltype(A)}(0)
    @inbounds for i in 1:length(allinds)
        if allinds[i] != v && allds[i] <= r
            push!(inds, allinds[i])
            push!(ds, allds[i])
        end
    end
    SparseVector(length(V), inds, ds)
end

function inball{S}(V::Vector{S}, dist::PreMetric, DS::TreeDistanceDS, v::Int, r, forwards::Bool = true)
    inds = inrange(DS.tree, V[v], r, true)
    inds = deleteat!(inds, searchsortedfirst(inds, v))
    SparseVector(length(V), inds, forwards ? colwise(dist, V[v], V[inds]) : colwise(dist, V[inds], V[v]))
end

function inball{S}(V::Vector{S}, dist::ChoppedPreMetric, DS::TreeDistanceDS, v::Int, r, forwards::Bool = true)
    allinds = inrange(DS.tree, V[v], r, true)
    allinds = deleteat!(allinds, searchsortedfirst(allinds, v))
    allds = forwards ? colwise(dist, V[v], V[allinds]) : colwise(dist, V[allinds], V[v])
    inds = Vector{Int}(0)
    ds = similar(allds, 0)
    @inbounds for i in 1:length(allinds)
        if allds[i] <= r
            push!(inds, allinds[i])
            push!(ds, allds[i])
        end
    end
    SparseVector(length(V), inds, ds)
end

for ff in [x -> Symbol("inball", x), x -> Symbol("inball", x, "!")]
    @eval $(ff("F"))(NN::MetricNN, args...) = $(ff(""))(NN, args...)
    @eval $(ff("B"))(NN::MetricNN, args...) = $(ff(""))(NN, args...)
end


text:
const VERSION_ID_REGEX = r"^[0-9a-zA-Z\._]{32}$"

struct S3Path{A<:AbstractS3PathConfig} <: AbstractPath
    segments::Tuple{Vararg{String}}
    root::String
    drive::String
    isdirectory::Bool
    version::Union{String,Nothing}
    config::A

    # Inner constructor performs no data checking and is only meant for direct use by
    # deserialization.
    function S3Path{A}(
        segments, root, drive, isdirectory, version, config::A
    ) where {A<:AbstractS3PathConfig}
        return new(segments, root, drive, isdirectory, version, config)
    end
end

function S3Path(
    segments, root, drive, isdirectory, version::AbstractString, config::A
) where {A<:AbstractS3PathConfig}
    # Validate the `version` string provided is valid. Having this check during construction
    # allows us to fail early instead of having to wait to make an API call to fail.
    if !occursin(VERSION_ID_REGEX, version)
        throw(ArgumentError("`version` string is invalid: $(repr(version))"))
    end

    return S3Path{A}(segments, root, drive, isdirectory, version, config)
end

function S3Path(
    segments, root, drive, isdirectory, version::Nothing, config::A
) where {A<:AbstractS3PathConfig}
    return S3Path{A}(segments, root, drive, isdirectory, version, config)
end

"""
    S3Path()
    S3Path(str; version::$(AbstractS3Version)=nothing, config::$(AbstractS3PathConfig)=nothing)

Construct a new AWS S3 path type which should be of the form
`"s3://<bucket>/prefix/to/my/object"`.

NOTES:

- Directories are required to have a trailing `/` due to how S3
  distinguishes files from folders, as internally they're just
  keys to objects.
- Objects `p"s3://bucket/a"` and `p"s3://bucket/a/b"` can co-exist.
  If both of these objects exist listing the keys for `p"s3://bucket/a"` returns
  `[p"s3://bucket/a"]` while `p"s3://bucket/a/"` returns `[p"s3://bucket/a/b"]`.
- The drive property will return `"s3://<bucket>"`
- On top of the standard path properties (e.g., `segments`, `root`, `drive`,
  `separator`), `S3Path`s also support `bucket` and `key` properties for your
  convenience.
- If `version` argument is `nothing`, will return latest version of object. Version
  can be provided via either kwarg `version` or as suffix `"?versionId=<object_version>"`
  of `str`, e.g., `"s3://<bucket>/prefix/to/my/object?versionId=<object_version>"`.
- If `config` is left at its default value of `nothing`, then the
  latest `global_aws_config()` will be used in any operations involving the
  path. To "freeze" the config at construction time, explicitly pass an
  `AbstractAWSConfig` to the `config` keyword argument.
"""
S3Path() = S3Path((), "/", "", true, nothing, nothing)

# below definition needed by FilePathsBase
S3Path{A}() where {A<:AbstractS3PathConfig} = S3Path()

function S3Path(
    bucket::AbstractString,
    key::AbstractString;
    isdirectory::Bool=endswith(key, "/"),
    version::AbstractS3Version=nothing,
    config::AbstractS3PathConfig=nothing,
)
    return S3Path(
        Tuple(filter!(!isempty, split(key, "/"))),
        "/",
        strip(startswith(bucket, "s3://") ? bucket : "s3://$bucket", '/'),
        isdirectory,
        version,
        config,
    )
end

function S3Path(
    bucket::AbstractString,
    key::AbstractPath;
    isdirectory::Bool=false,
    version::AbstractS3Version=nothing,
    config::AbstractS3PathConfig=nothing,
)
    return S3Path(
        key.segments, "/", normalize_bucket_name(bucket), isdirectory, version, config
    )
end

# To avoid a breaking change.
function S3Path(
    str::AbstractString;
    isdirectory::Union{Bool,Nothing}=nothing,
    version::AbstractS3Version=nothing,
    config::AbstractS3PathConfig=nothing,
)
    result = tryparse(S3Path, str; config=config)
    result !== nothing || throw(ArgumentError("Invalid s3 path string: $str"))
    ver = if version !== nothing
        if result.version !== nothing && result.version != version
            throw(ArgumentError("Conflicting object versions in `version` and `str`"))
        end
        version
    else
        result.version
    end

    # Replace the parsed isdirectory field with an explicit passed in argument.
    is_dir = isdirectory === nothing ? result.isdirectory : isdirectory

    # Warning: We need to use the full constructor because reconstructing with the bucket
    # and key results in inconsistent `root` fields.
    return S3Path(result.segments, result.root, result.drive, is_dir, ver, result.config)
end

function Base.tryparse(
    ::Type{<:S3Path}, str::AbstractString; config::AbstractS3PathConfig=nothing
)
    uri = URI(str)
    uri.scheme == "s3" || return nothing

    drive = "s3://$(uri.host)"
    root = isempty(uri.path) ? "" : "/"
    isdirectory = isempty(uri.path) || endswith(uri.path, '/')
    path = Tuple(split(uri.path, '/'; keepempty=false))
    version = get(queryparams(uri), "versionId", nothing)

    return S3Path(path, root, drive, isdirectory, version, config)
end

function normalize_bucket_name(bucket)
    return strip(startswith(bucket, "s3://") ? bucket : "s3://$bucket", '/')
end

function Base.print(io::IO, fp::S3Path)
    print(io, fp.anchor, fp.key)
    fp.version !== nothing && print(io, "?versionId=", fp.version)
    return nothing
end

function Base.:(==)(a::S3Path, b::S3Path)
    return (
        a.segments == b.segments &&
        a.root == b.root &&
        a.drive == b.drive &&
        a.isdirectory == b.isdirectory &&
        a.version == b.version
    )
end

function Base.getproperty(fp::S3Path, attr::Symbol)
    if attr === :anchor
        return fp.drive * fp.root
    elseif attr === :separator
        return "/"
    elseif attr === :bucket
        return split(fp.drive, "//")[2]
    elseif attr === :key
        if isempty(fp.segments)
            return ""
        end

        return join(fp.segments, '/') * (fp.isdirectory ? "/" : "")
    else
        return getfield(fp, attr)
    end
end

# We need to special case join and parents so that we propagate
# directories correctly (see type docstring for details)
function FilePathsBase.join(prefix::S3Path, pieces::AbstractString...)
    isempty(pieces) && return prefix

    segments = String[prefix.segments...]
    isdirectory = endswith(last(pieces), "/")

    for p in pieces
        append!(segments, filter!(!isempty, split(p, "/")))
    end

    return S3Path(
        tuple(segments...),
        "/",
        prefix.drive,
        isdirectory,
        nothing, # Version is per-object, so we should not propagate it from the prefix
        prefix.config,
    )
end

function FilePathsBase.parents(fp::S3Path)
    if hasparent(fp)
        return map(0:(length(fp.segments) - 1)) do i
            S3Path(fp.segments[1:i], fp.root, fp.drive, true, nothing, fp.config)
        end
    elseif fp.segments == tuple(".") || isempty(fp.segments)
        return [fp]
    else
        return [isempty(fp.root) ? Path(fp, tuple(".")) : Path(fp, ())]
    end
end

"""
    get_config(fp::S3Path)

Returns the AWS configuration used by the path `fp`.  This can be stored within the path itself, but if not
it will be fetched with `AWS.global_aws_config()`.
"""
get_config(fp::S3Path) = @something(fp.config, global_aws_config())

function FilePathsBase.exists(fp::S3Path)
    return s3_exists(get_config(fp), fp.bucket, fp.key; version=fp.version)
end

Base.isfile(fp::S3Path) = !fp.isdirectory && exists(fp)
function Base.isdir(fp::S3Path)
    fp.isdirectory || return false
    if isempty(fp.segments)  # special handling of buckets themselves
        # may not be super efficient for those with billions of buckets, but really our best option
        fp.bucket ∈ s3_list_buckets(get_config(fp))
    else
        exists(fp)
    end
end

function FilePathsBase.walkpath(fp::S3Path; kwargs...)
    # Select objects with that prefix
    objects = s3_list_objects(get_config(fp), fp.bucket, fp.key; delimiter="")
    root = joinpath(fp, "/")

    # Construct a new Channel using a recursive internal `_walkpath!` function
    return Channel(; ctype=typeof(fp), csize=128) do chnl
        _walkpath!(root, root, Iterators.Stateful(objects), chnl; kwargs...)
    end
end

function _walkpath!(
    root::S3Path, prefix::S3Path, objects, chnl; topdown=true, onerror=throw, kwargs...
)
    @assert root.isdirectory
    @assert prefix.isdirectory

    while true
        try
            # Start by inspecting the next element
            obj = Base.peek(objects)

            # Early exit condition if we've exhausted the iterator or just the current prefix.
            obj === nothing && return nothing

            # Extract the non-root part of the key
            k = chop(obj["Key"]; head=length(root.key), tail=0)

            fp = joinpath(root, k)
            _parents = parents(fp)

            # If the filepath matches our prefix then pop it off and continue
            # Cause we would have already processed it before recursing
            child = if prefix.segments == fp.segments
                popfirst!(objects)
                continue
                # If the filpath is a direct descendant of our prefix then check if it
                # is a directory too
            elseif last(_parents) == prefix
                popfirst!(objects)
                # If our current path is a prefix for the next path then we can assume that
                # the current path should be a directory without needing to call `isdir`
                next = Base.peek(objects)
                is_dir =
                    (next !== nothing && startswith(next["Key"], fp.key)) || isdir(fp)
                # Reconstruct our next object and explicitly specify whether it is a
                # directory.
                S3Path(
                    fp.bucket,
                    fp.key;
                    isdirectory=is_dir,
                    config=fp.config,
                    version=fp.version,
                )
                # If our filepath is a distance descendant of the prefix then start filling in
                # the intermediate paths
            elseif prefix in _parents
                i = findfirst(==(prefix), _parents)
                _parents[i + 1]
                # Otherwise we've established that the current filepath isn't a descendant
                # of the prefix and we should exit
            else
                return nothing
            end

            # If we aren't dealing with the root and we're doing topdown iteration then
            # insert the child into the results channel
            !isempty(k) && topdown && put!(chnl, child)

            # Apply our recursive call for the children as necessary
            # NOTE: We're relying on the `isdirectory` field rather than calling `isdir`
            # which will call out to AWS as a fallback.
            if child.isdirectory
                _walkpath!(
                    root, child, objects, chnl; topdown=topdown, onerror=onerror, kwargs...
                )
            end

            # If we aren't dealing with the root and we're doing bottom up iteration then
            # insert the child ion the result channel here
            !isempty(k) && !topdown && put!(chnl, child)
        catch e
            isa(e, Base.IOError) ? onerror(e) : rethrow()
        end
    end
end

"""
    stat(fp::S3Path)

Return the status struct for the S3 path analogously to `stat` for local directories.

Note that this cannot be used on a directory.  This is because S3 is a pure key-value store and internally does
not have a concept of directories.  In some cases, a directory may actually be an empty file, in which case
you should use `s3_get_meta`.
"""
function Base.stat(fp::S3Path)
    # Currently AWSS3 would require a s3_get_acl call to fetch
    # ownership and permission settings
    m = Mode(; user=(READ + WRITE), group=(READ + WRITE), other=(READ + WRITE))
    u = FilePathsBase.User()
    g = FilePathsBase.Group()
    blksize = 4096
    blocks = 0
    s = 0
    last_modified = DateTime(0)

    if isfile(fp)
        resp = s3_get_meta(get_config(fp), fp.bucket, fp.key; version=fp.version)

        # Example: "Thu, 03 Jan 2019 21:09:17 GMT"
        last_modified = DateTime(
            resp["Last-Modified"][1:(end - 4)], dateformat"e, d u Y H:M:S"
        )
        s = parse(Int, resp["Content-Length"])
        blocks = ceil(Int, s / 4096)
    end

    return Status(0, 0, m, 0, u, g, 0, s, blksize, blocks, last_modified, last_modified)
end

"""
    diskusage(fp::S3Path)

Compute the *total* size of all contents of a directory.  Note that there is no direct functionality
for this in the S3 API so it may be slow.
"""
function FilePathsBase.diskusage(fp::S3Path)
    return if isfile(fp)
        stat(fp).size
    else
        s3_directory_stat(get_config(fp), fp.bucket, fp.key)[1]
    end
end

"""
    lastmodified(fp::S3Path)

Returns a `DateTime` corresponding to the latest time at which the object (or, in the case of a
directory, any contained object) was modified.
"""
function lastmodified(fp::S3Path)
    return if isfile(fp)
        stat(fp).mtime
    else
        s3_directory_stat(get_config(fp), fp.bucket, fp.key)[2]
    end
end

# Need API for accessing object ACL permissions for this to work
FilePathsBase.isexecutable(fp::S3Path) = false
Base.isreadable(fp::S3Path) = true
Base.iswritable(fp::S3Path) = true
Base.ismount(fp::S3Path) = false

"""
    mkdir(fp::S3Path; recursive=false, exist_ok=false)

Create an empty directory at the S3 path `fp`.  If `recursive`, this will create any previously non-existent
directories which would contain `fp`.  An error will be thrown if an object exists at `fp` unless `exist_ok`.

Note that empty directories in S3 are actually 0-byte objects with the naming convention of a directory.

This will *not* create a bucket.
"""
function Base.mkdir(fp::S3Path; recursive=false, exist_ok=false)
    fp.isdirectory || throw(ArgumentError("S3Path folders must end with '/': $fp"))

    if exists(fp)
        !exist_ok && error("$fp already exists.")
    else
        if hasparent(fp) && !exists(parent(fp))
            if recursive
                # don't try to create buckets this way, minio at least really doesn't like it
                isempty(parent(fp).segments) ||
                    mkdir(parent(fp); recursive=recursive, exist_ok=exist_ok)
            else
                error(
                    "The parent of $fp does not exist. " *
                    "Pass `recursive=true` to create it.",
                )
            end
        end

        write(fp, "")
    end

    return fp
end

function Base.rm(fp::S3Path; recursive=false, kwargs...)
    if isdir(fp)
        files = readpath(fp)

        if recursive
            for f in files
                rm(f; recursive=recursive, kwargs...)
            end
        elseif length(files) > 0
            error("S3 path $fp is not empty. Use `recursive=true` to delete.")
        end
    end

    @debug "delete: $fp"
    return s3_delete(get_config(fp), fp.bucket, fp.key; version=fp.version)
end

# We need to special case sync with S3Paths because of how directories
# are handled again.
# NOTE: This method signature only makes sense with FilePathsBase 0.6.2, but
# 1) It'd be odd for other packages to restrict FilePathsBase to a patch release
# 2) Seems cleaner to have it fallback and error rather than having
# slightly inconsistent handling of edge cases between the two versions.
function FilePathsBase.sync(
    f::Function, src::AbstractPath, dst::S3Path; delete=false, overwrite=true
)
    # Throw an error if the source path doesn't exist at all
    exists(src) || throw(ArgumentError("Unable to sync from non-existent $src"))

    # If the top level source is just a file then try to just sync that
    # without calling walkpath
    if isfile(src)
        # If the destination exists then we should make sure it is a file and check
        # if we should copy the source over.
        if exists(dst)
            isfile(dst) || throw(ArgumentError("Unable to sync file $src to non-file $dst"))
            if overwrite && f(src, dst)
                cp(src, dst; force=true)
            end
        else
            cp(src, dst)
        end
    elseif isdir(src)
        if exists(dst)
            isdir(dst) ||
                throw(ArgumentError("Unable to sync directory $src to non-directory $dst"))
            # Create an index of all of the source files
            src_paths = collect(walkpath(src))

            #! format: off
            # https://github.com/domluna/JuliaFormatter.jl/issues/458
            index = Dict(
                Tuple(setdiff(p.segments, src.segments)) => i
                for (i, p) in enumerate(src_paths)
            )
            #! format: on

            for dst_path in walkpath(dst)
                k = Tuple(setdiff(dst_path.segments, dst.segments))

                if haskey(index, k)
                    src_path = src_paths[pop!(index, k)]
                    if overwrite && f(src_path, dst_path)
                        cp(src_path, dst_path; force=true)
                    end
                elseif delete
                    rm(dst_path; recursive=true)
                end
            end

            # Finally, copy over files that don't exist at the destination
            # But we need to iterate through it in a way that respects the original
            # walkpath order otherwise we may end up trying to copy a file before its parents.
            index_pairs = collect(pairs(index))
            index_pairs = index_pairs[sortperm(index_pairs; by=last)]
            for (seg, i) in index_pairs
                new_dst = S3Path(
                    tuple(dst.segments..., seg...),
                    dst.root,
                    dst.drive,
                    isdir(src_paths[i]),
                    nothing,
                    dst.config,
                )

                cp(src_paths[i], new_dst; force=true)
            end
        else
            cp(src, dst)
        end
    else
        throw(ArgumentError("$src is neither a file or directory."))
    end
end

# for some reason, sometimes we get back a `Pair`
# other times a `AbstractDict`.
function _pair_or_dict_get(p::Pair, k)
    first(p) == k || return nothing
    return last(p)
end
_pair_or_dict_get(d::AbstractDict, k) = get(d, k, nothing)

function _retrieve_prefixes!(results, objects, prefix_key, chop_head)
    objects === nothing && return nothing

    rm_key = s -> chop(s; head=chop_head, tail=0)

    for p in objects
        prefix = _pair_or_dict_get(p, prefix_key)

        if prefix !== nothing
            push!(results, rm_key(prefix))
        end
    end

    return nothing
end

function _readdir_add_results!(results, response, key_length)
    sizehint!(results, length(results) + parse(Int, response["KeyCount"]))

    common_prefixes = get(response, "CommonPrefixes", nothing)
    _retrieve_prefixes!(results, common_prefixes, "Prefix", key_length)

    contents = get(response, "Contents", nothing)
    _retrieve_prefixes!(results, contents, "Key", key_length)

    return get(response, "NextContinuationToken", nothing)
end

function Base.readdir(fp::S3Path; join=false, sort=true)
    if isdir(fp)
        k = fp.key
        key_length = length(k)
        results = String[]
        token = ""
        while token !== nothing
            response = @repeat 4 try
                params = Dict("delimiter" => "/", "prefix" => k)

                if !isempty(token)
                    params["continuation-token"] = token
                end
                parse(S3.list_objects_v2(fp.bucket, params; aws_config=get_config(fp)))
            catch e
                #! format: off
                @delay_retry if ecode(e) in ["NoSuchBucket"] end
                #! format: on
            end
            token = _readdir_add_results!(results, response, key_length)
        end

        # Filter out any empty object names which are valid in S3
        filter!(!isempty, results)

        # Sort results if sort=true
        sort && sort!(results)

        # Return results, possibly joined with the root path if join=true
        return join ? joinpath.(fp, results) : results
    else
        throw(ArgumentError("\"$fp\" is not a directory"))
    end
end

"""
    read(fp::S3Path; byte_range=nothing)

Fetch data from the S3 path as a `Vector{UInt8}`.  A subset of the object can be specified with
`byte_range` which should be a contiguous integer range, e.g. `1:4`.
"""
function Base.read(fp::S3Path; byte_range=nothing)
    return Vector{UInt8}(
        s3_get(
            get_config(fp),
            fp.bucket,
            fp.key;
            raw=true,
            byte_range=byte_range,
            version=fp.version,
        ),
    )
end

function Base.write(fp::S3Path, content::String; kwargs...)
    return Base.write(fp, Vector{UInt8}(content); kwargs...)
end

function Base.write(
    fp::S3Path,
    content::Vector{UInt8};
    part_size_mb=50,
    multipart::Bool=true,
    other_kwargs...,
)
    # avoid HTTPClientError('An HTTP Client raised an unhandled exception: string longer than 2147483647 bytes')
    MAX_HTTP_BYTES = 2147483647
    if fp.version !== nothing
        throw(ArgumentError("Can't write to a specific object version ($(fp.version))"))
    end

    if !multipart || length(content) < MAX_HTTP_BYTES
        return s3_put(get_config(fp), fp.bucket, fp.key, content)
    else
        io = IOBuffer(content)
        return s3_multipart_upload(
            get_config(fp), fp.bucket, fp.key, io, part_size_mb; other_kwargs...
        )
    end
end

function FilePathsBase.mktmpdir(parent::S3Path)
    fp = parent / string(uuid4(), "/")
    return mkdir(fp)
end


text:
# from https://github.com/queryverse/NamedTupleUtilities.jl

module NamedTupleUtilities
export remove, rename
"""
    remove(a::NamedTuple, v::Val{n})

Remove a field `n` from the `a` if it is in `a`.

```jldoctest
julia> remove((a=1,b=2,c=3),Val(:c))
(a = 1, b = 2)
```
"""
@generated function remove(a::NamedTuple{an}, ::Val{bn}) where {an, bn}
    names = ((i for i in an if i != bn)...,)
    types = Tuple{(fieldtype(a ,n) for n in names)...}
    vals = Expr[:(getfield(a, $(QuoteNode(n)))) for n in names]
    return :(NamedTuple{$names,$types}(($(vals...),)))
end

"""
    range(a::NamedTuple, b::Val{n}, c::Val{n})

Return a NamedTuple which retains the fields from `b` to `c` in `a`.
If `b` is not in `a`, then it will return empty NamedTuple.
If `c` is not in `a`, then it will return everything started with `b`.

```jldoctest
julia> range((a=1,b=2,c=3),Val(:a),Val(:b))
(a = 1, b = 2)
```
"""
@generated function Base.range(a::NamedTuple{an}, ::Val{bn}, ::Val{cn}) where {an, bn, cn}
    rangeStarted = false
    names = Symbol[]
    for n in an
        if n == bn
            rangeStarted = true
        end
        if rangeStarted
            push!(names, n)
        end
        if n == cn
            rangeStarted = false
            break
        end
    end
    types = Tuple{(fieldtype(a, n) for n in names)...}
    vals = Expr[:(getfield(a, $(QuoteNode(n)))) for n in names]
    return :(NamedTuple{$(names...,),$types}(($(vals...),)))
end

"""
    rename(a::NamedTuple, b::Val{n}, c::Val{n})

Return a NamedTuple derived from `a` in which the the field from `b` is renamed to `c`.
If `b` is not in `a`, then it will return the original NamedTuple.
If `c` is in `a`, then `ERROR: duplicate field name in NamedTuple: "c" is not unique` will occur.

```jldoctest
julia> rename((a = 1, b = 2, c = 3),Val(:a),Val(:d))
(d = 1, b = 2, c = 3)

julia> rename((a = 1, b = 2, c = 3),Val(:m),Val(:d))
(a = 1, b = 2, c = 3)

julia> rename((a = 1, b = 2, c = 3),Val(:a),Val(:c))
ERROR: duplicate field name in NamedTuple: "c" is not unique
```
"""
@generated function rename(a::NamedTuple{an}, ::Val{bn}, ::Val{cn}) where {an, bn, cn}
    names = Symbol[]
    typesArray = DataType[]
    vals = Expr[]
    for n in an
        if n == bn
            push!(names, cn)
        else
            push!(names, n)
        end
        push!(typesArray, fieldtype(a, n))
        push!(vals, :(getfield(a, $(QuoteNode(n)))))
    end
    types = Tuple{typesArray...}
    return :(NamedTuple{$(names...,),$types}(($(vals...),)))
end

"""
    startswith(a::NamedTuple, b::Val{n})

Return a NamedTuple which retains the fields with names started with `b` in `a`.

```jldoctest
julia> startswith((abc=1,bcd=2,cde=3),Val(:a))
(abc = 1,)
```
"""
@generated function Base.startswith(a::NamedTuple{an}, ::Val{bn}) where {an, bn}
    names = ((i for i in an if startswith(String(i), String(bn)))...,)
    types = Tuple{(fieldtype(a ,n) for n in names)...}
    vals = Expr[:(getfield(a, $(QuoteNode(n)))) for n in names]
    return :(NamedTuple{$names,$types}(($(vals...),)))
end

"""
    endswith(a::NamedTuple, b::Val{n})

Return a NamedTuple which retains the fields with names ended with `b` in `a`.

```jldoctest
julia> endswith((abc=1,bcd=2,cde=3),Val(:d))
(bcd = 2,)
```
"""
@generated function Base.endswith(a::NamedTuple{an}, ::Val{bn}) where {an, bn}
    names = ((i for i in an if endswith(String(i), String(bn)))...,)
    types = Tuple{(fieldtype(a ,n) for n in names)...}
    vals = Expr[:(getfield(a, $(QuoteNode(n)))) for n in names]
    return :(NamedTuple{$names,$types}(($(vals...),)))
end

"""
    occursin(a::NamedTuple, b::Val{n})

Return a NamedTuple which retains the fields with names contains `b` as a substring.

```jldoctest
julia> occursin((abc=1,bcd=2,cde=3),Val(:d))
(bcd = 2, cde = 3)
```
"""
@generated function Base.occursin(a::NamedTuple{an}, ::Val{bn}) where {an, bn}
    names = ((i for i in an if occursin(String(bn), String(i)))...,)
    types = Tuple{(fieldtype(a ,n) for n in names)...}
    vals = Expr[:(getfield(a, $(QuoteNode(n)))) for n in names]
    return :(NamedTuple{$names,$types}(($(vals...),)))
end

end # module


text:
using Polylogarithms
using SpecialFunctions
using Test
using DataFrames, CSV
import Base.MathConstants: π, pi, ℯ, e, γ, eulergamma, catalan, φ, golden
include("test_defs.jl")

@testset "Harmonic numbers" begin

    @testset "    throws errors" begin
        @test_throws DomainError stieltjes(-1)
        @test_throws DomainError stieltjes(11)
        @test_throws MethodError stieltjes(1.5)
    end

    @testset "    types" begin
        @test typeof(stieltjes(1)) == Float64
    end

    @testset "    values" begin
        @test stieltjes(0) ≅ γ
        function gen_euler_calc(n::Integer)
            m = 10000000 # just choose a big value for the moment
            total = -log(m)^(n+1) / (n+1)
            for k=1:m
                total += log(k)^n / k
            end
            return total
        end
        @test abs(stieltjes(1) - gen_euler_calc(1)) < 1.0e-5 # this series is a little crude as a test
        @test abs(stieltjes(2) - gen_euler_calc(2)) < 1.0e-4 # this series is a little crude as a test
        @test abs(stieltjes(3) - gen_euler_calc(3)) < 1.0e-3 # this series is a little crude as a test
    end

end


text:
using Documenter, UnixMmap
include(joinpath(@__DIR__, "sys_override.jl"))

doctest = "--fix"  in ARGS ? :fix :
          "--test" in ARGS ? true : false

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

makedocs(
    sitename = "UnixMmap.jl",
    authors  = "Justin Willmert",
    doctest  = doctest,
    doctestfilters = Regex[
        r"Ptr{0x[0-9a-f]+}",
        r"[0-9\.]+ seconds \(.*\)"
    ],
    format   = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical  = "https://jmert.github.io/UnixMmap.jl",
        assets     = String[],
        mathengine = Documenter.MathJax3(),
    ),
    pages = [
        "Home" => "index.md",
        "API Reference" => [
            "Public" => "lib/public.md",
        ]
    ],
    repo = "https://github.com/jmert/UnixMmap.jl/blob/{commit}{path}#L{line}",
)

deploydocs(
    repo = "github.com/jmert/UnixMmap.jl",
    push_preview = true,
)


text:
"""
    hahn(x, p::Integer, α, β, N::Integer)

Evaluate the Hahn polynomial with parameters `α`, `β`, `N` of degree `p` at `x`
using the three term recursion [Öffner, Zweidimensionale klassische und diskrete
orthogonale Polynome, Chapter 5].
"""
function hahn(x, p::Integer, α, β, N)
    T = typeof( (α+β+2)*x/(N*(α+1)) )
    a = one(T)
    b = 1 - (α+β+2)*x/(N*(α+1))

    if p <= 0
        return a
    elseif p == 1
        return b
    elseif typeof(N) <: Integer && p > N
        throw(ArgumentError("p==$p must not be larger than N==$N."))
    end

    for n in 2:p
        An = T( (n+α+β)*(n+α)*(N-n+1) ) / ( (2n-1+α+β)*(2n+α+β) )
        Cn = T( (n-1)*(n+α+β+N)*(n-1+β) ) / ( (2n-2+α+β)*(2n-1+α+β) )
        a, b = b, ( (An+Cn-x)*b - Cn*a ) / An
    end

    b
end