text:
@testset "mean" begin

@testset "create_mean and apply_mean! 1D" begin
    N = SVec(10,)
    M = prod(N)

    Fu = rand(N...)
    Gv = similar(Fu)
    gv = zeros(M)

    Mask = similar(Fu)  # masking array

    K = length(N)
    KM = K * M
    F = rand(ComplexF64, N..., K)
    f = F[:]
    G = similar(F)
    g = zeros(ComplexF64, KM)

    for isfwd = (true,false), isbloch = (true,false)  # (backward,forward difference), (Bloch,symmetric)
        ∆l = rand.(tuple(N...))
        ∆l′⁻¹ = rand.(tuple(N...))
        Mdiag = spzeros(KM,KM)

        for nw = 1:K
            sub_on = Vector{Int}(undef, K)
            sub_off = Vector{Int}(undef, K)

            Nw = N[nw]
            ∆w = ∆l[nw]
            ∆w′⁻¹ = ∆l′⁻¹[nw]

            Mws = spzeros(M,M)

            for ind_on = 1:M
                sub_on .= CartesianIndices(N.data)[ind_on].I  # subscripts of diagonal entry
                indw_on = sub_on[nw]
                ∆w_on = ∆w[indw_on]
                ∆w′_on⁻¹ = ∆w′⁻¹[indw_on]
                Mws[ind_on,ind_on] = 0.5 * ∆w_on * ∆w′_on⁻¹  # diagonal entries

                # Calculate the column index of the off-diagonal entry in the row `ind`.
                sub_off .= sub_on  # subscripts of off-diagonal entry
                if isfwd  # forward difference
                    if sub_off[nw] == Nw
                        sub_off[nw] = 1
                    else
                        sub_off[nw] += 1
                    end
                else  # backward difference
                    if sub_off[nw] == 1
                        if isbloch
                            sub_off[nw] = Nw
                        else  # symmetry boundary
                            # Leave sub_off[nw] at 1
                        end
                    else
                        sub_off[nw] -= 1
                    end
                end

                indw_off = sub_off[nw]
                ∆w_off = ∆w[indw_off]
                ind_off = LinearIndices(N.data)[sub_off...]  # index of off-diagonal entry
                if !isfwd && !isbloch && sub_on[nw]==1  # bacward difference && symmetry boundary at negative end
                    Mws[ind_on, ind_off] = ∆w_off * ∆w′_on⁻¹  # double diagonal entry
                else
                    Mws[ind_on, ind_off] = 0.5 * ∆w_off * ∆w′_on⁻¹  # off-diagonal entry
                end
            end

            if isfwd && !isbloch  # forward difference and symmetry boundary
                # Initialize the masking array.
                Mask .= 1
                Mask[Base.setindex(axes(Mask), 1, nw)...] = 0  # not .= 0 unlike code for 2D and 3D

                # The input fields at the symmetry boundary should be zero, so apply the
                # mask to the input field.
                Mws = Mws * spdiagm(0=>Mask[:])
            end

            # Test create_m.
            @test create_m̂(nw, isfwd, N, ∆w, ∆w′⁻¹, isbloch) == Mws

            # Test apply_m̂!.
            fu = Fu[:]
            mul!(gv, Mws, fu)
            apply_m̂!(Gv, Fu, Val(:(=)), nw, isfwd, ∆w, ∆w′⁻¹, isbloch)
            @test Gv[:] ≈ gv

            # Construct Mdiag, Msup, Msub.
            Is = 1+(nw-1)*M:nw*M
            Mdiag[Is,Is] .= Mws
        end  # for nw
        isfwdK = fill(isfwd, K)
        isblochK = fill(isbloch, K)
        @test create_mean(isfwdK, ∆l, ∆l′⁻¹, isblochK, order_cmpfirst=false) == Mdiag

        # Test apply_mean!.
        mul!(g, Mdiag, f)
        apply_mean!(G, F, Val(:(=)), isfwdK, ∆l, ∆l′⁻¹, isblochK)
        @test G[:] ≈ g
    end  # isfwd = ..., isbloch = ...
end  # @testset "create_mean and apply_mean! 1D"

@testset "create_mean and apply_mean! 2D" begin
    N = SVec(8,9)
    M = prod(N)

    Fu = rand(N...)
    Gv = similar(Fu)
    gv = zeros(M)

    Mask = similar(Fu)  # masking array

    K = length(N)
    KM = K * M
    F = rand(ComplexF64, N..., K)
    f = F[:]
    G = similar(F)
    g = zeros(ComplexF64, KM)

    for isfwd = (true,false), isbloch = (true,false)  # (backward,forward difference), (Bloch,symmetric)
        ∆l = rand.(tuple(N...))
        ∆l′⁻¹ = rand.(tuple(N...))
        Mdiag = spzeros(KM,KM)

        for nw = 1:K
            sub_on = Vector{Int}(undef, K)
            sub_off = Vector{Int}(undef, K)

            Nw = N[nw]
            ∆w = ∆l[nw]
            ∆w′⁻¹ = ∆l′⁻¹[nw]

            Mws = spzeros(M,M)

            for ind_on = 1:M
                sub_on .= CartesianIndices(N.data)[ind_on].I  # subscripts of diagonal entry
                indw_on = sub_on[nw]
                ∆w_on = ∆w[indw_on]
                ∆w′_on⁻¹ = ∆w′⁻¹[indw_on]
                Mws[ind_on,ind_on] = 0.5 * ∆w_on * ∆w′_on⁻¹  # diagonal entries

                # Calculate the column index of the off-diagonal entry in the row `ind`.
                sub_off .= sub_on  # subscripts of off-diagonal entry
                if isfwd  # forward difference
                    if sub_off[nw] == Nw
                        sub_off[nw] = 1
                    else
                        sub_off[nw] += 1
                    end
                else  # backward difference
                    if sub_off[nw] == 1
                        if isbloch
                            sub_off[nw] = Nw
                        else  # symmetry boundary
                            # Leave sub_off[nw] at 1
                        end
                    else
                        sub_off[nw] -= 1
                    end
                end

                indw_off = sub_off[nw]
                ∆w_off = ∆w[indw_off]
                ind_off = LinearIndices(N.data)[sub_off...]  # index of off-diagonal entry
                if !isfwd && !isbloch && sub_on[nw]==1  # bacward difference && symmetry boundary at negative end
                    Mws[ind_on, ind_off] = ∆w_off * ∆w′_on⁻¹  # double diagonal entry
                else
                    Mws[ind_on, ind_off] = 0.5 * ∆w_off * ∆w′_on⁻¹  # off-diagonal entry
                end
            end

            if isfwd && !isbloch  # forward difference and symmetry boundary
                # Initialize the masking array.
                Mask .= 1
                Mask[Base.setindex(axes(Mask), 1, nw)...] .= 0

                # The input fields at the symmetry boundary should be zero, so apply the
                # mask to the input field.
                Mws = Mws * spdiagm(0=>Mask[:])
            end

            # Test create_m.
            @test create_m̂(nw, isfwd, N, ∆w, ∆w′⁻¹, isbloch) == Mws

            # Test apply_m̂!.
            fu = Fu[:]
            mul!(gv, Mws, fu)
            apply_m̂!(Gv, Fu, Val(:(=)), nw, isfwd, ∆w, ∆w′⁻¹, isbloch)
            @test Gv[:] ≈ gv

            # Construct Mdiag, Msup, Msub.
            Is = 1+(nw-1)*M:nw*M
            Mdiag[Is,Is] .= Mws
        end  # for nw
        isfwdK = fill(isfwd, K)
        isblochK = fill(isbloch, K)
        @test create_mean(isfwdK, ∆l, ∆l′⁻¹, isblochK, order_cmpfirst=false) == Mdiag

        # Test apply_mean!.
        mul!(g, Mdiag, f)
        apply_mean!(G, F, Val(:(=)), isfwdK, ∆l, ∆l′⁻¹, isblochK)
        @test G[:] ≈ g
    end  # isfwd = ..., isbloch = ...
end  # @testset "create_mean and apply_mean! 2D"

@testset "create_mean and apply_mean! 3D" begin
    N = SVec(8,9,10)
    M = prod(N)

    Fu = rand(N...)
    Gv = similar(Fu)
    gv = zeros(M)

    Mask = similar(Fu)  # masking array

    K = length(N)
    KM = K * M
    F = rand(ComplexF64, N..., K)
    f = F[:]
    G = similar(F)
    g = zeros(ComplexF64, KM)

    for isfwd = (true,false), isbloch = (true,false)  # (backward,forward difference), (Bloch,symmetric)
        ∆l = rand.(tuple(N...))
        ∆l′⁻¹ = rand.(tuple(N...))
        Mdiag = spzeros(KM,KM)

        for nw = 1:K
            sub_on = Vector{Int}(undef, K)
            sub_off = Vector{Int}(undef, K)

            Nw = N[nw]
            ∆w = ∆l[nw]
            ∆w′⁻¹ = ∆l′⁻¹[nw]

            Mws = spzeros(M,M)

            for ind_on = 1:M
                sub_on .= CartesianIndices(N.data)[ind_on].I  # subscripts of diagonal entry
                indw_on = sub_on[nw]
                ∆w_on = ∆w[indw_on]
                ∆w′_on⁻¹ = ∆w′⁻¹[indw_on]
                Mws[ind_on,ind_on] = 0.5 * ∆w_on * ∆w′_on⁻¹  # diagonal entries

                # Calculate the column index of the off-diagonal entry in the row `ind`.
                sub_off .= sub_on  # subscripts of off-diagonal entry
                if isfwd  # forward difference
                    if sub_off[nw] == Nw
                        sub_off[nw] = 1
                    else
                        sub_off[nw] += 1
                    end
                else  # backward difference
                    if sub_off[nw] == 1
                        if isbloch
                            sub_off[nw] = Nw
                        else  # symmetry boundary
                            # Leave sub_off[nw] at 1
                        end
                    else
                        sub_off[nw] -= 1
                    end
                end

                indw_off = sub_off[nw]
                ∆w_off = ∆w[indw_off]
                ind_off = LinearIndices(N.data)[sub_off...]  # index of off-diagonal entry
                if !isfwd && !isbloch && sub_on[nw]==1  # bacward difference && symmetry boundary at negative end
                    Mws[ind_on, ind_off] = ∆w_off * ∆w′_on⁻¹  # double diagonal entry
                else
                    Mws[ind_on, ind_off] = 0.5 * ∆w_off * ∆w′_on⁻¹  # off-diagonal entry
                end
            end

            if isfwd && !isbloch  # forward difference and symmetry boundary
                # Initialize the masking array.
                Mask .= 1
                Mask[Base.setindex(axes(Mask), 1, nw)...] .= 0

                # The input fields at the symmetry boundary should be zero, so apply the
                # mask to the input field.
                Mws = Mws * spdiagm(0=>Mask[:])
            end

            # Test create_m.
            @test create_m̂(nw, isfwd, N, ∆w, ∆w′⁻¹, isbloch) == Mws

            # Test apply_m̂!.
            fu = Fu[:]
            mul!(gv, Mws, fu)
            apply_m̂!(Gv, Fu, Val(:(=)), nw, isfwd, ∆w, ∆w′⁻¹, isbloch)
            @test Gv[:] ≈ gv

            # Construct Mdiag, Msup, Msub.
            Is = 1+(nw-1)*M:nw*M
            Mdiag[Is,Is] .= Mws
        end  # for nw
        isfwdK = fill(isfwd, K)
        isblochK = fill(isbloch, K)
        @test create_mean(isfwdK, ∆l, ∆l′⁻¹, isblochK, order_cmpfirst=false) == Mdiag

        # Test apply_mean!.
        mul!(g, Mdiag, f)
        apply_mean!(G, F, Val(:(=)), isfwdK, ∆l, ∆l′⁻¹, isblochK)
        @test G[:] ≈ g
    end  # isfwd = ..., isbloch = ...
end  # @testset "create_mean and apply_mean! 3D"

end  # @testset "mean"


text:
# wordcount.jl
#
# Implementation of parallelized "word-count" of a text, inspired by the
# Hadoop WordCount example. Uses @spawn and fetch() to parallelize
# the "map" task. Reduce is currently done single-threaded.
#
# To run in parallel on a string stored in variable `text`:
#  julia -p <N>
#  julia> require("<julia_dir>/examples/wordcount.jl")
#  julia> ...(define text)...
#  julia> counts=parallel_wordcount(text)
#
# Or to run on a group of files, writing results to an output file:
#  julia -p <N>
#  julia> require("<julia_dir>/examples/wordcount.jl")
#  julia> wordcount_files("/tmp/output.txt", "/tmp/input1.txt","/tmp/input2.txt",...)

# "Map" function.
# Takes a string. Returns a Dict with the number of times each word
# appears in that string.
function wordcount(text)
    words=split(text,[' ','\n','\t','-','.',',',':',';'],false)
    counts=Dict()
    for w = words
        counts[w]=get(counts,w,0)+1
    end
    return counts
end

# "Reduce" function.
# Takes a collection of Dicts in the format returned by wordcount()
# Returns a Dict in which words that appear in multiple inputs
# have their totals added together.
function wcreduce(wcs)
    counts=Dict()
    for c in wcs, (k,v) in c
        counts[k] = get(counts,k,0)+v
    end
    return counts
end

