text:
"""
    FourierDerivativeOperator{T}

A derivative operator on a periodic grid with real scalar type `T` computing the
first derivative using a spectral Fourier expansion via real discrete Fourier
transforms.
"""
struct FourierDerivativeOperator{T<:Real, Grid, RFFT, BRFFT} <: AbstractPeriodicDerivativeOperator{T}
    jac::T
    Δx::T
    grid_compute::Grid   # N-1 nodes, including the left and excluding the right boundary
    grid_evaluate::Grid #  N  nodes, including both boundaries
    tmp::Vector{Complex{T}}
    rfft_plan::RFFT
    brfft_plan::BRFFT

    function FourierDerivativeOperator(jac::T, Δx::T, grid_compute::Grid, grid_evaluate::Grid,
                                        tmp::Vector{Complex{T}}, rfft_plan::RFFT, brfft_plan::BRFFT) where {T<:Real, Grid, RFFT, BRFFT}
        @argcheck length(brfft_plan) == length(tmp) DimensionMismatch
        @argcheck length(brfft_plan) == (length(rfft_plan)÷2)+1 DimensionMismatch
        @argcheck length(grid_compute) == length(rfft_plan) DimensionMismatch
        @argcheck length(grid_compute) == length(grid_evaluate)-1 DimensionMismatch
        @argcheck first(grid_compute) == first(grid_evaluate)
        @argcheck step(grid_compute) ≈ step(grid_evaluate)
        @argcheck last(grid_compute) < last(grid_evaluate)

        new{T, Grid, RFFT, BRFFT}(jac, Δx, grid_compute, grid_evaluate, tmp, rfft_plan, brfft_plan)
    end
end

"""
    FourierDerivativeOperator(xmin::T, xmax::T, N::Integer) where {T<:Real}

Construct the `FourierDerivativeOperator` on a uniform grid between `xmin` and
`xmax` using `N` nodes and `N÷2+1` complex Fourier modes.
"""
function FourierDerivativeOperator(xmin::T, xmax::T, N::Integer) where {T<:Real}
    @argcheck N >= 1

    jac = 2*T(π) / (xmax - xmin) / N # / N because of brfft instead of BRFFT
    Δx = (xmax - xmin) / N
    grid_evaluate = range(xmin, stop=xmax, length=N+1) # two boundary nodes
    grid_compute = range(xmin, stop=grid_evaluate[end-1], length=N)
    u = zero.(grid_compute)
    rfft_plan = plan_rfft(u)
    uhat = rfft_plan*u
    brfft_plan = plan_brfft(uhat, N)

    FourierDerivativeOperator(jac, Δx, grid_compute, grid_evaluate, uhat, rfft_plan, brfft_plan)
end

"""
    fourier_derivative_operator(xmin::Real, xmax::Real, N::Integer)
    fourier_derivative_operator(; xmin::Real, xmax::Real, N::Integer)

Construct the `FourierDerivativeOperator` on a uniform grid between `xmin` and
`xmax` using `N` nodes and `N÷2+1` complex Fourier modes.
"""
function fourier_derivative_operator(xmin::Real, xmax::Real, N::Integer)
    FourierDerivativeOperator(promote(xmin, xmax)..., N)
end

function fourier_derivative_operator(; xmin::Real, xmax::Real, N::Integer)
    FourierDerivativeOperator(promote(xmin, xmax)..., N)
end

derivative_order(D::FourierDerivativeOperator) = 1
LinearAlgebra.issymmetric(D::FourierDerivativeOperator) = false

function Base.show(io::IO, D::FourierDerivativeOperator)
    if get(io, :compact, false)
        print(io, "Periodic 1st derivative Fourier operator")
    else
        grid = D.grid_evaluate
        print(io, "Periodic 1st derivative Fourier operator {T=", eltype(D), "} \n")
        print(io, "on a grid in [", first(grid), ", ", last(grid),
                    "] using ", length(D.rfft_plan), " nodes and ",
                    length(D.brfft_plan), " modes")
    end
end

function mul!(dest::AbstractVector{T}, D::FourierDerivativeOperator, u::AbstractVector{T}) where {T}
    @unpack jac, tmp, rfft_plan, brfft_plan = D
    N, _ = size(D)
    @boundscheck begin
        @argcheck N == length(u)
        @argcheck N == length(dest)
    end

    mul!(tmp, rfft_plan, u)
    @inbounds @simd for j in Base.OneTo(length(tmp)-1)
        tmp[j] *= (j-1)*im * jac
    end
    # see e.g. Steven G. Johnson (2011) Notes on FFT based differentiation
    if iseven(N)
        @inbounds tmp[end] = zero(eltype(tmp))
    else
        @inbounds tmp[end] *= (length(tmp)-1)*im * jac
    end
    mul!(dest, brfft_plan, tmp)
end

# TODO there is no 5 argument mul! in FFTW.jl...
# function mul!(dest::AbstractVector{T}, D::FourierDerivativeOperator, u::AbstractVector{T}, α, β) where {T}
#     @unpack jac, tmp, rfft_plan, brfft_plan = D
#     N, _ = size(D)
#     @boundscheck begin
#         @argcheck N == length(u)
#         @argcheck N == length(dest)
#     end

#     mul!(tmp, rfft_plan, u)
#     @inbounds @simd for j in Base.OneTo(length(tmp)-1)
#         tmp[j] *= (j-1)*im * jac
#     end
#     if iseven(N)
#         @inbounds tmp[end] = zero(eltype(tmp))
#     else
#         @inbounds tmp[end] *= (length(tmp)-1)*im * jac
#     end
#     mul!(dest, brfft_plan, tmp, α, β)
# end

function integrate(func, u::AbstractVector, D::FourierDerivativeOperator)
    @boundscheck begin
        length(u) == length(grid(D))
    end
    @unpack Δx = D

    @inbounds res = sum(func, u)

    Δx * res
end

function mass_matrix(D::FourierDerivativeOperator)
    @unpack Δx = D

    Δx * I
end

"""
    fourier_derivative_matrix(N, xmin::Real=0.0, xmax::Real=2π)

Compute the Fourier derivative matrix with respect to the corresponding nodal
basis using `N` nodes, see
Kopriva (2009) Implementing Spectral Methods for PDEs, Algorithm 18.
"""
function fourier_derivative_matrix(N, xmin::Real=0.0, xmax::Real=2π)
    T = promote_type(typeof(xmin), typeof(xmax))
    jac_2 = T(π) / (xmax - xmin)
    D = Array{T}(undef, N, N)
    @inbounds for j in 1:N, i in 1:N
        j == i && continue
        D[i,j] = (-1)^(i+j) * cot((i-j)*T(π)/N) * jac_2
        D[i,i] -= D[i,j]
    end
    D
end

function _coef_end(coef, D1::FourierDerivativeOperator{T}) where {T}
    if isodd(size(D1, 1))
        return coef
    else
        return ntuple(i -> isodd(i) ? coef[i] : zero(T), length(coef))
    end
end
# fallback for periodic FD operators etc.
_coef_end(coef, D1) = coef

struct FourierPolynomialDerivativeOperator{T<:Real, Grid, RFFT, BRFFT, N} <: AbstractPeriodicDerivativeOperator{T}
    D1::FourierDerivativeOperator{T,Grid,RFFT,BRFFT}
    coef::NTuple{N,T}
    coef_end::NTuple{N,T}

    function FourierPolynomialDerivativeOperator(D1::FourierDerivativeOperator{T,Grid,RFFT,BRFFT}, coef::NTuple{N,T}) where {T<:Real, Grid, RFFT, BRFFT, N}
        coef_end = _coef_end(coef, D1)
        new{T,Grid,RFFT,BRFFT,N}(D1, coef, coef_end)
    end
end

function FourierPolynomialDerivativeOperator(D1::FourierDerivativeOperator)
    T = eltype(D1)
    FourierPolynomialDerivativeOperator(D1, (zero(T), one(T)))
end

Base.size(poly::FourierPolynomialDerivativeOperator) = size(poly.D1)
function LinearAlgebra.issymmetric(poly::FourierPolynomialDerivativeOperator)
    @unpack coef = poly
    all(iszero, coef[idx] for idx in eachindex(coef) if iseven(idx))
end
grid(poly::FourierPolynomialDerivativeOperator) = grid(poly.D1)

function Base.show(io::IO, poly::FourierPolynomialDerivativeOperator)
    if get(io, :compact, false)
        print(io, "Fourier polynomial")
    else
        print(io, "Fourier polynomial with coefficients\n")
        print(io, poly.coef)
        print(io, "\nof the operator:\n")
        print(io, poly.D1)
    end
end

function Base.:*(D1::FourierDerivativeOperator, D2::FourierDerivativeOperator)
    T = eltype(D1)
    @argcheck T == eltype(D2) ArgumentError
    @argcheck D1.jac == D2.jac ArgumentError
    @argcheck D1.Δx == D2.Δx ArgumentError
    @argcheck D1.grid_compute == D2.grid_compute DimensionMismatch
    @argcheck D1.grid_evaluate == D2.grid_evaluate DimensionMismatch

    FourierPolynomialDerivativeOperator(D1, (zero(T), zero(T), one(T)))
end

function Base.literal_pow(::typeof(^), D1::FourierDerivativeOperator, ::Val{P}) where {P}
    T = eltype(D1)
    coef = Base.setindex( ntuple(_->zero(T), Val{P+1}()), one(T), P+1)
    FourierPolynomialDerivativeOperator(D1, coef)
end

function Base.:*(factor::Union{Real,Integer}, poly::FourierPolynomialDerivativeOperator)
    @unpack coef = poly
    for idx in 1:length(coef)
        coef = Base.setindex(coef, factor*coef[idx], idx)
    end

    FourierPolynomialDerivativeOperator(poly.D1, coef)
end

function Base.:*(poly::FourierPolynomialDerivativeOperator, factor::Union{Real,Integer})
    factor * poly
end

function Base.:*(D::FourierDerivativeOperator, factor::Union{Real,Integer})
    FourierPolynomialDerivativeOperator(D) * factor
end

function Base.:*(factor::Union{Real,Integer}, D::FourierDerivativeOperator)
    D * factor
end

