text:
using QuantumStatistics, Parameters, Printf
using Roots, Polylogarithms
using Plots
# using PackageCompiler
# --sysimage sys_plots.so
# create_sysimage(:Plots, sysimage_path="sys_plots.so", precompile_execution_file="precompile_plots.jl")
include("parameter.jl")
# gr()
# Plots.GRBackend()

@with_kw struct Para
    n::Int = 0 # external Matsubara frequency
    Qsize::Int = 17
    Tsize::Int = 16
    extQ::Vector{SVector{3,Float64}} = [@SVector [q, 0.0, 0.0] for q in LinRange(0.0, 32.0 * kF, Qsize)]
    extT::Vector{SVector{3,Float64}} = [@SVector [T, 0.0, 0.0] for T in LinRange(10, 610, Tsize)]
    fixQ::Float64 = 10* kF
end

function chemical_potential(beta)
    f(β,μ)=real(polylog(3/2,-exp(β*μ)))+4/3/π^0.5*(β)^(3/2)
    g(μ)=f(beta,μ)
    return find_zero(g, (-10000,1))
end

function chi0vsq(beta)
    @unpack n, Qsize, extQ = Para()

    μ=chemical_potential(beta)*EF
    println("μ=$μ, reduced β=$beta")

    pdata = Matrix{Float64}(undef, Qsize, 2)
    for (idx, q) in enumerate(extQ)
        q = q[1]
        p, err = TwoPoint.LindhardΩnFiniteTemperature(dim, q, n, μ, kF, β, me, spin)
        @printf("%10.6f  %10.6f ± %10.6f\n", q / kF, p, err)
        pdata[idx,1] = q/kF
        # pdata[idx,1] = (kF/q)^2
        pdata[idx,2] = -p
        # pdata[idx,1] = (kF/q)^4
        pdata[idx,2] = -p-2/3/π^2*EF^(3/2)/q^2
    end
    return pdata
end

function chi0vsT()
    @unpack n, Tsize, extT, fixQ = Para()
    pdata = Matrix{Float64}(undef, Tsize, 2)
    for (idx, T) in enumerate(extT)
        T = T[1]
        μ=chemical_potential(1.0/T)*EF
        println("μ=$μ, reduced T=$T, fixed extMom q=$fixQ")
        β_in = 1.0/T/kF^2

        p, err = TwoPoint.LindhardΩnFiniteTemperature(dim, fixQ, n, μ, kF, β_in, me, spin)
        @printf("%10.6f  %10.6f ± %10.6f\n", T, p, err)
        # pdata[idx,1] = T
        pdata[idx,1] = 1.0/T
        pdata[idx,2] = -p
    end
    return pdata, fixQ
end
# labels = [string(pdata[i,1]) for i in 1:Tsize]

# @unpack Tsize, extT = Para()
# for (idx, T) in enumerate(extT)
#     T = T[1]
#     pdata = chi0vsq(1.0/T)
#     plot!(pdata[:,1], pdata[:,2], xlabel="q/k_F", ylabel="χ_0-(1/q^2 term)", label="T=$T")
#     # plt = plot(pdata[:,1], pdata[:,2], xlim=(0, 0.1), ylim=(0,) xlabel="(k_F/q)^2", ylabel="χ_0", label="β=$beta", title = "ideal polarization")
#     # plt = plot(pdata[:,1], pdata[:,2], xlabel="(k_F/q)^4", ylabel="χ_0-(1/q^2 term)", label="β=$beta", title = "ideal polarization")
# end
# savefig("chi0vsq_T_v1.pdf")
# show(plt)

pdata, fixQ = chi0vsT()
qint = round(Int, fixQ/kF)
# plt = plot(pdata[:,1], pdata[:,2], xlim=(0,0.2), ylim=(0,0.0017), xlabel="√β", ylabel="χ_0(q=$qint k_F)", title = "ideal polarization", legend=:none)
plt = plot(pdata[:,1], pdata[:,2], xlim=(0,0.04), ylim=(0,0.0017), xlabel="β", ylabel="χ_0(q=$qint k_F)", title = "ideal polarization", legend=:none)
savefig("chi0vsβ_q$qint.pdf")
display(plt)


text:
lb_param_ub = [
0.95    0.99   	1; 	    # beta
0.0	    2      	5;		# gam
0    	0.9   	1;	    # rho1
0       0.02 	0.1;    # sigma1
0    	0.7     1;      # rho2
0	    0.01  	0.1;    # sigma2
6/24    8/24	9/24	# nss
]

ext: jl

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
@testset "Orders" begin

    @testset "blub" begin
        @test 1 + 1 == 2
    end

end


text:
_nnz(Q::AbstractMatrix) = length(Q)
_nnz(Q::SparseMatrixCSC) = nnz(Q)
function write(io::IO,
    A::AbstractMatrix,       # the constraint matrix
    collb::Vector,  # vector of variable lower bounds
    colub::Vector,  # vector of variable upper bounds
    c::Vector,      # vector containing variable objective coefficients
    rowlb::Vector,  # constraint lower bounds
    rowub::Vector,  # constraint upper bounds
    sense::Symbol,           # model sense
    colcat::Vector{Symbol},  # constraint types
    sos::Vector{SOS},        # SOS information
    Q::AbstractMatrix,      #  Quadratic objectives 0.5 * x' Q x
    modelname::AbstractString="MPSWriter_jl",  # MPS model name
    colnames::Vector{String} = ["V$i" for i in 1:length(c)],
    rownames::Vector{String} = ["C$i" for i in 1:length(rowub)]
)
    for (i, cname) in enumerate(colnames)
        if !verifyname(cname)
            colnames[i] = correctname(cname)
        end
    end
    for (i, rname) in enumerate(rownames)
        if !verifyname(rname)
            rownames[i] = correctname(rname)
        end
    end
    if _nnz(Q) > 0
        error("LP writer does not support Quadratic objective")
    end
    if sense != :Max && sense != :Min
        error("sense must be either :Min or :Max. Currently sense =$(sense).")
    end

    if sense == :Max
        println(io,"Maximize")
    else
        println(io,"Minimize")
    end

    print_objective!(io, c, colnames)
    print_constraints!(io, A, rowlb, rowub, colnames, rownames)
    print_sos!(io, sos, colnames)
    print_bounds!(io, collb, colub, colnames)
    print_category!(io, colcat, colnames)

    println(io, "End")
end

function print_objective!(io, c, colnames)
    print(io, "obj: ")
    sp_c = sparse(c)
    is_first = true
    for (col, val) in zip(sp_c.nzind, sp_c.nzval)
        print_variable_coefficient!(io, val, colnames[col], is_first)
        is_first = false
    end
    println(io, "")
end

function print_variable_coefficient!(io, val, name, is_first)
    if is_first
        print_shortest(io, val)
    else
        print(io, val < 0?" - ":" + ")
        print_shortest(io, abs(val))
    end
    print(io, " $(name)")
end

function getrowsense{T1 <: Real, T2<: Real}(rowlb::Vector{T1}, rowub::Vector{T2})
    @assert length(rowlb) == length(rowub)
    row_sense = Array{Symbol}(length(rowub))
    hasranged = false
    for r=1:length(rowlb)
        @assert rowlb[r] <= rowub[r]
    	if (rowlb[r] == -Inf && rowub[r] != Inf) || (rowlb[r] == typemin(eltype(rowlb)) && rowub[r] != typemax(eltype(rowub)))
    		row_sense[r] = :(<=) # LE constraint
    	elseif (rowlb[r] != -Inf && rowub[r] == Inf)  || (rowlb[r] != typemin(eltype(rowlb)) && rowub[r] == typemax(eltype(rowub)))
    		row_sense[r] = :(>=) # GE constraint
    	elseif rowlb[r] == rowub[r]
    		row_sense[r] = :(==) # Eq constraint
        elseif (rowlb[r] == -Inf && rowub[r] == Inf)
            error("Cannot have a constraint with no bounds")
    	else
            row_sense[r] = :ranged
            hasranged = true
    	end
    end
    row_sense, hasranged
end

function print_constraints!(io, A, rowlb, rowub, colnames, rownames)
    println(io, "Subject To")
    row_sense, hasranged = getrowsense(rowlb, rowub)
    if hasranged
        error("LP Writer does not support ranged constraints")
    end

    sA = sparse(A)' # so we can iterate over rows instead of columns

    cols, vals = rowvals(sA), nonzeros(sA)
    for row in 1:length(rowlb)
        is_first = true
        print(io, "$(rownames[row]): ")
        for j in nzrange(sA, row)
            print_variable_coefficient!(io, vals[j], colnames[cols[j]], is_first)
            is_first = false
        end
        if row_sense[row] == :(>=)
            println(io, " $(row_sense[row]) $(rowlb[row])")
        else
            println(io, " $(row_sense[row]) $(rowub[row])")
        end
    end
end

function print_bounds!(io, collb, colub, colnames)
    println(io, "Bounds")
    for (i, (lb, ub, name)) in enumerate(zip(collb, colub, colnames))
        if lb == -Inf && ub == +Inf
            println(io, "$(name) free")
        else
            if lb == -Inf
                print(io, "-inf")
            else
                print_shortest(io, lb)
            end
             print(io, " <= $(name) <= ")
            if ub == Inf
                print(io, "+inf")
            else
                print_shortest(io, ub)
            end
            println(io)
        end
    end
end

function print_category!(io, colcat, colnames)
    if any(colcat .== :SemiCont) || any(colcat .== :SemiInt)
        error("The LP file writer does not currently support semicontinuous or semi-integer variables")
    end
    println(io, "General")# Integer
    for (cat, name) in zip(colcat, colnames)
        if cat == :Int
            println(io, "$(name)")
        end
    end
    println(io, "Binary")
    for (cat, name) in zip(colcat, colnames)
        if cat == :Bin
            println(io, "$(name)")
        end
    end
end

function print_sos!(io, sosconstraints::Vector{SOS}, colnames)
    for (i, sos) in enumerate(sosconstraints)
        print_sos!(io, string("csos", i), sos, colnames)
    end
end
function print_sos!(io, rowname::String, sos::SOS, colnames)
    # format
    # c5: S1:: col_1:1 y:2 z:3
    print(io, rowname, ": ")
    print(io, "S", string(sos[1]), "::")
    for (index, weight) in zip(sos[2], sos[3])
        print(io, " ", colnames[index], ":")
        print_shortest(io, weight)
    end
    println(io)
end


text:
import SoftSurfaceDynamics

msh = SoftSurfaceDynamics.unitsphere(3)

normals = SoftSurfaceDynamics.normals(msh.vertices,msh.faces)
curvatures = SoftSurfaceDynamics.meancurvatures(msh.vertices,msh.faces)
vareas = SoftSurfaceDynamics.vertexareas(msh.vertices,msh.faces)

volume = SoftSurfaceDynamics.surfacevolume(msh.vertices,msh.faces)


text:
export interactive_poincaresos

