text:
@testset "Zero" begin
    @test iszero(LogFixPoint16(0))
    @test iszero(LogFixPoint16(-0))
    @test LogFixPoint16(0) == LogFixPoint16(-0)

    @test LogFixPoint16(0) == zero(LogFixPoint16)
    @test 0x0000 == reinterpret(UInt16,zero(LogFixPoint16))
    @test 0x0000 == reinterpret(UInt16,-zero(LogFixPoint16))

    @test iszero(zero(LogFixPoint16))
end

@testset "NaR" begin
    @test reinterpret(UInt16,nan(LogFixPoint16)) == 0x8000
    @test isnan(nan(LogFixPoint16))
end

@testset "signbit" begin
    @test signbit(LogFixPoint16(-1))
    @test signbit(LogFixPoint16(-10))
    @test ~signbit(LogFixPoint16(0))
    @test ~signbit(LogFixPoint16(1))
    @test ~signbit(floatmax(LogFixPoint16))
    @test ~signbit(floatmin(LogFixPoint16))
end

@testset "floatmin/max" begin
    @test floatmin(LogFixPoint16) < floatmax(LogFixPoint16)
    @test zero(LogFixPoint16) < floatmin(LogFixPoint16)
    @test prevfloat(floatmin(LogFixPoint16)) == zero(LogFixPoint16)
    @test isnan(nextfloat(floatmax(LogFixPoint16)))
end

@testset "one" begin
    @test one(LogFixPoint16) == LogFixPoint16(1)
    @test reinterpret(UInt16,one(LogFixPoint16)) == 0x4000
end

@testset "negation" begin
    @test -(-(one(LogFixPoint16))) == one(LogFixPoint16)
    @test -zero(LogFixPoint16) == zero(LogFixPoint16)
    @test -nan(LogFixPoint16) != nan(LogFixPoint16)
    @test -LogFixPoint16(1) == LogFixPoint16(-1)
end

@testset "inverse" begin
    @test inv(one(LogFixPoint16)) == one(LogFixPoint16)
    @test isnan(inv(zero(LogFixPoint16)))
    @test inv(LogFixPoint16(2)) == LogFixPoint16(0.5)
    @test inv(LogFixPoint16(-2)) == LogFixPoint16(-0.5)
    @test inv(LogFixPoint16(4)) == LogFixPoint16(0.25)
    @test isnan(inv(nan(LogFixPoint16)))
end

@testset "conversion" begin
    @test Float32(LogFixPoint16(1)) == 1f0
    @test Float32(LogFixPoint16(1f0)) == 1f0
    @test Float32(LogFixPoint16(1.0)) == 1f0

    @test Float64(LogFixPoint16(1)) == 1.0
    @test Float64(LogFixPoint16(1f0)) == 1.0
    @test Float64(LogFixPoint16(1.0)) == 1.0

    @test Float32(LogFixPoint16(0)) == 0f0
    @test Float32(LogFixPoint16(-0)) == 0f0
    @test Float32(LogFixPoint16(-1)) == -1f0
    @test Float32(LogFixPoint16(-8)) == -8f0
    @test Float32(LogFixPoint16(-0.125)) == -0.125f0

    @test Int(LogFixPoint16(4)) == 4
    @test Int(LogFixPoint16(-2)) == -2

    @test isnan(LogFixPoint16(NaN32))
    @test isnan(LogFixPoint16(NaN))
    @test isnan(LogFixPoint16(Inf))
    @test isnan(LogFixPoint16(-Inf))
end

@testset "conversion to Float32 all" begin
    for hex in 0x0000:0xffff
        lf = reinterpret(LogFixPoint16,hex)
        lf2f2lf = LogFixPoint16(Float32(lf))
        @test reinterpret(UInt16,lf) == reinterpret(UInt16,lf2f2lf)
        # @test lf == lf2f2lf   # fails for NaR == NaR
    end
end

@testset "multiplication" begin
    @test LogFixPoint16(1)*LogFixPoint16(2) == LogFixPoint16(2)
    @test LogFixPoint16(123)*inv(LogFixPoint16(123)) == one(LogFixPoint16)
    @test Float32(LogFixPoint16(-8)*LogFixPoint16(8)) == -64f0
    @test Float32(LogFixPoint16(-8)*LogFixPoint16(-0.125)) == 1f0
    @test iszero(LogFixPoint16(0)*LogFixPoint16(12.1234))
    @test iszero(LogFixPoint16(0)*LogFixPoint16(-12.1234))
end

@testset "multiplication all" begin
    for _ in 1:1000000
        lf1,lf2 = rand(UInt16,2)
        lf1 = reinterpret(LogFixPoint16,lf1)
        lf2 = reinterpret(LogFixPoint16,lf2)

        result = LogFixPoint16(Float32(lf1)*Float32(lf2))
        result2 = lf1*lf2

        # result might yield NaN in Float32 arithmetic (overflow)
        # but in LogFixPoint16-arithmetic yield ±floatmax (no-overflow)
        # in case Float32 yields NaN test that LogFixPoint16 is either
        # NaR or ±floatmax
        if isnan(result)
            if isnan(result2)
                @test true
            else
                @test (result2 == floatmax(LogFixPoint16)) || (result2 == -floatmax(LogFixPoint16))
            end
        else
            @test result == result2
        end
    end
end

@testset "division" begin
    @test one(LogFixPoint16) / one(LogFixPoint16) == one(LogFixPoint16)
    @test LogFixPoint16(-0.1273) / LogFixPoint16(0.1273) == -one(LogFixPoint16)
    @test LogFixPoint16(4) / LogFixPoint16(2) == LogFixPoint16(2)
    @test isnan(LogFixPoint16(3)/zero(LogFixPoint16))
    @test LogFixPoint16(3.123)/LogFixPoint16(0.1232) ==
            LogFixPoint16(3.123)*inv(LogFixPoint16(0.1232))
end

@testset "division all" begin
    for _ in 1:1000000
        lf1,lf2 = rand(UInt16,2)
        lf1 = reinterpret(LogFixPoint16,lf1)
        lf2 = reinterpret(LogFixPoint16,lf2)

        result = LogFixPoint16(Float32(lf1)/Float32(lf2))
        result2 = lf1/lf2

        # result might yield NaN in Float32 arithmetic (overflow)
        # but in LogFixPoint16-arithmetic yield ±floatmax (no-overflow)
        # in case Float32 yields NaN test that LogFixPoint16 is either
        # NaR or ±floatmax
        if isnan(result)
            if isnan(result2)
                @test true
            else
                @test (result2 == floatmax(LogFixPoint16)) || (result2 == -floatmax(LogFixPoint16))
            end
        else
            @test result == result2
        end
    end
end

@testset "sqrt" begin
    @test sqrt(LogFixPoint16(4)) == LogFixPoint16(2)
    @test_skip sqrt(LogFixPoint16(12.56))^2 == prevfloat(LogFixPoint16(12.56))
    @test isnan(sqrt(LogFixPoint16(-1)))
    @test iszero(sqrt(LogFixPoint16(0)))
end

@testset "addition" begin
    a,b,c,d,e = LogFixPoint16.([-8,-4,-2,0,1])

    @test a+b == b+a
    @test a+d == a
    @test e+e == -c
    @test -(c+c) == -b
    @test d+e == e+d
end

@testset "addition all" begin
    for _ in 1:100000
        lf1ui,lf2ui = rand(UInt16,2)

        # abs to only test the addTable
        lf1 = abs(reinterpret(LogFixPoint16,lf1ui))
        lf2 = abs(reinterpret(LogFixPoint16,lf2ui))

        result = LogFixPoint16(Float32(lf1)+Float32(lf2))
        @test_skip reinterpret(UInt16,result) == reinterpret(UInt16,lf1+lf2)
    end
end

@testset "max diff resolvable" begin
    @test LogFixPoint16s.addTable[LogFixPoint16s.max_diff_resolvable[]] == 1
    if LogFixPoint16s.nfrac[] < 11
        @test LogFixPoint16s.addTable[LogFixPoint16s.max_diff_resolvable[]+1] == 0
    end
end

@testset "subtraction all" begin
    for _ in 1:100000
        lf1ui,lf2ui = rand(UInt16,2)

        # abs to only test the addTable
        lf1 = abs(reinterpret(LogFixPoint16,lf1ui))
        lf2 = abs(reinterpret(LogFixPoint16,lf2ui))

        result = LogFixPoint16(Float32(lf1)-Float32(lf2))
        @test_skip reinterpret(UInt16,result) == reinterpret(UInt16,lf1-lf2)
    end
end