# Splits input string into nprocs() equal-sized chunks (last one rounds up),
# and @spawns wordcount() for each chunk to run in parallel. Then fetch()s
# results and performs wcreduce().
function parallel_wordcount(text)
    lines=split(text,'\n',false)
    np=nprocs()
    unitsize=ceil(length(lines)/np)
    wcounts={}
    rrefs={}
    # spawn procs
    for i=1:np
        first=unitsize*(i-1)+1
        last=unitsize*i
        if last>length(lines)
            last=length(lines)
        end
        subtext=join(lines[int(first):int(last)],"\n")
        push!(rrefs, @spawn wordcount( subtext ) )
    end
    # fetch results
    while length(rrefs)>0
        push!(wcounts,fetch(pop!(rrefs)))
    end
    # reduce
    count=wcreduce(wcounts)
    return count
end

# Takes the name of a result file, and a list of input file names.
# Combines the contents of all files, then performs a parallel_wordcount
# on the resulting string. Writes the results to result_file.
function wordcount_files(result_file,inputs...)
    text = ""
    for file in inputs
        open(file) do f
            text *= readall(f)
        end
    end
    wc = parallel_wordcount(text)
    open(result_file,"w") do f
        for (k,v) in wc
            println(f, k,"=",v)
        end
    end
end

ext: jl

max_stars_repo_name: Aerlinger/julia

max_issues_repo_name: Aerlinger/julia

max_forks_repo_name: Aerlinger/julia


text:
function f(x)
a=x*5
b=a+3
c=a-4
d=b/c
end

max_stars_repo_path: Module 2/Chapter04/F(x).jl

max_issues_repo_path: Module 2/Chapter04/F(x).jl

max_forks_repo_path: Module 2/Chapter04/F(x).jl


text:
using Manifolds, Manopt, Test, Dates

@testset "proximal maps" begin
    #
    # prox_TV
    p = [1.0, 0.0, 0.0]
    q = [0.0, 1.0, 0.0]
    M = Sphere(2)
    N = PowerManifold(M, NestedPowerRepresentation(), 2)
    @test_throws ErrorException prox_distance(M, 1.0, p, q, 3)
    @test_throws ErrorException prox_distance!(M, p, 1.0, p, q, 3)
    @test distance(
        M, prox_distance(M, distance(M, p, q) / 2, p, q, 1), shortest_geodesic(M, p, q, 0.5)
    ) ≈ 0
    t = similar(p)
    prox_distance!(M, t, distance(M, p, q) / 2, p, q, 1)
    @test t == prox_distance(M, distance(M, p, q) / 2, p, q, 1)
    (r, s) = prox_TV(M, π / 4, (p, q))
    X = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
    prox_TV!(M, X, π / 4, (p, q))
    @test norm(r - s) < eps(Float64)
    @test norm(X[1] - s) < eps(Float64)
    @test norm(X[2] - r) < eps(Float64)
    # i.e. they are moved together
    @test distance(M, r, s) < eps(Float64)
    (t, u) = prox_TV(M, π / 8, (p, q))
    @test_throws ErrorException prox_TV(M, π, (p, q), 3)
    @test_throws ErrorException prox_TV!(M, [p, q], π, (p, q), 3)
    # they cross correlate
    @test (
        abs(t[1] - u[2]) < eps(Float64) &&
        abs(t[2] - u[1]) < eps(Float64) &&
        abs(t[3] - u[3]) < eps(Float64)
    )
    @test distance(M, t, u) == π / 4 # and have moved half their distance
    #
    (v, w) = prox_TV(M, 1.0, (p, q), 2)
    vC, wC = shortest_geodesic(M, p, q, [1 / 3, 2 / 3])
    @test distance(M, v, vC) ≈ 0
    @test distance(M, w, wC) ≈ 0
    P = [similar(p), similar(q)]
    prox_TV!(M, P, 1.0, (p, q), 2)
    @test P == [v, w]
    # prox_TV on Power
    T = prox_TV(N, π / 8, [p, q])
    @test distance(N, T, [t, u]) ≈ 0
    # parallelprox_TV
    N2 = PowerManifold(M, NestedPowerRepresentation(), 3)
    r = geodesic(M, p, q, 0.5)
    s, t = prox_TV(M, π / 16, (r, q))
    u, v = prox_TV(M, π / 16, (p, r))
    y = prox_parallel_TV(N2, π / 16, [[p, r, q], [p, r, q]])
    yM = [similar.(ye) for ye in y]
    prox_parallel_TV!(N2, yM, π / 16, [[p, r, q], [p, r, q]])
    @test y == yM
    @test distance(N2, y[1], [p, s, t]) ≈ 0 # even indices in first comp
    @test distance(N2, y[2], [u, v, q]) ≈ 0 # odd in second
    # dimensions of x have to fit, here they don't
    @test_throws ErrorException prox_parallel_TV(N2, π / 16, [[p, r, q]])
    @test_throws ErrorException prox_parallel_TV!(N2, yM, π / 16, [[p, r, q]])
    # prox_TV2
    p2, r2, q2 = prox_TV2(M, 1.0, (p, r, q))
    y = [similar(p) for _ in 1:3]
    prox_TV2!(M, y, 1.0, (p, r, q))
    @test y ≈ [p2, r2, q2]
    sum(distance.(Ref(M), [p, r, q], [p2, r2, q2])) ≈ 0
    @test_throws ErrorException prox_TV2(M, 1.0, (p, r, q), 2) # since prox_TV is only defined for p=1
    distance(
        PowerManifold(M, NestedPowerRepresentation(), 3),
        [p2, r2, q2],
        prox_TV2(PowerManifold(M, NestedPowerRepresentation(), 3), 1.0, [p, r, q]),
    ) ≈ 0
    # Circle
    M2 = Circle()
    N2 = PowerManifold(M2, 3)
    pS, rS, qS = [-0.5, 0.1, 0.5]
    d = sum([pS, rS, qS] .* [1.0, -2.0, 1.0])
    m = min(0.3, abs(Manopt.sym_rem(d) / 6))
    s = sign(Manopt.sym_rem(d))
    pSc, rSc, qSc = Manopt.sym_rem.([pS, rS, qS] .- m .* s .* [1.0, -2.0, 1.0])
    pSr, rSr, qSr = prox_TV2(M2, 0.3, (pS, rS, qS))
    @test sum(distance.(Ref(M2), [pSc, rSc, qSc], [pSr, rSr, qSr])) ≈ 0
    # p=2
    t = 0.3 * Manopt.sym_rem(d) / (1 + 0.3 * 6.0)
    @test sum(
        distance.(
            Ref(M2),
            [prox_TV2(M2, 0.3, (pS, rS, qS), 2)...],
            [pS, rS, qS] .- t .* [1.0, -2.0, 1.0],
        ),
    ) ≈ 0
    @test prox_TV2(N2, 0.3, [pS, rS, qS]) == [pSr, rSr, qSr]
    # others fail
    @test_throws ErrorException prox_TV2(M2, 0.3, (pS, rS, qS), 3)
    # Rn
    M3 = Euclidean(1)
    pR, rR, qR = [pS, rS, qS]
    m = min.(Ref(0.3), abs.([pR, rR, qR] .* [1.0, -2.0, 1.0]) / 6)
    s = sign(d)  # we can reuse d
    pRc, rRc, qRc = [pR, rR, qR] .- m .* s .* [1.0, -2.0, 1.0]
    pRr, rRr, qRr = prox_TV2(M3, 0.3, (pR, rR, qR))
    @test sum(distance.(Ref(M3), [pRc, rRc, qRc], [pRr, rRr, qRr])) ≈ 0
    # p=2
    t = 0.3 * d / (1 + 0.3 * 6.0)
    @test sum(
        distance.(
            Ref(M3),
            [prox_TV2(M3, 0.3, (pR, rR, qR), 2)...],
            [pR, rR, qR] .- t .* [1.0, -2.0, 1.0],
        ),
    ) ≈ 0
    # others fail
    @test_throws ErrorException prox_TV2(M3, 0.3, (pR, rR, qR), 3)
    #
    # collaborative integer tests
    #
    @test_throws ErrorException prox_TV2(M3, 0.3, (pS, rS, qS), 3)
    ξR, ηR, νR = [pS, rS, qS]
    N3 = PowerManifold(M3, 3)
    P = [pR rR qR]
    Ξ = [ξR ηR νR]
    Θ = similar(Ξ)
    @test project_collaborative_TV(N3, 0.0, P, Ξ, 1, 1) == Ξ
    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 1, 1)
    @test Θ == Ξ
    @test project_collaborative_TV(N3, 0.0, P, Ξ, 1.0, 1) == Ξ
    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 1.0, 1)
    @test Θ == Ξ
    @test project_collaborative_TV(N3, 0.0, P, Ξ, 1, 1.0) == Ξ
    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 1, 1.0)
    @test Θ == Ξ
    @test project_collaborative_TV(N3, 0.0, P, Ξ, 1.0, 1.0) == Ξ
    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 1.0, 1.0)
    @test Θ == Ξ

    @test project_collaborative_TV(N3, 0.0, P, Ξ, 2, 1) == Ξ
    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 2, 1)
    @test Θ == Ξ
    @test norm(N3, P, project_collaborative_TV(N3, 0.0, P, Ξ, 2, Inf)) ≈ norm(Ξ)
    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 2, Inf)
    @test norm(N3, P, Θ) ≈ norm(Ξ)
    @test sum(abs.(project_collaborative_TV(N3, 0.0, P, Ξ, 1, Inf))) ≈ 1.0
    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 1, Inf)
    @test sum(abs.(Θ)) ≈ 1.0
    @test norm(N3, P, project_collaborative_TV(N3, 0.0, P, Ξ, Inf, Inf)) ≈ norm(Ξ)
    project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, Inf, Inf)
    @test norm(N3, P, Θ) ≈ norm(Ξ)
    @test_throws ErrorException project_collaborative_TV(N3, 0.0, P, Ξ, 3, 3)
    @test_throws ErrorException project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 3, 3)
    @test_throws ErrorException project_collaborative_TV(N3, 0.0, P, Ξ, 3, 1)
    @test_throws ErrorException project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 3, 1)
    @test_throws ErrorException project_collaborative_TV(N3, 0.0, P, Ξ, 3, Inf)
    @test_throws ErrorException project_collaborative_TV!(N3, Θ, 0.0, P, Ξ, 3, Inf)

    @testset "Multivariate project collaborative TV" begin
        S = Sphere(2)
        M = PowerManifold(S, NestedPowerRepresentation(), 2, 2, 2)
        p = [zeros(3) for i in [1, 2], j in [1, 2], k in [1, 2]]
        p[1, 1, 1] = [1.0, 0.0, 0.0]
        p[1, 2, 1] = 1 / sqrt(2) .* [1.0, 1.0, 0.0]
        p[2, 1, 1] = 1 / sqrt(2) .* [1.0, 0.0, 1.0]
        p[2, 2, 1] = [0.0, 1.0, 0.0]
        p[:, :, 2] = deepcopy(p[:, :, 1])
        X = zero_vector(M, p)
        X[1, 1, 1] .= [0.0, 0.5, 0.5]
        Y = zero_vector(M, p)
        @test norm(project_collaborative_TV(M, 1, p, X, 2, 1)) ≈ 0
        project_collaborative_TV!(M, Y, 1, p, X, 2, 1)
        @test norm(Y) ≈ 0
        @test norm(project_collaborative_TV(M, 0.5, p, X, 2, 1)) ≈ (norm(X[1, 1, 1]) - 0.5)
        project_collaborative_TV!(M, Y, 0.5, p, X, 2, 1)
        @test norm(Y) ≈ (norm(X[1, 1, 1]) - 0.5)
        Nf = PowerManifold(S, NestedPowerRepresentation(), 2, 2, 1)
        @test_throws ErrorException project_collaborative_TV(Nf, 1, p, X, 2, 1)
        @test_throws ErrorException project_collaborative_TV!(Nf, Y, 1, p, X, 2, 1)
    end
end


text:
import FourierTools

# Keeps track of the number of spectra in this session.

let spectrumIndex = 0
    global spectrumCounter() = (spectrumIndex += 1)
end

