text:
_smallwerr(w) = throw(DomainError(w,"window size `w` must be greater than zero"))

function findnextextrema(cmp, x::AbstractVector, i::Int, w::Int, strict::Bool)
    w > 0 || _smallwerr(w)
    strict && i - w < firstindex(x) && (i = w + 1) # First peak can't be closer than w
    maxlast = strict ? lastindex(x) - w : lastindex(x)

    while i ≤ maxlast
        peak = true
        xi = x[i]
        if ismissing(xi) || isnan(xi)
            i = something(findnext(p -> !isnan(p) & !ismissing(p), x, i+1), lastindex(x)+1)
        else
            for j in max(firstindex(x),i-w):min(i+w,lastindex(x))
                j === i && continue
                xj = x[j]
                if coalesce(cmp(xi, xj), strict) || (strict && isnan(xj))
                    peak &= false
                    i = j > i ? j : i + 1
                    break # No reason to continue checking if the rest of the elements qualify
                elseif xi === xj # plateau
                    k = findnext(p -> xi !== p, x, j+1)
                    # @show k i j
                    if strict
                        if !isnothing(k) && i-1 ≥ firstindex(x)
                            xi1 = x[i-1]
                            xk = x[k]
                            # @show xi1 xi xk i j k
                            if !ismissing(xi1) && !isnan(xi1) && cmp(xi1, xi) &&
                                !ismissing(xk) && !isnan(xk) && cmp(xk, xi)
                                # plateau begins and ends with elements less than the current
                                return i
                            else
                                peak &= false
                                i = k > i ? k : i+1
                                break
                            end
                        else
                            peak &= false
                            i = isnothing(k) ? lastindex(x)+1 :
                                                    k > i ? k : i + 1
                            break
                        end
                    else
                        xk = isnothing(k) ? NaN : x[k]
                        xi1 = i-1 ≥ firstindex(x) ? x[i-1] : NaN
                        if something(k, lastindex(x)+1) ≥ i+w &&
                            (ismissing(xi1) || isnan(xi1) || cmp(xi1, xi)) &&
                            (ismissing(xk) || isnan(xk) || cmp(xk, xi))
                            return i
                        elseif !isnothing(k)
                            if k > i
                                peak &= false
                                i = k
                                break
                            end
                        else
                            peak &= false
                            i = lastindex(x)+1
                            break
                        end
                    end
                end
            end
            peak && return i
        end
    end

    return lastindex(x) + 1
end

"""
    findnextmaxima(x, i[, w=1; strict=true])

Find the index of the next maxima in `x` after or including `i`, where the maxima `i` is
either the maximum of `x[i-w:i+w]` or the first index of a plateau. Returns `lastindex(x) +
1` if no maxima occur after `i`.

If `strict` is `true`, all elements in `x[i-w:i+w]` or the bounds of a plateau must
exist and must not be `missing` or `NaN`. For `strict == false`, a maxima is the
maximum or first element of all existing and non-NaN or missing elements in `x[i-w:i+w]` or
the bounds of a plateau.

See also: [`argmaxima`](@ref)

# Examples
```jldoctest
julia> findnextmaxima([0,2,0,1,1,0], 2)
2

julia> findnextmaxima([0,2,0,1,1,0], 3)
4

```
"""
findnextmaxima(x, i, w=1; strict=true) = findnextextrema(<, x, i, w, strict)

"""
    argmaxima(x[, w=1; strict=true])

Find the indices of the local maxima of `x` where each maxima `i` is either the maximum of
`x[i-w:i+w]` or the first index of a plateau.

If `strict` is `true`, all elements of `x[i-w:i+w]` or the bounds of a plateau must
exist and must not be `missing` or `NaN`. For `strict == false`, a maxima is the
maximum or first element of all existing and non-NaN or missing elements in `x[i-w:i+w]` or
the bounds of a plateau.

See also: [`findmaxima`](@ref), [`findnextmaxima`](@ref)

# Examples
```jldoctest
julia> argmaxima([0,2,0,1,1,0])
2-element Vector{Int64}:
 2
 4

julia> argmaxima([2,0,1,1])
Int64[]

julia> argmaxima([2,0,1,1]; strict=false)
2-element Vector{Int64}:
 1
 3
```
"""
function argmaxima(
    x::AbstractVector{T}, w::Int=1; strict::Bool=true
) where T
    w > 0 || _smallwerr(w)
    pks = Int[]

    i = strict ? w + firstindex(x) : firstindex(x)
    i = findnextmaxima(x, i, w; strict=strict)
    while i ≤ lastindex(x)
        push!(pks, i)
        i = findnextmaxima(x, i+w+1, w; strict=strict)
    end

    return pks
end

"""
    maxima(x[, w=1; strict=true])

Find the values of the local maxima of `x` where each maxima `i` is either the maximum of
`x[i-w:i+w]` or the first index of a plateau.

See also: [`argmaxima`](@ref), [`findnextmaxima`](@ref)
"""
function maxima(
    x::AbstractVector{T}, w::Int=1; strict::Bool=true
) where T
    idxs = argmaxima(x, w; strict=strict)
    return x[idxs]
end

"""
    findmaxima(x[, w=1; strict=true]) -> (idxs, vals)

Find the indices and values of the local maxima of `x` where each maxima `i` is either the maximum of
`x[i-w:i+w]` or the first index of a plateau.

See also: [`argmaxima`](@ref), [`findnextmaxima`](@ref)
"""
function findmaxima(x, w::Int=1; strict::Bool=true)
    idxs = argmaxima(x, w; strict=strict)
    return (idxs, x[idxs])
end

"""
    findnextminima(x, i[, w=1, strict=true])

Find the index of the next minima in `x` after or including `i`, where the minima `i` is
either the minimum of `x[i-w:i+w]` or the first index of a plateau. Returns `lastindex(x) +
1` if no minima occur after `i`.

If `strict` is `true`, all elements in `x[i-w:i+w]` or the bounds of a plateau must
exist and must not be `missing` or `NaN`. For `strict == false`, a minima is the
minimum or first element of all existing and non-NaN or missing elements in `x[i-w:i+w]` or
the bounds of a plateau.

See also: [`argminima`](@ref)

# Examples
```jldoctest
julia> findnextminima([3,2,3,1,1,3], 2)
2

julia> findnextminima([3,2,3,1,1,3], 3)
4

```
"""
findnextminima(x, i, w=1; strict=true) = findnextextrema(>, x, i, w, strict)

"""
    argminima(x[, w=1; strict=false])

Find the indices of the local minima of `x` where each minima `i` is either the minimum of
`x[i-w:i+w]` or the first index of a plateau.

If `strict` is `true`, all elements of `x[i-w:i+w]` or the bounds of a plateau must
exist and must not be `missing` or `NaN`. For `strict == false`, a minima is the
minimum or first element of all existing and non-NaN or missing elements in `x[i-w:i+w]` or
the bounds of a plateau.

See also: [`findminima`](@ref), [`findnextminima`](@ref)

# Examples
```jldoctest
julia> argminima([3,2,3,1,1,3])
2-element Vector{Int64}:
 2
 4

julia> argminima([2,3,1,1])
Int64[]

julia> argminima([2,3,1,1]; strict=false)
2-element Vector{Int64}:
 1
 3
```
"""
function argminima(
    x::AbstractVector{T}, w::Int=1; strict::Bool=true
) where T
    w > 0 || _smallwerr(w)
    pks = Int[]

    i = strict ? w + firstindex(x) : firstindex(x)
    i = findnextminima(x, i, w; strict=strict)
    while i ≤ lastindex(x)
        push!(pks, i)
        i = findnextminima(x, i+w+1, w; strict=strict)
    end

    return pks
end

"""
    minima(x[, w=1; strict=true])

Find the values of the local minima of `x` where each minima `i` is either the minimum of
`x[i-w:i+w]` or the first index of a plateau.

See also: [`argminima`](@ref), [`findnextminima`](@ref)
"""
function minima(
    x::AbstractVector{T}, w::Int=1; strict::Bool=true
) where T
    idxs = argminima(x, w; strict=strict)
    return x[idxs]
end

"""
    findminima(x[, w=1; strict=true]) -> (idxs, vals)

Find the indices and values of the local minima of `x` where each minima `i` is either the minimum of
`x[i-w:i+w]` or the first index of a plateau.

See also: [`argminima`](@ref), [`findnextminima`](@ref)
"""
function findminima(x, w::Int=1; strict::Bool=true)
    idxs = argminima(x, w; strict=strict)
    return (idxs, x[idxs])
end

max_stars_repo_name: halleysfifthinc/Peaks

max_issues_repo_name: halleysfifthinc/Peaks

max_forks_repo_name: halleysfifthinc/Peaks


text:
# Load Julia packages (libraries) needed.

using StatisticalRethinking
using StructuralCausalModels

ProjDir = @__DIR__
cd(ProjDir) #do

# ### snippet 5.1

df1 = DataFrame!(CSV.File(rel_path("..", "data", "WaffleDivorce.csv"), delim=';'));

df = DataFrame()
df[!, :A] = df1[:, :MedianAgeMarriage]
df[!, :M] = df1[:, :Marriage]
df[!, :D] = df1[:, :Divorce]

d = OrderedDict(
  :M => [:A],
  :D => [:A]
);

dag = DAG("waffles", d, df=df);
show(dag)

fname = ProjDir * "/AMD_3.dot"
to_graphviz(dag, fname)
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

display(dag.s); println()

shipley_test_dag_3 = shipley_test(dag)
if !isnothing(shipley_test_dag_3)
  display(shipley_test_dag_3)
end
println()

f = [:M]; s = [:D]; sel = vcat(f, s)
c = [:A]

e = d_separation(dag, f, s)
println("d_separation($(dag.name), $f, $s) = $e")

e = d_separation(dag, f, s, c=c)
println("d_separation($(dag.name), $f, $s, c=$c) = $e\n")

bs = basis_set(dag)
display(bs)

adjustmentsets = adjustment_sets(dag, f[1], s[1])
println("Adjustment sets:")
adjustmentsets |> display

#end

ext: jl


text:
## Parquet

@info "Parquet.jl is available - including functionality to read / write Parquet files"

using .Parquet

"""
For strings (plus potentially additional data types), the binary representation in Parquet must be specified when reading the file.
A general mapping can be given with the keyword argument `map_logical_types`.
Alternatively, the string columns can be given

"""
function read_table(::ParquetFormat, filename; kwargs...)
    parfile = Parquet.File(filename; kwargs...)
    return RecordCursor(parfile)
end

function write_table!(::ParquetFormat, filename, table; kwargs...)
    _checktable(table)
    write_parquet(filename, table; kwargs...)
    nothing
end


text:
# Create Roadway
roadway_opts = MergingRoadwayOptions(lane_width=0.26)
roadway = build_roadway(roadway_opts)

 # Create Players
players = Vector{Player{T}}(undef, p)
player_opts = HighwayPlayerOptions{model.ni[1],model.mi[1],T}()
players[1] = Player(model, roadway.lane[5])
players[2] = Player(model, roadway.lane[3])

# Create Scenario
sce = Scenario(model, roadway, players)

# Create players
T = Float64
p = 2
model = DoubleIntegratorGame(p=p)
n = model.n
m = model.m

# Define the horizon of the problem
N = 20 # N time steps
dt = 0.1 # each step lasts 0.1 second
probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