@testset "addition with zero" begin
    lf2 = zero(LogFixPoint16)
    for hex in 0x0000:0xffff  # only positive numbers
        lf1 = reinterpret(LogFixPoint16,hex)
        @test hex == reinterpret(UInt16,lf1+lf2)
        @test hex == reinterpret(UInt16,lf2+lf1)
    end
end

@testset "subtraction with zero" begin
    lf2 = zero(LogFixPoint16)
    for hex in 0x0000:0xffff  # only positive numbers
        lf1 = reinterpret(LogFixPoint16,hex)
        @test hex == reinterpret(UInt16,lf1-lf2)
    end
end

@testset "subtraction" begin
    a,b,c,d,e = LogFixPoint16.([-8,-4,-2,0,1])

    @test a-b == -(b-a)
    @test c-a == -(a-c)
    @test d-a == -a
    @test a-a == d
end

@testset "nextfloat prevfloat" begin
    a = floatmin(LogFixPoint16)
    b = floatmax(LogFixPoint16)
    c = LogFixPoint16(123)

    @test nextfloat(prevfloat(c)) == c
    @test prevfloat(nextfloat(c)) == c
    @test prevfloat(c) < c
    @test prevfloat(-c) > -c
    @test nextfloat(c) > c
    @test nextfloat(-c) < -c

    @test iszero(prevfloat(a))
    @test iszero(nextfloat(-a))
    @test isnan(nextfloat(b))
    @test isnan(prevfloat(-b))

    @test nextfloat(nextfloat(c)) == nextfloat(c,2)
    @test prevfloat(prevfloat(c)) == prevfloat(c,2)
end

@testset "bitstring" begin
    lf = one(LogFixPoint16)
    ni = LogFixPoint16s.nint[]
    nf = LogFixPoint16s.nfrac[]
    @test bitstring(lf) == "0100000000000000"
    @test bitstring(lf,:split) == "0 1"*repeat("0",ni-1)*" "*repeat("0",nf)
end

@testset "power2" begin
    @test power2(LogFixPoint16(2)) == 2*LogFixPoint16(2)
    @test power2(zero(LogFixPoint16)) == zero(LogFixPoint16)
    @test isequal(power2(nan(LogFixPoint16)),nan(LogFixPoint16))
    @test power2(-one(LogFixPoint16)) == one(LogFixPoint16)
    @test power2(LogFixPoint16(0.5)) == LogFixPoint16(0.25)
end


text:
"""
Caluclates and returns the admittance laplacian of the powergrid pg.
The notation was choosen in accordance to: https://en.wikipedia.org/wiki/Nodal_admittance_matrix
Inputs:
        pg: Powergrid, a graph contaning nodes and lines
"""
function admittance_laplacian(pg::PowerGrid)
    Y =  zeros(Complex, length(pg.nodes), length(pg.nodes))
    y_sum = zeros(Complex, length(pg.nodes))
    for i in 1:length(pg.lines)
        From = pg.lines[i].from
        To = pg.lines[i].to
        Y[From,To] = -1 * pg.lines[i].Y
        y_sum[From] += pg.lines[i].Y
    end
    for j in 1:length(pg.nodes)
        Y[j,j] = pg.nodes[j].Y_n + y_sum[j]
    end
    return Y
end

"""
A Simple Power Flow Calculation for PowerDynamics.
The calculation is only correct for nodes where the phase θ,
voltage magnitude u are know. For powergrids contaning eg. SlackAlgebraic the
calculation can not be performed because the phases θ are unknown.
A Power-flow study would need to be performed before hand.
Inputs:
    pg: PowerGrid, a graph containg nodes and lines
    op: State, Operationpoint (or any other valid point) of the powergrid pg
Outputs:
    P: Vector, Active Power at the nodes
    Q: Vector, Reactive Power at the nodes
"""
function SimplePowerFlow(pg::PowerGrid, op::State)
    #if SlackAlgebraic ∈ pg.nodes .|> typeof
    #    throw(OperationPointError("Found a SlackAlgebraic in the powergrid."))
    #end
    Y = admittance_laplacian(pg)
    G = real(Y)
    B = imag(Y)

    P = zeros(Complex, length(pg.nodes))
    Q = zeros(Complex, length(pg.nodes))

    θ_idx = findall(:θ .∈ symbolsof.(pg.nodes))
    for i in θ_idx
        Vi = op[i, :v] # Voltage Magnitude of node i
        for k in θ_idx
            Vk = op[k,:v] # Volatge Magnitude of node k
            θik = op[i, :θ] - op[k, :θ] # Phase difference between node i and k
            P[i] += Vk * (G[i,k] * sin(θik) + B[i,k] * sin(θik))
            Q[i] += Vk * (G[i,k] * cos(θik) - B[i,k] * cos(θik))
        end
        P[i] *= Vi
        Q[i] *= Vi
    end
    return P, Q
end

ext: jl

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
using Distributions, JuMP, Ipopt, Random
include("KPNormal.jl")
using KP

function rob(cs, muhat, vs, r)
    const sqrt_vmin = sqrt(minimum(vs))
    function f(Gamma)
        xs, lam = KP.x_l2reg(cs, muhat, vs, Gamma)
#        println(Gamma, "\t", Gamma * sqrt_vmin * v_norm(xs, vs) - r)
        return Gamma * sqrt_vmin * v_norm(xs, vs) - r
    end
    println(f(.01), "\t", f(100))
    Gammastar = fzero(f, .01, 100)
    return KP.x_l2reg(cs, muhat, vs, Gammastar)
end

function rob_direct(cs, muhat, vs, r)
    m = Model(solver=IpoptSolver())
    @variable(m, xs[1:n] >= 0)
    @constraint(m, sum(cs[i]* xs[i] for i = 1:n) <= n)
    @NLobjective(m, Max, sum(muhat[i] * xs[i]/n for i = 1:n) - r/n * sqrt( sum(xs[i]^2 for i = 1:n)))
    solve(m)
end

v_norm(xs, vs) = sqrt(sum(xs.^2 ./ vs))

n =  1000
Random.seed!(8675309)
thetas = ones(n)
vs = ones(n)
cs = ones(n)

#From above
theta_l, v_l, c_l = 0.0, .01,  10.001
theta_m, v_m, c_m = 1.0, .75, 10.001
theta_h, v_h, c_h = 0.3,  2, 10.001

thetas[1:3:n] = theta_l
vs[1:3:n] = v_l
cs[1:3:n] = c_l

thetas[2:3:n] = theta_m
vs[2:3:n] = v_m
cs[2:3:n] = c_m

thetas[3:3:n] = theta_h
vs[3:3:n] = v_h
cs[3:3:n] = c_h

muhat = thetas + randn(n) ./ sqrt(vs)
sqrt_vmin = sqrt(minimum(vs))

tic()
xs, lam = rob(cs, muhat, vs, quantile(Normal(0, 1), .99))
toc()

tic()
xs2 = rob_direct(cs, muhat, vs, quantile(Normal(0, 1), .99))
toc()


text:
using Test
using LinearAlgebra

@testset "Moreau envelope" begin

@testset "Box indicator" for R in [Float32, Float64]

    f = IndBox(-1, 1)

    for g in [
        MoreauEnvelope(f),
        MoreauEnvelope(f, R(0.01))
    ]

        predicates_test(g)

        @test ProximalOperators.is_smooth(g) == true
        @test ProximalOperators.is_quadratic(g) == false
        @test ProximalOperators.is_set(g) == false

        x = R[1.0, 2.0, 3.0, 4.0, 5.0]

        grad_g_x, g_x = gradient(g, x)

        y, g_y = prox_test(g, x, R(1/2))
        grad_g_y, _ = gradient(g, y)

        @test y + grad_g_y / 2 ≈ x
        @test g(y) ≈ g_y
    end

end

@testset "L2 norm" for R in [Float32, Float64]

    for (g, h) in [
        (MoreauEnvelope(NormL2()), HuberLoss()),
        (MoreauEnvelope(NormL2(R(1)), R(0.1)), HuberLoss(R(0.1), R(1)/R(0.1)))
    ]

        predicates_test(g)

        @test ProximalOperators.is_smooth(g) == true
        @test ProximalOperators.is_quadratic(g) == false
        @test ProximalOperators.is_set(g) == false

        x = R[1.0, 2.0, 3.0, 4.0, 5.0]

        @test g(x) ≈ h(x)

        grad_g_x, g_x = gradient(g, x)
        grad_h_x, h_x = gradient(h, x)

        @test g_x ≈ g(x)
        @test h_x ≈ h(x)
        @test all(grad_g_x .≈ grad_h_x)

        y, g_y = prox_test(g, x, R(1/2))
        grad_g_y, _ = gradient(g, y)

        @test y + grad_g_y / 2 ≈ x
        @test g(y) ≈ g_y

    end