function Base.:*(poly::FourierPolynomialDerivativeOperator, scaling::UniformScaling)
    scaling.λ * poly
end

function Base.:*(scaling::UniformScaling, poly::FourierPolynomialDerivativeOperator)
    poly * scaling
end

function Base.:*(D::FourierDerivativeOperator, scaling::UniformScaling)
    scaling * FourierPolynomialDerivativeOperator(D)
end

function Base.:*(scaling::UniformScaling, D::FourierDerivativeOperator)
    FourierPolynomialDerivativeOperator(D) * scaling
end

function Base.:*(poly1::FourierPolynomialDerivativeOperator, poly2::FourierPolynomialDerivativeOperator)
    T = eltype(poly1.D1)
    @argcheck T == eltype(poly2.D1) ArgumentError
    @argcheck poly1.D1.jac == poly2.D1.jac ArgumentError
    @argcheck poly1.D1.Δx == poly2.D1.Δx ArgumentError
    @argcheck poly1.D1.grid_compute == poly2.D1.grid_compute DimensionMismatch
    @argcheck poly1.D1.grid_evaluate == poly2.D1.grid_evaluate DimensionMismatch

    coef = mul_poly(poly1.coef, poly2.coef)
    FourierPolynomialDerivativeOperator(poly1.D1, coef)
end

function Base.:*(poly1::FourierPolynomialDerivativeOperator, D2::FourierDerivativeOperator)
    poly1 * FourierPolynomialDerivativeOperator(D2)
end

function Base.:*(D1::FourierDerivativeOperator, poly2::FourierPolynomialDerivativeOperator)
    FourierPolynomialDerivativeOperator(D1) * poly2
end

function Base.:+(poly1::FourierPolynomialDerivativeOperator, poly2::FourierPolynomialDerivativeOperator)
    T = eltype(poly1.D1)
    @argcheck T == eltype(poly2.D1) ArgumentError
    @argcheck poly1.D1.jac == poly2.D1.jac ArgumentError
    @argcheck poly1.D1.Δx == poly2.D1.Δx ArgumentError
    @argcheck poly1.D1.grid_compute == poly2.D1.grid_compute DimensionMismatch
    @argcheck poly1.D1.grid_evaluate == poly2.D1.grid_evaluate DimensionMismatch

    coef = add_poly(poly1.coef, poly2.coef)
    FourierPolynomialDerivativeOperator(poly1.D1, coef)
end

function Base.:+(D1::FourierDerivativeOperator, poly2::FourierPolynomialDerivativeOperator)
    FourierPolynomialDerivativeOperator(D1) + poly2
end

function Base.:+(poly1::FourierPolynomialDerivativeOperator, D2::FourierDerivativeOperator)
    poly1 + FourierPolynomialDerivativeOperator(D2)
end

function Base.:+(poly::FourierPolynomialDerivativeOperator, scaling::UniformScaling)
    @unpack coef = poly
    coef = Base.setindex(coef, coef[1] + scaling.λ, 1)

    FourierPolynomialDerivativeOperator(poly.D1, coef)
end

function Base.:+(scaling::UniformScaling, poly::FourierPolynomialDerivativeOperator)
    poly + scaling
end

function Base.:+(D::FourierDerivativeOperator, scaling::UniformScaling)
    FourierPolynomialDerivativeOperator(D) + scaling
end

function Base.:+(scaling::UniformScaling, D::FourierDerivativeOperator)
    D + scaling
end

function Base.:-(poly1::FourierPolynomialDerivativeOperator, poly2::FourierPolynomialDerivativeOperator)
    T = eltype(poly1.D1)
    @argcheck T == eltype(poly2.D1) ArgumentError
    @argcheck poly1.D1.jac == poly2.D1.jac ArgumentError
    @argcheck poly1.D1.Δx == poly2.D1.Δx ArgumentError
    @argcheck poly1.D1.grid_compute == poly2.D1.grid_compute DimensionMismatch
    @argcheck poly1.D1.grid_evaluate == poly2.D1.grid_evaluate DimensionMismatch

    coef = subtract_poly(poly1.coef, poly2.coef)
    FourierPolynomialDerivativeOperator(poly1.D1, coef)
end

function Base.:-(D1::FourierDerivativeOperator, poly2::FourierPolynomialDerivativeOperator)
    FourierPolynomialDerivativeOperator(D1) - poly2
end

function Base.:-(poly1::FourierPolynomialDerivativeOperator, D2::FourierDerivativeOperator)
    poly1 - FourierPolynomialDerivativeOperator(D2)
end

function Base.:-(poly::FourierPolynomialDerivativeOperator, scaling::UniformScaling)
    @unpack coef = poly
    coef = Base.setindex(coef, coef[1] - scaling.λ, 1)

    FourierPolynomialDerivativeOperator(poly.D1, coef)
end

function Base.:-(scaling::UniformScaling, poly::FourierPolynomialDerivativeOperator)
    @unpack coef = poly
    coef = Base.setindex(coef, scaling.λ - coef[1], 1)
    for idx in 2:length(coef)
        coef = Base.setindex(coef, -coef[idx], idx)
    end

    FourierPolynomialDerivativeOperator(poly.D1, coef)
end

function Base.:-(D::FourierDerivativeOperator, scaling::UniformScaling)
    FourierPolynomialDerivativeOperator(D) - scaling
end

function Base.:-(scaling::UniformScaling, D::FourierDerivativeOperator)
    scaling - D
end

function mul!(dest::AbstractVector{T}, poly::FourierPolynomialDerivativeOperator, u::AbstractVector{T}) where {T}
    @unpack D1, coef, coef_end = poly
    @unpack jac, tmp, rfft_plan, brfft_plan = D1
    N, _ = size(D1)
    @boundscheck begin
        @argcheck N == length(u)
        @argcheck N == length(dest)
    end

    mul!(tmp, rfft_plan, u)
    @inbounds @simd for j in Base.OneTo(length(tmp)-1)
        # *N ) / N: brfft instead of irfft
        tmp[j] *= evalpoly((j-1)*im * jac*N, coef) / N
    end
    # see e.g. Steven G. Johnson (2011) Notes on FFT based differentiation
    @inbounds tmp[end] *= evalpoly((length(tmp)-1)*im * jac*N, coef_end) / N
    mul!(dest, brfft_plan, tmp)
end

function LinearAlgebra.ldiv!(dest::AbstractVector{T}, rat::FourierPolynomialDerivativeOperator, u::AbstractVector{T}) where {T}
    @unpack D1, coef, coef_end = rat
    @unpack jac, tmp, rfft_plan, brfft_plan = D1
    N, _ = size(D1)
    @boundscheck begin
        @argcheck N == length(u)
        @argcheck N == length(dest)
    end

    mul!(tmp, rfft_plan, u)
    @inbounds @simd for j in Base.OneTo(length(tmp)-1)
        # *N ) / N: brfft instead of irfft
        tmp[j] /= (evalpoly((j-1)*im * jac*N, coef) * N)
    end
    # see e.g. Steven G. Johnson (2011) Notes on FFT based differentiation
    @inbounds tmp[end] /= (evalpoly((length(tmp)-1)*im * jac*N, coef_end) * N)
    mul!(dest, brfft_plan, tmp)
end

struct FourierRationalDerivativeOperator{T<:Real, Grid, RFFT, BRFFT, Nnum, Nden} <: AbstractPeriodicDerivativeOperator{T}
    D1::FourierDerivativeOperator{T,Grid,RFFT,BRFFT}
    num_coef::NTuple{Nnum,T}
    num_coef_end::NTuple{Nnum,T}
    den_coef::NTuple{Nden,T}
    den_coef_end::NTuple{Nden,T}

    function FourierRationalDerivativeOperator(D1::FourierDerivativeOperator{T,Grid,RFFT,BRFFT}, num_coef::NTuple{Nnum,T}, den_coef::NTuple{Nden,T}) where {T<:Real, Grid, RFFT, BRFFT, Nnum, Nden}
        num_coef_end = _coef_end(num_coef, D1)
        den_coef_end = _coef_end(den_coef, D1)
        new{T,Grid,RFFT,BRFFT,Nnum,Nden}(D1, num_coef, num_coef_end, den_coef, den_coef_end)
    end
end

Base.size(rat::FourierRationalDerivativeOperator) = size(rat.D1)
grid(rat::FourierRationalDerivativeOperator) = grid(rat.D1)

function Base.show(io::IO, rat::FourierRationalDerivativeOperator)
    if get(io, :compact, false)
        print(io, "Rational Fourier operator")
    else
        print(io, "Rational Fourier operator with coefficients\n")
        print(io, rat.num_coef)
        print(io, "\nand\n")
        print(io, rat.den_coef)
        print(io, "\nof the operator:\n")
        print(io, rat.D1)
    end
end

function LinearAlgebra.issymmetric(rat::FourierRationalDerivativeOperator)
    @unpack num_coef, den_coef = rat
    num_is_even = all(iszero, num_coef[idx] for idx in eachindex(num_coef) if iseven(idx))
    den_is_even = all(iszero, den_coef[idx] for idx in eachindex(den_coef) if iseven(idx))

    num_is_even == den_is_even
end

function FourierRationalDerivativeOperator(num::FourierPolynomialDerivativeOperator)
    T = eltype(num)
    FourierRationalDerivativeOperator(num.D1, num.coef, (one(T),))
end

function FourierRationalDerivativeOperator(D::FourierDerivativeOperator)
    FourierRationalDerivativeOperator(FourierPolynomialDerivativeOperator(D))
end

function Base.:/(num::FourierPolynomialDerivativeOperator, den::FourierPolynomialDerivativeOperator)
    @argcheck num.D1.jac == den.D1.jac ArgumentError
    @argcheck num.D1.Δx == den.D1.Δx ArgumentError
    @argcheck num.D1.grid_compute == den.D1.grid_compute DimensionMismatch
    @argcheck num.D1.grid_evaluate == den.D1.grid_evaluate DimensionMismatch

    FourierRationalDerivativeOperator(num.D1, num.coef, den.coef)
end

function Base.:/(num::FourierDerivativeOperator, den::FourierPolynomialDerivativeOperator)
    FourierPolynomialDerivativeOperator(num) / den