"""
    Spectrum{T<:Real} <: AbstractVector{T}

    Spectrum(energy::EnergyScale, data::Vector{<:Real}, props::Dict{Symbol,Any})

Construct a structure to hold spectrum data (energy scale, counts and metadata).

See [`NeXLSpectrum.EnergyScale`](@ref) or [`NeXLSpectrum.LinearEnergyScale`](@ref)

Example:

    julia> spec = Spectrum(LinearEnergyScale(0.0,10.0),
                     collect(1:1024),
                     Dict{Symbol,Any}(:BeamEnergy=>10.0e3, :LiveTime=>30.0))
                                                                                ** 1024.0
                                                                           *********
                                                                   *****************
                                                            ************************
                                                     *******************************
                                              **************************************
                                       *********************************************
                                ****************************************************
                         ***********************************************************
                  ******************************************************************
           *************************************************************************
    ******************************************************************************** 10.0 keV
    Spectrum[3062][10.0 keV, Unknown, 525000.0 counts]

Spectra are usually loaded from disk using:

    s1 = loadspectrum(joinpath(path,"ADM-6005a_1.msa")) # Auto-detects the file type (supports ISO/EMSA, Bruker, ASPEX, ??? files)
      *                                                                                128860.0
      *
      *
      *
      *
      *
      *     *
      *     *
      *     *
      * ** **
      * ** *****        **  **
    **************************************************************************************************** 20.0 keV
    Spectrum{Float64}[ADM-6005a_1, -484.20818 + 5.01716⋅ch eV, 4096 ch, 20.0 keV, Unknown, 6.81e6 counts]

`Spectrum` implements indexing using various different mechanisms.  If `spec` is a `Spectrum` then

    spec[123] # will return the number of counts in channel 123
	spec[123:222] # will return a Vector of counts from channel 123:222
    spec[134.] # will return the number of counts in the channel at energy 134.0 eV
    spec[134.0:270.0] # will return a Vector of counts for channels with energies between 134.0 eV and 270.0 eV
    spec[n"Ca L3-M5"] # Counts in the channel at the Ca L3-M5 characteristic X-ray energy
    spec[:Comment] # will return the meta-data item named :Comment

Basic spectrum math is supported using the operators *, /, ÷, +, and -.  If `s1`, `s2` and `s3` are `Spectrum` objects
then:

    2*s1 # A Spectrum containing twice as many counts in each channel
    2.0*s1 # A Spectrum containing twice as many counts in each channel
    s1/2.0 # A spectrum containing half as many counts in each channel
    s1÷2 # A spectrum containing half as many counts in each channel (Only works for Spectrum{<:Integer})
    s1 + s2 # A Spectrum containing channel-by-channel sum of counts in s1 and s1 and common properties.
    s1 - s2 # The channel-by-channel difference

Spectrum metadata is identified by a `Symbol`. These `Symbol`s are used within NeXLSpectrum. You can create others
to associate other data items with a `Spectrum`.

    :BeamEnergy    # In eV
	:Elevation     # In radians
	:TakeOffAngle  # In radians (Detector position)
    :Azimuthal     # In radians (Detector position)
	:WorkingDistance # In cm (not mm!!!!)
    :LiveTime      # In seconds
    :RealTime      # In seconds
    :DeadFraction  # Fractional
    :ProbeCurrent  # In nano-amps
    :Name          # A `String`
    :Owner         # A `String`
    :Sample        # Description of the sample
    :StagePosition # A Dict{Symbol,Real} with entries :X, :Y, :Z, :R, :T, B: in cm and degrees
    :Comment       # A `String`
    :Composition   # A `Material` (known composition, not measured)
    :Elements      # A collection of elements in the material like `[ n"Fe", n"Ca", n"Si" ]`
    :Detector      # A Detector like a BasicEDS or another EDSDetector
    :Filename      # Source filename
    :Coating       # A `Film` or `Film[]` (eg. 10 nm of C|Au etc.)
	:AcquisitionTime # Date and time of acquisition (`DateTime` struct)
	:Signature     # Dict{Element,Real} with the "particle signature"

Spectrum Image items:

    :ImageRotation # Rotation of the primary scan direction from the :X axis towards the :Y axis in radians

Less common items:

	:ImageMag	   # Magnification (assuming a 3.5\" image) of the first image
	:ImageZoom     # Additional zoom for second image in a two image TIFF
	:Operator      # Analyst in ASPEX TIFF files
    :Image1, :Image2 ... # Images associated with the spectrum
    :BrukerThroughtput # Nominal throughtput setting on a Bruker detector
    :DetectorSerialNumber # EDS detector serial number
    :DetectorModel # Vendor model name
    :DetectorThickness # Thickness of detector active area
    :DeadLayerThickness # Thickness of Si dead layer on the entrance surface of the detector
    :Window        # Window construction details
    :DetectorSolidAngle # Collection solid angle of the X-ray detector
    :ChamberPressure # Vacuum presure in the sample chamber
    :ChamberAtmosphere # Nominally the composition of the residual gas in the chamber

XRF related items:

    :BeamEnergy  # Accelarating voltage within X-ray tube (eV)
    :XRFTubeAnode    # Element from which the X-ray tube is constructed
    :ProbeCurrent  # Electron current in the X-ray tube
    :XRFTubeIncidentAngle # Incident angle of electron beam in tube
    :XRFTubeTakeOffAngle # Take-off angle from tube
    :XRFExcitationAngle # Angle of incidence of the X-ray beam on the sample
    :XRFDetectionAngle # Angle of the detector relative to the sample
    :XRFExcitationPathLength # Distance from X-ray source to sample
    :XRFDetectionPathLength # Distance from the sample to the X-ray detector
    :XRFSampleTilt    #  Additional tilt of the sample
    :XRFTubeWindow   # Construction of the tube window

Not all spectra will define all properties.  Algorithms can define the `NeXLCore.minproperties(ty::Type)` method to specify
which properties are required by an algorithm of `ty::Type`.  Then `hasminrequired` and `requiredbutmissing` methods
will determine whether a `Spectrum` or `Dict{Symbol,Any}` is suitable for an algorithm.
"""
struct Spectrum{T<:Real} <: AbstractVector{T}
    energy::EnergyScale
    counts::Vector{T}
    properties::Dict{Symbol,Any}
    hash::UInt  # Stays fixed as underlying data changes
    function Spectrum(energy::EnergyScale, data::AbstractVector{<:Real}, props::Dict{Symbol,<:Any})
        res = new{typeof(data[1])}(
            energy,
            data,
            convert(Dict{Symbol, Any}, props),
            hash(energy, hash(data, hash(props))),
        )
        res.properties[:Name] = get(res.properties, :Name, "Spectrum[$(spectrumCounter())]")
        res
    end
end

Base.hash(spec::Spectrum, h::UInt) = hash(spec.hash, h)
Base.isequal(spec1::Spectrum, spec2::Spectrum) =
    (spec1===spec2) || (
    (spec1.hash == spec2.hash) &&
    isequal(spec1.energy, spec2.energy) &&
    isequal(spec1.properties, spec2.properties) &&
    isequal(spec1.counts, spec2.counts))
Base.isless(s1::Spectrum, s2::Spectrum) =
    isequal(s1[:Name], s2[:Name]) ? isless(s1.hash, s2.hash) : isless(s1[:Name], s2[:Name])
# Make it act like an AbstractVector
Base.eltype(::Spectrum{T}) where { T <: Real } = T
Base.length(spec::Spectrum) = length(spec.counts)
Base.ndims(spec::Spectrum) = 1
Base.size(spec::Spectrum) = size(spec.counts)
Base.axes(spec::Spectrum) = Base.axes(spec.counts)
Base.eachindex(spec::Spectrum) = eachindex(spec.counts)
Base.strides(spec::Spectrum) = strides(spec.counts)
Base.firstindex(spec::Spectrum) = firstindex(spec.counts)
Base.lastindex(spec::Spectrum) = lastindex(spec.counts)
Base.IndexStyle(::Type{<:Spectrum{T}}) where {T <: Real} = IndexStyle(Vector{T})
Base.IteratorSize(::Type{<:Spectrum{T}}) where {T <: Real} =
    Base.IteratorSize(Vector{T})
Base.IteratorEltype(::Type{<:Spectrum{T}}) where {T <: Real} =
    Base.IteratorEltype(Vector{T})
Base.iterate(spec::Spectrum, args...) = iterate(spec.counts, args...)
# Integer indices
Base.@propagate_inbounds Base.getindex(spec::Spectrum, I...) =
    getindex(spec.counts, I...)
Base.get(spec::Spectrum, idx::Int, def = zero(eltype(spec.counts))) =
    get(spec.counts, idx, def)
# AbstractFloat indices
Base.getindex(spec::Spectrum, energy::AbstractFloat) = getindex(spec.counts, channel(energy, spec))
Base.getindex(spec::Spectrum, sr::StepRangeLen{<:AbstractFloat}) =
    spec.counts[channel(sr[1], spec):channel(sr[end], spec)]
# CharXRay indices
Base.getindex(spec::Spectrum, cxr::CharXRay) = spec.counts[channel(energy(cxr), spec)]
# Symbol indices (ie Properties)
Base.get(spec::Spectrum, sym::Symbol, def::Any = missing) = get(spec.properties, sym, def)
Base.getindex(spec::Spectrum, sym::Symbol)::Any = getindex(spec.properties, sym)
# setindex!(...) for channel data
Base.setindex!(spec::Spectrum{T}, val::Real, idx...) where {T <: Real} =
    setindex!(spec.counts, convert(T, val), idx...)
# setindex!(...) for spectrum properties
Base.setindex!(spec::Spectrum, val::Real, sym::Symbol) =
    setindex!(spec.properties, val, sym)
Base.setindex!(spec::Spectrum, val::Any, sym::Symbol) =
    setindex!(spec.properties, val, sym)

Base.copy(spec::Spectrum) = Spectrum(spec.energy, copy(spec.counts), copy(spec.properties))
Base.merge!(spec::Spectrum, props::Dict{Symbol,Any}) = merge!(spec.properties, props)
properties(spec::Spectrum) = spec.properties
NeXLCore.name(spec::Spectrum) = spec[:Name]
Base.convert(::Type{Spectrum{T}}, spec::Spectrum{T}) where { T <: Real } = spec
function Base.convert(::Type{Spectrum{T}}, spec::Spectrum{U})::Spectrum{T} where {T <: Real, U <: Real }
    return Spectrum(spec.energy, T.(spec.counts), copy(spec.properties))
end
function Base.similar(spec::Spectrum{T}, ::Type{U}) where { T <: Real, U <: Real }
    return Spectrum(spec.energy, similar(spec.counts, U), copy(spec.properties))
end

# Spectrum math simply performs the operation on the channel data but doesn't change any properties except the name.
Base.:*(a::Real, s::Spectrum) = property!(Spectrum(s.energy, a*s.counts, copy(s.properties)), :Name, "$a⋅$(s[:Name])")
Base.:*(s::Spectrum, a::Real) = property!(Spectrum(s.energy, a*s.counts, copy(s.properties)), :Name, "$a⋅$(s[:Name])")
Base.:/(s::Spectrum, a::Real) = property!(Spectrum(s.energy, s.counts/a, copy(s.properties)), :Name, "$(s[:Name])/$a")
Base.:÷(s::Spectrum{T}, a::Integer) where {T <: Integer} = property!(Spectrum(s.energy, s.counts .÷ a, copy(s.properties)), :Name, "$(s[:Name])÷$a")
# Operations with two Spectrum arguments combine the common properties and rename.
function Base.:+(s1::Spectrum, s2::Spectrum)
    @assert isequal(s1.energy,s2.energy) "The energy calibration must be the same to add two spectra."
    r = Base.OneTo(min(length(s1.counts),length(s2.counts)))
    property!(Spectrum(s1.energy, view(s1.counts,r) + view(s2.counts,r), commonproperties(s1,s2)), :Name, "$(s1[:Name]) + $(s2[:Name])")
end
function Base.:-(s1::Spectrum, s2::Spectrum)
    @assert isequal(s1.energy,s2.energy) "The energy calibration must be the same to subtract one spectrum from another."
    r = Base.OneTo(min(length(s1.counts),length(s2.counts)))
    property!(Spectrum(s1.energy, view(s1.counts,r) - view(s2.counts,r), commonproperties(s1,s2)), :Name, "$(s1[:Name]) - $(s2[:Name])")
end

"""
    offset(s::Spectrum, dcounts::Real)

Returns a `Spectrum` like `s` but with dcounts added to each channel.
"""
offset(s::Spectrum, dcounts::Real) = property!(Spectrum(s.energy, s.counts .+ dcounts, copy(s.properties)), :Name, "$(s[:Name]) offset by $dcounts")

"""
    property!(spec::Spectrum, sym::Symbol, val::Any)

Useful to broadcast properties over many spectra.
"""
function property!(spec::Spectrum, sym::Symbol, val::Any)
    spec[sym]=val
    spec
end

"""
    rangeofenergies(ch::Integer, spec::Spectrum)

Returns the low and high energy extremes for the channels `ch`.
"""
rangeofenergies(ch::Integer, spec::Spectrum) =
    (energy(ch, spec.energy), energy(ch + 1, spec.energy))

"""
    NeXLCore.hasminrequired(ty::Type, spec::Spectrum)

Does `spec` have the necessary property items for the algorithm `ty`?
"""
NeXLCore.hasminrequired(ty::Type, spec::Spectrum) = hasminrequired(ty, spec.properties)

"""
    NeXLCore.requiredbutmissing(ty::Type, spec::Spectrum)

Which properties are missing from `spec` but are required for the algorithm `ty`?
"""
NeXLCore.requiredbutmissing(ty::Type, spec::Spectrum) =
    requiredbutmissing(ty, spec.properties)

maxproperty(specs::AbstractArray{<:Spectrum}, prop::Symbol) =
    maximum(spec -> spec[prop], specs)
minproperty(specs::AbstractArray{<:Spectrum}, prop::Symbol) =
    minimum(spec -> spec[prop], specs)
sameproperty(specs::AbstractArray{<:Spectrum}, prop::Symbol) =
    all(spec -> spec[prop] == specs[1][prop], specs) ? #
    specs[1][prop] : #
    error("The property $prop is not equal for all these spectra.")

function Base.show(io::IO, ::MIME"text/plain", spec::Spectrum{<:Real})
    comp = haskey(spec, :Composition) ? name(spec[:Composition]) : "Unknown"
    e0 =
        haskey(spec, :BeamEnergy) ? #
        "$(round(spec[:BeamEnergy]/1000.0,sigdigits=3))" :  #
        "Unknown"
    cnts = "$(round(sum(spec.counts),sigdigits=3))"
    if !(get(io, :compact, false) || haskey(io, :SHOWN_SET))
        textplot(io, spec, size = (12, 100))
        print(io, "\n")
    end
    print(
        io,
        "Spectrum{$(eltype(spec.counts))}[$(spec[:Name]), $(spec.energy), $(length(spec.counts)) ch, $e0 keV, $comp, $cnts counts]",
    )
end

Base.show(io::IO, spec::Spectrum) = show(io, "text/plain", spec)