end

end


text:
module AlgoritmiStruttureDati

include("./laboratorio3/Laboratorio3.jl")

export laboratorio3

end # module


text:
"""
    PeriodicDerivativeCoefficients

The coefficients of a derivative operator on a periodic grid.
"""
struct PeriodicDerivativeCoefficients{T,LowerOffset,UpperOffset,ExecutionMode,SourceOfCoefficients} <: AbstractDerivativeCoefficients{T}
    # coefficients defining the operator and its action
    lower_coef::SVector{LowerOffset, T}
    central_coef::T
    upper_coef::SVector{UpperOffset, T}
    mode::ExecutionMode
    # corresponding orders etc.
    derivative_order::Int
    accuracy_order  ::Int
    symmetric       ::Bool
    source_of_coefficients::SourceOfCoefficients

    function PeriodicDerivativeCoefficients(
            lower_coef::SVector{LowerOffset, T}, central_coef::T, upper_coef::SVector{UpperOffset, T},
            mode::ExecutionMode, derivative_order::Int, accuracy_order::Int,
            source_of_coefficients::SourceOfCoefficients) where {T,LowerOffset,UpperOffset,ExecutionMode,SourceOfCoefficients}
        symmetric = LowerOffset == UpperOffset
        if symmetric
            @inbounds for i in Base.OneTo(LowerOffset)
                symmetric = symmetric && lower_coef[i] == upper_coef[i]
            end
        end
        new{T,LowerOffset,UpperOffset,ExecutionMode,SourceOfCoefficients}(
            lower_coef, central_coef, upper_coef, mode,
            derivative_order, accuracy_order, symmetric, source_of_coefficients)
    end
end

# Compute `α*D*u + β*dest` and store the result in `dest`.
function mul!(dest::AbstractVector, coefficients::PeriodicDerivativeCoefficients{T,LowerOffset,UpperOffset}, u::AbstractVector, α, β) where {T,LowerOffset,UpperOffset}
    @boundscheck begin
        @argcheck length(dest) == length(u) DimensionMismatch
        @argcheck length(u) > LowerOffset+UpperOffset DimensionMismatch
    end

    @unpack lower_coef, central_coef, upper_coef, mode = coefficients
    convolve_periodic_boundary_coefficients!(dest, lower_coef, central_coef, upper_coef, u, α, β)
    convolve_interior_coefficients!(dest, lower_coef, central_coef, upper_coef, u, α, β, static_length(lower_coef), static_length(upper_coef), mode)
end

# Compute `α*D*u` and store the result in `dest`.
function mul!(dest::AbstractVector, coefficients::PeriodicDerivativeCoefficients{T,LowerOffset,UpperOffset}, u::AbstractVector, α) where {T,LowerOffset,UpperOffset}
    @boundscheck begin
        @argcheck length(dest) == length(u) DimensionMismatch
        @argcheck length(u) > LowerOffset+UpperOffset DimensionMismatch
    end

    @unpack lower_coef, central_coef, upper_coef, mode = coefficients
    convolve_periodic_boundary_coefficients!(dest, lower_coef, central_coef, upper_coef, u, α)
    convolve_interior_coefficients!(dest, lower_coef, central_coef, upper_coef, u, α, static_length(lower_coef), static_length(upper_coef), mode)
end

@generated function convolve_periodic_boundary_coefficients!(dest::AbstractVector, lower_coef::SVector{LowerOffset}, central_coef, upper_coef::SVector{UpperOffset},
                                                             u::AbstractVector, α, β) where {LowerOffset, UpperOffset}
    ex_lower = :( nothing )
    for i in 1:LowerOffset
        ex = :( lower_coef[$LowerOffset]*u[end-$(LowerOffset-i)] )
        for j in LowerOffset-1:-1:1
            if i-j < 1
                ex = :( $ex + lower_coef[$j]*u[end-$(j-i)] )
            else
                ex = :( $(ex) + lower_coef[$j]*u[$(i-j)] )
            end
        end
        ex = :( $ex + central_coef*u[$i] )
        for j in 1:UpperOffset
            ex = :( $ex + upper_coef[$j]*u[$(i+j)] )
        end
        ex_lower = quote
            $ex_lower
            @inbounds dest[$i] = β*dest[$i] + α*$ex
        end
    end

    ex_upper = :( nothing )
    for i in (UpperOffset-1):-1:0
        ex = :( lower_coef[$LowerOffset]*u[end-$(i+LowerOffset)] )
        for j in LowerOffset-1:-1:1
            ex = :( $ex + lower_coef[$j]*u[end-$(j+i)] )
        end
        ex = :( $ex + central_coef*u[end-$i] )
        for j in 1:UpperOffset
            if i-j < 0
                ex = :( $ex + upper_coef[$j]*u[$(j-i)] )
            else
                ex = :( $ex + upper_coef[$j]*u[end-$(i-j)] )
            end
        end
        ex_upper = quote
            $ex_upper
            @inbounds dest[end-$i] = β*dest[end-$i] + α*$ex
        end
    end

    quote
        $ex_lower
        $ex_upper
    end
end

@generated function convolve_periodic_boundary_coefficients!(dest::AbstractVector, lower_coef::SVector{LowerOffset}, central_coef, upper_coef::SVector{UpperOffset},
                                                             u::AbstractVector, α) where {LowerOffset, UpperOffset}
    ex_lower = :( nothing )
    for i in 1:LowerOffset
        ex = :( lower_coef[$LowerOffset]*u[end-$(LowerOffset-i)] )
        for j in LowerOffset-1:-1:1
            if i-j < 1
                ex = :( $ex + lower_coef[$j]*u[end-$(j-i)] )
            else
                ex = :( $(ex) + lower_coef[$j]*u[$(i-j)] )
            end
        end
        ex = :( $ex + central_coef*u[$i] )
        for j in 1:UpperOffset
            ex = :( $ex + upper_coef[$j]*u[$(i+j)] )
        end
        ex_lower = quote
            $ex_lower
            @inbounds dest[$i] = α*$ex
        end
    end

    ex_upper = :( nothing )
    for i in (UpperOffset-1):-1:0
        if LowerOffset > 0
            ex = :( lower_coef[$LowerOffset]*u[end-$(i+LowerOffset)] )
            for j in LowerOffset-1:-1:1
                ex = :( $ex + lower_coef[$j]*u[end-$(j+i)] )
            end
            ex = :( $ex + central_coef*u[end-$i] )
        else
            ex = :( central_coef*u[end-$i] )
        end
        for j in 1:UpperOffset
            if i-j < 0
                ex = :( $ex + upper_coef[$j]*u[$(j-i)] )
            else
                ex = :( $ex + upper_coef[$j]*u[end-$(i-j)] )
            end
        end
        ex_upper = quote
            $ex_upper
            @inbounds dest[end-$i] = α*$ex
        end
    end

    quote
        $ex_lower
        $ex_upper
    end
end

"""
    BeljaddLeFlochMishraParés2017()

Coefficients of the periodic operators given in
- Beljadid, LeFloch, Mishra, Parés (2017)
  Schemes with Well-Controlled Dissipation. Hyperbolic Systems in
    Nonconservative Form.
  Communications in Computational Physics 21.4, pp. 913-946.
"""
struct BeljaddLeFlochMishraParés2017 <: SourceOfCoefficients end

function Base.show(io::IO, source::BeljaddLeFlochMishraParés2017)
    if get(io, :compact, false)
        summary(io, source)
    else
        print(io,
            "Beljadid, LeFloch, Mishra, Parés (2017) \n",
            "  Schemes with Well-Controlled Dissipation. Hyperbolic Systems in \n",
            "    Nonconservative Form. \n",
            "  Communications in Computational Physics 21.4, pp. 913-946.")
    end
end

