text:
abstract type UAVDynamicsModel end
abstract type UAVState end
abstract type UAVAction end

"""
Represents state of simple 2D multirotor, with position and velocity in each direction.
"""
struct MultiRotorUAVState <: UAVState
    x::Float64
    y::Float64
    xdot::Float64
    ydot::Float64
end

function get_position(state::MultiRotorUAVState)
    return Point(state.x, state.y)
end

function get_speed(state::MultiRotorUAVState)
    return sqrt(state.xdot^2 + state.ydot^2)
end

function rel_state(source::MultiRotorUAVState, target::MultiRotorUAVState)
    return MultiRotorUAVState(source.x - target.x, source.y - target.y, source.xdot, source.ydot)
end
"""
Represents control action for simple 2D multirotor, with acceleration in each direction.
"""
struct MultiRotorUAVAction <: UAVAction
    xddot::Float64
    yddot::Float64
end

"""
Simulates physical dynamics of UAV, mapping state and control to next state
"""
struct MultiRotorUAVDynamicsModel <: UAVDynamicsModel
    timestep::Float64
    noise::Distributions.ZeroMeanDiagNormal
    params::Parameters
end

"""
    MultiRotorUAVDynamicsModel(t::Float64, params::Parameters, sig::Float64)

Define a dynamics model where t is the timestep duration, params is the collection of system parameters,
and sig is the standard deviation along each axis (sig_xy) is vector of standard deviations for each corresponding axis
"""
function MultiRotorUAVDynamicsModel(params::Parameters)
    # Generate diagonal covariance matrix
    noise = Distributions.MvNormal([params.scale_params.ACC_NOISE_STD, params.scale_params.ACC_NOISE_STD])
    return MultiRotorUAVDynamicsModel(params.time_params.MDP_TIMESTEP, noise, params)
end

function MultiRotorUAVDynamicsModel(sig_xy::StaticVector{2,Float64}, params::Parameters)
    # Generate diagonal covariance matrix
    noise = Distributions.MvNormal([sig_xy[1], sig_xy[2]])
    return MultiRotorUAVDynamicsModel(t, noise, params)
end

"""
    get_uav_dynamics_actions(model::MultiRotorUAVDynamicsModel)

Compiles a vector of all multirotor acceleration actions within limits, based on the resolution parameters
"""
function get_uav_dynamics_actions(::Type{MultiRotorUAVAction}, params::Parameters)

    multirotor_actions = Vector{MultiRotorUAVAction}(undef, 0)
    acc_vals = range(-params.scale_params.ACCELERATION_LIM, stop=params.scale_params.ACCELERATION_LIM, length=params.scale_params.ACCELERATION_NUMVALS)

    for xddot in acc_vals
        for yddot in acc_vals
            push!(multirotor_actions, MultiRotorUAVAction(xddot, yddot))
        end
    end

    return multirotor_actions
end