"""
    interactive_poincaresos(cds, plane, idxs, complete; kwargs...)
Open an interactive application for exploring a Poincaré surface of section (PSOS)
of the continuous dynamical system `cds`. Return an observable containing the
latest initial state created by `complete`, as well as its color.

The `plane` can only be the `Tuple` type accepted by [`poincaresos`](@ref),
i.e. `(i, r)` for the `i`th variable crossing the value `r`. `idxs` gives the two
indices of the variables to be displayed, since the PSOS plot is always a 2D scatterplot.
I.e. `idxs = (1, 2)` will plot the 1st versus 2nd variable of the PSOS. It follows
that `plane[1] ∉ idxs` must be true.

`complete` is a three-argument **function** that completes the new initial state
during interactive use, see below.

The function returns: an observable containing the latest initial `state`
and the `scene` that is plotted. The scatter plot is `scene.children[2]`.

## Keyword Arguments
* `direction, rootkw` : Same use as in [`poincaresos`](@ref).
* `tfinal` : A 2-element tuple for the range of values for the total integration time
  (chosen interactively).
* `Ttr` : A 2-element tuple for the range of values for the transient integration time
  (chosen interactively).
* `markersizes = (-4, -1)` : A 2-element tuple for the range of the marker sizes
  (which scale exponentially: the actual size is `10.0^markersize`).
* `color` : A *function* of the system's initial condition, that returns a color to
  plot the new points with. A random color is chosen by default. Notice that for
  type stability reasons this function must return an instance of
  `RGBf0(red, green, blue)`.
* `labels = ("u₁" , "u₂")` : Axis labels.
* `diffeq...` : Any extra keyword arguments are passed into `init` of DiffEq.

## Interaction
The application is a standard AbstractPlotting scatterplot, which shows the PSOS of the system,
initially using the system's `u0`. Two sliders control the final evolution time and
the size of the marker points.

Upon clicking within the bounds of the scatter plot your click is transformed into
a new initial condition, which is further evolved and its PSOS is computed and then
plotted into the scatter plot.

Your click is transformed into a full `D`-dimensional initial condition through
the function `complete`. The first two arguments of the function are the positions
of the click on the PSOS. The third argument is the value of the variable the PSOS
is defined on. To be more exact, this is how the function is called:
```julia
x, y = mouseclick; z = plane[2]
newstate = complete(x, y, z)
```
The `complete` function can throw an error for ill-conditioned `x, y, z`.
This will be properly handled instead of breaking the application.
This `newstate` is also given to the function `color` that
gets a new color for the new points.
"""
function interactive_poincaresos(ds::ContinuousDynamicalSystem{IIP, S, D}, plane, idxs, complete;
                         # PSOS kwargs:
                         direction = -1,
                         Ttr = (0.0, 1000.0),
                         tfinal = (1000.0, 10.0^4),
                         rootkw = (xrtol = 1e-6, atol = 1e-6),
                         # AbstractPlotting kwargs:
                         color = _randomcolor, resolution = (750, 750),
                         makiekwargs = (),
                         markersizes = (-4, -1),
                         labels = ("u₁" , "u₂"),
                         # DiffEq kwargs:
                         diffeq...) where {IIP, S, D}

    @assert typeof(plane) <: Tuple
    @assert length(idxs) == 2
    @assert eltype(idxs) == Int
    @assert plane[1] ∉ idxs
    u0 = get_state(ds)

    # This is the low-level call of poincaresos:
    ChaosTools._check_plane(plane, D)
    integ = integrator(ds, u0; diffeq...)
    planecrossing = PlaneCrossing(plane, direction > 0)
    i = SVector{2, Int}(idxs)

    # Time sliders:
    ui_tf, tf = AbstractPlotting.textslider(
        range(tfinal[1], tfinal[2], length = 1000), "tfinal", start=tfinal[1]
    )
    ui_ttr, Ttr = AbstractPlotting.textslider(
        range(Ttr[1], Ttr[2], length = 100), "Ttr", start=Ttr[1]
    )

    # Initial Section
    data = poincaresos(integ, planecrossing, tf[], Ttr[], i, rootkw)
    length(data) == 0 && error(ChaosTools.PSOS_ERROR)

    # Plot the first trajectory on the section:
    ui_ms, mss = AbstractPlotting.textslider(range(markersizes[1], markersizes[2];
    length = 100), "size", start = markersizes[2])
    ms = lift(a -> 10.0^a, mss)

    scene = Scene(resolution = (1500, 1000))
    positions_node = Observable(data)
    colors = (c = color(u0); [c for i in 1:length(data)])
    colors_node = Observable(colors)
    scplot = scatter(positions_node, color = colors_node, markersize = ms)
    scplot[Axis][:names][:axisnames] = labels

    laststate = Observable((u0, color(u0)))

    # Interactive clicking on the psos:
    on(events(scplot).mousebuttons) do buttons
        if (ispressed(scplot, Mouse.left) && !ispressed(scplot, Keyboard.space) &&
            AbstractPlotting.is_mouseinside(scplot))

            pos = mouseposition(scplot)

            x, y = pos; z = plane[2] # third variable comes from plane

            newstate = try
               complete(x, y, z)
            catch err
               @error "Could not get state, got error: " exception=err
               return
            end

            reinit!(integ, newstate)

            data = poincaresos(integ, planecrossing, tf[], Ttr[], i, rootkw)

            positions = positions_node[]; colors = colors_node[]
            append!(positions, data)
            c = color(newstate)
            append!(colors, fill(c, length(data)))

            # Notify the signals
            positions_node[] = positions; colors_node[] = colors

            # Update last state
            laststate[] = (newstate, c)

            # AbstractPlotting.scatter!(scplot, data; makiekwargs..., color = color(newstate))
            # display(scene)
        end
    end

    # Button to print current state:
    statebutton = AbstractPlotting.button(Theme(raw = true, camera = campixel!), "latest state")
    on(statebutton[end][:clicks]) do c
        u0, col = laststate[]
        println("Latest state: ")
        println(u0)
        println("with color: $(col)")
    end

    AbstractPlotting.hbox(AbstractPlotting.vbox(ui_ms, ui_tf, ui_ttr, statebutton), scplot, parent=scene)
    display(scene)
    return laststate, scene
end

_randomcolor(args...) = RGBf0(rand(Float32), rand(Float32), rand(Float32))

# TODO: Better estimate of marker size and last Tfinal


text:
let hyperspectrumIndex = 0
    global hyperspectrumCounter() = (hyperspectrumIndex += 1)
end

"""
    HyperSpectrum(arr::Array{T<:Real}, energy::EnergyScale, props::Array{Symbol, Any})

    HyperSpectrum(energy::EnergyScale, props::Dict{Symbol,Any}, arr::Array{<:Real}; axisnames = ( "X", "Y", "Z", "A", "B", "C" ), fov = ( 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 ), livetime=)
    HyperSpectrum(energy::EnergyScale, props::Dict{Symbol,Any}, arr::AxisArray)
    HyperSpectrum(energy::EnergyScale, props::Dict{Symbol,Any}, dims::NTuple{<:Integer}, depth::Int, type::Type{Real}; axisnames = ( "X", "Y", "Z", "A", "B", "C" ), fov = ( 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 )

The HyperSpectrum struct represents a multi-dimensional array of Spectrum objects.  The dimension of a HyperSpectrum
may be 1 for a traverse or a line-scan, 2 for a spectrum image or >2 for time-series of spectrum images or
multi-slice spectrum images.

The first constructor is used to create a HyperSpectrum from a raw Array of data.  The second to construct a HyperSpectrum
from another HyperSpectrum or an AxisArray.  The third from a description of the intended contents.

  * `axisnames`: A list of the names by which the axis can be referred
  * `fov`: The full width of the dimension in mm.

HyperSpectra differ from Array{Spectrum} in that the spectra in a HyperSpectrum must share properties like
:ProbeCurrent and :BeamEnergy.  However, each pixel can have it's own livetime.  HyperSpectrum objects can refer
to line-scans (1D), spectrum images (2D), slice-and-view (3D), time sequenced images (3D), or higher dimension spectrum images.

Internally, HyperSpectrum reinterpretes an Array{T<:Real, N+1} as an Array{Spectrum{T<:Real},N-1}.

HyperSpectrum objects can be read from a RPL/RAW file (using `readrplraw(filenamebase::AbstractString)`) but
can be constructed from any Array{<:Real}.

HyperSpectrum objects can be indexed using the standard Julia array idioms including a single integer index or a
CartesianIndex.  For example, to iterate over every spectrum in a HyperSpectrum

    % Construct a 21 row × 19 column spectrum image with 2048 channels of [0,255].
    hs = HyperSpectrum{es, props, (21,19), 2048, UInt8}
    for idx in eachindex(hs)
        spec = hs[idx]   % get a Spectrum representing the 2048 channels of data at idx
        spec[22] = 1     % Set the 22nd channel to 1
    end
    % Indices into a HyperSpectrum are (row, column)
"""
struct HyperSpectrum{T<:Real, N, NP} <: AbstractArray{Spectrum{T}, N}
    counts::AxisArray{T, NP, Array{T,NP}} # Organized as (ch, r, c) or (ch, y, x)
    energy::EnergyScale
    livetime::Array{Float64, N} # A matrix containing the pixel-by-pixel livetime
    properties::Dict{Symbol,Any}
    stagemap::Type{<:StageMapping} # maps pixel coordinates into stage coordinates
    hash::UInt

    function HyperSpectrum(energy::EnergyScale, props::Dict{Symbol,Any}, arr::Array{<:Real};
        axisnames = ( "Y", "X", "Z", "A", "B", "C" ), #
        fov = [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 ], #
        offset = 0.0 * collect(fov), #
        stagemap::Type{<:StageMapping}=DefaultStageMapping, #
        livetime=fill(get(props, :LiveTime, 1.0), size(arr)[2:end]...)
    )
        axes = Axis[ Axis{:Channel}(1:size(arr,1)) ]
        for i in 1:ndims(arr)-1
            left = offset[i]-fov[i]/2.0
            push!(axes, Axis{Symbol(axisnames[i])}(left:fov[i]/(size(arr,i+1)-1):(left+fov[i])))
        end
        haskey(props, :Name) || (props[:Name] = "HS$(hyperspectrumCounter())")
        counts = AxisArray(arr, axes...)
        new{eltype(arr), ndims(arr) - 1, ndims(arr)}(
            counts,
            energy,
            livetime,
            props,
            stagemap,
            hash(counts, hash(energy, hash(livetime, hash(props, hash(stagemap)))))
        )
    end
    function HyperSpectrum(energy::EnergyScale, props::Dict{Symbol,Any}, arr::AxisArray; stagemap::Type{<:StageMapping}=DefaultStageMapping, livetime=get(props, :LiveTime, 1.0)  * ones(Float64, size(arr[2:end])...))
        haskey(props, :Name) || (props[:Name] = "HS$(hyperspectrumCounter())")
        new{eltype(arr), ndims(arr) - 1, ndims(arr)}(
            arr,
            energy,
            livetime,
            props,
            stagemap,
            hash(arr, hash(energy, hash(livetime, hash(props, hash(stagemap)))))
        )
    end
    function HyperSpectrum(
        energy::EnergyScale,
        props::Dict{Symbol,Any},
        dims::Tuple,
        depth::Int,
        type::Type{<:Real};
        axisnames = ( "Y", "X", "Z", "A", "B", "C" ),
        fov = ( 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 ),
        offset = ( 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ), #
        stagemap::Type{<:StageMapping}=DefaultStageMapping
    )
        axes = Axis[ Axis{:Channel}(1:depth) ]
        for i in 1:length(dims)
            left = offset[i]-fov[i]/2.0
            push!(axes, Axis{Symbol(axisnames[i])}(left:fov[i]/(dims[i]-1):(left+fov[i])))
        end
        haskey(props, :Name) || (props[:Name] = "HS$(hyperspectrumCounter())")
        counts = AxisArray(zeros(type, depth, dims...), axes...)
        livetime = fill(Float64(get(props,:LiveTime,1.0)), dims...)
        new{type, length(dims), length(dims)+1}(
            counts,
            energy,
            livetime,
            props,
            stagemap,
            hash(counts, hash(energy, hash(livetime, hash(props, hash(stagemap)))))
        )
    end