end

function Base.:/(num::FourierPolynomialDerivativeOperator, den::FourierDerivativeOperator)
    num / FourierPolynomialDerivativeOperator(den)
end

function Base.inv(rat::FourierRationalDerivativeOperator)
    FourierRationalDerivativeOperator(rat.D1, rat.den_coef, rat.num_coef)
end

function Base.inv(den::Union{FourierDerivativeOperator,FourierPolynomialDerivativeOperator})
    inv(FourierRationalDerivativeOperator(den))
end

function Base.:+(rat1::FourierRationalDerivativeOperator, rat2::FourierRationalDerivativeOperator)
    T = eltype(rat1)
    @argcheck T == eltype(rat2) ArgumentError
    @argcheck rat1.D1.jac == rat2.D1.jac ArgumentError
    @argcheck rat1.D1.Δx == rat2.D1.Δx ArgumentError
    @argcheck rat1.D1.grid_compute == rat2.D1.grid_compute DimensionMismatch
    @argcheck rat1.D1.grid_evaluate == rat2.D1.grid_evaluate DimensionMismatch

    num_coef = add_poly(mul_poly(rat1.num_coef, rat2.den_coef), mul_poly(rat1.den_coef, rat2.num_coef))
    den_coef = mul_poly(rat1.den_coef, rat2.den_coef)
    FourierRationalDerivativeOperator(rat1.D1, num_coef, den_coef)
end

function Base.:+(rat1::FourierRationalDerivativeOperator, rat2::Union{FourierDerivativeOperator,FourierPolynomialDerivativeOperator})
    rat1 + FourierRationalDerivativeOperator(rat2)
end

function Base.:+(rat1::Union{FourierDerivativeOperator,FourierPolynomialDerivativeOperator}, rat2::FourierRationalDerivativeOperator)
    FourierRationalDerivativeOperator(rat1) + rat2
end

function Base.:-(rat1::FourierRationalDerivativeOperator, rat2::FourierRationalDerivativeOperator)
    T = eltype(rat1)
    @argcheck T == eltype(rat2) ArgumentError
    @argcheck rat1.D1.jac == rat2.D1.jac ArgumentError
    @argcheck rat1.D1.Δx == rat2.D1.Δx ArgumentError
    @argcheck rat1.D1.grid_compute == rat2.D1.grid_compute DimensionMismatch
    @argcheck rat1.D1.grid_evaluate == rat2.D1.grid_evaluate DimensionMismatch

    num_coef = subtract_poly(mul_poly(rat1.num_coef, rat2.den_coef), mul_poly(rat1.den_coef, rat2.num_coef))
    den_coef = mul_poly(rat1.den_coef, rat2.den_coef)
    FourierRationalDerivativeOperator(rat1.D1, num_coef, den_coef)
end

function Base.:-(rat1::FourierRationalDerivativeOperator, rat2::Union{FourierDerivativeOperator,FourierPolynomialDerivativeOperator})
    rat1 - FourierRationalDerivativeOperator(rat2)
end

function Base.:-(rat1::Union{FourierDerivativeOperator,FourierPolynomialDerivativeOperator}, rat2::FourierRationalDerivativeOperator)
    FourierRationalDerivativeOperator(rat1) - rat2
end

function Base.:*(rat1::FourierRationalDerivativeOperator, rat2::FourierRationalDerivativeOperator)
    T = eltype(rat1)
    @argcheck T == eltype(rat2) ArgumentError
    @argcheck rat1.D1.jac == rat2.D1.jac ArgumentError
    @argcheck rat1.D1.Δx == rat2.D1.Δx ArgumentError
    @argcheck rat1.D1.grid_compute == rat2.D1.grid_compute DimensionMismatch
    @argcheck rat1.D1.grid_evaluate == rat2.D1.grid_evaluate DimensionMismatch

    num_coef = mul_poly(rat1.num_coef, rat2.num_coef)
    den_coef = mul_poly(rat1.den_coef, rat2.den_coef)
    FourierRationalDerivativeOperator(rat1.D1, num_coef, den_coef)
end

function Base.:*(rat1::FourierRationalDerivativeOperator, rat2::Union{FourierDerivativeOperator,FourierPolynomialDerivativeOperator})
    rat1 * FourierRationalDerivativeOperator(rat2)
end

function Base.:*(rat1::Union{FourierDerivativeOperator,FourierPolynomialDerivativeOperator}, rat2::FourierRationalDerivativeOperator)
    FourierRationalDerivativeOperator(rat1) * rat2
end

function Base.:/(rat1::FourierRationalDerivativeOperator, rat2::FourierRationalDerivativeOperator)
    T = eltype(rat1)
    @argcheck T == eltype(rat2) ArgumentError
    @argcheck rat1.D1.jac == rat2.D1.jac ArgumentError
    @argcheck rat1.D1.Δx == rat2.D1.Δx ArgumentError
    @argcheck rat1.D1.grid_compute == rat2.D1.grid_compute DimensionMismatch
    @argcheck rat1.D1.grid_evaluate == rat2.D1.grid_evaluate DimensionMismatch

    num_coef = mul_poly(rat1.num_coef, rat2.den_coef)
    den_coef = mul_poly(rat1.den_coef, rat2.num_coef)
    FourierRationalDerivativeOperator(rat1.D1, num_coef, den_coef)
end

function Base.:/(rat1::FourierRationalDerivativeOperator, rat2::Union{FourierDerivativeOperator,FourierPolynomialDerivativeOperator})
    rat1 / FourierRationalDerivativeOperator(rat2)
end

function Base.:/(rat1::Union{FourierDerivativeOperator,FourierPolynomialDerivativeOperator}, rat2::FourierRationalDerivativeOperator)
    FourierRationalDerivativeOperator(rat1) / rat2
end

function mul!(dest::AbstractVector{T}, rat::FourierRationalDerivativeOperator, u::AbstractVector{T}) where {T}
    @unpack D1, num_coef, num_coef_end, den_coef, den_coef_end = rat
    @unpack jac, tmp, rfft_plan, brfft_plan = D1
    N, _ = size(D1)
    @boundscheck begin
        @argcheck N == length(u)
        @argcheck N == length(dest)
    end

    mul!(tmp, rfft_plan, u)
    @inbounds @simd for j in Base.OneTo(length(tmp)-1)
        # *N ) / N: brfft instead of irfft
        tmp[j] *= evalpoly((j-1)*im * jac*N, num_coef) / (N * evalpoly((j-1)*im * jac*N, den_coef))
    end
    # see e.g. Steven G. Johnson (2011) Notes on FFT based differentiation
    @inbounds tmp[end] *= evalpoly((length(tmp)-1)*im * jac*N, num_coef_end) / (N * evalpoly((length(tmp)-1)*im * jac*N, den_coef_end))
    mul!(dest, brfft_plan, tmp)
end

function LinearAlgebra.ldiv!(dest::AbstractVector{T}, rat::FourierRationalDerivativeOperator, u::AbstractVector{T}) where {T}
    @unpack D1, num_coef, num_coef_end, den_coef, den_coef_end = rat
    @unpack jac, tmp, rfft_plan, brfft_plan = D1
    N, _ = size(D1)
    @boundscheck begin
        @argcheck N == length(u)
        @argcheck N == length(dest)
    end

    mul!(tmp, rfft_plan, u)
    @inbounds @simd for j in Base.OneTo(length(tmp)-1)
        # *N ) / N: brfft instead of irfft
        tmp[j] *= evalpoly((j-1)*im * jac*N, den_coef) / (N * evalpoly((j-1)*im * jac*N, num_coef))
    end
    # see e.g. Steven G. Johnson (2011) Notes on FFT based differentiation
    @inbounds tmp[end] *= evalpoly((N-1)*im * jac*N, den_coef_end) / (N * evalpoly((N-1)*im * jac*N, num_coef_end))
    mul!(dest, brfft_plan, tmp)
end

function Base.:\(rat::Union{FourierRationalDerivativeOperator,FourierPolynomialDerivativeOperator}, u::AbstractVector{T}) where {T}
    dest = similar(u)
    ldiv!(dest, rat, u)
end

struct PeriodicDerivativeOperatorQuotient{T<:Real, numDtype<:Union{PeriodicDerivativeOperator{T},FourierDerivativeOperator{T}}, denDtype<:Union{PeriodicDerivativeOperator{T},FourierDerivativeOperator{T}}, Nnum, Nden, RFFT, IRFFT} <: AbstractPeriodicDerivativeOperator{T}
    num_D::numDtype
    den_D::denDtype
    num_coef::NTuple{Nnum,T}
    den_coef::NTuple{Nden,T}
    num_coef_end::NTuple{Nnum,T}
    den_coef_end::NTuple{Nden,T}
    tmp::Vector{Complex{T}}
    num_eigval::Vector{Complex{T}}
    den_eigval::Vector{Complex{T}}
    rfft_plan::RFFT
    irfft_plan::IRFFT

    function PeriodicDerivativeOperatorQuotient(num_D::numDtype, den_D::denDtype, num_coef::NTuple{Nnum,T}, den_coef::NTuple{Nden,T}, tmp::Vector{Complex{T}}, num_eigval::Vector{Complex{T}}, den_eigval::Vector{Complex{T}}, rfft_plan::RFFT, irfft_plan::IRFFT) where {T<:Real, numDtype<:Union{PeriodicDerivativeOperator{T},FourierDerivativeOperator{T}}, denDtype<:Union{PeriodicDerivativeOperator{T},FourierDerivativeOperator{T}}, Nnum, Nden, RFFT, IRFFT}
        @argcheck length(irfft_plan) == length(tmp) DimensionMismatch
        @argcheck length(irfft_plan) == length(num_eigval) DimensionMismatch
        @argcheck length(irfft_plan) == length(den_eigval) DimensionMismatch
        @argcheck length(irfft_plan) == (length(rfft_plan)÷2)+1 DimensionMismatch
        @argcheck length(grid(num_D)) == length(rfft_plan) DimensionMismatch
        @argcheck length(grid(den_D)) == length(rfft_plan) DimensionMismatch

        num_coef_end = _coef_end(num_coef, num_D)
        den_coef_end = _coef_end(den_coef, den_D)

        new{T,numDtype,denDtype,Nnum,Nden,RFFT,IRFFT}(num_D, den_D, num_coef, den_coef, num_coef_end, den_coef_end, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan)
    end