"""
    periodic_central_derivative_coefficients(derivative_order, accuracy_order, T=Float64, mode=FastMode())

Create the `PeriodicDerivativeCoefficients` approximating the `derivative_order`-th
derivative with an order of accuracy `accuracy_order` and scalar type `T`.
The evaluation of the derivative can be parallized using threads by chosing
`mode=ThreadedMode())`.
"""
function periodic_central_derivative_coefficients(derivative_order, accuracy_order, T=Float64, mode=FastMode())
    @argcheck derivative_order > 0
    @argcheck accuracy_order > 0

    if isodd(accuracy_order)
        @warn("Central derivative operators support only even orders of accuracy.")
    end

    if derivative_order == 1
        # Exact evaluation of the coefficients, see
        # Beljadid, LeFloch, Mishra, Parés (2017)
        # Schemes with Well-Controlled Dissipation. Hyperbolic Systems in Nonconservative Form.
        # Communications in Computational Physics 21.4, pp. 913-946.
        n = (accuracy_order+1) ÷ 2
        coef = Vector{T}(undef, n)
        for j in 1:n
            tmp = one(Rational{Int128})
            if iseven(j)
                tmp *= -1
            end
            for i in 1:j
                tmp *= (n+1-i) // (n+i)
            end
            coef[j] =  tmp / j
        end
        upper_coef = SVector{n,T}(coef)
        central_coef = zero(T)
        lower_coef = -upper_coef
        source = BeljaddLeFlochMishraParés2017()
    elseif derivative_order == 2
        # Exact evaluation of the coefficients, see
        # Beljadid, LeFloch, Mishra, Parés (2017)
        # Schemes with Well-Controlled Dissipation. Hyperbolic Systems in Nonconservative Form.
        # Communications in Computational Physics 21.4, pp. 913-946.
        n = (accuracy_order+1) ÷ 2
        coef = Vector{T}(undef, n)
        for j in 1:n
            tmp = one(Rational{Int128})
            if iseven(j)
                tmp *= -1
            end
            for i in 1:j
                tmp *= (n+1-i) // (n+i)
            end
            coef[j] =  2*tmp / j^2
        end
        upper_coef = SVector{n,T}(coef)
        central_coef = -2*sum(upper_coef)
        lower_coef = upper_coef
        source = BeljaddLeFlochMishraParés2017()
    elseif derivative_order == 3
        # Exact evaluation of the coefficients, see
        # Beljadid, LeFloch, Mishra, Parés (2017)
        # Schemes with Well-Controlled Dissipation. Hyperbolic Systems in Nonconservative Form.
        # Communications in Computational Physics 21.4, pp. 913-946.
        n = (accuracy_order+1) ÷ 2
        coef = Vector{T}(undef, n)
        for j in 1:n
            tmp = one(Rational{Int128})
            if isodd(j)
                tmp *= -1
            end
            for i in 1:j
                tmp *= (n+1-i) // (n+i)
            end
            fac = zero(Rational{Int128})
            for k in 1:n
                k == j && continue
                fac += 1 // k^2
            end
            coef[j] =  6*tmp*fac/ j
        end
        upper_coef = SVector{n,T}(coef)
        central_coef = zero(T)
        lower_coef = -upper_coef
        source = BeljaddLeFlochMishraParés2017()
    else
        throw(ArgumentError("Derivative order $derivative_order not implemented yet."))
    end

    PeriodicDerivativeCoefficients(lower_coef, central_coef, upper_coef, mode, derivative_order, accuracy_order, source)
end

"""
    Fornberg1998()

Coefficients of the periodic operators given in
- Fornberg (1998)
  Calculation of Weights in Finite Difference Formulas.
  SIAM Rev. 40.3, pp. 685-691.
"""
struct Fornberg1998 <: SourceOfCoefficients end

function Base.show(io::IO, source::Fornberg1998)
    if get(io, :compact, false)
        summary(io, source)
    else
        print(io,
            "Fornberg (1998) \n",
            "  Calculation of Weights in Finite Difference Formulas. \n",
            "  SIAM Rev. 40.3, pp. 685-691.")
    end
end

"""
    fornberg(x::Vector{T}, m::Int) where {T}

Calculate the weights of a finite difference approximation of the `m`th derivative
with maximal order of accuracy at `0` using the nodes `x`, see
Fornberg (1998)
Calculation of Weights in Finite Difference Formulas
SIAM Rev. 40.3, pp. 685-691.
"""
function fornberg(xx::Vector{T}, m::Int) where {T}
    x = sort(xx)
    @argcheck x[1] <= 0
    @argcheck x[end] >= 0
    @argcheck length(unique(x)) == length(x)
    @argcheck zero(T) in x
    @argcheck m >= 1

    z = zero(T)
    n = length(x) - 1
    c = fill(zero(T), length(x), m+1)
    c1 = one(T)
    c4 = x[1] - z
    c[1,1] = one(T)
    for i in 1:n
        mn = min(i,m)
        c2 = one(T)
        c5 = c4
        c4 = x[i+1] - z
        for j in 0:i-1
            c3 = x[i+1] - x[j+1]
            c2 = c2 * c3
            if j == i-1
                for k in mn:-1:1
                    c[i+1,k+1] = c1 * (k*c[i,k]-c5*c[i,k+1]) / c2
                end
                c[i+1,1] = -c1*c5*c[i,1] / c2
            end
            for k in mn:-1:1
                c[j+1,k+1] = (c4*c[j+1,k+1]-k*c[j+1,k]) / c3
            end
            c[j+1,1] = c4*c[j+1,1] / c3
        end
        c1 = c2
    end

    c[:,end]
end

"""
    periodic_derivative_coefficients(derivative_order, accuracy_order,
                                     left_offset=-(accuracy_order+1)÷2,
                                     T=Float64, mode=FastMode())

Create the `PeriodicDerivativeCoefficients` approximating the `derivative_order`-th
derivative with an order of accuracy `accuracy_order` and scalar type `T` where
the leftmost grid point used is determined by `left_offset`.
The evaluation of the derivative can be parallized using threads by chosing
mode=ThreadedMode()`.
"""
function periodic_derivative_coefficients(derivative_order, accuracy_order,
                                          left_offset::Int=-(accuracy_order+1)÷2,
                                          T=Float64, mode=FastMode())
    @argcheck -accuracy_order <= left_offset <= 0

    x = Rational{Int128}.(left_offset:left_offset+accuracy_order)
    c = fornberg(x, derivative_order)

    z = zero(eltype(x))
    lower_idx = x .< z
    central_idx = first(findall(xx->xx==z, x))
    upper_idx = x .> z

    LowerOffset = sum(lower_idx)
    UpperOffset = sum(upper_idx)

    lower_coef = SVector{LowerOffset, T}(T.(reverse(c[lower_idx])))
    central_coef = T(c[central_idx])
    upper_coef = SVector{UpperOffset, T}(T.(c[upper_idx]))

    source = Fornberg1998()

    PeriodicDerivativeCoefficients(lower_coef, central_coef, upper_coef, mode, derivative_order, accuracy_order, source)
end

"""
    Holoborodko2008()

Coefficients of the periodic operators given in
- Holoborodko (2008)
  Smooth Noise Robust Differentiators.
  http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
"""
struct Holoborodko2008 <: SourceOfCoefficients end

function Base.show(io::IO, source::Holoborodko2008)
    if get(io, :compact, false)
        summary(io, source)
    else
        print(io,
            "  Holoborodko (2008) \n",
            "  Smooth Noise Robust Differentiators. \n",
            "  http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/")
    end
end