function textplot(io::IO, spec::Spectrum; size = (16, 80))
    (rows, cols) = size
    # how much to plot
    e0_eV = haskey(spec, :BeamEnergy) ? min(spec[:BeamEnergy], energy(length(spec), spec)) : energy(length(spec), spec)
    maxCh = min(channel(convert(Float64, e0_eV), spec), length(spec))
    step, max = maxCh ÷ cols, maximum(spec.counts)
    maxes = [rows * (maximum(spec.counts[(i-1)*step+1:i*step]) / max) for i in 1:cols]
    for r in 1:rows
        print(io, join(map(i -> r ≥ rows - maxes[i] ? '*' : ' ', 1:cols)))
        if r == 1
            println(io, " $(round(max,digits=0))")
        elseif r == rows
            print(io, " $(round(0.001*e0_eV,digits=2)) keV")
        else
            println(io)
        end
    end
end

textplot(spec::Spectrum; size = (16, 80)) = textplot(stdout, spec, size = size)

"""
    NeXLUncertainties.asa(::Type{DataFrame}, spec::Spectrum; properties::Bool = false)

Converts the spectrum energy and counts data into a DataFrame.
"""
NeXLUncertainties.asa(
    ::Type{DataFrame},
    spec::Spectrum;
    properties::Bool = false,
)::DataFrame =
    properties ?
    DataFrame(
        Keys = [keys(spec.properties)...],
        Values = repr.([values(spec.properties)...]),
    ) : #
    DataFrame(E = energyscale(spec), I = counts(spec))

"""
    apply(spec::Spectrum, det::EDSDetector)::Spectrum

Applys the specified detector to this spectrum by ensuring that the energy scales match and spec[:Detector] = det.
Creates a copy of the original spectrum unless the detector is already det.
"""

function apply(spec::Spectrum, det::EDSDetector)::Spectrum
    if (haskey(spec, :Detector)) && (spec[:Detector] == det)
        return spec
    else
        props = copy(spec.properties)
        props[:Detector] = det
        return Spectrum(det.scale, spec.counts, props)
    end
end

"""
	matching(spec::Spectrum, resMnKa::Float64, lld::Int=1, minByFam::Dict{Shell,Element} = Dict())::BasicEDS

Build an EDSDetector to match the channel count and energy scale in this spectrum.
"""
matching(spec::Spectrum, resMnKa::Float64, lld::Int = 1)::BasicEDS =
    BasicEDS(length(spec), spec.energy, MnKaResolution(resMnKa), lld)

matching(
    spec::Spectrum,
    resMnKa::Float64,
    lld::Int,
    minByFam::Dict{Shell,Element},
)::BasicEDS = BasicEDS(length(spec), spec.energy, MnKaResolution(resMnKa), lld, minByFam)
"""
    matches(spec::Spectrum, det::Detector, tol::Float64 = 1.0)::Bool
    matches(spec1::Spectrum, spec2::Spectrum, tol::Float64 = 1.0)::Bool

Does the calibration of the Spectrum (approximately) match the calibration of the Detector
or the other Spectrum?
"""
matches(spec::Spectrum, det::Detector, tol::Float64 = 1.0)::Bool = false
function matches(spec::Spectrum, det::EDSDetector, tol::Float64 = 1.0)::Bool
    return abs(energy(1, spec) - energy(1, det)) < tol * channelwidth(1, det) &&
           abs(energy(length(spec), spec) - energy(length(spec), det)) <
           tol * channelwidth(length(spec), det)
end
function matches(spec1::Spectrum, spec2::Spectrum, tol::Float64 = 1.0)::Bool
    return abs(energy(1, spec1) - energy(1, spec2)) < tol * channelwidth(1, spec1) &&
           abs(energy(length(spec1), spec1) - energy(length(spec2), spec2)) <
           tol * channelwidth(length(spec1), spec1)
end

"""
    Base.keys(spec::Spectrum)

Return the defined properties as a set of `Symbol`s.
"""
Base.keys(spec::Spectrum) = keys(spec.properties)

"""
    Base.haskey(spec::Spectrum, sym::Symbol)

Is the specified key defined?
"""
Base.haskey(spec::Spectrum, sym::Symbol) = haskey(spec.properties, sym)

"""
    elms(spec::Spectrum, withcoating = false)::Set{Element}

Returns a set of the elements associated with this spectrum. `withcoating` determines whether the coating
elements are also added.
"""
function NeXLCore.elms(spec::Spectrum, withcoating = false)::Set{Element}
    res = Set{Element}()
    if haskey(spec, :Elements)
        union!(res, spec[:Elements])
    end
    if haskey(spec, :Composition)
        union!(res, keys(spec[:Composition]))
    end
    if haskey(spec, :Signature)
        union!(res, keys(spec[:Signature]))
    end
    if withcoating && haskey(spec, :Coating)
        union!(res, keys(material(spec[:Coating])))
    end
    return res
end

"""
    dose(spec::Spectrum, def=missing)
    dose(spec::Spectrum)

The probe dose in nano-amp seconds
"""
function dose(spec::Spectrum, def::Union{Float64,Missing})::Union{Float64,Missing}
    res = get(spec.properties, :LiveTime, missing) * get(spec, :ProbeCurrent, missing)
    return isequal(res, missing) ? def : res
end

function dose(spec::Spectrum)
    if haskey(spec.properties, :LiveTime) && haskey(spec.properties, :ProbeCurrent)
        return get(spec.properties, :LiveTime, missing) * get(spec, :ProbeCurrent, missing)
    else
        @error "One or more properties necessary to calculate the dose is missing.\n  You must provide both :ProbeCurrent and :LiveTime."
    end
end

"""
    NeXLCore.energy(ch::Int, spec::Spectrum)::Float64

What energy is associated with `ch` in `spec`?
"""
NeXLCore.energy(ch::Int, spec::Spectrum)::Float64 = energy(ch, spec.energy)

"""
    coordinate(spec::Spectrum)
    coordinate(hss::Hyperspectrum, idx::NTuple{N, <Integer})
    coordinate(hss::HyperSpectrum, ci::CartesianIndex)

Returns the stage coordinate corresponding most closely to the :StagePosition (compensated
for the image offset in a `HyperSpectrum`) Defaults to :X=>0, :Y=>0

Nominally coordinates are:  :X, :Y, :Z, :R, :T, :B

Use with: `image2stage(...)` and `stage2image(...)` to calculate the stage coordinate of image pixels.
"""
function coordinate(spec::Spectrum)
    return get(spec, :StagePosition, Dict{Symbol,Float64}(:X=>0.0, :Y=>0.0))
end

"""
    channelwidth(ch::Int, spec::Spectrum)::Float64
    channelwidth(ch::Int, spec::HyperSpectrum)::Float64

Returns the width of the `ch` channel in eV.
"""
channelwidth(ch::Int, spec::Spectrum) = energy(ch + 1, spec) - energy(ch, spec)

channel(eV::Float64, spec::Spectrum)::Int = channel(eV, spec.energy)

"""
    counts(spec::Spectrum, ::Type{T}, applyLLD=false)::Vector{T} where {T<:Number}
    counts(spec::Spectrum, channels::AbstractUnitRange{<:Integer}, ::Type{T}, applyLLD=false)::Vector{T} where {T<:Number}

Creates a copy of the spectrum counts data as the specified Number type. If the spectrum has a :Detector
property then the detector's lld (low-level discriminator) and applyLLD=true then the lld is applied to the result
by setting all channels less-than-or-equal to det.lld to zero. This method throws an error if the counts data
cannot be converted without loss to the type `T`.
"""
function counts(
    spec::Spectrum,
    ::Type{T} = Float64,
    applyLLD = false,
)::Array{T} where {T<:Number}
    res = T.(spec.counts)
    if applyLLD && haskey(spec, :Detector)
        fill!(view(res, 1:lld(spec[:Detector])), zero(T))
    end
    return res
end
function counts(
    spec::Spectrum,
    channels::AbstractRange{<:Integer},
    ::Type{T} = Float64,
    applyLLD = false,
)::Array{T} where {T<:Real}
    if (first(channels) >= 1) && (last(channels) <= length(spec))
        res = T.(spec.counts[channels])
    else
        res = zeros(T, length(channels))
        r = max(first(channels), 1):min(last(channels), length(spec))
        res[first(r)-first(channels)+1:last(r)-first(channels)+1] = spec.counts[r]
    end
    if applyLLD && haskey(spec, :Detector) && (lld(spec) <= first(channels))
        fill!(view(res, 1:lld(spec)-first(channels)+1), zero(T))
    end
    return res
end
function counts(
    spec::Spectrum,
    channels::AbstractVector{<:Integer},
    ::Type{T} = Float64,
    applyLLD = false,
)::Array{T} where {T<:Real}
    res = T.(spec.counts[channels])
    if applyLLD && haskey(spec, :Detector) && (lld(spec) <= first(channels))
        fill!(view(res, 1:lld(spec)-first(channels)+1), zero(T))
    end
    return res
end

"""
    lld(spec::Spectrum)

Gets the low-level discriminator associated with this spectrum if there is one.
"""
lld(spec::Spectrum) =
    haskey(spec.properties, :Detector) ? lld(spec.properties[:Detector]) : 1

"""
	Base.findmax(spec::Spectrum, chs::AbstractRange{<:Integer})
    Base.findmax(spec::Spectrum)

Returns the (maximum intensity, channel index) over the specified range of channels
"""
function Base.findmax(spec::Spectrum, chs::AbstractRange{<:Integer})
    max = findmax(spec.counts[chs])
    return (max[1] + first(chs) - 1, max[2])
end
function Base.findmax(spec::Spectrum)
    last = min(
        haskey(spec, :BeamEnergy) ? channel(spec[:BeamEnergy], spec) : length(spec.counts),
        length(spec.counts),
    )
    return findmax(spec.counts[lld(spec):last])
end

"""
    integrate(spec::Spectrum, channels::AbstractUnitRange{<:Integer})
    integrate(spec::Spectrum, energyRange::StepRangeLen{Float64})

Sums all the counts in the specified channels.  No background correction.

    integrate(spec::Spectrum, back1::AbstractUnitRange{<:Integer}, peak::AbstractUnitRange{<:Integer}, back2::AbstractUnitRange{<:Integer})::UncertainValue
    integrate(spec::Spectrum, back1::StepRangeLen{Float64}, peak::StepRangeLen{Float64}, back2::StepRangeLen{Float64})::UncertainValue

Sums all the counts in the specified channels with background correction using the background intervals.

    integrate(spec::Spectrum)

Total integral of all counts from the LLD to the beam energy

"""
integrate(spec::Spectrum, channels::AbstractUnitRange{<:Integer}) =
    sum(spec.counts[channels])

integrate(spec::Spectrum, energyRange::StepRangeLen{Float64}) =
    sum(spec.counts[channel(energyRange[1], spec):channel(energyRange[end], spec)])

function integrate(
    spec::Spectrum,
    back1::AbstractUnitRange{<:Integer},
    peak::AbstractUnitRange{<:Integer},
    back2::AbstractUnitRange{<:Integer},
)::UncertainValue
    bL, bH = mean(spec.counts[back1]), mean(spec.counts[back2])
    cL, cH, cP = mean(back1), mean(back2), mean(peak)
    iP, w, t = sum(spec.counts[peak]), length(peak), ((cP - cL) / (cH - cL))
    return uv(
        iP - w * (bL * (1.0 - t) + bH * t),
        sqrt(iP + (sqrt(bL) * w * (1.0 - t))^2 + (sqrt(bH) * w * t)^2),
    )
end

function integrate(
    spec::Spectrum,
    back1::StepRangeLen{Float64},
    peak::StepRangeLen{Float64},
    back2::StepRangeLen{Float64},
)::UncertainValue
    b1i = channel(back1[1], spec):channel(back1[end], spec)
    b2i = channel(back2[1], spec):channel(back2[end], spec)
    pi = channel(peak[1], spec):channel(peak[end], spec)
    integrate(spec, b1i, pi, b2i)
end

function integrate(spec::Spectrum)
    last = min(
        haskey(spec, :BeamEnergy) ? channel(spec[:BeamEnergy], spec) : length(spec.counts),
        length(spec.counts),
    )
    return integrate(spec, lld(spec):last)
end

"""
    kratio(unk::Spectrum, std::Spectrum, back1::AbstractUnitRange{<:Integer}, peak::AbstractUnitRange{<:Integer}, back2::AbstractUnitRange{<:Integer})::UncertainValue
    kratio(unk::Spectrum, std::Spectrum, back1::StepRangeLen{Float64}, peak::StepRangeLen{Float64}, back2::StepRangeLen{Float64})::UncertainValue

The k-ratio of unk relative to std corrected for dose.  Requires that `unk` and `std` have the properties
`:LiveTime` and `:ProbeCurrent` defined.
"""
function kratio(
    unk::Spectrum,
    std::Spectrum,
    back1::AbstractUnitRange{<:Integer},
    peak::AbstractUnitRange{<:Integer},
    back2::AbstractUnitRange{<:Integer},
)
    iunk, istd = integrate(unk, back1, peak, back2) / dose(unk),
    integrate(std, back1, peak, back2) / dose(std)
    f = NeXLUncertainties.value(iunk) / NeXLUncertainties.value(istd)
    return uv(
        f,
        abs(f) * sqrt(
            (σ(iunk) / NeXLUncertainties.value(iunk))^2 +
            (σ(istd) / NeXLUncertainties.value(istd))^2,
        ),
    )
end

kratio(
    unk::Spectrum,
    std::Spectrum,
    back1::StepRangeLen,
    peak::StepRangeLen,
    back2::StepRangeLen,
) = kratio(
    unk,
    std, #
    channel(back1[1], unk):channel(back1[end], unk), #
    channel(peak[1], unk):channel(peak[end], unk), #
    channel(back2[1], unk):channel(back2[end], unk),
)