end

Base.hash(hs::HyperSpectrum, h::UInt) = hash(hs.hash, h)

function Base.show(io::IO, ::MIME"text/plain", hss::HyperSpectrum{T,N,NP}) where {T<:Real, N, NP}
    sz = join(repr.(size(hss)), " × ")
    print(
        io,
        "$sz HyperSpectrum{$T,$(axisnames(hss))}[$(name(hss)), $(hss.energy), $(depth(hss)) ch]",
    )
end

function Base.show(io::IO, ::MIME"text/html", hss::HyperSpectrum{T,N,NP}) where {T<:Real, N, NP}
    sz = join(repr.(size(hss)), " × ")
    print(
        io,
        "<p>$sz HyperSpectrum{$T,$(axisnames(hss))}[$(name(hss))), $(hss.energy), $(depth(hss)) ch]</p>",
    )
end

function Base.show(io::IO, hss::HyperSpectrum{T,N,NP}) where {T<:Real, N, NP}
    sz = join(repr.(size(hss)), " × ")
    print(
        io,
        "$sz HyperSpectrum{$T,$(axisnames(hss))}[$(name(hss)), $(hss.energy), $(depth(hss)) ch]",
    )
end

Base.eltype(::HyperSpectrum{T,N,NP}) where {T<:Real,N,NP} = Spectrum{T}
Base.size(hss::HyperSpectrum) = size(hss.counts)[2:end]
Base.length(hss::HyperSpectrum) = prod(size(hss))
Base.ndims(::HyperSpectrum{T,N,NP}) where {T<:Real,N,NP} = N
Base.axes(hss::HyperSpectrum) = Base.axes(counts(hss))[2:end]
Base.firstindex(hss::HyperSpectrum) = 1
Base.lastindex(hss::HyperSpectrum) = length(hss)
Base.eachindex(hss::HyperSpectrum) = Base.OneTo(length(hss))
Base.strides(hss::HyperSpectrum) = strides(counts(hss))[2:end]

# Implement iterators over Spectrum{T}
Base.IteratorSize(::Type{<:HyperSpectrum{T, N, NP}}) where {T <: Real, N, NP} =
    Base.HasShape{NP}()
Base.IndexStyle(::Type{<:HyperSpectrum{T, N, NP}}) where {T <: Real, N, NP} = IndexCartesian()
Base.IteratorEltype(::Type{<:HyperSpectrum{T, N, NP}}) where {T <: Real, N, NP} = Base.HasEltype()
Base.iterate(hss::HyperSpectrum{T, N, NP}, state=1) where {T <: Real, N, NP} =
    state > length(hss) ? nothing : ( getindex(hss, state), state + 1 )

NeXLCore.name(hss::HyperSpectrum) = get(hss.properties, :Name, "HyperSpectrum")
AxisArrays.axisnames(hss::HyperSpectrum) = AxisArrays.axisnames(hss.counts)[2:end]
AxisArrays.axisvalues(hss::HyperSpectrum) = AxisArrays.axisvalues(hss.counts)[2:end]
AxisArrays.axes(hss::HyperSpectrum) = AxisArrays.axes(hss.counts)[2:end]
axisname(hss::HyperSpectrum, ax::Int) = AxisArrays.axisnames(hss.counts)[ax+1]
axisvalue(hss::HyperSpectrum, ax::Int, j::Int) = AxisArrays.axisvalues(hss.counts)[ax+1][j]
axisrange(hss::HyperSpectrum, ax::Int) = AxisArrays.axisvalues(hss.counts)[ax+1]
Base.CartesianIndices(hss::HyperSpectrum{T, N, NP}) where {T<:Real, N, NP} = CartesianIndices(size(hss))
Base.haskey(hss::HyperSpectrum, sym::Symbol) = haskey(hss.properties, sym)

depth(hss::HyperSpectrum) = size(counts(hss), 1)