"""
    periodic_derivative_coefficients(source::Holoborodko2008, derivative_order, accuracy_order;
                                     T=Float64, mode=FastMode(),
                                     stencil_width=accuracy_order+3)

Create the `PeriodicDerivativeCoefficients` approximating the `derivative_order`-th
derivative with an order of accuracy `accuracy_order` and scalar type `T` given
by [`Holoborodko2008`](@ref).
The evaluation of the derivative can be parallized using threads by chosing
mode=ThreadedMode()`.
"""
function periodic_derivative_coefficients(source::Holoborodko2008, derivative_order, accuracy_order;
                                          T=Float64, mode=FastMode(),
                                          stencil_width=accuracy_order+3)
    method_exists = true
    if derivative_order == 1
        if accuracy_order == 2
            if stencil_width == 5
                lower_coef = SVector{2, T}([
                    -2//8, -1//8
                ])
                upper_coef = -lower_coef
                central_coef = T(0)
            elseif stencil_width == 7
                lower_coef = SVector{3, T}([
                    -5//32, -4//32, -1//32
                ])
                upper_coef = -lower_coef
                central_coef = T(0)
            elseif stencil_width == 9
                lower_coef = SVector{4, T}([
                    -14//128, -14//128, -6//128, -1//128
                ])
                upper_coef = -lower_coef
                central_coef = T(0)
            elseif stencil_width == 11
                lower_coef = SVector{5, T}([
                    -42//512, -48//512, -27//512, -8//512, -1//512
                ])
                upper_coef = -lower_coef
                central_coef = T(0)
            else
                method_exists = false
            end
        elseif accuracy_order == 4
            if stencil_width == 7
                lower_coef = SVector{3, T}([
                    -39//96, -12//96, 5//96
                ])
                upper_coef = -lower_coef
                central_coef = T(0)
            elseif stencil_width == 9
                lower_coef = SVector{4, T}([
                    -27//96, -16//96, 1//96, 2//96
                ])
                upper_coef = -lower_coef
                central_coef = T(0)
            elseif stencil_width == 11
                lower_coef = SVector{5, T}([
                    -322//1536, -256//1536, -39//1536, 32//1536, 11//1536
                ])
                upper_coef = -lower_coef
                central_coef = T(0)
            else
                method_exists = false
            end
        else
            method_exists = false
        end
    elseif derivative_order == 2
        if accuracy_order == 2
            if stencil_width == 5
                lower_coef = SVector{2, T}([
                    0, 1//4
                ])
                upper_coef = lower_coef
                central_coef = T(-2//4)
            elseif stencil_width == 7
                lower_coef = SVector{3, T}([
                    -1//16, 2//16, 1//16
                ])
                upper_coef = lower_coef
                central_coef = T(-4//16)
            elseif stencil_width == 9
                lower_coef = SVector{4, T}([
                    -4//64, 4//64, 4//64, 1//64
                ])
                upper_coef = lower_coef
                central_coef = T(-10//64)
            else
                method_exists = false
            end
        elseif accuracy_order == 4
            if stencil_width == 7
                lower_coef = SVector{3, T}([
                    1//12, 5//12, -1//12
                ])
                upper_coef = lower_coef
                central_coef = T(-10//12)
            elseif stencil_width == 9
                lower_coef = SVector{4, T}([
                    -12//192, 52//192, 12//192, -7//192
                ])
                upper_coef = lower_coef
                central_coef = T(-90//192)
            else
                method_exists = false
            end
        else
            method_exists = false
        end
    else
        method_exists = false
    end
    if method_exists == false
        throw(ArgumentError("Method with derivative_order=$derivative_order, accuracy_order=$accuracy_order, stencil_width=$stencil_width not implemented/derived."))
    end

    PeriodicDerivativeCoefficients(lower_coef, central_coef, upper_coef, mode, derivative_order, accuracy_order, source)
end

"""
    PeriodicDerivativeOperator

A derivative operator on a uniform periodic grid.
See [`periodic_derivative_operator`](@ref) and
[`periodic_central_derivative_operator`](@ref).
"""
struct PeriodicDerivativeOperator{T,LowerOffset,UpperOffset,ExecutionMode,SourceOfCoefficients,Grid} <: AbstractPeriodicDerivativeOperator{T}
    coefficients::PeriodicDerivativeCoefficients{T,LowerOffset,UpperOffset,ExecutionMode,SourceOfCoefficients}
    grid_compute::Grid   # N   nodes, including the left and excluding the right boundary
    grid_evaluate::Grid #  N+1 nodes, including both boundaries
    Δx::T
    factor::T

    function PeriodicDerivativeOperator(coefficients::PeriodicDerivativeCoefficients{T,LowerOffset,UpperOffset,ExecutionMode,SourceOfCoefficients},
                                        grid_evaluate::Grid) where {T,LowerOffset,UpperOffset,ExecutionMode,SourceOfCoefficients,Grid}
        @argcheck length(grid_evaluate) > LowerOffset+UpperOffset DimensionMismatch
        grid_compute = range(grid_evaluate[1], stop=grid_evaluate[end-1], length=length(grid_evaluate)-1)

        Δx = step(grid_evaluate)
        factor = inv(Δx^coefficients.derivative_order)
        new{T,LowerOffset,UpperOffset,ExecutionMode,SourceOfCoefficients,Grid}(coefficients, grid_compute, grid_evaluate, Δx, factor)
    end
end

@inline source_of_coefficients(D::PeriodicDerivativeOperator) = source_of_coefficients(D.coefficients)

function Base.show(io::IO, D::PeriodicDerivativeOperator{T,LowerOffset,UpperOffset}) where {T,LowerOffset,UpperOffset}
    x = D.grid_evaluate
    if derivative_order(D) == 1
        print(io, "Periodic first-derivative operator")
    elseif  derivative_order(D) == 2
        print(io, "Periodic second-derivative operator")
    elseif  derivative_order(D) == 3
        print(io, "Periodic third-derivative operator")
    else
        print(io, "Periodic ", derivative_order(D),
              "-derivative operator")
    end
    print(io, " of order ", accuracy_order(D))
    if get(io, :compact, false) == false
        print(io, " on a grid in [", first(x), ", ", last(x),
                "] using ", size(D, 1), " nodes, \n")
        print(io, "stencils with ", LowerOffset, " nodes to the left, ",
                UpperOffset, " nodes to the right, and coefficients")
    end
    print(io, " of ", source_of_coefficients(D))
end

# Compute `α*D*u + β*dest` and store the result in `dest`.
Base.@propagate_inbounds function mul!(dest::AbstractVector, D::PeriodicDerivativeOperator, u::AbstractVector, α, β)
    @boundscheck begin
        @argcheck size(D, 2) == length(u) DimensionMismatch
        @argcheck size(D, 1) == length(dest) DimensionMismatch
    end
    @inbounds mul!(dest, D.coefficients, u, α*D.factor, β)
end

# Compute `α*D*u` and store the result in `dest`.
Base.@propagate_inbounds function mul!(dest::AbstractVector, D::PeriodicDerivativeOperator, u::AbstractVector, α)
    @boundscheck begin
        @argcheck size(D, 2) == length(u) DimensionMismatch
        @argcheck size(D, 1) == length(dest) DimensionMismatch
    end
    @inbounds mul!(dest, D.coefficients, u, α*D.factor)
end

"""
    integrate(func, u, D::PeriodicDerivativeOperator)

Map the function `func` to the coefficients `u` and integrate with respect to
the quadrature rule associated with the periodic derivative operator `D`.
"""
function integrate(func, u::AbstractVector, D::PeriodicDerivativeOperator)
    @boundscheck begin
        length(u) == length(grid(D))
    end
    @unpack Δx = D

    @inbounds res = sum(func, u)

    Δx * res
end

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

    Δx * I
end

"""
    periodic_central_derivative_operator(derivative_order, accuracy_order,
                                         xmin, xmax, N, mode=FastMode())

Create a [`PeriodicDerivativeOperator`](@ref) approximating the `derivative_order`-th
derivative on a uniform grid between `xmin` and `xmax` with `N` grid points up
to order of accuracy `accuracy_order`.
The evaluation of the derivative can be parallized using threads by chosing
`mode=ThreadedMode()`.
"""
function periodic_central_derivative_operator(derivative_order, accuracy_order,
                                              xmin, xmax, N, mode=FastMode())
    if mode === Val(:serial) || mode === Val(:threads)
        # TODO: deprecated in v0.5
        Base.depwarn("Providing the argument `parallel` is deprecated." *
                     "Use `mode` instead.", :periodic_central_derivative_operator)
        mode = _parallel_to_mode(mode)
    end
    grid_evaluate = range(xmin, stop=xmax, length=N+1) # N includes two identical boundary nodes
    coefficients = periodic_central_derivative_coefficients(derivative_order, accuracy_order, eltype(grid_evaluate), mode)
    PeriodicDerivativeOperator(coefficients, grid_evaluate)
end

"""
    periodic_derivative_operator(derivative_order, accuracy_order,
                                 xmin, xmax, N,
                                 left_offset=-(accuracy_order+1)÷2,
                                 mode=FastMode())
    periodic_derivative_operator(; derivative_order, accuracy_order,
                                 xmin, xmax, N,
                                 left_offset=-(accuracy_order+1)÷2,
                                 mode=FastMode())

Create a [`PeriodicDerivativeOperator`](@ref) approximating the `derivative_order`-th
derivative on a uniform grid between `xmin` and `xmax` with `N` grid points up
to order of accuracy `accuracy_order` where the leftmost grid point used is
determined by `left_offset`.
The evaluation of the derivative can be parallized using threads by chosing
`mode=ThreadedMode())`.

## Examples

```jldoctest
julia> periodic_derivative_operator(derivative_order=1, accuracy_order=2,
                                    xmin=0.0, xmax=1.0, N=11)
Periodic first-derivative operator of order 2 on a grid in [0.0, 1.0] using 11 nodes,
stencils with 1 nodes to the left, 1 nodes to the right, and coefficients of Fornberg (1998)
  Calculation of Weights in Finite Difference Formulas.
  SIAM Rev. 40.3, pp. 685-691.
```
"""
function periodic_derivative_operator(derivative_order::Integer, accuracy_order,
                                      xmin, xmax, N, left_offset::Int=-(accuracy_order+1)÷2,
                                      mode=FastMode())
    if mode === Val(:serial) || mode === Val(:threads)
        # TODO: deprecated in v0.5
        Base.depwarn("Providing the argument `parallel` is deprecated." *
                     "Use `mode` instead.", :periodic_derivative_operator)
        mode = _parallel_to_mode(mode)
    end
    grid = range(xmin, stop=xmax, length=N+1) # N includes two identical boundary nodes
    coefficients = periodic_derivative_coefficients(derivative_order, accuracy_order, left_offset, eltype(grid), mode)
    PeriodicDerivativeOperator(coefficients, grid)
end

@inline function periodic_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N,
                                              mode, left_offset::Int=-(accuracy_order+1)÷2)
    if mode === Val(:serial) || mode === Val(:threads)
        # TODO: deprecated in v0.5
        Base.depwarn("Providing the argument `parallel` is deprecated." *
                     "Use `mode` instead.", :periodic_derivative_operator)
        mode = _parallel_to_mode(mode)
    end
    periodic_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N, left_offset, mode)
end

function periodic_derivative_operator(; derivative_order, accuracy_order,
                                      xmin, xmax, N, left_offset::Int=-(accuracy_order+1)÷2,
                                      mode=FastMode())
    if mode === Val(:serial) || mode === Val(:threads)
        # TODO: deprecated in v0.5
        Base.depwarn("Providing the argument `parallel` is deprecated." *
                     "Use `mode` instead.", :periodic_derivative_operator)
        mode = _parallel_to_mode(mode)
    end
    periodic_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N,
                                 left_offset, mode)
end

"""
    periodic_derivative_operator(source::Holoborodko2008,
                                 derivative_order, accuracy_order,
                                 xmin, xmax, N; mode=FastMode(), kwargs...)
    periodic_derivative_operator(source::Holoborodko2008;
                                 derivative_order, accuracy_order,
                                 xmin, xmax, N, mode=FastMode(), kwargs...)

Create a `PeriodicDerivativeOperator` approximating the `derivative_order`-th
derivative on a uniform grid between `xmin` and `xmax` with `N` grid points up
to order of accuracy `accuracy_order` where the leftmost grid point used is
determined by `left_offset`.
The evaluation of the derivative can be parallized using threads by chosing
`mode=ThreadedMode()`.

## Examples

```jldoctest
julia> periodic_derivative_operator(Holoborodko2008(), derivative_order=1, accuracy_order=2,
                                    xmin=0.0, xmax=1.0, N=11)
Periodic first-derivative operator of order 2 on a grid in [0.0, 1.0] using 11 nodes,
stencils with 2 nodes to the left, 2 nodes to the right, and coefficients of   Holoborodko (2008)
  Smooth Noise Robust Differentiators.
  http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
```
"""
function periodic_derivative_operator(source::Holoborodko2008, derivative_order, accuracy_order,
                                      xmin, xmax, N; mode=FastMode(), parallel=nothing, kwargs...)
    if parallel !== nothing
        # TODO: deprecated in v0.5
        Base.depwarn("Providing the keyword argument `parallel` is deprecated." *
                     "Use `mode` instead.", :periodic_derivative_operator)
        mode = _parallel_to_mode(parallel)
    end
    grid = range(xmin, stop=xmax, length=N+1) # N includes two identical boundary nodes
    coefficients = periodic_derivative_coefficients(source, derivative_order, accuracy_order;
                                                    kwargs..., T=eltype(grid), mode=mode)
    PeriodicDerivativeOperator(coefficients, grid)
end

function periodic_derivative_operator(source::Holoborodko2008; derivative_order, accuracy_order,
                                      xmin, xmax, N, mode=FastMode(), parallel=nothing, kwargs...)
    if parallel !== nothing
        # TODO: deprecated in v0.5
        Base.depwarn("Providing the keyword argument `parallel` is deprecated." *
                     "Use `mode` instead.", :periodic_derivative_operator)
        mode = _parallel_to_mode(parallel)
    end
    periodic_derivative_operator(source, derivative_order, accuracy_order, xmin, xmax, N;
                                 mode=mode, kwargs...)
end

"""
    periodic_derivative_operator(derivative_order, accuracy_order, grid,
                                 left_offset=-(accuracy_order+1)÷2, mode=FastMode())

Create a `PeriodicDerivativeOperator` approximating the `derivative_order`-th
derivative on thr uniform `grid` up to order of accuracy `accuracy_order` where
the leftmost grid point used is determined by `left_offset`.
The evaluation of the derivative can be parallized using threads by chosing
`mode=ThreadedMode())`.
"""
function periodic_derivative_operator(derivative_order, accuracy_order,
                                      grid::Union{LinRange,StepRangeLen},
                                      left_offset::Int=-(accuracy_order+1)÷2,
                                      mode=FastMode())
    if mode === Val(:serial) || mode === Val(:threads)
        # TODO: deprecated in v0.5
        Base.depwarn("Providing the argument `parallel` is deprecated." *
                     "Use `mode` instead.", :periodic_derivative_operator)
        mode = _parallel_to_mode(mode)
    end
    coefficients = periodic_derivative_coefficients(derivative_order, accuracy_order,
        left_offset, eltype(grid), mode)
    PeriodicDerivativeOperator(coefficients, grid)
end

@inline function periodic_derivative_operator(derivative_order, accuracy_order,
                                              grid::Union{LinRange,StepRangeLen},
                                              mode::AbstractExecutionMode,
                                              left_offset::Int=-(accuracy_order+1)÷2)
    if mode === Val(:serial) || mode === Val(:threads)
        # TODO: deprecated in v0.5
        Base.depwarn("Providing the argument `parallel` is deprecated." *
                     "Use `mode` instead.", :periodic_derivative_operator)
        mode = _parallel_to_mode(mode)
    end
    periodic_derivative_operator(derivative_order, accuracy_order, grid, left_offset, mode)
end

"""
    PeriodicDissipationOperator

A dissipation operator on a periodic finite difference grid.
See [`dissipation_operator`](@ref).
"""
struct PeriodicDissipationOperator{T,LowerOffset,UpperOffset,ExecutionMode,SourceOfCoefficients,Grid} <: AbstractPeriodicDerivativeOperator{T}
    factor::T
    Di::PeriodicDerivativeOperator{T,LowerOffset,UpperOffset,ExecutionMode,SourceOfCoefficients,Grid}

    function PeriodicDissipationOperator(factor::T, Di::PeriodicDerivativeOperator{T,LowerOffset,UpperOffset,ExecutionMode,SourceOfCoefficients,Grid}) where {T,LowerOffset,UpperOffset,ExecutionMode,SourceOfCoefficients,Grid}
        # check validity
        @argcheck iseven(derivative_order(Di))
        @argcheck (-1)^(1+derivative_order(Di)÷2)*factor >= 0
        new{T,LowerOffset,UpperOffset,ExecutionMode,SourceOfCoefficients,Grid}(factor, Di)
    end
end

grid(Di::PeriodicDissipationOperator) = grid(Di.Di)
derivative_order(Di::PeriodicDissipationOperator) = derivative_order(Di.Di)
accuracy_order(Di::PeriodicDissipationOperator) = accuracy_order(Di.Di)
source_of_coefficients(Di::PeriodicDissipationOperator) = MattssonSvärdNordström2004()

function Base.show(io::IO, Di::PeriodicDissipationOperator)
    if  derivative_order(Di) == 2
        print(io, "Periodic second-derivative dissipation operator")
    else
        print(io, "Periodic ", derivative_order(Di),
              "-derivative dissipation operator")
    end
    print(io, " of order ", accuracy_order(Di))
    if get(io, :compact, false) == false
        print(io, " on a grid in [", first(Di.Di.grid_evaluate), ", ", last(Di.Di.grid_evaluate),
                "] using ", size(Di, 1), " nodes \n")
        print(io, "and coefficients")
    end
    print(io, " of ", source_of_coefficients(Di))
end

# Compute `α*D*u + β*dest` and store the result in `dest`.
Base.@propagate_inbounds function mul!(dest::AbstractVector, Di::PeriodicDissipationOperator,
                                       u::AbstractVector, α, β)
    mul!(dest, Di.Di, u, Di.factor*α, β)
end

# Compute `α*D*u` and store the result in `dest`.
Base.@propagate_inbounds function mul!(dest::AbstractVector, Di::PeriodicDissipationOperator,
                                       u::AbstractVector, α)
    mul!(dest, Di.Di, u, Di.factor*α)
end

"""
    dissipation_operator(D::PeriodicDerivativeOperator;
                         strength=one(eltype(D)),
                         order=accuracy_order(D),
                         mode=D.coefficients.mode)

Create a negative semidefinite `DissipationOperator` using undivided differences
approximating a `order`-th derivative with strength `strength` adapted to the
derivative operator `D`.
The evaluation of the derivative can be parallized using threads by chosing
`mode=ThreadedMode()`.
"""
function dissipation_operator(D::PeriodicDerivativeOperator;
                              strength=eltype(D)(1),
                              order=accuracy_order(D),
                              mode=D.coefficients.mode,
                              parallel=nothing)
    if parallel !== nothing
        # TODO: deprecated in v0.5
        Base.depwarn("Providing the keyword argument `parallel` is deprecated." *
                     "Use `mode` instead.", :dissipation_operator)
        mode = _parallel_to_mode(parallel)
    end
    @argcheck iseven(order) ArgumentError("Dissipation operators require even derivatives.")
    # account for the grid spacing
    factor = strength * (-1)^(1 + order÷2) * D.Δx^order
    x = D.grid_evaluate
    Di = periodic_derivative_operator(order, order, first(x), last(x), size(D, 1), mode)
    PeriodicDissipationOperator(factor, Di)
end

struct PeriodicRationalDerivativeOperator{T<:Real, Dtype<:PeriodicDerivativeOperator{T}, Nnum, Nden, RFFT, IRFFT} <: AbstractPeriodicDerivativeOperator{T}
    D::Dtype
    num_coef::NTuple{Nnum,T}
    den_coef::NTuple{Nden,T}
    tmp::Vector{Complex{T}}
    eigval::Vector{Complex{T}}
    rfft_plan::RFFT
    irfft_plan::IRFFT

    function PeriodicRationalDerivativeOperator(D::Dtype, num_coef::NTuple{Nnum,T}, den_coef::NTuple{Nden,T}, tmp::Vector{Complex{T}}, eigval::Vector{Complex{T}}, rfft_plan::RFFT, irfft_plan::IRFFT) where {T<:Real, Dtype<:PeriodicDerivativeOperator{T}, Nnum, Nden, RFFT, IRFFT}
        @argcheck length(irfft_plan) == length(tmp) DimensionMismatch
        @argcheck length(irfft_plan) == length(eigval) DimensionMismatch
        @argcheck length(irfft_plan) == (length(rfft_plan)÷2)+1 DimensionMismatch
        @argcheck length(grid(D)) == length(rfft_plan) DimensionMismatch

        new{T,Dtype,Nnum,Nden,RFFT,IRFFT}(D, num_coef, den_coef, tmp, eigval, rfft_plan, irfft_plan)
    end
end

function _eigvals!(eigval::Vector{Complex{T}}, D::PeriodicDerivativeOperator{T}) where {T<:Real}
    @unpack factor = D
    @unpack lower_coef, central_coef, upper_coef = D.coefficients
    N = length(grid(D))

    @inbounds for idx in 1:length(eigval)
        tmp_eigval = zero(T)
        for j in Base.OneTo(length(lower_coef))
            tmp_eigval += factor * lower_coef[j] * (exp(2π * im * (idx-1) * (-j) / N))
        end
        tmp_eigval += factor * central_coef
        for j in Base.OneTo(length(upper_coef))
            tmp_eigval += factor * upper_coef[j] * (exp(2π * im * (idx-1) * (+j) / N))
        end
        eigval[idx] = tmp_eigval
    end

    eigval
end

function PeriodicRationalDerivativeOperator(D::PeriodicDerivativeOperator{T}) where {T<:Real}
    x = grid(D)
    u = zero.(x)
    rfft_plan = plan_rfft(u)
    tmp = rfft_plan * u
    irfft_plan = plan_irfft(tmp, length(x))

    eigval = similar(tmp)
    _eigvals!(eigval, D)

    PeriodicRationalDerivativeOperator(D, (zero(T), one(T)), (one(T),), tmp, eigval, rfft_plan, irfft_plan)
end

function PeriodicRationalDerivativeOperator(rat::PeriodicRationalDerivativeOperator, num_coef, den_coef)

    PeriodicRationalDerivativeOperator(rat.D, num_coef, den_coef, rat.tmp, rat.eigval, rat.rfft_plan, rat.irfft_plan)
end

function PeriodicRationalDerivativeOperator(D::PeriodicDerivativeOperator, num_coef, den_coef)

    PeriodicRationalDerivativeOperator(PeriodicRationalDerivativeOperator(D), num_coef, den_coef)
end

Base.size(rat::PeriodicRationalDerivativeOperator) = size(rat.D)
grid(rat::PeriodicRationalDerivativeOperator) = grid(rat.D)

function Base.show(io::IO, rat::PeriodicRationalDerivativeOperator)
    if get(io, :compact, false)
        print(io, "Rational periodic operator")
    else
        print(io, "Rational periodic 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.D)
    end
end

function mul_poly(coef1, coef2)
    T = promote_type(eltype(coef1), eltype(coef2))

    coef = ntuple(idx->zero(T), (length(coef1)-1) + (length(coef2)-1) + 1)
    for idx1 in 0:length(coef1)-1, idx2 in 0:length(coef2)-1
        coef = Base.setindex(coef, coef[idx1+idx2+1] + coef1[idx1+1] * coef2[idx2+1], idx1+idx2+1)
    end

    coef
end

function add_poly(coef1, coef2)
    T = promote_type(eltype(coef1), eltype(coef2))

    coef = ntuple(idx->zero(T), max(length(coef1), length(coef2)))
    for idx in 1:length(coef1)
        coef = Base.setindex(coef, coef[idx] + coef1[idx], idx)
    end
    for idx in 1:length(coef2)
        coef = Base.setindex(coef, coef[idx] + coef2[idx], idx)
    end

    coef
end

function subtract_poly(coef1, coef2)
    T = promote_type(eltype(coef1), eltype(coef2))

    coef = ntuple(idx->zero(T), max(length(coef1), length(coef2)))
    for idx in 1:length(coef1)
        coef = Base.setindex(coef, coef[idx] + coef1[idx], idx)
    end
    for idx in 1:length(coef2)
        coef = Base.setindex(coef, coef[idx] - coef2[idx], idx)
    end

    coef
end

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

function Base.:*(D1::PeriodicDerivativeOperator, D2::PeriodicDerivativeOperator)
    T = eltype(D1)
    @argcheck D1 == D2 ArgumentError

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

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

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

    PeriodicRationalDerivativeOperator(rat, num_coef, rat.den_coef)
end

function Base.:*(rat::PeriodicRationalDerivativeOperator, factor::Union{Real,Integer})
    factor * rat
end

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

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

function Base.:*(rat::PeriodicRationalDerivativeOperator, scaling::UniformScaling)
    scaling.λ * rat
end

function Base.:*(scaling::UniformScaling, rat::PeriodicRationalDerivativeOperator)
    rat * scaling
end

function Base.:*(D::PeriodicDerivativeOperator, scaling::UniformScaling)
    scaling * PeriodicRationalDerivativeOperator(D)
end

function Base.:*(scaling::UniformScaling, D::PeriodicDerivativeOperator)
    PeriodicRationalDerivativeOperator(D) * scaling
end

function Base.:+(rat1::PeriodicRationalDerivativeOperator, rat2::PeriodicRationalDerivativeOperator)
    @argcheck rat1.D == rat2.D ArgumentError

    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)
    PeriodicRationalDerivativeOperator(rat1, num_coef, den_coef)
end

function Base.:+(rat1::PeriodicRationalDerivativeOperator, rat2::PeriodicDerivativeOperator)
    rat1 + PeriodicRationalDerivativeOperator(rat2)
end

function Base.:+(rat1::PeriodicDerivativeOperator, rat2::PeriodicRationalDerivativeOperator)
    PeriodicRationalDerivativeOperator(rat1) + rat2
end

function Base.:+(rat::PeriodicRationalDerivativeOperator, scaling::UniformScaling)
    @unpack num_coef, den_coef = rat
    num_coef = add_poly(num_coef, den_coef)

    PeriodicRationalDerivativeOperator(rat, num_coef, den_coef)
end

function Base.:+(scaling::UniformScaling, rat::PeriodicRationalDerivativeOperator)
    rat + scaling
end

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

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

function Base.:-(rat1::PeriodicRationalDerivativeOperator, rat2::PeriodicRationalDerivativeOperator)
    @argcheck rat1.D == rat2.D ArgumentError

    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)
    PeriodicRationalDerivativeOperator(rat1, num_coef, den_coef)
end

function Base.:-(rat1::PeriodicRationalDerivativeOperator, rat2::PeriodicDerivativeOperator)
    rat1 - PeriodicRationalDerivativeOperator(rat2)
end

function Base.:-(rat1::PeriodicDerivativeOperator, rat2::PeriodicRationalDerivativeOperator)
    PeriodicRationalDerivativeOperator(rat1) - rat2
end

function Base.:-(rat::PeriodicRationalDerivativeOperator, scaling::UniformScaling)
    @unpack num_coef, den_coef = rat
    num_coef = subtract_poly(num_coef, den_coef)

    PeriodicRationalDerivativeOperator(rat, num_coef, den_coef)
end

function Base.:-(scaling::UniformScaling, rat::PeriodicRationalDerivativeOperator)
    @unpack num_coef, den_coef = rat
    num_coef = subtract_poly(den_coef, num_coef)

    PeriodicRationalDerivativeOperator(rat, num_coef, den_coef)
end

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

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

function Base.:*(rat1::PeriodicRationalDerivativeOperator, rat2::PeriodicRationalDerivativeOperator)
    @argcheck rat1.D == rat2.D ArgumentError

    num_coef = mul_poly(rat1.num_coef, rat2.num_coef)
    den_coef = mul_poly(rat1.den_coef, rat2.den_coef)
    PeriodicRationalDerivativeOperator(rat1, num_coef, den_coef)
end

function Base.:*(rat1::PeriodicRationalDerivativeOperator, rat2::PeriodicDerivativeOperator)
    rat1 * PeriodicRationalDerivativeOperator(rat2)
end

function Base.:*(rat1::PeriodicDerivativeOperator, rat2::PeriodicRationalDerivativeOperator)
    PeriodicRationalDerivativeOperator(rat1) * rat2
end

function Base.:/(rat1::PeriodicRationalDerivativeOperator, rat2::PeriodicRationalDerivativeOperator)
    @argcheck rat1.D == rat2.D ArgumentError

    num_coef = mul_poly(rat1.num_coef, rat2.den_coef)
    den_coef = mul_poly(rat1.den_coef, rat2.num_coef)
    PeriodicRationalDerivativeOperator(rat1, num_coef, den_coef)
end

function Base.:/(rat1::PeriodicRationalDerivativeOperator, rat2::PeriodicDerivativeOperator)
    rat1 / PeriodicRationalDerivativeOperator(rat2)
end

function Base.:/(rat1::PeriodicDerivativeOperator, rat2::PeriodicRationalDerivativeOperator)
    PeriodicRationalDerivativeOperator(rat1) / rat2
end

function mul!(dest::AbstractVector{T}, rat::PeriodicRationalDerivativeOperator, u::AbstractVector{T}) where {T}
    @unpack D, num_coef, den_coef, tmp, eigval, rfft_plan, irfft_plan = rat
    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))
        tmp[j] *= evalpoly(eigval[j], num_coef) / evalpoly(eigval[j], den_coef)
    end
    mul!(dest, irfft_plan, tmp)