"""
    energyscale(spec::Spectrum)
    energyscale(spec::Spectrum, chs::AbstractRange{<:Integer})

Returns an array with the bin-by-bin energies
"""
energyscale(spec::Spectrum) = energyscale(spec.energy, eachindex(spec))
energyscale(spec::Spectrum, chs::AbstractRange{<:Integer}) = energyscale(spec.energy, chs)

"""
    simpleEDS(spec::Spectrum, fwhmatmnka::Float64)

Build a `EDSDetector` object for this spectrum with the specified FWHM at Mn Kα.
"""
simpleEDS(spec::Spectrum, fwhmatmnka::Float64) =
    BasicEDS(length(spec), spec.energy, MnKaResolution(fwhmatmnka))

"""
    subsample(spec::Spectrum, frac::Float64)

Subsample the counts data in a spectrum according to a statistically valid algorithm.  Returns
`spec` if frac>=1.0.
"""
function subsample(spec::Spectrum, frac::Float64)::Spectrum
    @assert frac > 0.0 "frac must be larger than zero."
    @assert frac <= 1.0 "frac must be less than or equal to 1.0."
    @assert haskey(spec, :LiveTime) "Please specify a :LiveTime in subsample"
    ss(n, f) = n > 0 ? sum(rand() <= f ? 1 : 0 for _ in 1:n) : 0
    frac = max(0.0, min(frac, 1.0))
    if frac < 1.0
        props = deepcopy(spec.properties)
        props[:LiveTime] = frac * spec[:LiveTime] # Must have
        if haskey(spec, :RealTime)
            props[:RealTime] = frac * spec[:RealTime] # Might have
        end
        return Spectrum(spec.energy, map(n -> ss(floor(Int, n), frac), spec.counts), props)
    else
        @warn "Not actually subsampling the spectrum because $frac > 1.0."
        return spec
    end
end
"""
    subdivide(spec::Spectrum, n::Int)::Vector{Spectrum}

Splits the event data in one spectrum into n spectra by assigning each event
to a pseudo-random choice of one of the n result spectra.  Produces n spectra
that act as though the original spectrum was collected in n time intervals
of LiveTime/n.  This is quite slow because it needs to call rand() for each
count in the spectrum (not just each channel).
"""
function subdivide(spec::Spectrum{T}, n::Int)::Vector{Spectrum{T}} where {T <: Real}
    res = zeros(Int, n, length(spec.counts))
    for ch in eachindex(spec.counts)
        # Assign each event to one and only one detector
        si = rand(1:n, floor(Int, spec[ch]))
        for i in 1:n
            res[i, ch] = count(e -> e == i, si)
        end
    end
    return map(1:n) do i
        props = deepcopy(spec.properties)
        props[:Name] = "Sub[$(spec[:Name]),$(i) of $(n)]"
        props[:LiveTime] = spec[:LiveTime] / n # Must have
        if haskey(spec, :RealTime)
            props[:RealTime] = spec[:RealTime] / n # Might have
        end
        Spectrum(spec.energy, res[i, :], props)
    end
end

"""
    estimatebackground(data::AbstractArray{Float64}, channel::Int, width::Int=5, order::Int=2)

Returns the tangent to the a quadratic fit to the counts data centered at channel with width
"""
function estimatebackground(
    data::AbstractArray{Float64},
    channel::Int,
    width::Int = 5,
    order::Int = 2,
)::ImmutablePolynomial
    minCh, maxCh = max(1, channel - width), min(length(data), channel + width)
    if maxCh - minCh >= order
        fr = fit(
            ImmutablePolynomial,
            (minCh-channel):(maxCh-channel),
            data[minCh:maxCh],
            order,
        )
        return ImmutablePolynomial([fr(0), derivative(fr)(0)]) # Linear
    else
        return ImmutablePolynomial([mean(data[minCh:maxCh]), 0.0])
    end
end

"""
    modelBackground(spec::Spectrum, chs::AbstractUnitRange{<:Integer}, ash::AtomicSubShell)

spec: A spectrum containing a peak centered on chs
chs:  A range of channels containing a peak
ash:  The edge (as an AtomicSubShell)

A simple model for modeling the background under a characteristic x-ray peak. The model
fits a line to low and high energy background regions around firsr(chs) and last(chs). If
the low energy line extended out to the edge energy is larger than the high energy line
at the same energy, then a negative going edge is fit between the two. Otherwise a line
is fit between the low energy side and the high energy side. This model only works when
there are no peak interference over the range chs.

    modelBackground(spec::Spectrum, chs::AbstractUnitRange{<:Integer})
    modelBackground(spec::Spectrum, chs::AbstractUnitRange{<:Integer}, ash::AtomicSubShell)

`spec`: A spectrum containing a peak centered on chs
`chs`:  A range of channels containing a peak
`ash`:  The largest edge within the range of channels `chs` associated with the characteristic peak

A simple model for modeling the background under a characteristic x-ray peak. The model
fits a line between the  low and high energy background regions around first(chs) and last(chs).
This model only works when there are no peak interference over the range chs.
"""
function modelBackground(
    spec::Spectrum,
    chs::AbstractUnitRange{<:Integer},
    ash::AtomicSubShell,
)
    cnts, ec = counts(spec), channel(energy(ash), spec)
    bl, bh = estimatebackground(cnts, first(chs), 5), estimatebackground(cnts, last(chs), 5)
    # bh = ch-> mean(cnts[last(chs):min(length(cnts),last(chs)+5)])
    if (ec < last(chs)) &&
       (bl(ec - first(chs)) > bh(ec - last(chs))) &&
       (energy(ash) < 2.0e3)  && ( shell(ash.subshell)==KShell )
        res = zeros(Float64, length(chs))
        res[1:ec-first(chs)] .= (bl(y - 1) for y in 1:ec-first(chs))
        res[ec-first(chs)+1] = 0.5 * (bl(ec - first(chs) + 1) + bh(ec - last(chs)))
        res[ec-first(chs)+2:last(chs)-first(chs)+1] .=
            (bh(y - last(chs)) for y in ec+1:last(chs))
    else
        s = (bh(0) - bl(0)) / length(chs)
        back = ImmutablePolynomial([bl(0), s])
        res = back.(collect(0:length(chs)-1))
    end
    return res
end

function modelBackground(spec::Spectrum, chs::AbstractUnitRange{<:Integer})
    cxs = counts(spec)
    bl = estimatebackground(cxs, first(chs), 5)
    bh = estimatebackground(cxs, last(chs), 5)
    s = (bh(0) - bl(0)) / length(chs)
    back = ImmutablePolynomial([bl(0), s])
    return back.(collect(0:length(chs)-1))
end

"""
    extractcharacteristic(spec::Spectrum, lowBack::AbstractUnitRange{<:Integer}, highBack::AbstractUnitRange{<:Integer})::Vector{Float64}

Extract the characteristic intensity for the peak located within chs with an edge at ash.
"""
function extractcharacteristic(
    spec::Spectrum,
    chs::AbstractUnitRange{<:Integer},
    ash::AtomicSubShell,
)::Vector{Float64}
    return counts(spec, chs, Float64) - modelBackground(spec, chs, ash)
end

"""
    peak(spec::Spectrum, chs::AbstractUnitRange{<:Integer}, ash::AtomicSubShell)::Float64

Estimates the peak intensity for the characteristic X-ray in the specified range of channels.
"""
peak(spec::Spectrum, chs::AbstractUnitRange{<:Integer})::Float64 =
    return sum(counts(spec, chs, Float64)) - background(spec, chs)

"""
    background(spec::Spectrum, chs::AbstractUnitRange{<:Integer}, ash::AtomicSubShell)::Float64

Estimates the background intensity for the characteristic X-ray in the specified range of channels.
"""
background(spec::Spectrum, chs::AbstractUnitRange{<:Integer})::Float64 =
    sum(modelBackground(spec, chs))

"""
    peaktobackground(spec::Spectrum, chs::AbstractUnitRange{<:Integer}, ash::AtomicSubShell)::Float64

Estimates the peak-to-background ratio for the characteristic X-ray intensity in the specified range of channels
which encompass the specified AtomicSubShell.
"""
function peaktobackground(
    spec::Spectrum,
    chs::AbstractUnitRange{<:Integer},
    ash::AtomicSubShell,
)::Float64
    back = sum(modelBackground(spec, chs, ash))
    return (sum(counts(spec, chs, Float64)) - back) / back
end

"""
    estkratio(unk::Spectrum, std::Spectrum, chs::AbstractUnitRange{<:Integer})

Estimates the k-ratio from niave models of peak and background intensity.  Only works if the peak is not interfered.
"""
estkratio(unk::Spectrum, std::Spectrum, chs::AbstractUnitRange{<:Integer}) =
    peak(unk, chs) * dose(std) / (peak(std, chs) * dose(unk))

"""
    NeXLUncertainties.asa(::Type{DataFrame}, spec::AbstractArray{Spectrum})::DataFrame

Returns a DataFrame that summarizes the list of spectra.
"""
function NeXLUncertainties.asa(
    ::Type{DataFrame},
    specs::AbstractArray{<:Spectrum},
)::DataFrame
    _asname(comp) = ismissing(comp) ? missing : name(comp)
    unf, unl, uns = Union{Float64,Missing}, Union{Film,Nothing}, Union{String,Missing}
    nme, e0, pc, lt, rt, coat, integ, comp =
        String[], unf[], unf[], unf[], unf[], unl[], Float64[], uns[]
    for spec in specs
        push!(nme, spec[:Name])
        push!(e0, get(spec, :BeamEnergy, missing))
        push!(pc, get(spec, :ProbeCurrent, missing))
        push!(lt, get(spec, :LiveTime, missing))
        push!(rt, get(spec, :RealTime, missing))
        push!(coat, get(spec, :Coating, nothing))
        push!(integ, integrate(spec))
        push!(comp, _asname(get(spec, :Composition, missing)))
    end
    return DataFrame(
        Name = nme,
        BeamEnergy = e0,
        ProbeCurrent = pc,
        LiveTime = lt,
        RealTime = rt,
        Coating = coat,
        Integral = integ,
        Material = comp,
    )
end

"""
    NeXLUncertainties.asa(::Type{DataFrame}, spec::AbstractDict{Element, Spectrum})::DataFrame

Returns a DataFrame that summarizes a dictionary of standard spectra.
"""
function NeXLUncertainties.asa(
    ::Type{DataFrame},
    stds::AbstractDict{Element,Spectrum},
)::DataFrame
    _asname(comp) = ismissing(comp) ? missing : name(comp)
    unf, unl, uns = Union{Float64,Missing}, Union{Film,Nothing}, Union{String,Missing}
    elm, zs, mfs, nme, e0, pc, lt, rt, coat, integ, comp = String[],
    Int[],
    Float64[],
    String[],
    unf[],
    unf[],
    unf[],
    unf[],
    unl[],
    Float64[],
    uns[]
    for el in sort(collect(keys(stds)))
        spec = stds[el]
        push!(elm, el.symbol)
        push!(zs, el.number)
        push!(nme, spec[:Name])
        push!(mfs, ismissing(spec[:Composition]) ? missing : spec[:Composition][el])
        push!(e0, get(spec, :BeamEnergy, missing))
        push!(pc, get(spec, :ProbeCurrent, missing))
        push!(lt, get(spec, :LiveTime, missing))
        push!(rt, get(spec, :RealTime, missing))
        push!(coat, get(spec, :Coating, nothing))
        push!(integ, integrate(spec))
        push!(comp, _asname(get(spec, :Composition, missing)))
    end
    return DataFrame(
        Element = elm,
        Z = zs,
        Name = nme,
        Material = comp,
        MassFrac = mfs,
        BeamEnergy = e0,
        ProbeCurrent = pc,
        LiveTime = lt,
        RealTime = rt,
        Coating = coat,
        Integral = integ,
    )
end

"""
    details(io, spec::Spectrum)
    details(spec::Spectrum)

Outputs a description of the data in the spectrum.
"""
function details(io::IO, spec::Spectrum)
    println(io, "           Name:   $(spec[:Name])")
    println(io, "    Calibration:   $(spec.energy)")
    println(io, "    Beam energy:   $(get(spec, :BeamEnergy, missing)/1000.0) keV")
    println(io, "  Probe current:   $(get(spec, :ProbeCurrent, missing)) nA")
    println(io, "      Live time:   $(get(spec, :LiveTime, missing)) s")
    println(io, "        Coating:   $(get(spec,:Coating, "None"))")
    println(io, "       Detector:   $(get(spec, :Detector, missing))")
    println(io, "        Comment:   $(get(spec, :Comment, missing))")
    println(io, "       Integral:   $(integrate(spec)) counts")
    comp = get(spec, :Composition, missing)
    if !ismissing(comp)
        println(io, "    Composition:   $(comp)")
        det = get(spec, :Detector, missing)
        if !ismissing(det)
            coating = get(spec, :Coating, missing)
            comp2 = collect(keys(comp))
            if !ismissing(coating)
                append!(comp2, keys(coating.material))
            end
            for elm1 in keys(comp)
                for ext1 in extents(elm1, det, 1.0e-4)
                    print(io, "            ROI:   $(elm1.symbol)[$(ext1)]")
                    intersects = []
                    for elm2 in comp2
                        if elm2 ≠ elm1
                            for ext2 in extents(elm2, det, 1.0e-4)
                                if length(intersect(ext1, ext2)) > 0
                                    push!(intersects, "$(elm2.symbol)[$(ext2)]")
                                end
                            end
                        end
                    end
                    if length(intersects) > 0
                        println(io, " intersects $(join(intersects,", "))")
                    else
                        p, b = peak(spec, ext1), background(spec, ext1)
                        σ = p / sqrt(b)
                        println(
                            io,
                            " = $(round(Int,p)) counts over $(round(Int,b)) counts - σ = $(round(Int,σ))",
                        )
                    end
                end
            end
        end
    end
    return nothing