end

Base.size(quot::PeriodicDerivativeOperatorQuotient) = size(quot.num_D)
grid(quot::PeriodicDerivativeOperatorQuotient) = grid(quot.num_D)

function Base.show(io::IO, quot::PeriodicDerivativeOperatorQuotient)
    if get(io, :compact, false)
        print(io, "Rational operator (", quot.num_D, ", ", quot.den_D, ")")
    else
        print(io, "Quotient of the polynomial with coefficients\n")
        print(io, quot.num_coef)
        print(io, "\nof the operator:\n")
        print(io, quot.num_D)
        print(io, "and the polynomial with coefficients\n")
        print(io, quot.den_coef)
        print(io, "\nof the operator:\n")
        print(io, quot.den_D)
    end
end

function _eigvals!(eigval::Vector{Complex{T}}, D::FourierDerivativeOperator{T}) where {T<:Real}
    N = length(grid(D))
    jac = D.jac * N
    @inbounds for idx in 1:(length(eigval))
        eigval[idx] = (idx-1)*im * jac
    end

    eigval
end

function Base.://(num::FourierDerivativeOperator, den::PeriodicRationalDerivativeOperator)
    T = eltype(num)
    @argcheck T == eltype(den) ArgumentError
    @argcheck grid(num) == grid(den) ArgumentError
    @argcheck den.den_coef == (one(T),) ArgumentError

    @unpack tmp, rfft_plan = num
    irfft_plan = plan_irfft(tmp, length(grid(num)))

    num_eigval = similar(tmp)
    _eigvals!(num_eigval, num)

    den_eigval = similar(tmp)
    _eigvals!(den_eigval, den.D)

    PeriodicDerivativeOperatorQuotient(num, den.D, (zero(T), one(T)), den.num_coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan)
end

function Base.://(num::FourierPolynomialDerivativeOperator, den::PeriodicRationalDerivativeOperator)
    T = eltype(num)
    @argcheck T == eltype(den) ArgumentError
    @argcheck grid(num) == grid(den) ArgumentError
    @argcheck den.den_coef == (one(T),) ArgumentError

    @unpack tmp, rfft_plan = num.D1
    irfft_plan = plan_irfft(tmp, length(grid(num)))

    num_eigval = similar(tmp)
    _eigvals!(num_eigval, num.D1)

    den_eigval = similar(tmp)
    _eigvals!(den_eigval, den.D)

    PeriodicDerivativeOperatorQuotient(num.D1, den.D, num.coef, den.num_coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan)
end

function Base.://(num::PeriodicDerivativeOperator, den::PeriodicRationalDerivativeOperator)
    T = eltype(num)
    @argcheck T == eltype(den) ArgumentError
    @argcheck grid(num) == grid(den) ArgumentError
    @argcheck den.den_coef == (one(T),) ArgumentError

    x = grid(num)
    u = zero.(x)
    rfft_plan = plan_rfft(u)
    tmp = rfft_plan * u
    irfft_plan = plan_irfft(tmp, length(grid(num)))

    num_eigval = similar(tmp)
    _eigvals!(num_eigval, num)

    den_eigval = similar(tmp)
    _eigvals!(den_eigval, den.D)

    PeriodicDerivativeOperatorQuotient(num, den.D, (zero(T), one(T)), den.num_coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan)
end

function Base.://(num::PeriodicRationalDerivativeOperator, den::PeriodicRationalDerivativeOperator)
    T = eltype(num)
    @argcheck T == eltype(den) ArgumentError
    @argcheck grid(num) == grid(den) ArgumentError
    @argcheck num.den_coef == (one(T),) ArgumentError
    @argcheck den.den_coef == (one(T),) ArgumentError

    x = grid(num)
    u = zero.(x)
    rfft_plan = plan_rfft(u)
    tmp = rfft_plan * u
    irfft_plan = plan_irfft(tmp, length(grid(num)))

    num_eigval = similar(tmp)
    _eigvals!(num_eigval, num.D)

    den_eigval = similar(tmp)
    _eigvals!(den_eigval, den.D)

    PeriodicDerivativeOperatorQuotient(num.D, den.D, num.num_coef, den.num_coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan)
end

function Base.://(num::FourierDerivativeOperator, den::FourierPolynomialDerivativeOperator)
    T = eltype(num)
    @argcheck T == eltype(den) ArgumentError
    @argcheck grid(num) == grid(den) ArgumentError

    @unpack tmp, rfft_plan = den.D1
    irfft_plan = plan_irfft(tmp, length(grid(num)))

    num_eigval = similar(tmp)
    _eigvals!(num_eigval, num)

    den_eigval = similar(tmp)
    _eigvals!(den_eigval, den.D1)

    PeriodicDerivativeOperatorQuotient(num, den.D1, (zero(T), one(T)), den.coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan)
end

function Base.://(num::FourierPolynomialDerivativeOperator, den::FourierPolynomialDerivativeOperator)
    T = eltype(num)
    @argcheck T == eltype(den) ArgumentError
    @argcheck grid(num) == grid(den) ArgumentError

    @unpack tmp, rfft_plan = den.D1
    irfft_plan = plan_irfft(tmp, length(grid(num)))

    num_eigval = similar(tmp)
    _eigvals!(num_eigval, num.D1)

    den_eigval = similar(tmp)
    _eigvals!(den_eigval, den.D1)

    PeriodicDerivativeOperatorQuotient(num.D1, den.D1, num.coef, den.coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan)
end

function Base.://(num::PeriodicDerivativeOperator, den::FourierPolynomialDerivativeOperator)
    T = eltype(num)
    @argcheck T == eltype(den) ArgumentError
    @argcheck grid(num) == grid(den) ArgumentError

    @unpack tmp, rfft_plan = den.D1
    irfft_plan = plan_irfft(tmp, length(grid(num)))

    num_eigval = similar(tmp)
    _eigvals!(num_eigval, num)

    den_eigval = similar(tmp)
    _eigvals!(den_eigval, den.D1)

    PeriodicDerivativeOperatorQuotient(num, den.D1, (zero(T), one(T)), den.coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan)
end

function Base.://(num::PeriodicRationalDerivativeOperator, den::FourierPolynomialDerivativeOperator)
    T = eltype(num)
    @argcheck T == eltype(den) ArgumentError
    @argcheck grid(num) == grid(den) ArgumentError
    @argcheck num.den_coef == (one(T),) ArgumentError

    @unpack tmp, rfft_plan = den.D1
    irfft_plan = plan_irfft(tmp, length(grid(num)))

    num_eigval = similar(tmp)
    _eigvals!(num_eigval, num.D)

    den_eigval = similar(tmp)
    _eigvals!(den_eigval, den.D1)

    PeriodicDerivativeOperatorQuotient(num.D, den.D1, num.num_coef, den.coef, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan)
end

function mul!(dest::AbstractVector{T}, quot::PeriodicDerivativeOperatorQuotient, u::AbstractVector{T}) where {T}
    @unpack num_D, num_coef, den_coef, num_coef_end, den_coef_end, tmp, num_eigval, den_eigval, rfft_plan, irfft_plan = quot
    N, _ = size(num_D)
    @boundscheck begin
        @argcheck N == length(u)
        @argcheck N == length(dest)
    end

    mul!(tmp, rfft_plan, u)
    @inbounds @simd for j in Base.OneTo(length(tmp)-1)
        tmp[j] *= evalpoly(num_eigval[j], num_coef) / evalpoly(den_eigval[j], den_coef)
    end
    tmp[end] *= evalpoly(num_eigval[end], num_coef_end) / evalpoly(den_eigval[end], den_coef_end)
    mul!(dest, irfft_plan, tmp)
end

"""
    ConstantFilter(D::FourierDerivativeOperator, filter)

Create a modal filter with constant parameters adapted to the Fourier
derivative operator `D` with parameters given by the filter function `filter`.
"""
function ConstantFilter(D::FourierDerivativeOperator{T}, filter) where {T}
    Np1 = length(D.brfft_plan)
    coefficients = Array{T}(undef, Np1)
    set_filter_coefficients!(coefficients, filter)
    tmp = copy(D.tmp)
    modal2nodal = plan_irfft(D.tmp, length(D.rfft_plan))
    nodal2modal = D.rfft_plan
    ConstantFilter(coefficients, nodal2modal, modal2nodal, tmp, filter)
end

abstract type AbstractFourierViscosity{T} <: AbstractPeriodicDerivativeOperator{T} end

@inline source_of_coefficients(Di::AbstractFourierViscosity) = (Di.source_of_coefficients)

LinearAlgebra.issymmetric(Di::AbstractFourierViscosity) = true
grid(Di::AbstractFourierViscosity) = grid(Di.D)

function mul!(dest::AbstractVector{T}, Di::AbstractFourierViscosity{T}, u::AbstractVector{T}) where {T}
    @unpack coefficients, D = Di
    @unpack tmp, rfft_plan, brfft_plan = D
    N = size(D, 1)
    @boundscheck begin
        @argcheck N == length(u)
        @argcheck N == length(dest)
        @argcheck length(tmp) == length(coefficients)
    end

    mul!(tmp, rfft_plan, u)
    @inbounds @simd for j in Base.OneTo(length(tmp))
        tmp[j] *= coefficients[j]
    end
    mul!(dest, brfft_plan, tmp)

    nothing
end

"""
    FourierConstantViscosity

Fourier viscosity operator with constant coefficients for the periodic 1st
derivative Fourier operator.
"""
struct FourierConstantViscosity{T<:Real, Grid, RFFT, BRFFT} <: AbstractFourierViscosity{T}
    coefficients::Vector{T}
    D::FourierDerivativeOperator{T,Grid,RFFT,BRFFT}
    parameters
    source_of_coefficients

    function FourierConstantViscosity(D::FourierDerivativeOperator{T,Grid,RFFT,BRFFT}, parameters, source_of_coefficients) where {T<:Real, Grid, RFFT, BRFFT}
        # precompute coefficients
        N = size(D,1)
        jac = N * D.jac # *N: brfft instead of irfft
        coefficients = Vector{T}(undef, length(D.brfft_plan))
        set_filter_coefficients!(coefficients, jac, N, parameters, source_of_coefficients)
        new{T,Grid,RFFT,BRFFT}(coefficients, D, parameters, source_of_coefficients)
    end