"""
    generate_start_state(model::MultiRotorUAVDynamicsModel, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

Generate a MultiRotorUAVState with a random location inside the grid and at rest
"""
function generate_start_state(model::MultiRotorUAVDynamicsModel, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    x = rand(rng, Uniform(-model.params.scale_params.XY_LIM, model.params.scale_params.XY_LIM))
    y = rand(rng, Uniform(-model.params.scale_params.XY_LIM, model.params.scale_params.XY_LIM))

    xdot = 0.
    ydot = 0.

    return MultiRotorUAVState(x, y, xdot, ydot)
end

function generate_start_state(::Type{MultiRotorUAVState}, params::Parameters, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    x = rand(rng, Uniform(-params.scale_params.XY_LIM, params.scale_params.XY_LIM))
    y = rand(rng, Uniform(-params.scale_params.XY_LIM, params.scale_params.XY_LIM))

    xdot = 0.
    ydot = 0.

    return MultiRotorUAVState(x, y, xdot, ydot)
end

"""
    get_state_at_rest(model::MultiRotorUAVDynamicsModel, p::Point)

Given a Point instance, generate the corresponding MultiRotorUAVState at rest at that position
"""
function get_state_at_rest(::Type{MultiRotorUAVState}, p::Point)
    return MultiRotorUAVState(p.x, p.y, 0.0, 0.0)
end

"""
get_relative_state_to_goal(model::MultiRotorUAVDynamicsModel, goal_pos::Point, state::MultiRotorUAVState)

Given some goal position, get the MultiRotorUAVState with relative position and own velocity
"""
function get_relative_state_to_goal(::Type{MultiRotorUAVState}, goal_pos::Point, state::MultiRotorUAVState)
    return MultiRotorUAVState(state.x - goal_pos.x, state.y - goal_pos.y, state.xdot, state.ydot)
end

"""
    apply_controls(model::MultiRotorUAVDynamicsModel, state::MultiRotorUAVState, xddot::Float64, yddot::Float64)

Propagate the MultiRotorUAVState through the dynamics model (without noise)
"""
function apply_controls(model::MultiRotorUAVDynamicsModel, state::MultiRotorUAVState, xddot::Float64, yddot::Float64)

    # Update position and velocity exactly
    xdot = clamp(state.xdot + xddot*model.timestep, -model.params.scale_params.XYDOT_LIM, model.params.scale_params.XYDOT_LIM)
    ydot = clamp(state.ydot + yddot*model.timestep, -model.params.scale_params.XYDOT_LIM, model.params.scale_params.XYDOT_LIM)

    # Get true effective accelerations
    true_xddot = (xdot - state.xdot)/model.timestep
    true_yddot = (ydot - state.ydot)/model.timestep

    x = state.x + state.xdot*model.timestep + 0.5*true_xddot*(model.timestep^2)
    x = clamp(x, -model.params.scale_params.XY_LIM, model.params.scale_params.XY_LIM)

    y = state.y + state.ydot*model.timestep + 0.5*true_yddot*(model.timestep^2)
    y = clamp(y, -model.params.scale_params.XY_LIM, model.params.scale_params.XY_LIM)

    return MultiRotorUAVState(x, y, xdot, ydot)
end

function apply_controls(params::Parameters, state::MultiRotorUAVState, xddot::Float64, yddot::Float64)

    timestep = params.time_params.MDP_TIMESTEP

    # Update position and velocity exactly
    xdot = clamp(state.xdot + xddot*timestep, -params.scale_params.XYDOT_LIM, params.scale_params.XYDOT_LIM)
    ydot = clamp(state.ydot + yddot*timestep, -params.scale_params.XYDOT_LIM, params.scale_params.XYDOT_LIM)

    # Get true effective accelerations
    true_xddot = (xdot - state.xdot)/timestep
    true_yddot = (ydot - state.ydot)/timestep

    x = state.x + state.xdot*timestep + 0.5*true_xddot*(timestep^2)
    x = clamp(x, -params.scale_params.XY_LIM, params.scale_params.XY_LIM)

    y = state.y + state.ydot*timestep + 0.5*true_yddot*(timestep^2)
    y = clamp(y, -params.scale_params.XY_LIM, params.scale_params.XY_LIM)

    return MultiRotorUAVState(x, y, xdot, ydot)
end

"""
    next_state(model::MultiRotorUAVDynamicsModel, state::MultiRotorUAVState, action::MultiRotorUAVAction, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

Generate the next state by propagating the action noisily with the dynamics
"""
function next_state(model::MultiRotorUAVDynamicsModel, state::MultiRotorUAVState, action::MultiRotorUAVAction, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    # Sample acceleration noise vector
    noiseval = rand(rng,model.noise)

    # Augment actual acceleration with noise
    xddot = action.xddot + noiseval[1]
    yddot = action.yddot + noiseval[2]

    return apply_controls(model, state, xddot, yddot)
end

function next_state(params::Parameters, state::MultiRotorUAVState, action::MultiRotorUAVAction, rng::RNG=Random.GLOBAL_RNG) where {RNG <: AbstractRNG}

    # Sample acceleration noise vector
    noise = Distributions.MvNormal([params.scale_params.ACC_NOISE_STD, params.scale_params.ACC_NOISE_STD])
    noiseval = rand(rng, noise)

    # Augment actual acceleration with noise
    xddot = action.xddot + noiseval[1]
    yddot = action.yddot + noiseval[2]

    return apply_controls(params, state, xddot, yddot)
end

"""
    function dynamics_cost(model::MultiRotorUAVDynamicsModel, state::MultiRotorUAVState, next_state::MultiRotorUAVState)

Compute the energy/dynamics cost of going from one state to another, based on the cost parameters in the model.
N.B - In the most general case, cost should depend on the action too.
"""
function dynamics_cost(model::MultiRotorUAVDynamicsModel, state::MultiRotorUAVState, next_state::MultiRotorUAVState)
    old_point = Point(state.x, state.y)
    new_point = Point(next_state.x, next_state.y)
    dyn_dist = point_dist(old_point, new_point)

    cost = 0.0

    # Compute cost due to flying or due to hovering
    if dyn_dist < model.timestep*model.params.scale_params.EPSILON &&
        sqrt(next_state.xdot^2 + next_state.ydot^2) < model.timestep*model.params.scale_params.EPSILON
        cost += model.params.cost_params.HOVER_COEFFICIENT*model.timestep
    else
        cost += model.params.cost_params.FLIGHT_COEFFICIENT*dyn_dist
    end

    return cost
end

function dynamics_cost(params::Parameters, state::MultiRotorUAVState, next_state::MultiRotorUAVState)
    old_point = Point(state.x, state.y)
    new_point = Point(next_state.x, next_state.y)
    dyn_dist = point_dist(old_point, new_point)

    cost = 0.0

    # Compute cost due to flying or due to hovering
    if dyn_dist < params.time_params.MDP_TIMESTEP*params.scale_params.EPSILON &&
        sqrt(next_state.xdot^2 + next_state.ydot^2) < params.time_params.MDP_TIMESTEP*params.scale_params.EPSILON
        cost += params.cost_params.HOVER_COEFFICIENT*params.time_params.MDP_TIMESTEP
    else
        cost += params.cost_params.FLIGHT_COEFFICIENT*dyn_dist
    end

    return cost
end


text:
module Colors

using FixedPointNumbers
using Reexport
using Printf

@reexport using ColorTypes

# TODO: why these types are defined here? Can they move to ColorTypes.jl?
AbstractAGray{C<:AbstractGray,T} = AlphaColor{C,T,2}
AbstractGrayA{C<:AbstractGray,T} = ColorAlpha{C,T,2}

import Base: ==, +, -, *, /
import Base: convert, eltype, hex, isless, range, show, typemin, typemax

# Additional exports, not exported by ColorTypes
export weighted_color_mean,
       hex, @colorant_str,
       protanopic, deuteranopic, tritanopic,
       distinguishable_colors, whitebalance,
       colordiff, DE_2000, DE_94, DE_JPC79, DE_CMC, DE_BFD, DE_AB, DE_DIN99, DE_DIN99d, DE_DIN99o,
       MSC, sequential_palette, diverging_palette, colormap,
       colormatch, CIE1931_CMF, CIE1964_CMF, CIE1931J_CMF, CIE1931JV_CMF

# Early utilities
include("utilities.jl")

# Include other module components
include("conversions.jl")
include("promotions.jl")
include("algorithms.jl")
include("parse.jl")
include("differences.jl")
include("colormaps.jl")
include("display.jl")
include("colormatch.jl")

end # module

ext: jl


text:
abstract type AbstractObjectDetection{T <: Real} end

@enum DetectorType HUMAN MACHINE

struct ObjectDetectionCore
    class_name::String
    confidence::Union{Float32, Nothing}
    image_width::Union{Int, Nothing}
    image_height::Union{Int, Nothing}
    detector_type::DetectorType
    detector_name::String
end

function ObjectDetectionCore(class_name::String, detector_type::DetectorType, detector_name::String)
    ObjectDetectionCore(class_name, nothing, nothing, nothing, detector_type, detector_name)
end

function ObjectDetectionCore(class_name::String, confidence::Float32, detector_type::DetectorType, detector_name::String)
    ObjectDetectionCore(class_name, confidence, nothing, nothing, detector_type, detector_name)
end

class_name(detection::AbstractObjectDetection) = detection.core.class_name
confidence(detection::AbstractObjectDetection) = detection.core.confidence
image_width(detection::AbstractObjectDetection) = detection.core.image_width
image_height(detection::AbstractObjectDetection) = detection.core.image_height
detector_type(detection::AbstractObjectDetection) = detection.core.detector_type
detector_name(detection::AbstractObjectDetection) = detection.core.detector_name

function centroid(detection::AbstractObjectDetection{T}) :: Point2{T} where T
    error("No implementation for $(typeof(detection))")
end

function bounding_box_detection(detection::AbstractObjectDetection{T}) :: BoundingBoxDetection{T} where T
    error("No implementation for $(typeof(detection))")
end

function bounding_box(detection::AbstractObjectDetection{T}) :: Rect2{T} where T
    return bounding_box_detection(detection).rect
end

function iou(d1::AbstractObjectDetection{T}, d2::AbstractObjectDetection{T}) where T
    bb1 = bounding_box(d1)
    bb2 = bounding_box(d2)
    left_1, top_1 = bb1.origin; right_1, bottom_1 = bb1.origin + bb1.widths
    left_2, top_2 = bb2.origin; right_2, bottom_2 = bb2.origin + bb2.widths

    if bottom_1 <= top_2 || bottom_2 <= top_1 || right_1 <= left_2 || right_2 <= left_1
        return 0.0
    end

    area_1 = (bottom_1 - top_1) * (right_1 - left_1)
    area_2 = (bottom_2 - top_2) * (right_2 - left_2)

    a = min(right_1, right_2)
    b = max(top_1, top_2)
    c = max(left_1, left_2)
    d = min(bottom_1, bottom_2)
    intersection = (d - b) * (a - c)
    union = area_1 + area_2 - intersection

    return intersection / union
end


text:
using LazySets.Arrays: _rationalize, append_zeros, prepend_zeros

for N in [Float64, Float32, Rational{Int}]
    v1 = N[0, 4, 0]
    v2 = sparsevec(N[2], N[4], 3)
    v3 = SingleEntryVector(2, 3, N(4))
    @test v1 == v2 == v3
    v1a = append_zeros(v1, 2)
    v1p = prepend_zeros(v1, 2)
    @test v1a == append_zeros(v2, 2) == append_zeros(v3, 2) == N[0, 4, 0, 0, 0]
    @test v1p == prepend_zeros(v2, 2) == prepend_zeros(v3, 2) == N[0, 0, 0, 4, 0]
end

for N in [Float32, Float64]
    x = [N(1), N(2), N(3)]
    out = [1//1, 2//1, 3//1]
    @test _rationalize(x) == out
    v = _rationalize(BigInt, x)
    @test v[1] isa Rational{BigInt}
    @test _rationalize(x, tol=2*eps(N)) == out
end


text:
# _Constant--DecisionRef
Base.:(+)(lhs::_Constant, rhs::DecisionRef) = DAE(_VAE(convert(Float64, lhs)), _DAE(0.0, rhs => +one(Float64)), _KAE(0.0))
Base.:(-)(lhs::_Constant, rhs::DecisionRef) = DAE(_VAE(convert(Float64, lhs)), _DAE(0.0, rhs => -one(Float64)), _KAE(0.0))
Base.:(*)(lhs::_Constant, rhs::DecisionRef) = DAE(_VAE(0.0), _DAE(0.0, rhs => lhs), _KAE(0.0))

# _Constant--KnownRef
Base.:(+)(lhs::_Constant, rhs::KnownRef) = DAE(_VAE(convert(Float64, lhs)), _DAE(0.0), _KAE(0.0, rhs => +one(Float64)))
Base.:(-)(lhs::_Constant, rhs::KnownRef) = DAE(_VAE(convert(Float64, lhs)), _DAE(0.0), _KAE(0.0, rhs => -one(Float64)))
Base.:(*)(lhs::_Constant, rhs::KnownRef) = DAE(_VAE(0.0), _DAE(0.0), _KAE(0.0, rhs => lhs))

# _Constant--DecisionAffExpr{C}
Base.:(+)(lhs::_Constant, rhs::DecisionAffExpr{C}) where C =
    DecisionAffExpr(convert(C,lhs) + rhs.variables,
                    copy(rhs.decisions),
                    copy(rhs.knowns))
Base.:(-)(lhs::_Constant, rhs::DecisionAffExpr{C}) where C =
    DecisionAffExpr(convert(C,lhs) - rhs.variables,
                    -rhs.decisions,
                    -rhs.knowns)
Base.:(*)(lhs::_Constant, rhs::DecisionAffExpr{C}) where C =
    DecisionAffExpr(convert(C,lhs) * rhs.variables,
                    convert(C,lhs) * rhs.decisions,
                    convert(C,lhs) * rhs.knowns)

# _Constant--DecisionQuadExpr{C}
Base.:(+)(lhs::_Constant, rhs::DecisionQuadExpr{C}) where C =
    DecisionQuadExpr(convert(C,lhs) + rhs.variables,
                     copy(rhs.decisions),
                     copy(rhs.knowns),
                     copy(rhs.cross_terms),
                     copy(rhs.known_variable_terms),
                     copy(rhs.known_decision_terms))
Base.:(-)(lhs::_Constant, rhs::DecisionQuadExpr{C}) where C =
    DecisionQuadExpr(convert(C,lhs) - rhs.variables,
                     copy(rhs.decisions),
                     copy(rhs.knowns),
                     copy(rhs.cross_terms),
                     copy(rhs.known_variable_terms),
                     copy(rhs.known_decision_terms))
Base.:(*)(lhs::_Constant, rhs::DecisionQuadExpr{C}) where C =
    DecisionQuadExpr(convert(C,lhs) * rhs.variables,
                     convert(C,lhs) * rhs.decisions,
                     convert(C,lhs) * rhs.knowns,
                     _map_cross_terms(c -> convert(C,lhs) * c,
                                      copy(rhs.cross_terms)),
                     _map_cross_terms(c -> convert(C,lhs) * c,
                                      copy(rhs.known_variable_terms)),
                     _map_cross_terms(c -> convert(C,lhs) * c,
                                      copy(rhs.known_decision_terms)))

#=
    VariableRef
=#

# VariableRef--DecisionRef
Base.:(+)(lhs::VariableRef, rhs::DecisionRef) = DAE(_VAE(0.0, lhs => 1.0), _DAE(0.0, rhs =>  1.0), _KAE(0.0))
Base.:(-)(lhs::VariableRef, rhs::DecisionRef) = DAE(_VAE(0.0, lhs => 1.0), _DAE(0.0, rhs => -1.0), _KAE(0.0))
function Base.:(*)(lhs::VariableRef, rhs::DecisionRef)
    result = zero(DQE)
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# VariableRef--KnownRef
Base.:(+)(lhs::VariableRef, rhs::KnownRef) = DAE(_VAE(0.0, lhs => 1.0), _DAE(0.0), _KAE(0.0, rhs =>  1.0))
Base.:(-)(lhs::VariableRef, rhs::KnownRef) = DAE(_VAE(0.0, lhs => 1.0), _DAE(0.0), _KAE(0.0, rhs => -1.0))
function Base.:(*)(lhs::VariableRef, rhs::KnownRef)
    result = zero(DQE)
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# VariableRef--_DecisionAffExpr{C}
Base.:(+)(lhs::VariableRef, rhs::_DecisionAffExpr{C}) where C = DecisionAffExpr{C}(_VariableAffExpr{C}(zero(C), lhs => 1.), copy(rhs), zero(_KnownAffExpr{C}))
Base.:(-)(lhs::VariableRef, rhs::_DecisionAffExpr{C}) where C = DecisionAffExpr{C}(_VariableAffExpr{C}(zero(C), lhs => 1.), -rhs, zero(_KnownAffExpr{C}))
function Base.:(*)(lhs::VariableRef, rhs::_DecisionAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# VariableRef--_DecisionQuadExpr{C}
function Base.:(+)(lhs::VariableRef, rhs::_DecisionQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expressoin(result, rhs)
    return result
end
function Base.:(-)(lhs::VariableRef, rhs::_DecisionQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

# VariableRef--_KnownAffExpr{C}
Base.:(+)(lhs::VariableRef, rhs::_KnownAffExpr{C}) where C = DecisionAffExpr{C}(_VariableAffExpr{C}(zero(C), lhs => 1.), zero(_DecisionAffExpr{C}), copy(rhs))
Base.:(-)(lhs::VariableRef, rhs::_KnownAffExpr{C}) where C = DecisionAffExpr{C}(_VariableAffExpr{C}(zero(C), lhs => 1.), zero(_DecisionAffExpr{C}), -rhs)
function Base.:(*)(lhs::VariableRef, rhs::_KnownAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# VariableRef--_KnownQuadExpr{C}
function Base.:(+)(lhs::VariableRef, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, rhs)
    return result
end
function Base.:(-)(lhs::VariableRef, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

# VariableRef--DecisionAffExpr{C}
Base.:(+)(lhs::VariableRef, rhs::DecisionAffExpr{C}) where C = DecisionAffExpr{C}(lhs + rhs.variables, copy(rhs.decisions), copy(rhs.knowns))
Base.:(-)(lhs::VariableRef, rhs::DecisionAffExpr{C}) where C = DecisionAffExpr{C}(lhs - rhs.variables, -rhs.decisions, -rhs.knowns)
function Base.:(*)(lhs::VariableRef, rhs::DecisionAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end
Base.:/(lhs::VariableRef, rhs::DecisionAffExpr) = error("Cannot divide a variable by an affine expression")

# VariableRef--DecisionQuadExpr{C}
function Base.:(+)(lhs::VariableRef, rhs::DecisionQuadExpr{C}) where C
    result = copy(rhs)
    JuMP.add_to_expression!(result, lhs)
    return result
end
function Base.:(-)(lhs::VariableRef, rhs::DecisionQuadExpr{C}) where C
    result = -rhs
    JuMP.add_to_expression!(result, lhs)
    return result
end
Base.:(*)(lhs::VariableRef, rhs::DecisionQuadExpr) = error("Cannot multiply a quadratic expression by a variable")

#=
    DecisionRef
=#

# DecisionRef
Base.:(-)(lhs::DecisionRef) = DAE(_VAE(0.0), _DAE(0.0, lhs => -1.0), _KAE(0.0))

# DecisionRef--_Constant
Base.:(+)(lhs::DecisionRef, rhs::_Constant) = (+)(rhs, lhs)
Base.:(-)(lhs::DecisionRef, rhs::_Constant) = (+)(-rhs, lhs)
Base.:(*)(lhs::DecisionRef, rhs::_Constant) = (*)(rhs, lhs)
Base.:(/)(lhs::DecisionRef, rhs::_Constant) = (*)(1.0 / rhs, lhs)

# DecisionRef--VariableRef
Base.:(+)(lhs::DecisionRef, rhs::VariableRef) = (+)(rhs, lhs)
Base.:(-)(lhs::DecisionRef, rhs::VariableRef) = (+)(-rhs, lhs)
Base.:(*)(lhs::DecisionRef, rhs::VariableRef) = (*)(rhs, lhs)

# DecisionRef--DecisionRef
Base.:(+)(lhs::DecisionRef, rhs::DecisionRef) = DAE(_VAE(0.0), _DAE(0.0, lhs => 1.0, rhs => +1.0), _KAE(0.0))
Base.:(-)(lhs::DecisionRef, rhs::DecisionRef) = DAE(_VAE(0.0), _DAE(0.0, lhs => 1.0, rhs => -1.0), _KAE(0.0))
function Base.:(*)(lhs::DecisionRef, rhs::DecisionRef)
    result = zero(DQE)
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# DecisionRef--KnownRef
Base.:(+)(lhs::DecisionRef, rhs::KnownRef) = DAE(_VAE(0.0), _DAE(0.0, lhs => 1.0), _KAE(0.0, rhs => +1.0))
Base.:(-)(lhs::DecisionRef, rhs::KnownRef) = DAE(_VAE(0.0), _DAE(0.0, lhs => 1.0), _KAE(0.0, rhs => -1.0))
function Base.:(*)(lhs::DecisionRef, rhs::KnownRef)
    result = zero(DQE)
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# DecisionRef--_VariableAffExpr{C}
Base.:(+)(lhs::DecisionRef, rhs::_VariableAffExpr{C}) where C = DecisionAffExpr{C}(rhs, _DecisionAffExpr{C}(zero(C), lhs => 1.), zero(_KnownAffExpr{C}))
Base.:(-)(lhs::DecisionRef, rhs::_VariableAffExpr{C}) where C = DecisionAffExpr{C}(-rhs, _DecisionAffExpr{C}(zero(C), lhs => 1.), zero(_KnownAffExpr{C}))
function Base.:(*)(lhs::DecisionRef, rhs::_VariableAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# DecisionRef--_VariableQuadExpr{C}
function Base.:(+)(lhs::DecisionRef, rhs::_VariableQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, rhs)
    return result
end
function Base.:(-)(lhs::DecisionRef, rhs::_VariableQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

# DecisionRef--_DecisionAffExpr{C}
Base.:(+)(lhs::DecisionRef, rhs::_DecisionAffExpr{C}) where C = (+)(_DecisionAffExpr{C}(zero(C), lhs => 1.0),  rhs)
Base.:(-)(lhs::DecisionRef, rhs::_DecisionAffExpr{C}) where C = (+)(_DecisionAffExpr{C}(zero(C), lhs => 1.0), -rhs)
function Base.:(*)(lhs::DecisionRef, rhs::_DecisionAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# DecisionRef--_DecisionQuadExpr{C}
function Base.:(+)(lhs::DecisionRef, rhs::_DecisionQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expressoin(result, rhs)
    return result
end
function Base.:(-)(lhs::DecisionRef, rhs::_DecisionQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expressoin(result, lhs)
    JuMP.add_to_expressoin(result, -1., rhs)
    return result
end

# DecisionRef--_KnownAffExpr{C}
Base.:(+)(lhs::DecisionRef, rhs::_KnownAffExpr{C}) where C = DecisionAffExpr{C}(zero(_VariableAffExpr{C}), _DecisionAffExpr{C}(zero(C), lhs => 1), copy(rhs))
Base.:(-)(lhs::DecisionRef, rhs::_KnownAffExpr{C}) where C = DecisionAffExpr{C}(zero(_VariableAffExpr{C}), _DecisionAffExpr{C}(zero(C), lhs => 1.), -rhs)
function Base.:(*)(lhs::DecisionRef, rhs::_KnownAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# DecisionRef--_KnownQuadExpr{C}
function Base.:(+)(lhs::DecisionRef, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, rhs)
    return result
end
function Base.:(-)(lhs::DecisionRef, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

# DecisionRef--DecisionAffExpr{C}
Base.:(+)(lhs::DecisionRef, rhs::DecisionAffExpr{C}) where C = DecisionAffExpr{C}(copy(rhs.variables), lhs+rhs.decisions, copy(rhs.knowns))
Base.:(-)(lhs::DecisionRef, rhs::DecisionAffExpr{C}) where C = DecisionAffExpr{C}(-rhs.variables, lhs-rhs.decisions, -rhs.knowns)
function Base.:(*)(lhs::DecisionRef, rhs::DecisionAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end
Base.:/(lhs::DecisionRef, rhs::DecisionAffExpr) = error("Cannot divide a decision by an affine expression")

# DecisionRef--DecisionQuadExpr{C}
function Base.:(+)(lhs::DecisionRef, rhs::DecisionQuadExpr{C}) where C
    result = copy(rhs)
    JuMP.add_to_expression!(result, lhs)
    return result
end
function Base.:(-)(lhs::DecisionRef, rhs::DecisionQuadExpr{C}) where C
    result = -rhs
    JuMP.add_to_expression!(result, lhs)
    return result
end
Base.:(*)(lhs::DecisionRef, rhs::DecisionQuadExpr) = error("Cannot multiply a quadratic expression by a variable")

#=
    KnownRef
=#

# KnownRef
Base.:(-)(lhs::KnownRef) = DAE(_VAE(0.0), _DAE(0.0), _KAE(0.0, lhs => -1.0))

# KnownRef--_Constant
Base.:(+)(lhs::KnownRef, rhs::_Constant) = (+)(rhs, lhs)
Base.:(-)(lhs::KnownRef, rhs::_Constant) = (+)(-rhs, lhs)
Base.:(*)(lhs::KnownRef, rhs::_Constant) = (*)(rhs, lhs)
Base.:(/)(lhs::KnownRef, rhs::_Constant) = (*)(1.0 / rhs, lhs)

# KnownRef--VariableRef
Base.:(+)(lhs::KnownRef, rhs::VariableRef) = (+)(rhs, lhs)
Base.:(-)(lhs::KnownRef, rhs::VariableRef) = (+)(-rhs, lhs)
Base.:(*)(lhs::KnownRef, rhs::VariableRef) = (*)(rhs, lhs)

# KnownRef--DecisionRef
Base.:(+)(lhs::KnownRef, rhs::DecisionRef) = (+)(rhs, lhs)
Base.:(-)(lhs::KnownRef, rhs::DecisionRef) = (+)(-rhs, lhs)
Base.:(*)(lhs::KnownRef, rhs::DecisionRef) = (*)(rhs, lhs)

# KnownRef--KnownRef
Base.:(+)(lhs::KnownRef, rhs::KnownRef) = DAE(_VAE(0.0), _DAE(0.0), _KAE(0.0, lhs => 1.0, rhs => +1.0))
Base.:(-)(lhs::KnownRef, rhs::KnownRef) = DAE(_VAE(0.0), _DAE(0.0), _KAE(0.0, lhs => 1.0, rhs => -1.0))
function Base.:(*)(lhs::KnownRef, rhs::KnownRef)
    result = zero(DQE)
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# KnownRef--_VariableAffExpr{C}
Base.:(+)(lhs::KnownRef, rhs::_VariableAffExpr{C}) where C = DecisionAffExpr{C}(copy(rhs), zero(_DecisionAffExpr{C}), _KnownAffExpr{C}(zero(C), lhs => 1.))
Base.:(-)(lhs::KnownRef, rhs::_VariableAffExpr{C}) where C = DecisionAffExpr{C}(-rhs, zero(_DecisionAffExpr{C}), _KnownAffExpr{C}(zero(C), lhs => 1.))
function Base.:(*)(lhs::KnownRef, rhs::_VariableAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# KnownRef--_VariableQuadExpr{C}
function Base.:(+)(lhs::KnownRef, rhs::_VariableQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, rhs)
    return result
end
function Base.:(-)(lhs::KnownRef, rhs::_VariableQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

# KnownRef--_DecisionAffExpr{C}
Base.:(+)(lhs::KnownRef, rhs::_DecisionAffExpr{C}) where C = DecisionAffExpr{C}(zero(_VariableAffExpr{C}), rhs, _KnownAffExpr{C}(zero(C), lhs => 1.))
Base.:(-)(lhs::KnownRef, rhs::_DecisionAffExpr{C}) where C = DecisionAffExpr{C}(zero(_VariableAffExpr{C}), -rhs, _KnownAffExpr{C}(zero(C), lhs => 1.))
function Base.:(*)(lhs::KnownRef, rhs::_DecisionAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# KnownRef--_DecisionQuadExpr{C}
function Base.:(+)(lhs::KnownRef, rhs::_DecisionQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expressoin(result, rhs)
    return result
end
function Base.:(-)(lhs::KnownRef, rhs::_DecisionQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expressoin(result, lhs)
    JuMP.add_to_expressoin(result, -1., rhs)
    return result
end

# KnownRef--_KnownAffExpr{C}
Base.:(+)(lhs::KnownRef, rhs::_KnownAffExpr{C}) where C = (+)(_KnownAffExpr{C}(zero(C), lhs => 1.0),  rhs)
Base.:(-)(lhs::KnownRef, rhs::_KnownAffExpr{C}) where C = (+)(_KnownAffExpr{C}(zero(C), lhs => 1.0), -rhs)
function Base.:(*)(lhs::KnownRef, rhs::_KnownAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# KnownRef--_KnownQuadExpr{C}
function Base.:(+)(lhs::KnownRef, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, rhs)
    return result
end
function Base.:(-)(lhs::KnownRef, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

# KnownRef--DecisionAffExpr{C}
Base.:(+)(lhs::KnownRef, rhs::DecisionAffExpr{C}) where C = DecisionAffExpr{C}(copy(rhs.variables), copy(rhs.decisions), lhs+rhs.knowns)
Base.:(-)(lhs::KnownRef, rhs::DecisionAffExpr{C}) where C = DecisionAffExpr{C}(-rhs.variables, -rhs.decisions, lhs-rhs.knowns)
function Base.:(*)(lhs::KnownRef, rhs::DecisionAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end
Base.:/(lhs::KnownRef, rhs::DecisionAffExpr) = error("Cannot divide a known decision by an affine expression")

# KnownRef--DecisionQuadExpr{C}
function Base.:(+)(lhs::KnownRef, rhs::DecisionQuadExpr{C}) where C
    result = copy(rhs)
    JuMP.add_to_expression!(result, lhs)
    return result
end
function Base.:(-)(lhs::KnownRef, rhs::DecisionQuadExpr{C}) where C
    result = -rhs
    JuMP.add_to_expression!(result, lhs)
    return result
end
Base.:(*)(lhs::KnownRef, rhs::DecisionQuadExpr) = error("Cannot multiply a quadratic expression by a variable")

#=
    _VariableAffExpr{C}
=#

# _VariableAffExpr--DecisionRef/KnownRef
Base.:(+)(lhs::_VariableAffExpr, rhs::Union{DecisionRef, KnownRef}) = (+)(rhs, lhs)
Base.:(-)(lhs::_VariableAffExpr, rhs::Union{DecisionRef, KnownRef}) = (+)(-rhs, lhs)
Base.:(*)(lhs::_VariableAffExpr, rhs::Union{DecisionRef, KnownRef}) = (*)(rhs, lhs)

# _VariableAffExpr{C}--_DecisionAffExpr{C}
Base.:(+)(lhs::_VariableAffExpr{C}, rhs::_DecisionAffExpr{C}) where C = DecisionAffExpr{C}(copy(lhs), copy(rhs), zero(_KnownAffExpr{C}))
Base.:(-)(lhs::_VariableAffExpr{C}, rhs::_DecisionAffExpr{C}) where C = DecisionAffExpr{C}(copy(lhs), -rhs, zero(_KnownAffExpr{C}))
function Base.:(*)(lhs::_VariableAffExpr{C}, rhs::_DecisionAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# _VariableAffExpr{C}--_DecisionQuadExpr{C}
function Base.:(+)(lhs::_VariableAffExpr{C}, rhs::_DecisionQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expressoin(result, rhs)
    return result
end
function Base.:(-)(lhs::_VariableAffExpr{C}, rhs::_DecisionQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

# _VariableAffExpr{C}--_KnownAffExpr{C}
Base.:(+)(lhs::_VariableAffExpr{C}, rhs::_KnownAffExpr{C}) where C = DecisionAffExpr{C}(copy(lhs), zero(_DecisionAffExpr{C}), copy(rhs))
Base.:(-)(lhs::_VariableAffExpr{C}, rhs::_KnownAffExpr{C}) where C = DecisionAffExpr{C}(copy(lhs), zero(_DecisionAffExpr{C}), -rhs)
function Base.:(*)(lhs::_VariableAffExpr{C}, rhs::_KnownAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# _VariableAffExpr{C}--_KnownQuadExpr{C}
function Base.:(+)(lhs::_VariableAffExpr{C}, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, rhs)
    return result
end
function Base.:(-)(lhs::_VariableAffExpr{C}, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end
Base.:(*)(lhs::_VariableAffExpr, rhs::_KnownQuadExpr) = error("Cannot multiply a quadratic expression by an affine expression")

# VariableRef{C}--DecisionAffExpr{C}
Base.:(+)(lhs::_VariableAffExpr{C}, rhs::DecisionAffExpr{C}) where C = DecisionAffExpr{C}(lhs+rhs.variables, copy(rhs.decisions), copy(rhs.knowns))
Base.:(-)(lhs::_VariableAffExpr{C}, rhs::DecisionAffExpr{C}) where C = DecisionAffExpr{C}(lhs-rhs.variables, -rhs.decisions, -rhs.knowns)
function Base.:(*)(lhs::_VariableAffExpr{C}, rhs::DecisionAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# _VariableAffExpr{C}--DecisionQuadExpr{C}
function Base.:(+)(lhs::_VariableAffExpr{C}, rhs::DecisionQuadExpr{C}) where C
    result = copy(rhs)
    JuMP.add_to_expression!(result, lhs)
    return result
end
function Base.:(-)(lhs::_VariableAffExpr{C}, rhs::DecisionQuadExpr{C}) where C
    result = -rhs
    JuMP.add_to_expression!(result, lhs)
    return result
end
Base.:(*)(lhs::_VariableAffExpr, rhs::DecisionQuadExpr) = error("Cannot multiply a quadratic expression by an affine expression")

#=
    _DecisionAffExpr{C}
=#

# _DecisionAffExpr--VariableRef/DecisionRef/KnownRef
Base.:(+)(lhs::_DecisionAffExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = (+)(rhs, lhs)
Base.:(-)(lhs::_DecisionAffExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = (+)(-rhs, lhs)
Base.:(*)(lhs::_DecisionAffExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = (*)(rhs, lhs)

# _DecisionAffExpr--_VariableAffExpr
Base.:(+)(lhs::_DecisionAffExpr, rhs::_VariableAffExpr) = (+)(rhs, lhs)
Base.:(-)(lhs::_DecisionAffExpr, rhs::_VariableAffExpr) = (+)(-rhs, lhs)
Base.:(*)(lhs::_DecisionAffExpr, rhs::_VariableAffExpr) = (*)(rhs, lhs)

# _DecisionAffExpr{C}--_VariableQuadExpr{C}
function Base.:(+)(lhs::_DecisionAffExpr{C}, rhs::_VariableQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, rhs)
    return result
end
function Base.:(-)(lhs::_DecisionAffExpr{C}, rhs::_VariableQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

# _DecisionAffExpr--_KnownAffExpr
Base.:(+)(lhs::_DecisionAffExpr, rhs::_KnownAffExpr) = DecisionAffExpr{C}(zero(_VariableAffExpr{C}), copy(lhs), copy(rhs))
Base.:(-)(lhs::_DecisionAffExpr, rhs::_KnownAffExpr) = DecisionAffExpr{C}(zero(_VariableAffExpr{C}), copy(lhs), -rhs)
function Base.:(*)(lhs::_DecisionAffExpr{C}, rhs::_KnownAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# _DecisionAffExpr{C}--_KnownQuadExpr{C}
function Base.:(+)(lhs::_DecisionAffExpr{C}, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, rhs)
    return result
end
function Base.:(-)(lhs::_DecisionAffExpr{C}, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

# _DecisionAffExpr{C}--DecisionAffExpr{C}
Base.:(+)(lhs::_DecisionAffExpr{C}, rhs::DecisionAffExpr{C}) where C = DecisionAffExpr{C}(copy(rhs.variables), lhs+rhs.decisions, copy(rhs.knowns))
Base.:(-)(lhs::_DecisionAffExpr{C}, rhs::DecisionAffExpr{C}) where C = DecisionAffExpr{C}(-rhs.variables, lhs-rhs.decisions, -rhs.knowns)
function Base.:(*)(lhs::_DecisionAffExpr{C}, rhs::DecisionAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# _DecisionAffExpr{C}--DecisionQuadExpr{C}
function Base.:(+)(lhs::_DecisionAffExpr{C}, rhs::DecisionQuadExpr{C}) where C
    result = copy(rhs)
    JuMP.add_to_expression!(result, lhs)
    return result
end
function Base.:(-)(lhs::_DecisionAffExpr{C}, rhs::DecisionQuadExpr{C}) where C
    result = -rhs
    JuMP.add_to_expression!(result, lhs)
    return result
end
Base.:(*)(lhs::_DecisionAffExpr, rhs::DecisionQuadExpr) = error("Cannot multiply a quadratic expression by an affine expression")

#=
    _KnownAffExpr{C}
=#

# _KnownAffExpr--VariableRef/DecisionRef/KnownRef
Base.:(+)(lhs::_KnownAffExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = (+)(rhs, lhs)
Base.:(-)(lhs::_KnownAffExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = (+)(-rhs, lhs)
Base.:(*)(lhs::_KnownAffExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = (*)(rhs, lhs)

# _KnownAffExpr--_VariableAffExpr
Base.:(+)(lhs::_KnownAffExpr, rhs::_VariableAffExpr) = (+)(rhs, lhs)
Base.:(-)(lhs::_KnownAffExpr, rhs::_VariableAffExpr) = (+)(-rhs, lhs)
Base.:(*)(lhs::_KnownAffExpr, rhs::_VariableAffExpr) = (*)(rhs, lhs)

# _KnownAffExpr{C}--_VariableQuadExpr{C}
function Base.:(+)(lhs::_KnownAffExpr{C}, rhs::_VariableQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, rhs)
    return result
end
function Base.:(-)(lhs::_KnownAffExpr{C}, rhs::_VariableQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

# _KnownAffExpr--_DecisionAffExpr
Base.:(+)(lhs::_KnownAffExpr, rhs::_DecisionAffExpr) = (+)(rhs, lhs)
Base.:(-)(lhs::_KnownAffExpr, rhs::_DecisionAffExpr) = (+)(-rhs, lhs)
Base.:(*)(lhs::_KnownAffExpr, rhs::_DecisionAffExpr) = (*)(rhs, lhs)

# _KnownAffExpr{C}--_DecisionQuadExpr{C}
function Base.:(+)(lhs::_KnownAffExpr{C}, rhs::_DecisionQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expressoin(result, rhs)
    return result
end
function Base.:(-)(lhs::_KnownAffExpr{C}, rhs::_DecisionQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expressoin(result, lhs)
    JuMP.add_to_expressoin(result, -1., rhs)
    return result
end

# _KnownAffExpr{C}--DecisionAffExpr{C}
Base.:(+)(lhs::_KnownAffExpr{C}, rhs::DecisionAffExpr{C}) where C = DecisionAffExpr{C}(copy(rhs.variables), copy(rhs.decisions), lhs+rhs.knowns)
Base.:(-)(lhs::_KnownAffExpr{C}, rhs::DecisionAffExpr{C}) where C = DecisionAffExpr{C}(-rhs.variables, -rhs.decisions, lhs-rhs.knowns)
function Base.:(*)(lhs::_KnownAffExpr{C}, rhs::DecisionAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# _KnownAffExpr{C}--DecisionQuadExpr{C}
function Base.:(+)(lhs::_KnownAffExpr{C}, rhs::DecisionQuadExpr{C}) where C
    result = copy(rhs)
    JuMP.add_to_expression!(result, lhs)
    return result
end
function Base.:(-)(lhs::_KnownAffExpr{C}, rhs::DecisionQuadExpr{C}) where C
    result = -rhs
    JuMP.add_to_expression!(result, lhs)
    return result
end

#=
    DecisionAffExpr{C}
=#

Base.:(-)(lhs::DecisionAffExpr{C}) where C = DecisionAffExpr{C}(-lhs.variables, -lhs.decisions, -lhs.knowns)

# DecisionAffExpr--_Constant
Base.:(+)(lhs::DecisionAffExpr, rhs::_Constant) = (+)(rhs, lhs)
Base.:(-)(lhs::DecisionAffExpr, rhs::_Constant) = (+)(-rhs, lhs)
Base.:(*)(lhs::DecisionAffExpr, rhs::_Constant) = (*)(rhs, lhs)
Base.:/(lhs::DecisionAffExpr, rhs::_Constant) = map_coefficients(c -> c/rhs, lhs)
function Base.:^(lhs::DecisionAffExpr{C}, rhs::Integer) where C
    if rhs == 2
        return lhs*lhs
    elseif rhs == 1
        return convert(DecisionQuadExpr{C}, lhs)
    elseif rhs == 0
        return one(GenericQuadExpr{C})
    else
        error("Only exponents of 0, 1, or 2 are currently supported.")
    end
end
Base.:^(lhs::DecisionAffExpr, rhs::_Constant) = error("Only exponents of 0, 1, or 2 are currently supported.")

# DecisionAffExpr--VariableRef/DecisionRef/KnownRef
Base.:(+)(lhs::DecisionAffExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = (+)(rhs, lhs)
Base.:(-)(lhs::DecisionAffExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = (+)(-rhs, lhs)
Base.:(*)(lhs::DecisionAffExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = (*)(rhs, lhs)
Base.:/(lhs::DecisionAffExpr, rhs::VariableRef) = error("Cannot divide affine expression by a variable")
Base.:/(lhs::DecisionAffExpr, rhs::DecisionRef) = error("Cannot divide affine expression by a decision")
Base.:/(lhs::DecisionAffExpr, rhs::KnownRef) = error("Cannot divide affine expression by a known decision s")

# DecisionAffExpr--_VariableAffExpr/_DecisionAffExpr/_KnownAffExpr
Base.:(+)(lhs::DecisionAffExpr, rhs::Union{_VariableAffExpr, _DecisionAffExpr, _KnownAffExpr}) = (+)(rhs, lhs)
Base.:(-)(lhs::DecisionAffExpr, rhs::Union{_VariableAffExpr, _DecisionAffExpr, _KnownAffExpr}) = (+)(-rhs, lhs)
Base.:(*)(lhs::DecisionAffExpr, rhs::Union{_VariableAffExpr, _DecisionAffExpr, _KnownAffExpr}) = (*)(rhs, lhs)

# DecisionAffExpr{C}--DecisionAffExpr{C}
Base.:(+)(lhs::DecisionAffExpr{C}, rhs::DecisionAffExpr{C}) where C =
    DecisionAffExpr{C}(lhs.variables+rhs.variables,
                       lhs.decisions+rhs.decisions,
                       lhs.knowns+rhs.knowns)
Base.:(-)(lhs::DecisionAffExpr{C}, rhs::DecisionAffExpr{C}) where C =
    DecisionAffExpr{C}(lhs.variables-rhs.variables,
                       lhs.decisions-rhs.decisions,
                       lhs.knowns-rhs.knowns)
function Base.:(*)(lhs::DecisionAffExpr{C}, rhs::DecisionAffExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs, rhs)
    return result
end

# DecisionAffExpr{C}--_VariableQuadExpr{C}
function Base.:(+)(lhs::DecisionAffExpr{C}, rhs::_VariableQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, rhs)
    return result
end
function Base.:(-)(lhs::DecisionAffExpr{C}, rhs::_VariableQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

# DecisionAffExpr{C}--_DecisionQuadExpr{C}
function Base.:(+)(lhs::DecisionAffExpr{C}, rhs::_DecisionQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, rhs)
    return result
end
function Base.:(-)(lhs::DecisionAffExpr{C}, rhs::_DecisionQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

# DecisionAffExpr{C}--_KnownQuadExpr{C}
function Base.:(+)(lhs::DecisionAffExpr{C}, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, rhs)
    return result
end
function Base.:(-)(lhs::DecisionAffExpr{C}, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

Base.:(*)(lhs::DecisionAffExpr, rhs::Union{_VariableQuadExpr, _DecisionQuadExpr, _KnownQuadExpr}) = error("Cannot multiply a quadratic expression by an affine expression")

# DecisionAffExpr{C}--DecisionQuadExpr{C}
function Base.:(+)(lhs::DecisionAffExpr{C}, rhs::DecisionQuadExpr{C}) where C
    result = copy(rhs)
    JuMP.add_to_expression!(result, lhs)
    return result
end
function Base.:(-)(lhs::DecisionAffExpr{C}, rhs::DecisionQuadExpr{C}) where C
    result = -rhs
    JuMP.add_to_expression!(result, lhs)
    return result
end
Base.:(*)(lhs::DecisionAffExpr, rhs::DecisionQuadExpr) = error("Cannot multiply a quadratic expression by an affine expression")

#=
    _VariableQuadExpr{C}
=#

# _VariableQuadExpr--DecisionRef/KnownRef
Base.:(+)(lhs::_VariableQuadExpr, rhs::Union{DecisionRef, KnownRef}) = (+)(rhs, lhs)
Base.:(-)(lhs::_VariableQuadExpr, rhs::Union{DecisionRef, KnownRef}) = (+)(-rhs, lhs)

# _VariableQuadExpr--_DecisionAffExpr/_KnownAffExpr/DecisionAffExpr
Base.:(+)(lhs::_VariableQuadExpr, rhs::Union{_DecisionAffExpr, _KnownAffExpr, DecisionAffExpr}) = (+)(rhs, lhs)
Base.:(-)(lhs::_VariableQuadExpr, rhs::Union{_DecisionAffExpr, _KnownAffExpr, DecisionAffExpr}) = (+)(-rhs, lhs)
Base.:(*)(lhs::_VariableQuadExpr, rhs::DecisionAffExpr) = error("Cannot multiply a quadratic expression by an affine expression")

# _VariableQuadExpr{C}--_DecisionQuadExpr{C}
function Base.:(+)(lhs::_VariableQuadExpr{C}, rhs::_DecisionQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expressoin(result, rhs)
    return result
end
function Base.:(-)(lhs::_VariableQuadExpr{C}, rhs::_DecisionQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

# _VariableQuadExpr{C}--_KnownQuadExpr{C}
function Base.:(+)(lhs::_VariableQuadExpr{C}, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, rhs)
    return result
end
function Base.:(-)(lhs::_VariableQuadExpr{C}, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

# _VariableQuadExpr{C}--DecisionQuadExpr{C}
function Base.:(+)(lhs::_VariableQuadExpr{C}, rhs::DecisionQuadExpr{C}) where C
    result = copy(rhs)
    JuMP.add_to_expression!(result, lhs)
    return result
end
function Base.:(-)(lhs::_VariableQuadExpr{C}, rhs::DecisionQuadExpr{C}) where C
    result = -rhs
    JuMP.add_to_expression!(result, lhs)
    return result
end

#=
    _DecisionQuadExpr{C}
=#

# _DecisionQuadExpr--VariableRef/DecisionRef/KnownRef
Base.:(+)(lhs::_DecisionQuadExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = (+)(rhs, lhs)
Base.:(-)(lhs::_DecisionQuadExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = (+)(-rhs, lhs)

# _DecisionQuadExpr--_VariableAffExpr
Base.:(+)(lhs::_DecisionQuadExpr, rhs::Union{_VariableAffExpr, _DecisionAffExpr, _KnownAffExpr, DecisionAffExpr}) = (+)(rhs, lhs)
Base.:(-)(lhs::_DecisionQuadExpr, rhs::Union{_VariableAffExpr, _DecisionAffExpr, _KnownAffExpr, DecisionAffExpr}) = (+)(-rhs, lhs)
Base.:(*)(lhs::_DecisionQuadExpr, rhs::DecisionAffExpr) = error("Cannot multiply a quadratic expression by an affine expression")

# _DecisionQuadExpr{C}--_VariableQuadExpr{C}
Base.:(+)(lhs::_DecisionQuadExpr, rhs::_VariableQuadExpr) = (+)(rhs, lhs)
Base.:(-)(lhs::_DecisionQuadExpr, rhs::_VariableQuadExpr) = (+)(-rhs, lhs)

# _DecisionQuadExpr{C}--_KnownQuadExpr{C}
function Base.:(+)(lhs::_DecisionQuadExpr{C}, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, rhs)
    return result
end
function Base.:(-)(lhs::_DecisionQuadExpr{C}, rhs::_KnownQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, lhs)
    JuMP.add_to_expression!(result, -1., rhs)
    return result
end

# _DecisionQuadExpr{C}--DecisionQuadExpr{C}
function Base.:(+)(lhs::_DecisionQuadExpr{C}, rhs::DecisionQuadExpr{C}) where C
    result = copy(rhs)
    JuMP.add_to_expression!(result, lhs)
    return result
end
function Base.:(-)(lhs::_DecisionQuadExpr{C}, rhs::DecisionQuadExpr{C}) where C
    result = -rhs
    JuMP.add_to_expression!(result, lhs)
    return result
end

#=
    _KnownQuadExpr{C}
=#

# _KnownQuadExpr--VariableRef/DecisionRef/KnownRef
Base.:(+)(lhs::_KnownQuadExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = (+)(rhs, lhs)
Base.:(-)(lhs::_KnownQuadExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = (+)(-rhs, lhs)

# _KnownQuadExpr--_VariableAffExpr/_DecisionAffExpr/_KnownAffExpr/DecisionAffExpr
Base.:(+)(lhs::_KnownQuadExpr, rhs::Union{_VariableAffExpr, _DecisionAffExpr, _KnownAffExpr, DecisionAffExpr}) = (+)(rhs, lhs)
Base.:(-)(lhs::_KnownQuadExpr, rhs::Union{_VariableAffExpr, _DecisionAffExpr, _KnownAffExpr, DecisionAffExpr}) = (+)(-rhs, lhs)
Base.:(*)(lhs::_KnownQuadExpr, rhs::DecisionAffExpr) = error("Cannot multiply a quadratic expression by an affine expression")

# _KnownQuadExpr--_VariableQuadExpr/_DecisionQuadExpr
Base.:(+)(lhs::_KnownQuadExpr, rhs::Union{_VariableQuadExpr, _DecisionQuadExpr}) = (+)(rhs, lhs)
Base.:(-)(lhs::_KnownQuadExpr, rhs::Union{_VariableQuadExpr, _DecisionQuadExpr}) = (+)(-rhs, lhs)

# _KnownQuadExpr{C}--DecisionQuadExpr{C}
Base.:(+)(lhs::_KnownQuadExpr, rhs::DecisionQuadExpr) = (+)(rhs, lhs)
Base.:(-)(lhs::_KnownQuadExpr, rhs::DecisionQuadExpr) = (+)(-rhs, lhs)

#=
    DecisionQuadExpr{C}
=#

function Base.:(-)(lhs::DecisionQuadExpr{C}) where C
    result = zero(DecisionQuadExpr{C})
    JuMP.add_to_expression!(result, -1., lhs.variables)
    JuMP.add_to_expression!(result, -1., lhs.decisions)
    JuMP.add_to_expression!(result, -1., lhs.knowns)
    for (vars, coef) in lhs.cross_terms
        JuMP._add_or_set!(result.cross_terms, vars, -coef)
    end
    for (vars, coef) in lhs.known_variable_terms
        JuMP._add_or_set!(result.known_variable_terms, vars, -coef)
    end
    for (vars, coef) in lhs.known_decision_terms
        JuMP._add_or_set!(result.known_decision_terms, vars, -coef)
    end
    return result
end

# DecisionQuadExpr--_Constant
Base.:(+)(lhs::DecisionQuadExpr, rhs::_Constant) = (+)(rhs, lhs)
Base.:(-)(lhs::DecisionQuadExpr, rhs::_Constant) = (+)(-rhs, lhs)
Base.:(*)(lhs::DecisionQuadExpr, rhs::_Constant) = (*)(rhs, lhs)
Base.:(/)(lhs::DecisionQuadExpr, rhs::_Constant) = (*)(inv(rhs), lhs)

# DecisionQuadExpr--VariableRef/DecisionRef/KnownRef
Base.:(+)(lhs::DecisionQuadExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = (+)(rhs, lhs)
Base.:(-)(lhs::DecisionQuadExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = (+)(-rhs, lhs)
Base.:(/)(lhs::DecisionQuadExpr, rhs::Union{VariableRef, DecisionRef, KnownRef}) = error("Cannot divide a quadratic expression by a variable")

# DecisionQuadExpr--_VariableAffExpr/_DecisionAffExpr/_KnownAffExpr
Base.:(+)(lhs::DecisionQuadExpr, rhs::Union{_VariableAffExpr, _DecisionAffExpr, _KnownAffExpr, DecisionAffExpr}) = (+)(rhs, lhs)
Base.:(-)(lhs::DecisionQuadExpr, rhs::Union{_VariableAffExpr, _DecisionAffExpr, _KnownAffExpr, DecisionAffExpr}) = (+)(-rhs, lhs)
Base.:(*)(lhs::DecisionQuadExpr, rhs::Union{_VariableAffExpr, _DecisionAffExpr, _KnownAffExpr, DecisionAffExpr}) = error("Cannot multiply a quadratic expression by an affine expression")
Base.:(/)(lhs::DecisionQuadExpr, rhs::Union{GenericAffExpr, DecisionAffExpr}) = error("Cannot divide a quadratic expression by an affine expression")

# DecisionQuadExpr{C}--_VariableQuadExpr/_DecisionQuadExpr/_KnownQuadExpr
Base.:(+)(lhs::DecisionQuadExpr, rhs::Union{_VariableQuadExpr, _DecisionQuadExpr, _KnownQuadExpr}) = (+)(rhs, lhs)
Base.:(-)(lhs::DecisionQuadExpr, rhs::Union{_VariableQuadExpr, _DecisionQuadExpr, _KnownQuadExpr}) = (+)(-rhs, lhs)

# DecisionQuadExpr{C}--DecisionQuadExpr{C}
function Base.:(+)(lhs::DecisionQuadExpr{C}, rhs::DecisionQuadExpr{C}) where C
    result = copy(lhs)
    JuMP.add_to_expression!(result, rhs.variables)
    JuMP.add_to_expression!(result, rhs.decisions)
    JuMP.add_to_expression!(result, rhs.knowns)
    add_cross_terms!(result_terms, terms) = begin
        for (vars, coef) in terms
            JuMP._add_or_set!(result_terms, vars, coef)
        end
    end
    add_cross_terms!(result.cross_terms, rhs.cross_terms)
    add_cross_terms!(result.known_variable_terms, rhs.known_variable_terms)
    add_cross_terms!(result.known_decision_terms, rhs.known_decision_terms)
    return result
end
function Base.:(-)(lhs::DecisionQuadExpr{C}, rhs::DecisionQuadExpr{C}) where C
    result = -rhs
    JuMP.add_to_expression!(result, lhs.variables)
    JuMP.add_to_expression!(result, lhs.decisions)
    JuMP.add_to_expression!(result, lhs.knowns)
    add_cross_terms!(result_terms, terms) = begin
        for (vars, coef) in terms
            JuMP._add_or_set!(result_terms, vars, -coef)
        end
    end
    add_cross_terms!(result.cross_terms, lhs.cross_terms)
    add_cross_terms!(result.known_variable_terms, lhs.known_variable_terms)
    add_cross_terms!(result.known_decision_terms, lhs.known_decision_terms)
    return result
end
Base.:(*)(lhs::DecisionQuadExpr, rhs::DecisionQuadExpr) = error("Cannot multiply a quadratic expression by a quadratic expression")


text:
function my_Gale_Shap{T<:Int64}(prop_prefs::AbstractArray{T, 2}, resp_prefs::AbstractArray{T, 2}, caps::Vector{Int})
    m = size(prop_prefs, 2)
    n = size(resp_prefs, 2)
    L = sum(caps)
    prop_matched = zeros(Int64, m)
    resp_matched = zeros(Int64, L)
    prop_pool = collect(1:m)

    indptr = Array(Int, n+1)
    indptr[1] = 1
    for i in 1:n
        indptr[i+1] = indptr[i] + caps[i]
    end

    while length(prop_pool) != 0
        i = pop!(prop_pool)

        for j in prop_prefs[:, i]

            if j == 0
                prop_matched[i] = 0
                break

                elseif j == 1
                count = 0
                k = 1
                l = caps[1]

                for t in k:l
                    if findnext(resp_prefs[:, j], resp_matched[t], 1) > findnext(resp_prefs[:, j], i, 1)
                        if resp_matched[t] == 0
                            resp_matched[t] = i
                            prop_matched[i] = j
                            count = 1
                            break

                        else
                            push!(prop_pool, resp_matched[t])
                            resp_matched[t] = i
                            prop_matched[i] = j
                            count = 1
                            break

                        end
                    else
                        resp_matched[t] = resp_matched[t]
                        prop_matched[i] = prop_matched[i]
                    end

                end
                if count == 1
                    break
                end

            else
                count = 0
                k = sum(caps[1:j-1]) + 1
                l = sum(caps[1:j])
                for t in k:l
                    if findnext(resp_prefs[:, j], resp_matched[t], 1) > findnext(resp_prefs[:, j], i, 1)
                        if resp_matched[t] == 0
                            resp_matched[t] = i
                            prop_matched[i] = j
                            count = 1
                            break

                        else
                            push!(prop_pool, resp_matched[t])
                            resp_matched[t] = i
                            prop_matched[i] = j
                            count = 1
                            break

                        end
                    else
                        resp_matched[t] = resp_matched[t]
                        prop_matched[i] = prop_matched[i]
                    end

                end
                if count == 1
                    break
                end
            end
        end

    end
    return prop_matched, resp_matched, indptr
end

function my_Gale_Shap{T<:Int64}(prop_prefs::AbstractArray{T, 2}, resp_prefs::AbstractArray{T, 2})
    caps = ones(Int, size(resp_prefs, 2))
    prop_matched, resp_matched, indptr = my_Gale_Shap(prop_prefs, resp_prefs, caps)
    return prop_matched, resp_matched
end

max_stars_repo_name: SUZUKITAISHI/matching

max_issues_repo_name: SUZUKITAISHI/matching

max_forks_repo_name: SUZUKITAISHI/matching


text:
#========================== Scoring functions =================================#
struct Score
    val::Int
end
makescore(mg::Int, eg::Int) = Int32(mg + (eg << 16))
scoreMG(s::Integer) = Int(Int32(unsigned(s) & 0x7fff) - Int32(unsigned(s) & 0x8000))
scoreEG(s::Integer) = Int((s + 0x8000) >> 16)

#============================ Piece square tables =============================#

# Accesses the piece square tables, applying any scaling according to the game phase.
function psqt(piece::Piece, sqr::Integer)
    pt = type(piece)
    col = color(piece)
    factor = 1
    if col == BLACK
        sqr = 65 - sqr
        factor = -1
    end
    @inbounds file = FILE_TO_QSIDE_MAP[fileof(sqr)]
    rank = rankof(sqr)
    @inbounds psqt_vec = (PSQT[pt.val][rank][file] + PVALS[pt.val]) * factor
end

#========================= Late Move Reduction tables =========================#

function init_reduction_table()
    lmrtable = zeros(Int, (64, 64))
    for depth in 1:64
        for played in 1:64
            @inbounds lmrtable[depth, played] = floor(Int, 0.6 + log(depth) * log(played) / 2.0)
        end
    end
    lmrtable
end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
using Printf

"""
    prettytime(t)

Convert a floating point value `t` representing an amount of time in seconds to a more
human-friendly formatted string with three decimal places. Depending on the value of `t`
the string will be formatted to show `t` in nanoseconds (ns), microseconds (μs),
milliseconds (ms), seconds (s), minutes (min), hours (hr), or days (day).
"""
function prettytime(t)
    # Modified from: https://github.com/JuliaCI/BenchmarkTools.jl/blob/master/src/trials.jl
    iszero(t) && return "0.000 s"

    if t < 1e-6
        value, units = t * 1e9, "ns"
    elseif t < 1e-3
        value, units = t * 1e6, "μs"
    elseif t < 1
        value, units = t * 1e3, "ms"
    elseif t < minute
        value, units = t, "s"
    elseif t < hour
        value, units = t / minute, "min"
    elseif t < day
        value, units = t / hour, "hr"
    else
        value, units = t / day, "day"
    end

    return @sprintf("%.3f", value) * " " * units
end


text:
x = randn(10000000)
y = randn(10000000)

function inner(x, y)
    result = 0.0
    for i=1:length(x)
        result += x[i] * y[i]
    end
    result
end

# gc_disable()
# tic()
@time for i=1:50 result = sum(x .* y); end
# timing = toc();
# println("Array operations took ", (timing / 50) * 1000, " ms")
# with 1 processor: Array operations took 89.43973046 ms - avg: 90 ms
# with 9 processors: elapsed time: 5.358992903 seconds (4000026400 bytes allocated, 34.95% gc time)

@time for i=1:50 result = inner(x, y) end
# elapsed time: 1.126870943 seconds (800 bytes allocated)

x = randn(1, 10000000)
y = randn(10000000, 1)
#tic()
@time for i=1:50 result = x * y end
# elapsed time: 2.401834172 seconds (4800 bytes allocated)
#timing = toc();
#println("BLAS took ", (timing / 50) * 1000, "ms")


text:
# before running this code:
# start a Julia REPL with: julia -p N with N > 4
# Otherwise you get a BoundsError in the call workers()[1:4]
arr = drand((100,100), workers()[1:4], [1,4])
# 100x100 DArray{Float64,2,Array{Float64,2}}:
#  0.815847   0.354301   0.212216   …  0.649562   0.0508848  0.803529
#  0.214374   0.132324   0.500233      0.147626   0.887912   0.0595262
#  0.396211   0.77957    0.290449      0.856027   0.558378   0.635368
#  0.398797   0.924516   0.473312      0.754448   0.310684   0.214169
#  0.472551   0.291316   0.166627      0.358318   0.711451   0.0217645
#  0.211057   0.0705473  0.0590814  …  0.981098   0.834826   0.289676
#  0.457238   0.666109   0.505208      0.173396   0.589662   0.0842167
#  0.2103     0.793901   0.0772641     0.600755   0.511281   0.482069
#  0.0666974  0.271477   0.47484       0.344734   0.400214   0.529893
#  0.0929895  0.334918   0.317601      0.0641198  0.93372    0.119353
#  ⋮                                ⋱
#  0.480727   0.205939   0.510688   …  0.472838   0.293027   0.7352
#  0.146304   0.431916   0.386177      0.441221   0.750997   0.126209
#  0.454018   0.935533   0.369919      0.483729   0.832018   0.369996
#  0.374881   0.542864   0.642512      0.781415   0.903669   0.669261
#  0.181886   0.602693   0.50041       0.750187   0.553714   0.43828
#  0.699433   0.904601   0.351032   …  0.764433   0.0109924  0.110523
#  0.95344    0.693372   0.25179       0.965437   0.0590245  0.698218
#  0.518332   0.241622   0.476134      0.201053   0.253833   0.882931
#  0.183282   0.134111   0.40739       0.184264   0.104197   0.795514
#  0.539012   0.453164   0.0306549     0.926738   0.326411   0.284542

arr.dims # (100, 100)
arr.pmap # on which workers ? 4-element Array{Int64,1}: 2 3 4 5
arr.indexes # data division:
# 1x4 Array{(UnitRange{Int64},UnitRange{Int64}),2}:
#  (1:100,1:25)  (1:100,26:50)  (1:100,51:75)  (1:100,76:100)
arr.cuts # data division:
# 2-element Array{Array{Int64,1},1}:
#  [1,101]
#  [1,26,51,76,101]
arr.chunks # references on the workers:
# 1x4 Array{RemoteRef,2}:
#  RemoteRef(2,1,11164)  RemoteRef(3,1,11165)  …  RemoteRef(5,1,11167)

da = @parallel [2i for i = 1:10]
# 10-element DArray{Int64,1,Array{Int64,1}}:
#   2
#   4
#   6
#   8
#  10
#  12
#  14
#  16
#  18
#  20

DArray((10,10)) do I
    println(I)
    return zeros(length(I[1]),length(I[2]))
end
#         From worker 2:  (1:5,1:3)
#         From worker 8:  (1:5,9:10)
#         From worker 4:  (1:5,4:5)
#         From worker 3:  (6:10,1:3)
#         From worker 5:  (6:10,4:5)
#         From worker 7:  (6:10,6:8)
#         From worker 6:  (1:5,6:8)
#         From worker 9:  (6:10,9:10)
# 10x10 DArray{Float64,2,Array{Float64,2}}:
#  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
#  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
#  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
#  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
#  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
#  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
#  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
#  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
#  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
#  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

max_stars_repo_path: Module 1/Chapter08/distrib_arrays.jl

max_issues_repo_path: Module 1/Chapter08/distrib_arrays.jl

max_forks_repo_path: Module 1/Chapter08/distrib_arrays.jl


text:
# # Dataframes.jl: a Perfect Sidekick for Your Next Data Science Project
# ## Bogumił Kamiński
# ### Prepared for https://odsc.com/europe/, 10:15-11:45 AM GMT, June 9, 2021

# The tutorial is loosely based on:
# https://cran.r-project.org/web/packages/data.table/vignettes/datatable-intro.html

# ## Intended usage of this tutorial

# Before you start make sure that you have
# [Julia](https://julialang.org/) installed
# and ensure that your system can find the `julia` executable.

# The tutorial was tested under Julia 1.6.1.

# Next execute the following:
# 1. Clone GitHub repository https://github.com/bkamins/ODSC-EUROPE-2021.git containing the tutorial into a local folder
# 2. Execute `julia --project` in the project folder to start Julia REPL
# 3. pesss `]` to enter package manager mode and execute `instantiate`; press backspace
# 4. Follow the code contained in this file

# ## Load the packages we are going to use in this tutorial

using Arrow
using Chain
using CSV
using DataFrames
using Dates
using GLM
using HTTP
using Plots
using Statistics

# ## Basic operations on data frames

# Create a samle data frame

df = DataFrame(id = ["b","b","b","a","a","c"],
               a = 1:6, b = 7:12, c = 13:18)

#-

typeof(df)

# check its dimensions

nrow(df)

#-

ncol(df)

# You can index into a data frame just like you would index a matrix

df[1:2, 2:end] # subset rows and columns

#-

df[1, :] # get one row

#-

df[:, 1] # get one column

#-

df[1, 1] # get a cell

# The extra functionality is that:

# a) you can select columns using their names

df.id # a most common way to do the selection - treat column name as a field

#-

df[:, "id"] # column names can be passed as strings

#-

df[:, :id] # or as Symbols, which are slightly faster

# b) you can use non-standard: `Not`, `Between`, `Cols`, `All`, and `Regex` selectors

df[Not(1), Between("a", "c")]

# The indexing is much more powerful than what we have shown.
# You can find all the rules here: https://dataframes.juliadata.org/stable/lib/indexing/

# In order to get summary statistics of a data frame use the describe function

describe(df)

# ## Reading and writing CSV and Apache Arrow files

# Fetch from the Internet and load a file that we will use in this tutorial

input = "https://raw.githubusercontent.com/Rdatatable" *
        "/data.table/master/vignettes/flights14.csv"

#-

flights = CSV.read(HTTP.get(input).body, DataFrame)

# R users will likely ask if it is possible to use piping to perform
# the steps shown above. Here is how you can do it:

@chain input begin
    HTTP.get
    _.body
    CSV.read(DataFrame)
end

# by default `@chain` passes the value of the previous operation as a first
# argument to the next operation; unless you explicitly use `_` in which
# case you can put it anywhere in the expression.

# Let us inspect how `@chain` rewrites our code:

@macroexpand @chain input begin
    HTTP.get
    _.body
    CSV.read(DataFrame)
end

# The data frame

# Writing a CSV file is also easy

CSV.write("flights14.csv", flights)

# And now save the file in Apache Arrow format

Arrow.write("flights14.arrow", flights)

# Now read back the files

tmp1 = CSV.read("flights14.csv", DataFrame)

#-

tmp2 = Arrow.Table("flights14.arrow") |> DataFrame

# Check if all three tables hold the same data

flights == tmp1 == tmp2

# ## Filtering

# Assume we are only interested in fligts from "EWR" to "PHL"

# We use indexing

flights[(flights.origin .== "EWR") .& (flights.dest .== "PHL"), :]

# or specialized functions

# `filter` works row-wise

filter(row -> row.origin == "EWR" && row.dest == "PHL", flights)

# `subset` takes whole columns

subset(flights, :origin => x -> x .== "EWR", :dest => x -> x .== "PHL")

#-

subset(flights, :origin => ByRow(==("EWR")), :dest => ByRow(==("PHL")))

# you can also set an index on `:origin` and `:dest` columns for fast lookup
# if you expect to do it many times

flights_idx = groupby(flights, [:origin, :dest])

#-

flights_idx[("EWR", "PHL")]

# The detailed comparison of `filter` and `subset` design is discussed in
# the following post: https://bkamins.github.io/julialang/2021/05/07/subset.html.
# Here let me just comment that working on whole columns is useful when you
# want to define filtering conditions based on group aggregates.
# For instance for every `:origin`-`:dest` pair pick the flight that had highest `:air_time`

subset(flights_idx, :air_time => x -> x .== maximum(x))

# Why do we have more rows than `:origin`-`:dest` pairs? This question naturally
# leads us to the next section.

# ## Aggregation using split-apply-combine

combine(flights_idx) do sdf
    max_air_time = maximum(sdf.air_time)
    return count(sdf.air_time .== max_air_time)
end

# Let us keep only the cases when we have more than one entry per group

@chain flights_idx begin
    combine(_) do sdf
        max_air_time = maximum(sdf.air_time)
        return count(sdf.air_time .== max_air_time)
    end
    filter(:x1 => >(1), _)
end

# Simple aggregations are supported quite conveniently.
# For each `:month` pair find number of flights and average `:dep_delay`

@chain flights begin
    groupby(:month)
    combine(nrow, :dep_delay => mean)
end

# DataFrames.jl ensures a reasonably good performance of such aggregations.
# See https://h2oai.github.io/db-benchmark/ for a benchmark

# ## Sorting

# Assume we want to sort the data frame that we have just obtained by `:dep_delay_mean`

@chain flights begin
    groupby(:month)
    combine(nrow, :dep_delay => mean)
    sort(:dep_delay_mean)
end

# You can sort by multiple columns and can specify complex ordering rules, e.g.

sort(df, [:id, order(:a, rev=true)])

# ## Joining

# Assume we have a dictionary table with month names:

months = DataFrame(month=1:10,
                   month_name=["Jan", "Feb", "Mar", "Apr", "May",
                               "Jun", "Jul", "Aug", "Sep", "Oct"])

# We want to add the column `:month_name` to our original table

leftjoin(flights, months, on=:month)

# Let us add it to our earlier analysis

@chain flights begin
    leftjoin(months, on=:month)
    groupby(:month_name)
    combine(nrow, :dep_delay => mean)
    sort(:dep_delay_mean)
end

# All standard joins are supported: `innerjoin`, `leftjoin`, `rightjoin`,
# `outerjoin`, `semijoin`, `antijoin`, `crossjoin`

# For joins DataFrames.jl also ensures a reasonably good performance.
# See https://h2oai.github.io/db-benchmark/ for a benchmark

# ## Reshaping

# Assume we want to see a cross tabultaion of number of flights per carrier
# in consecutive months.
# We already know how to get this information if long format.

@chain flights begin
    groupby([:month, :carrier], sort=true)
    combine(nrow)
end

# If we wanted to get carriers as consecutive columns we can
# reshape the data frame into wide format

@chain flights begin
    groupby([:month, :carrier], sort=true)
    combine(nrow)
    unstack(:month, :carrier, :nrow)
end

# In order to go from wide to long format use the `stack` function.

# ## Transforming

# We already know the combine function which combines rows by aggregating them.
# DataFrames.jl also provides `select` and `transform` functions that always
# retain the number and order of rows when doing transformations

# Assume we wanted to add a `:total_delay` column to our data frame that
# is a sum of `:dep_delay` and `:arr_delay`.

# The simplest way to do it is just

flights.total_delay = flights.dep_delay + flights.arr_delay

#-

flights

# A declarative way to do the same is (note `!` which signals in-place operation)

transform!(flights, [:dep_delay, :arr_delay] => (+) => :total_delay)

# The general syntax of all transformations in DataFrames.jl is
# `source_columns => transformation_function => target_column_names`

# Let us make one more example: convert `:year`, `:month`, and `:day` columns into a date:

select(flights, [:year, :month, :day], [:year, :month, :day] => ByRow(Date))

# Note that `select`, as opposed to `transform`, only keeps the columns that are specified.
# Also observe that DataFrames.jl tries to auto-generate a meaningful target column
# name in case it is not passed by the user.

# ## Integration with the ecosystem: plotting data, building basic predictive models

# Before concluding the tutorial let us build a linear model and plot explaining
# average `:total_delay` by `:month`

# Note that `@aside` inside `@chain` allows us to fork the execution of the pipe.

@chain flights begin
    groupby(:month, sort=true)
    combine(:total_delay => mean => :mean_delay)
    @aside lm(@formula(mean_delay~month), _) |> display
    plot(_.month, _.mean_delay, label=nothing,
         xlabel="month", ylabel="mean delay")
end

# ## Why do I call DataFrames.jl a sidekick?

# Consider the following data frame (example supplied by Nils Gudat):

df_long = DataFrame(id=1:10^7,
                    trials=repeat(1:10, 10^6),
                    heads=repeat(1:2, 5*10^6))

# We want to expand this table to hold `:id` and `:toss` columns
# where `:toss` is a `Bool` variable indicating if head of tail was tossed
# You can do it e.g. like this

function expand(t, h)
    @assert 0 <= h <= t
    x = falses(t)
    x[1:h] .= true
    return x
end

@chain df_long begin
    combine(:id => ByRow(identity) => :id,
            [:trials, :heads] => ByRow(expand) => :toss)
    @aside println(first(_, 5))
    flatten(:toss)
end

# However, it is quite slow

@time @chain df_long begin
    combine(:id => ByRow(identity) => :id,
            [:trials, :heads] => ByRow(expand) => :toss)
    flatten(:toss)
end;

# The beauty of Julia is in that case you can write
# a user defined function that will be fast and fully generic

function fast_expand(id, t, h)
    @assert length(id) == length(t) == length(h)
    @assert all(x -> 0 <= x[1] <= x[2], zip(h, t))
    st = sum(t)
    x = falses(st)
    idx = similar(id, st)
    i = 1
    for (idv, nt, ns) in zip(id, t, h)
       x[i:i+ns-1] .= true
       idx[i:i+nt-1] .= Ref(idv)
       i += nt
    end
    return DataFrame("id" => idx, "toss" => x, copycols=false)
end

# additionally we can learn to use the `eachcol` function

fast_expand(eachcol(df_long)...)

# let us check if the resutls match

fast_expand(eachcol(df_long)...) == @chain df_long begin
    combine(:id => ByRow(identity) => :id,
            [:trials, :heads] => ByRow(expand) => :toss)
    flatten(:toss)
end

# and the timing

@time fast_expand(eachcol(df_long)...);

# ## Concluding remarks

# Before I finish let me mention two packages.

# DataFramesMeta.jl: https://github.com/JuliaData/DataFramesMeta.jl
# is a very nice package that allows one to use R-like references to variable
# names in expressions

# Literate.jl can be used to convert this file into Jupyter Notebook.
# Just run the following commands in your OS prompt to see the result:
# ```
# julia --project -e "using Literate; Literate.notebook(\"flights.jl\", execute=false)"
# julia --project -e "using IJulia; notebook(dir=pwd())"
# ```


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Level-set
# ------------------------------------------------------------
const LVConstraint = CI{AffineDecisionFunction{Float64}, MOI.LessThan{Float64}}
@with_kw mutable struct LVData{T <: AbstractFloat}
    Q̃::T = 1e10
    objective::AffineDecisionFunction{T} = zero(AffineDecisionFunction{T})
    constraint::LVConstraint = LVConstraint(0)
    objective_cached::Bool = false
    incumbent::Int = 1
    major_iterations::Int = 0
    minor_iterations::Int = 0
    levelindex::Int = -1
end

@with_kw mutable struct LVParameters{T <: AbstractFloat}
    λ::T = 0.5
end

"""
    LevelSet

Functor object for using level-set regularization in an L-shaped algorithm. Create by supplying an [`LV`](@ref) object through `regularize` in `LShaped.Optimizer` or by setting the [`Regularizer`](@ref) attribute.

...
# Parameters
- `λ::AbstractFloat = 0.5`: Controls the level position L = (1-λ)*θ + λ*Q̃, a convex combination of the current lower and upper bound.
- `penaltyterm::PenaltyTerm = Quadratic`: Specify penaltyterm variant ([`Quadratic`](@ref), [`InfNorm`](@ref), [`ManhattanNorm`][@ref])
...
"""
struct LevelSet{T <: AbstractFloat, A <: AbstractVector, PT <: AbstractPenaltyTerm} <: AbstractRegularization
    data::LVData{T}
    parameters::LVParameters{T}

    decisions::DecisionMap
    projection_targets::Vector{MOI.VariableIndex}
    ξ::Vector{Decision{T}}

    Q̃_history::A
    levels::A
    incumbents::Vector{Int}

    penaltyterm::PT

    function LevelSet(decisions::DecisionMap, ξ₀::AbstractVector, penaltyterm::AbstractPenaltyTerm; kw...)
        T = promote_type(eltype(ξ₀), Float32)
        A = Vector{T}
        ξ = map(ξ₀) do val
            KnownDecision(val, T)
        end
        PT = typeof(penaltyterm)
        return new{T, A, PT}(LVData{T}(),
                             LVParameters{T}(;kw...),
                             decisions,
                             Vector{MOI.VariableIndex}(undef, length(ξ₀)),
                             ξ,
                             A(),
                             A(),
                             Vector{Int}(),
                             penaltyterm)
    end
end

function initialize_regularization!(lshaped::AbstractLShaped, lv::LevelSet)
    # Add projection targets
    add_projection_targets!(lv, lshaped.master)
    # Sense-correct initial objective
    sense = MOI.get(lshaped.master, MOI.ObjectiveSense())
    coeff = sense == MOI.MIN_SENSE ? 1.0 : -1.0
    lv.data.Q̃ *= coeff
    # Run penalty initialization to check support
    initialize_penaltyterm!(lv.penaltyterm,
                            lshaped.master,
                            1.0,
                            all_decisions(lv.decisions),
                            lv.projection_targets)
    # Delete penalty-term
    remove_penalty!(lv.penaltyterm, lshaped.master)
    return nothing
end

function restore_regularized_master!(lshaped::AbstractLShaped, lv::LevelSet)
    # Delete level-set constraint
    if !iszero(lv.data.constraint.value)
        MOI.delete(lshaped.master, lv.data.constraint)
        lv.data.constraint = LVConstraint(0)
    end
    # Delete penalty-term
    remove_penalty!(lv.penaltyterm, lshaped.master)
    # Delete projection targets
    for var in lv.projection_targets
        MOI.delete(lshaped.master, var)
        StochasticPrograms.remove_decision!(lv.decisions, var)
    end
    empty!(lv.projection_targets)
    return nothing
end

function filter_variables!(lv::LevelSet, list::Vector{MOI.VariableIndex})
    # Filter projection targets
    filter!(vi -> !(vi in lv.projection_targets), list)
    # Filter any auxiliary penaltyterm variables
    remove_penalty_variables!(lv.penaltyterm, list)
    return nothing
end

function filter_constraints!(lv::LevelSet, list::Vector{<:CI})
    # Filter any auxiliary penaltyterm constraints
    remove_penalty_constraints!(lv.penaltyterm, list)
    # Filter level-set constraint
    i = something(findfirst(isequal(lv.data.constraint), list), 0)
    if !iszero(i)
        MOI.deleteat!(list, i)
    end
    return nothing
end

function log_regularization!(lshaped::AbstractLShaped, lv::LevelSet)
    @unpack Q̃,incumbent = lv.data
    push!(lv.Q̃_history, Q̃)
    push!(lv.incumbents, incumbent)
    return nothing
end

function log_regularization!(lshaped::AbstractLShaped, t::Integer, lv::LevelSet)
    @unpack Q̃,incumbent = lv.data
    lv.Q̃_history[t] = Q̃
    lv.incumbents[t] = incumbent
    return nothing
end

function take_step!(lshaped::AbstractLShaped, lv::LevelSet)
    @unpack Q = lshaped.data
    @unpack τ = lshaped.parameters
    @unpack Q̃ = lv.data
    sense = MOI.get(lshaped.master, MOI.ObjectiveSense())
    coeff = sense == MOI.MIN_SENSE ? 1.0 : -1.0
    if coeff*Q <= coeff*Q̃ - τ
        x = current_decision(lshaped)
        for i in eachindex(lv.ξ)
            lv.ξ[i].value = x[i]
        end
        lv.data.Q̃ = Q
        lv.data.incumbent = timestamp(lshaped)
        lv.data.major_iterations += 1
    else
        lv.data.minor_iterations += 1
    end
    if count(active_model_objectives(lshaped)) != num_thetas(lshaped)
        # Only project if all master variables have been added
        return nothing
    end
    if !lv.data.objective_cached
        F = MOI.get(lshaped.master, MOI.ObjectiveFunctionType())
        lv.data.objective = MOI.get(lshaped.master, MOI.ObjectiveFunction{F}())
        lv.data.objective_cached = true
    end
    @unpack Q̃, objective = lv.data
    @unpack λ = lv.parameters
    # Calculate new θ
    # ---------------- #
    # Delete level-set constraint
    if !iszero(lv.data.constraint.value)
        MOI.delete(lshaped.master, lv.data.constraint)
        lv.data.constraint = LVConstraint(0)
    end
    # Delete penalty-term
    remove_penalty!(lv.penaltyterm, lshaped.master)
    # Re-add objective
    F = typeof(objective)
    MOI.set(lshaped.master, MOI.ObjectiveFunction{F}(), objective)
    # Solve unregularized master
    status = solve_master!(lshaped)
    if !(status ∈ AcceptableTermination)
        # Early termination
        return nothing
    end
    # Update model objective
    θ = MOI.get(lshaped.master, MOI.ObjectiveValue())
    lshaped.data.θ = θ
    # Calculate new level
    L = (1-λ)*θ + λ*Q̃
    push!(lv.levels, L)
    # Reformulate projection problem
    F = MOI.ScalarAffineFunction{Float64}
    MOI.set(lshaped.master, MOI.ObjectiveFunction{F}(), zero(F))
    initialize_penaltyterm!(lv.penaltyterm,
                            lshaped.master,
                            1.0,
                            all_decisions(lv.decisions),
                            lv.projection_targets)
    # Add level constraint
    sense = MOI.get(lshaped.master, MOI.ObjectiveSense())
    coeff = sense == MOI.MIN_SENSE ? 1.0 : -1.0
    lv.data.constraint =
        MOI.add_constraint(lshaped.master, coeff*objective,
                           MOI.LessThan(coeff*L))
    return nothing
end

function process_cut!(lshaped::AbstractLShaped, cut::AnyOptimalityCut, lv::LevelSet)
    if lshaped.execution isa AsynchronousExecution
        # Shift level by cut value to prevent infeasibility
        increment = cut(lshaped.x)
        if !iszero(lv.data.constraint.value)
            L = MOI.get(lshaped.master, MOI.ConstraintSet(), lv.data.constraint)
            MOI.set(lshaped.master, MOI.ConstraintSet(), lv.data.constraint, MOIU.shift_constant(L, -increment))
        end
    end
    return nothing
end

function str(::LevelSet)
    return "L-shaped using level sets"
end

# API
# ------------------------------------------------------------
"""
    LV

Factory object for [`LevelSet`](@ref). Pass to `regularize` in `LShaped.Optimizer` or set the [`Regularizer`](@ref) attribute. Equivalent factory calls: `LV`, `WithLV`, `LevelSet`, `WithLevelSets`. See ?LevelSet for parameter descriptions.

"""
mutable struct LV <: AbstractRegularizer
    penaltyterm::AbstractPenaltyTerm
    parameters::LVParameters{Float64}
end
LV(; penaltyterm = Quadratic(), kw...) = LV(penaltyterm, LVParameters(; kw...))
WithLV(; penaltyterm = Quadratic(), kw...) = LV(penaltyterm, LVParameters(; kw...))
LevelSet(; penaltyterm = Quadratic(), kw...) = LV(penaltyterm, LVParameters(; kw...))
WithLevelSets(; penaltyterm = Quadratic(), kw...) = LV(penaltyterm, LVParameters(; kw...))

function MOI.get(lv::LV, ::RegularizationPenaltyTerm)
    return lv.penaltyterm
end

function MOI.set(lv::LV, ::RegularizationPenaltyTerm, penaltyterm::AbstractPenaltyTerm)
    return lv.penaltyterm = penaltyterm
end

function (lv::LV)(decisions::DecisionMap, x::AbstractVector)
    return LevelSet(decisions, x, lv.penaltyterm; type2dict(lv.parameters)...)
end

function str(::LV)
    return "L-shaped using level sets"
end


text: compute(n::Int)::Int = sum(powermod.(1:n, 1:n, 10_000_000_000)) % 10_000_000_000

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
import Polyhedra # needed for some algorithms
using Reachability: constrained_dimensions

include("models/bouncing_ball.jl")
include("models/thermostat.jl")

models = [bouncing_ball, thermostat]

for model in models
    system, options = model()

    # --- reachability algorithms ---

    # default algorithm
    sol = solve(system, options)

    # --- BFFPSV18

    opC = BFFPSV18(:δ=>0.1)

    # concrete discrete-post operator
    sol = solve(system, options, opC, ConcreteDiscretePost())

    # lazy discrete-post operator
    sol = solve(system, options, opC, LazyDiscretePost())

    # overapproximating discrete-post operator
    sol = solve(system, options, opC, ApproximatingDiscretePost())

    # --- BFFPS19

    opC = BFFPS19(:δ=>0.1)

    sol = solve(system, options, opC, DecomposedDiscretePost())

    # --- GLGM06

    opC = GLGM06(:δ=>0.1)

    # concrete discrete-post operator
    if model == bouncing_ball  # currently this crashes for the thermostat
        opD = ConcreteDiscretePost(:check_invariant_intersection => true)
        sol = solve(system, options, opC, opD);
    end

    # lazy discrete-post operator
    opD = LazyDiscretePost(:check_invariant_intersection => true)
    sol = solve(system, options, opC, opD);

    # --- model analysis ---

    # constrained_dimensions
    HS = system.s
    if model == bouncing_ball
        @test constrained_dimensions(HS) == Dict{Int,Vector{Int}}(1=>[1, 2])
    elseif model == thermostat
        @test constrained_dimensions(HS) == Dict{Int,Vector{Int}}(1=>[1], 2=>[1])

        opC = BFFPS19(:δ=>0.1)

        prop_poly = HPolyhedron([HalfSpace([1.], 17.5)])
        property = BadStatesProperty(prop_poly)

        options[:property]= property
        options[:mode]= "check"

        sol = solve(system, options, opC, DecomposedDiscretePost())
        @test sol.satisfied == true
    end
end


text:
"""
VolumeSlices

TODO add function signatures
TODO add descripton

## Attributes
$(ATTRIBUTES)
"""
@recipe(VolumeSlices, x, y, z, volume) do scene
    Attributes(
        colormap = theme(scene, :colormap),
        colorrange = nothing,
        alpha = 0.1,
        contour = Attributes(),
        heatmap = Attributes(),
    )
end

convert_arguments(::Type{<: VolumeSlices}, x, y, z, volume) = convert_arguments(Contour, x, y, z, volume)
function plot!(vs::VolumeSlices)
    @extract vs (x, y, z, volume)
    replace_automatic!(vs, :colorrange) do
        map(extrema, volume)
    end
    keys = (:colormap, :alpha, :colorrange)
    contour!(vs, vs[:contour], x, y, z, volume)
    planes = (:xy, :xz, :yz)
    hattributes = vs[:heatmap]
    sliders = map(zip(planes, (x, y, z))) do plane_r
        plane, r = plane_r
        idx = Node(1)
        vs[plane] = idx
        hmap = heatmap!(vs, hattributes, x, y, zeros(length(x[]), length(y[]))).plots[end]
        on(idx) do i
            transform!(hmap, (plane, r[][i]))
            indices = ntuple(Val(3)) do j
                planes[j] == plane ? i : (:)
            end
            hmap[3][] = view(volume[], indices...)
        end
        idx
    end
    plot!(scene, vs, rest)
end


text:
module SimplePCHIP
#
# Simple PCHIP implementation following Fritsch and Carlson, SIAM J. NUMER. ANAL. 17 (1980) 238-246.
# Derivative calculated in a fashion similar to SciPy's PchipInterpolate
#

export interpolate

struct _pchip
    N :: Int
    xs :: Vector{Float64}
    ys :: Vector{Float64}
    ds :: Vector{Float64}
end

ϕ(t) = 3t^2 - 2t^3
ψ(t) = t^3 - t^2

function _interp(pchip :: _pchip, x :: Number)
    i = _pchip_index(pchip, x)
    x1, x2 = pchip.xs[i:i+1]
    y1, y2 = pchip.ys[i:i+1]
    d1, d2 = pchip.ds[i:i+1]
    h = x2 - x1

    (y1 * ϕ((x2-x)/h)
     + y2 * ϕ((x-x1)/h)
     - d1*h * ψ((x2-x)/h)
     + d2*h * ψ((x-x1)/h))
end

function _pchip_index(pchip :: _pchip, x)
    N = pchip.N
    if N < 200  # Approximate performance cross-over on my old intel i7-3517U
        i = _pchip_index_linear_search(pchip, x)
    else
        i = _pchip_index_bisectional_search(pchip, x)
    end
    if i == N
        # Treat right endpoint as part of rightmost interval
        @assert (x ≈ pchip.xs[N])
        i = N-1
    end
    i
end

function _pchip_index_linear_search(pchip :: _pchip, x)
    xmin = pchip.xs[1]
    @assert (x >= xmin)

    i = 1
    N = pchip.N
    while i < N  &&  x >= pchip.xs[i+1]
        i = i + 1
    end
    i
end

function _pchip_index_bisectional_search(pchip :: _pchip, x)
    N = pchip.N
    imin, imax = 1, N
    xmin = pchip.xs[imin]
    xmax = pchip.xs[imax]
    @assert (x >= xmin && x <= xmax)

    i = imin + div(imax - imin + 1, 2)
    while imin < imax
        if x < pchip.xs[i]
            imax = i - 1
        elseif x >= pchip.xs[i+1]
            imin = i + 1
        else
            break
        end
        i = imin + div(imax - imin + 1, 2)
    end
    i
end

"Similar to how SciPy's PCHIP does it"
function _initial_ds_scipy(xs, ys)
    h(i) = xs[i+1]-xs[i]
    Δ(i) = (ys[i+1]-ys[i]) / h(i)

    N = length(xs)
    ds = similar(xs)
    if N == 2
        ds[:] .= Δ(1)
    else
        Δl = Δ(1)
        hl = h(1)
        for i ∈ 2:N-1
            Δr = Δ(i)
            hr = h(i)
            if sign(Δl) != sign(Δr) || Δl ≈ 0.0 || Δr ≈ 0.0
                ds[i] = 0.0
            else
                wl = 2hl + hr
                wr = hl + 2hr
                ds[i] = (wl + wr) / (wl/Δl + wr/Δr)
            end
            Δl = Δr
            hl = hr
        end
        ds[1] = _edge_derivative(h(1), h(2), Δ(1), Δ(2))
        ds[N] = _edge_derivative(h(N-1), h(N-2), Δ(N-1), Δ(N-2))
    end
    ds
end

function _edge_derivative(h1, h2, Δ1, Δ2)
    d = ((2h1 + h2)*Δ1 - h2*Δ2) / (h1 + h2)
    if sign(d) != sign(Δ1)
        d = 0.0
    elseif sign(Δ1) != sign(Δ2)  &&  abs(d) > abs(3Δ1)
        d = 3Δ1
    end
    d
end

function interpolate(xs, ys)
    xs_ = [x for x ∈ xs]
    ys_ = [y for y ∈ ys]
    _assert_xs_ys(xs_, ys_)
    ds = _initial_ds_scipy(xs_, ys_)
    pchip = _pchip(length(xs_), xs_, ys_, ds)

    x -> _interp(pchip, x)
end

function _assert_xs_ys(xs, ys)
    N = length(xs)
    @assert (N > 1)
    @assert (N == length(ys))
    assert_monotonic_increase(xs)
end

function assert_monotonic_increase(xs)
    foldl((a,b) -> (@assert (a < b); b), xs)
end

end  # module

max_stars_repo_name: genkuroki/SimplePCHIP

max_issues_repo_name: genkuroki/SimplePCHIP

max_forks_repo_name: genkuroki/SimplePCHIP


text:
function versioninfo(io::IO=stdout)
    println("HSA Runtime ($(functional(:hsa) ? "ready" : "MISSING"))")
    if functional(:hsa)
        println("- Path: $libhsaruntime_path")
        println("- Version: $(libhsaruntime_version)")
        #println("- Initialized: $(repr(HSA_REFCOUNT[] > 0))")
    end
    println("ld.lld ($(functional(:lld) ? "ready" : "MISSING"))")
    if functional(:lld)
        println("- Path: $ld_lld_path")
    end
    println("ROCm-Device-Libs ($(functional(:device_libs) ? "ready" : "MISSING"))")
    if functional(:device_libs)
        println("- Path: $device_libs_path")
        # TODO: println("- Version: $(device_libs_version)")
        println("- Downloaded: $(repr(device_libs_downloaded))")
    end
    println("HIP Runtime ($(functional(:hip) ? "ready" : "MISSING"))")
    if functional(:hip)
        println("- Path: $libhip_path")
        # TODO: println("- Version: $(libhip_version)")
    end
    println("rocBLAS ($(functional(:rocblas) ? "ready" : "MISSING"))")
    if functional(:rocblas)
        println("- Path: $(Libdl.dlpath(librocblas))")
    end
    println("rocSOLVER ($(functional(:rocsolver) ? "ready" : "MISSING"))")
    if functional(:rocsolver)
        println("- Path: $(Libdl.dlpath(librocsolver))")
    end
    println("rocALUTION ($(functional(:rocalution) ? "ready" : "MISSING"))")
    if functional(:rocalution)
        println("- Path: $(Libdl.dlpath(librocalution))")
    end
    println("rocSPARSE ($(functional(:rocsparse) ? "ready" : "MISSING"))")
    if functional(:rocsparse)
        println("- Path: $(Libdl.dlpath(librocsparse))")
    end
    println("rocRAND ($(functional(:rocrand) ? "ready" : "MISSING"))")
    if functional(:rocrand)
        println("- Path: $(Libdl.dlpath(librocrand))")
    end
    println("rocFFT ($(functional(:rocfft) ? "ready" : "MISSING"))")
    if functional(:rocfft)
        println("- Path: $(Libdl.dlpath(librocfft))")
    end
    println("MIOpen ($(functional(:miopen) ? "ready" : "MISSING"))")
    if functional(:miopen)
        println("- Path: $(Libdl.dlpath(libmiopen))")
    end

    if functional(:hsa)
        println("HSA Agents ($(length(agents()))):")
        for agent in agents()
            println("- ", repr(agent))
        end
    end
end

"""
    functional() -> Bool

Returns `true` if AMDGPU is nominally functional; "functional" currently means
that HSA, HIP, lld, and device libraries are available (although it does not
imply that usages of these components will be successful).

Packages may use the result of this query to determine whether it is safe to:
- Use AMDGPU to compile code
- Query devices, queues, and other runtime state
- Launch compiled kernels on a device
- Wait on launched kernels to complete
- Utilize external ROCm libraries (rocBLAS et. al)

If the full compilation and launch pipeline is desired, then this query should
be sufficient for most packages and applications. This query combines
sub-queries of multiple components; a failing sub-query will propagate to a
`false` return value. For more fine-grained queries, use `functional(::Symbol)`.

This query should never throw.
"""
function functional()
    return functional(:hsa) &&
           functional(:hip) &&
           functional(:lld) &&
           functional(:device_libs)
end
"""
    functional(component::Symbol) -> Bool

Returns `true` if the ROCm component `component` is configured and expected to
function correctly. Available `component` values are:

- `:hsa`         - Queries ROCR library availability, initialization status, and GPU device availability
- `:hip`         - Queries HIP library availability
- `:lld`         - Queries `ld.lld` tool availability
- `:device_libs` - Queries ROCm device libraries availability
- `:rocblas`     - Queries rocBLAS library availability
- `:rocsolver`   - Queries rocSOLVER library availability
- `:rocalution`  - Queries rocALUTION library availability
- `:rocsparse`   - Queries rocSPARSE library availability
- `:rocrand`     - Queries rocRAND library availability
- `:rocfft`      - Queries rocFFT library availability
- `:miopen`      - Queries MIOpen library availability
- `:all`         - Queries all above components

This query should never throw for valid `component` values.
"""
function functional(component::Symbol)
    if component == :hsa
        return hsa_configured &&
               HSA_REFCOUNT[] > 0 &&
               length(get_agents()) > 0
    elseif component == :hip
        return hip_configured
    elseif component == :lld
        return lld_configured
    elseif component == :device_libs
        return device_libs_configured
    elseif component == :rocblas
        return librocblas !== nothing
    elseif component == :rocsolver
        return librocsolver !== nothing
    elseif component == :rocalution
        return librocalution !== nothing
    elseif component == :rocsparse
        return librocsparse !== nothing
    elseif component == :rocrand
        return librocrand !== nothing
    elseif component == :rocfft
        return librocfft !== nothing
    elseif component == :miopen
        return libmiopen !== nothing
    elseif component == :all
        for component in (:hsa, :hip, :lld, :device_libs, :rocblas, :rocsolver,
                          :rocalution, :rocsparse, :rocrand, :rocfft, :miopen)
            functional(component) || return false
        end
        return true
    else
        throw(ArgumentError("Unknown component $(repr(component))"))
    end
end

function has_rocm_gpu()
    if !functional(:hsa)
        return false
    else
        return length(get_agents(:gpu)) > 0
    end
end

function print_build_diagnostics()
    println("Diagnostics:")
    println("-- deps/build.log")
    println(String(read(joinpath(@__DIR__, "..", "deps", "build.log"))))
    println("-- deps/ext.jl")
    println(String(read(joinpath(@__DIR__, "..", "deps", "ext.jl"))))
    println("-- permissions")
    run(`ls -lah /dev/kfd`)
    run(`ls -lah /dev/dri`)
    run(`id`)
end


text:
using FRGRealTime.MathFuns

function massflow2(k, m2, T, lambda, Nf=4)
    (1 / 6) * (Nf + 2) * lambda * TF.Fb1(k, m2, T)
end

#zero momentum lambda flow
function lambdaflow0(k, m2, T, lambda, Nf=4)
    -(1 / 6) * lambda^2 * (Nf + 8) * TF.Fb2(k, m2, T)
end

function RelambdaflowEu(
    p0,
    q0,
    k,
    m2,
    lambda,
    T,
    RelambdaEu,
    Nf = 4,
)
    -(1 / 3) * (TF.ReFb2Eu(p0, k, m2, T) * RelambdaEu^2) -
    (1 / 3) * (TF.ReFb2Eu(p0, k, m2, T) * RelambdaEu^2) -
    (1 / 6) * lambda^2 * (Nf + 4) * TF.Fb2(k, m2, T)
end

function RelambdaflowEutest(
    p0,
    q0,
    k,
    m2,
    lambda,
    T,
    RelambdaEu,
    Nf = 4,
)
    -(1 / 3) * (ReFb2Eutest(p0, k, m2, T) * RelambdaEu^2)
end

function ReFb2Eutest(p0, k, m2, T)
    -1 / 12 * k^4 / (pi^2 * Eb(k, m2)^2 * (p0^2 + 4 * Eb(k, m2)^2))
end

max_stars_repo_path: Applications/O(N) t-chanel Scheme2 Eu/flow.jl