# Define the objective of each player
# We use a LQR cost
Q = [Diagonal(3*[0.0, 1.0, 2.0, 1.0]) for i=1:p] # Quadratic state cost
R = [Diagonal(0.3*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
# Desrired state
xf = [SVector{model.ni[1],T}([1.0,+0.35,0.5,0]),
      SVector{model.ni[2],T}([1.0,-0.05,0.5,0]),
      ]
# Desired control
uf = [zeros(SVector{model.mi[i],T}) for i=1:p]
# Objectives of the game
game_obj = Algames.GameObjective(Q,R,xf,uf,N,model)
radius = 0.2*ones(p)
μ = 15.0*ones(p)
add_collision_cost!(game_obj, radius, μ)

# Define the constraints that each player must respect
game_con = Algames.GameConstraintValues(probsize)
game_con = GameConstraintValues(N, sce)

# Define the initial state of the system
x0 = SVector{model.n,T}([
    0.2,   0.2,
   -0.45, -0.05,
    0.4,   0.5,
    0.1,   0.0,
    ])

# Define the Options of the solver
opts = Options()
opts.ls_iter = 15
opts.outer_iter = 20
opts.inner_iter = 20
opts.ρ_0 = 1e0
opts.reg_0 = 1e-8
opts.α_dual = 0.0
opts.λ_max = 1.0*1e7
opts.ϵ_dyn = 1e-6
opts.ϵ_sta = 1e-6
opts.ϵ_con = 1e-6
opts.ϵ_opt = 1e-6
opts.regularize = true
# Define the game problem
prob = Algames.GameProblem(N,dt,x0,model,opts,game_obj,game_con)

# Solve the problem
newton_solve!(prob)

plot!(prob.model, prob.pdtraj.pr)
plot!(prob.stats)

function mpc_loop(prob::Algames.GameProblem, x0::SVx, sce::Scenario, vis::Vector{Visualizer},
    mpc_opts::MPCOptions{T}, mpc_stats::MPCStatistics{T}, traj_stats::TrajStatistics{T};
    seed::Vector{Int}=[1,6]) where {T,SVx}
    # We select the seed to make sure both players start with a different local nash equilibrium.
    probsize = prob.probsize
    p = probsize.p
    pz = probsize.pz

    # Initialization
    reset!(mpc_stats)
    reset!(traj_stats)
    probs = [deepcopy(prob) for i=1:p]

    prob.opts.shift = 2^10
    prob.opts.f_init = x -> (rand(x) .- 0.5)
    prob.opts.amplitude_init = 3.0

    # Set tolerances
    prob.opts.ϵ_opt = 3e-3
    prob.opts.ϵ_dyn = 3e-3
    prob.opts.ϵ_con = 3e-3
    prob.opts.ϵ_sta = 3e-3
    prob.opts.inner_print = false
    prob.opts.outer_print = false

    for i = 1:p
        probs[i].opts.seed = seed[1+(i-1)%p]
        probs[i].x0 = x0

        newton_solve!(probs[i])
        probs[i].opts.shift = 0
        probs[i].opts.amplitude_init = 1e-8
        plot!(probs[i].model, probs[i].pdtraj.pr)
    end

    dt = zeros(p)
    for j = 1:mpc_opts.M
        println("$(round(mpc_stats.t,digits=2))/$(round(mpc_opts.h, digits=1))")
        trajs = deepcopy([probs[i].pdtraj.pr for i=1:p])
        # Solve
        for i = 1:p
            set_state!(vis[i], probs[i].model, sce, probs[i].x0)
            set_line_traj!(vis[i], probs[i].model, sce.player, probs[i].pdtraj.pr)
            dt[i] = @elapsed newton_solve!(probs[i])
        end
        # Record
        z0 = probs[1].pdtraj.pr[1]
        z = Algames.KnotPoint(Algames.state(z0), Algames.control(z0), minimum(dt))
        record!(traj_stats, z.dt, trajs)
        @show traj_stats.iter
        @show length(traj_stats.traj[1])

        record!(mpc_stats, z)
        if mpc_stats.t > mpc_opts.h
            break
        end

        # Update solver
        X = []
        S = []
        for i = 1:p
            x,s = first_order_rollout(probs[i].pdtraj.pr, minimum(dt))
            push!(X, x)
            push!(S, s)
        end
        x = zeros(n)
        for i = 1:p
            x[pz[i]] = X[i][pz[i]]
        end
        x = SVector{n,T}(x)
        for i = 1:p
            probs[i].x0 = x
            probs[i].opts.shift = minimum(S)
        end
    end
    return nothing
end

# # Initialize visualizers
# vis = [Visualizer() for i=1:p]
# open.(vis)

# # Visualize trajectories
# for i = 1:p
#     set_scenario!(vis[i], sce)
#     set_env!(vis[i], VehicleState(1.0, 0.0, 0.0, 0.0))
#     set_camera_birdseye!(vis[i], height=8.0)
# end

mpc_opts = MPCOptions(M=100, h=4.0)
mpc_stats = MPCStatistics()
traj_stats = TrajStatistics(probsize)

mpc_loop(prob, x0, sce, vis, mpc_opts, mpc_stats, traj_stats, seed=[5,13])

for t = 1:traj_stats.iter
    for i = 1:p
        set_line_traj!(vis[i], model, sce.player, traj_stats.traj[t][i])
        sleep(traj_stats.dt[t])
    end
end

# for i = 1:15
#     @show i
#     mpc_loop(prob, x0, sce, vis, mpc_opts, mpc_stats, seed=[i,i])
#     sleep(0.5)
# end

# init_traj, sol_traj = get_figure(prob, sce, vis_2, vis_1)


text:
__precompile__()

"""
Functions to learn a metamotif model of nucleosome positions by nested sampling
"""

module nnlearn
    using BenchmarkTools, BioSequences, DataFrames, Distributed, Distributions, ProgressMeter, Serialization, UnicodePlots, JLD2
    import ProgressMeter: AbstractProgress
    import Printf: @sprintf
    import StatsFuns: logaddexp, logsumexp #both are needed as logsumexp for two terms is deprecated
    import Random: rand, seed!, shuffle!
    import Distances: euclidean

    mutable struct Model_Record #record struct to associate a log_Li with a saved, calculated model
        path::String
        log_Li::AbstractFloat
    end

    function read_fa_wms_tr(path::String)
        wms=Vector{Matrix{AbstractFloat}}()
        wm=zeros(1,4)
        f=open(path)
        for line in eachline(f)
            prefix=line[1:2]
            prefix == "01" && (wm=transpose([parse(AbstractFloat,i) for i in split(line)[2:end]]))
            prefix != "01" && prefix != "NA" && prefix != "PO" && prefix != "//" && (wm=vcat(wm, transpose([parse(AbstractFloat,i) for i in split(line)[2:end]])))
            prefix == "//" && push!(wms, wm)
        end
        return wms
    end

    #wm_samples are in decimal probability space, not log space
    function assemble_source_priors(no_sources::Integer, wm_samples::Vector{Matrix{AbstractFloat}}, prior_wt::AbstractFloat, uninform_length_range::UnitRange{Integer}, unspecified_wm_prior=Dirichlet(ones(4)/4)) #estimate a dirichlet prior on wm_samples inputs; unless otherwise specified, fill non-wm_samples-specified priors with uninformative priors of random lengths in the specified range
        source_priors = Vector{Vector{Dirichlet{AbstractFloat}}}()
        for source in 1:no_sources
            if source <= length(wm_samples)
                push!(source_priors, estimate_dirichlet_prior_on_wm(wm_samples[source], prior_wt))
            else
                uninformative_prior = Vector{Dirichlet{AbstractFloat}}()
                for pos in 1:rand(uninform_length_range)
                    push!(uninformative_prior, unspecified_wm_prior)
                end
                push!(source_priors, uninformative_prior)
            end
        end
        return source_priors
    end

                function estimate_dirichlet_prior_on_wm(wm::Matrix{AbstractFloat}, wt::AbstractFloat)
                    for i in 1:size(wm)[1]
                        !(isprobvec(wm[i,:])) && throw(DomainError("Bad weight vec supplied to estimate_dirichlet_prior_on_wm! $(wm[i,:])"))
                    end
                    prior = Vector{Dirichlet{AbstractFloat}}()
                    for position in 1:size(wm)[1]
                        normvec=wm[position,:]
                        zero_idxs=findall(isequal(0.),wm[position,:])
                        normvec[zero_idxs].+=10^-99
                        push!(prior, Dirichlet(normvec.*wt))
                    end
                    return prior
                end

    function cluster_mix_prior!(df::DataFrame, wms::Vector{Matrix{AbstractFloat}})
        mix=falses(size(df,1),length(wms))
        for (o, row) in enumerate(eachrow(df))
            row.cluster != 0 && (mix[o,row.cluster]=true)
        end

        represented_sources=unique(df.cluster)
        wms=wms[represented_sources]
        return mix[:,represented_sources]
    end

    #subfuncs to handle sums of log probabilities that may include -Inf (ie p=0), returning -Inf in this case rather than NaNs
    function lps(adjuvants::AbstractArray)
        prob = sum(adjuvants) ; isnan(prob) ? - Inf : prob
    end

    function lps(base, adjuvants...)::AbstractFloat
        prob = base+sum(adjuvants) ; isnan(prob) ? -Inf : prob
    end

    function infocenter_wms_trim(wm::Matrix{AbstractFloat}, trimsize::Integer)
        !(size(wm,2)==4) && throw(DomainError("Bad wm! 2nd dimension should be size 4"))
        infovec=get_pwm_info(wm, logsw=false)
        maxval, maxidx=findmax(infovec)
        upstream_extension=Int(floor((trimsize-1)/2))
        downstream_extension=Int(ceil((trimsize-1)/2))
        1+upstream_extension+downstream_extension > size(wm,1) && throw(DomainError("Src too short for trim! $upstream_extension $downstream_extension"))
        return wm[max(1,maxidx-upstream_extension):min(maxidx+downstream_extension,size(wm,1)),:]
    end

    function filter_priors(target_src_no::Integer, target_src_size::Integer, prior_wms::Vector{Matrix{AbstractFloat}}, prior_mix::BitMatrix)
        wms=Vector{Matrix{AbstractFloat}}(undef, target_src_no)
        freqsort_idxs=sortperm([sum(prior_mix[:,s]) for s in 1:length(prior_wms)])
        for i in 1:target_src_no
            target_src_idx=freqsort_idxs[i]
            wms[i]=infocenter_wms_trim(prior_wms[target_src_idx], target_src_size)
        end
        return wms
    end

    function combine_filter_priors(target_src_no::Integer, target_src_size::Integer, prior_wms::Tuple{Vector{Matrix{AbstractFloat}},Vector{Matrix{AbstractFloat}}}, prior_mix::Tuple{BitMatrix,BitMatrix})
        wms=Vector{Matrix{AbstractFloat}}(undef, target_src_no)
        cat_wms=vcat(prior_wms[1],prior_wms[2])
        first_freq=[sum(prior_mix[1][:,s]) for s in 1:length(prior_wms[1])]
        second_freq=[sum(prior_mix[2][:,s]) for s in 1:length(prior_wms[2])]
        freqsort_idxs=sortperm(vcat(first_freq,second_freq))
        for i in 1:target_src_no
            target_src_idx=freqsort_idxs[i]
            wms[i]=infocenter_wms_trim(cat_wms[target_src_idx], target_src_size)
        end
        return wms
    end

    include("ICA_PWM_model.jl")
    include("Bayes_IPM_ensemble.jl")
    include("nested_sampler.jl")
    include("model_display.jl")
    include("ns_progressmeter.jl")
    include("ensemble_utils.jl")
    include("permute_control.jl")
    include("permute_instructions.jl")

    #include("performance.jl")
end # module


text:
export network, NetworkObservable

using Sockets

import Base: ==
import Base: show

"""
    NetworkObservable{D, Address, Port, S}()

NetworkObservable listens for the messages of type `D` from remote server with specified `Address` and `Port` parameters.

# See also: [`network`](@ref), [`Subscribable`](@ref)
"""
@subscribable struct NetworkObservable{D, A, P, S} <: Subscribable{D} end

function on_subscribe!(observable::NetworkObservable{D, Address, Port, S}, actor) where { D, Address, Port, S }
    clientside = Sockets.connect(Address, Port)
    @async begin
        try
            while isopen(clientside)
                count = read(clientside, Int)
                if count === 0
                    complete!(actor)
                    return nothing
                elseif count === -1
                    error!(actor, ErrorException("NetworkObservableError"))
                    return nothing
                else
                    next!(actor, read(clientside, D))
                end
            end
        catch err
            if !(err isa EOFError)
                error!(actor, err)
            end
        end
    end
    return NetworkObservableSubscritpion(clientside)
end

function on_subscribe!(observable::NetworkObservable{Vector{D}, Address, Port, S}, actor) where { D, Address, Port, S }
    clientside = Sockets.connect(Address, Port)
    @async begin
        try
            buffer = Vector{D}(undef, S)
            while isopen(clientside)
                count = read(clientside, Int)
                if count === 0
                    complete!(actor)
                    return nothing
                elseif count === -1
                    error!(actor, ErrorException("NetworkObservableError"))
                    return nothing
                else
                    unsafe_read(clientside, pointer(buffer), count * sizeof(D))
                    next!(actor, unsafe_wrap(Vector{D}, pointer(buffer), count))
                end
            end
        catch err
            if !(err isa EOFError)
                error!(actor, err)
            end
        end
    end
    return NetworkObservableSubscritpion(clientside)
end

struct NetworkObservableSubscritpion <: Teardown
    clientside :: TCPSocket
end

as_teardown(::Type{<:NetworkObservableSubscritpion}) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::NetworkObservableSubscritpion)
    close(subscription.clientside)
end

"""
    network(::Type{D}, port::Int)             where D
    network(::Type{D}, address::A, port::Int) where { D, A <: IPAddr }

    network(::Type{Vector{D}}, port::Int, buffer_size::Int)             where D
    network(::Type{Vector{D}}, address::A, port::Int, buffer_size::Int) where { D, A <: IPAddr }

Creation operator for the `NetworkObservable` that emits messages from the server with specified `address` and `port` arguments.

See also: [`NetworkObservable`](@ref), [`subscribe!`](@ref)
"""
network(::Type{D}, port::Int)             where D                  = network(D, Sockets.localhost, port)
network(::Type{D}, address::A, port::Int) where { D, A <: IPAddr } = NetworkObservable{D, address, port, 0}()

network(::Type{Vector{D}}, port::Int)             where D                  = error("Specify maximum buffer size for input data")
network(::Type{Vector{D}}, address::A, port::Int) where { D, A <: IPAddr } = error("Specify maximum buffer size for input data")

network(::Type{Vector{D}}, port::Int, buffer_size::Int)             where D                  = network(Vector{D}, Sockets.localhost, port, buffer_size)
network(::Type{Vector{D}}, address::A, port::Int, buffer_size::Int) where { D, A <: IPAddr } = NetworkObservable{Vector{D}, address, port, buffer_size}()

Base.:(==)(::NetworkObservable{D1, A1, P1}, ::NetworkObservable{D2, A2, P2}) where { D1, A1, P1 } where { D2, A2, P2 } = D1 == D2 && A1 == A2 && P1 == P2

Base.show(io::IO, observable::NetworkObservable{D, A, P}) where { D, A, P } = print(io, "NetworkObservable($D, address = $A, port = $P)")


text:
# you will need this to use Tropical numbers in CUDA.
using CUDA
CUDA.allowscalar(false)
using TropicalGEMM: XTranspose, NativeTypes, Tropical, TropicalTypes, CountingTropical
using LinearAlgebra
using Test

const CTranspose{T} = Transpose{T, <:StridedCuVecOrMat{T}}
for TT in [:(Tropical{<:NativeTypes}), :TropicalTypes]
    for RT in [TT, :Real]
        for (TA, CTA) in [(:CuMatrix, :CuMatrix), (:CTranspose, :(Transpose{<:Any, <:StridedCuVecOrMat}))]
            for (TB, CTB) in [(:CuMatrix, :CuMatrix), (:CTranspose, :(Transpose{<:Any, <:StridedCuVecOrMat}))]
                @eval function LinearAlgebra.mul!(o::CuMatrix{T}, a::$TA{T}, b::$TB{T}, α::$RT, β::$RT) where {T<:$TT}
                    CUDA.CUBLAS.gemm_dispatch!(o, a, b, α, β)
                end
            end
        end
    end
end

@testset "cuda patch" begin
    for T in [Tropical{Float64}, CountingTropical{Float64,Float64}]
        a = T.(CUDA.randn(4, 4))
        b = T.(CUDA.randn(4))
        for A in [transpose(a), a, transpose(b)]
            for B in [transpose(a), a, b]
                if !(size(A) == (1,4) && size(B) == (4,))
                    res0 = Array(A) * Array(B)
                    res1 = A * B
                    res2 = mul!(CUDA.zeros(T, size(res0)...), A, B, true, false)
                    @test Array(res1) ≈ res0
                    @test Array(res2) ≈ res0
                end
            end
        end
    end
end


text:
# Functions for querying network properties.

# DEPRECIATED after v9.0
function params(network)
    Base.depwarn("`params` is depreciated, please use `ModelingToolkit.parameters` for all system and subsystem parameters, or `reactionparams` for all parameters within system and subsystem `Reaction`s.", :params, force=true)
    parameters(network)
end

function numparams(network)
    Base.depwarn("`numparams` is depreciated, please use `length(ModelingToolkit.parameters)` for the total number of parameters across all systems and subsystems, or `numreactionparams` for the number of parameters within system and subsystem `Reaction`s.", :params, force=true)
    length(parameters(network))
end

"""
    merge(network1::ReactionSystem, network2::ReactionSystem)

Create a new network merging `network1` and `network2`.

Notes:
- Duplicate reactions between the two networks are not filtered out.
- [`Reaction`](@ref)s are not deepcopied to minimize allocations, so the new
  network will share underlying data arrays.
- Subsystems are not deepcopied between the two networks and will hence be
  shared.
- Returns the merged network.
"""
function Base.merge(network1::ReactionSystem, network2::ReactionSystem)
    Base.depwarn("`merge(sys1::ReactionSystem, sys2::ReactionNetwork)` is depreciated, please use `ModelingToolkit.extend` instead.", :merge, force=true)
    network = make_empty_network()
    merge!(network, network1)
    merge!(network, network2)
    network
end

######### Accessors: #########

function filter_nonrxsys(network)
    systems   = get_systems(network)
    rxsystems = ReactionSystem[]
    for sys in systems
        (sys isa ReactionSystem) && push!(rxsystems, sys)
    end
    rxsystems
end

function species(network::ReactionSystem, sts)
    [MT.renamespace(network,st) for st in sts]
end

"""
    species(network)

Given a [`ReactionSystem`](@ref), return a vector of all species defined in the system
and any subsystems that are of type `ReactionSystem`. To get the variables in the system
and all subsystems, including non-`ReactionSystem` subsystems, uses `states(network)`.

Notes:
- If `ModelingToolkit.get_systems(network)` is non-empty will allocate.
"""
function species(network)
    sts = get_states(network)
    systems = filter_nonrxsys(network)
    isempty(systems) && return sts
    unique([sts; reduce(vcat, map(sys -> species(sys,species(sys)), systems))])
end

"""
    reactionparams(network)

Given a [`ReactionSystem`](@ref), return a vector of all parameters defined
within the system and any subsystems that are of type `ReactionSystem`. To get
the parameters in the system and all subsystems, including non-`ReactionSystem`
subsystems, use `parameters(network)`.

Notes:
- If `ModelingToolkit.get_systems(network)` is non-empty will allocate.
"""
function reactionparams(network)
    ps = get_ps(network)
    systems = filter_nonrxsys(network)
    isempty(systems) && return ps
    unique([ps; reduce(vcat, map(sys -> species(sys,reactionparams(sys)), systems))])
end

function namespace_reactions(network::ReactionSystem)
    rxs = reactions(network)
    isempty(rxs) && return Reaction[]
    map(rx -> namespace_equation(rx, network), rxs)
end

"""
    reactions(network)

Given a [`ReactionSystem`](@ref), return a vector of all `Reactions` in the system.

Notes:
- If `ModelingToolkit.get_systems(network)` is not empty, will allocate.
"""
function reactions(network)
    rxs = get_eqs(network)
    systems = filter_nonrxsys(network)
    isempty(systems) && (return rxs)
    [rxs; reduce(vcat, namespace_reactions.(systems); init=Reaction[])]
end

"""
    speciesmap(network)

Given a [`ReactionSystem`](@ref), return a Dictionary mapping species that
participate in `Reaction`s to their index within [`species(network)`](@ref).
"""
function speciesmap(network)
    Dict(S => i for (i,S) in enumerate(species(network)))
end

"""
    paramsmap(network)

Given a [`ReactionSystem`](@ref), return a Dictionary mapping from all
parameters that appear within the system to their index within
`parameters(network)`.
"""
function paramsmap(network)
    Dict(p => i for (i,p) in enumerate(parameters(network)))
end

"""
    reactionparamsmap(network)

Given a [`ReactionSystem`](@ref), return a Dictionary mapping from parameters that
appear within `Reaction`s to their index within [`reactionparams(network)`](@ref).
"""
function reactionparamsmap(network)
    Dict(p => i for (i,p) in enumerate(reactionparams(network)))
end

"""
    numspecies(network)

Return the total number of species within the given [`ReactionSystem`](@ref) and
subsystems that are `ReactionSystem`s.

Notes
- If there are no subsystems this will be fast.
- As this calls [`species`](@ref), it can be slow and will allocate if there are
  any subsystems.
"""
function numspecies(network)
    length(species(network))
end

"""
    numreactions(network)

Return the total number of reactions within the given [`ReactionSystem`](@ref)
and subsystems that are `ReactionSystem`s.
"""
function numreactions(network)
    nr = length(get_eqs(network))
    for sys in get_systems(network)
        (sys isa ReactionSystem) && (nr += numreactions(sys))
    end
    nr
end

"""
    numreactionparams(network)

Return the total number of parameters within the given [`ReactionSystem`](@ref)
and subsystems that are `ReactionSystem`s.

Notes
- If there are no subsystems this will be fast.
- As this calls [`reactionparams`](@ref), it can be slow and will allocate if
  there are any subsystems.
"""
function numreactionparams(network)
    length(reactionparams(network))
end

"""
    dependents(rx, network)

Given a [`Reaction`](@ref) and a [`ReactionSystem`](@ref), return a vector of
`ModelingToolkit.Num`s corresponding to species the *reaction rate
law* depends on. E.g., for

`k*W, 2X + 3Y --> 5Z + W`

the returned vector would be `[W(t),X(t),Y(t)]`.

Notes:
- Allocates
- Does not check for dependents within any subsystems.
"""
function dependents(rx, network)
    if rx.rate isa Number
        return rx.substrates
    else
        rvars = get_variables(rx.rate, species(network))
        return union!(rvars, rx.substrates)
    end
end

"""
    dependents(rx, network)

See documentation for [`dependents`](@ref).
"""
function dependants(rx, network)
    dependents(rx, network)
end

"""
    reactionrates(network)

Given a [`ReactionSystem`](@ref), returns a vector of the symbolic reaction
rates for each reaction.
"""
function reactionrates(rn)
    [r.rate for r in reactions(rn)]
end

"""
    substoichmat(rn; sparse=false, smap=speciesmap(rn))

Returns the substrate stoichiometry matrix, ``S``, with ``S_{i j}`` the
stoichiometric coefficient of the ith substrate within the jth reaction.

Note:
- Set sparse=true for a sparse matrix representation
"""
function substoichmat(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem; smap=speciesmap(rn))
    Is=Int[];  Js=Int[];  Vs=Int[];
    for (k,rx) in enumerate(reactions(rn))
        stoich = rx.substoich
        for (i,sub) in enumerate(rx.substrates)
            push!(Js, k)
            push!(Is, smap[sub])
            push!(Vs, stoich[i])
        end
    end
    sparse(Is,Js,Vs,numspecies(rn),numreactions(rn))
end
function substoichmat(::Type{Matrix{Int}},rn::ReactionSystem; smap=speciesmap(rn))
    smat = zeros(Int, numspecies(rn), numreactions(rn))
    for (k,rx) in enumerate(reactions(rn))
        stoich = rx.substoich
        for (i,sub) in enumerate(rx.substrates)
            smat[smap[sub],k] = stoich[i]
        end
    end
    smat
end
function substoichmat(rn::ReactionSystem; sparse::Bool=false, smap=speciesmap(rn))
    isempty(get_systems(rn)) || error("substoichmat does not currently support subsystems.")
	sparse ? substoichmat(SparseMatrixCSC{Int,Int}, rn; smap=smap) : substoichmat(Matrix{Int}, rn; smap=smap)
end

"""
    prodstoichmat(rn; sparse=false, smap=speciesmap(rn))

Returns the product stoichiometry matrix, ``P``, with ``P_{i j}`` the
stoichiometric coefficient of the ith product within the jth reaction.

Note:
- Set sparse=true for a sparse matrix representation
"""
function prodstoichmat(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem; smap=speciesmap(rn))
    Is=Int[];  Js=Int[];  Vs=Int[];
    for (k,rx) in enumerate(reactions(rn))
        stoich = rx.prodstoich
        for (i,prod) in enumerate(rx.products)
			push!(Js, k)
			push!(Is, smap[prod])
			push!(Vs, stoich[i])
        end
    end
    sparse(Is,Js,Vs,numspecies(rn),numreactions(rn))
end
function prodstoichmat(::Type{Matrix{Int}},rn::ReactionSystem; smap=speciesmap(rn))
    pmat = zeros(Int, numspecies(rn), numreactions(rn))
    for (k,rx) in enumerate(reactions(rn))
        stoich = rx.prodstoich
        for (i,prod) in enumerate(rx.products)
            pmat[smap[prod],k] = stoich[i]
        end
    end
    pmat
end
function prodstoichmat(rn::ReactionSystem; sparse=false, smap=speciesmap(rn))
    isempty(get_systems(rn)) || error("prodstoichmat does not currently support subsystems.")
	sparse ? prodstoichmat(SparseMatrixCSC{Int,Int}, rn; smap=smap) : prodstoichmat(Matrix{Int}, rn; smap=smap)
end

"""
    netstoichmat(rn, sparse=false; smap=speciesmap(rn))

Returns the net stoichiometry matrix, ``N``, with ``N_{i j}`` the net
stoichiometric coefficient of the ith species within the jth reaction.

Note:
- Set sparse=true for a sparse matrix representation
"""
function netstoichmat(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem; smap=speciesmap(rn))
    Is=Int[];  Js=Int[];  Vs=Int[];
    for (k,rx) in pairs(reactions(rn))
        for (spec,coef) in rx.netstoich
			push!(Js, k)
			push!(Is, smap[spec])
			push!(Vs, coef)
        end
    end
    sparse(Is,Js,Vs,numspecies(rn),numreactions(rn))
end
function netstoichmat(::Type{Matrix{Int}},rn::ReactionSystem; smap=speciesmap(rn))
    nmat = zeros(Int,numspecies(rn),numreactions(rn))
    for (k,rx) in pairs(reactions(rn))
        for (spec,coef) in rx.netstoich
            nmat[smap[spec],k] = coef
        end
    end
    nmat
end
function netstoichmat(rn::ReactionSystem; sparse=false, smap=speciesmap(rn))
    isempty(get_systems(rn)) || error("netstoichmat does not currently support subsystems.")
	sparse ? netstoichmat(SparseMatrixCSC{Int,Int}, rn; smap=smap) : netstoichmat(Matrix{Int}, rn; smap=smap)
end

"""
    setdefaults!(rn, newdefs)

Sets the default (initial) values of parameters and species in the
`ReactionSystem`, `rn`.

For example,
```julia
sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end β ν
setdefaults!(sir, [:S => 999.0, :I => 1.0, :R => 1.0, :β => 1e-4, :ν => .01])

# or
@parameter β ν
@variables t S(t) I(t) R(t)
setdefaults!(sir, [S => 999.0, I => 1.0, R => 0.0, β => 1e-4, ν => .01])
```
gives initial/default values to each of `S`, `I` and `β`

Notes:
- Can not be used to set default values for species, variables or parameters of
  subsystems or constraint systems. Either set defaults for those systems
  directly, or [`flatten`](@ref) to collate them into one system before setting
  defaults.
- Defaults can be specified in any iterable container of symbols to value pairs
  or symbolics to value pairs.
"""
function setdefaults!(rn, newdefs)
    defs = eltype(newdefs) <: Pair{Symbol} ? symmap_to_varmap(rn,newdefs) : newdefs
    rndefs = MT.get_defaults(rn)
    for (var,val) in defs
        rndefs[var] = value(val)
    end
    nothing
end

function __unpacksys(rn)
    ex = :(begin end)
    for key in keys(get_var_to_name(rn))
        var = MT.getproperty(rn, key, namespace=false)
        push!(ex.args, :($key = $var))
    end
    ex
end

"""
    @unpacksys sys::ModelingToolkit.AbstractSystem

Loads all species, variables, parameters, and observables defined in `sys` as
variables within the calling module.

For example,
```julia
sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end β ν
@unpacksys sir
```
will load the symbolic variables, `S`, `I`, `R`, `ν` and `β`.

Notes:
- Can not be used to load species, variables, or parameters of subsystems or
  constraints. Either call `@unpacksys` on those systems directly, or
  [`flatten`](@ref) to collate them into one system before calling.
- Note that this places symbolic variables within the calling module's scope, so
  calling from a function defined in a script or the REPL will still result in
  the symbolic variables being defined in the `Main` module.
"""
macro unpacksys(rn)
    quote
    ex = Catalyst.__unpacksys($(esc(rn)))
    Base.eval($(__module__), ex)
    end
end

# convert symbol of the form :sys.a.b.c to a symbolic a.b.c
function _symbol_to_var(sys, sym)
    if hasproperty(sys, sym)
        var = getproperty(sys, sym, namespace=false)
    else
        strs = split(String(sym), "₊")   # need to check if this should be split of not!!!
        if length(strs) > 1
            var = getproperty(sys, Symbol(strs[1]), namespace=false)
            for str in view(strs, 2:length(strs))
                var = getproperty(var, Symbol(str), namespace=true)
            end
        else
            throw(ArgumentError("System $(nameof(sys)): variable $sym does not exist"))
        end
    end
    var
end

"""
    symmap_to_varmap(sys, symmap)

Given a system and map of `Symbol`s to values, generates a map from
corresponding symbolic variables/parameters to the values that can be used to
pass initial conditions and parameter mappings.

For example,
```julia
sir = @reaction_network sir begin
    β, S + I --> 2I
    ν, I --> R
end β ν
subsys = @reaction_network subsys begin
    k, A --> B
end k
@named sys = compose(sir, [subsys])
```
gives
```
Model sys with 3 equations
States (5):
  S(t)
  I(t)
  R(t)
  subsys₊A(t)
  subsys₊B(t)
Parameters (3):
  β
  ν
  subsys₊k
```
to specify initial condition and parameter mappings from *symbols* we can use
```julia
symmap = [:S => 1.0, :I => 1.0, :R => 1.0, :subsys₊A => 1.0, :subsys₊B => 1.0]
u0map  = symmap_to_varmap(sys, symmap)
pmap   = symmap_to_varmap(sys, [:β => 1.0, :ν => 1.0, :subsys₊k => 1.0])
```
`u0map` and `pmap` can then be used as input to various problem types.

Notes:
- Any `Symbol`, `sym`, within `symmap` must be a valid field of `sys`. i.e.
  `sys.sym` must be defined.
"""
function symmap_to_varmap(sys, symmap::Tuple)
    if all(p -> p isa Pair{Symbol}, symmap)
        return ((_symbol_to_var(sys,sym) => val for (sym,val) in symmap)...,)
    else  # if not all entries map a symbol to value pass through
        return symmap
    end
end

symmap_to_varmap(sys, symmap::AbstractArray{Pair{Symbol,T}}) where {T} =
    [_symbol_to_var(sys,sym) => val for (sym,val) in symmap]

symmap_to_varmap(sys, symmap::Dict{Symbol,T}) where {T} =
    Dict(_symbol_to_var(sys,sym) => val for (sym,val) in symmap)

# don't permute any other types and let varmap_to_vars handle erroring
symmap_to_varmap(sys, symmap) = symmap
#error("symmap_to_varmap requires a Dict, AbstractArray or Tuple to map Symbols to values.")

######################## reaction complexes and reaction rates ###############################
"""
$(TYPEDEF)
One reaction complex element

# Fields
$(FIELDS)
"""
struct ReactionComplexElement{T}
    """The integer id of the species representing this element."""
    speciesid::Int
    """The stoichiometric coefficient of this species."""
    speciesstoich::T
end

"""
$(TYPEDEF)
One reaction complex.

# Fields
$(FIELDS)
"""
struct ReactionComplex{V<:Integer} <: AbstractVector{ReactionComplexElement{V}}
    """The integer ids of all species participating in this complex."""
    speciesids::Vector{Int}
    """The stoichiometric coefficients of all species participating in this complex."""
    speciesstoichs::Vector{V}
end

function (==)(a::ReactionComplex{V},b::ReactionComplex{V}) where {V <: Integer}
    (a.speciesids == b.speciesids) &&
    (a.speciesstoichs == b.speciesstoichs)
end
hash(rc::ReactionComplex,h::UInt) = Base.hash(rc.speciesids,Base.hash(rc.speciesstoichs,h))
Base.size(rc::ReactionComplex) = size(rc.speciesids)
Base.length(rc::ReactionComplex) = length(rc.speciesids)
Base.getindex(rc::ReactionComplex, i...) =
        ReactionComplexElement(getindex(rc.speciesids, i...), getindex(rc.speciesstoichs, i...))
Base.setindex!(rc::ReactionComplex, t::ReactionComplexElement, i...) =
    (setindex!(rc.speciesids, t.speciesid, i...); setindex!(rc.speciesstoichs, t.speciesstoich, i...); rc)
Base.isless(a::ReactionComplexElement, b::ReactionComplexElement) = isless(a.speciesid, b.speciesid)
Base.Sort.defalg(::ReactionComplex{T}) where {T <: Integer} = Base.DEFAULT_UNSTABLE

"""
    reactioncomplexmap(rn::ReactionSystem; smap=speciesmap(rn))

Find each [`ReactionComplex`](@ref) within the specified system, constructing a
mapping from the complex to vectors that indicate which reactions it appears in
as substrates and products.

Notes:
- Each [`ReactionComplex`](@ref) is mapped to a vector of pairs, with each pair
  having the form `reactionidx => ± 1`, where `-1` indicates the complex appears
  as a substrate and `+1` as a product in the reaction with integer label
  `reactionidx`.
"""
function reactioncomplexmap(rn::ReactionSystem; smap=speciesmap(rn))
    isempty(get_systems(rn)) || error("reactioncomplexmap does not currently support subsystems.")
    rxs = reactions(rn)
    numreactions(rn) > 0 || error("There must be at least one reaction to find reaction complexes.")
    complextorxsmap = OrderedDict{ReactionComplex{eltype(rxs[1].substoich)},Vector{Pair{Int,Int}}}()
    for (i,rx) in enumerate(rxs)
        reactantids = isempty(rx.substrates) ? Vector{Int}() : [smap[sub] for sub in rx.substrates]
        subrc = sort!(ReactionComplex(reactantids, copy(rx.substoich)))
        if haskey(complextorxsmap, subrc)
            push!(complextorxsmap[subrc], i => -1)
        else
            complextorxsmap[subrc] = [i => -1]
        end

        productids = isempty(rx.products) ? Vector{Int}() : [smap[prod] for prod in rx.products]
        prodrc = sort!(ReactionComplex(productids, copy(rx.prodstoich)))
        if haskey(complextorxsmap, prodrc)
            push!(complextorxsmap[prodrc], i => 1)
        else
            complextorxsmap[prodrc] = [i => 1]
        end
    end
    complextorxsmap
end

function reactioncomplexes(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem, smap, complextorxsmap)
    complexes = collect(keys(complextorxsmap))
    Is=Int[];  Js=Int[];  Vs=Int[];
	for (i,c) in enumerate(complexes)
        for (j,σ) in complextorxsmap[c]
			push!(Is, i)
			push!(Js, j)
			push!(Vs, σ)
        end
    end
	B = sparse(Is,Js,Vs,length(complexes),numreactions(rn))
    complexes,B
end
function reactioncomplexes(::Type{Matrix{Int}}, rn::ReactionSystem, smap, complextorxsmap)
    complexes = collect(keys(complextorxsmap))
    B = zeros(Int, length(complexes), numreactions(rn));
    for (i,c) in enumerate(complexes)
        for (j,σ) in complextorxsmap[c]
            B[i,j] = σ
        end
    end
    complexes,B
end

@doc raw"""
    reactioncomplexes(network::ReactionSystem; sparse=false, smap=speciesmap(rn),
                      complextorxsmap=reactioncomplexmap(rn; smap=smap))

Calculate the reaction complexes and complex incidence matrix for the given
[`ReactionSystem`](@ref).

Notes:
- returns a pair of a vector of [`ReactionComplex`](@ref)s and the complex
  incidence matrix.
- An empty [`ReactionComplex`](@ref) denotes the null (∅) state (from reactions
  like ∅ -> A or A -> ∅).
- The complex incidence matrix, ``B``, is number of complexes by number of
  reactions with
```math
B_{i j} = \begin{cases}
-1, &\text{if the i'th complex is the substrate of the j'th reaction},\\
1, &\text{if the i'th complex is the product of the j'th reaction},\\
0, &\text{otherwise.}
\end{cases}
```
- Set sparse=true for a sparse matrix representation of the incidence matrix
"""
function reactioncomplexes(rn::ReactionSystem; sparse=false, smap=speciesmap(rn),
                           complextorxsmap=reactioncomplexmap(rn; smap=smap))
    isempty(get_systems(rn)) || error("reactioncomplexes does not currently support subsystems.")
	sparse ? reactioncomplexes(SparseMatrixCSC{Int,Int}, rn, smap, complextorxsmap) :
             reactioncomplexes(Matrix{Int}, rn, smap, complextorxsmap)
end

function complexstoichmat(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem, rcs)
    Is=Int[];  Js=Int[];  Vs=Int[];
    for (i,rc) in enumerate(rcs)
        for rcel in rc
			push!(Is, rcel.speciesid)
			push!(Js, i)
			push!(Vs, rcel.speciesstoich)
        end
    end
    Z = sparse(Is,Js,Vs, numspecies(rn), length(rcs))
end
function complexstoichmat(::Type{Matrix{Int}}, rn::ReactionSystem, rcs)
    Z=zeros(Int, numspecies(rn), length(rcs))
    for (i,rc) in enumerate(rcs)
        for rcel in rc
            Z[rcel.speciesid,i] = rcel.speciesstoich
        end
    end
    Z
end

"""
    complexstoichmat(network::ReactionSystem; sparse=false, rcs=keys(reactioncomplexmap(rn)))

Given a [`ReactionSystem`](@ref) and vector of reaction complexes, return a
matrix with positive entries of size number of species by number of complexes,
where the non-zero positive entries in the kth column denote stoichiometric
coefficients of the species participating in the kth reaction complex.

Notes:
- `rcs` correspond to an iterable of the `ReactionComplexes`, i.e.
  `rcs=keys(reactioncomplexmap(rn))` or `reactioncomplexes(rn)[1]`.
- Set sparse=true for a sparse matrix representation
"""
function complexstoichmat(rn::ReactionSystem; sparse=false, rcs=keys(reactioncomplexmap(rn)))
    isempty(get_systems(rn)) || error("complexstoichmat does not currently support subsystems.")
    sparse ? complexstoichmat(SparseMatrixCSC{Int,Int}, rn, rcs) :
             complexstoichmat(Matrix{Int}, rn, rcs)
end

function complexoutgoingmat(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem, B)
    n = size(B,2)
	rows = rowvals(B)
	vals = nonzeros(B)
    Is = Int[]; Js = Int[]; Vs = Int[]
    sizehint!(Is, div(length(vals),2))
    sizehint!(Js, div(length(vals),2))
    sizehint!(Vs, div(length(vals),2))
	for j = 1:n
	   for i in nzrange(B, j)
	      if vals[i] != one(eltype(vals))
            push!(Is, rows[i])
            push!(Js, j)
            push!(Vs, vals[i])
          end
	   end
	end
    sparse(Is,Js,Vs,size(B,1),size(B,2))
end
function complexoutgoingmat(::Type{Matrix{Int}}, rn::ReactionSystem, B)
    Δ = copy(B)
    for (I,b) in pairs(Δ)
        (b == 1) && (Δ[I] = 0)
    end
    Δ
end

@doc raw"""
    complexoutgoingmat(network; sparse=false, B=reactioncomplexes(rn)[2])

Given a [`ReactionSystem`](@ref) and complex incidence matrix, ``B``, return a
matrix of size num of complexes by num of reactions that identifies substrate
complexes.

Notes:
- The complex outgoing matrix, ``\Delta``, is defined by
```math
\Delta_{i j} = \begin{cases}
    = 0,    &\text{if } B_{i j} = 1, \\
    = B_{i j}, &\text{otherwise.}
\end{cases}
```
- Set sparse=true for a sparse matrix representation
"""
function complexoutgoingmat(rn::ReactionSystem; sparse=false, B=reactioncomplexes(rn,sparse=sparse)[2])
    isempty(get_systems(rn)) || error("complexoutgoingmat does not currently support subsystems.")
    sparse ? complexoutgoingmat(SparseMatrixCSC{Int,Int}, rn, B) :
             complexoutgoingmat(Matrix{Int}, rn, B)
end

"""
    incidencematgraph(incidencemat)

Given an incidence matrix of a reaction-network, construct a directed simple
graph where nodes correspond to reaction complexes and directed edges to
reactions converting between two complexes.

For example,
```julia
sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end β ν
rcs,incidencemat = reactioncomplexes(sir)
incidencegraph   = incidencematgraph(incidencemat)
```
"""
function incidencematgraph(incidencemat::Matrix{Int})
    @assert all(∈([-1,0,1]) ,incidencemat)
    n = size(incidencemat,1)  # no. of nodes/complexes
    graph = Graphs.DiGraph(n)
    for col in eachcol(incidencemat)
        src = 0; dst = 0;
        for i in eachindex(col)
                (col[i] == -1) && (src = i)
                (col[i] == 1) && (dst = i)
                (src != 0) && (dst != 0) && break
        end
        Graphs.add_edge!(graph, src, dst)
    end
    return graph
end
function incidencematgraph(incidencemat::SparseMatrixCSC{Int,Int})
    @assert all(∈([-1,0,1]) ,incidencemat)
    m,n = size(incidencemat)
    graph = Graphs.DiGraph(m)
    rows = rowvals(incidencemat)
    vals = nonzeros(incidencemat)
    for j = 1:n
        inds=nzrange(incidencemat, j)
        row = rows[inds];
        val = vals[inds];
        if val[1] == -1
            Graphs.add_edge!(graph, row[1], row[2])
        else
            Graphs.add_edge!(graph, row[2], row[1])
        end
    end
    return graph
end

"""
    linkageclasses(incidencegraph)

Given the incidence graph of a reaction network, return a vector of the
connected components of the graph (i.e. sub-groups of reaction complexes that
are connected in the incidence graph).

For example, continuing the example from [`incidencematgraph`](@ref)
```julia
julia> linkageclasses(incidencegraph)
2-element Vector{Vector{Int64}}:
 [1, 2]
 [3, 4]
```
"""
function linkageclasses(incidencegraph)
    Graphs.connected_components(incidencegraph)
end

@doc raw"""
    deficiency(netstoich_mat, incidence_graph, linkage_classes)

Calculate the deficiency of a reaction network.

Here the deficiency, ``\delta``, of a network with ``n`` reaction complexes,
``\ell`` linkage classes and a rank ``s`` stoichiometric matrix is

```math
\delta = n - \ell - s
```

For example,
```julia
sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end β ν
rcs,incidencemat = reactioncomplexes(sir)
incidence_graph  = incidencematgraph(incidencemat)
linkage_classes   = linkageclasses(incidence_graph)
netstoich_mat    = netstoichmat(sir)
δ = deficiency(netstoich_mat, incidence_graph, linkage_classes)
```
"""
function deficiency(ns, ig, lc)
    Graphs.nv(ig) - length(lc) - AA.rank(AA.matrix(AA.ZZ,ns))
end

function subnetworkmapping(linkageclass, allrxs, complextorxmap, p)
    rxinds  = sort!(collect(Set(rxidx for rcidx in linkageclass for rxidx in complextorxmap[rcidx])))
    rxs     = allrxs[rxinds]
    specset = Set(substrate for rx in rxs for substrate in rx.substrates)
    for rx in rxs
        for product in rx.products
            push!(specset, product)
        end
    end
    specs = collect(specset)
    newps = Vector{eltype(p)}()
    for rx in rxs
        Symbolics.get_variables!(newps, rx.rate, p)
    end
    rxs, specs, newps   # reactions and species involved in reactions of subnetwork
end

"""
    subnetworks(network, linkage_classes ; rxs = reactions(network),
                  complextorxmap = collect(values(reactioncomplexmap(network))),
                  p = parameters(network))

Find subnetworks corresponding to the each linkage class of reaction network

For example, continuing the example from [`deficiency`](@ref)
```julia
   subnets = subnetworks(sir, linkage_classes)
```
"""
function subnetworks(rs::ReactionSystem, lcs::AbstractVector;
                  rxs = reactions(rs),
                  complextorxmap = [map(first,rcmap) for rcmap in values(reactioncomplexmap(rs))],
                  p = parameters(rs))
    isempty(get_systems(rs)) || error("subnetworks does not currently support subsystems.")
    t = get_iv(rs)
    subreac = Vector{ReactionSystem}()
    for i in 1:length(lcs)
        reacs,specs,newps = subnetworkmapping(lcs[i], rxs, complextorxmap, p)
        newname = Symbol(nameof(rs), "_", i)
        push!(subreac, ReactionSystem(reacs, t, specs, newps; name=newname))
    end
    subreac
end

"""
    linkagedeficiencies(subnetworks::AbstractVector, linkage_classes::AbstractVector)

Calculates the deficiency of each sub-reaction network defined by a collection
of linkage_classes.

For example, continuing the example from [`deficiency`](@ref)
```julia
subnets = subnetworks(sir, linkage_classes)
linkage_deficiencies = linkagedeficiency(subnets, linkage_classes)
```
"""
function linkagedeficiencies(subnets, lcs)
    δ = zeros(Int,length(lcs))
    for (i,subnet) in enumerate(subnets)
        ns_sub = netstoichmat(subnet)
        δ[i] = length(lcs[i]) - 1 - AA.rank(AA.matrix(AA.ZZ, ns_sub))
    end
    δ
end

"""
    isreversible(incidencegraph)

Given an incidence graph of the reaction network, returns if the network is reversible or not.
For example, continuing the example from [`linkagedeficiencies`](@ref)
```julia
isreversible(incidence_graph)
```
"""
function isreversible(ig::Graphs.SimpleDiGraph)
    Graphs.reverse(ig) == ig
end

"""
    isweaklyreversible(subnetworks)

Given the subnetworks corresponding to the each linkage class of reaction network,
determines if the reaction network is weakly reversible or not.
For example, continuing the example from [`isreversible`](@ref)
```julia
isweaklyreversible(subnets)
```
"""
function isweaklyreversible(subnets::Vector{ReactionSystem})
    igs = [incidencematgraph(reactioncomplexes(subrs)[2]) for subrs in subnets]
    all([Graphs.is_strongly_connected(ig) for ig in igs])
end

################################################################################################
######################## conservation laws ###############################

"""
    conservationlaws(netstoichmat::AbstractMatrix)::Matrix

Given the net stoichiometry matrix of a reaction system, computes a matrix of
conservation laws, each represented as a row in the output.
"""
function conservationlaws(nsm::AbstractMatrix)

    # We basically have to compute the left null space of the matrix
    # over the integers. We do this using AbstractAlgebra's integer (ZZ) interface.
    N   = AA.nullspace(AA.matrix(AA.ZZ, nsm'))[2]

    # to save allocations we manually take the adjoint when converting back
    # to a Julia integer matrix from the Nemo matrix.
    ret = [convert(Int,N[i,j]) for j=1:size(N,2), i=1:size(N,1)]

    # If all coefficients for a conservation law are negative
    # we might as well flip them to become positive
    for retrow in eachrow(ret)
        all(r -> r <= 0, retrow) && (retrow .*= -1)
    end

    ret
end

"""
    conservedquantities(state, cons_laws)

Compute conserved quantities for a system with the given conservation laws.
"""
conservedquantities(state, cons_laws) = cons_laws * state

######################## reaction network operators #######################

"""
    ==(rx1::Reaction, rx2::Reaction)

Tests whether two [`Reaction`](@ref)s are identical.

Notes:
- Ignores the order in which stoichiometry components are listed.
- *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be
    considered different than `(A+1)^2`.
"""
function (==)(rx1::Reaction, rx2::Reaction)
    isequal(rx1.rate, rx2.rate) || return false
    issetequal(zip(rx1.substrates,rx1.substoich), zip(rx2.substrates,rx2.substoich)) || return false
    issetequal(zip(rx1.products,rx1.prodstoich), zip(rx2.products,rx2.prodstoich)) || return false
    issetequal(rx1.netstoich, rx2.netstoich) || return false
    rx1.only_use_rate == rx2.only_use_rate
end

function hash(rx::Reaction, h::UInt)
    h = Base.hash(rx.rate, h)
    h = Base.hash(rx.substrates, h)
    h = Base.hash(rx.products, h)
    h = Base.hash(rx.prodstoich, h)
    h = Base.hash(rx.substoich, h)
    h = Base.hash(rx.netstoich, h)
    Base.hash(rx.only_use_rate, h)
end

"""
    isequal_ignore_names(rn1::ReactionSystem, rn2::ReactionSystem)

Tests whether the underlying species, parameters and reactions are the same in
the two [`ReactionSystem`](@ref)s. Ignores the names of the systems in testing
equality.

Notes:
- *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be
    considered different than `(A+1)^2`.
- Does not include `defaults` in determining equality.
"""
function isequal_ignore_names(rn1::ReactionSystem, rn2::ReactionSystem)
    isequal(get_iv(rn1), get_iv(rn2)) || return false
    issetequal(get_states(rn1), get_states(rn2)) || return false
    issetequal(get_ps(rn1), get_ps(rn2)) || return false
    issetequal(MT.get_observed(rn1), MT.get_observed(rn2)) || return false
    issetequal(get_eqs(rn1), get_eqs(rn2)) || return false
    (get_constraints(rn1) == get_constraints(rn2)) || return false

    # subsystems
    (length(get_systems(rn1)) == length(get_systems(rn2))) || return false
    issetequal(get_systems(rn1), get_systems(rn2)) || return false

    true
end

"""
    ==(rn1::ReactionSystem, rn2::ReactionSystem)

Tests whether the underlying species, parameters and reactions are the same in
the two [`ReactionSystem`](@ref)s. Requires the systems to have the same names
too.

Notes:
- *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be
    considered different than `(A+1)^2`.
- Does not include `defaults` in determining equality.
"""
function (==)(rn1::ReactionSystem, rn2::ReactionSystem)
    (nameof(rn1) == nameof(rn2)) || return false
    isequal_ignore_names(rn1,rn2)
end

######################## functions to extend a network ####################

"""
    make_empty_network(; iv=DEFAULT_IV, name=gensym(:ReactionSystem))

Construct an empty [`ReactionSystem`](@ref). `iv` is the independent variable,
usually time, and `name` is the name to give the `ReactionSystem`.
"""
function make_empty_network(; iv=DEFAULT_IV, name=gensym(:ReactionSystem))
    ReactionSystem(Reaction[], iv, [], []; name=name)
end

"""
    addspecies!(network::ReactionSystem, s::Symbolic; disablechecks=false)

Given a [`ReactionSystem`](@ref), add the species corresponding to the variable
`s` to the network (if it is not already defined). Returns the integer id of the
species within the system.

Notes:
- `disablechecks` will disable checking for whether the passed in variable is
  already defined, which is useful when adding many new variables to the system.
  *Do not disable checks* unless you are sure the passed in variable is a new
  variable, as this will potentially leave the system in an undefined state.
"""
function addspecies!(network::ReactionSystem, s::Symbolic; disablechecks=false)

    # we don't check subsystems since we will add it to the top-level system...
    curidx = disablechecks ? nothing : findfirst(S -> isequal(S, s), get_states(network))
    if curidx === nothing
        push!(get_states(network), s)
        MT.process_variables!(get_var_to_name(network), get_defaults(network), [s])
        return length(get_states(network))
    else
        return curidx
    end
end

"""
    addspecies!(network::ReactionSystem, s::Num; disablechecks=false)

Given a [`ReactionSystem`](@ref), add the species corresponding to the
variable `s` to the network (if it is not already defined). Returns the
integer id of the species within the system.

- `disablechecks` will disable checking for whether the passed in variable is
  already defined, which is useful when adding many new variables to the system.
  *Do not disable checks* unless you are sure the passed in variable is a new
  variable, as this will potentially leave the system in an undefined state.
"""
function addspecies!(network::ReactionSystem, s::Num; disablechecks=false)
    addspecies!(network, value(s), disablechecks=disablechecks)
end

"""
    reorder_states!(rn, neworder)

Given a [`ReactionSystem`](@ref) and a vector `neworder`, orders the states of `rn` accordingly to `neworder`.

Notes:
- Currently only supports `ReactionSystem`s without constraints or subsystems.
"""
function reorder_states!(rn, neworder)
   (get_constraints(rn) === nothing) && isempty(get_systems(rn)) ||
           error("Reordering of states is only supported for systems without constraints or subsystems.")
    permute!(get_states(rn), neworder)
end

"""
    addparam!(network::ReactionSystem, p::Symbolic; disablechecks=false)

Given a [`ReactionSystem`](@ref), add the parameter corresponding to the
variable `p` to the network (if it is not already defined). Returns the integer
id of the parameter within the system.

- `disablechecks` will disable checking for whether the passed in variable is
  already defined, which is useful when adding many new variables to the system.
  *Do not disable checks* unless you are sure the passed in variable is a new
  variable, as this will potentially leave the system in an undefined state.
"""
function addparam!(network::ReactionSystem, p::Symbolic; disablechecks=false)
    # we don't check subsystems since we will add it to the top-level system...
    if istree(p) && !(operation(p) isa Sym)
        error("If the passed in parameter is an expression, it must correspond to an underlying Variable.")
    end
    curidx = disablechecks ? nothing : findfirst(S -> isequal(S, p), get_ps(network))
    if curidx === nothing
        push!(get_ps(network), p)
        MT.process_variables!(get_var_to_name(network), get_defaults(network), [p])
        return length(get_ps(network))
    else
        return curidx
    end
end

"""
    addparam!(network::ReactionSystem, p::Num; disablechecks=false)

Given a [`ReactionSystem`](@ref), add the parameter corresponding to the
variable `p` to the network (if it is not already defined). Returns the
integer id of the parameter within the system.

- `disablechecks` will disable checking for whether the passed in variable is
  already defined, which is useful when adding many new variables to the system.
  *Do not disable checks* unless you are sure the passed in variable is a new
  variable, as this will potentially leave the system in an undefined state.
"""
function addparam!(network::ReactionSystem, p::Num; disablechecks=false)
    addparam!(network, value(p); disablechecks=disablechecks)
end

"""
    addreaction!(network::ReactionSystem, rx::Reaction)

Add the passed in reaction to the [`ReactionSystem`](@ref). Returns the
integer id of `rx` in the list of `Reaction`s within `network`.

Notes:
- Any new species or parameters used in `rx` should be separately added to
    `network` using [`addspecies!`](@ref) and [`addparam!`](@ref).
"""
function addreaction!(network::ReactionSystem, rx::Reaction)
    push!(get_eqs(network), rx)
    length(get_eqs(network))
end

"""
    merge!(network1::ReactionSystem, network2::ReactionSystem)

Merge `network2` into `network1`.

Notes:
- Duplicate reactions between the two networks are not filtered out.
- [`Reaction`](@ref)s are not deepcopied to minimize allocations, so both
  networks will share underlying data arrays.
- Subsystems are not deepcopied between the two networks and will hence be
  shared.
- Returns `network1`.
"""
function Base.merge!(network1::ReactionSystem, network2::ReactionSystem)
    ((get_constraints(network1) === nothing) && (get_constraints(network2) === nothing)) ||
        error("merge! does not currently support ReactionSystems with constraints, consider ModelingToolkit.extend instead.")
    isequal(get_iv(network1), get_iv(network2)) ||
        error("Reaction networks must have the same independent variable to be mergable.")
    append!(get_eqs(network1), get_eqs(network2))
    union!(get_states(network1), get_states(network2))
    union!(get_ps(network1), get_ps(network2))
    append!(get_observed(network1), get_observed(network2))
    append!(get_systems(network1), get_systems(network2))
    merge!(get_defaults(network1), get_defaults(network2))
    network1
end

###############################   units   #####################################

"""
    validate(rx::Reaction; info::String = "")

Check that all substrates and products within the given [`Reaction`](@ref) have
the same units, and that the units of the reaction's rate expression are
internally consistent (i.e. if the rate involves sums, each term in the sum has
the same units).

"""
function validate(rx::Reaction; info::String = "")
    validated = MT._validate([rx.rate], [string(rx, ": rate")], info = info)

    subunits = isempty(rx.substrates) ? nothing : get_unit(rx.substrates[1])
    for i in 2:length(rx.substrates)
        if get_unit(rx.substrates[i]) != subunits
            validated = false
            @warn(string("In ", rx, " the substrates have differing units."))
        end
    end

    produnits = isempty(rx.products) ? nothing : get_unit(rx.products[1])
    for i in 2:length(rx.products)
        if get_unit(rx.products[i]) != produnits
            validated = false
            @warn(string("In ", rx, " the products have differing units."))
        end
    end

    if (subunits !== nothing) && (produnits !== nothing) && (subunits != produnits)
        validated = false
        @warn(string("in ", rx, " the substrate units are not consistent with the product units."))
    end

    validated
end

"""
    validate(rs::ReactionSystem, info::String="")

Check that all species in the [`ReactionSystem`](@ref) have the same units, and
that the rate laws of all reactions reduce to units of (species units) / (time
units).

Notes:
- Does not check subsystems too.
"""
function validate(rs::ReactionSystem, info::String="")
    specs = get_states(rs)

    # if there are no species we don't check units on the system
    isempty(specs) && return true

    specunits = get_unit(specs[1])
    validated = true
    for spec in specs
        if get_unit(spec) != specunits
            validated = false
            @warn(string("Species are expected to have units of ", specunits, " however, species ", spec, " has units ", get_unit(spec), "."))
        end
    end
    timeunits = get_unit(get_iv(rs))

    # no units for species, time or parameters then assume validated
    (specunits in (MT.unitless,nothing)) && (timeunits in (MT.unitless,nothing)) &&
        MT.all_dimensionless(get_ps(rs)) && return true

    rateunits = specunits / timeunits
    for rx in get_eqs(rs)
        rxunits = get_unit(rx.rate)
        for (i,sub) in enumerate(rx.substrates)
            rxunits *= get_unit(sub) ^ rx.substoich[i]
        end

        if rxunits != rateunits
            validated = false
            @warn(string("Reaction rate laws are expected to have units of ", rateunits, " however, ", rx, " has units of ", rxunits, "."))
        end
    end

    validated
end


text:
module FIRLSFilterDesign
"""
https://cnx.org/contents/6x7LNQOp@7/Linear-Phase-Fir-Filter-Design-By-Least-Squares
https://www.dsprelated.com/Plishowarticle/808.php
https://eeweb.engineering.nyu.edu/iselesni/EL713/zoom/linphase.pdf
"""
# using LinearAlgebra
using LinearAlgebra
import LinearAlgebra: I, Diagonal, UniformScaling

include("utils.jl")

export firls_design

abstract type FIR end
struct FIR_I <: FIR end
struct FIR_II <: FIR end
struct FIR_III <: FIR end
struct FIR_IV <: FIR end

@doc """
    infer_fir_type(is_odd, is_antisymmetric)

Determines the type of FIR filter to be designed, based on:
* Whether the number of filter coefficients is odd (`is_odd`)
* Whether the filter should be antisymmetric (`is_antisymmetric`)

The result is a `fir_type`, which can be:
1. `FIR_I`, when filter length is odd and the filter is not antisymmetric
2. `FIR_II`, when filter length is even and the filter is not antisymmetric
3. `FIR_III`, when filter length is odd and the filter is antisymmetric
4. `FIR_IV`, when filter length is even and the filter is antisymmetric

# Arguments
- `is_odd::Bool`
- `is_antisymmetric::Bool`

# Outputs
- `fir_type::FIR`
"""
function infer_fir_type(is_odd, is_antisymmetric)
    if !is_antisymmetric
        fir_type = is_odd ? FIR_I : FIR_II
    else
        fir_type = is_odd ? FIR_III : FIR_IV
    end
    fir_type()
end

function validate_inputs(filter_order, bands_D, D, fs)
    @assert filter_order >= 0 "Filter order cannot be negative."
    @assert fs > 0 "Sampling frequency should be larger than 0 Hz"
    validate_bands_D(bands_D, D, fs)
end
function validate_inputs(filter_order, bands_DW, D, W, fs)
    @assert filter_order >= 0 "Filter order cannot be negative."
    @assert fs > 0 "Sampling frequency should be larger than 0 Hz"
    validate_bands_D(bands_DW, D, fs)
    validate_bands_W(bands_DW, W, fs)
end
function validate_inputs(filter_order, bands_D, D, bands_W, W, fs)
    @assert filter_order >= 0 "Filter order cannot be negative."
    @assert fs > 0 "Sampling frequency should be larger than 0 Hz"
    validate_bands_D(bands_D, D, fs)
    validate_bands_W(bands_W, W, fs)
end

function validate_bands_D(bands_D, D, fs)
    validate_bands(bands_D, fs)
    @assert size(bands_D,2) == size(D,2) == 2 "Frequency bands and desired response should be N x 2 matrices"
    @assert size(bands_D,1) == size(D,1) "Frequency bands and desired response should be N x 2 matrices"
end

function validate_bands_W(bands_W, W, fs)
    validate_bands(bands_W, fs)
    @assert size(bands_W,2) == size(W,2) == 2 "Frequency band and weight matrix should be N x 2."
    @assert size(bands_W,1) == size(W,1) "Number of frequency bands should match number of weights."
end

function validate_bands(fbands::Matrix{T}, fs) where T
    @assert fbands[1] == 0 "Frequency bands should start at 0."
    @assert fbands[end] == fs * (1//2) "Frequency bands should end at fs/2"
    @views @assert all((fbands[2:end,1] .- fbands[1:end-1,2]) .== zero(T)) "Frequency bands should cover entire [0,fs/2] interval, without gaps or overlaps."
end

@doc """
    get_flength_M(filter_order)

Determines the length of the FIR filter and the number of amplitude coefficients needed, based on its order.

# Arguments
- `filter_order::Integer` : the order of the FIR filter.

# Outputs
- `filter_length::Integer` : the number of filter coefficients
- `M::Integer` : the number of unique amplitude coefficients needed to form the filter is equal to ``M+1``
"""
function get_flength_M(filter_order)
    filter_length = filter_order + 1
    M = isodd(filter_length) ? (filter_length-1)÷2 : filter_length÷2-1
    return filter_length, M
end

to_matrix_simple(v::Vector) = hcat(v,v)
to_matrix_simple(m::Matrix) = m
knotpoints_to_matrix(v::Vector) = @views hcat(v[1:end-1], v[2:end])

@doc """
    firls_design(filter_order::Integer, bands_DW::Matrix, D::Matrix, W::Matrix, antisymmetric::Bool; fs::Real = 1, solver::Function = \\)

Designs a linear-phase FIR filter.

# Arguments
- `filter_order::Integer`   : the order of the FIR filter.
- `bands_DW::Matrix`        : a matrix of size (N,2) which contains rows of sequential frequency bands, spanning [0, fs/2].
- `D::Matrix`               : a matrix of size (N,2) which contains rows of amplitude responses for the frequency bands in `bands_DW`. The first and second columns indicate the amplitude response at the lower and upper bound of the frequency bands, interpolated linearly in between.
- `W::Matrix`               : a matrix of size (N,2) which contains rows of weighting coefficients for the frequency bands in `bands_DW`. The first and second columns indicate the weighting at the lower and upper bound of the frequency bands, interpolated linearly in between.
- `antisymmetric::Bool`     : a Boolean that signifies whether the filter coefficients will be anti-symmetric, as used in type III and IV FIR filters.
- `fs::Real`                : the sampling frequency.
- `solver::Function`        : the function that is called to solve the equation ``Qa = b``, with the function call: `solver(Q,b)` which returns `a`.

# Outputs
- `h` : a vector of linear-phase FIR filter coefficients.
"""
function firls_design(filter_order::Integer, bands_DW::Matrix, D::Matrix, W::Matrix, antisymmetric::Bool; fs::Real = 1, solver::Function = \)
    validate_inputs(filter_order, bands_DW, D, W, fs)
    filter_length, M = get_flength_M(filter_order)
    fir_type = infer_fir_type(isodd(filter_length), antisymmetric)
    a = calc_amplitude_coeff(M, bands_DW, D, W, solver, fir_type)
    h = _to_impulse_response(a, fir_type)
end
@doc """
    firls_design(filter_order::Integer, bands_DW::Matrix, D::Union{Vector,Matrix}, W::Union{Vector,Matrix}, antisymmetric::Bool; fs::Real = 1, solver::Function = \\)

# Arguments
- `filter_order::Integer`   : the order of the FIR filter.
- `bands_DW::Matrix`        : a matrix of size `(N,2)` which contains rows of sequential frequency bands, spanning [0, fs/2].
- `D::Union{Vector,Matrix}` : a matrix of size `(N,2)`, or a vector of size `(N,)`, which amplitude responses for the frequency bands in `bands_DW`. In the case of a matrix, the first and second columns indicate the amplitude response at the lower and upper bound of the frequency bands, interpolated linearly in between. In the case of a vector the elements the amplitude response is piecewise constant.
- `W::Union{Vector,Matrix}` : a matrix of size `(N,2)`, or a vector of size `(N,)`, which contains weighting function values for the frequency bands in `bands_DW`. In the case of a matrix, the first and second columns indicate the weighting function values at the lower and upper bound of the frequency bands, interpolated linearly in between. In the case of a vector the elements the weighting function is piecewise constant.
- `antisymmetric::Bool`     : a Boolean that signifies whether the filter coefficients will be anti-symmetric, as used in type III and IV FIR filters.
- `fs::Real`                : the sampling frequency.
- `solver::Function`        : the function that is called to solve the equation ``Qa = b``, with the function call: `solver(Q,b)` which returns `a`.

# Outputs
- `h` : a vector of linear-phase FIR filter coefficients.
"""
function firls_design(filter_order::Integer, bands_DW::Matrix, D::Union{Vector,Matrix}, W::Union{Vector,Matrix}, antisymmetric::Bool; fs::Real = 1, solver::Function = \)
    firls_design(filter_order, bands_DW, to_matrix_simple(D), to_matrix_simple(W), antisymmetric; fs = fs, solver = solver)
end

@doc """
    firls_design(filter_order::Integer, knotpoints_DW::Vector, D::Vector, W::Vector, antisymmetric::Bool; fs::Real = 1, solver::Function = \\)

# Arguments
- `filter_order::Integer`   : the order of the FIR filter.
- `knotpoints_DW::Vector`   : a vector of size `(N,)` which contains frequency knotpoints, spanning [0, fs/2].
- `D::Vector`               : a vector of size `(N,)` which contains amplitude response values for the frequency knotpoints in `knotpoints_DW`.
- `W::Vector`               : a vector of size `(N,)` which contains weighting function values for the frequency knotpoints in `knotpoints_DW`.
- `antisymmetric::Bool`     : a Boolean that signifies whether the filter coefficients will be anti-symmetric, as used in type III and IV FIR filters.
- `fs::Real`                : the sampling frequency.
- `solver::Function`        : the function that is called to solve the equation ``Qa = b``, with the function call: `solver(Q,b)` which returns `a`.

# Outputs
- `h` : a vector of linear-phase FIR filter coefficients.
"""
function firls_design(filter_order::Integer, knotpoints_DW::Vector, D::Vector, W::Vector, antisymmetric::Bool; fs::Real = 1, solver::Function = \)
    firls_design(filter_order, knotpoints_to_matrix(knotpoints_DW), knotpoints_to_matrix(D), knotpoints_to_matrix(W), antisymmetric; fs = fs, solver = solver)
end

@doc """
    firls_design(filter_order::Integer, bands_DW::Matrix, D::Union{Vector,Matrix}, antisymmetric::Bool; fs::Real = 1, solver::Function = \\)

# Arguments
- `filter_order::Integer`   : the order of the FIR filter.
- `bands_DW::Matrix`        : a matrix of size `(N,2)` which contains rows of sequential frequency bands, spanning [0, fs/2].
- `D::Union{Vector,Matrix}` : a matrix of size `(N,2)`, or a vector of size `(N,)`, which amplitude responses for the frequency bands in `bands_DW`. In the case of a matrix, the first and second columns indicate the amplitude response at the lower and upper bound of the frequency bands, interpolated linearly in between. In the case of a vector the elements the amplitude response is piecewise constant.
- `antisymmetric::Bool`     : a Boolean that signifies whether the filter coefficients will be anti-symmetric, as used in type III and IV FIR filters.
- `fs::Real`                : the sampling frequency.
- `solver::Function`        : the function that is called to solve the equation ``Qa = b``, with the function call: `solver(Q,b)` which returns `a`.

# Outputs
- `h` : a vector of linear-phase FIR filter coefficients.
"""
function firls_design(filter_order::Integer, bands_DW::Matrix, D::Union{Vector,Matrix}, antisymmetric::Bool; fs::Real = 1, solver::Function = \)
    D = to_matrix_simple(D)
    validate_inputs(filter_order, bands_DW, D, fs)
    filter_length, M = get_flength_M(filter_order)
    fir_type = infer_fir_type(isodd(filter_length), antisymmetric)
    a = calc_amplitude_coeff(M, bands_DW, D, solver, fir_type)
    h = _to_impulse_response(a, fir_type)
end

@doc """
    firls_design(filter_order::Integer, knotpoints_D::Vector, D::Vector, antisymmetric::Bool; fs::Real = 1, solver::Function = \\)

# Arguments
- `filter_order::Integer`   : the order of the FIR filter.
- `knotpoints_D::Vector`    : a vector of size `(N,)` which contains frequency knotpoints, spanning [0, fs/2].
- `D::Vector`               : a vector of size `(N,)` which contains amplitude response values for the frequency knotpoints in `knotpoints_D::Vector`.
- `antisymmetric::Bool`     : a Boolean that signifies whether the filter coefficients will be anti-symmetric, as used in type III and IV FIR filters.
- `fs::Real`                : the sampling frequency.
- `solver::Function`        : the function that is called to solve the equation ``Qa = b``, with the function call: `solver(Q,b)` which returns `a`.

# Outputs
- `h` : a vector of linear-phase FIR filter coefficients.
"""
function firls_design(filter_order::Integer, knotpoints_D::Vector, D::Vector, antisymmetric::Bool; fs::Real = 1, solver::Function = \)
    firls_design(filter_order, knotpoints_to_matrix(knotpoints_D), knotpoints_to_matrix(D), antisymmetric; fs = fs, solver = solver)
end

function calc_amplitude_coeff(M, bands_DW, D, W, solver, fir_type)
    Q = get_Q(M, bands_DW, W, fir_type)
    b = get_b(M, bands_DW, D, W, fir_type)
    a = solve(Q, b, solver, fir_type)
end
function calc_amplitude_coeff(M, bands_DW, D, solver, fir_type)
    Q = get_Q(M, fir_type)
    b = get_b(M, bands_DW, D, ones(size(D)), fir_type)
    a = solve(Q, b, solver, fir_type)
end

solve(Q, b, solver, fir_type) = solver(Q, b)
function solve(Q, b, solver, fir_type::FIR_III)
    a = zeros(length(b))
    if length(a) > 1
        @views a[2:end] .= solver(Q[2:end,2:end], b[2:end])
    end
    return a
end
function solve(Q::UniformScaling, b, solver, fir_type::FIR_III)
    a = zeros(length(b))
    if length(a) > 1
        @views a[2:end] .= solver(Q, b[2:end])
    end
    return a
end

"""
    get_Q(M, f, W, fir_type)

Constructs the matrix ``Q`` used in the equation ``Qa = b``, based on a set of weights.

# Arguments
- `M::Integer`      : indicator of the amount of elements needed.
- `f::Matrix`       : a matrix of size `(N,2)` which contains rows of sequential frequency bands, spanning [0, fs/2].
- `W::Matrix`       : a matrix of size `(N,2)` which contains rows of weighting coefficients for the frequency bands in `f`. The first and second columns indicate the weighting at the lower and upper bound of the frequency bands, interpolated linearly in between.
- `fir_type::FIR`   : indicates the type of FIR filter.

# Outputs
- `Q::Matrix` : the matrix ``Q`` used in the equation ``Qa = b``.
"""
function get_Q(M, f, W, fir_type)
    q = get_q(M, f, W, fir_type)
    Q1, Q2 = q_to_Q1Q2(q, M, fir_type)
    Q = Q1Q2_to_Q(Q1, Q2, fir_type)
end
"""
    get_Q(M, fir_type)

Constructs the matrix ``Q`` used in the equation ``Qa = b``, when there are no weights. Which results in ``Q`` being the identity matrix.

# Arguments
- `M::Integer`      : indicator of the amount of elements needed.
- `fir_type::FIR`   : indicates the type of FIR filter.

# Outputs
- `Q::Matrix` : the matrix ``Q`` used in the equation ``Qa = b``.
"""
get_Q(M, fir_type) = I
"""
    get_Q(M, fir_type::FIR_I)

Constructs the matrix ``Q`` used in the equation ``Qa = b``, when there are no weights and the FIR filter is of type I.

# Arguments
- `M::Integer`      : indicator of the amount of elements needed.
- `fir_type::FIR_I`   : indicates the type of FIR filter is I.

# Outputs
- `Q::Matrix` : the matrix ``Q`` used in the equation ``Qa = b``.
"""
function get_Q(M, fir_type::FIR_I)
    v_diag = fill(1., M+1)
    v_diag[1] *= 2
    Q = Diagonal(v_diag)
end

function q_to_Q1Q2(q, M, fir_type::Union{FIR_I,FIR_III})
    Q1 = @views to_toeplitz(q[1:M+1])
    Q2 = @views to_hankel(q[1:M+1], q[M+1:end])
    return Q1, Q2
end
function q_to_Q1Q2(q, M, fir_type::Union{FIR_II,FIR_IV})
    Q1 = @views to_toeplitz(q[1:M+1])
    Q2 = @views to_hankel(q[1+1:M+1+1], q[M+1+1:end])
    return Q1, Q2
end

Q1Q2_to_Q(Q1, Q2, fir_type::Union{FIR_I,FIR_II}) = Q1 .+ Q2
Q1Q2_to_Q(Q1, Q2, fir_type::Union{FIR_III,FIR_IV}) = Q1 .- Q2

"""
    get_q(M, f, W, fir_type)

Finds the vector ``q`` which is used to populate the matrix ``Q``.
...
# Arguments
- `M::Integer`      : indicator of the amount of elements needed.
- `f::Matrix`       : a matrix of size `(N,2)` which contains rows of sequential frequency bands, spanning [0, fs/2].
- `W::Matrix`       : a matrix of size `(N,2)` which contains rows of weighting coefficients for the frequency bands in `f`. The first and second columns indicate the weighting at the lower and upper bound of the frequency bands, interpolated linearly in between.
- `fir_type::FIR`   : indicates the type of FIR filter.

# Outputs
- `q_out::Vector`   : a vector of q-values that are used to fill in the Q-matrix.
...
"""
function get_q(M, f, W, fir_type)
    a, b, α, β, γ, k = constants_q(f, W)
    _αn, _βn⁻¹, _γ⁻² = copy(α), copy(β), copy(γ)
    q_out, _qn = allocate_q(M, fir_type), zeros(size(f))
    q_out[1] = qn!(_qn, k, f, a, b)
    for idx = 2:length(q_out)
        _αn .= α; _βn⁻¹ .= β; _γ⁻² .= γ
        n = idx2n_q(idx, fir_type)
        q_out[idx] = qn!(_qn, n, k, _αn, _βn⁻¹, _γ⁻²)
    end
    return q_out
end

function constants_q(f, W)
    fs = 2f[end,2]; k = 2/fs
    a = @views @. (W[:,2] - W[:,1]) / (f[:,2] - f[:,1])
    b = @views @. W[:,1] - f[:,1]*a
    α = @. π*k*f
    β = @. (1/(π*k)) * (a*f + b)
    γ = @. (1/(π*k)^2) * a
    return a, b, α, β, γ, k
end

allocate_q(M, fir_type::Union{FIR_I,FIR_III}) = q_out = zeros(2M+1)
allocate_q(M, fir_type::Union{FIR_II,FIR_IV}) = q_out = zeros(2M+1+1)

idx2n_q(idx, fir_type::Union{FIR_I,FIR_III}) = idx-1.
idx2n_q(idx, fir_type::Union{FIR_II,FIR_IV}) = idx-1.

function qn!(_qn, k, f, a, b)
    # Fallback for when n = 0
    # _qn = (af²/2 + bf)
    #_qn = 1/2 * (af² + 2bf)
    # bn = k * (_qn[:,2] - _qn[:,1])
    @. _qn = b*f
    @. _qn *= 2
    @. _qn += a*f^2
    nan2zero!(_qn)
    bn = @views k/2 * (sum(_qn[:,2]) - sum(_qn[:,1]))
end

function qn!(_qn, n, k, _αn, _βn⁻¹, _γn⁻²)
    _αn .*= n
    @. _βn⁻¹ *= 1/n
    @. _γn⁻² *= 1/n^2
    @. _qn = _βn⁻¹ * sin(_αn) + _γn⁻² * cos(_αn)
    nan2zero!(_qn)
    bn = @views k * (sum(_qn[:,2]) - sum(_qn[:,1]))
end

"""
    get_b(M, f, D, W, fir_type)

Finds the vector ``b`` used in the equation ``Qa = b``.
...
# Arguments
- `M::Integer`      : Size of the `b` vector is `M+1` .
- `f::Matrix`       : a matrix of size `(N,2)` which contains rows of sequential frequency bands, spanning the interval [0, fs/2].
- `D::Matrix`       : a matrix of size `(N,2)` which contains rows of desired frequency response values for the frequency bands in `f`. The first and second columns indicate the desired response at the lower and upper bound of the frequency bands, interpolated linearly in between.
- `W::Matrix`       : a matrix of size `(N,2)` which contains rows of weighting coefficients for the frequency bands in `f`. The first and second columns indicate the weighting at the lower and upper bound of the frequency bands, interpolated linearly in between.
- `fir_type::FIR`   : indicates the type of FIR filter.

# Outputs
- `b_out::Vector`   : a vector of size `(M+1,)`, the b-vector used in the equation ``Qa = b``.
...
"""
function get_b(M, f, D, W, fir_type)
    a, b, c, d, α, β, γ, δ, k = constants_b(f, D, W)
    _αn, _βn², _δn = copy(α), copy(β), copy(δ)
    b_out, _bn = zeros(M+1), zeros(size(f))
    b_out[1] = bn_n0!(_bn, idx2n_b(1, fir_type), k, f, a, b, c, d, _αn, _βn², γ, _δn, fir_type)
    for idx = 2:length(b_out)
        _αn .= α; _βn² .= β; _δn .= δ
        n = idx2n_b(idx, fir_type)
        b_out[idx] = bn!(_bn, n, k, _αn, _βn², γ, _δn, fir_type)
    end
    return b_out
end

@doc raw"""
    constants_b(f, D, W)

Calculates data that is reused at every evaluation of [`bn!`](@ref).

...
# Arguments
- `f::Matrix` : a matrix of size `(N,2)` which contains rows of sequential frequency bands, spanning the interval [0, fs/2].
- `D::Matrix` : a matrix of size `(N,2)` which contains rows of desired frequency response values for the frequency bands in `f`. The first and second columns indicate the desired response at the lower and upper bound of the frequency bands, interpolated linearly in between.
- `W::Matrix` : a matrix of size `(N,2)` which contains rows of weighting coefficients for the frequency bands in `f`. The first and second columns indicate the weighting at the lower and upper bound of the frequency bands, interpolated linearly in between.
...
"""
function constants_b(f, D, W)
    fs = 2f[end,2]; k = 2/fs
    a = @views @. (D[:,2] - D[:,1]) / (f[:,2] - f[:,1])
    b = @views @. D[:,1] - f[:,1]*a
    c = @views @. (W[:,2] - W[:,1]) / (f[:,2] - f[:,1])
    d = @views @. W[:,1] - f[:,1]*c
    α = @. k*π*f
    β = @. (k*π)^2 * (a*c*f^2 + (a*d+b*c)*f + b*d)
    γ = @. -2a*c
    δ = @. k*π * (2a*c*f + a*d + b*c)
    return a, b, c, d, α, β, γ, δ, k
end

@doc raw"""
    bn_n0!(_bn, n, k, f, a, b, c, d, _αn, _βn², γ, _δn, fir_type)

Dispatches to the correct function to calculate the first element of the `b` vector, based on the type of FIR filter.
Needed because for type I FIR filters the value of `n` at the first iteration is `0`.
"""
function bn_n0!(_bn, n, k, f, a, b, c, d, _αn, _βn², γ, _δn, fir_type)
    bn!(_bn, n, k, _αn, _βn², γ, _δn, fir_type)
end
function bn_n0!(_bn, n, k, f, a, b, c, d, _αn, _βn², γ, _δn, fir_type::FIR_I)
    bn!(_bn, k, f, a, b, c, d)
end

@doc raw"""
    idx2n_b(idx, fir_type::Union{FIR_I,FIR_III})

Determine the value of ``n`` based on the index in the b-vector.
For type I and III FIR filters the following holds: ``n = i - 1``.
Based on page 10 and 12 of [this](https://eeweb.engineering.nyu.edu/iselesni/EL713/zoom/linphase.pdf).

...
# Arguments
- `i::Integer` : index in the b-vector.
- `fir_type::Union{FIR_I,FIR_III}`
# Outputs
- `n::Real`
...
"""
idx2n_b(i, fir_type::Union{FIR_I,FIR_III}) = n = i-1.
@doc raw"""
    idx2n_b(idx, fir_type::Union{FIR_II,FIR_IV})

Determine the value of ``n`` based on the index in the b-vector, Where ``n`` is used in
For type II and IV FIR filters the following holds: ``n = i - 1/2``.
Based on page 11 and 13 of [this](https://eeweb.engineering.nyu.edu/iselesni/EL713/zoom/linphase.pdf).

...
# Arguments
- `i::Integer` : index in the b-vector.
- `fir_type::Union{FIR_II,FIR_IV}`

# Outputs
- `n::Real`
...
"""
idx2n_b(i, fir_type::Union{FIR_II,FIR_IV}) = n = i-0.5

@doc raw"""
    bn!(_bn, n, k, _αn, _βn², γ, _δn, fir_type)

Calculates the elements of the b-vector, which are equal to:
```math
b[i] = \frac{2}{f_s} \int_0^{f_s/2} W(f) D(f) cos(\pi \frac{2}{f_s} n f) df, \quad i = 1, 2, \cdots, M+1
```

Both ``D(f)`` and ``W(f)`` are piecewise linear functions, and ``n`` is calculated by [`idx2n_b`](@ref). Using ``k = 2/f_s``, this integral becomes:
```math
b[i] = k \sum_{j=1}^J \int_{F_{j,1}}^{F_{j,2}} g_j(f,n) df = k \sum_{j=1}^J \int_{F_{j,1}}^{F_{j,2}} \big(c_j f+d_j\big) \big(a_j f+b_j\big) \cos(\pi k n f) df, \quad i = 1, 2, \cdots, M+1
```
Where:
* ``F_{j,1}`` and ``F_{j,2}`` are the lower and upper bound of the ``j^{th}`` frequency band.
* ``a_j`` and ``b_j`` are the parameters of the linear function that describes the desired frequency response in the ``j^{th}`` frequency band.
* ``c_j`` and ``d_j`` are the parameters of the linear function that describes the error weighting function in the ``j^{th}`` frequency band.

The antiderivative of ``g(f,n)`` is equal to:
```math
G_j(f,n) = \frac{1}{\pi^3 k^3 n^3} \bigg(\sin\big(\alpha(f) n\big)\Big(\beta_j(f) n^2 + \gamma_j(f)\Big) + \delta_j(f) n \cos\big(\alpha(f) n\big)\bigg) + constants
```
Where:
*  ``\alpha(f) = \pi kf``,
* ``\beta(f) = \pi^2 k^2 \big(acf^2 + (ad+bc)f + bd\big)``,
* ``\gamma(f) = -2ac``,
* ``\delta(f) = \pi k\big(2acf + ad + bc\big)``.
Note that the subscripts have been dropped here for clarity.

...
# Arguments
- `_bn::Vector`   : a vector of size `(J,)` that is used to store the intermediate values.
- `n::Integer`    : integer denoting the current cosine mode.
- `k::Real`       : equal to `2/fs`.
- `_αn::Vector`   : a vector of size `(J,)` that holds the values of ``\alpha n``.
- `_βn²::Vector`  : a vector of size `(J,)` that holds the values of ``\beta n^2``.
- `γ::Vector`     : a vector of size `(J,)` that holds the values of ``\gamma``.
- `_δn::Vector`   : a vector of size `(J,)` that holds the values of ``\delta n``.
- `fir_type::FIR` : indicates the type of FIR filter.

# Outputs
- `bn`            : ``n^{th}`` element in the ``b`` vector.
...
"""
function bn!(_bn, n, k, _αn, _βn², γ, _δn, fir_type)
    # if antisymmetric is false: _bn = 1/(π³k³n³) * (sin(αn) * (βn² + γ) + δncos(αn))
    # if antisymmetric is true:  _bn = 1/(π³k³n³) * (-cos(αn) * (βn² + γ) + δnsin(αn))
    #                 equal to:  _bn = 1/(π³k³n³) * (sin(αn-π/2) * (βn² + γ) + δncos(αn-π/2))
    # see _update_trig_arg_b
    #  bn = k * (_bn[:,2] - _bn[:,1])
    _update_trig_arg_b!(_αn, n, fir_type)
    @. _βn² *= n^2
    @. _δn *= n
    @. _bn = sin(_αn)*(_βn² + γ) + _δn*cos(_αn)
    nan2zero!(_bn)
    bn = @views k/(π*n*k)^3 * (sum(_bn[:,2]) - sum(_bn[:,1]))
end

@doc raw"""
    bn!(_bn, k, f, a, b, c, d)

Special case for when ``n = 0``, since then the integral is simplified:
```math
g_j(f,0) = (c_j f+d_j) (a_j f+b_j) \cos(\pi k f 0) = (c_j f+d_j) (a_j f+b_j)
```
And the antiderivative becomes:
```math
G_j(f,0) = a_j c_j \frac{f^3}{3}  + (a_j d_j + b_j c_j)\frac{f^2}{2} + b_j d_j f
```

...
# Arguments
- `_bn::Vector` : a vector of size `(J,)` that is used to store the intermdediate values.
- `k::Real` : equal to ``2/f_s``.
- `f::Matrix` : a matrix of size `(J,2)` which contains rows of sequential frequency bands, spanning [0, fs/2].
- `a::Vector` : a vector of size `(J,)` with the ``a_j`` values in the equation ``a_j f + b_j`` that equates to the linear function that describes the desired frequency response in the ``j^{th}`` frequency band.
- `b::Vector` : a vector of size `(J,)` with the ``b_j`` values in the equation ``a_j f + b_j`` that equates to the linear function that describes the desired frequency response in the ``j^{th}`` frequency band.
- `c::Vector` : a vector of size `(J,)` with the ``c_j`` values in the equation ``c_j f + d_j`` that equates to the linear function that describes the error weighting function in the ``j^{th}`` frequency band.
- `d::Vector` : a vector of size `(J,)` with the ``d_j`` values in the equation ``c_j f + d_j`` that equates to the linear function that describes the error weighting function in the ``j^{th}`` frequency band.
...
"""
function bn!(_bn, k, f, a, b, c, d)
    # _bn = acf³/3 + (ad+bc)f²/2 + bdf
    # bn = k * (_bn[:,2] - _bn[:,1])
    @. _bn = b*d*f
    @. _bn *= 2
    @. _bn += (a*d+b*c)*f^2
    @. _bn *= 3/2
    @. _bn += a*c*f^3
    nan2zero!(_bn)
    bn = @views k/3 * (sum(_bn[:,2]) - sum(_bn[:,1]))
end

@doc raw"""
    _update_trig_arg_b!(_αn, n, fir_type::Union{FIR_I,FIR_II})

Updates the argument of the trigonometric functions in [`bn!`](@ref) by multiplying with the current `n`.

...
# Arguments
- `_αn::Vector` : a vector of size `(J,)`
- `n::Real`
- `fir_type{Union{FIR_I,FIR_II}}`
...
"""
_update_trig_arg_b!(_αn, n, fir_type::Union{FIR_I,FIR_II}) = _αn .*= n
@doc raw"""
    _update_trig_arg_b!(_αn, n, fir_type::Union{FIR_III,FIR_IV})

Updates the argument of the trigonometric functions in [`bn!`](@ref) by multiplying with ``n`` and subtracting ``\pi``.
The subtraction of ``\pi`` is necessary because when the filter is antisymmetric (type III and IV FIR filters), the filter response is a sum of sines instead of cosines and ``\sin(x) = \cos(x - \pi/2)`` (see page 12 and 13 of [this](https://eeweb.engineering.nyu.edu/iselesni/EL713/zoom/linphase.pdf)).

...
# Arguments
- `_αn::Vector` : a vector of size `(J,)`
- `n::Real`
- `fir_type::Union{FIR_III,FIR_IV}`
...
"""
function _update_trig_arg_b!(_αn, n, fir_type::Union{FIR_III,FIR_IV})
    _αn .*= n
    _αn .-= π/2
end

@doc raw"""
    _to_impulse_response(a, fir_type)

Creates a linear phase FIR filter based on `fir_type` and the coefficients in vector ``a``, which was obtained by solving the linear equation ``Qa = b``.

...
# Arguments
-`a`::Vector      : a vector of size `(M+1,)` with coefficients.
- `fir_type::FIR` : indicates the type of FIR filter.

# Outputs
-`h`::Vector      : a vector of size (filter_order+1,) with the filter coefficients.
...

For type I FIR filters:
```math
    h = \bigg[ a[M+1] \quad a[M] \quad \cdots \quad a[2] \quad a[1] \quad a[2] \quad \cdots \quad a[M+1] \bigg]^T
```
"""
function _to_impulse_response(a, fir_type::FIR_I)
    filter_length = 2length(a) - 1
    N_a, h = length(a), zeros(eltype(a), filter_length)
    h[1:N_a-1] .= @view(a[end:-1:2])
    h[N_a] = 2a[1]
    h[N_a+1:end] .= @view(a[2:end])
    return h
end
@doc raw"""
For type II FIR filters:
```math
    h = \bigg[ a[M+1] \quad a[M] \quad \cdots \quad a[2] \quad a[1] \quad a[1] \quad a[2] \quad \cdots \quad a[M+1] \bigg]^T
```
"""
function _to_impulse_response(a, fir_type::FIR_II)
    filter_length = 2length(a)
    N_a, h = length(a), zeros(eltype(a), filter_length)
    h[1:N_a] .= @view(a[end:-1:1])
    h[N_a+1:end] .= a
    return h
end
@doc raw"""
For type III FIR filters:
```math
    h = \bigg[ a[M+1] \quad a[M] \quad \cdots \quad a[2] \quad 0 \quad -a[2] \quad \cdots \quad -a[M+1] \bigg]^T
```
"""
function _to_impulse_response(a, fir_type::FIR_III)
    filter_length = 2length(a)-1
    N_a, h = length(a), zeros(eltype(a), filter_length)
    h[1:N_a-1] .= @view(a[end:-1:2])
    h[N_a+1:end] .-= @view(a[2:end])
    return h
end
@doc raw"""
For type IV FIR filters:
```math
    h = \bigg[ a[M+1] \quad a[M] \quad \cdots \quad a[2] \quad 0 \quad 0 \quad -a[2] \quad \cdots \quad -a[M+1] \bigg]^T
```
"""
function _to_impulse_response(a, fir_type::FIR_IV)
    filter_length = 2length(a)
    N_a, h = length(a), zeros(eltype(a), filter_length)
    h[1:N_a] .= @view(a[end:-1:1])
    h[N_a+1:end] .-= a
    return h
end

end