end

function Base.show(io::IO, Di::FourierConstantViscosity)
    if get(io, :compact, false)
        print(io, "Fourier viscosity operator with constant coefficients of ")
    else
        grid = Di.D.grid_evaluate
        print(io, "Fourier viscosity operator with constant coefficients for the periodic 1st\n")
        print(io, "derivative Fourier operator {T=", eltype(Di), "} on a grid in [",
                first(grid), ", ", last(grid), "]\n")
        print(io, "using ", length(Di.D.rfft_plan), " nodes and ",
                length(Di.D.brfft_plan), " modes with coefficients of ")
    end
    print(io, Di.source_of_coefficients)
end

function dissipation_operator(source_of_coefficients, D::FourierDerivativeOperator; kwargs...)
    parameters = get_parameters(source_of_coefficients, D; kwargs...)
    FourierConstantViscosity(D, parameters, source_of_coefficients)
end

"""
    Tadmor1989()

Coefficients of the Fourier spectral viscosity given in
- Tadmor (1989)
  Convergence of Spectral Methods for Nonlinear Conservation Laws.
  SIAM Journal on Numerical Analysis 26.1, pp. 30-44.
"""
struct Tadmor1989 <: SourceOfCoefficients end

function Base.show(io::IO, source::Tadmor1989)
    if get(io, :compact, false)
        summary(io, source)
    else
        print(io,
            "Tadmor (1989) \n",
            "  Convergence of Spectral Methods for Nonlinear Conservation Laws. \n",
            "  SIAM Journal on Numerical Analysis 26.1, pp. 30-44. \n")
    end
end

function set_filter_coefficients!(coefficients::AbstractVector{T},
                                  jac::T, N::Int,
                                  parameters, source::Tadmor1989) where {T<:Real}
    @unpack strength, cutoff = parameters
    @argcheck cutoff >= 1
    fill!(coefficients, zero(T))
    jac = jac^2 / N # ^2: 2nd derivative; /N: brfft instead of irfft
    @simd for j in cutoff:length(coefficients)
        @inbounds coefficients[j] = -strength * (j-1)^2 * jac
    end
end

"""
    MadayTadmor1989()

Coefficients of the Fourier spectral viscosity given in
- Maday, Tadmor (1989)
  Analysis of the Spectral Vanishing Viscosity Method for Periodic Conservation
    Laws.
  SIAM Journal on Numerical Analysis 26.4, pp. 854-870.
"""
struct MadayTadmor1989 <: SourceOfCoefficients end

function Base.show(io::IO, source::MadayTadmor1989)
    if get(io, :compact, false)
        summary(io, source)
    else
        print(io,
            "Maday, Tadmor (1989) \n",
            "  Analysis of the Spectral Vanishing Viscosity Method for Periodic Conservation\n",
            "    Laws. \n",
            "  SIAM Journal on Numerical Analysis 26.4, pp. 854-870.")
    end
end

function set_filter_coefficients!(coefficients::AbstractVector{T},
                                  jac::T, N::Int,
                                  parameters, source::MadayTadmor1989) where {T<:Real}
    @unpack strength, cutoff = parameters
    @argcheck cutoff >= 1
    fill!(coefficients, zero(T))
    jac = jac^2 / N # ^2: 2nd derivative; /N: brfft instead of irfft
    @inbounds @simd for j in cutoff:min(2cutoff,length(coefficients))
        coefficients[j] = -strength * (j-1)^2 * jac * (j-cutoff)/cutoff
    end
    @inbounds @simd for j in 2cutoff:length(coefficients)
        coefficients[j] = -strength * (j-1)^2 * jac
    end
end

"""
    TadmorWaagan2012Standard()

Coefficients of the Fourier spectral viscosity given in
- Tadmor, Waagan (2012)
  Adaptive Spectral Viscosity for Hyperbolic Conservation Laws.
  SIAM Journal on Scientific Computing 34.2, pp. A993-A1009.
"""
struct TadmorWaagan2012Standard <: SourceOfCoefficients end

function Base.show(io::IO, source::TadmorWaagan2012Standard)
    if get(io, :compact, false)
        summary(io, source)
    else
        print(io,
            "Tadmor, Waagan (2012) \n",
            "  Adaptive Spectral Viscosity for Hyperbolic Conservation Laws. \n",
            "  SIAM Journal on Scientific Computing 34.2, pp. A993-A1009.")
    end
end

function set_filter_coefficients!(coefficients::AbstractVector{T},
                                  jac::T, N::Int,
                                  parameters, source::TadmorWaagan2012Standard) where {T<:Real}
    @unpack strength, cutoff = parameters
    @argcheck cutoff >= 1
    fill!(coefficients, zero(T))
    jac = jac^2 / N # ^2: 2nd derivative; /N: brfft instead of irfft
    @inbounds @simd for j in cutoff+1:length(coefficients)
        coefficients[j] = -strength * (j-1)^2 * jac * exp(-((length(coefficients)-j)/(j-cutoff))^2)
    end
end

"""
    TadmorWaagan2012Convergent()

Coefficients of the Fourier spectral viscosity given in
- Tadmor, Waagan (2012)
  Adaptive Spectral Viscosity for Hyperbolic Conservation Laws.
  SIAM Journal on Scientific Computing 34.2, pp. A993-A1009.

See also
- Schochet (1990)
  The Rate of Convergence of Spectral-Viscosity Methods for Periodic Scalar
    Conservation Laws.
  SIAM Journal on Numerical Analysis 27.5, pp. 1142-1159.
"""
struct TadmorWaagan2012Convergent <: SourceOfCoefficients end

function Base.show(io::IO, source::TadmorWaagan2012Convergent)
    if get(io, :compact, false)
        summary(io, source)
    else
        print(io,
            "Tadmor, Waagan (2012) \n",
            "  Adaptive Spectral Viscosity for Hyperbolic Conservation Laws. \n",
            "  SIAM Journal on Scientific Computing 34.2, pp. A993-A1009.")
    end
end

function set_filter_coefficients!(coefficients::AbstractVector{T},
                                  jac::T, N::Int,
                                  parameters, source::TadmorWaagan2012Convergent) where {T<:Real}
    @unpack strength, cutoff = parameters
    @argcheck cutoff >= 1
    fill!(coefficients, zero(T))
    jac = jac^2 / N # ^2: 2nd derivative; /N: brfft instead of irfft
    @inbounds @simd for j in cutoff:min(2cutoff,length(coefficients))
        coefficients[j] = -strength * (j-1)^2 * jac * exp(-((2cutoff-j)/(j-cutoff))^2)
    end
    @inbounds @simd for j in 2cutoff:length(coefficients)
        coefficients[j] = -strength * (j-1)^2 * jac
    end
end

const FourierSpectralViscosityCoefficients = Union{Tadmor1989,MadayTadmor1989,TadmorWaagan2012Standard,TadmorWaagan2012Convergent}

function get_parameters(source_of_coefficients::FourierSpectralViscosityCoefficients,
                        D::FourierDerivativeOperator;
                        strength::Real=eltype(D)(1)/size(D,2),
                        cutoff::Int=1+round(Int, sqrt(size(D,2))), #+1: 1 based indexing
                        kwargs...)
    @argcheck cutoff >= 1

    Dict(:strength=>strength, :cutoff=>cutoff)
end

"""
    Tadmor1993()

Coefficients of the Fourier super spectral viscosity given in
- Tadmor (1993)
  Super Viscosity and Spectral Approximations of Nonlinear Conservation Laws.
  Numerical Methods for Fluid Dynamics IV, pp. 69-82.
"""
struct Tadmor1993 <: SourceOfCoefficients end

function Base.show(io::IO, source::Tadmor1993)
    if get(io, :compact, false)
        summary(io, source)
    else
        print(io,
            "Tadmor (1993) \n",
            "  Super Viscosity and Spectral Approximations of Nonlinear Conservation Laws. \n",
            "  Numerical Methods for Fluid Dynamics IV, pp. 69-82.")
    end
end

function set_filter_coefficients!(coefficients::AbstractVector{T},
                                  jac::T, N::Int,
                                  parameters, source::Tadmor1993) where {T<:Real}
    @unpack order, strength, cutoff = parameters
    @argcheck order >= 1
    @argcheck cutoff >= 1
    fill!(coefficients, zero(T))
    jac = jac^2order / N # ^2order: order-th derivative; /N: brfft instead of irfft
    @inbounds @simd for j in cutoff:length(coefficients)
        coefficients[j] = -strength * (j-1)^2order * jac * (1 - (cutoff/j)^2order)
    end
end

const FourierSuperSpectralViscosityCoefficients = Tadmor1993

function get_parameters(source::FourierSuperSpectralViscosityCoefficients,
                        D::FourierDerivativeOperator;
                        order::Int=1,
                        strength=eltype(D)(1)/size(D,2)^(2order-1),
                        cutoff::Int=1+round(Int, size(D,2)^(1-1/2order)), #+1: 1 based indexing
                        kwargs...)
    @argcheck order >= 1
    @argcheck cutoff >= 1

    Dict(:order=>order, :strength=>strength, :cutoff=>cutoff)
end


text:
### A Pluto.jl notebook ###
# v0.18.1

using Markdown
using InteractiveUtils

# ╔═╡ 58c3936e-c10a-40bd-b8f9-a15804c27e59
using DataFrames

# ╔═╡ 636f1f75-2055-4cd3-95bc-3e4f1885f79c
html"""
<style>
	main {
		margin: 0 auto;
		max-width: 2000px;
    	padding-left: max(160px, 10%);
    	padding-right: max(160px, 10%);
	}
</style>
"""

# ╔═╡ 2a70a1c7-054f-43d1-84ac-0c1b737233cb
df = DataFrame(a=1:3, b=4:6)

# ╔═╡ 8993cf19-ae9b-4add-95bc-78207dd6ec9c
md" #### `getindex` (used in df on the right hand side of =)."