end

details(spec::Spectrum) = details(stdout, spec)

"""
    commonproperties(specs::AbstractArray{Spectrum})
    commonproperties(props1::Dict{Symbol,Any}, props2::Dict{Symbol,Any})

Return the properties that are held in common by all the spectra.
"""
function commonproperties(props1::Dict{Symbol,Any}, props2::Dict{Symbol,Any})
    res = Dict{Symbol,Any}()
    for (key, sp1) in props1
        if isequal(sp1, get(props2, key, missing))
            res[key] = sp1
        end
    end
    return res
end
function commonproperties(specs::Vararg{<:Spectrum})
    props = specs[1].properties
    for sp2 in specs[2:end]
        props = commonproperties(props, sp2.properties)
    end
    return props
end
commonproperties(specs::AbstractArray{<:Spectrum}) = commonproperties(specs...)

"""
    maxspectrum(specs::AbstractArray{Spectrum{T}})::Spectrum{T} where T<:Real

Compute the max-pixel spectrum for the specified spectra.
"""
function maxspectrum(specs::AbstractArray{Spectrum{T}})::Spectrum{T} where {T<:Real}
    res = zeros(T, maximum(length.(specs)))
    for spec in specs
        for i in eachindex(spec)
            res[i] = max(res[i], spec[i])
        end
    end
    props = commonproperties(specs)
    props[:Name] = "MaxSpectrum[$(length(specs)) spectra]"
    return Spectrum(spec1.energy, res, props)
end

"""
    minspectrum(specs::AbstractArray{Spectrum{T}})::Spectrum{T} where T<:Real

Compute the min-pixel spectrum for the specified spectra.
"""
function minspectrum(specs::AbstractArray{Spectrum{T}})::Spectrum{T} where {T<:Real}
    res = zeros(T, maximum(length.(specs)))
    fill!(res, typemax(T))
    for spec in specs
        for i in eachindex(spec)
            res[i] = min(res[i], spec[i])
        end
    end
    props = commonproperties(specs)
    props[:Name] = "MinSpectrum[$(length(specs)) spectra]"
    return Spectrum(spec1.energy, res, props)
end

"""
    Base.sum(specs::AbstractArray{Spectrum{T}}; restype::Type{<:Real}=T, applylld=false, name=Nothing|String)::Spectrum{T}

Computes the sum spectrum over an `AbstractArray{Spectrum}` where the :ProbeCurrent and :LiveTime will be maintained
in a way that maintains the sum of the individual doses.  This function assumes (but does not check) that the energy
scales are equivalent for all the spectra.  The resultant energy scale is the scale of the first spectrum.  Other than
:ProbeCurrent, :LiveTime and :RealTime, only the properties that the spectra hold in common will be maintained.
"""
function Base.sum(
    specs::AbstractArray{Spectrum{T}},
    restype::Type{<:Real} = T;
    applylld = false,
    name = nothing,
)::Spectrum where {T<:Real}
    cxs = zeros(restype, maximum(length.(specs)))
    for spec in specs
        cxs .+= counts(spec, eachindex(cxs), restype, applylld)
    end
    props = commonproperties(specs)
    props[:ProbeCurrent] = mean(sp[:ProbeCurrent] for sp in specs)
    props[:LiveTime] = sum(dose.(specs)) / props[:ProbeCurrent]
    rt = sum(get(sp, :RealTime, NaN64) for sp in specs)
    if !isnan(rt)
        props[:RealTime] = rt
    end
    props[:Name] = something(name, "Sum[$(length(specs)) spectra]")
    return Spectrum(specs[1].energy, cxs, props)
end

"""
    χ²(s1::Spectrum{T}, s2::Spectrum{T}, chs)::T where {T<:Real}
    χ²(specs::AbstractArray{Spectrum{T}}, chs)::Matrix{T}

Computes the dose corrected χ² metric of similarity between `s1` and `s2`.
Or computes the matrix of χ² for the array of spectra.
"""
function χ²(s1::Spectrum{T}, s2::Spectrum{T}, chs)::T where {T<:Real}
    k1, k2 = 1.0 / dose(s1), 1.0 / dose(s2)
    d = sum((k1 * s1[ch] - k2 * s2[ch])^2 for ch in chs)
    s = sum(sqrt(k1^2 * max(one(T), s1[ch]) + k2^2 * max(one(T), s2[ch])) for ch in chs)
    return d / s
end

function χ²(specs::AbstractVector{Spectrum{T}}, chs)::Matrix{T} where {T<:Real}
    χ2s = zeros(Float64, (length(specs), length(specs)))
    for i in eachindex(specs), j in i+1:length(specs)
        χ2s[i, j] = χ²(specs[i], specs[j], chs)
        χ2s[j, i] = χ2s[i, j]
    end
    return χ2s
end

"""
    measure_dissimilarity(specs::AbstractArray{Spectrum{T}}, chs)::Vector{Float64}
    measure_dissimilarity(specs::AbstractArray{Spectrum}, minE::Float64=100.0)::Vector{Float64}
    measure_dissimilarity(specs::AbstractArray{Spectrum{T}}, det::Detector, mat::Material)::Vector{Float64}

Returns a vector of χ²s which measure how different the i-th `Spectrum` is from the other spectra.
The first version covers all the channels between minE and the nominal beam energy. The second version
only considers those channels representing peaks in a spectrum from the `Material` on the `Detector`.
The statistic returned is the χ² for each spectrum relative to the mean spectrum.  It should be nominally
a vector of 1's for spectra that differ only by count statistics.
"""
function measure_dissimilarity(
    specs::AbstractArray{<:Spectrum},
    chs,
)::Vector{Float64}
    s = sum(specs)
    return [χ²(spec, s, chs) for spec in specs]
end

function measure_dissimilarity(
    specs::AbstractArray{<:Spectrum},
    minE::Float64 = 100.0,
)::Vector{Float64}
    e0 = maximum(spec[:BeamEnergy] for spec in specs)
    chs =
        minimum(
            channel(minE, spec) for spec in specs
        ):maximum(channel(e0, spec) for spec in specs)
    return measure_dissimilarity(specs, chs)
end

function measure_dissimilarity(
    specs::AbstractArray{<:Spectrum},
    det::Detector,
    mat::Material,
)::Vector{Float64}
    function mrg(inp::Vector{UnitRange{Int}})::Vector{UnitRange{Int}}
        simp = sort(inp)
        st, res = simp[1], UnitRange{Int}[]
        for r in simp[2:end]
            if isempty(intersect(st, r))
                push!(res, st)
                st = r
            else
                st = first(st):max(last(r), last(st))
            end
        end
        push!(res, st)
        return res
    end
    e0 = maximum(spec[:BeamEnergy] for spec in specs)
    # Figure out the contiguous ROIs and the channels in the ROIs
    rois = mrg(
        mapreduce(
            elm -> extents(characteristic(elm, alltransitions, 0.01, e0), det, 0.001),
            append!,
            keys(mat),
        ),
    )
    chs = mapreduce(collect, append!, rois)
    return measure_dissimilarity(specs, chs)
end

"""
    findsimilar(specs::AbstractArray{Spectrum{T}}; tol = 1.8, minspecs=3)::Vector{Spectrum{T}}

Filters a collection of spectra for the ones most similar to the average by
removing the least similar spectrum sequentially until all the remaining spectra are within
tol times the mean(χ²).  This is useful for finding which of a set of replicate spectra are
sufficiently similar to each other.
"""
function findsimilar(
    specs::AbstractArray{Spectrum{T}};
    tol = 1.8,
    minspecs = 3,
)::Vector{Spectrum{T}} where {T<:Real}
    if length(specs) > minspecs
        σs = measure_dissimilarity(specs)
        (fmσ, fmi) = findmax(σs)
        rem = filter(σ -> σ ≠ fmσ, σs)
        # Now perform this recursively until all are within tol or we hit minspecs
        if fmσ > tol * Statistics.mean(rem)
            return findsimilar(
                filter(s -> s != specs[fmi], specs),
                tol = tol,
                minspecs = minspecs,
            )
        end
    end
    return specs
end

"""
    duane_hunt(spec::Spectrum, def = spec[:BeamEnergy])

Estimates the Duane-Hunt limit (the energy at which the continuum goes to effectively zero.)
"""
function duane_hunt(spec::Spectrum, def = spec[:BeamEnergy])
    if channel(spec[:BeamEnergy], spec) < length(spec)
        este0 = let # Sometimes the D-H can be much below `spec[:BeamEnergy]`
            ed, cd = energyscale(spec), counts(spec)
            lim = min(10.0, max(1.0, 1.0e-5 * maximum(cd)))
            len = min(10, max(5, length(spec) ÷ 400))
            # this handles pulse pile-up above the D-H
            ed[findlast(map(i -> mean(cd[i:i+len]), eachindex(cd[1:end-len])) .> lim)]
        end
        xdata, ydata = let
            chs = let # Range of channels above/below `este0`
                ch0 = channel(este0, spec)
                max(1,(19*ch0)÷20):min((21*ch0)÷20,length(spec))
            end
            energyscale(spec)[chs], counts(spec, chs)
        end
        esti0 = mean(ydata[1:min(10,length(ydata))]) * xdata[1] / (este0 - xdata[1])
        model(xs, p) = map(x -> p[1] * bremsstrahlung(Kramers1923, x, p[2], n"H"), xs)
        curve_fit(model, xdata, ydata, [esti0, este0]).param[2]
    else
        def
    end
end

"""
    uv(spec::Spectrum, chs::AbstractRange{<:Integer}=eachindex(spec))::Vector{UncertainValue}

Converts the count's data in a spectrum into an Vector{UncertainValue} assuming
count statistics can be approximated by C ± √C.
"""
function NeXLUncertainties.uv(spec::Spectrum, chs::AbstractRange{<:Integer}=eachindex(spec))::Vector{UncertainValue}
    c = counts(spec, chs)
    return uv.(c, sqrt.(max.(c,one(eltype(spec)))))
end

"""
    sigma(spec::Spectrum, specs::AbstractArray{<:Spectrum}, chs::AbstractRange{<:Integer})::Vector{Float64}

Computes on a channel-by-channel basis how much `spec` spectrum deviates from the mean of the
other spectra in `specs`.  The result is expressed in terms of the standard deviation expected
from count statistics alone.   Assuming `spec` varies only by count statistics we expect
the result values have a mean 0.0 and a standard deviation of 1.0.
"""
function sigma(spec::Spectrum, specs::AbstractArray{<:Spectrum}, chs::AbstractRange{<:Integer})::Vector{Float64}
    function doseaverage(specs, chs)
        t = [ uv(spec, chs)/dose(spec) for spec in specs ]
        return map(j->mean(collect(t[i][j] for i in eachindex(t))), eachindex(t[1]))
    end
    function delta(spec, specs, chs)
        minus(uv1,uv2) = uv(value(uv1)-value(uv2),sqrt(variance(uv1)+variance(uv2)))
        return minus.(uv(spec, chs), dose(spec)*doseaverage(filter(s->s!=spec, specs), chs))
    end
    return map(v->value(v)/σ(v), delta(spec, specs, chs))
end

"""
    dosenormalize(spectrum::Spectrum{T}, dose=60.0)::Spectrum{T} where { T <: Real }
    dosenormalize(spectrum::Spectrum{T}, dose=60.0)::Spectrum{Float64} where { T <: Integer }

Compute a spectrum which is `spectrum` rescaled to a live time times probe current equal to `dose`.
Useful for setting spectra on an equivalent acquisition duration scale.
"""
function dosenormalize(spectrum::Spectrum{T}, dose=60.0)::Spectrum{T} where { T <: AbstractFloat }
    res = copy(spectrum)
    scale = dose / NeXLSpectrum.dose(res)
    res.counts .*= scale
    res[:LiveTime] *= scale
    res[:Name] = "N[$(spectrum[:Name]), $dose nA⋅s]"
    return res
end
function dosenormalize(spectrum::Spectrum{T}, dose=60.0)::Spectrum{Float64} where { T <: Integer }
    scale = dose / NeXLSpectrum.dose(spectrum)
    newProps = copy(spectrum.properties)
    newProps[:LiveTime] *= scale
    newProps[:Name] = "N[$(spectrum[:Name]), $dose nA⋅s]"
    return Spectrum(spectrum.energy, Float64.(spectrum.counts) * scale, newProps)
end

"""
    recalibrate(s::Spectrum{T}, es::LinearEnergyScale)

Allows changing the energy scale on a spectrum from one LinearEnergyScale to another as though the spectrum were
measured on a different detector.  The algorith uses a FFT-base scheme to rescale and shift the spectral data.
Ths scheme allows for fractional shifts of offset and fractional changes in the width.  It is limited in that
the change in width must produce an integral number of channels in the resulting spectrum.  The algorithm
maintains the total spectrum integral so the new spectrum can be used for quantitative purposes.
Plotting one spectrum over the other should maintain peak position but is likely to change the channel counts.
"""
function recalibrate(s::Spectrum{T}, es::LinearEnergyScale) where { T <: Real }
    (!(s.energy isa LinearEnergyScale)) && error("The rescale(...) function requires that the spectrum has a LinearEnergyScale.")
    # Resample the spectrum onto a new channel width (offset remains same)
    cxs, oldsum = counts(s), sum(s.counts)
    oldlen, newlen = length(cxs), round(Int, (s.energy.width/es.width)*length(cxs))
    if newlen ≠ length(cxs)
        cxs = FourierTools.resample(cxs, newlen)
    end
    newWidth = s.energy.width*(oldlen/newlen)
    # Shift the first channel from s.energy.offset to es.offset
    shft = (s.energy.offset - es.offset) / newWidth
    if abs(shft) > 0.01
        cxs = FourierTools.shift(cxs, shft)
    end
    return Spectrum(es, cxs*(oldsum/sum(cxs)), copy(s.properties))