function matches(spec1::HyperSpectrum, spec2::HyperSpectrum, 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

"""
    coordinate(hss::HyperSpectrum, idx::Tuple{<:Int})

Computes the stage coordinate centering the pixel specified by `idx` using the `StageMapping` `hss.stagemap`.
"""
function coordinate(hss::HyperSpectrum, idx::Tuple)
    av, an = axisvalues(hss.counts), axisnames(hss.counts)
    img_coord = Dict(  an[i+1] => av[i+1][ii] for (i,ii) in enumerate(idx) )
    stage_coord = get(hss.properties, :StagePosition, Dict(sy=>zero(typeof(av[2][1])) for sy in (:X, :Y) ))
    th = get(hss.properties, :ImageRotation, 0.0)
    return image2stage(hss.stagemap, stage_coord, img_coord, th)
end
coordinate(hss::HyperSpectrum, ci::CartesianIndex) = #
    coordinate(hss, ci.I)

"""
    dose(hss::HyperSpectrum) # Average dose per pixel
    dose(hss::HyperSpectrum, idx...) # Dose for the `idx` pixel

Returns the product of the probe current and the live-time on a per pixel basis.
"""
dose(hss::HyperSpectrum) = hss[:ProbeCurrent] * mean(hss.livetime)
dose(hss::HyperSpectrum, idx...)  = hss[:ProbeCurrent] * hss.livetime[idx...]

"""
    livetime!(hss::HyperSpectrum, lt::AbstractFloat, idx...)
    livetime!(hss::HyperSpectrum{T,N}, lt::AbstractFloat) # All pixels to lt

Set the livetime on a per pixel basis.
"""
function livetime!(hss::HyperSpectrum, lt::AbstractFloat, idx...)
    hss.livetime[idx...] = lt
end
function livetime!(hss::HyperSpectrum, lt::AbstractFloat)
    hss.livetime .= lt
end

livetime(hss::HyperSpectrum, idx...) = hss.livetime[idx...]

function Base.getindex(
    hss::HyperSpectrum{T,N,NP},
    ci::CartesianIndex
)::Spectrum{T} where { T<: Real, N, NP}
    props = copy(hss.properties)
    props[:Cartesian] = ci
    props[:Name] = "$(name(hss))[$(join(["$i" for i in ci.I],","))]"
    props[:StageCoordinate]=coordinate(hss, ci)
    props[:LiveTime] = hss.livetime[ci]
    return Spectrum(hss.energy, counts(hss)[:, ci], props)
end
function Base.getindex(
    hss::HyperSpectrum{T, N, NP},
    idx::Vararg{Int,N},
)::Spectrum{T} where {T<:Real, N, NP}
    getindex(hss, CartesianIndex(idx...))
end
function Base.getindex(hss::HyperSpectrum{T, N, NP}, idx::Int)::Spectrum{T} where {T<:Real, N, NP}
    getindex(hss, CartesianIndices(hss)[idx])
end
function Base.getindex(
    hss::HyperSpectrum{T, N, NP},
    vidx::Vector
)::Vector{Spectrum{T}} where { T<: Real, N, NP}
    return [ hss[i] for i in vidx ]
end
function Base.getindex(
    hss::HyperSpectrum{T, N, NP},
    rng::AbstractRange{<:Integer},
)::Vector{Spectrum{T}} where { T<: Real, N, NP}
    return [ hss[i] for i in rng ]
end
function Base.getindex(
    hss::HyperSpectrum{T, N, NP},
    idx...
)::HyperSpectrum{T, N, NP} where {T<:Real, N, NP}
    props = copy(hss.properties)
    props[:Name] = "$(name(hss))[$idx]"
    fov = map(1:length(idx)) do i
        ax = AxisArrays.axes(hss,i)
        ax[last(idx[i])]-ax[first(idx[i])]
    end
    offset = map(1:length(idx)) do i
        ax = AxisArrays.axes(hss,i)
        0.5*(ax[last(idx[i])]+ax[first(idx[i])])
    end
    return HyperSpectrum(hss.energy, props, counts(hss)[:, idx...], livetime=hss.livetime[idx...], fov = fov, offset = offset)
end

function Base.getindex(
    hss::HyperSpectrum,
    cxr::CharXRay
)
    return roiimage(hss, cxr)
end
function Base.getindex(
    hss::HyperSpectrum,
    cxrs::Vector{CharXRay}
)
    return roiimages(hss, cxrs)
end
function Base.getindex(
    hss::HyperSpectrum,
    cxr1::CharXRay,
    cxr2::CharXRay
)
    colorize(hss, [cxr1, cxr2], :All)
end
function Base.getindex(
    hss::HyperSpectrum,
    cxr1::CharXRay,
    cxr2::CharXRay,
    cxr3::CharXRay
)
    colorize(hss, [cxr1, cxr2, cxr3], :Each)
end

function Base.get(hss::HyperSpectrum, sym::Symbol, def)
    get(hss.properties, sym, def)
end

"""
    region(hss::HyperSpectrum{T, N}, ranges::AbstractRange...)::HyperSpectrum where {T<:Real,N}

Creates a view of a HyperSpectrum to represent the range of pixels within the `hss` HyperSpectrum.  Does
not copy the data or properties so any modifications to the region are also made to `hss`.
"""
function region(
    hss::HyperSpectrum{T, N, NP},
    ranges::AbstractRange...,
)::HyperSpectrum{T, N, NP} where {T<:Real, N, NP}
    HyperSpectrum(hss.energy, hss.properties, counts(hss)[:, ranges...], livetime=hss.livetime[ranges...])
end

function Base.getindex(hss::HyperSpectrum, sy::Symbol)
    @assert sy!=:LiveTime "Use the livetime(...) function to get the :LiveTime property for hyperspectra."
    getindex(hss.properties, sy)
end
function Base.setindex!(hss::HyperSpectrum, val::Any, sy::Symbol)
    @assert sy!=:LiveTime "Use the livetime!(...) function to set the :LiveTime property for hyperspectra."
    isequal(sy, :Detector) && @assert isequal(hss.energy, val.scale) "The energy calibration in $val must match the HyperSpectrum's: $(hss.energy)."
    setindex!(hss.properties, val, sy)
end

NeXLCore.energy(ch::Integer, hss::HyperSpectrum) = energy(ch, hss.energy)
channel(energy::Real, hss::HyperSpectrum) = channel(energy, hss.energy)
channelwidth(ch::Int, hss::HyperSpectrum) = energy(ch + 1, hss) - energy(ch, hss)
rangeofenergies(ch::Integer, hss::HyperSpectrum) =
    (energy(ch, hss.energy), energy(ch + 1, hss.energy))
properties(hss::HyperSpectrum)::Dict{Symbol,Any} = hss.properties

"""
	matching(spec::HyperSpectrum, resMnKa::Float64, lld::Int=1)::BasicEDS

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

matching(
    spec::HyperSpectrum,
    resMnKa::Float64,
    lld::Int,
    minByFam::Dict{Shell,Element},
)::BasicEDS = BasicEDS(depth(spec), spec.energy, MnKaResolution(resMnKa), lld, minByFam)

"""
    counts(hss::HyperSpectrum{T, N, NP})::Array{T,NP}

Creates type-friendly view of the counts data array.  Use of this function helps to avoid performance
penalties associated with boxing/unboxing the counts data.

    counts(hss::HyperSpectrum{T, N, NP}, ci::CartesianIndex)::Vector{T}

Access the counts data associated with the pixel `ci`.

    counts(hss::HyperSpectrum{T, N, NP}, ci::CartesianIndex, ch::Int)::T

Access the counts data at the pixel represented by `ci` and the channel represented by `ch`.
"""
counts(hss::HyperSpectrum{T, N, NP}) where {T<:Real, N, NP} = parent(hss.counts)
counts(hss::HyperSpectrum{T, N, NP}, ci::CartesianIndex, ch::Int) where {T<:Real, N, NP} = counts(hss)[ch, ci]
counts(hss::HyperSpectrum{T, N, NP}, ci::CartesianIndex) where {T<:Real, N, NP} = view(counts(hss)[:, ci])

"""
    compress(hss::HyperSpectrum)

Returns a HyperSpectrum with smaller or equal storage space to `hss` without losing or truncating any counts
(note: AbstractFloat compresses to Float32 with loss of precision).  Can change the storage type and/or
reduce the depth of hss.  If there is nothing that can be done to reduce the size of `hss` then `compress(hss)===hss`.

Example:

     hs = compress(hs) # Replace `hs` with a version that uses less memory (if possible).
"""
function compress(hss::HyperSpectrum{T, N, NP}) where {T<:Integer, N, NP}
    (minval, maxval) = extrema(counts(hss))
    tyopts = minval < 0 ? (Int8, Int16, Int32, Int64) : (UInt8, UInt16, UInt32, UInt64)
    restype = tyopts[findfirst(tyopts) do ty
        minval >= typemin(ty) && maxval <= typemax(ty)
    end]
    @assert sizeof(restype) <= sizeof(T) "$restype"
    return compact(restype, hss)
end
function compress(hss::HyperSpectrum{T, N, NP}) where {T<:AbstractFloat, N, NP}
    (minval, maxval) = extrema(counts(hss))
    restype=T
    if sizeof(T) < sizeof(Float32) &&
        minval >= typemin(Float32) &&
        maxval <= typemax(Float32)
        restype=Float32
    end
    return compact(restype, hss)
end

"""
    compact(S::Type{<:Real}, hss::HyperSpectrum)

Converts the counts data type to S.  Also eliminates any channel planes on the end that
are entirely less than or equal to zero.
"""
function compact(S::Type{<:Real}, hss::HyperSpectrum)
    data = counts(hss)
    itr = Iterators.reverse(Base.axes(data,1))
    maxi = findfirst(itr) do ch # First non-negative channel
        !all((<=)(zero(eltype(data))), @view data[ch, Base.axes(data)[2:end]...])
    end
    @assert !isnothing(maxi) "There are no non-negative counts in this hyperspectrum."
    maxch = itr[maxi]
    if S != eltype(data) || maxch < size(data,1)
        sdata = maxch == size(data,1) ? S.(data) : S.(data[1:maxch, Base.axes(data)[2:end]...])
        return HyperSpectrum(hss.energy, copy(hss.properties), sdata, livetime=hss.livetime)
    else
        return hss
    end
end

"""
   plane(hss::HyperSpectrum, chs::AbstractUnitRange{<:Integer}, normalize=false)

Sums a contiguous range of data channels into an Array. The dimension of the result is
one less than the dimension of the HyperSpectrum and is stored as a Float64 to ensure that not information is lost.
"""
function plane(
    hss::HyperSpectrum{T, N, NP},
    chs::AbstractUnitRange{<:Integer},
    normalize = false,
) where {T<:Real, N, NP}
    data = counts(hss)
    res = map(CartesianIndices(hss)) do ci
        sum(view(data, chs, ci), init=zero(Float64))
    end
    return normalize ? res ./= maximum(res) : res
end
function plane(
    hss::HyperSpectrum{T, N, NP},
    cxr::CharXRay,
    normalize = false,
) where {T<:Real, N, NP}
    plane(hss, fwhmroi(hss, cxr), normalize)
end

"""
   plane(hss::HyperSpectrum, ch::Int, normalize=false)

Extracts a single channel plane from a HyperSpectrum. The dimension of the result is
one less than the dimension of the HyperSpectrum.
"""
function plane(hss::HyperSpectrum, ch::Int, normalize = false)
    res = counts(hss)[ch, CartesianIndices(hss)]
    return normalize ? res / maximum(res) : res
end

"""
    maxpixel(hss::HyperSpectrum, filt=ci->true)
    maxpixel(hss::HyperSpectrum, cis::CartesianIndices, filt=ci->true)
    maxpixel(hss::HyperSpectrum, mask::BitArray)
    minpixel(hss::HyperSpectrum, filt=ci->true)
    minpixel(hss::HyperSpectrum, cis::CartesianIndices, filt=ci->true)
    minpixel(hss::HyperSpectrum, mask::BitArray)

Compute Bright's Max-Pixel derived signal for the entire HyperSpectrum or a rectanglar sub-region.
"""
maxpixel(hss::HyperSpectrum, mask::BitArray) =
    maxpixel(hss, CartesianIndices(hss), ci -> mask[ci])
function maxpixel(
    hss::HyperSpectrum{T, N, NP},
    cis::CartesianIndices,
    filt::Function = ci -> true,
)::Spectrum{T} where {T<:Real, N, NP}
    data, res = counts(hss), zeros(T, depth(hss))
    for ci in filter(ci -> filt(ci), cis), ch in eachindex(res)
        res[ch] = max(res[ch], data[ch, ci])
    end
    props = copy(hss.properties)
    props[:LiveTime] = mean(hss.livetime)
    props[:Name] = "MaxPixel[$(props[:Name])]"
    return Spectrum(hss.energy, res, props)
end
function maxpixel(
    hss::HyperSpectrum{T, N, NP},
)::Spectrum{T} where {T<:Real, N, NP}
    data, res  = counts(hss), fill(zero(T), depth(hss))
    for ci in CartesianIndices(hss), ch in eachindex(res)
        res[ch] = max(res[ch], data[ch, ci])
    end
    props = copy(hss.properties)
    props[:LiveTime] = mean(hss.livetime)
    props[:Name] = "MaxPixel[$(props[:Name])]"
    return Spectrum(hss.energy, res, props)
end

function minpixel(
    hss::HyperSpectrum{T, N, NP},
)::Spectrum{T} where {T<:Real, N, NP}
    data, res  = counts(hss), fill(typemax(T), depth(hss))
    for ci in CartesianIndices(hss), ch in eachindex(res)
        res[ch] = min(res[ch], data[ch, ci])
    end
    props = copy(hss.properties)
    props[:LiveTime] = mean(hss.livetime)
    props[:Name] = "MaxPixel[$(props[:Name])]"
    return Spectrum(hss.energy, res, props)
end
minpixel(hss::HyperSpectrum, mask::BitArray) =
    minpixel(hss, CartesianIndices(hss), ci -> mask[ci])
function minpixel(
    hss::HyperSpectrum{T, N, NP},
    cis::CartesianIndices,
    filt::Function = ci -> true,
) where {T<:Real, N, NP}
    data, res  = counts(hss), fill(typemax(T), depth(hss))
    for ci in filter(ci -> filt(ci), cis), ch in eachindex(res)
        res[ch] = min(res[ch], data[ch, ci])
    end
    props = copy(hss.properties)
    props[:LiveTime] = mean(hss.livetime)
    props[:Name] = "MaxPixel[$(props[:Name])]"
    return Spectrum(hss.energy, res, props)
end
avgpixel(hss::HyperSpectrum, mask::BitArray) =
    avgpixel(hss, CartesianIndices(hss), ci -> mask[ci])
function avgpixel(
    hss::HyperSpectrum{T, N, NP},
    cis::CartesianIndices,
    filt::Function = ci -> true,
) where {T<:Real, N, NP}
    data = counts(hss)
    res = zeros(T isa Int ? Int64 : Float64, depth(hss))
    for ci in filter(ci -> filt(ci), cis)
        res .+= @view data[:, ci]
    end
    res /= prod(size(hss))
    lt = mean(hss.livetime)
    props = copy(hss.properties)
    props[:Name] = something(name, "Average[$(props[:Name])]")
    return Spectrum(hss.energy, res, props)
end
function avgpixel(
    hss::HyperSpectrum{T, N, NP}
) where {T<:Real, N, NP}
    data = counts(hss)
    res = zeros(T isa Int ? Int64 : Float64, depth(hss))
    for ci in CartesianIndices(hss)
        res .+= @view data[:, ci]
    end
    res /= prod(size(hss))
    lt = mean(hss.livetime)
    props = copy(hss.properties)
    props[:Name] = something(name, "Average[$(props[:Name])]")
    return Spectrum(hss.energy, res, props)
end

"""
    sumcounts(hss::HyperSpectrum, cis::CartesianIndices = CartesianIndices(hss))

An array containing the number of counts at each pixel.
"""
function sumcounts(hss::HyperSpectrum, cis::CartesianIndices = CartesianIndices(hss))
    data = counts(hss)
    map(cis) do ci
        sum(view(data, :, ci), init=zero(Float64))
    end
end

"""
    indexofmaxpixel(hss::HyperSpectrum, ch::Int) # at channel `ch`
    indexofmaxpixel(hss::HyperSpectrum) # all channels
    indexofmaxpixel(hss::HyperSpectrum, ch::Int, cis::CartesianIndices)
    indexofmaxpixel(hss::HyperSpectrum, cis::CartesianIndices)

Find the indices producing the maximum value in data[ch] or data[:] within 'cis' or full
spatial dimensions.
"""
function indexofmaxpixel(hss::HyperSpectrum, ch::Int, cis::CartesianIndices)
    data = counts(hss)
    maxidx, max = cis[1], data[ch, cis[1]]
    for idx in cis
        if data[ch, idx] > max
            maxidx = idx
            max = data[ch, idx]
        end
    end
    return maxidx
end

indexofmaxpixel(hss::HyperSpectrum, ch::Int) =
    indexofmaxpixel(hss, ch, CartesianIndices(hss))

indexofmaxpixel(hss::HyperSpectrum) = indexofmaxpixel(hss, CartesianIndices(hss))

function indexofmaxpixel(hss::HyperSpectrum{T, N, NP}, cis::CartesianIndices) where {T<:Real, N, NP}
    res, cix = zeros(T, depth(hss)), CartesianIndex[CartesianIndices(hss)[1] for _ = 1:depth(hss)]
    data = counts(hss)
    for ci in cis, i in filter(i -> data[i, ci] > res[i], eachindex(res))
        res[i] = data[i, ci]
        cix[i] = ci
    end
    return cix
end

"""
    Base.sum(hss::HyperSpectrum{T, N}) where {T<:Real, N}
    Base.sum(hss::HyperSpectrum{T, N}, mask::Union{BitArray, Array{Bool}}) where {T<:Real, N}
    Base.sum(hss::HyperSpectrum{T, N}, func::Function) where {T<:Real, N}
    Base.sum(hss::HyperSpectrum{T, N}, cis::CartesianIndices) where {T<:Real, N}

Compute a sum spectrum for all or a subset of the pixels in `hss`.

Where func(hss::HyperSpectrum, ci::CartesianIndex)::Bool
"""
function Base.sum(
    hss::HyperSpectrum{T, N, NP},
    mask::Union{BitArray,Array{Bool}};
    name = nothing
) where {T<:Real, N, NP}
    @assert size(mask) == size(hss) "Mask size[$(size(mask))] ≠ Hyperspectrum size[$(size(hss))]"
    data = counts(hss)
    res, lt = zeros(T isa Int ? Int64 : Float64, depth(hss)), 0.0
    for ci in filter(ci->mask[ci], CartesianIndices(hss))
        res .+= data[:, ci]
        lt += hss.livetime[ci]
    end
    props = copy(hss.properties)
    props[:Name] = something(name, "MaskedSum[$(props[:Name])]")
    return Spectrum(hss.energy, res, props)
end
function Base.sum(
    hss::HyperSpectrum{T, N, NP},
    filt::Function;
    name = nothing
) where {T<:Real, N, NP}
    data = counts(hss)
    res, lt = zeros(T isa Int ? Int64 : Float64, depth(hss)), 0.0
    for ci in filter(ci->filt(hss,ci), CartesianIndices(hss))
        res .+= data[:, ci]
        lt += hss.livetime[ci]
    end
    props = copy(hss.properties)
    props[:Name] = something(name, "FilteredSum[$(props[:Name])]")
    return Spectrum(hss.energy, res, props)
end
function Base.sum(hss::HyperSpectrum{T, N, NP}; name = nothing) where {T<:Real, N, NP}
    data = counts(hss)
    res = zeros(T isa Int ? Int64 : Float64, depth(hss))
    for ci in CartesianIndices(hss)
        res .+= @view data[:, ci]
    end
    lt = sum(hss.livetime)
    props = copy(hss.properties)
    props[:Name] = something(name, "Sum[$(props[:Name])]")
    return Spectrum(hss.energy, res, props)
end
function Base.sum(hss::HyperSpectrum{T, N, NP}, cis::CartesianIndices; name = nothing) where {T<:Real, N, NP}
    data = counts(hss)
    res, lt = zeros(T isa Int ? Int64 : Float64, depth(hss)), 0.0
    for ci in cis
        res .+= data[:, ci]
        lt += hss.livetime[ci]
    end
    props = copy(hss.properties)
    props[:Name] = something(name, "IndexedSum[$(props[:Name])]")
    return Spectrum(hss.energy, res, props)
end

using Images

"""
    roiimages(hss::HyperSpectrum, achs::AbstractVector{<:AbstractUnitRange{<:Integer}})

Create an array of Gray images representing the intensity in each range of channels in
in `achs`.  They are normalized such the the most intense pixel in any of them defines white.
(Accounts for differences in :LiveTime between pixels.)
"""
function roiimages(hss::HyperSpectrum, achs::AbstractVector{<:AbstractUnitRange{<:Integer}})
    data = counts(hss)
    ps = map(achs) do chs
        map(CartesianIndices(hss)) do ci
            sum(view(data, chs, ci), init=zero(Float64)) / hss.livetime[ci]
        end
    end
    maxval = maximum(map(p -> maximum(p), ps))
    return map(p -> Gray.(convert.(N0f8, p / maxval)), ps)
end

"""
    roiimage(hss::HyperSpectrum, chs::AbstractUnitRange{<:Integer})

Create a count map from the specified contiguous range of channels.
(Accounts for differences in :LiveTime between pixels.)
"""
function roiimage(hss::HyperSpectrum, chs::AbstractUnitRange{<:Integer})
    data = counts(hss)
    tmp = map(CartesianIndices(hss)) do ci
        sum(view(data, chs, ci), init=zero(Float64)) / hss.livetime[ci]
    end
    Gray.(N0f8.(tmp ./= maximum(tmp)))
end

function fwhmroi(hss::Union{Spectrum, HyperSpectrum}, cxr::CharXRay, k=1.0)
    ecxr = energy(cxr)
    r = haskey(hss.properties, :Detector) ? resolution(ecxr, hss[:Detector]) : resolution(ecxr, MnKaResolution(130.0))
    return channel(ecxr-0.5*r*k, hss.energy):channel(ecxr+0.5*r*k, hss.energy)
end

"""
    roiimage(hss::HyperSpectrum, cxr::CharXRay)

Create a count map for the specified characteristic X-ray.  By default, integrates
for one FWHM at `cxr`.  If hss[:Detector] is an `EDSDetector`, the FWHM is taken from it.
Otherwise, a FWHM of 130 eV at Mn Kα is assumed.
"""
roiimage(hss::HyperSpectrum, cxr::CharXRay) = roiimage(hss, fwhmroi(hss, cxr))

"""
    roiimages(hss::HyperSpectrum, cxrs::AbstractVector{CharXRay})

Create an array of Gray images representing the intensity in each of the CharXRay lines
in `cxrs`.  They are normalized such the the most intense pixel in any of them defines white.
By default, integrates for one FWHM at `cxr`.  If hss[:Detector] is an `EDSDetector`, the FWHM
is taken from it.   Otherwise, a FWHM of 130 eV at Mn Kα is assumed.
"""
roiimages(hss::HyperSpectrum, cxrs::AbstractVector{CharXRay}) = roiimages(hss, map(cxr->fwhmroi(hss, cxr), cxrs))

"""
    colorize(hss::HyperSpectrum, cxrs::AbstractVector{CharXRay}, normalize=:All)

Create RGB colorized images from up to three `CharXRay` which define channels over which
the count signal is integrated. `normalize=:All` puts the intensities on a common scale
using the `roiimages(...)` method.  Otherwise each image is scaled independently based on
the brightest pixel.
"""
function NeXLCore.colorize(hss::HyperSpectrum, cxrs::AbstractVector{CharXRay}, normalize=:All)
    if normalize==:All
        imgs = roiimages(hss, cxrs)
        colorview(RGB, imgs[1], length(imgs) > 1 ? imgs[2] : zeroarray, length(imgs) > 2 ? imgs[3] : zeroarray )
    else
        # Normalize relative to max of each ROIs independently
        colorview(RGB, roiimage(hss,cxrs[1]), length(cxrs)>1 ? roiimage(hss,cxrs[2]) : zeroarray, length(cxrs) > 2 ? roiimage(hss, cxrs[3]) : zeroarray)
    end
end

"""
    block(hss::HyperSpectrum{T,N,NP}, steps::NTuple{N, Int})::HyperSpectrum{T,N,NP} where {T<:Real, N, NP}
    block(hss::HyperSpectrum{T,N,NP}, step::Int) where {T<:Real, N, NP}

Reduce the size of a HyperSpectrum by summing together blocks of adjacent pixels.  For example, `steps=(4,4)` would
sum together blocks of 16 spectra in `hss` to form a single pixel in the resulting `Hyperspectrum`.
"""
function block(hss::HyperSpectrum{T,N,NP}, step::Int) where {T<:Real, N, NP}
    block(hss, ntuple(_ -> step, N))
end
function block(hss::HyperSpectrum{T,N,NP}, steps::NTuple{N, Int})::HyperSpectrum where {T<:Real, N, NP}
    dims = size(hss) .÷ steps
    fov = [last(ax.val)-first(ax.val) for ax in AxisArrays.axes(hss.counts)[2:end]]
    offset = [first(ax.val) for ax in AxisArrays.axes(hss.counts)[2:end]]
    et = eltype(hss.counts) isa Integer ? widen(eltype(hss.counts)) : eltype(hss.counts)
    res = HyperSpectrum(hss.energy, copy(hss.properties), dims, depth(hss), et,
        axisnames=axisnames(hss.counts)[2:end], fov=fov, offset=offset, stagemap=hss.stagemap)
    cx, rcx = counts(hss), counts(res) # Force the type to eliminate boxing
    ii = ntuple(_ -> 1, length(steps))
    tmp = zeros(et, depth(hss)) # accumulate it here
    for ci in CartesianIndices(dims)
        tmp .= zero(et)
        lt = zero(eltype(res.livetime))
        for inner in CartesianIndices(steps)
            idx = (ci.I .- ii) .* steps .+ inner.I
            tmp .+= @view cx[:, idx...] # @view is important!!!
            lt += hss.livetime[idx...]
        end
        rcx[:, ci] .= tmp
        res.livetime[ci] = lt
    end
    res
end

"""
    linescan(hss::HyperSpectrum{T,2,3}, ci1::CartesianIndex{2}, ci2::CartesianIndex{2}, width::Int=1)

Extract pixels from `hss` along the line from `ci1` to `ci2` as a 1D HyperSpectrum.  The `width` argument integrates
the linescan along a line perpendicular to the primary axis. Only works on 2D SpectrumImages and for odd values of `width`.
The algorithm does double count the occasional pixel but the length of each perpendicular is maintained at `width`.
The `AxisArrays.axes(...)` scaling is maintained so lengths on the linescan can be compared to lengths on the original map.
"""
function linescan(hss::HyperSpectrum{T,2,3}, ci1::CartesianIndex{2}, ci2::CartesianIndex{2}, width::Int=1) where {T<:Real}
    @assert width >= 1 "Only positive widths are supported in linescan."
    @assert ci1 in CartesianIndices(hss) "$ci1 is not inside $(CartesianIndices(hss))"
    @assert ci2 in CartesianIndices(hss) "$ci2 is not inside $(CartesianIndices(hss))"
    function perpray(c1, c2, w) # Computes a list of CI's starting at c1 towards c2 of length w skipping c1
        i, res = 0, CartesianIndex[]
        drawray(c1, c2) do ci
            if (i>0) && (i<=w)
                push!(res, CartesianIndex(ci))
            end
            (i+=1)<=w
        end
        return res
    end
    width % 2 == 0 && @info "Even widths are rounded down to the nearest odd width."
    fov = [ sqrt(sum(map(1:2) do ii
        a = AxisArrays.axes(hss, ii)
        (a[ci2.I[ii]]-a[ci1.I[ii]])^2
    end)) ]
    len = sum(abs.(ci2.I .- ci1.I)) + 1
    hs2 = HyperSpectrum(hss.energy, copy(hss.properties), (len, ), depth(hss), T,
        axisnames=[ "Linescan[$(ci1.I), $(ci2.I), $width]" ], fov = fov, offset = 0.0*fov)
    cxh, cxr, idx = counts(hss), counts(hs2), 0
    drawline(ci1, ci2, true) do ci
        idx+=1
        cxr[:, idx] .= @view cxh[:, ci...]
        hs2.livetime[idx] = hss.livetime[ci...]
    end
    if (width-1) ÷ 2 > 0
        bounds = CartesianIndices(hss)
        for rot in ( [ 0 -1; 1 0 ], [ 0 1; -1 0 ])
            perp = rot * collect(ci2.I .- ci1.I)
            r1 = perpray(ci1, CartesianIndex((perp  .+ ci1.I)...), (width-1) ÷ 2)
            r2 = perpray(ci2, CartesianIndex((perp  .+ ci2.I)...), (width-1) ÷ 2)
            for (ci1r, ci2r) in zip(r1, r2)
                idx = 0
                drawline(ci1r, ci2r, true) do ci
                    idx+=1
                    if all(i->ci[i] in bounds.indices[i], 1:2)
                        cxr[:, idx] .+= @view cxh[:, ci...]
                        hs2.livetime[idx] += hss.livetime[ci...]
                    end
                end
            end
        end
    end
    return hs2
end

function gridize(img,step=16)
    for ci in CartesianIndices(img)
        if any(i->i%step==0, ci.I)
            img[ci] = 1-img[ci]
        end
    end
    return img
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
module BuildExecutable
export build_executable
@windows_only using WinRPM
# Builds an executable that doesn't require any julia source code.
# The user needs to provide a julia script that contains a function main(),
# taking no argument, which will be run when executing the
# produced executable.

# Note on packages:
# Even if the script contains using statements, exported functions
# will not be available in main(). Full qualification of names is
# required. It's suggested to replace using statements with import
# statements to produce a consistent result between running main() in
# the REPL and running the executable.

type Executable
    name
    filename
    buildfile
    targetfile
    libjulia
end

function Executable(exename, targetdir, debug)
    if debug
        exename = exename * "-debug"
    end
    filename = exename
    @windows_only filename = filename * ".exe"
    buildfile = abspath(joinpath(JULIA_HOME, filename))
    targetfile = targetdir == nothing ? buildfile : joinpath(targetdir, filename)
    libjulia = debug ? "-ljulia-debug" : "-ljulia"

    Executable(exename, filename, buildfile, targetfile, libjulia)
end

type SysFile
    buildpath
    buildfile
	inference
    inference0
end

function SysFile(exename, debug=false)
    buildpath = abspath(dirname(Libdl.dlpath(debug ? "libjulia-debug" : "libjulia")))
    buildfile = joinpath(buildpath, "lib"*exename)
    inference = joinpath(buildpath, "inference")
    inference0 = joinpath(buildpath, "inference0")
    SysFile(buildpath, buildfile, inference, inference0)
end

function build_executable(exename, script_file, targetdir=nothing, cpu_target="native";
                          force=false, debug=false)
    julia = abspath(joinpath(JULIA_HOME, debug ? "julia-debug" : "julia"))
    if !isfile(julia * @windows? ".exe" : "")
        println("ERROR: file '$(julia)' not found.")
        return 1
    end
    build_sysimg = abspath(dirname(@__FILE__), "build_sysimg.jl")
	if !isfile(build_sysimg)
		build_sysimg = abspath(JULIA_HOME, "..", "..", "contrib", "build_sysimg.jl")
		if !isfile(build_sysimg)
			println("ERROR: build_sysimg.jl not found.")
			return 1
		end
	end

    if targetdir != nothing
        patchelf = find_patchelf()
        if patchelf == nothing && !(OS_NAME == :Windows)
            println("ERROR: Using the 'targetdir' option requires the 'patchelf' utility. Please install it.")
            return 1
        end
    end

    if !isfile(script_file)
        println("ERROR: $(script_file) not found.")
        return 1
    end

    tmpdir = mktempdir()
    cfile = joinpath(tmpdir, "start_func.c")
    userimgjl = joinpath(tmpdir, "userimg.jl")
    script_file = abspath(script_file)

    if targetdir != nothing
        targetdir = abspath(targetdir)
        if !isdir(targetdir)
            println("ERROR: targetdir is not a directory.")
            return 1
        end
    end

    exe_file = Executable(exename, targetdir, debug)
    sys = SysFile(exename, debug)

    if !force
        for f in [cfile, userimgjl, "$(sys.buildfile).$(Libdl.dlext)", "$(sys.buildfile).ji", exe_file.buildfile]
            if isfile(f)
                println("ERROR: File '$(f)' already exists. Delete it or use --force.")
                return 1
            end
        end

        if targetdir != nothing && !isempty(readdir(targetdir))
            println("ERROR: targetdir is not an empty diectory. Delete all contained files or use --force.")
            return 1
        end
    end

    emit_cmain(cfile, exename, targetdir != nothing)
    emit_userimgjl(userimgjl, script_file)

    empty_cmd_str = ``
    println("running: $(julia) $(build_sysimg) $(sys.buildfile) $(cpu_target) $(userimgjl) --force" * (debug ? " --debug" : ""))
    run(`$(julia) $(build_sysimg) $(sys.buildfile) $(cpu_target) $(userimgjl) --force $(debug ? "--debug" : empty_cmd_str)`)
    println()

    gcc = find_system_gcc()
    win_arg = ``
    # This argument is needed for the gcc, see issue #9973
    @windows_only win_arg = Base.WORD_SIZE==32 ? `-D_WIN32_WINNT=0x0502 -march=pentium4` : `-D_WIN32_WINNT=0x0502`
    incs = get_includes()
    ENV2 = deepcopy(ENV)
    @windows_only begin
        if contains(gcc, "WinRPM")
            # This should not bee necessary, it is done due to WinRPM's gcc's
            # include paths is not correct see WinRPM.jl issue #38
            ENV2["PATH"] *= ";" * dirname(gcc)
            push!(incs, "-I"*abspath(joinpath(dirname(gcc),"..","include")))
        end
    end

    println("running: $gcc -g $win_arg $(join(incs, " ")) $(cfile) -o $(exe_file.buildfile) -Wl,-rpath,$(sys.buildpath) -L$(sys.buildpath) $(exe_file.libjulia) -l$(exename)")
    cmd = setenv(`$gcc -g $win_arg $(incs) $(cfile) -o $(exe_file.buildfile) -Wl,-rpath,$(sys.buildpath) -Wl,-rpath,$(sys.buildpath*"/julia") -L$(sys.buildpath) $(exe_file.libjulia) -l$(exename)`, ENV2)
    run(cmd)
    println()

    println("running: rm -rf $(tmpdir) $(sys.buildfile).o $(sys.inference).o $(sys.inference).ji $(sys.inference0).o $(sys.inference0).ji")
    map(f-> rm(f, recursive=true), [tmpdir, sys.buildfile*".o", sys.inference*".o", sys.inference*".ji", sys.inference0*".o", sys.inference0*".ji"])
    println()

    if targetdir != nothing
        # Move created files to target directory
        for file in [exe_file.buildfile, sys.buildfile * ".$(Libdl.dlext)", sys.buildfile * ".ji"]
            mv(file, joinpath(targetdir, basename(file)), remove_destination=force)
        end

        # Copy needed shared libraries to the target directory
        tmp = ".*\.$(Libdl.dlext).*"
        paths = [sys.buildpath]
        VERSION>v"0.5.0-dev+5537" && is_unix() && push!(paths, sys.buildpath*"/julia")
        for path in paths
            shlibs = filter(Regex(tmp),readdir(path))
            for shlib in shlibs
                cp(joinpath(path, shlib), joinpath(targetdir, shlib), remove_destination=force)
            end
        end

        @unix_only begin
            # Fix rpath in executable and shared libraries
            # old implementation for fixing rpath in shared libraries
            #=
            shlibs = filter(Regex(tmp),readdir(targetdir))
            push!(shlibs, exe_file.filename)
            for shlib in shlibs
                rpath = readall(`$(patchelf) --print-rpath $(joinpath(targetdir, shlib))`)[1:end-1]
                # For debug purpose
                #println("shlib=$shlib\nrpath=$rpath")
                if Base.samefile(rpath, sys.buildpath)
                    run(`$(patchelf) --set-rpath $(targetdir) $(joinpath(targetdir, shlib))`)
                end
            end
            =#
            # New implementation
            shlib = exe_file.filename
            @linux_only begin
                run(`$(patchelf) --set-rpath \$ORIGIN/ $(joinpath(targetdir, shlib))`)
            end
            @osx_only begin
                # For debug purpose
                #println(readall(`otool -L $(joinpath(targetdir, shlib))`)[1:end-1])
                #println("sys.buildfile=",sys.buildfile)
                run(`$(patchelf) -rpath $(sys.buildpath) @executable_path/ $(joinpath(targetdir, shlib))`)
                run(`$(patchelf) -change $(sys.buildfile).$(Libdl.dlext) @executable_path/$(basename(sys.buildfile)).$(Libdl.dlext) $(joinpath(targetdir, shlib))`)
                #println(readall(`otool -L $(joinpath(targetdir, shlib))`)[1:end-1])
            end
        end
    end

    println("$(exe_file.targetfile) successfully created.")
    return 0
end

function find_patchelf()
    @linux_only for patchelf in [joinpath(JULIA_HOME, "patchelf"), "patchelf"]
        try
            if success(`$(patchelf) --version`)
                return patchelf
            end
        end
    end
    @osx_only "install_name_tool"
end

function get_includes()
    ret = []

    # binary install
    incpath = abspath(joinpath(JULIA_HOME, "..", "include", "julia"))
    push!(ret, "-I$(incpath)")

    # Git checkout
    julia_root = abspath(joinpath(JULIA_HOME, "..", ".."))
    push!(ret, "-I$(julia_root)src")
    push!(ret, "-I$(julia_root)src/support")
    push!(ret, "-I$(julia_root)usr/include")

    ret
end

function emit_cmain(cfile, exename, relocation)
    if relocation
        sysji = joinpath("lib"*exename)
    else
        sysji = joinpath(dirname(Libdl.dlpath("libjulia")), "lib"*exename)
    end
    sysji = escape_string(sysji)
    if VERSION > v"0.5.0-dev+4397"
        arr = "jl_alloc_vec_any"
        str = "jl_string_type"
    else
        arr = "jl_alloc_cell_1d"
        str = "jl_utf8_string_type"
    end
    f = open(cfile, "w")
    write( f, """
        #include <julia.h>
        #include <stdlib.h>
        #include <stdio.h>
        #include <assert.h>
        #include <string.h>
        #if defined(_WIN32) || defined(_WIN64)
        #include <malloc.h>
        #endif

        void failed_warning(void) {
            if (jl_base_module == NULL) { // image not loaded!
                char *julia_home = getenv("JULIA_HOME");
                if (julia_home) {
                    fprintf(stderr,
                            "\\nJulia init failed, "
                            "a possible reason is you set an envrionment variable named 'JULIA_HOME', "
                            "please unset it and retry.\\n");
                }
            }
        }

        int main(int argc, char *argv[])
        {
            char sysji[] = "$(sysji).ji";
            char *sysji_env = getenv("JULIA_SYSIMAGE");
            char mainfunc[] = "main()";

            assert(atexit(&failed_warning) == 0);

            jl_init_with_image(NULL, sysji_env == NULL ? sysji : sysji_env);

            // set Base.ARGS, not Core.ARGS
            if (jl_base_module != NULL) {
                jl_array_t *args = (jl_array_t*)jl_get_global(jl_base_module, jl_symbol("ARGS"));
                if (args == NULL) {
                    args = $arr(0);
                    jl_set_const(jl_base_module, jl_symbol("ARGS"), (jl_value_t*)args);
                }
                assert(jl_array_len(args) == 0);
                jl_array_grow_end(args, argc - 1);
                int i;
                for (i=1; i < argc; i++) {
                    jl_value_t *s = (jl_value_t*)jl_cstr_to_string(argv[i]);
                    jl_set_typeof(s,$str);
                    jl_arrayset(args, s, i - 1);
                }
            }

            // call main
            jl_eval_string(mainfunc);

            int ret = 0;
            if (jl_exception_occurred())
            {
                jl_show(jl_stderr_obj(), jl_exception_occurred());
                jl_printf(jl_stderr_stream(), "\\n");
                ret = 1;
            }

            jl_atexit_hook(ret);
            return ret;
        }
        """
    )
    close(f)
end

function emit_userimgjl(userimgjl, script_file)
    open(userimgjl, "w") do f
        write( f, "include(\"$(escape_string(script_file))\")")
    end
end

function find_system_gcc()
    # On Windows, check to see if WinRPM is installed, and if so, see if gcc is installed
    @windows_only try
        winrpmgcc = joinpath(WinRPM.installdir,"usr","$(Sys.ARCH)-w64-mingw32",
            "sys-root","mingw","bin","gcc.exe")
        if success(`$winrpmgcc --version`)
            return winrpmgcc
        end
    end

    # See if `gcc` exists
    @unix_only try
        if success(`gcc -v`)
            return "gcc"
        end
    end

    error( "GCC not found on system: " * @windows? "GCC can be installed via `Pkg.add(\"WinRPM\"); WinRPM.install(\"gcc\")`" : "" )
end

end # module


text:
# Try converting a string into a number
function tryNumberify(tagValue::String)
    outValue = tagValue
    try outValue = float(tagValue)
    catch
        splitTagValue = split(tagValue, " ")
        try outValue = [float(i) for i in splitTagValue]
        end
    end
    return outValue
end

# Parse the header files (i.e. d3proc, reco, visu_pars)
function parseHeaderFile(pathToFile::String)
    file_IO = open(pathToFile)
    hdrText = readstring(file_IO)
    close(file_IO)

    hdrElements = replace(hdrText, "\n", "")
    hdrElements = split(hdrElements, "##\$")
    deleteat!(hdrElements,1)

    tagNames=[replace(i, r"(.*?)=(.*)", s"\1") for i in hdrElements] # Get tag names/fields
    tagValues=[replace(i, r"(.*?)=(.*)", s"\2") for i in hdrElements] # Get values for each tag name/field
    tagValues=[replace(i, r"\([^\)]+\)\s?", "") for i in tagValues] # Remove parenthesis
    tagValues=[replace(i, r"\$.*", "") for i in tagValues] # Remove the '$$ @vis' endings
    tagValues=[replace(i, r"[<>]", "") for i in tagValues] # Remove '<' and '>'
    myDict = Dict(zip(tagNames,[tryNumberify(i) for i in tagValues]))
    return myDict
end

# MAIN FUNCTION - loads image data from paravision's d3seq
function load2dseq(pathTo2dseq::String)
    baseDir = dirname(pathTo2dseq)

    filesToLookFor = ["d3proc", "reco", "visu_pars"]
    hdrDict = Dict(zip(filesToLookFor,[Dict() for i in filesToLookFor]))
    for file in filesToLookFor
        if isfile(joinpath(baseDir,file))
           hdrDict[file] = parseHeaderFile(joinpath(baseDir,file))
        end
    end

    if length(hdrDict["visu_pars"]) > 0
        visuPars = hdrDict["visu_pars"]
        coreSize = visuPars["VisuCoreSize"]
        coreFrames = visuPars["VisuCoreFrameCount"]
        byteType = visuPars["VisuCoreWordType"]
        slopes = visuPars["VisuCoreDataSlope"]
        offsets = visuPars["VisuCoreDataOffs"]
    elseif length(hdrDict["reco"]) > 0
        visuPars = hdrDict["reco"]
        coreSize = visuPars["RECO_size"]
        slopes = 1 ./ visuPars["RECO_map_slope"] # RECOslope is 1/visuparsSlope
        offsets = visuPars["RECO_map_offset"]
        byteType = visuPars["RECO_wordtype"]
        # Infer number of frames by number of slopes
        coreFrames = length(slopes)
    else
        print("Need visu_pars or reco file to read data")
    end

    totalSize = [coreSize; coreFrames]
    imgData = read(pathTo2dseq, byteDict[byteType], tuple(Int.(totalSize)...))

    # Using a shortcut that assumes that all frames have same slope, so verify that first
    @assert minimum(slopes) == maximum(slopes)
    imgData = imgData .* slopes[1] .+ offsets[1]

    return(imgData, hdrDict)
end

byteDict = Dict("_8BIT_UNSGN_INT" => Int8, "_16BIT_SGN_INT" => Int16, "_32BIT_SGN_INT" => Int32, "_32BIT_FLOAT" => Float32)


text:
include("../src/Thirian/fractional_shift.jl")

# include("../src/FractionalDelayFilters.jl")
# import .FractionalDelayFilters

using LinearAlgebra
using FFTW
import PyPlot

import Utilities # https://gitlab.com/RoyCCWang/utilities

PyPlot.close("all")
fig_num = 1

PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

# sine for debug.
f = sin
P = LinRange(-pi, pi, 3000)
U = P
f_U = f.(U)

## filter.
s = real.(f_U)

N_thirian = 15
b_thirian = Vector{Float64}(undef, N_thirian)

τ = 0.31
applydelaythirian!(b_thirian, s, τ)

## visualize.

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(P, real.(f_U), label = "f")
PyPlot.plot(P, s, label = "s", "--")

PyPlot.legend()
PyPlot.xlabel("ppm")
PyPlot.ylabel("real")
PyPlot.title("reversible fractional delay filtering")

x_ranges = [ P ]
s_itp, d_s_itp, d2_s_itp = Utilities.setupcubicitp(real.(f_U), x_ranges, 1.0)
s_func = xx->s_itp([xx])

S = FiniteSignalType(real.(f_U), P[1], P[end])

# no shifts.
j = 14 # test this time index.
println("no shift: j ", j)
s0 = real.(f_U)
println("zero: evalsignal( S, j ) - s0[j] = ", evalsignal( S, j ) - s0[j])
println("zero: s0[j] - s_func(P[j] = ", s0[j] - s_func(P[j]))
println()

# integer shifts.
d = 19
applyintegershift!(S, d)

j = d+4 # test this time index.
println("shift: j ", j)
evalsignal( S, j ) # evaluates s0[j-d].
println("zero: s0[j-d] - evalsignal( S, j ) = ", s0[j-d] - evalsignal( S, j ))
println()

# fractional shift, τ ∈ [0, 0.5].
ΔP = P[2]-P[1]

d = 2 # this is the integer delay (positive).
applyintegershift!(S, d)

τ = 0.319 # fractional delay (positive).

resetsignal!(S)
applydelaythirian!(b_thirian, S.s, τ)
evals_S = collect( evalsignal(S, j) for j = 1:length(S.s) )

evals_s_func = collect( s_func( P[j]-ΔP*τ ) for j = 1:length(P) )

println("norm difference between cubic itp and fractional delay: ", norm(evals_s_func - evals_S))
#[evals_S evals_s_func]

# visualize.
PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(P, real.(f_U), label = "f")
PyPlot.plot(P, evals_S, label = "evals_S", "--")
PyPlot.plot(P, evals_s_func, label = "evals_s_func", "--")

PyPlot.legend()
PyPlot.xlabel("ppm")
PyPlot.ylabel("real")
PyPlot.title("cubic itp vs. fractional delay filtering")

### test reversible property of thirian delay filter. see https://lcondat.github.io/publis/condat_icip08_rotations.pdf
τ = 0.21
resetsignal!(S)
applydelaythirian!(b_thirian, S.s, τ)

tmp_s = copy(S.s)

S.s[:] = reverse(S.s)
applydelaythirian!(b_thirian, S.s, τ)
S.s[:] = reverse(S.s)

#reverse!(tmp_s)
applyadvancethirian!(b_thirian, tmp_s, τ) # bake reverse!() into this function.
#reverse!(tmp_s)
println("test applyadvancethirian!: ", norm(tmp_s-S.s))

## timing.
# import BenchmarkTools
# BenchmarkTools.@btime applydelaythirian!(b_thirian, tmp_s, τ)
# BenchmarkTools.@btime applyadvancethirian!(b_thirian, tmp_s, τ)

evals_S = collect( evalsignal(S, j) for j = 1:length(S.s) )
println("reversible test: ", norm(evals_S - s0))
println()

### arb shift, test reversible, frac part ∈ [0, 0.5]
resetsignal!(S)

shift = 1.23
applyshift!(S, shift, b_thirian)

shift = -1.23
applyshift!(S, shift, b_thirian)
evals_S = collect( evalsignal(S, j) for j = 1:length(S.s) )
println("reversible test: ", norm(evals_S - s0))
println()

### arb shift, test reversible, frac part ∈ [0.5, 1]
resetsignal!(S)

shift = 1.9
applyshift!(S, shift, b_thirian)

shift = -1.9
applyshift!(S, shift, b_thirian)
evals_S = collect( evalsignal(S, j) for j = 1:length(S.s) )
println("reversible test: ", norm(evals_S - s0))
println()

######## get shift coordinates from Δp, which is in units of P.
Δp = 2.0
a = P[1]
b = P[end]

resetsignal!(S)

shift = Δp*length(S.s)/(b-a)
applyshift!(S, shift, b_thirian)

# visualize.
evals_S = collect( evalsignal(S, j) for j = 1:length(S.s) )
evals_s_func = collect( s_func( P[j]-ΔP*τ ) for j = 1:length(P) )

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(P, real.(f_U), label = "f")
PyPlot.plot(P, evals_S, label = "evals_S", "--")
#PyPlot.plot(P, evals_s_func, label = "evals_s_func", "--")

PyPlot.legend()
PyPlot.xlabel("ppm")
PyPlot.ylabel("real")
PyPlot.title("fractional shift with Δp = $(Δp)")

max_stars_repo_name: RoyCCWang/FractionalDelayFilters

max_issues_repo_name: RoyCCWang/FractionalDelayFilters

max_forks_repo_name: RoyCCWang/FractionalDelayFilters


text:
@testset "$TEST $G" begin

# construct graph
g = DG(3)
add_edge!(g,1,2)
add_edge!(g,2,3)

# source and sink terminals
source, target = 1, 3

# default capacity
capacity_matrix = Erdos.DefaultCapacity(g)

# state variables
flow_matrix = zeros(typeof(signed(V(0))), 3, 3)

TREE = zeros(Int, 3)
TREE[source] = 1
TREE[target] = 2

PARENT = zeros(Int, 3)

A = [source,target]

residual_graph = complete(g)

path = Erdos.find_path!(residual_graph, source, target, flow_matrix, capacity_matrix, PARENT, TREE, A)

@test path == [1,2,3]

end # testset


text:
using CUDA
using Oceananigans

using ConvergenceTests.ForcedFlowFixedSlip: setup_and_run_xy, setup_and_run_xz

arch = CUDA.has_cuda() ? GPU() : CPU()

# Run 4 simulations:
Nx = [16, 32, 64, 128]
stop_time = 0.01
h = π / maximum(Nx)

Δt = 0.01 .* h^2 # Initial guess

# Modify so simulations stop on a dime
stop_iteration = round(Int, stop_time / Δt)
Δt = stop_time / stop_iteration

# Run (x, y) and (x, z) simulations
for N in Nx
    setup_and_run_xy(architecture=arch, Nx=N, Δt=Δt, stop_iteration=stop_iteration)
    setup_and_run_xz(architecture=arch, Nx=N, Δt=Δt, stop_iteration=stop_iteration)
end


text:
using Test
using ExtensibleEffects
using DataTypesBasic
using TypeClasses
using Suppressor
splitln(str) = split(strip(str), "\n")
import ExtensibleEffects: autorun

# Callable
# --------

# callable is the only standard type which needs a custom handler
# hence a perfect example of how such custom handlers are lovely integrated into the API

myflatmap = @syntax_flatmap begin
  a = Callable(x -> x+2)
  b = Callable(x -> a + x)
  @pure a, b
end
@test myflatmap(3) == (5, 8)

myeff = Callable(function(args...; kwargs...)
  @runhandlers CallableHandler(args...; kwargs...) @syntax_eff begin
    a = Callable(x -> x+2)
    b = Callable(x -> a + x)
    @pure a, b
  end
end)
@test myeff(3) == (5, 8)

# we provide a custom macro for this common callable pattern which does exactly the same

myeff2 = @runcallable @syntax_eff begin
  a = Callable(x -> x+2)
  b = Callable(x -> a + x)
  @pure a, b
end
@test myeff2(3) == (5, 8)

myeff3 = @runcallable @syntax_eff begin
  v = [1,3,4]
  a = Callable(x -> x+v)
  o = isodd(a) ? Option(100) : Option()
  b = Callable(x -> x + a + o)
  @pure v, a, o, b
end

@test myeff3(1) == [Const(nothing), Const(nothing), Identity((4,5,100,106))]
@test myeff3(2) == [Identity((1,3,100,105)), Identity((3,5,100,107)), Const(nothing)]

max_stars_repo_name: schlichtanders/ExtensibleEffects

max_issues_repo_name: schlichtanders/ExtensibleEffects

max_forks_repo_name: schlichtanders/ExtensibleEffects


text:
"""
    dare(A, B, Q, R)

Compute `X`, the solution to the discrete-time algebraic Riccati equation,
defined as A'XA - X - (A'XB)(B'XB + R)^-1(B'XA) + Q = 0, where Q>=0.