# ╔═╡ d2959067-58d9-487f-af6f-2ce2980d63ae
let
	copyOfa = df[:, :a]
	copyOfa .= 0
	copyOfa
end

# ╔═╡ 2434c37c-a7d4-4e1a-985d-7824885fbab4
df

# ╔═╡ f4ef304f-2ca5-4e95-86b0-d83865d74aeb
md" ##### Note that column `:a` has not changed."

# ╔═╡ fb181af4-b916-4aac-b997-1b89dbb99fb5
md" #### `setindex` (used in df on left hand side of =)."

# ╔═╡ 2b49dcd4-4054-4300-b75a-c109b40fc815
begin
	df[:, :a] .= 0
	df
end

# ╔═╡ 957207cf-47ce-4f2d-ba3b-53b609401c8d
md" ##### Now df has changed!"

# ╔═╡ a373977d-c27b-4478-8beb-58b15b7995a9
begin
	df .= 0
	df
end

# ╔═╡ 6111833d-6c40-4de1-9692-93a276a5730a
md" ##### In `getindex` context a user can ask for the following behaviors:

1. get a vector “as is”, without copying; this is achieved with df[!, :a] and df.a;
2. get a vector with copying; that is achieved with df[:, :a];
3. get a view of the vector; this is achieved with view(df, :, :a) or view(df, !, :a)
"

# ╔═╡ 331bde26-9892-4045-8355-85f666514487
md""" ##### In `setindex!` + broadcasting assignment context the user might ask for the following:

1. set an existing vector in-place, this is done by: df[:, :a] = vec and df[:, :a] .= something

2. replace an existing vector without copying: df[!, :a] = vec or df.a = vec

3. replace an existing vector with copying: df[!, :a] .= vec and df[!, :a] .= something

4. add a new vector without copying: df[!, :new] = vec or df.new = vec

5. add a new vector with copying: df[:, :new] = vec or df[:, :new] .= something

or

6. df[!, :new] .= something

(this last rule breaks things a bit as ! and : behave here the in same way, but it was added for user convenience)
"""

# ╔═╡ 9e510dbd-2e9c-40b3-a5ca-4ab4aa02d019
begin
	a = [7; 8; 9]
	df[:, :a] = a
	df
end

# ╔═╡ 2d15edc9-1f27-4d89-8824-7f6f27385b67
let
	b = [-1; -2; -3]
	df[!, :a] = b
	df
end

# ╔═╡ 5e631910-7a4f-4651-b605-c12966468757
begin
	b = 10:12
	df
end

# ╔═╡ 2b33122b-4e96-42b4-a2f3-83b6f231fd1c
df

# ╔═╡ 9d8a3a01-7c1c-4505-b7d1-42849d44ce3d

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"

[compat]
DataFrames = "~1.3.2"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.8.0-beta1"
manifest_format = "2.0"
project_hash = "c5a2b0816e3f46eedf015440dc2fa38edb4deba5"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
version = "1.1.1"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.42.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
version = "0.5.0+0"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "ae02104e835f219b8930c7664b8012c93475c340"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.3.2"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Downloads]]
deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
version = "1.6.0"

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
version = "0.6.3"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
version = "7.81.0+0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
version = "1.10.2+0"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
version = "2.28.0+0"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
version = "2022.2.1"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
version = "1.2.0"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
version = "0.3.17+2"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
version = "1.8.0"

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "db3a23166af8aebf4db5ef87ac5b00d36eb771e2"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.0"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "dfb54c4e414caa595a1f2ed759b160f5a3ddcba5"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.3.1"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
version = "0.7.0"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
version = "1.0.0"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
git-tree-sha1 = "5ce79ce186cc678bbb5c5681ca3379d1ddae11a1"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.7.0"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
version = "1.10.0"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
version = "1.2.12+1"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
version = "5.0.1+0"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
version = "1.41.0+1"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
version = "16.2.1+1"
"""

# ╔═╡ Cell order:
# ╠═636f1f75-2055-4cd3-95bc-3e4f1885f79c
# ╠═58c3936e-c10a-40bd-b8f9-a15804c27e59
# ╠═2a70a1c7-054f-43d1-84ac-0c1b737233cb
# ╟─8993cf19-ae9b-4add-95bc-78207dd6ec9c
# ╠═d2959067-58d9-487f-af6f-2ce2980d63ae
# ╠═2434c37c-a7d4-4e1a-985d-7824885fbab4
# ╟─f4ef304f-2ca5-4e95-86b0-d83865d74aeb
# ╟─fb181af4-b916-4aac-b997-1b89dbb99fb5
# ╠═2b49dcd4-4054-4300-b75a-c109b40fc815
# ╟─957207cf-47ce-4f2d-ba3b-53b609401c8d
# ╠═a373977d-c27b-4478-8beb-58b15b7995a9
# ╟─6111833d-6c40-4de1-9692-93a276a5730a
# ╟─331bde26-9892-4045-8355-85f666514487
# ╠═9e510dbd-2e9c-40b3-a5ca-4ab4aa02d019
# ╠═2d15edc9-1f27-4d89-8824-7f6f27385b67
# ╠═5e631910-7a4f-4651-b605-c12966468757
# ╠═2b33122b-4e96-42b4-a2f3-83b6f231fd1c
# ╠═9d8a3a01-7c1c-4505-b7d1-42849d44ce3d
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl


text:
module ErrorFree

import Base:fma

export eftRecip, eftSquare, eftCube,
       eftSqrt,  eftRecipSqrt, # see comment for eftSqrt
       eftSum2,  eftSum2inOrder, eftDiff2, eftDiff2inOrder,
       eftProd2, eftDiv2,
       eftSum3,  eftSum3as2, eftSum3inOrder, eftSum3inOrderAs2,
       eftProd3, eftProd3as2,
       eftFMA,   eftFMS

#>   one operand error free transformations

@inline function eftSquare(a::Float64)
    x = a * a
    y = fma(a, a, -x)
    x,y
end

function eftCube(a::Float64)
    p = a*a; e = fma(a, a, -p)
    x = p*a; p = fma(p, a, -x)
    y = e*a
    x,y
end

@inline function eftRecip(a::Float64)
     x = one(Float64)/a
     y = -(fma(x,a,-1.0)/a)
     x,y
end

#=
   While not strictly an error-free transformation it is quite reliable and recommended for use.
   Augmented precision square roots, 2-D norms and discussion on correctly reounding sqrt(x^2 + y^2)
   by Nicolas Brisebarre, Mioara Joldes, Erik Martin-Dorel, Hean-Michel Muller, Peter Kornerup
=#
@inline function eftSqrt(a::Float64)
     x = sqrt(a)
     t = fma(x,-x,a)
     y = t / (x*2.0)
     x,y
end

@inline function eftRecipSqrt(a::Float64)
     r = 1.0/a
     x = sqrt(r)
     t = fma(x,-x,r)
     y = t / (x*2.0)
     x,y
end

#>   two operand error free transformations

"""
 (a,b) ↦ (x,y)\\
x⊕y ≖ a⊕b and x⊕y≖x
"""
@inline function eftSum2{T<:Float64}(a::T, b::T)
  x = a + b
  t = x - a
  y = (a - (x - t)) + (b - t)
  x,y
end

"""
 (a,b) ↦ (x,y)\\
x⊕y ≖ a⊕b and x⊕y≖x\\
 *presumes* |a| ≥ |b|
"""
@inline function eftSum2inOrder{T<:Float64}(a::T, b::T)
  x = a + b
  y = b - (x - a)
  x,y
end

"""
 (a,b) ↦ (x,y)\\
x⊕y ≖ a⊝b and x⊕y≖x
"""
@inline function eftDiff2{T<:Float64}(a::T, b::T)
  x = a - b
  t = x - a
  y = (a - (x - t)) - (b + t)
  x,y
end

"""
 (a,b) ↦ (x,y)\\
x⊕y ≖ a⊝b and x⊕y≖x\\
*presumes* |a| ≥ |b|
"""
@inline function eftDiff2inOrder{T<:Float64}(a::T, b::T)
  x = a - b
  y = (a - x) - b
  x,y
end

@inline function eftProd2{T<:Float64}(a::T, b::T)
    x = a * b
    y = fma(a, b, -x)
    x,y
end

@inline function eftDiv2{T<:Float64}(a::T,b::T)
     x = a/b
     y = -(fma(x,b,-a)/b)
     x,y
end

#>   three operand error free transformations

function eftSum3{T<:Float64}(a::T,b::T,c::T)
    s,t = eftSum2(b, c)
    x,u = eftSum2(a, s)
    y,z = eftSum2(u, t)
    x,y = eftSum2inOrder(x, y)
    x,y,z
end

function eftSum3inOrder{T<:Float64}(a::T,b::T,c::T)
    s,t = eftSum2inOrder(b, c)
    x,u = eftSum2inOrder(a, s)
    y,z = eftSum2inOrder(u, t)
    x,y = eftSum2inOrder(x, y)
    x,y,z
end

@inline function eftSum3as2{T<:Float64}(a::T,b::T,c::T)
    s,t = eftSum2(b, c)
    x,u = eftSum2(a, s)
    y = u+t
    x,y = eftSum2inOrder(x, y)
    x,y
end

@inline function eftSum3inOrderAs2{T<:Float64}(a::T,b::T,c::T)
    s,t = eftSum2inOrder(b, c)
    x,u = eftSum2inOrder(a, s)
    y = u+t
    x,y = eftSum2inOrder(x, y)
    x,y
end

function eftProd3{T<:Float64}(a::T, b::T, c::T)
    p,e = eftProd2(a,b)
    x,p = eftProd2(p,c)
    y,z = eftProd2(e,c)
    x,y,z
end

@inline function eftProd3as2{T<:Float64}(a::T, b::T, c::T)
    p,e = eftProd2(a,b)
    x,p = eftProd2(p,c)
    y = e*c
    x,y
end

end # module ErrorFree


text:
"""
	coverage(ls)

Return the coverage of the given link stream.
By definition, the coverage of a link stream is 1.
"""
coverage(ls::Union{LinkStream,DirectedLinkStream})=1.0

"""
	coverage(s)