end

function LinearAlgebra.ldiv!(dest::AbstractVector{T}, rat::PeriodicRationalDerivativeOperator, u::AbstractVector{T}) where {T}
    @unpack D, num_coef, den_coef, tmp, eigval, rfft_plan, irfft_plan = rat
    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))
        tmp[j] *= evalpoly(eigval[j], den_coef) / evalpoly(eigval[j], num_coef)
    end
    mul!(dest, irfft_plan, tmp)
end

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


text:
function maxpartsum(a)
    sum    = 0
    mins   = 0
    answer = 0
    for i in 1:lastindex(a)
        sum = sum+a[i]
        answer = max(answer, sum - mins )
        mins = min(sum, mins)
    end
    return answer
end


text:
"""
    edit_distance(G₁, G₂;
           insert_cost::Function=v->1.0,
           delete_cost::Function=u->1.0,
           subst_cost::Function=(u,v)->0.5,
           heuristic::Function=DefaultEditHeuristic)

Computes the edit distance between graphs G₁ and G₂.

Returns the minimum edit cost and edit path to transform graph
G₁ into graph G₂. An edit path consists of a sequence of pairs
of vertices (u,v) ∈ [0,|G₁|] × [0,|G₂|] representing vertex
operations:

- (0,v): insertion of vertex v ∈ G₂
- (u,0): deletion of vertex u ∈ G₁
- (u>0,v>0): substitution of vertex u ∈ G₁ by vertex v ∈ G₂

By default, the algorithm uses constant operation costs. The
user can provide classical Minkowski costs computed from vertex
labels μ₁ (for G₁) and μ₂ (for G₂) in order to further guide the
search, for example:
```
edit_distance(G₁, G₂, subst_cost=MinkowskiCost(μ₁, μ₂))
```
A custom heuristic can be provided to the A* search in case the
default heuristic is not satisfactory.
Performance tips:
-----------------
- Given two graphs ``|G₁| < |G₂|``, `edit_distance(G₁, G₂)` is faster to
compute than `edit_distance(G₂, G₁)`. Consider swapping the arguments
if involved costs are `symmetric`.
- The use of simple Minkowski costs can improve performance considerably.
- Exploit vertex attributes when designing operation costs.

For further details, please refer to:

RIESEN, K., 2015. Structural Pattern Recognition with Graph Edit
Distance: Approximation Algorithms and Applications. (Chapter 2)

Author: Júlio Hoffimann Mendes (juliohm@stanford.edu)
"""
function edit_distance(G₁::AGraphOrDiGraph, G₂::AGraphOrDiGraph;
                       insert_cost::Function=v->1.0,
                       delete_cost::Function=u->1.0,
                       subst_cost::Function=(u,v)->0.5,
                       heuristic::Function=DefaultEditHeuristic)

  # A* search heuristic
  h(λ) = heuristic(λ, G₁, G₂)

  # initialize open set
  OPEN = PriorityQueue{Vector{Tuple}, Float64}()
  for v in 1:nv(G₂)
    enqueue!(OPEN, [(1,v)], subst_cost(1,v) + h([(1,v)]))
  end
  enqueue!(OPEN, [(1,0)], delete_cost(1) + h([(1,0)]))

  while true
    # minimum (partial) edit path
    λ, cost = peek(OPEN)
    dequeue!(OPEN)

    if is_complete_path(λ, G₁, G₂)
      return cost, λ
    else
      k, _ = λ[end]
      vs = setdiff(1:nv(G₂), [v for (u,v) in λ])

      if k < nv(G₁) # there are still vertices to process in G₁?
        for v in vs
          λ⁺ = [λ; (k+1,v)]
          enqueue!(OPEN, λ⁺, cost + subst_cost(k+1,v) + h(λ⁺) - h(λ))
        end
        λ⁺ = [λ; (k+1,0)]
        enqueue!(OPEN, λ⁺, cost + delete_cost(k+1) + h(λ⁺) - h(λ))
      else
        # add remaining vertices of G₂ to the path
        λ⁺ = [λ; [(0,v) for v in vs]]
        total_insert_cost = sum(insert_cost, vs)
        enqueue!(OPEN, λ⁺, cost + total_insert_cost + h(λ⁺) - h(λ))
      end
    end
  end