end

"""
    shift(s::Spectrum, ev::AbstractFloat)::Spectrum

Shift the entire spectrum along the energy axis by a specified number of ev by modifying the
counts data.  This function can shift by a fractional number of channels.
"""
function shift(s::Spectrum, ev::AbstractFloat)::Spectrum
    return Spectrum(s.energy, FourierTools.shift(counts(s), ev/s.energy.width), copy(s.properties))
end


text:
module WeaklyHard

#= Uses for all includes etc =#
using LinearAlgebra

#= Types and constructors =#
include("types.jl")

#= Includes =#
include("alphabet.jl")          # Defines the functions on the Words and Characters
include("constraints.jl")       # Defines constraints, dominance relations, equivalences, and satisfaction functions
include("dominant_set.jl")      # Defines functions to generate the dominant subset of a constraint set.

include("datastructures.jl")    # Defines the datastructures necessary to improve execution time of the automaton construction.
include("automaton_single.jl")  # Defines the automata generation functions for single constraints
include("automaton_set.jl")     # Defines the automata generation functions for sets of constraints

include("automaton_util.jl")    # Defines the utility functions applied to the automaton
include("sequences.jl")         # Defines functions that generate sequences of a specific automata construction.
include("printing.jl")          # Defines the printing statements for the existing structs

end # module


text:
#!/usr/bin/env bash
    #=
    exec julia -tauto --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

using CSV, DataFrames, StatsPlots, ProgressMeter

include(joinpath(dirname(dirname(@__FILE__)), "src", "CodingTheory.jl"))
using .CodingTheory

stop_at = parse(BigInt, ARGS[1])

# You might want to use your code to search for parameters where q^n divided by sphere_packing_bound(q,n,d) is an integer. I think this is essentially how Golay started his search - it would be interesting to replicate some of that search.
function integer_search(stop_at::Integer)::Array{Array{Number, 1}}
    A = []
    upper_bound = 10
    increment_bound = 10
	processed = []
	i = 0
    p = Progress(stop_at, 1) # minimum update interval: 1 second

    while true
        for q in 1:upper_bound, n in 1:upper_bound, d in 1:upper_bound
			# skip configurations that have already been processed
			# arelessthan(upper_bound - increment_bound, q, n, d) && continue
			(q, n, d) ∈ processed && continue
            # note that we actually don't want to filter out non-linear codes
			# distance shouldn't be larger than the block length; filter trivial distances of one; we filter out combinations when all are equal; filter out codes that are generalised hamming codes; filter out even distances
			# if d < n && ! isone(q) #&& ! CodingTheory.allequal(q, n, d) && ! ishammingperfect(q, n, d) && ! iseven(d) && ! isprimepower(q)
			q, n, d = big(q), big(n), big(d)
				hb = hamming_bound(q, n, d, no_round)
                sb = singleton_bound(q, n, d, no_round)

                # if isinteger(hb) #&& ! isone(hb) && ! isone(sb) && hb ≤ sb
					push!(processed, (q, n, d))
					# i += 1; println("$i:\t$q, $n, $d")
                    push!(A, [q, n, d, hb, sb])
                    isequal(length(A), stop_at) && return A
					next!(p)
                # end
            # end
        end
        upper_bound += increment_bound
    end
end

function make_integer_csv(stop_at::Integer)
    A = integer_search(stop_at)
    D = DataFrame(q = Number[], n = Number[], d = Number[], hamming_bound = Number[], singleton_bound = Number[])

    for i in A
        push!(D, i)
    end

	data_file_desktop = joinpath(homedir(), "Desktop", "bound_integers_$(stop_at).csv")
	data_file_other = joinpath(dirname(dirname(@__FILE__)), "other", "bound_integers_$(stop_at).csv")

    CSV.write(data_file_desktop, D)
    CSV.write(data_file_other, D)
	println("\nWrote data to $(data_file_other).")
end

function bound_comparison(stop_at::Integer)::Tuple{Int, Array{Array{Number, 1}}}
    A = Array{AbstractFloat}[]
	processed = []
    upper_bound = 10
    increment_bound = 10
	i = 0

    while true
        @showprogress for q in 1:upper_bound, n in 1:upper_bound, d in 1:upper_bound
			# skip configurations that have already been processed
			# arelessthan(upper_bound - increment_bound, q, n, d) && continue
			(q, n, d) ∈ processed && continue
			# remove q non prime powers
            if isprimepower(big(q))
                hb = hamming_bound(q, n, d)
                sb = singleton_bound(q, n, d)
				# filter out the more trivial bounds that are one (or if distance is one); filter out distances that are more than or equal to the block length, as they don't make sense in our context; filter out even distances
				# if d < n || 1 ∉ (hb, sb, d)
				if d < n && ! isone(d) && ! CodingTheory.allequal(q, n, d) && ! iszero(mod(d, 2))
					push!(processed, (q, n, d))
					# i += 1; println("$i:\t$q, $n, $d")
                    comparison = hb > sb ? 1 : (sb > hb ? -1 : 0)
                    push!(A, [q, n, d, hb, sb, comparison])
                    isequal(length(A), stop_at) && return stop_at, A
                end
            end
        end
        upper_bound += increment_bound
    end
end

function plot_bound_comparison(stop_at::Integer)
	number_of_bounds, A = bound_comparison(stop_at)
	D = DataFrame(q = Number[], n = Number[], d = Number[], hamming_bound = Number[], singleton_bound = Number[], comparison = Number[])

	for i in A
	    push!(D, i)
	end

	data_file_desktop = joinpath(homedir(), "Desktop", "bound_comparison_$(stop_at).csv")
	data_file_other = joinpath(dirname(dirname(@__FILE__)), "other", "bound_comparison_$(stop_at).csv")

	CSV.write(data_file_desktop, D)
	CSV.write(data_file_other, D)
	println("Data written to $(data_file_desktop) and $(data_file_other)")

	`Rscript --vanilla "~"/projects/CodingTheory.jl/other/regression-tree.R $data_file_other $stop_at`
	println("Regression tree saved at $(joinpath(dirname(dirname(@__FILE__)), "other", "Rplots_$(stop_at).pdf"))")
end

# make_integer_csv constructs a csv file which looks for hamming bounds (given certain conditions) that are integers before rounding
make_integer_csv(stop_at) # takes a command line argument

# plot_bound_comparison constructs a dataframe which compares the hamming bound to the singleton bound and plots a regression tree using r
# plot_bound_comparison(stop_at) # takes a command line argument


text:
#лежит внутри DataGroup и хранит историю операций
mutable struct UndoRedoTool
    history::Vector{String} #храним операции пользователя
    sourse #DataGroup #источник данных
    state::Int64 #в каком состоянии находимся (0 - не было изменений)
    UndoRedoTool() = new([],[],0)
end

#отмена для большой истории правок
function undo!(AllObj::Dict,baseIP::IPv4,port::Union{String,Int64})
    if !haskey(AllObj,"history") || AllObj["state"]==0
        return AllObj
    end
    #определяем, над кем была последняя операция
    manualEvent = JSON.parse(AllObj["history"][AllObj["state"]])
    #понижаем статус
    AllObj["state"] -= 1

    #сопоставляем внешние названия и внутренние данные
    gp_name = manualEvent["chName"]

    #откатываем последнего товарища
    undo!(AllObj["dataStorage"][gp_name],baseIP, port)

    return AllObj
end

#возвращение для большой истории правок
function redo!(AllObj::Dict,baseIP::IPv4,port::Union{String,Int64})
    if !haskey(AllObj,"history") ||  AllObj["state"] == size(AllObj["history"],1)
        return AllObj
    end
    #повышаем статус
    AllObj["state"] += 1
    #определяем, над кем была операция
    manualEvent = JSON.parse(AllObj["history"][AllObj["state"]])
    #сопоставляем внешние названия и внутренние данные
    gp_name = manualEvent["chName"]

    #откатываем последнего товарища
    redo!(AllObj["dataStorage"][gp_name],baseIP, port)

    return AllObj
end

function Base.joinpath(a::Vector{T}) where T
    Path = a[1]
    for i=2:length(a)
        Path = joinpath(Path,a[i])
    end
    Path
end

#перевод информации о действии пользователя из словаря в строку
#данные берутся из JSON
#action - словарь, собранный из JSON с данными о команде
function commandToString(action::Dict,resStr = "")
    for key in keys(action)
        if isa(action[key],Dict)
            #"=>" значит переход на следующий уровень вложенности словаря
            resStr = commandToString(action[key],resStr)
        else
            if resStr == ""
                resStr = key*": "*string(action[key])
            else
                resStr = resStr*"; "*key*": "*string(action[key])
            end
        end
    end
    resStr
end


text:
using Luna

a = 50e-6
b = 10e-6
gas = :Ar
pres = 5
L = 15e-2/10

τfwhm = 30e-15
λ0 = 800e-9
energy = 5e-6

grid = Grid.RealGrid(L, λ0, (160e-9, 3000e-9), 1e-12)

modes = collect(RectModes.RectMode(a, b, gas, pres, :Ag, n=n, m=m) for m in 1:3 for n in 1:6)
nmodes = length(modes)

densityfun = let dens0=PhysData.density(gas, pres)
    z -> dens0
end

ionpot = PhysData.ionisation_potential(gas)
ionrate = Ionisation.ionrate_fun!_ADK(ionpot)

responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
             #Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot))

inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)

Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, modes, :x; full=true)

linop = LinearOps.make_const_linop(grid, modes, λ0)
statsfun = Stats.default(grid, Eω, modes, linop, transform; gas=gas, windows=((150e-9, 300e-9),))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)

Luna.run(Eω, grid, linop, transform, FT, output)

##
Plotting.pygui(true)
Plotting.stats(output)
Plotting.prop_2D(output)
Plotting.prop_2D(output, bandpass=(150e-9, 300e-9))
Plotting.time_1D(output, [5e-2, 9.8e-2])
Plotting.time_1D(output, [5e-2, 9e-2], modes=:sum, bandpass=(150e-9, 300e-9))
Plotting.spec_1D(output, [5e-2, 9.8e-2])
Plotting.spec_1D(output, [5e-2, 9.8e-2], modes=:sum)


text:
"""
Test for https://github.com/JuliaPolyhedra/Polyhedra.jl/issues/301
"""
function issue301test(lib::Polyhedra.Library)
    p = polyhedron(vrep([[3, 0], [0, 3], [0, 0], [1, 1]]), lib)
    removevredundancy!(p; verbose=0, ztol=1e-7)
    @test npoints(p) == 3
    removehredundancy!(p; verbose=0, ztol=1e-7)
    @test nhalfspaces(p) == 3
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
module ImageMorphology

using ImageCore
using ImageCore: GenericGrayImage
using LinearAlgebra
using TiledIteration: EdgeIterator, SplitAxis, SplitAxes
using Requires

include("convexhull.jl")
include("connected.jl")
include("clearborder.jl")
include("dilation_and_erosion.jl")
include("thinning.jl")
include("imfill.jl")
include("maxtree.jl")

include("feature_transform.jl")
using .FeatureTransform

include("deprecations.jl")

export
    dilate,
    dilate!,
    erode,
    erode!,
    opening,
    closing,
    tophat,
    bothat,
    morphogradient,
    morpholaplace,

    # connected.jl
    label_components,
    label_components!,
    component_boxes,
    component_lengths,
    component_indices,
    component_subscripts,
    component_centroids,

    # convexhull.jl
    convexhull,

    # thinning.jl
    thinning,
    GuoAlgo,
    imfill,

    # maxtree.jl
    MaxTree,
    areas, boundingboxes, diameters,
    area_opening, area_opening!, area_closing, area_closing!,
    diameter_opening, diameter_opening!, diameter_closing, diameter_closing!,
    local_maxima!, local_maxima, local_minima!, local_minima,

    #feature_transform.jl
    feature_transform,
    distance_transform,

    clearborder

function __init__()
    @require ImageMetadata = "bc367c6b-8a6b-528e-b4bd-a4b897500b49" begin
        # morphological operations for ImageMeta
        dilate(img::ImageMetadata.ImageMeta, region=coords_spatial(img)) = ImageMetadata.shareproperties(img, dilate!(copy(ImageMetadata.arraydata(img)), region))
        erode(img::ImageMetadata.ImageMeta, region=coords_spatial(img)) = ImageMetadata.shareproperties(img, erode!(copy(ImageMetadata.arraydata(img)), region))
    end
end

end # module


text:
@testset "Vessels" begin
    ## Need to have range in a container for broadcasting to work correctly
    @test all(rand_t .∈ Ref(1:length(Ω)))
    @test all(pref_t .∈ Ref(1:length(Ω)))
    # @test
    # @test all(pref_hist.weights[1] .== pref_hist.weights[end]
    @test q_const[(l = 5, t = 0)] == Catchability(0.2)
    @test q_const[(l = 200, t = 1)] == Catchability(0.2)
    # @test q_const[1, 1] == Catchability(0.2)
    # @test q_const[5, 15] == Catchability(0.2)
    @test q_spat[(l = 500, t = 100)].catchability > 0
    # @test q_sptemp[(l = 400, t = 1)].catchability > 0
    # @test q_vary[35].catchability == 1.0
    @test q_hab[10].catchability == q_hab.base_catchability * hab2[2][10] * 1.2

    @test c1 isa Catch
    @test any(P.P .!= 0.005)

    @test vessels(fleet) == fleet.vessels
    @test c1.catch_biomass ≥ 0

    @test sum(getfield.(c2, :catch_biomass)) > 0
    @test sum(P.P) < 1

    fleet2 = Fleet(fleet.vessels, fleet.total_effort, [1, 2, 3, 4])
    eff2 = FisherySim.order_effort(fleet2)
    @test issorted(eff2)