Return the coverage of the given stream graph:
```math
cov \\left( S \\right) = \\frac{\\left| W \\right|}{\\left| T \\times V \\right|}
```
Note: For stream graphs with no node and/or no duration, this
function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function coverage(s::Union{StreamGraph,DirectedStreamGraph})
	if ((length(s.V)!=0) && (duration(s)!=0))
		sum([duration(n) for (k,n) in s.W])/(duration(s)*length(s.V))
	else
		0.0
	end
end


text:
"""
Calculation of differential cross sections
from structure functions.

Based on the Rxsecnc.h C++ code provided
by R. Aggarwal.
"""

using QCDNUM

# These globals will later be set by user
ZMass = 91.1876
WMass = 80.398
AlphaEM = 7.297352570e-03
GFermi = 1.16637e-05
TopMass = 171.2
BottomMass = 4.20

Vub = 41.2e-3
Vcb = 3.93e-3

Sin2ThetaW = 0.23127
Sin2ThetaC = 0.05
vu = 0.19164
vd = -0.34582
ve = -0.03746

au = 0.5
ad = -0.5
ae = -0.5

# Variables for Rxsecnc_xq2
sqrt_s = 318.1 # Should be configurable for calculation of y
Lepcharge = 1 # Should be configurable for calculation of Y

export dd_xsecnc_xq2_i
export f2_lo
export xf3_lo
export fl_lo

"""
    f2_lo(x, q2)

Calculate the f2_lo structure function term.
To be run after the evolution of PDFs with QCDNUM.
"""
function f2_lo(x::Float64, q2::Float64)::Float64

    # For weights
    pz = q2 / ((ZMass^2 + q2) * (4*Sin2ThetaW * (1 - Sin2ThetaW)))

    Au = 4.0/9.0 -2*pz*(2.0/3.0)*vu*ve + pz^2*(ve^2 + ae^2)*(vu^2 + au^2)

    Ad = 1.0/9.0 -2*pz*(-1.0/3.0)*vd*ve + pz^2*(ve^2 + ae^2)*(vd^2 + ad^2)

    # As in HERAPDF (top set to 0)
    weights = [0., Ad, Au, Ad, Au, Ad, 0., Ad, Au, Ad, Au, Ad, 0.]

    # Structure function calculation
    output = QCDNUM.zmstfun(2, weights, [x], [q2], 1, 0)

    output[1]
end

"""
    xf3_lo(x, q2)

Calculate the xf3_lo structure function term.
To be run after the evolution of PDFs with QCDNUM.
"""
function xf3_lo(x::Float64, q2::Float64)::Float64

    # For weights
    pz = q2 / ((ZMass^2 + q2) * (4*Sin2ThetaW * (1 - Sin2ThetaW)))

    Bu = -2*(2.0/3.0)*au*ae*pz + 4*au*ae*vu*ve*pz^2

    Bd = -2*(-1.0/3.0)*ad*ae*pz + 4*ad*ae*vd*ve*pz^2

    # weights = [-Bu, -Bd, -Bu, -Bd, -Bu, -Bd, 0.,  Bd, Bu, Bd, Bu, Bd, Bu]

    # As in HERAPDF (top set to 0)
    weights = [0., -Bd, -Bu, -Bd, -Bu, -Bd, 0.,  Bd, Bu, Bd, Bu, Bd, 0.]

    # Structure function calculation
    output = QCDNUM.zmstfun(3, weights, [x], [q2], 1, 0)

    output[1]
end

"""
    fl_lo(x, q2)

Calculate the fl_lo structure function term.
To be run after the evolution of PDFs with QCDNUM.
"""
function fl_lo(x::Float64, q2::Float64)::Float64

    # For weights
    pz = q2 / ((ZMass^2 + q2) * (4*Sin2ThetaW * (1 - Sin2ThetaW)))

    Au = 4.0/9.0 -2*pz*(2.0/3.0)*vu*ve + pz^2*(ve^2 + ae^2)*(vu^2 + au^2)

    Ad = 1.0/9.0 -2*pz*(-1.0/3.0)*vd*ve + pz^2(ve^2 + ae^2)*(vd^2 + ad^2)

    # As in HERAPDF (top set to 0)
    weights = [0., Ad, Au, Ad, Au, Ad, 0.,  Ad, Au, Ad, Au, Ad, 0.]

    output = QCDNUM.zmstfun(1, weights, [x], [q2], 1, 0)

    output[1]
end

"""
    rxsecnc_xq2_i(x, q2)

Reduced cross section for single x, q2.
"""
function rxsecnc_xq2_i(x::Float64, q2::Float64)::Float64

    rxsec = -1.0 #?
    y = 0.04  #?
    y = q2 / sqrt_s / sqrt_s / x

    Y_plus = 1 + (1 - y)^2
    Y_minus = 1 - (1 - y)^2

    F2 = f2_lo(x, q2)
    xF3 = xf3_lo(x, q2)
    FL = fl_lo(x, q2)

    if (Lepcharge == 1)

        rxsec =  F2 - (Y_minus / Y_plus)*xF3 - (y^2 / Y_plus) * FL

    elseif (Lepcharge == -1)

        rxsec =  F2 + (Y_minus / Y_plus)*xF3 - (y^2 / Y_plus) * FL

    end

    rxsec
end

"""
    rxsecnc_xq2(x, q2)

Reduced cross section for all bins.
"""
function rxsecnc_xq2(x_bin_cen::Array{Float64}, q2_bin_cen::Array{Float64})::Array{Float64}

    n_bins = length(x_bin_cen)

    # Check same number of x and q2 vals
    if (n_bins != length(q2_bin_cen))

        throw(DimensionMismatch("x_bin_cen and q2_bin_cen must have same length."))

    end

    # Calculate for each bin
    rxsec = zeros(n_bins)

    for i = 1:n_bins

        rxsec[i] = rxsecnc_xq2_i(x_bin_cen[i], q2_bin_cen[i])

    end

    rxsec
end

"""
    nc_propagator(q2, x)
"""
function nc_propagator(q2::Float64, x::Float64)::Float64

    y = q2 / sqrt_s / sqrt_s / x
    Yplus = 1 + (1 - y)^2
    alpha = AlphaEM
    conversion_factor = 0.3894e9;  # GeV^2 -> pb^2

    conversion_factor*2*π*Yplus*alpha^2 / (x * q2^2)
end

"""
    dd_xsecnc_xq2_i(x, q2)

Double differential cross section for single
x and q2 values.
NB: modifications needed to include pol and order.
"""
function dd_xsecnc_xq2_i(x::Float64, q2::Float64)

    dd_xsec = -1.0

    dd_xsec = nc_propagator(q2, x) * rxsecnc_xq2_i(x, q2)

    dd_xsec
end

"""
    dd_xsecnc_xq2(x_bin_cen, q2_bin_cen)

Double differential cross section for all x and
q2 bins.
NB: modifications needed to include pol and order.
"""
function dd_xsecnc_xq2(x_bin_cen::Array{Float64},
                       q2_bin_cen::Array{Float64})::Array{Float64}

    n_bins = length(x_bin_cen)

    # Check same number of x and q2 vals
    if (n_bins != length(q2_bin_cen))

        throw(DimensionMismatch("x_bin_cen and q2_bin_cen must have same length."))

    end

    xsec = zeros(n_bins)

    for i = 1:n_bins

        xsec[i] = dd_xsecnc_xq2_i(x_bin_cen[i], q2_bin_cen[i])

    end

  xsec
end


text:
InputFolder = './Images/Zebrafish/';
OutputFolder = './Results/Images/Zebrafish/';

@iA = '*.tif';

@fxg_gDestripe [iA] > [D];
params.Angle = 0;
params.Scale = 0.25;
params.TopRad = 40;
params.OpenRad = 80;
params.SharpAmount = 0.25;
/endf

/show iA >;
/show D >;
/keep D > jp2;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
"""
computes the eigenvalues of the polynomial matrix P
"""
function eig(P::AbstractArray{Pol{T}, N}, tol=0) where {T<:Number, N}

    e, d, Ar, Br = colNullities(P, tol)
    infMult = isempty(d) ? 0 : sum(collect(1:length(d)).*d)

    e1, d1, Ar, Br = rowNullities(Ar, Br, tol)
    e2, _, _, _ = rowNullities(P, tol)
    (e1 == e2 && isempty(d1)) || return @error "could not compute the eigenvalues. Try using a higher tolerance"
    eigs = eigvals(Ar, Br)
    append!(eigs, [Inf for _ in 1:infMult])

    return eigs
end

"""
computes the column indices of the polynomial matrix P
"""
function colIndices(P, tol=0)
    e, _, _, _ = colNullities(P, tol)
    _, _, Ar, Br = rowNullities(P, tol)
    e1, _, _, _ = colNullities(Ar, Br, tol)
    e1 == e || return @error "could not compute the column indices. Try using a higher tolerance"
    ϵ = null2indices(e)
    k = max(deg(P), 1)
    return [i-k+1 for i in ϵ]
end

"""
computes the row indices of the polynomial matrix P
"""
function rowIndices(P, tol=0)
    e, _, _, _ = rowNullities(P, tol)
    _, _, Ar, Br = colNullities(P, tol)
    e1, _, _, _ = rowNullities(Ar, Br, tol)
    e == e1 || return @error "could not compute the row indices. Try using a higher tolerance"

    return null2indices(e)
end

"""
computes how many times ∞ is an eigenvalue of the polynomial matrix P. That is, how many
times 0 is an eigenvalue of the reversal of P
"""
function infMultiplicities(P, tol=0)
    _, d, _, _ = colNullities(P, tol)
    _, d1, _, _ = rowNullities(P, tol)
    d == d1 || return @error "could not compute the infinite eigenvalues. Try using a higher tolerance"
    return inf2indices(d)
end

"""
computes the multiplicity of the eigenvalue λ of P.
"""
function eigMultiplicities(λ::Number, P, tol=0)
    e = eigNullities(λ, P, tol)
    return inf2indices(e)
end