end

function is_complete_path(λ, G₁, G₂)
  us = Set(); vs = Set()
  for (u,v) in λ
    push!(us, u)
    push!(vs, v)
  end
  delete!(us, 0)
  delete!(vs, 0)

  return length(us) == nv(G₁) && length(vs) == nv(G₂)
end

function DefaultEditHeuristic(λ, G₁::AGraphOrDiGraph, G₂::AGraphOrDiGraph)
  vs = Set()
  for (u,v) in λ
    push!(vs, v)
  end
  delete!(vs, 0)

  nv(G₂) - length(vs)
end

#-------------------------
# Edit path cost functions
#-------------------------

"""
For labels μ₁ on the vertices of graph G₁ and labels μ₂ on the vertices
of graph G₂, compute the p-norm cost of substituting vertex u ∈ G₁ by
vertex v ∈ G₂.
"""
function MinkowskiCost(μ₁::AbstractVector, μ₂::AbstractVector; p::Real=1)
  (u,v) -> norm(μ₁[u] - μ₂[v], p)
end

"""
Similar to MinkowskiCost, but ensures costs smaller than 2τ.
"""
function BoundedMinkowskiCost(μ₁::AbstractVector, μ₂::AbstractVector; p::Real=1, τ::Real=1)
  (u,v) -> 1 / (1/(2τ) + exp(-norm(μ₁[u] - μ₂[v], p)))
end