This version relaxes the requirement that R>0, rather
R = [R1 0;0 R2] where R1 >0, R2 < 0. This is a usuful formulation for
H_∞ synthesis.

Algorithm taken from:
Laub, "A Schur Method for Solving Algebraic Riccati Equations."
Relaxes the positive definiteness of R, useful for robust control.

http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf

Implementation stolen from: ControlSystems.jl
"""
function dare(A, B, Q, R)
    if !issemiposdef(Q)
        error("Q must be positive-semidefinite.");
    end
    n = size(A, 1);

    E = [
        Matrix{Float64}(I, n, n) B*(R\B');
        zeros(size(A)) A'
    ];
    F = [
        A zeros(size(A));
        -Q Matrix{Float64}(I, n, n)
    ];

    QZ = schur(F, E);
    QZ = ordschur(QZ, abs.(QZ.alpha./QZ.beta) .< 1);

    return QZ.Z[(n+1):end, 1:n]/QZ.Z[1:n, 1:n];
end

issemiposdef(A) = ishermitian(A) && minimum(real.(eigvals(A))) >= 0
issemiposdef(A::UniformScaling) = real(A.λ) >= 0
issquare(A) = size(A)[1] == size(A)[2]


text:
module TypeSynonyms

using OffsetArrays

const Float = Float64
export Float

# OffsetArray{ComplexF64}
const ComplexOffsetMatrix = OffsetArray{Complex{Float},2, Array{Complex{Float}, 2}}
export ComplexOffsetMatrix

# OffsetVector{ComplexF64}
const ComplexOffsetVector = OffsetArray{Complex{Float},1,Array{Complex{Float}, 1}}
export ComplexOffsetVector

# OffsetVector{Int}
const IntOffsetVector = OffsetArray{Int, 1, Array{Int, 1}}
export IntOffsetVector

const RationalComplex = Complex{Rational{Int}}
export RationalComplex

const RComplex2IntDict = Dict{RationalComplex, Int}
export RComplex2IntDict

const RComplex2OffsetMatrix = Dict{RationalComplex,ComplexOffsetMatrix}
export RComplex2OffsetMatrix

const Variable = Int
export Variable

const Coefficient = ComplexF64
export Coefficient

const Pole = RationalComplex
export Pole

const Power = Int
export Power

end # module TypeSynonims


text:
#include("SparseDeconvolution.jl")
#results = SparseDeconvolution.run_demo()
#SparseDeconvolution.show_results(results...)
#
#results = SparseDeconvolution.run_demo_JuMP()
#SparseDeconvolution.show_results(results...)

#include("LineSpectraEstimation.jl")
#results = LineSpectraEstimation.run_demo()
###results = LineSpectraEstimation.run_demo_Convex()
#LineSpectraEstimation.show_results(results...)

#include("MatrixDecomposition.jl")
#results = MatrixDecomposition.run_demo()
#MatrixDecomposition.show_results(results...)

include("DNN.jl")
results = DNN.run_demo()
DNN.show_results(results...)

#include("TotalVariation.jl")
#results = TotalVariation.run_demo()
#TotalVariation.show_results(results...)

#include("AudioDeclipping.jl")
#results = AudioDeclipping.run_demo()
#AudioDeclipping.show_results(results...)
#AudioDeclipping.save_wav(results...)


text:
"""
    ewma_canonical(series, λ) -> Float64

