text:
#=
    Functions in this file test functionality relating to variables in MOI.

### Functionality currently tested
    - add_variables
    - add_variable
    - deleting variables
    - get/set VariableName
    - is_valid for VariableIndex
    - get VariableIndex by name
    - NumberOfVariables

### Functionality not yet tested
    - VariablePrimalStart
    - VariablePrimal
    - ListOfVariableIndices
=#

"""
    add_variable(model::MOI.ModelLike, config::Config)

Test adding a single variable.
"""
function add_variable(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    v = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1
end
unittests["add_variable"] = add_variable

"""
    add_variables(model::MOI.ModelLike, config::Config)

Test adding multiple variables.
"""
function add_variables(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    v = MOI.add_variables(model, 2)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2
end
unittests["add_variables"] = add_variables

"""
    delete_variable(model::MOI.ModelLike, config::Config)

Tess adding, and then deleting, a single variable.
"""
function delete_variable(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    v = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1
    MOI.delete(model, v)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
end
unittests["delete_variable"] = delete_variable

"""
    delete_variables(model::MOI.ModelLike, config::Config)

Test adding, and then deleting, multiple variables.
"""
function delete_variables(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    v = MOI.add_variables(model, 2)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2
    MOI.delete(model, v)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    v = MOI.add_variables(model, 2)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2
    MOI.delete(model, v[1])
    @test MOI.get(model, MOI.NumberOfVariables()) == 1
    @test_throws MOI.InvalidIndex{MOI.VariableIndex} MOI.delete(model, v[1])
    try
        MOI.delete(model, v[1])
    catch err
        @test err.index == v[1]
    end
    @test !MOI.is_valid(model, v[1])
    @test MOI.is_valid(model, v[2])
end
unittests["delete_variables"] = delete_variable

"""
    delete_nonnegative_variables(model::MOI.ModelLike, config::Config)

Test adding, and then deleting, nonnegative variables.
"""
function delete_nonnegative_variables(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    v, cv = MOI.add_constrained_variables(model, MOI.Nonnegatives(2))
    @test MOI.get(model, MOI.NumberOfVariables()) == 2
    MOI.delete(model, v)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    @test !MOI.is_valid(model, v[1])
    @test_throws MOI.InvalidIndex(v[1]) MOI.delete(model, v[1])
    @test !MOI.is_valid(model, v[2])
    @test_throws MOI.InvalidIndex(v[2]) MOI.delete(model, v[2])
    @test !MOI.is_valid(model, cv)
    v, cv = MOI.add_constrained_variables(model, MOI.Nonnegatives(1))
    @test MOI.get(model, MOI.NumberOfVariables()) == 1
    MOI.delete(model, v[1])
    @test !MOI.is_valid(model, v[1])
    @test_throws MOI.InvalidIndex(v[1]) MOI.delete(model, v[1])
    @test !MOI.is_valid(model, cv)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
end
unittests["delete_nonnegative_variables"] = delete_nonnegative_variables

"""
    update_dimension_nonnegative_variables(model::MOI.ModelLike, config::Config)

Test adding, and then deleting one by one, nonnegative variables.
"""
function update_dimension_nonnegative_variables(
    model::MOI.ModelLike,
    config::Config,
)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    v, cv = MOI.add_constrained_variables(model, MOI.Nonnegatives(2))
    @test MOI.get(model, MOI.NumberOfVariables()) == 2
    MOI.delete(model, v[1])
    @test !MOI.is_valid(model, v[1])
    @test_throws MOI.InvalidIndex(v[1]) MOI.delete(model, v[1])
    @test MOI.is_valid(model, cv)
    @test MOI.is_valid(model, v[2])
    @test MOI.get(model, MOI.ConstraintFunction(), cv) ==
          MOI.VectorOfVariables([v[2]])
    @test MOI.get(model, MOI.ConstraintSet(), cv) == MOI.Nonnegatives(1)
    MOI.delete(model, v[2])
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    @test !MOI.is_valid(model, v[1])
    @test_throws MOI.InvalidIndex(v[1]) MOI.delete(model, v[1])
    @test !MOI.is_valid(model, v[2])
    @test_throws MOI.InvalidIndex(v[2]) MOI.delete(model, v[2])
    @test !MOI.is_valid(model, cv)
end
unittests["update_dimension_nonnegative_variables"] =
    update_dimension_nonnegative_variables

"""
    delete_soc_variables(model::MOI.ModelLike, config::Config)

Test adding, and then deleting, second-order cone variables.
"""
function delete_soc_variables(model::MOI.ModelLike, config::Config)
    MOI.supports_add_constrained_variables(model, MOI.SecondOrderCone) || return
    MOI.empty!(model)
    @test MOI.is_empty(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    v, cv = MOI.add_constrained_variables(model, MOI.SecondOrderCone(3))
    @test MOI.get(model, MOI.NumberOfVariables()) == 3
    MOI.delete(model, v)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    @test !MOI.is_valid(model, v[1])
    @test_throws MOI.InvalidIndex(v[1]) MOI.delete(model, v[1])
    @test !MOI.is_valid(model, v[2])
    @test_throws MOI.InvalidIndex(v[2]) MOI.delete(model, v[2])
    @test !MOI.is_valid(model, cv)
    v, cv = MOI.add_constrained_variables(model, MOI.SecondOrderCone(3))
    @test MOI.get(model, MOI.NumberOfVariables()) == 3
    @test_throws MOI.DeleteNotAllowed MOI.delete(model, v[1])
end
unittests["delete_soc_variables"] = delete_soc_variables

"""
    getvariable(model::MOI.ModelLike, config::Config)

Test getting variables by name.
"""
function getvariable(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    MOIU.loadfromstring!(
        model,
        """
    variables: x
    minobjective: 2.0x
    x >= 1.0
    x <= 2.0
""",
    )
    @test MOI.get(model, MOI.VariableIndex, "y") === nothing
    x = MOI.get(model, MOI.VariableIndex, "x")
    @test MOI.is_valid(model, x)
end
unittests["getvariable"] = getvariable

"""
    variablenames(model::MOI.ModelLike, config::Config)

Test getting and setting variable names.
"""
function variablenames(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    v = MOI.add_variable(model)
    @test MOI.get(model, MOI.VariableName(), v) == ""
    @test MOI.supports(model, MOI.VariableName(), typeof(v))
    MOI.set(model, MOI.VariableName(), v, "x")
    @test MOI.get(model, MOI.VariableName(), v) == "x"
    MOI.set(model, MOI.VariableName(), v, "y")
    @test MOI.get(model, MOI.VariableName(), v) == "y"
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    @test MOI.get(model, MOI.VariableName(), x) == "x"
end
unittests["variablenames"] = variablenames

"""
    solve_with_upperbound(model::MOI.ModelLike, config::Config)

Test setting the upper bound of a variable, confirm that it solves correctly,
and if `config.duals=true`, check that the dual is computed correctly.
"""
function solve_with_upperbound(model::MOI.ModelLike, config::Config)
    atol, rtol = config.atol, config.rtol
    MOI.empty!(model)
    @test MOI.is_empty(model)
    MOIU.loadfromstring!(
        model,
        """
    variables: x
    maxobjective: 2.0x
    x <= 1.0
    x >= 0.0
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    c1 = MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}}(x.value)
    # We test this after the creation of every `VariableIndex` constraint
    # to ensure a good coverage of corner cases.
    @test c1.value == x.value
    c2 =
        MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}}(x.value)
    @test c2.value == x.value
    return test_model_solution(
        model,
        config;
        objective_value = 2.0,
        variable_primal = [(x, 1.0)],
        constraint_primal = [(c1, 1.0), (c2, 1.0)],
        constraint_dual = [(c1, -2.0), (c2, 0.0)],
    )
end
unittests["solve_with_upperbound"] = solve_with_upperbound

"""
    solve_with_lowerbound(model::MOI.ModelLike, config::Config)

Test setting the lower bound of a variable, confirm that it solves correctly,
and if `config.duals=true`, check that the dual is computed correctly.
"""
function solve_with_lowerbound(model::MOI.ModelLike, config::Config)
    atol, rtol = config.atol, config.rtol
    MOI.empty!(model)
    @test MOI.is_empty(model)
    MOIU.loadfromstring!(
        model,
        """
    variables: x
    minobjective: 2.0x
    x >= 1.0
    x <= 2.0
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    c1 =
        MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}}(x.value)
    @test c1.value == x.value
    c2 = MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}}(x.value)
    @test c2.value == x.value
    return test_model_solution(
        model,
        config;
        objective_value = 2.0,
        variable_primal = [(x, 1.0)],
        constraint_primal = [(c1, 1.0), (c2, 1.0)],
        constraint_dual = [(c1, 2.0), (c2, 0.0)],
    )
end
unittests["solve_with_lowerbound"] = solve_with_lowerbound

"""
    solve_integer_edge_cases(model::MOI.ModelLike, config::Config)

Test a variety of edge cases related to binary and integer variables.
"""
function solve_integer_edge_cases(model::MOI.ModelLike, config::Config)
    @testset "integer with lower bound" begin
        MOI.empty!(model)
        @test MOI.is_empty(model)
        MOIU.loadfromstring!(
            model,
            """
    variables: x
    minobjective: 2.0x
    x >= 1.5
    x in Integer()
""",
        )
        x = MOI.get(model, MOI.VariableIndex, "x")
        test_model_solution(
            model,
            config;
            objective_value = 4.0,
            variable_primal = [(x, 2.0)],
        )
    end
    @testset "integer with upper bound" begin
        MOI.empty!(model)
        @test MOI.is_empty(model)
        MOIU.loadfromstring!(
            model,
            """
    variables: x
    minobjective: -2.0x
    x <= 1.5
    x in Integer()
""",
        )
        x = MOI.get(model, MOI.VariableIndex, "x")
        test_model_solution(
            model,
            config;
            objective_value = -2.0,
            variable_primal = [(x, 1.0)],
        )
    end
    @testset "binary with upper" begin
        MOI.empty!(model)
        @test MOI.is_empty(model)
        MOIU.loadfromstring!(
            model,
            """
    variables: x
    minobjective: -2.0x
    x <= 2.0
    x in ZeroOne()
""",
        )
        x = MOI.get(model, MOI.VariableIndex, "x")
        test_model_solution(
            model,
            config;
            objective_value = -2.0,
            variable_primal = [(x, 1.0)],
        )
    end
    @testset "binary with 0 upper" begin
        MOI.empty!(model)
        @test MOI.is_empty(model)
        MOIU.loadfromstring!(
            model,
            """
    variables: x
    minobjective: 1.0x
    x <= 0.0
    x in ZeroOne()
""",
        )
        x = MOI.get(model, MOI.VariableIndex, "x")
        test_model_solution(
            model,
            config;
            objective_value = 0.0,
            variable_primal = [(x, 0.0)],
        )
    end
end
unittests["solve_integer_edge_cases"] = solve_integer_edge_cases


text:
using Rankers.Storages # store_with
using LOLTools.MatchV4: MatchV4, http_action
using LOLToolsDataDragon.DataDragon: en_US, ko_KR

api_key = get(ENV, "RIOT_TOKEN", "")
region = "kr"
matchId = 3999591966 # 3999564765
timelines = MatchV4.timelines(api_key, region, matchId; action=store_with(region, http_action, debug=true));
frame = timelines.frames[end]
frame_par = frame.participantFrames[9]
@info frame_par
@info frame.events
@info :timestamp frame.timestamp
@info :frameInterval timelines.frameInterval
using Dates
@info Time(Nanosecond(Millisecond(frame.timestamp - timelines.frameInterval)))

ext: jl


text:
# plots the contours of the OLS objective function,
# with the true and estimated coefficients
# generate data
using Plots, Econometrics
n = 20
x = [ones(n) randn(n)]
beta = randn(2)
y = x*beta+randn(n)
# ols estimate and objective function
betahat = x\y
s = (a,b)-> (1/n)*sum((y - x*[a,b]).^2)
# plot the contours and the points
closeall()
b1 = range(-4.0,stop=4.0,length=100)
b2 = range(-4.0,stop=4.0,length=100)
p = contour(b1,b2,(b1,b2)->s(b1,b2),fill=true, c=:viridis)
scatter!([beta[1]],[beta[2]], markersize=10, label="true")
scatter!([betahat[1]],[betahat[2]], markersize=10, label="estimated")
# do OLS by simulated annealing, and plot the path
lb = -4*ones(2,1)
ub = -lb
thetastart = fill(-3.0, 2, 1)
obj = theta-> (1/n)*sum((y - x*theta).^2)[1,1]
results = samin(obj, thetastart, lb, ub, verbosity=2)
path = results[4]
plot!(path[:,4], path[:,5])

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
using ADCME
using PyPlot
using ADCMEKit

n = 200
pml = 30
C = 1000.0
NT = 1000
Δt = 1.0/NT
x0 = LinRange(0, 1, n+1)
h = 1/n
xE = Array((0:n+2pml)*h .- pml*h)
xH = (xE[2:end]+xE[1:end-1])/2
N = n + 2pml + 1

σE = zeros(N)
for i = 1:pml
    d = i*h
    σE[pml + n + 1 + i] = C* (d/(pml*h))^3
    σE[pml+1-i] = C* (d/(pml*h))^3
end

σH = zeros(N-1)
for i = 1:pml
    d = (i-1/2)*h
    σH[pml + n + i] = C* (d/(pml*h))^3
    σH[pml+1-i] = C* (d/(pml*h))^3
end

function ricker(dt = 0.002, f0 = 5.0)
    nw = 2/(f0*dt)
    nc = floor(Int, nw/2)
    t = dt*collect(-nc:1:nc)
    b = (π*f0*t).^2
    w = @. (1 - 2b)*exp(-b)
end
R = ricker()
if length(R)<NT+1
    R = [R;zeros(NT+1-length(R))]
end
R = R[1:NT+1]
# tn = ( 0:NT ) *Δt
# R = @. exp( -20(tn-0.3)^2)
# error()
R_ = constant(R)

cH = ones(length(xH)) * 2.
cH[pml + Int(round(0.5÷h)): pml + Int(round(0.75÷h))] .= 1.0
cH[pml + Int(round(0.75÷h)): end] .= 1.

mask = ones(length(xH))
mask[pml+1:end-pml] .= 0.0
c_ = Variable(ones(N-2pml))

x0 = collect(((pml:N-pml -1) .- pml)*h)
# c_ = squeeze(fc(x0, [20,20,1])) + 1.0
# c_ = 2tanh(c_) + 2.5
cH = scatter_update(constant(cH.*mask), pml+1:N-pml, c_)
# cH = constant(cH)
cE = (cH[1:end-1]+cH[2:end])/2
Z = zeros(N)
Z[pml + pml÷2] = 1.0
Z = Z[2:end-1]

function condition(i, E_arr, H_arr)
    i<=NT+1
end

function body(i, E_arr, H_arr)
    E = read(E_arr, i-1)
    H = read(H_arr, i-1)
    ΔH = cH * (E[2:end]-E[1:end-1])/h - σH*H
    H += ΔH * Δt
    ΔE = cE * (H[2:end]-H[1:end-1])/h - σE[2:end-1]*E[2:end-1] + R_[i] * Z
    E = scatter_add(E, 2:N-1, ΔE * Δt)
    # E = scatter_update(E, N÷2, R_[i])
    i+1, write(E_arr, i, E), write(H_arr, i, H)
end

E_arr = TensorArray(NT+1)
H_arr = TensorArray(NT+1)

E_arr = write(E_arr, 1, zeros(N))
H_arr = write(H_arr, 1, zeros(N-1))

i = constant(2, dtype = Int32)

_, E, H = while_loop(condition, body, [i, E_arr, H_arr])

E = stack(E); E = set_shape(E, (NT+1, N))
H = stack(H)

loss = sum((E[:,pml+pml÷2] - E_[:, pml+pml÷2])^2)

sess = Session(); init(sess)
BFGS!(sess, loss)


text:
# associate names with matchers so that debugging is easier

set_name(x::Any, s::Symbol) = x
set_name(m::Matcher, s::Symbol) = (m.name = s; m)

set_names(x) = x
function set_names(node::Expr)
    if node.head == :(=) && length(node.args) == 2 && isa(node.args[1], Symbol)
        node.args[2] = Expr(:call, :set_name, node.args[2], QuoteNode(node.args[1]))
    end
    node.args = map(set_names, node.args)
    node
end

# this should be applied to a begin/end block, and sets the names of all
# matchers within that block to match the variables they are asigned to.

# so, for example
#   foo = Dot()
# will set the name field of Dot() to :foo.

# since this clutters the code with extra calls (to set_name()) it should only
# be used around code that is evaluated once, to construct the grammar.

macro with_names(block)
    esc(set_names(block))
end


text:
#Unary Operators
"""
Identity: `z=x`
"""
UnaryOps.IDENTITY
juliaop(::typeof(UnaryOps.IDENTITY)) = identity
UnaryOps.UnaryOp(::typeof(identity)) = UnaryOps.IDENTITY
"""
Additive Inverse: `z=-x`
"""
UnaryOps.AINV
juliaop(::typeof(UnaryOps.AINV)) = -
UnaryOps.UnaryOp(::typeof(-)) = UnaryOps.AINV
"""
Logical Negation

`z=¬x::Bool`

`Real`:  `z=¬(x::ℝ ≠ 0)`
"""
UnaryOps.LNOT
juliaop(::typeof(UnaryOps.LNOT)) = !
UnaryOps.UnaryOp(::typeof(!)) = UnaryOps.LNOT
"""
Multiplicative Inverse: `z=1/x`
"""
UnaryOps.MINV
juliaop(::typeof(UnaryOps.MINV)) = inv
UnaryOps.UnaryOp(::typeof(inv)) = UnaryOps.MINV
"""
One: `z=one(x)`
"""
UnaryOps.ONE
juliaop(::typeof(UnaryOps.ONE)) = one
UnaryOps.UnaryOp(::typeof(one)) = UnaryOps.ONE
"""
Absolute Value: `z=|x|`
"""
UnaryOps.ABS
juliaop(::typeof(UnaryOps.ABS)) = abs
UnaryOps.UnaryOp(::typeof(abs)) = UnaryOps.ABS
"""
Bitwise Negation: `z=¬x`
"""
UnaryOps.BNOT
juliaop(::typeof(UnaryOps.BNOT)) = ~
UnaryOps.UnaryOp(::typeof(~)) = UnaryOps.BNOT
"""
Square Root: `z=√(x)`
"""
UnaryOps.SQRT
juliaop(::typeof(UnaryOps.SQRT)) = sqrt
UnaryOps.UnaryOp(::typeof(sqrt)) = UnaryOps.SQRT
"""
Natural Logarithm: `z=logₑ(x)`
"""
UnaryOps.LOG
juliaop(::typeof(UnaryOps.LOG)) = log
UnaryOps.UnaryOp(::typeof(log)) = UnaryOps.LOG
"""
Natural Base Exponential: `z=eˣ`
"""
UnaryOps.EXP
juliaop(::typeof(UnaryOps.EXP)) = exp
UnaryOps.UnaryOp(::typeof(exp)) = UnaryOps.EXP
"""
Log Base 2: `z=log₂(x)`
"""
UnaryOps.LOG2
juliaop(::typeof(UnaryOps.LOG2)) = log2
UnaryOps.UnaryOp(::typeof(log2)) = UnaryOps.LOG2
"""
Sine: `z=sin(x)`
"""
UnaryOps.SIN
juliaop(::typeof(UnaryOps.SIN)) = sin
UnaryOps.UnaryOp(::typeof(sin)) = UnaryOps.SIN
"""
Cosine: `z=cos(x)`
"""
UnaryOps.COS
juliaop(::typeof(UnaryOps.COS)) = cos
UnaryOps.UnaryOp(::typeof(cos)) = UnaryOps.COS
"""
Tangent: `z=tan(x)`
"""
UnaryOps.TAN
juliaop(::typeof(UnaryOps.TAN)) = tan
UnaryOps.UnaryOp(::typeof(tan)) = UnaryOps.TAN
"""
Inverse Cosine: `z=cos⁻¹(x)`
"""
UnaryOps.ACOS
juliaop(::typeof(UnaryOps.ACOS)) = acos
UnaryOps.UnaryOp(::typeof(acos)) = UnaryOps.ACOS
"""
Inverse Sine: `z=sin⁻¹(x)`
"""
UnaryOps.ASIN
juliaop(::typeof(UnaryOps.ASIN)) = asin
UnaryOps.UnaryOp(::typeof(asin)) = UnaryOps.ASIN
"""
Inverse Tangent: `z=tan⁻¹(x)`
"""
UnaryOps.ATAN
juliaop(::typeof(UnaryOps.ATAN)) = atan
UnaryOps.UnaryOp(::typeof(atan)) = UnaryOps.ATAN
"""
Hyperbolic Sine: `z=sinh(x)`
"""
UnaryOps.SINH
juliaop(::typeof(UnaryOps.SINH)) = sinh
UnaryOps.UnaryOp(::typeof(sinh)) = UnaryOps.SINH
"""
Hyperbolic Cosine: `z=cosh(x)`
"""
UnaryOps.COSH
juliaop(::typeof(UnaryOps.COSH)) = cosh
UnaryOps.UnaryOp(::typeof(cosh)) = UnaryOps.COSH
"""
Hyperbolic Tangent: `z=tanh(x)`
"""
UnaryOps.TANH
juliaop(::typeof(UnaryOps.TANH)) = tanh
UnaryOps.UnaryOp(::typeof(tanh)) = UnaryOps.TANH
"""
Inverse Hyperbolic Sine: `z=sinh⁻¹(x)`
"""
UnaryOps.ASINH
juliaop(::typeof(UnaryOps.ASINH)) = asinh
UnaryOps.UnaryOp(::typeof(asinh)) = UnaryOps.ASINH
"""
Inverse Hyperbolic Cosine: `z=cosh⁻¹(x)`
"""
UnaryOps.ACOSH
juliaop(::typeof(UnaryOps.ACOSH)) = acosh
UnaryOps.UnaryOp(::typeof(acosh)) = UnaryOps.ACOSH
"""
Inverse Hyperbolic Tangent: `z=tanh⁻¹(x)`
"""
UnaryOps.ATANH
juliaop(::typeof(UnaryOps.ATANH)) = atanh
UnaryOps.UnaryOp(::typeof(atanh)) = UnaryOps.ATANH
"""
Sign Function: `z=signum(x)`
"""
UnaryOps.SIGNUM
juliaop(::typeof(UnaryOps.SIGNUM)) = sign
UnaryOps.UnaryOp(::typeof(sign)) = UnaryOps.SIGNUM
"""
Ceiling Function: `z=⌈x⌉`
"""
UnaryOps.CEIL
juliaop(::typeof(UnaryOps.CEIL)) = ceil
UnaryOps.UnaryOp(::typeof(ceil)) = UnaryOps.CEIL
"""
Floor Function: `z=⌊x⌋`
"""
UnaryOps.FLOOR
juliaop(::typeof(UnaryOps.FLOOR)) = floor
UnaryOps.UnaryOp(::typeof(floor)) = UnaryOps.FLOOR
"""
Round to nearest: `z=round(x)`
"""
UnaryOps.ROUND
juliaop(::typeof(UnaryOps.ROUND)) = round
UnaryOps.UnaryOp(::typeof(round)) = UnaryOps.ROUND
"""
Truncate: `z=trunc(x)`
"""
UnaryOps.TRUNC
juliaop(::typeof(UnaryOps.TRUNC)) = trunc
UnaryOps.UnaryOp(::typeof(trunc)) = UnaryOps.TRUNC
"""
Base-2 Exponential: `z=2ˣ`
"""
UnaryOps.EXP2
juliaop(::typeof(UnaryOps.EXP2)) = exp2
UnaryOps.UnaryOp(::typeof(exp2)) = UnaryOps.EXP2
"""
Natural Exponential - 1: `z=eˣ - 1`
"""
UnaryOps.EXPM1
juliaop(::typeof(UnaryOps.EXPM1)) = expm1
UnaryOps.UnaryOp(::typeof(expm1)) = UnaryOps.EXPM1
"""
Log Base 10: `z=log₁₀(x)`
"""
UnaryOps.LOG10
juliaop(::typeof(UnaryOps.LOG10)) = log10
UnaryOps.UnaryOp(::typeof(log10)) = UnaryOps.LOG10
"""
Natural Log of x + 1: `z=logₑ(x + 1)`
"""
UnaryOps.LOG1P
juliaop(::typeof(UnaryOps.LOG1P)) = log1p
UnaryOps.UnaryOp(::typeof(log1p)) = UnaryOps.LOG1P
"""
Log of Gamma Function: `z=log(|Γ(x)|)`
"""
UnaryOps.LGAMMA
juliaop(::typeof(UnaryOps.LGAMMA)) = lgamma
UnaryOps.UnaryOp(::typeof(lgamma)) = UnaryOps.LGAMMA
"""
Gamma Function: `z=Γ(x)`
"""
UnaryOps.TGAMMA
juliaop(::typeof(UnaryOps.TGAMMA)) = gamma
UnaryOps.UnaryOp(::typeof(gamma)) = UnaryOps.TGAMMA
"""
Error Function: `z=erf(x)`
"""
UnaryOps.ERF
juliaop(::typeof(UnaryOps.ERF)) = erf
UnaryOps.UnaryOp(::typeof(erf)) = UnaryOps.ERF
"""
Complimentary Error Function: `z=erfc(x)`
"""
UnaryOps.ERFC
juliaop(::typeof(UnaryOps.ERFC)) = erfc
UnaryOps.UnaryOp(::typeof(erfc)) = UnaryOps.ERFC

#There is no exact equivalent here, since Julia's frexp returns (frexpx, frexpe).
"""
Normalized Exponent: `z=frexpe(x)`
"""
UnaryOps.FREXPE
function frexpe end
juliaop(::typeof(UnaryOps.FREXPE)) = frexpe
UnaryOps.UnaryOp(::typeof(frexpe)) = UnaryOps.FREXPE
"""
Normalized Fraction: `z=frexpx(x)`
"""
UnaryOps.FREXPX
function frexpx end
juliaop(::typeof(UnaryOps.FREXPX)) = frexpx
UnaryOps.UnaryOp(::typeof(frexpx)) = UnaryOps.frexpx

"""
Complex Conjugate: `z=x̄`
"""
UnaryOps.CONJ
juliaop(::typeof(UnaryOps.CONJ)) = conj
UnaryOps.UnaryOp(::typeof(conj)) = UnaryOps.CONJ
"""
Real Part: `z=real(x)`
"""
UnaryOps.CREAL
juliaop(::typeof(UnaryOps.CREAL)) = real
UnaryOps.UnaryOp(::typeof(real)) = UnaryOps.CREAL
"""
Imaginary Part: `z=imag(x)`
"""
UnaryOps.CIMAG
juliaop(::typeof(UnaryOps.CIMAG)) = imag
UnaryOps.UnaryOp(::typeof(imag)) = UnaryOps.CIMAG
"""
Angle: `z=carg(x)`
"""
UnaryOps.CARG
juliaop(::typeof(UnaryOps.CARG)) = angle
UnaryOps.UnaryOp(::typeof(angle)) = UnaryOps.CARG
"""
isinf: `z=(x == ±∞)`
"""
UnaryOps.ISINF
juliaop(::typeof(UnaryOps.ISINF)) = isinf
UnaryOps.UnaryOp(::typeof(isinf)) = UnaryOps.ISINF
"""
isnan: `z=(x == NaN)`
"""
UnaryOps.ISNAN
juliaop(::typeof(UnaryOps.ISNAN)) = isnan
UnaryOps.UnaryOp(::typeof(isnan)) = UnaryOps.ISNAN
"""
isfinite: `z=isfinite(x)`
"""
UnaryOps.ISFINITE
juliaop(::typeof(UnaryOps.ISFINITE)) = isfinite
UnaryOps.UnaryOp(::typeof(isfinite)) = UnaryOps.ISFINITE
"""
0-based Row Index: `z=i`
"""
UnaryOps.POSITIONI #No Julia version since it's 0-based.
"""
1-Based Row Index: `z=i + 1`
"""
UnaryOps.POSITIONI1
function positioni end
juliaop(::typeof(UnaryOps.POSITIONI1)) = positioni
UnaryOps.UnaryOp(::typeof(positioni)) = UnaryOps.POSITIONI1
"""
0-Based Column Index: `z=j`
"""
UnaryOps.POSITIONJ #No Julia version since it's 0-based.
"""
1-Based Column Index: `z=j + 1`
"""
UnaryOps.POSITIONJ1
function positionj end
juliaop(::typeof(UnaryOps.POSITIONJ1)) = positionj
UnaryOps.UnaryOp(::typeof(positionj)) = UnaryOps.POSITIONJ1

#Binary Operators
"""
First argument: `f(x::T,y::T)::T = x`
"""
BinaryOps.FIRST
juliaop(::typeof(BinaryOps.FIRST)) = first
BinaryOps.BinaryOp(::typeof(first)) = BinaryOps.FIRST
"""
Second argument: `f(x::T,y::T)::T = y`
"""
BinaryOps.SECOND
function second end
juliaop(::typeof(BinaryOps.SECOND)) = second
BinaryOps.BinaryOp(::typeof(second)) = BinaryOps.SECOND
"""
Power: `f(x::T,y::T)::T = xʸ`
"""
BinaryOps.POW
juliaop(::typeof(BinaryOps.POW)) = ^
BinaryOps.BinaryOp(::typeof(^)) = BinaryOps.POW
"""
Addition: `f(x::T,y::T)::T = x + y`
"""
BinaryOps.PLUS
juliaop(::typeof(BinaryOps.PLUS)) = +
BinaryOps.BinaryOp(::typeof(+)) = BinaryOps.PLUS
"""
Subtraction: `f(x::T,y::T)::T = x - y`
"""
BinaryOps.MINUS
juliaop(::typeof(BinaryOps.MINUS)) = -
BinaryOps.BinaryOp(::typeof(-)) = BinaryOps.MINUS
"""
Multiplication: `f(x::T,y::T)::T = xy`
"""
BinaryOps.TIMES
juliaop(::typeof(BinaryOps.TIMES)) = *
BinaryOps.BinaryOp(::typeof(*)) = BinaryOps.TIMES
"""
Division: `f(x::T,y::T)::T = x / y`
"""
BinaryOps.DIV
juliaop(::typeof(BinaryOps.DIV)) = /
BinaryOps.BinaryOp(::typeof(/)) = BinaryOps.DIV
"""
Reverse Subtraction: `f(x::T,y::T)::T = y - x`
"""
BinaryOps.RMINUS
function rminus end
juliaop(::typeof(BinaryOps.RMINUS)) = rminus
BinaryOps.BinaryOp(::typeof(rminus)) = BinaryOps.RMINUS
"""
Reverse Division: `f(x::T,y::T)::T = y / x`
"""
BinaryOps.RDIV
juliaop(::typeof(BinaryOps.RDIV)) = \
BinaryOps.BinaryOp(::typeof(\)) = BinaryOps.RDIV
"""
One when both x and y exist: `f(x::T,y::T)::T = 1`
"""
BinaryOps.PAIR
function pair end
juliaop(::typeof(BinaryOps.PAIR)) = pair
BinaryOps.BinaryOp(::typeof(pair)) = BinaryOps.PAIR
"""
Pick x or y arbitrarily: `f(x::T,y::T)::T = x or y`
"""
BinaryOps.ANY
#This is sort of incorrect
juliaop(::typeof(BinaryOps.ANY)) = any
BinaryOps.BinaryOp(::typeof(any)) = BinaryOps.ANY
"""
Equal: `f(x::T,y::T)::T = x == y``
"""
BinaryOps.ISEQ
function iseq end
juliaop(::typeof(BinaryOps.ISEQ)) = iseq
BinaryOps.BinaryOp(::typeof(iseq)) = BinaryOps.ISEQ
"""
Not Equal: `f(x::T,y::T)::T = x ≠ y`
"""
BinaryOps.ISNE
function isne end
juliaop(::typeof(BinaryOps.ISNE)) = isne
BinaryOps.BinaryOp(::typeof(isne)) = BinaryOps.ISNE
"""
Greater Than: `f(x::ℝ,y::ℝ)::ℝ = x > y`
"""
BinaryOps.ISGT
function isgt end
juliaop(::typeof(BinaryOps.ISGT)) = isgt
BinaryOps.BinaryOp(::typeof(isgt)) = BinaryOps.ISGT
"""
Less Than: `f(x::ℝ,y::ℝ)::ℝ = x < y`
"""
BinaryOps.ISLT
function islt end
juliaop(::typeof(BinaryOps.ISLT)) = islt
BinaryOps.BinaryOp(::typeof(islt)) = BinaryOps.ISLT
"""
Greater Than or Equal: `f(x::ℝ,y::ℝ)::ℝ = x ≥ y`
"""
BinaryOps.ISGE
function isge end
juliaop(::typeof(BinaryOps.ISGE)) = isge
BinaryOps.BinaryOp(::typeof(isge)) = BinaryOps.ISGE
"""
Less Than or Equal: `f(x::ℝ,y::ℝ)::ℝ = x ≤ y`
"""
BinaryOps.ISLE
function isle end
juliaop(::typeof(BinaryOps.ISLE)) = isle
BinaryOps.BinaryOp(::typeof(isle)) = BinaryOps.ISLE
"""
Minimum: `f(x::ℝ,y::ℝ)::ℝ = min(x, y)`
"""
BinaryOps.MIN
juliaop(::typeof(BinaryOps.MIN)) = min
BinaryOps.BinaryOp(::typeof(min)) = BinaryOps.MIN
"""
Maximum: `f(x::ℝ,y::ℝ)::ℝ = max(x, y)`
"""
BinaryOps.MAX
juliaop(::typeof(BinaryOps.MAX)) = max
BinaryOps.BinaryOp(::typeof(max)) = BinaryOps.MAX
"""
Logical OR: `f(x::ℝ,y::ℝ)::ℝ = (x ≠ 0) ∨ (y ≠ 0)`
"""
BinaryOps.LOR
function ∨ end
juliaop(::typeof(BinaryOps.LOR)) = ∨
BinaryOps.BinaryOp(::typeof(∨)) = BinaryOps.LOR
"""
Logical AND: `f(x::ℝ,y::ℝ)::ℝ = (x ≠ 0) ∧ (y ≠ 0)`
"""
BinaryOps.LAND
function ∧ end
juliaop(::typeof(BinaryOps.LAND)) = ∧
BinaryOps.BinaryOp(::typeof(∧)) = BinaryOps.LAND
"""
Logical AND: `f(x::ℝ,y::ℝ)::ℝ = (x ≠ 0) ⊻ (y ≠ 0)`
"""
BinaryOps.LXOR
function lxor end
juliaop(::typeof(BinaryOps.LXOR)) = lxor
BinaryOps.BinaryOp(::typeof(lxor)) = BinaryOps.LXOR
"""
4-Quadrant Arc Tangent: `f(x::F, y::F)::F = tan⁻¹(y/x)`
"""
BinaryOps.ATAN2
juliaop(::typeof(BinaryOps.ATAN2)) = atan
BinaryOps.BinaryOp(::typeof(atan)) = BinaryOps.ATAN2
"""
Hypotenuse: `f(x::F, y::F)::F = √(x² + y²)`
"""
BinaryOps.HYPOT
juliaop(::typeof(BinaryOps.HYPOT)) = hypot
BinaryOps.BinaryOp(::typeof(hypot)) = BinaryOps.HYPOT
"""
Float remainder of x / y rounded towards zero.
"""
BinaryOps.FMOD
#Is this available?
function fmod end
juliaop(::typeof(BinaryOps.FMOD)) = fmod
BinaryOps.BinaryOp(::typeof(fmod)) = BinaryOps.FMOD
"""
Float remainder of x / y rounded towards nearest integral value.
"""
BinaryOps.REMAINDER
juliaop(::typeof(BinaryOps.REMAINDER)) = rem
BinaryOps.BinaryOp(::typeof(rem)) = BinaryOps.REMAINDER
"""
LDEXP: `f(x::F, y::F)::F = x × 2ⁿ`
"""
BinaryOps.LDEXP
juliaop(::typeof(BinaryOps.LDEXP)) = ldexp
BinaryOps.BinaryOp(::typeof(ldexp)) = BinaryOps.LDEXP
"""
Copysign: Value with magnitude of x and sign of y.
"""
BinaryOps.COPYSIGN
juliaop(::typeof(BinaryOps.COPYSIGN)) = copysign
BinaryOps.BinaryOp(::typeof(copysign)) = BinaryOps.COPYSIGN
"""
Bitwise OR: `f(x::ℤ, y::ℤ)::ℤ = x | y`
"""
BinaryOps.BOR
juliaop(::typeof(BinaryOps.BOR)) = |
BinaryOps.BinaryOp(::typeof(|)) = BinaryOps.BOR
"""
Bitwise AND: `f(x::ℤ, y::ℤ)::ℤ = x & y`
"""
BinaryOps.BAND
juliaop(::typeof(BinaryOps.BAND)) = &
BinaryOps.BinaryOp(::typeof(&)) = BinaryOps.BAND
"""
Bitwise XOR: `f(x::ℤ, y::ℤ)::ℤ = x ^ y`
"""
BinaryOps.BXOR
juliaop(::typeof(BinaryOps.BXOR)) = ⊻
BinaryOps.BinaryOp(::typeof(⊻)) = BinaryOps.BXOR
"""
Bitwise XNOR: : `f(x::ℤ, y::ℤ)::ℤ = ~(x ^ y)`
"""
BinaryOps.BXNOR
juliaop(::typeof(BinaryOps.BXNOR)) = !⊻
BinaryOps.BinaryOp(::typeof(!⊻)) = BinaryOps.BXNOR
"""
BGET: `f(x::ℤ, y::ℤ)::ℤ = get bit y of x.`
"""
BinaryOps.BGET
"""
BSET: `f(x::ℤ, y::ℤ)::ℤ = set bit y of x.`
"""
BinaryOps.BSET
"""
BCLR: `f(x::ℤ, y::ℤ)::ℤ = clear bit y of x.`
"""
BinaryOps.BCLR
"""
BSHIFT: `f(x::ℤ, y::Int8)::ℤ = bitshift(x, y)`
"""
BinaryOps.BSHIFT

"""
Equals: `f(x::T, y::T)::Bool = x == y`
"""
BinaryOps.EQ
juliaop(::typeof(BinaryOps.EQ)) = ==
BinaryOps.BinaryOp(::typeof(==)) = BinaryOps.EQ
"""
Not Equals: `f(x::T, y::T)::Bool = x ≠ y`
"""
BinaryOps.NE
juliaop(::typeof(BinaryOps.NE)) = !=
BinaryOps.BinaryOp(::typeof(!=)) = BinaryOps.NE
"""
Greater Than: `f(x::T, y::T)::Bool = x > y`
"""
BinaryOps.GT
juliaop(::typeof(BinaryOps.GT)) = >
BinaryOps.BinaryOp(::typeof(>)) = BinaryOps.GT
"""
Less Than: `f(x::T, y::T)::Bool = x < y`
"""
BinaryOps.LT
juliaop(::typeof(BinaryOps.LT)) = <
BinaryOps.BinaryOp(::typeof(<)) = BinaryOps.LT
"""
Greater Than or Equal: `f(x::T, y::T)::Bool = x ≥ y`
"""
BinaryOps.GE
juliaop(::typeof(BinaryOps.GE)) = >=
BinaryOps.BinaryOp(::typeof(>=)) = BinaryOps.GE
"""
Less Than or Equal: `f(x::T, y::T)::Bool = x ≤ y`
"""
BinaryOps.LE
juliaop(::typeof(BinaryOps.LE)) = <=
BinaryOps.BinaryOp(::typeof(<=)) = BinaryOps.LE
"""
Complex: `f(x::F, y::F)::Complex = x + y × i`
"""
BinaryOps.CMPLX
juliaop(::typeof(BinaryOps.CMPLX)) = complex
BinaryOps.BinaryOp(::typeof(complex)) = BinaryOps.CMPLX
"""
0-Based row index of a: `f(aᵢⱼ::T, bₖₗ::T)::Int64 = i`
"""
BinaryOps.FIRSTI
"""
1-Based row index of a: `f(aᵢⱼ::T, bₖₗ::T)::Int64 = i + 1`
"""
BinaryOps.FIRSTI1
function firsti end
juliaop(::typeof(BinaryOps.FIRSTI1)) = firsti
BinaryOps.BinaryOp(::typeof(firsti)) = BinaryOps.FIRSTI1
"""
0-Based column index of a: `f(aᵢⱼ::T, bₖₗ::T)::Int64 = j`
"""
BinaryOps.FIRSTJ
"""
1-Based column index of a: `f(aᵢⱼ::T, bₖₗ::T)::Int64 = j + 1`
"""
BinaryOps.FIRSTJ1
function firstj end
juliaop(::typeof(BinaryOps.FIRSTJ1)) = firstj
BinaryOps.BinaryOp(::typeof(firstj)) = BinaryOps.FIRSTJ1
"""
0-Based row index of b: `f(aᵢⱼ::T, bₖₗ::T)::Int64 = k`
"""
BinaryOps.SECONDI
"""
0-Based row index of b: `f(aᵢⱼ::T, bₖₗ::T)::Int64 = k + 1`
"""
BinaryOps.SECONDI1
function secondi end
juliaop(::typeof(BinaryOps.SECONDI1)) = secondi
BinaryOps.BinaryOp(::typeof(secondi)) = BinaryOps.SECONDI1
"""
0-Based column index of b: `f(aᵢⱼ::T, bₖₗ::T)::Int64 = l`
"""
BinaryOps.SECONDJ
"""
1-Based column index of b: `f(aᵢⱼ::T, bₖₗ::T)::Int64 = l + 1`
"""
BinaryOps.SECONDJ1
function secondj end
juliaop(::typeof(BinaryOps.SECONDJ1)) = secondj
BinaryOps.BinaryOp(::typeof(secondj)) = BinaryOps.SECONDJ1

#All binary ops will default to emul
defaultadd(f) = emul
# Default to eadd. We're limiting this to + and OR for now to enable easy graph unions.
for op ∈ [
    :+,
    :∨,
]
    funcquote = quote
        defaultadd(::typeof($op)) = eadd
    end
    @eval($funcquote)
end

#Monoid operators
"""
Minimum monoid: `f(x::ℝ, y::ℝ)::ℝ = min(x, y)`
* Identity: +∞
* Terminal: -∞
"""
Monoids.MIN_MONOID
op(::typeof(Monoids.MIN_MONOID)) = BinaryOps.MIN
Monoids.Monoid(::typeof(min)) = Monoids.MIN_MONOID
"""
Max monoid: `f(x::ℝ, y::ℝ)::ℝ = max(x, y)`
* Identity: -∞
* Terminal: +∞
"""
Monoids.MAX_MONOID
op(::typeof(Monoids.MAX_MONOID)) = BinaryOps.MAX
Monoids.Monoid(::typeof(max)) = Monoids.MAX_MONOID
"""
Plus monoid: `f(x::T, y::T)::T = x + y`
* Identity: 0
* Terminal: nothing
"""
Monoids.PLUS_MONOID
op(::typeof(Monoids.PLUS_MONOID)) = BinaryOps.PLUS
Monoids.Monoid(::typeof(+)) = Monoids.PLUS_MONOID
"""
Times monoid: `f(x::T, y::T)::T = xy`
* Identity: 1
* Terminal: 0 for non Floating-point numbers.
"""
Monoids.TIMES_MONOID
op(::typeof(Monoids.TIMES_MONOID)) = BinaryOps.TIMES
Monoids.Monoid(::typeof(*)) = Monoids.TIMES_MONOID
"""
Any monoid: `f(x::T, y::T)::T = x or y`
* Identity: any
* Terminal: any
"""
Monoids.ANY_MONOID
op(::typeof(Monoids.ANY_MONOID)) = BinaryOps.ANY
Monoids.Monoid(::typeof(any)) = Monoids.ANY_MONOID
"""
Logical OR monoid: `f(x::Bool, y::Bool)::Bool = x ∨ y`
* Identity: false
* Terminal: true
"""
Monoids.LOR_MONOID
op(::typeof(Monoids.LOR_MONOID)) = BinaryOps.LOR
Monoids.Monoid(::typeof(∨)) = Monoids.LOR_MONOID
"""
Logical AND monoid: `f(x::Bool, y::Bool)::Bool = x ∧ y`
* Identity: true
* Terminal: false
"""
Monoids.LAND_MONOID
op(::typeof(Monoids.LAND_MONOID)) = BinaryOps.LAND
Monoids.Monoid(::typeof(∧)) = Monoids.LAND_MONOID
"""
Logical XOR monoid: `f(x::Bool, y::Bool)::Bool = x ⊻ y`
* Identity: false
* Terminal: nothing
"""
Monoids.LXOR_MONOID
op(::typeof(Monoids.LXOR_MONOID)) = BinaryOps.LXOR
Monoids.Monoid(::typeof(lxor)) = Monoids.LXOR_MONOID
"""
Logical XNOR monoid: `f(x::Bool, y::Bool)::Bool = x == y`
* Identity: true
* Terminal: nothing
"""
Monoids.LXNOR_MONOID
#Don't care, this is ==.
"""
Boolean Equality `f(x::Bool, y::Bool)::Bool = x == y`.
"""
Monoids.EQ_MONOID
op(::typeof(Monoids.EQ_MONOID)) = BinaryOps.EQ
Monoids.Monoid(::typeof(==)) = Monoids.EQ_MONOID
"""
Bitwise OR monoid: `f(x::ℤ, y::ℤ)::ℤ = x|y`
* Identity: All bits `0`.* Terminal: All bits `1`.
"""
Monoids.BOR_MONOID
op(::typeof(Monoids.BOR_MONOID)) = BinaryOps.BOR
Monoids.Monoid(::typeof(|)) = Monoids.BOR_MONOID
"""
Bitwise AND monoid: `f(x::ℤ, y::ℤ)::ℤ = x&y`
* Identity: All bits `1`.
* Terminal: All bits `0`.
"""
Monoids.BAND_MONOID
op(::typeof(Monoids.BAND_MONOID)) = BinaryOps.BAND
Monoids.Monoid(::typeof(&)) = Monoids.BAND_MONOID
"""
Bitwise XOR monoid: `f(x::ℤ, y::ℤ)::ℤ = x^y`
* Identity: All bits `0`.
* Terminal: nothing
"""
Monoids.BXOR_MONOID
op(::typeof(Monoids.BXOR_MONOID)) = BinaryOps.BXOR
Monoids.Monoid(::typeof(⊻)) = Monoids.BXOR_MONOID
"""
Bitwise XNOR monoid: `f(x::ℤ, y::ℤ)::ℤ = ~(x^y)`
* Identity: All bits `1`.
* Terminal: nothing
"""
Monoids.BXNOR_MONOID
op(::typeof(Monoids.BXNOR_MONOID)) = BinaryOps.BXNOR
Monoids.Monoid(::typeof(!⊻)) = Monoids.BXNOR_MONOID

mulop(rig::AbstractSemiring) = nothing
addop(rig::AbstractSemiring) = nothing

# TODO:
# The definitions below are less than ideal, these should likely be fields of the semiring
# structs, not evaluated functions.
for oplus ∈ [(:max, "MAX"), (:min, "MIN"), (:+, "PLUS"), (:*, "TIMES"), (:any, "ANY")]
    for otimes ∈ [
        (:/, "DIV"),
        (:\, "RDIV"),
        (:first, "FIRST"),
        (:firsti, "FIRSTI1"),
        (:firstj, "FIRSTJ1"),
        (:iseq, "ISEQ"),
        (:isge, "ISGE"),
        (:isgt, "ISGT"),
        (:isle, "ISLE"),
        (:islt, "ISLT"),
        (:isne, "ISNE"),
        (:∧, "LAND"),
        (:∨, "LOR"),
        (:lxor, "LXOR"),
        (:max, "MAX"),
        (:min, "MIN"),
        (:-, "MINUS"),
        (:rminus, "RMINUS"),
        (:second, "SECOND"),
        (:secondi, "SECONDI1"),
        (:secondj, "SECONDJ1"),
        (:*, "TIMES"),
        (:+, "PLUS"),
        (:pair, "PAIR")
    ]
    rig = Symbol(oplus[2], "_", otimes[2])
    funcquote = quote
        Semirings.Semiring(::typeof($(oplus[1])), ::typeof($(otimes[1]))) = $rig
        addop(::typeof($rig)) = $(oplus[1])
        mulop(::typeof($rig)) = $(otimes[1])
    end
    @eval($funcquote)
    end
end

for otimes ∈ [(:(==), "EQ"), (:(>=), "GE"), (:(>), "GT"), (:(<=), "LE"), (:(<), "LT"), (:(!=), "NE")]
    oplus = (:any, "ANY")
    rig = Symbol(oplus[2], "_", otimes[2])
    funcquote = quote
        Semirings.Semiring(::typeof($(oplus[1])), ::typeof($(otimes[1]))) = $rig
        addop(::typeof($rig)) = $(oplus[1])
        mulop(::typeof($rig)) = $(otimes[1])
    end
    @eval($funcquote)
end
for oplus ∈ [(:(==), "EQ"), (:∧, "LAND"), (:∨, "LOR"), (:lxor, "LXOR")]
    for otimes ∈ [
        (:(==), "EQ"),
        (:(>=), "GE"),
        (:(>), "GT"),
        (:(<=), "LE"),
        (:(<), "LT"),
        (:(!=), "NE"),
        (:first, "FIRST"),
        (:second, "SECOND"),
        (:∧, "LAND"),
        (:∨, "LOR"),
        (:lxor, "LXOR"),
        (:pair, "PAIR")
    ]
    rig = Symbol(oplus[2], "_", otimes[2])
    funcquote = quote
        Semirings.Semiring(::typeof($(oplus[1])), ::typeof($(otimes[1]))) = $rig
        addop(::typeof($rig)) = $(oplus[1])
        mulop(::typeof($rig)) = $(otimes[1])
    end
    @eval($funcquote)
    end
end

#Select Ops

"""
    select(SuiteSparseGraphBLAS.TRIL, A, k=0)
    select(tril, A, k=0)

Select the entries on or below the `k`th diagonal of A.

See also: `LinearAlgebra.tril`
"""
TRIL
SelectOp(::typeof(LinearAlgebra.tril)) = TRIL
"""
    select(SuiteSparseGraphBLAS.TRIU, A, k=0)
    select(triu, A, k=0)

Select the entries on or above the `k`th diagonal of A.

See also: `LinearAlgebra.triu`
"""
TRIU
SelectOp(::typeof(LinearAlgebra.triu)) = TRIU
"""
    select(DIAG, A, k=0)

Select the entries on the `k`th diagonal of A.

See also: `LinearAlgebra.diag`
"""
DIAG
SelectOp(::typeof(LinearAlgebra.diag)) = DIAG
"""
    select(OFFDIAG, A, k=0)

Select the entries **not** on the `k`th diagonal of A.
"""
OFFDIAG
SelectOp(::typeof(offdiag)) = OFFDIAG
"""
    select(NONZERO, A)
    select(nonzeros, A)
Select all entries in A with nonzero value.
"""
NONZERO
SelectOp(::typeof(nonzeros)) = NONZERO

# I don't believe these should have Julia equivalents.
# Instead select(==, A, 0) will find EQ_ZERO internally.
"""
    select(EQ_ZERO, A)

Select all entries in A equal to zero.
"""
EQ_ZERO
"""
    select(GT_ZERO, A)

Select all entries in A greater than zero.
"""
GT_ZERO
"""
    select(GE_ZERO, A)

Select all entries in A greater than or equal to zero.
"""
GE_ZERO
"""
    select(LT_ZERO, A)

Select all entries in A less than zero.
"""
LT_ZERO
"""
    select(LE_ZERO, A)

Select all entries in A less than or equal to zero.
"""
LE_ZERO
"""
    select(NE, A, k)
    select(!=, A, k)
Select all entries not equal to `k`.
"""
NE
SelectOp(::typeof(!=)) = NE
"""
    select(EQ, A, k)
    select(==, A, k)
Select all entries equal to `k`.
"""
EQ
SelectOp(::typeof(==)) = EQ
"""
    select(GT, A, k)
    select(>, A, k)
Select all entries greater than `k`.
"""
GT
SelectOp(::typeof(>)) = GT
"""
    select(GE, A, k)
    select(>=, A, k)
Select all entries greater than or equal to `k`.
"""
GE
SelectOp(::typeof(>=)) = GE
"""
    select(LT, A, k)
    select(<, A, k)
Select all entries less than `k`.
"""
LT
SelectOp(::typeof(<)) = LT
"""
    select(LE, A, k)
    select(<=, A, k)
Select all entries less than or equal to `k`.
"""
LE
SelectOp(::typeof(<=)) = LE


text:
function ShiftingTransformProvider(arg0::TransformProvider, arg1::CartesianDerivativesFilter, arg2::AngularDerivativesFilter, arg3::jint, arg4::jdouble, arg5::jint, arg6::jdouble, arg7::jdouble)
    return ShiftingTransformProvider((TransformProvider, CartesianDerivativesFilter, AngularDerivativesFilter, jint, jdouble, jint, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function get_grid_points(obj::ShiftingTransformProvider)
    return jcall(obj, "getGridPoints", jint, ())
end

function get_raw_provider(obj::ShiftingTransformProvider)
    return jcall(obj, "getRawProvider", TransformProvider, ())
end

function get_step(obj::ShiftingTransformProvider)
    return jcall(obj, "getStep", jdouble, ())
end

function get_transform(obj::ShiftingTransformProvider, arg0::AbsoluteDate)
    return jcall(obj, "getTransform", Transform, (AbsoluteDate,), arg0)
end

function get_transform(obj::ShiftingTransformProvider, arg0::FieldAbsoluteDate)
    return jcall(obj, "getTransform", FieldTransform, (FieldAbsoluteDate,), arg0)
end


text:
using PermutedArrays
using Test
using Random

@testset "PermutedArrays.jl" begin

    function testiperm(v)
        n = length(v)
        @test v.perm[v.iperm] == 1:n
        @test v.iperm[v.perm] == 1:n
    end

    v =PermutedVector(collect(1:9));

    @testset "`permute!`" begin
        permute!(v, randperm(9));
        testiperm(v)
    end

    @testset "`swap!`" begin
        swap!(v, 2, 7);
        testiperm(v)
    end

    @testset "`deleteat!`" begin
        deleteat!(v, 3)
        testiperm(v)
    end

    @testset "`sortperm`" begin
        ℓ = sort(randn(100));
        x = PermutedVector(1:100);

        p = randperm(100);

        permute!(x,p);
        permute!(ℓ,p);

        @test sortperm(ℓ) == sortperm(x)
        testiperm(v)
    end

    @testset "`x .= x[p]`" begin
        ℓ = sort(randn(100));
        x = PermutedVector([1:100;]);

        p = randperm(100);

        x .= x[p]
        ℓ .= ℓ[p]

        @test sortperm(ℓ) == sortperm(x)
        testiperm(v)
    end

end

max_stars_repo_name: cscherrer/PermutedArrays

max_issues_repo_name: cscherrer/PermutedArrays

max_forks_repo_name: cscherrer/PermutedArrays


text:
using PhysicsCodesLabXD
using Test

@testset "PhysicsCodesLabXD.jl" begin
    # Write your tests here.
end


text:
# module TestODEModel
using CovidModels
using CovidModels.Data
using OrdinaryDiffEq
import OrdinaryDiffEq: ODEProblem
using Petri
using Test
using Plots
using StatsPlots
using Colors

function fitplot(sol, df)
    p = plot(sol, vars=[2, 4, 5, 6], label=["IA" "I" "R" "D"])
    @df df scatter!(p,
            0:length(italy.Deaths)-1,
            [:Cases, :Deaths, :Recoveries],
            label=["cases" "deaths" "recoveries"],
            marker=true,
            color=[RGB(1,0,0) RGB(0.5,0,0.5) RGB(0,1,0)]
            )
    return p
end

function makeplots_seir(sol, prefix)
    mkpath(dirname(prefix))
    p1 = plot(sol,vars=[1,2,3,4], xlabel="", ylabel="people", linewidth=3,title="Cities", legend=false)
    p2 = plot(sol,vars=[5,6,7,8], xlabel="", ylabel="people", linewidth=3, legend=false)
    p3 = plot(sol,vars=[9,10,11,12], xlabel="time", ylabel="people", linewidth=3, legend=false)
    p4 = plot(sol,vars=[2,6,10], xlabel="", linewidth=3, labels=["i1" "i2" "i3"], legend=true)
    p5 = plot(sol,vars=[3,7,11], xlabel="", linewidth=3,title="Populations", labels=["e1" "e2" "e3"], legend=true)
    p6 = plot(sol,vars=[4,8,12], xlabel="time", linewidth=3, labels=["r1" "r2" "r3"], legend=true)
    p = plot(p1, p5, p2, p4, p3, p6, layout=(3,2), linewidth=3, link=:both)
    savefig(p, "$(prefix)combined.pdf")
    p
end

function makeplots_seird(sol, prefix)
    mkpath(dirname(prefix))
    p1 = plot(sol,vars=[1,2,3,4,5], xlabel="", ylabel="people", linewidth=3,title="Cities", legend=false)
    p2 = plot(sol,vars=[6,7,8,9,10], xlabel="", ylabel="people", linewidth=3, legend=false)
    p3 = plot(sol,vars=[11,12,13,14,15], xlabel="time", ylabel="people", linewidth=3, legend=false)
    p4 = plot(sol,vars=[2,7,12], xlabel="", linewidth=3, labels=["i1" "i2" "i3"], legend=true)
    p5 = plot(sol,vars=[3,8,13], xlabel="", linewidth=3,title="Populations", labels=["e1" "e2" "e3"], legend=true)
    p6 = plot(sol,vars=[5,10,15], xlabel="time", linewidth=3, labels=["d1" "d2" "d3"], legend=true)
    p = plot(p1, p5, p2, p4, p3, p6, layout=(3,2), linewidth=3, link=:both)
    savefig(p, "$(prefix)combined.pdf")
    p
end

T = [
    # City 1 SEIR
    ([1, 2], [3, 2]),
    ([3], [2]), # E→I
    ([2], [4]), # I→R
    # outflow 1→2
    ([1], [5]), # S→S′
    ([2], [6]), # I→I′
    ([3], [7]), # E→E′
    # City 2 SEIR
    ([5, 6], [7, 6]),
    ([7], [6]), # E→I
    ([6], [8]), # I→R
    # outflow 2→3
    ([5], [9]), # S→S′
    ([6], [10]),# I→I′
    ([7], [11]),# E→E′
    # City 3 SEIR
    ([9, 10], [11, 10]),
    ([11], [10]), # E→I
    ([10], [12]) # I→R
]

m = Petri.Model(1:12, T, missing, missing)
nS = length(m.S)
β = ones(Float64, length(T))

tspan = (0.0,60.0)
prob = ODEProblem(fluxes(m), u₀(m, 1,0), tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

u0 = zeros(Float64, length(m.S))
u0[1]  = 10000
u0[5]  = 10000
u0[9]  = 10000
u0[2]  = 1
u0
βseir = [10/sum(u0), 1/2, 1/5]
βtravel = [1/2, 1/2, 1/2]/1000
β = vcat(βseir, βtravel, βseir, βtravel, βseir)
@show β
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
@show sol.u[end]

using Plots
makeplots_seir(sol, "img/estimation/seir")

T = [
    # City 1 SEIₐIRD
    ([1, 2], [3, 2]), # S+Iₐ → E+Iₐ
    ([3], [2]),       # E  → Iₐ
    ([2], [4]),       # Iₐ → I
    ([2], [5]),       # Iₐ → R
    ([4], [5]),       # I  → R
    ([4], [6]),       # I  → D
]

m = Petri.Model(1:6, T, missing, missing)
nS = length(m.S)
β = ones(Float64, length(T))

tspan = (0.0,60.0)
prob = ODEProblem(fluxes(m), u₀(m, 1,0), tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

u0 = zeros(Float64, length(m.S))
u0[1]  = 60000000
u0[2]  = 1
u0

βseir = [10/5sum(u0), 1/2, 1/5, 1/3, 1/4, 1/7]
# βtravel = [1/2, 1/2, 1/2]/1000
# β = vcat(βseir, βtravel, βseir, βtravel, βseir)
β = βseir
@show β
prob = ODEProblem(fluxes(m), u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
@show sol.u[end]

# makeplots_seir(sol, "img/estimation/seiird")

plot(sol, vars=[3, 4, 5, 6], label=["IA" "I" "R" "D"])
# end

using CSV
using Query
using DataFrames

df = CSV.read("../eucovid.csv")

italy = @from i in df begin
    @where i.Country == "Italy" && i.Cases > 0
    @select i
    @collect DataFrame
end

caseplot(italy)

function lossfunc(df::DataFrame, m, u0)
    tspan = (0.0, length(df.Deaths)-1)
    mean(x) = sum(x)/length(x)
    t = tspan[1]:tspan[end]
    f(sol) = begin
        deaths = [u[end] for u in sol.(t)]
        recovr = [u[end-1] for u in sol.(t)]
        infsym = [u[end-2] for u in sol.(t)]
        dt = mean((deaths .- df.Deaths).^2) / mean(df.Deaths)
        rt = mean((recovr .- df.Recoveries).^2) / mean(df.Recoveries)
        it = mean((infsym .- df.Cases).^2) / mean(df.Cases)
        return dt + rt + it
    end

    g(θ) = begin
        prob = ODEProblem(fluxes(m), u0, tspan, θ)
        sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
        neg = sum(abs.(θ[θ .< 0]))
        λ = 10000000
        # if final deaths > final recoveries then penalty
        # deathpenalty = sol.u[end][end] > sol.u[end][end-1] ? 1 : 0
        # iapenalty = sum(tanh.([u[3] - u[2] for u in sol.(t)]).+1)
        iapenalty = 0
        lossval, regval = f(sol), λ*(neg + iapenalty)
        # lossval, regval = f(sol) , λ*(neg + deathpenalty)
        println("loss:$lossval\tregval:$regval")
        lossval+regval
    end
    return g
end

u0_italy = [60000000.0, 0, 1, 0, 0, 0]
ℓ = lossfunc(italy, m, u0_italy)

ℓ(β)

using Optim

function calibrate(k)
    β = [40/5sum(u0), 11/2, 1/1, 2/13, 1/18, 1/35]/k
    prob = ODEProblem(fluxes(m), u0_italy, (0.0,95), β)
    sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
    # plot(sol, vars=[2, 4, 5, 6, 3], label=["IA" "I" "R" "D" "E"])
    # p = fitplot(sol, italy)
    p = plot(sol, vars=[4, 5, 6], label=["I" "R" "D"],
            color=[RGB(1,0,0) RGB(0,1,0) RGB(0.5,0,0.5)])
    @df df scatter!(p,
            0:length(italy.Deaths)-1,
            [:Cases, :Deaths, :Recoveries],
            label=["cases" "deaths" "recoveries"],
            marker=true,
            color=[RGB(1,0,0) RGB(0.5,0,0.5) RGB(0,1,0)]
            )
    return β, sol, p
end

calibrate(15.5)[end]

β = [40/5sum(u0), 11/2, 1/1, 2/13, 1/18, 1/35]/15
# β = [60/5sum(u0), 11/2, 1/1, 2/13, 1/24, 1/39]/2
prob = ODEProblem(fluxes(m), u0_italy, (0.0,64), β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
plot(sol, vars=[2, 4, 5, 6, 3], label=["IA" "I" "R" "D" "E"])
fitplot(sol, italy)

lb, ub = zeros(Float64, 6), 100*ones(Float64, 6)

res = optimize(ℓ, β)
# res = optimize(ℓ, lb, ub, β, Fminbox(NelderMead()))
@show res.minimizer
@show res.minimizer -β

# candidate minimizer
β̂ = [
 9.478916256341902e-9
 0.2915773795849941
 0.09434402982528367
 0.0061044262074341285
 0.016241417062038518
 0.02234216514573622
]

prob = ODEProblem(fluxes(m), u0_italy, (0.0,128), res.minimizer)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
fitplot(sol, italy)

prob = ODEProblem(fluxes(m), u0_italy, (0.0,70), res.minimizer)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
fitplot(sol, df)

# trying to add constraints to the parameter space
realparams(θ) = [θ[1:4]..., θ[5], θ[5]/5]
res = optimize(θ->ℓ(realparams(θ)), β[1:6])
@show res.minimizer
@show res.minimizer -β

prob = ODEProblem(fluxes(m), u0_italy, (0.0,128), realparams(res.minimizer))
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
fitplot(sol, df)


text:
#= Usage:

# Start Julia
julia --color=yes

# Install Reduce.jl package (only necessary once)
import Pkg; Pkg.add("Reduce");

# Load file
julia> using Revise
julia> includet("euler-manufactured.jl")

# Run methods to generate source terms. The source terms that need to be
# implemented can be found in the `source_...` variables of the quote'd output.
# If you want, you can also modify, e.g., the ini method to get different manufactured solutions.
julia> euler2d()
julia> euler3d()
=#

using Reduce
@force using Reduce.Algebra

# Original Reduce code (CompressibleEulerEquations 2D)
#=
clear(γ,f,A,ω,c,ini,rho,rho_v1,rho_v2,rho_v3,rho_e,v1,v2,p,x,y,t,u1,u2,u3,u4);

ini := c + A * sin(ω * (x + y - t));
rho := ini;
rho_v1 := ini;
rho_v2 := ini;
rho_e := ini^2;

v1 := rho_v1 / rho;
v2 := rho_v2 / rho;
p := (γ - 1) * (rho_e - 1/2 * rho * (v1^2 + v2^2));

source_rho    := df(rho, t)    + df(rho_v1, x)           + df(rho_v2, y);
source_rho_v1 := df(rho_v1, t) + df(rho * v1^2 + p, x)   + df(rho * v1 * v2, y);
source_rho_v2 := df(rho_v2, t) + df(rho * v1 * v2, x)    + df(rho * v2^2 + p, y);
source_rho_e  := df(rho_e, t)  + df((rho_e + p) * v1, x) + df((rho_e + p) * v2, y);
=#

function euler2d()
  quote
    ini = c + a * sin(ω * (x + y - t))
    rho = ini
    rho_v1 = ini
    rho_v2 = ini
    rho_e = ini^2

    v1 = rho_v1 / rho
    v2 = rho_v2 / rho
    p = (γ - 1) * (rho_e - 1/2 * rho * (v1^2 + v2^2))

    source_rho    = df(rho, t)    + df(rho_v1, x)           + df(rho_v2, y)
    source_rho_v1 = df(rho_v1, t) + df(rho * v1^2 + p, x)   + df(rho * v1 * v2, y)
    source_rho_v2 = df(rho_v2, t) + df(rho * v1 * v2, x)    + df(rho * v2^2 + p, y)
    source_rho_e  = df(rho_e, t)  + df((rho_e + p) * v1, x) + df((rho_e + p) * v2, y)
  end |> rcall
end

# Original Reduce code (CompressibleEulerEquations 3D)
#=
clear(γ,f,A,ω,c,a1,a2,a3,ini,rho,rho_v1,rho_v2,rho_v3,rho_e,v1,v2,v3,p,x,y,z,t);

ini := c + A * sin(ω * (x + y + z - t));
rho := ini;
rho_v1 := ini;
rho_v2 := ini;
rho_v3 := ini;
rho_e := ini^2;

v1 := rho_v1 / rho;
v2 := rho_v2 / rho;
v3 := rho_v3 / rho;
p := (γ - 1) * (rho_e - 1/2 * rho * (v1^2 + v2^2 + v3^2));

source_rho    := df(rho, t)    + df(rho_v1, x)           + df(rho_v2, y)           + df(rho_v3, z);
source_rho_v1 := df(rho_v1, t) + df(rho * v1^2 + p, x)   + df(rho * v1 * v2, y)    + df(rho * v1 * v3, z);
source_rho_v2 := df(rho_v2, t) + df(rho * v1 * v2, x)    + df(rho * v2^2 + p, y)   + df(rho * v2 * v3, z);
source_rho_v3 := df(rho_v3, t) + df(rho * v1 * v3, x)    + df(rho * v3 * v3, y)    + df(rho * v3^2 + p, z);
source_rho_e  := df(rho_e, t)  + df((rho_e + p) * v1, x) + df((rho_e + p) * v2, y) + df((rho_e + p) * v3, z);
=#

function euler3d()
  quote
    ini = c + a * sin(ω * (x + y + z - t))
    rho = ini
    rho_v1 = ini
    rho_v2 = ini
    rho_v3 = ini
    rho_e = ini^2

    v1 = rho_v1 / rho
    v2 = rho_v2 / rho
    v3 = rho_v3 / rho
    p = (γ - 1) * (rho_e - 1/2 * rho * (v1^2 + v2^2 + v3^2))

    source_rho    = df(rho, t)    + df(rho_v1, x)           + df(rho_v2, y)           + df(rho_v3, z)
    source_rho_v1 = df(rho_v1, t) + df(rho * v1^2 + p, x)   + df(rho * v1 * v2, y)    + df(rho * v1 * v3, z)
    source_rho_v2 = df(rho_v2, t) + df(rho * v1 * v2, x)    + df(rho * v2^2 + p, y)   + df(rho * v2 * v3, z)
    source_rho_v3 = df(rho_v3, t) + df(rho * v1 * v3, x)    + df(rho * v3 * v3, y)    + df(rho * v3^2 + p, z)
    source_rho_e  = df(rho_e, t)  + df((rho_e + p) * v1, x) + df((rho_e + p) * v2, y) + df((rho_e + p) * v3, z)
  end |> rcall
end


text:
module SimultaneousPerturbationStochasticApproximation
using Sobol, Random, Statistics, Dates, Printf, HypothesisTests

export SimpleHeuristic, SPSA, LearningRateUpdate, NoInitializer, RandomInitializer,
SobolInitializer, minimize!, minimize, ConvergenceWelchTest

mutable struct SimpleHeuristic{Ta,U,T}
    a::Ta
    A::T
    c::T
    c_min::T
    a_max::T
    α::T
    γ::T
    update::U
    ninit::Int
end
"""
    SimpleHeuristic(; elementwise = false, A = 500, c = 0., c_min = 1e-3,
                      a_max = 1e4, α = .602, γ = .101, ninit = 100,
                      update = nothing)

`update` could also be [`LearningRateUpdate`](@ref).
"""
function SimpleHeuristic(; elementwise = false, A = 500, c = 0., c_min = 1e-3,
                           a_max = 1e4, α = .602, γ = .101, ninit = 100,
                           update = nothing)
    a = elementwise ? zeros(1) : 0.
    SimpleHeuristic(a, float(A), c, c_min, a_max, α, γ, update, ninit)
end
function init!(rng, s::SimpleHeuristic{Ta}, f, θ, N, maxstepsize) where Ta
    K = floor(Int, s.ninit/3)
    nfevals = 0
    if s.c == 0
        y = [f(θ) for _ in 1:K]
        s.c = max(s.c_min, std(y))
        nfevals += K
    end
    if sum(s.a) == 0
        grad = [estimate_grad(rng, θ, s.c, f)[1] for _ in 1:K]
        nfevals += 2K
        a = min.(s.a_max, abs.(maxstepsize * (s.A + 1)^s.α ./ mean(grad)))
        s.a = Ta <: Number ? mean(a) : a
    end
    s, nfevals
end
update!(h::SimpleHeuristic, θ) = h
mutable struct LearningRateUpdate
    i::Int
    N::Int
    factor::Float64
    range::Vector{Float64}
    θs::Vector{Vector{Float64}}
end
"""
    LearningRateUpdate(lower, upper; pi_a = .7, N = 1000, factor = 1.5)
"""
function LearningRateUpdate(lower, upper; pi_a = .7, N = 1000, factor = 1.5)
    LearningRateUpdate(0, N, factor, pi_a * (upper .- lower), [])
end
function update!(h::SimpleHeuristic{<:AbstractVector,LearningRateUpdate}, θ)
    u = h.update
    if u.i % u.N == u.N - 1
        @views differences = hcat((u.θs[2:end] .- u.θs[1:end-1])...)
        for i in eachindex(h.a)
            θs = [x[i] for x in u.θs]
            if maximum(θs) - minimum(θs) > u.range[i]
                h.a[i] /= u.factor
            elseif pvalue(OneSampleTTest(@view(differences[i, :]))) < .05
                h.a[i] *= u.factor
            end
        end
        empty!(u.θs)
    end
    u.i += 1
    push!(u.θs, copy(θ))
    h
end
struct SPSA{H,I,C,T,Tl,Tu}
    heuristic::H
    initializer::I
    convergencetest::C
    lower::Tl
    upper::Tu
    maxstepsize::Vector{T}
    theta::Vector{T}
end
converged!(::Any, f, θ) = false
mutable struct ConvergenceWelchTest{T}
    i::Int
    N::Int
    K::Int
    m_old::T
    s_old::T
    N_positivetests::Int
    previoustests::Vector{Bool}
end
"""
    ConvergenceWelchTest(; N = 1000, K = 25, T = Float64, N_positivetests = 3)
"""
function ConvergenceWelchTest(; N = 1000, K = 25, T = Float64, N_positivetests = 3)
    ConvergenceWelchTest(0, N, K, zero(T), one(T), N_positivetests, Bool[])
end
function converged!(c::ConvergenceWelchTest, f, θ)
    if c.i % c.N == 0
        f_new = [f(θ) for _ in 1:c.K]
        m_new = mean(f_new)
        s_new = var(f_new)/c.K + eps()
        if length(c.previoustests) != 0
            xbar = c.m_old - m_new
            stderr = sqrt(c.s_old + s_new)
            df = stderr^4/((c.s_old^2 + s_new^2)/(c.K-1))
            p = pvalue(UnequalVarianceTTest(c.K, c.K, xbar, df, stderr, xbar/stderr, 0), tail = :right)
            push!(c.previoustests, p > .05)
            if length(c.previoustests) >= c.N_positivetests && (&)(c.previoustests[end-c.N_positivetests+1:end]...)
                return true
            end
        else
            push!(c.previoustests, false)
        end
        c.m_old = m_new
        c.s_old = s_new
    end
    c.i += 1
    false
end
"""
    NoInitializer()
"""
struct NoInitializer end
init!(rng, spsa::SPSA{<:Any, NoInitializer}, f) = [spsa.theta]
Base.length(::NoInitializer) = 1
fevals(::NoInitializer) = 0
"""
    RandomInitializer()
"""
struct RandomInitializer end
function init!(rng, spsa::SPSA{<:Any, RandomInitializer, <:Any, T}, f) where T
    [rand(rng, T, length(spsa.lower)) .* (spsa.upper .- spsa.lower) .+ spsa.lower]
end
Base.length(::RandomInitializer) = 1
fevals(::RandomInitializer) = 0
"""
    SobolInitializer(; restarts, N)

During initialization, the `SobolInitializer` evaluates the objective function
at `N` points and returns the `restarts` best ones as initial conditions for SPSA.
"""
Base.@kwdef struct SobolInitializer
    restarts::Int
    N::Int
end
function init!(rng, spsa::SPSA{<:Any,SobolInitializer, <:Any, T}, f) where T
    s = SobolSeq(spsa.lower, spsa.upper)
    res = []
    θs = Vector{T}[]
    for θ in Iterators.take(s, spsa.initializer.N)
        push!(res, f(θ))
        push!(θs, copy(θ))
    end
    sp = sortperm(res, rev = true)
    θs[sp[1:spsa.initializer.restarts]]
end
fevals(i::SobolInitializer) = i.N
Base.length(i::SobolInitializer) = i.restarts
"""
    SPSA(; lower, upper,
           heuristic = SimpleHeuristic(elementwise = true,
                                       update = LearningRateUpdate(lower, upper)),
           initializer = SobolInitializer(restarts = 1, N = 100),
           convergencetest = ConvergenceWelchTest(),
           pi_max = .1,
           init = lower)
"""
function SPSA(; lower, upper,
                heuristic = SimpleHeuristic(elementwise = true,
                                            update = LearningRateUpdate(lower, upper)),
                initializer = SobolInitializer(restarts = 1, N = 100),
                convergencetest = ConvergenceWelchTest(),
                pi_max = .1,
                init = lower)
    SPSA(heuristic, initializer, convergencetest, lower, upper,
         pi_max * (upper .- lower), float.(init))
end
struct SPSA_Iterator{A}
    spsa::A
    N::Int
end
Base.length(a::SPSA_Iterator) = a.N
function max_iterations(maxfevals, initializer)
    repetitions = length(initializer)
    evals = initevals = fevals(initializer)
    n = 0
    while evals < maxfevals
        n += 1
        evals += 2 * repetitions
    end
    n, initevals
end
function Base.iterate(a::SPSA_Iterator{<:SPSA{<:SimpleHeuristic}}, n = 1)
    n > a.N && return nothing
    x = a.spsa.heuristic
    A = x.A
    (n, x.a/(n + A)^x.α, x.c/n^x.γ), n + 1
end
_clamp!(x, lo::Nothing, hi::Nothing) = x
_clamp!(x, lo::Number, hi::Number) = clamp!(x, lo, hi)
function _clamp!(x, lo::AbstractVector, hi::AbstractVector)
    @inbounds for i in eachindex(x)
        x[i] = clamp(x[i], lo[i], hi[i])
    end
    x
end

function estimate_grad(rng, θ, c, f)
    d = length(θ)
    δ = 2 * (rand(rng, d) .> .5) .- 1
    θ⁺ = θ .+ c * δ
    θ⁻ = θ .- c * δ
    S⁺ = f(θ⁺)
    S⁻ = f(θ⁻)
    (S⁺ - S⁻)/(2*c) * δ, (S⁺ + S⁻)/2
end
"""
    minimize(f; callback = () -> nothing, maxfevals, verbose = true, kwargs...)

`kwargs` are passed to [`SPSA`](@ref).
"""
function minimize(f; callback = () -> nothing, maxfevals, verbose = true, kwargs...)
    spsa = SPSA(; kwargs...)
    minimize!(spsa, f, callback = callback, maxfevals = maxfevals, verbose = verbose)
end
"""
    minimize!(spsa::SPSA, f; kwargs...) = minimize!(Random.GLOBAL_RNG, spsa, f; kwargs...)
"""
minimize!(spsa::SPSA, f; kwargs...) = minimize!(Random.GLOBAL_RNG, spsa, f; kwargs...)
"""
    minimize!(rng::Random.AbstractRNG,
                   spsa::SPSA{<:Any,<:Any,<:Any,T},
                   f;
                   callback = () -> nothing,
                   maxfevals,
                   gamma_fhat = .9, # just for tracking the learning progress
                   verbose = true) where T
"""
function minimize!(rng::Random.AbstractRNG,
                   spsa::SPSA{<:Any,<:Any,<:Any,T},
                   f;
                   callback = () -> nothing,
                   maxfevals,
                   gamma_fhat = .9, # this is still a hack; should be estimated
                   verbose = true) where T
    start_time = last_time = now()
    verbose && @printf "%6.s %6.s %13.s %10.s %10.s\n" "round" "iter" "elapsed" "fevals" "f"
    results = []
    N, evals = max_iterations(maxfevals, spsa.initializer)
    initial_points = init!(rng, spsa, f)
    _, tune_evals = init!(rng, spsa.heuristic, f, initial_points[1], N, spsa.maxstepsize)
    evals += tune_evals
    restarts = length(initial_points)
    N -= floor(Int, tune_evals/(2*restarts))
    fhat_average = Inf
    for (r, θ) in enumerate(initial_points)
        spsa.theta .= θ # this is a hack to have spsa.theta always point
        θ = spsa.theta  # to the currently updated vector
        for (n, aₙ, cₙ) in SPSA_Iterator(spsa, N)
            grad, fhat = estimate_grad(rng, θ, cₙ, f)
            evals += 2
            fhat_average = isinf(fhat_average) ? fhat : gamma_fhat * fhat_average + (1 - gamma_fhat) * fhat
            @. θ -= clamp(aₙ * grad, -spsa.maxstepsize, spsa.maxstepsize)
            _clamp!(θ, spsa.lower, spsa.upper)
            callback()
            if verbose && (n % div(N, 50) == 0 || now() - last_time > Second(5))
                last_time = now()
                @printf "%4.g/%g %6.g %13.s %10.g %10.f\n" r restarts n round(last_time - start_time, Second) evals fhat_average
            end
            update!(spsa.heuristic, θ)
            converged!(spsa.convergencetest, f, θ) && break # these evals are not counted currently
        end
        push!(results, (x = copy(θ), f = fhat_average, evals = evals))
    end
    sort(results, by = x -> x.f)
end

end # module


text:
struct HuberAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int,Int}
    M::Real

    function HuberAtom(x::AbstractExpr, M::Real)
        if sign(x) == ComplexSign()
            error("Arguemt must be real")
        elseif M <= 0
            error("Huber parameter must by a positive scalar")
        end
        children = (x,)
        return new(:huber, hash((children, M)), children, x.size, M)
    end
end

function sign(x::HuberAtom)
    return Positive()
end

function monotonicity(x::HuberAtom)
    return (Nondecreasing() * sign(x.children[1]),)
end

function curvature(x::HuberAtom)
    return ConvexVexity()
end

function evaluate(x::HuberAtom)
    c = evaluate(x.children[1])
    for i in 1:length(c)
        if c[i] <= x.M
            c[i] = c[i]^2
        else
            c[i] = 2 * x.M * c[i] - x.M^2
        end
    end
    return c
end

function conic_form!(x::HuberAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, x)
        c = x.children[1]
        s = Variable(c.size)
        n = Variable(c.size)

        # objective given by s.^2 + 2 * M * |n|
        objective =
            conic_form!(square(s) + 2 * x.M * abs(n), unique_conic_forms)
        conic_form!(c == s + n, unique_conic_forms)

        cache_conic_form!(unique_conic_forms, x, objective)
    end
    return get_conic_form(unique_conic_forms, x)
end

huber(x::AbstractExpr, M::Real = 1.0) = HuberAtom(x, M)


text:
# check the Bernoulli numbers and polynomials
#    this is mostly just a check that code is correct
#    most errors are arising because of errors in the functions being used
using Polylogarithms
using PyPlot
using SpecialFunctions
import Base: parse
using Printf

Q = Polylogarithms.Q
c1 = Polylogarithms.c
c2 = Polylogarithms.c_crandall
b = Polylogarithms.b_crandall
f = Polylogarithms.f_crandall
g = Polylogarithms.g_crandall

# check f
@printf("Testing f\n")
for k=0:5
    @printf(" k = %d\n", k)
    @printf("   f(k,1) = %f, -H_k= %f\n", f(k,1), -harmonic(k) )
    @printf("   f(k,2) = %f, (H_k^2 + H_{k,2})/2= %f\n", f(k,2), (harmonic(k)^2 + harmonic(k,2))/2 )
end
@printf("\n")

# # check g against numerical derivates of Gamma at 1
@printf("Testing g\n")
function G(t::Integer, x::Float64; dx::Float64 = 1.0e-3)
    if t==0
        return 1
    elseif t==1
        return ( SpecialFunctions.gamma(x+dx) - SpecialFunctions.gamma(x) )/dx
    elseif t>1
        return ( G(t-1,x+dx) - G(t-1,x) ) /dx
    else
        throw(DomainError(t))
    end
end
for t=0:5
    @printf(" t = %d, ", t)
    @printf(" G^{t}(1) = %.7f,  g(t) = %.7f\n",  G(t,1.0), g(t) )
end
@printf("\n")

# check c
z = complex(-1.0/2.0)
μ = log(z)
ℒ = log(-μ)
for k=0:5
    @printf(" k = %d\n", k)
    for j=0:2
        @printf("  j = %d,  ", j)
        @printf("  c_explicit = %.12f + %.12fi, ", real(c1(k, j, ℒ)), imag(c1(k, j, ℒ)) )
        @printf("  c_crandall = %.12f + %.12fi, ", real(c2(k, j, ℒ)), imag(c2(k, j, ℒ)) )
        @printf("\n")
   end
   @printf("\n")
end

c2(1, 3, ℒ)


text:
using OrdinaryDiffEq
using Trixi
using Printf
using PrettyTables

######### EINSTELLUNGEN ############
N_vec = [3]
Nq_vec = [2, 4, 8, 16]
CFL = 0.1
# timespan for periodic solution
tspan = (0.0, 2.0)
mu = 0.0001
latex = false

surface_flux = FluxPlusDissipation(flux_chandrashekar, DissipationLocalLaxFriedrichs(max_abs_speed_naive)) # flux_lax_friedrichs
volume_flux  = flux_chandrashekar
volume_integral = VolumeIntegralFluxDifferencing(volume_flux)

# mapping as described in the worksheet
function mapping(xi_, eta_)

    xi = xi_
    eta = eta_

    x = xi + 0.15 * cos(0.5 * pi * xi) * cos((3 / 2) * pi * eta)
    y = eta + 0.15 * cos(2 * pi * xi) * cos(0.5 * pi * eta)

    return SVector(x, y)
end

###############################################################################
# semidiscretization of the compressible Euler equations

# equations = CompressibleEulerEquations2D(1.4)
equations = CompressibleEulerEquations2D(1.4, viscous=true, mu=0.00001)

# choose the initilial condition and source terms
# here we use the already in Trixi used ones
initial_condition =initial_condition_convergence_test #initial_condition_convergence_test_gassner #
source_terms = source_terms_convergence_test #source_terms_convergence_test_gassner #

# coordinates neccessary for periodic solution with periodic boundaries
coordinates_min = (0.0, 0.0)
coordinates_max = (2.0,  2.0)

####

    # Table preparation
kwargs = Dict{Symbol,Any}(
        :formatters => (
            ft_printf("%.5e", [2]),
            ft_printf("%.2f", [3])
        )
    )

if latex
    kwargs[:backend] = :latex
    kwargs[:highlighters] = LatexHighlighter(
            (data, i, j) -> j > 1,
            (data, i, j, str) -> "\$\\num{$str}\$"
        )
end

error = Vector{Float64}(undef, length(Nq_vec))
for N in N_vec

    basis = LobattoLegendreBasis(N)
    solver = DGSEM(basis, surface_flux, volume_integral)
    println(typeof(solver.volume_integral))
    # solver = DGSEM(polydeg=N, surface_flux=surface_flux, volume_integral=volume_integral)
    for j in eachindex(Nq_vec)
        Nq = Nq_vec[j]

        cells_per_dimension = (Nq, Nq)

        # mesh = CurvedMesh(cells_per_dimension, mapping)#(-1.0, -1.0), (1.0, 1.0))
        mesh = CurvedMesh(cells_per_dimension, coordinates_min, coordinates_max)
        boundary_conditions = (x_neg = boundary_condition_periodic,
                       x_pos = boundary_condition_periodic,
                       y_neg = boundary_condition_periodic,
                       y_pos = boundary_condition_periodic)

        semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions, source_terms=source_terms_convergence_test)

            ###############################################################################
            # ODE solvers, callbacks etc.

        ode = semidiscretize(semi, tspan)

        summary_callback = SummaryCallback()
        analysis_interval = 1000

        analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true)

        stepsize_callback = StepsizeCallback(cfl=CFL)

        callbacks = CallbackSet(
                                    analysis_callback,
                                    stepsize_callback)

            ###############################################################################
            # run the simulation

        sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
                        dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
                        save_everystep=false, callback=callbacks);
        summary_callback() # print the timer summary

        error[j] = maximum(abs.(sol.u[1] - sol.u[2])) # maximum(abs.(exact - sol))
    end

        # Compute EOC
    eoc = vcat("", [log(error[j + 1] / error[j]) /
            log(Nq_vec[j] / Nq_vec[j + 1]) for j in 1:length(Nq_vec) - 1])
    println("N = $N")

    pretty_table(hcat(Nq_vec, error, eoc), ["Nq", "Error", "EOC"]; kwargs...)

    mean = sum(eoc[2:end]) / (length(Nq_vec) - 1)
    @printf("Mean EOC: %.2f\n", mean)
end

max_stars_repo_path: src/WR2/MS5 EOC.jl

max_issues_repo_path: src/WR2/MS5 EOC.jl

max_forks_repo_path: src/WR2/MS5 EOC.jl


text:
export totalorder

# totalorder: Compute indices of total degree
# polynomial expansion. Output is an (ncoeff x d) matrix
# where 0 corresponds to the constant function
"""
$(TYPEDSIGNATURES)

Compute indices of total degree polynomial expansion.
"""
function totalorder(order_list::Array{Int64,1})

    if isempty(order_list)
        return zeros(Int64,0,0)
    end

    # Determine dimension and max_order
    d = size(order_list, 1)
    max_order = maximum(order_list)

    # Initialize multi_idx with zeros
    midxs_new = zeros(Int64, 1, d)
    multi_idxs = deepcopy(midxs_new)

    # Initialize midxs_old_set
    if sum(midxs_new) < max_order
        midxs_old_set = deepcopy(midxs_new)
    else
        midxs_old_set = zeros(Int64, 0, d)
    end

    # Generate higher order multi-indices
    for i = 1:max_order
        # Initialize empty set of starting multi-indices
        midxs_new_set = zeros(Int64, 0, d)

        # Extract each multi_idx in midxs_old_set
        for j = 1:size(midxs_old_set,1)
            midxs_old_j = deepcopy(midxs_old_set[j,:])
            # Expand index set along each direction
            for k = 1:d

                # If allowable, add new multi_idx
                if midxs_old_j[k] < order_list[k]
                    midx_new = deepcopy(midxs_old_j)
                    midx_new[k] += 1
                    multi_idxs = vcat(multi_idxs, reshape(midx_new,(1,d)))

                    # If boundary of orders isn't added, expand
                    # in the next iteration by adding to set
                    if sum(midx_new) < max_order
                        midxs_new_set = vcat(midxs_new_set, reshape(midx_new,(1,d)))
                    end
                end
            end
        end

        # Overwrite midxs_old_set for next iteration
        midxs_old_set = deepcopy(midxs_new_set)
    end

    # Remove duplicates in multi_idxs
    multi_idxstmp = unique(eachslice(multi_idxs; dims = 1))
    cardinal = size(multi_idxstmp,1)
    multi_idxs = zeros(Int64, cardinal,d)

    @inbounds for i=1:cardinal
        multi_idxs[i,:] = multi_idxstmp[i]
    end

    multi_idxs = sortslices(multi_idxs; dims = 1)

    return multi_idxs

end


text:
@require Juno begin

import Media

Juno.render(i::Juno.PlotPane, n::AbstractWidget) =
    Juno.render(i, WebIO.render(n))

Juno.render(i::Juno.Editor, n::AbstractWidget) =
    Juno.render(i, WebIO.render(n))

Juno.media(AbstractWidget, Media.Graphical)

end


text:
using MathematicalSystems, LazySets

export ContinuousSystem,
       DiscreteSystem,
       add_dimension,
       next_set

# name alises
const LCS = LinearContinuousSystem
const LDS = LinearDiscreteSystem
const CLCCS = ConstrainedLinearControlContinuousSystem
const CLCDS = ConstrainedLinearControlDiscreteSystem

import Base: *

*(M::AbstractMatrix, input::ConstantInput) =  ConstantInput(M * input.U)

# no input: x' = Ax, x(0) = X0
ContinuousSystem(A::AbstractMatrix, X0::LazySet) = IVP(LCS(A), X0)
DiscreteSystem(A::AbstractMatrix, X0::LazySet) = IVP(LDS(A), X0)

# with constant input: x' = Ax + u, x(0) = X0, u(t) = U
ContinuousSystem(A::AbstractMatrix, X0::LazySet, U::ConstantInput) = IVP(CLCCS(A, eye(A), nothing, U), X0)
ContinuousSystem(A::AbstractMatrix, X0::LazySet, U::LazySet) = ContinuousSystem(A, X0, ConstantInput(U))

DiscreteSystem(A::AbstractMatrix, X0::LazySet, U::ConstantInput) = IVP(CLCDS(A, eye(A), nothing, U), X0)
DiscreteSystem(A::AbstractMatrix, X0::LazySet, U::LazySet) = DiscreteSystem(A, X0, ConstantInput(U))

# with time-varying input: x' = Ax + u, x(0) = X0, u(t) = U(t)
ContinuousSystem(A::AbstractMatrix, X0::LazySet, U::VaryingInput) = IVP(CLCCS(A, eye(A), nothing, U), X0)
ContinuousSystem(A::AbstractMatrix, X0::LazySet, U::Vector{<:LazySet}) = ContinuousSystem(A, X0, VaryingInput(U))

DiscreteSystem(A::AbstractMatrix, X0::LazySet, U::VaryingInput) = IVP(CLCDS(A, eye(A), nothing, U), X0)
DiscreteSystem(A::AbstractMatrix, X0::LazySet, U::Vector{<:LazySet}) = DiscreteSystem(A, X0, VaryingInput(U))

# convenience functions
next_set(inputs::ConstantInput) = next(inputs, 1)[1]
next_set(inputs::AbstractInput, state::Int64) = next(inputs, state)[1]

"""
    add_dimension(A::AbstractMatrix, m=1)

Adds an extra zero row and column to a matrix.

### Input

- `A` -- matrix
- `m` -- (optional, default: `1`) the number of extra dimensions

### Examples

```jldoctest add_dimension_test
julia> A = [0.4 0.25; 0.46 -0.67]
2×2 Array{Float64,2}:
 0.4    0.25
 0.46  -0.67
julia> add_dimension(A)
3×3 Array{Float64,2}:
 0.4    0.25  0.0
 0.46  -0.67  0.0
 0.0    0.0   0.0
```

```jldoctest add_dimension_test
julia> add_dimension(A, 2)
4×4 Array{Float64,2}:
 0.4    0.25  0.0  0.0
 0.46  -0.67  0.0  0.0
 0.0    0.0   0.0  0.0
 0.0    0.0   0.0  0.0
```
"""
function add_dimension(A::AbstractMatrix, m=1)
    n = size(A, 1)
    return vcat(hcat(A, zeros(n, m)), zeros(m, n+m))
end

"""
    add_dimension(X::LazySet, m=1)::LazySet

Adds an extra dimension to a LazySet through a Cartesian product.

### Input

- `X` -- a lazy set
- `m` -- (optional, default: `1`) the number of extra dimensions

### Examples

```jldoctest add_dimension_set
julia> X = BallInf(ones(9), 0.5);

julia> dim(X)
9

julia> Xext = add_dimension(X);

julia> dim(Xext)
10

julia> X = ZeroSet(4);

julia> dim(add_dimension(X))
5

julia> typeof(X)
LazySets.ZeroSet{Float64}
```

More than one dimension can be added passing the second argument:

```jldoctest add_dimension_set
julia> Xext = add_dimension(BallInf(zeros(10), 0.1), 4);

julia> dim(Xext)
14
```

### Notes

In the special case that the given set is a zero set, instead of cartesian product
a new zero set with extended dimensions is returned.
"""
function add_dimension(X::LazySet, m=1)::LazySet
    return X * ZeroSet(m)
end

function add_dimension(X::ZeroSet, m=1)::ZeroSet
    return ZeroSet(dim(X)+m)
end

"""
    add_dimension(cs, m=1)

Adds an extra dimension to a continuous system.

### Input

- `cs` -- continuous system
- `m` -- (optional, default: `1`) the number of extra dimensions

### Examples

```jldoctest add_dimension_cont_sys
julia> using MathematicalSystems
julia> A = sprandn(3, 3, 0.5);
julia> X0 = BallInf(zeros(3), 1.0);
julia> s = ContinuousSystem(A, X0);
julia> sext = add_dimension(s);
julia> statedim(sext)
4
```

If there is an input set, it is also extended:

```jldoctest add_dimension_cont_sys
julia> U = Ball2(ones(3), 0.1);
julia> s = ContinuousSystem(A, X0, U);
julia> sext = add_dimension(s);
julia> statedim(sext)
4
julia> dim(next_set(inputset(sext)))
4
```

Extending a system with a varying input set:

If there is an input set, it is also extended:

```jldoctest add_dimension_cont_sys
julia> U = [Ball2(ones(3), 0.1 * i) for i in 1:3];
julia> s = ContinuousSystem(A, X0, U);
julia> sext = add_dimension(s);
julia> statedim(sext)
4
julia> dim(next_set(inputset(sext), 1))
4
```

Extending a varing input set with more than one extra dimension:

```jldoctest add_dimension_cont_sys
julia> sext = add_dimension(s, 7);
julia> statedim(sext)
10
julia> dim(next_set(inputset(sext), 1))
10
```
"""
function add_dimension(cs, m=1)
    Aext = add_dimension(cs.s.A, m)
    X0ext = add_dimension(cs.x0, m)
    if method_exists(inputset, Tuple{typeof(cs.s)})
        Uext = map(x -> add_dimension(x, m), inputset(cs))
        return ContinuousSystem(Aext, X0ext, Uext)
    else
        return ContinuousSystem(Aext, X0ext)
    end
end