end


text:
module TestFFIBooster

using LightGBM
using Test
using Random

# we don't want the LightGBM vom
redirect_stderr()

verbosity = "verbose=-1"

@testset "LGBM_BoosterCreate" begin

    mymat = [1. 2.; 3. 4.; 5. 6.]
    dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)

    booster = LightGBM.LGBM_BoosterCreate(dataset, verbosity)

    @test booster.handle != C_NULL

end

@testset "LGBM_BoosterCreateFromModelfile" begin

    booster = LightGBM.LGBM_BoosterCreateFromModelfile(joinpath(@__DIR__, "data", "test_tree"))
    @test booster.handle != C_NULL

end

@testset "LGBM_BoosterLoadModelFromString" begin

    load_str = read(joinpath(@__DIR__, "data", "test_tree"), String)
    booster = LightGBM.LGBM_BoosterLoadModelFromString(load_str)
    @test booster.handle != C_NULL

end

@testset "LGBM_BoosterFree" begin

    mymat = [1. 2.; 3. 4.; 5. 6.]
    dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)

    booster = LightGBM.LGBM_BoosterCreate(dataset, verbosity)

    # These tests exposed a double-free
    @test LightGBM.LGBM_BoosterFree(booster) == nothing
    @test booster.handle == C_NULL

end

@testset "LGBM_BoosterMerge" begin

    booster1 = LightGBM.LGBM_BoosterCreateFromModelfile(joinpath(@__DIR__, "data", "test_tree"))
    booster2 = LightGBM.LGBM_BoosterCreateFromModelfile(joinpath(@__DIR__, "data", "test_tree"))

    @test booster1.handle != booster2.handle

    # record the pointer values so we can verify they dont get mutated as part of the ops
    handle1 = UInt64(booster1.handle)
    handle2 = UInt64(booster2.handle)

    @test LightGBM.LGBM_BoosterMerge(booster1, booster2) == nothing

    @test booster1.handle != C_NULL
    @test booster2.handle != C_NULL

    @test UInt64(booster1.handle) == handle1
    @test UInt64(booster2.handle) == handle2

    @test booster1.handle != booster2.handle

end

@testset "LGBM_BoosterAddValidData" begin

    mymat = [1. 2.; 3. 4.; 5. 6.]
    dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)
    v_dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat .+ 1., verbosity)

    booster = LightGBM.LGBM_BoosterCreate(dataset, verbosity)

    @test LightGBM.LGBM_BoosterAddValidData(booster, v_dataset) == nothing
    @test v_dataset.handle in getfield.(booster.datasets, :handle)
    @test length(booster.datasets) == 2

end

@testset "LGBM_BoosterResetTrainingData" begin

    mymat = [1. 2.; 3. 4.; 5. 6.]
    dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)
    e_dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat .+ 1., verbosity)
    v_dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat .+ 2., verbosity)

    booster = LightGBM.LGBM_BoosterCreate(dataset, verbosity)

    @test LightGBM.LGBM_BoosterResetTrainingData(booster, e_dataset) == nothing
    @test length(booster.datasets) == 1
    @test first(booster.datasets).handle == e_dataset.handle

    # rewrite it but add another dataset
    booster = LightGBM.LGBM_BoosterCreate(dataset, verbosity)
    LightGBM.LGBM_BoosterAddValidData(booster, v_dataset)

    @test LightGBM.LGBM_BoosterResetTrainingData(booster, e_dataset) == nothing
    @test length(booster.datasets) == 2
    @test first(booster.datasets).handle == e_dataset.handle
    @test last(booster.datasets).handle == v_dataset.handle

end

@testset "LGBM_BoosterGetNumClasses" begin

    mymat = [1. 2.; 3. 4.; 5. 6.]
    dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)
    v_dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat .+ 1., verbosity)

    for n in 2:20

        booster = LightGBM.LGBM_BoosterCreate(dataset, "objective=multiclass num_class=$(n) $verbosity")

        @test LightGBM.LGBM_BoosterGetNumClasses(booster) == n

    end

end

@testset "LGBM_BoosterUpdateOneIter" begin

    mymat = [2. 1.; 4. 3.; 5. 6.; 3. 4.; 6. 5.; 2. 1.]
    labels = [1.f0, 0.f0, 1.f0, 1.f0, 0.f0, 0.f0]
    dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)
    LightGBM.LGBM_DatasetSetField(dataset, "label", labels)
    # default params won't allow this to learn anything from this useless data set (i.e. splitting completes)
    booster = LightGBM.LGBM_BoosterCreate(dataset, verbosity)

    finished = LightGBM.LGBM_BoosterUpdateOneIter(booster)
    @test finished == 1
    finished = LightGBM.LGBM_BoosterUpdateOneIter(booster)
    @test finished == 1 # just checking nothing silly happens if you try to continue

    # Feed the tree nuts data so it has a hard time learning
    mymat = randn(1000, 2)
    labels = randn(1000)
    dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)
    LightGBM.LGBM_DatasetSetField(dataset, "label", labels)

    booster = LightGBM.LGBM_BoosterCreate(dataset, verbosity)
    finished = LightGBM.LGBM_BoosterUpdateOneIter(booster)
    @test finished == 0

end

@testset "LGBM_BoosterRollbackOneIter" begin

    # I don't know how to test this, needs thought
    @test_broken false

end

@testset "LGBM_BoosterGetCurrentIteration" begin

    mymat = [2. 1.; 4. 3.; 5. 6.; 3. 4.; 6. 5.; 2. 1.]
    labels = [1.f0, 0.f0, 1.f0, 1.f0, 0.f0, 0.f0]
    dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)
    LightGBM.LGBM_DatasetSetField(dataset, "label", labels)
    # default params won't allow this to learn anything from this useless data set (i.e. splitting completes)
    booster = LightGBM.LGBM_BoosterCreate(dataset, verbosity)

    # check its finished
    @test LightGBM.LGBM_BoosterUpdateOneIter(booster) == 1

    for i in 1:10

        # Iteration number shouldn't increment cause we're not adding any more trees
        LightGBM.LGBM_BoosterUpdateOneIter(booster)
        @test LightGBM.LGBM_BoosterGetCurrentIteration(booster) == 1

    end

    # do again with nuts data
    mymat = randn(1000, 2)
    labels = randn(1000)
    dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)
    LightGBM.LGBM_DatasetSetField(dataset, "label", labels)

    # deliberately constrain the trees to ensure it can't converge by 10 gradient boosts by small chance
    booster = LightGBM.LGBM_BoosterCreate(dataset, "max_depth=1 num_leaf=2 $verbosity")

    for i in 1:10

        LightGBM.LGBM_BoosterUpdateOneIter(booster)
        @test LightGBM.LGBM_BoosterGetCurrentIteration(booster) == i

    end

end

@testset "LGBM_BoosterGetEvalCounts" begin

    # Unsure how to test
    @test_broken false

end

@testset "LGBM_BoosterGetEvalNames" begin

    mymat = [1. 2.; 3. 4.; 5. 6.]
    dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)
    booster = LightGBM.LGBM_BoosterCreate(dataset, "metric=binary_logloss $(verbosity)")
    @test LightGBM.LGBM_BoosterGetEvalNames(booster) |> first == "binary_logloss"

end

@testset "LGBM_BoosterGetFeatureNames" begin

    mymat = [1. 2.; 3. 4.; 5. 6.]
    dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)
    randfeatures = [randstring(2000), randstring(3000)]
    LightGBM.LGBM_DatasetSetFeatureNames(dataset, randfeatures)
    booster = LightGBM.LGBM_BoosterCreate(dataset, verbosity)
    @test LightGBM.LGBM_BoosterGetFeatureNames(booster) == randfeatures

end

@testset "LGBM_BoosterGetNumFeature" begin

    mymat = [1. 2.; 3. 4.; 5. 6.]
    dataset = LightGBM.LGBM_DatasetCreateFromMat(mymat, verbosity)
    booster = LightGBM.LGBM_BoosterCreate(dataset, verbosity)
    @test LightGBM.LGBM_BoosterGetNumFeature(booster) == 2

end

@testset "LGBM_BoosterGetEval" begin

    # Needs implementing
    @test_broken false

end

@testset "LGBM_BoosterGetNumPredict" begin

    # Needs implementing
    @test_broken false

end

@testset "LGBM_BoosterGetPredict" begin

    # Needs implementing
    @test_broken false

end

@testset "LGBM_BoosterCalcNumPredict" begin

    # Needs implementing
    @test_broken false

end

@testset "LGBM_BoosterPredictForMat" begin

    # Needs implementing
    @test_broken false

end

@testset "LGBM_BoosterSaveModel" begin

    booster = LightGBM.LGBM_BoosterCreateFromModelfile(joinpath(@__DIR__, "data", "test_tree"))
    model_path = joinpath(@__DIR__, "data", "test_model.txt")
    if isfile(model_path)
        rm(model_path, force = true)
    end
    LightGBM.LGBM_BoosterSaveModel(booster, 0, 0, 0, model_path )
    @test isfile(model_path)
    rm(model_path; force = true)

end

@testset "LGBM_BoosterSaveModelToString" begin

    booster = LightGBM.LGBM_BoosterCreateFromModelfile(joinpath(@__DIR__, "data", "test_tree"))
    string_repr = LightGBM.LGBM_BoosterSaveModelToString(booster)
    # so it turns out that the string save and file save aren't necesarily the same so..
    # check a bunch of expected substrings, etc
    @test occursin("version=v3", string_repr)
    @test occursin("num_leaves=1", string_repr)
    @test occursin("end of trees", string_repr)
    @test occursin("feature_importances:", string_repr)
    @test occursin("parameters:", string_repr)
    @test occursin("[convert_model: gbdt_prediction.cpp]", string_repr)
    @test occursin("Tree=0", string_repr)
    @test occursin("end of parameters", string_repr)

    # this is an additional test to check the presence and correctness of split and gain after saving model to string
    # due to LGBM v3.0.0 changing the LGBM_BoosterSaveModelToString API that now includes an additional parameter
    # for either split or gain parameters
    # this test is to show that the correct split and gain can be computed regardless of whether 0(split) or 1(gain)
    # `feature_importance_type` parameter is used to save model

    # load a model from file
    booster = LightGBM.LGBM_BoosterCreateFromModelfile(joinpath(@__DIR__, "data", "gain_test_booster"))
    # compute both its gain and split importance
    split_importance = LightGBM.LGBM_BoosterFeatureImportance(booster, 0, 0)
    gain_importance = LightGBM.LGBM_BoosterFeatureImportance(booster, 0, 1)
    # save the model twice, once with each of those gain and split parameters
    split_model_save = LightGBM.LGBM_BoosterSaveModelToString(booster, 0, 0, 0)
    gain_model_save = LightGBM.LGBM_BoosterSaveModelToString(booster, 0, 0, 1)
    # load each of the saved models to a new model
    split_load_model = LightGBM.LGBM_BoosterLoadModelFromString(split_model_save)
    gain_load_model = LightGBM.LGBM_BoosterLoadModelFromString(gain_model_save)
    # for EACH model, compute gain and split importances
    split_importance_from_split_loaded_model = LightGBM.LGBM_BoosterFeatureImportance(split_load_model, 0, 0)
    gain_importance_from_split_loaded_model= LightGBM.LGBM_BoosterFeatureImportance(split_load_model, 0, 1)
    gain_importance_from_gain_loaded_model = LightGBM.LGBM_BoosterFeatureImportance(gain_load_model, 0, 1)
    split_importance_from_gain_loaded_model = LightGBM.LGBM_BoosterFeatureImportance(gain_load_model, 0, 0)
    # and show that both values are the same as the originals for both models
    @test split_importance == split_importance_from_split_loaded_model
    @test split_importance == split_importance_from_gain_loaded_model
    @test gain_importance == gain_importance_from_split_loaded_model
    @test gain_importance == gain_importance_from_gain_loaded_model

end

@testset "LGBM_BoosterFeatureImportance" begin

    booster_path = joinpath(@__DIR__, "data", "gain_test_booster")
    booster = LightGBM.LGBM_BoosterCreateFromModelfile(booster_path)

    split_importance = LightGBM.LGBM_BoosterFeatureImportance(booster, 0, 0)
    gain_importance = LightGBM.LGBM_BoosterFeatureImportance(booster, 0, 1)

    split_sub_importance = LightGBM.LGBM_BoosterFeatureImportance(booster, 2, 0)
    gain_sub_importance = LightGBM.LGBM_BoosterFeatureImportance(booster, 2, 1)

    # this is copy/paste from python output without tweaks to display nums, it should be fine though
    # splits are actually ints
    expected_gain = [89.73980999, 65.49232054, 112.80447054, 107.81817985, 124.81229973]
    expected_split = [17, 11, 19, 20, 23]

    @test isapprox(gain_importance, expected_gain, atol=1e-4)
    @test split_importance == expected_split

    expected_sub_gain = [69.14388967, 52.85783052, 75.43093014, 56.17689991, 96.83372998]
    expected_sub_split = [12,  8, 13, 10, 17]

    @test isapprox(gain_sub_importance, expected_sub_gain, atol=1e-4)
    @test split_sub_importance == expected_sub_split

end

end # module