Canonical definition of Exponentially Weighted Moving Average.
"""
function ewma_canonical(series::Vector{Float64}, λ::Float64) :: Float64
    n = length(series)

    local sumprod::Float64 = 0.0
    local sumweights::Float64 = 0.0

    for (i, serie_element) in enumerate(series)
        weight = λ^(n-i)
        sumprod += serie_element * weight
        sumweights += weight
    end

    return sumprod / sumweights
end

"""
    ewma(returns, λ) -> Float64

Same as `ewma_canonical`, but squares returns and applies sqrt on the result.
Use this function to estimate the volatility of a return series with the EWMA model.
"""
function ewma(returns::Vector{Float64}, λ::Float64) :: Float64
    n = length(returns)

    local sumprod::Float64 = 0.0
    local sumweights::Float64 = 0.0

    for (i, return_element) in enumerate(returns)
        weight = λ^(n-i)
        sumprod += (return_element^2) * weight
        sumweights += weight
    end

    return sqrt(sumprod / sumweights)
end

"""
    ewma_cov(returns::Array{Float64, 2}, λ::Float64) -> Array{Float64, 2}

In the case of a matrix of returns, considers that each column is a time series of returns.
The result is the covariance matrix.
"""
function ewma_cov(returns::Array{Float64, 2}, λ::Float64) :: Array{Float64, 2}
    (rows, cols) = size(returns)
    covmatrix = Array{Float64, 2}(undef, cols, cols)

    @inbounds for j in 1:cols, i in 1:j
        local sumprod::Float64 = 0.0
        local sumweights::Float64 = 0.0

        for r in 1:rows
            weight = λ^(rows-r)
            sumprod += weight * returns[r, i] * returns[r, j]
            sumweights += weight
        end

        covmatrix[i,j] = sumprod / sumweights

        if i != j
            covmatrix[j,i] = covmatrix[i,j]
        end
    end

    return covmatrix
end

"""
    ewma(returns::Array{Float64, 2}, λ::Float64) -> Array{Float64, 2}