"""
computes the multiplicity of each eigenvalue in arr.
"""
function eigMultiplicities(arr, P, tol=0)
    d = minimum(abs.(diff(arr)))
    d < 1e-6*maximum(abs.(arr)) && @warn "Some eigenvalues are close to each others (d=$d). The multiplicities might be incorrect"
    nulls = Dict{T, Array{Int, 1}}()
    @inbounds for λ in arr
        e = eigMultiplicities(λ, P, tol)
        push!(nulls, λ => e)
    end
    return nulls
end


text:
RenderObject(data::Dict{Symbol}, program, bbs=Signal(AABB{Float32}(Vec3f0(0),Vec3f0(1))), main=nothing) = RenderObject(convert(Dict{Symbol,Any}, data), program, bbs, main)

function Base.show(io::IO, obj::RenderObject)
    println(io, "RenderObject with ID: ", obj.id)
end

Base.getindex(obj::RenderObject, symbol::Symbol)         = obj.uniforms[symbol]
Base.setindex!(obj::RenderObject, value, symbol::Symbol) = obj.uniforms[symbol] = value

Base.getindex(obj::RenderObject, symbol::Symbol, x::Function)     = getindex(obj, Val{symbol}(), x)
Base.getindex(obj::RenderObject, ::Val{:prerender}, x::Function)  = obj.prerenderfunctions[x]
Base.getindex(obj::RenderObject, ::Val{:postrender}, x::Function) = obj.postrenderfunctions[x]

Base.setindex!(obj::RenderObject, value, symbol::Symbol, x::Function)     = setindex!(obj, value, Val{symbol}(), x)
Base.setindex!(obj::RenderObject, value, ::Val{:prerender}, x::Function)  = obj.prerenderfunctions[x] = value
Base.setindex!(obj::RenderObject, value, ::Val{:postrender}, x::Function) = obj.postrenderfunctions[x] = value

"""
Function which sets an argument of a Context/RenderObject.
If multiple RenderObjects are supplied, it'll try to set the same argument in all
of them.
"""
function set_arg!(robj::RenderObject, sym, value)
    current_val = robj[sym]
    set_arg!(robj, sym, current_val, value)
    # GLVisualize relies on reactives event system no for rendering
    # so if a change should be visible there must be an event to indicate change
    Reactive.post_empty()
    nothing
end
function set_arg!(robj::Context, sym, value)
    set_arg!(robj.children, sym, value)
    nothing
end
function set_arg!(robj::Vector, sym, value)
    for elem in robj
        set_arg!(elem, sym, value)
    end
    nothing
end

function set_arg!(robj::RenderObject, sym, to_update::GPUArray, value)
    update!(to_update, value)
end
function set_arg!(robj::RenderObject, sym, to_update, value)
    robj[sym] = value
end
function set_arg!(robj::RenderObject, sym, to_update::Signal, value::Signal)
    robj[sym] = value
end
function set_arg!(robj::RenderObject, sym, to_update::Signal, value)
    push!(to_update, value)
end

"""
Represents standard sets of function applied before rendering
"""
immutable StandardPrerender
end

@compat function (::StandardPrerender)()
    glEnable(GL_DEPTH_TEST)
    glDepthMask(GL_TRUE)
    glDepthFunc(GL_LEQUAL)
    glDisable(GL_STENCIL_TEST)
    glStencilMask(0xff)
    glDisable(GL_CULL_FACE)
    enabletransparency()
end
immutable StandardPostrender
    vao::GLVertexArray
    primitive::GLenum
end
@compat function (sp::StandardPostrender)()
    render(sp.vao, sp.primitive)
end
immutable StandardPostrenderInstanced{T}
    main::T
    vao::GLVertexArray
    primitive::GLenum
end
@compat function (sp::StandardPostrenderInstanced)()
    renderinstanced(sp.vao, value(sp.main), sp.primitive)
end

immutable EmptyPrerender
end
@compat function (sp::EmptyPrerender)()
end
export EmptyPrerender
export prerendertype

function instanced_renderobject(data, program, bb=Signal(AABB(Vec3f0(0), Vec3f0(1))), primitive::GLenum=GL_TRIANGLES, main=nothing)
    pre = StandardPrerender()
    robj = RenderObject(data, program, pre, nothing, bb, main)
    robj.postrenderfunction = StandardPostrenderInstanced(main, robj.vertexarray, primitive)
    robj
end

function std_renderobject(data, program, bb=Signal(AABB(Vec3f0(0), Vec3f0(1))), primitive=GL_TRIANGLES, main=nothing)
    pre = StandardPrerender()
    robj = RenderObject(data, program, pre, nothing, bb, main)
    robj.postrenderfunction = StandardPostrender(robj.vertexarray, primitive)
    robj
end

prerendertype{Pre}(::Type{RenderObject{Pre}}) = Pre
prerendertype{Pre}(::RenderObject{Pre}) = Pre

extract_renderable(context::Vector{RenderObject}) = context
extract_renderable(context::RenderObject) = [context]
extract_renderable{T <: Composable}(context::Vector{T}) = map(extract_renderable, context)
function extract_renderable(context::Context)
    result = extract_renderable(context.children[1])
    for elem in context.children[2:end]
        push!(result, extract_renderable(elem)...)
    end
    result
end
transformation(c::RenderObject) = c[:model]
transformation(c::RenderObject, model) = (c[:model] = const_lift(*, model, c[:model]))
transform!(c::RenderObject, model) = (c[:model] = const_lift(*, model, c[:model]))

"""
Copy function for a context. We only need to copy the children
"""
function Base.copy{T}(c::GLAbstraction.Context{T})
    new_children = [copy(child) for child in c.children]
    Context{T}(new_children, c.boundingbox, c.transformation)
end

"""
Copy function for a RenderObject. We only copy the uniform dict
"""
function Base.copy{Pre}(robj::RenderObject{Pre})
    uniforms = Dict{Symbol, Any}([(k,v) for (k,v) in robj.uniforms])
    robj = RenderObject{Pre}(
        robj.main,
        uniforms,
        robj.vertexarray,
        robj.prerenderfunction,
        robj.postrenderfunction,
        robj.boundingbox,
    )
    Context(robj)
end

# """
# If you have an array of OptimizedPrograms, you only need to put PreRender in front.
# """
# type OptimizedProgram{PreRender}
#     program::GLProgram
#     uniforms::FixedDict
#     vertexarray::GLVertexArray
#     gl_parameters::PreRender
#     renderfunc::Callable
#     visible::Boolean
# end


text:
################################################################################
#
# Run inversion for the BG Compass model
#
################################################################################

### Module loading
using LinearAlgebra, Dates, Distributed, Images, Random
using JUDI, JUDI.SLIM_optim, JLD2

R = Float32
### Load synthetic data
base_path = dirname(pathof(JUDI))*"/../examples/twri/"
try
    @load string(base_path*"data/BGCompass_data_tti.jld") model_true dat fsrc
catch e
    @info "Data not found, modeling true data"
    include(base_path*"data/gen_data_bg_tti.jl")
    @load string(base_path*"data/BGCompass_data_tti.jld") model_true dat fsrc
end
# Different sigmas for smoothing
vvar = [20, 25, 30, 40]
# Global parameters
idx_w = 17
vare = 15
fevals = 20
batchsize = 25

# Bound projection
function proj_bounds(m, mmin, mmax)
    m[m[:] .< mmin] .= mmin
    m[m[:] .> mmax] .= mmax
    return m
end
mmin = 1f0/4.5f0^2
mmax = 1f0/1.4f0^2
ProjBound(x) = proj_bounds(x, mmin,  mmax)

for var=vvar
    for anis=["tt", "st"]
        ### Background model
        model0 = deepcopy(model_true)
        model0.m[:, idx_w+1:end] = R.(imfilter(model0.m[:, idx_w+1:end], Kernel.gaussian(var)))
        if anis == "st"
            model0.epsilon[:, idx_w+1:end] = R.(imfilter(model0.epsilon[:, idx_w+1:end], Kernel.gaussian(vare)))
            model0.delta[:, idx_w+1:end] = R.(imfilter(model0.delta[:, idx_w+1:end], Kernel.gaussian(vare)))
            model0.theta[:, idx_w+1:end] = R.(imfilter(model0.theta[:, idx_w+1:end], Kernel.gaussian(vare)))
        end
        m0 = model0.m

        ### Set objective functional
        mask = zeros(Float32, model0.n)
        mask[:, idx_w+1:end] .= 1

        # [TWRIdual]
        ε = zeros(Float32, fsrc.nsrc)

        ## Run
        # Very basic approximate inverse Hessian
        DZ = judiDepthScaling(model0)

        # Optimization parameters
        optwri = TWRIOptions(;grad_corr=false, comp_alpha=false, weight_fun=nothing, eps=ε, params=:m)

        # Objective function for minConf library
        function objective_function(x, fun="fwi")
            flush(Base.stdout)
            x = reshape(x, model0.n)
            # fwi function value and gradient
            i = randperm(dat.nsrc)[1:batchsize]
            # weight_fun_pars=weight_fun_pars
            if fun == "wri"
                fval, grad = twri_objective(model0, fsrc[i], dat[i], nothing; optionswri=optwri)
            else
                fval, grad = fwi_objective(model0, fsrc[i], dat[i])
            end
            grad = DZ*(vec(mask) .* vec(grad))
            grad = .1f0 .* grad ./ maximum(abs.(grad))  # scale for line search
            return fval, grad
        end

        wri_fun(x) = objective_function(x, "wri")
        fwi_fun(x) = objective_function(x, "fwi")

        # FWI with SPG
        x0 = vec(m0)
        options = spg_options(verbose=3, maxIter=fevals, memory=3, interp=0)

        sol = spg(wri_fun, x0, ProjBound, options)
        xwri, hwri = sol.x, sol.ϕ_trace
        @save string(base_path*"data/wriw_$(anis)_$(var).jld") xwri hwri

        x0 = vec(m0)
        sol = spg(fwi_fun, x0, ProjBound, options)
        xfwi, hfwi = sol.x, sol.ϕ_trace
        @save string(base_path*"data/fwiw_$(anis)_$(var).jld") xfwi hfwi
    end
end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
struct Point2D
    x::Float64
    y::Float64
end

max_stars_repo_name: JuliaTagBot/Cincinnati

max_issues_repo_name: JuliaTagBot/Cincinnati

max_forks_repo_name: JuliaTagBot/Cincinnati