In the case of a matrix of returns, considers that each column is a time series of returns.
The result is the correlation matrix.
"""
function ewma(returns::Array{Float64, 2}, λ::Float64) :: Array{Float64, 2}
    (rows, cols) = size(returns)
    covmatrix = ewma_cov(returns, λ)
    @assert size(covmatrix) == (cols, cols)

    # vols is the Volatility vector
    # We can get them from the sqrt of the diagonal of the cov matrix
    vols = sqrt.(LinearAlgebra.diag(covmatrix))

    # turns covmatrix into a correlation matrix
    @inbounds for j in 1:cols, i in 1:j
        covmatrix[i,j] = covmatrix[i,j] / (vols[i] * vols[j])

        if i != j
            covmatrix[j,i] = covmatrix[i,j]
        end
    end

    return covmatrix
end


text:
using Documenter, ReferenceTests

format = Documenter.HTML(edit_link = "master",
                         prettyurls = get(ENV, "CI", nothing) == "true")

makedocs(
    format = format,
    modules = [ReferenceTests],
    sitename = "ReferenceTests.jl",
    authors = "Christof Stocker",
    linkcheck = !("skiplinks" in ARGS),
    pages = Any[
        "Home" => "index.md",
    ]
)

deploydocs(repo = "github.com/johnnychen94/DemoCards.jl.git")


text:
using IHT
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
include("MendelIHT_test.jl")
include("MendelIHT_utilities_test.jl")

# To see code coverage, cd to the IHT folder and run the following
# julia -e 'Pkg.test("IHT",coverage=true)'
# @show get_summary(process_file("src/MendelIHT.jl"))
# @show get_summary(process_file("src/MendelIHT_utilities_test.jl"))


text:
using DistCtrl4DistMan: runExp
using Plots

DEP_params = Dict{String,Any}("maxDist" => (300e-6, 550e-6));
MAG_params = Dict{String,Any}("maxDist" => (50e-3, 75e-3));
ACU_params = Dict{String,Any}("maxDist" => (45e-3, 65e-3));

## DEP - two microparticles, 16x16 array of electrodes
runExp(platform=:DEP, N_iter=25, N_agnts=2, N_acts=(16, 16), convanalysis=false, saveplots=true, figfilename="DEParray_small.pdf", params=DEP_params);

## DEP - six microparticles, 16x16 array of electrodes
runExp(platform=:DEP, N_iter=25, N_agnts=6, N_acts=(16, 16),convanalysis=false, saveplots=true, figfilename="DEParray.pdf", params=DEP_params);

## MAG - five steel balls, 8x8 array of coils
runExp(platform=:MAG, N_iter=25, N_agnts=5, N_acts=(8, 8), convanalysis=false, saveplots=false, figfilename="MAGarray_8x8.pdf", plotActuatorCommands=true);

## MAG - five steel balls, 12x12 array of coils
runExp(platform=:MAG, N_iter=25, N_agnts=5, N_acts=(12, 12), convanalysis=false, saveplots=true, figfilename="MAGarray_12x12.pdf", plotActuatorCommands=true);

## ACU - six pressure points, 16x16 array of actuators
runExp(platform=:ACU, N_iter=25, N_agnts=6, N_acts=(16, 16),convanalysis=false, saveplots=true, figfilename="ACUarray.pdf");

## DEP - six microparticles, 16x16 array of electrodes, Algorithm 1 vs Algorithm 2
N_agnts=2;
N_iter=25;
N_acts = 16;

# run the experiment with Algorithm 2 (i.e. with different c_i and d_i sets)
params, exp_data_A2 = runExp(platform=:DEP, N_iter=N_iter, N_agnts=N_agnts, N_acts=(N_acts,N_acts),
                            convanalysis=true,
                            saveplots=false,
                            showplots=false,
                            params=DEP_params);

# Extract the randomly generated agents' positions and desdired DEP forces
agnts = exp_data_A2[1]["Agents"];
agent_positions = [agnts[i].pos for i in 1:length(agnts)];
Fdes = [agnts[i].Fdes .* agnts[i].Fdes_sc for i in 1:length(agnts)];

# run the experiment with Algorithm 1 (c_i = d_i)
DEP_params_A1 = copy(DEP_params);
DEP_params_A1["maxDist"] = (300e-6, 300e-6);
params_A1, exp_data_A1 = runExp(platform=:DEP, N_iter=N_iter, N_agnts=N_agnts, N_acts=(N_acts,N_acts),
                            convanalysis=true,
                            saveplots=false,
                            showplots=false,
                            agent_positions = agent_positions,
                            Fdes = Fdes,
                            params = DEP_params_A1);

# Compute the developed force by the phase shifts found by Alg. 2
Fdev_arr_A2 = Array{NTuple{params["force_dim"], Float64}, 1}();
controls_A2 = exp_data_A2[1]["ActuatorCommands"];
actuatorArray_A2 = exp_data_A2[1]["ActuatorArray"];
for agnt_pos in agent_positions
    Fdev = params["calcForce"](actuatorArray_A2, agnt_pos, controls_A2);
    push!(Fdev_arr_A2, Fdev);
end

# Compute the developed force by the phase shifts found by Alg. 1
Fdev_arr_A1 = Array{NTuple{params["force_dim"], Float64}, 1}();
controls_A1 = exp_data_A1[1]["ActuatorCommands"];
actuatorArray_A1 = exp_data_A1[1]["ActuatorArray"];
for agnt_pos in agent_positions
    Fdev2 = params["calcForce"](actuatorArray_A1, agnt_pos, controls_A1);
    push!(Fdev_arr_A1, Fdev2);
end

# Compare the developed forces in a plot
Plots.display(plot(actuatorArray_A2, agents=agnts,
                    dev_forces=Fdev_arr_A2,
                    aux_forces=Fdev_arr_A1,
                    showReqForces=true,
                    Fscale=params_A1["Fscale"]*1.5,
                    legend=nothing,
                    size=(500,500)))
##
savefig("DEP_Alg1_vs_Alg2.pdf")