text:
using DashBootstrapComponents, DashHtmlComponents

row = html_div([
    dbc_row(
        [
            dbc_col(html_div("One of two columns"), width = 4),
            dbc_col(html_div("One of two columns"), width = 4),
        ],
        justify = "start",
    ),
    dbc_row(
        [
            dbc_col(html_div("One of two columns"), width = 4),
            dbc_col(html_div("One of two columns"), width = 4),
        ],
        justify = "center",
    ),
    dbc_row(
        [
            dbc_col(html_div("One of two columns"), width = 4),
            dbc_col(html_div("One of two columns"), width = 4),
        ],
        justify = "end",
    ),
    dbc_row(
        [
            dbc_col(html_div("One of two columns"), width = 4),
            dbc_col(html_div("One of two columns"), width = 4),
        ],
        justify = "between",
    ),
    dbc_row(
        [
            dbc_col(html_div("One of two columns"), width = 4),
            dbc_col(html_div("One of two columns"), width = 4),
        ],
        justify = "around",
    ),
]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
immutable TSymbol
  s::AbstractString

  TSymbol(x) = new(string(x))
end

import Base.==, Base.hash, Base.string

function ==(ts::TSymbol, other::TSymbol)
    ts.s == other.s
end

const hashtsymbol_seed = UInt === UInt64 ? 0x8ee11137cfaae3ce : 0x2f28e1ce

function hash(ts::TSymbol, h::UInt)
    h = hash(hashtsymbol_seed, h)
    h $= hash(ts.s)
    return h
end

function string(ts::TSymbol)
    ts.s
end

ext: jl


text:
# This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license

# PriorityQueue
# -------------

"""
    PriorityQueue(K, V, [ord])

Construct a new [`PriorityQueue`](@ref), with keys of type
`K` and values/priorites of type `V`.
If an order is not given, the priority queue is min-ordered using
the default comparison for `V`.

A `PriorityQueue` acts like a `Dict`, mapping values to their
priorities, with the addition of a `dequeue!` function to remove the
lowest priority element.

```jldoctest
julia> a = PriorityQueue(["a","b","c"],[2,3,1],Base.Order.Forward)
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "c" => 1
  "b" => 3
  "a" => 2
```
"""
@compat type PriorityQueue{K,V,O<:Ordering} <: Associative{K,V}
    # Binary heap of (element, priority) pairs.
    xs::Array{Pair{K,V}, 1}
    o::O

    # Map elements to their index in xs
    index::Dict{K, Int}

    function (::Type{PriorityQueue{K,V,O}}){K,V,O<:Ordering}(o::O)
        new{K,V,O}(Vector{Pair{K,V}}(0), o, Dict{K, Int}())
    end

    (::Type{PriorityQueue{K,V,O}}){K,V,O<:Ordering}() = PriorityQueue{K,V,O}(Forward)

    function (::Type{PriorityQueue{K,V,O}}){K,V,O<:Ordering}(ks::AbstractArray{K},
                                                             vs::AbstractArray{V}, o::O)
        # TODO: maybe deprecate
        if length(ks) != length(vs)
            throw(ArgumentError("key and value arrays must have equal lengths"))
        end
        PriorityQueue{K,V,O}(zip(ks, vs), o)
    end

    function (::Type{PriorityQueue{K,V,O}}){K,V,O<:Ordering}(itr, o::O)
        xs = Vector{Pair{K,V}}(length(itr))
        index = Dict{K, Int}()
        for (i, (k, v)) in enumerate(itr)
            xs[i] = Pair{K,V}(k, v)
            if haskey(index, k)
                throw(ArgumentError("PriorityQueue keys must be unique"))
            end
            index[k] = i
        end
        pq = new{K,V,O}(xs, o, index)

        # heapify
        for i in heapparent(length(pq.xs)):-1:1
            percolate_down!(pq, i)
        end

        pq
    end
end

PriorityQueue(o::Ordering=Forward) = PriorityQueue{Any,Any,typeof(o)}(o)
PriorityQueue{K,V}(::Type{K}, ::Type{V}, o::Ordering=Forward) = PriorityQueue{K,V,typeof(o)}(o)

# TODO: maybe deprecate
PriorityQueue{K,V}(ks::AbstractArray{K}, vs::AbstractArray{V},
                   o::Ordering=Forward) = PriorityQueue{K,V,typeof(o)}(ks, vs, o)

PriorityQueue{K,V}(kvs::Associative{K,V}, o::Ordering=Forward) = PriorityQueue{K,V,typeof(o)}(kvs, o)

PriorityQueue{K,V}(a::AbstractArray{Tuple{K,V}}, o::Ordering=Forward) = PriorityQueue{K,V,typeof(o)}(a, o)

length(pq::PriorityQueue) = length(pq.xs)
isempty(pq::PriorityQueue) = isempty(pq.xs)
haskey(pq::PriorityQueue, key) = haskey(pq.index, key)

"""
    peek(pq)

Return the lowest priority key from a priority queue without removing that
key from the queue.
"""
peek(pq::PriorityQueue) = pq.xs[1]

function percolate_down!(pq::PriorityQueue, i::Integer)
    x = pq.xs[i]
    @inbounds while (l = heapleft(i)) <= length(pq)
        r = heapright(i)
        j = r > length(pq) || lt(pq.o, pq.xs[l].second, pq.xs[r].second) ? l : r
        if lt(pq.o, pq.xs[j].second, x.second)
            pq.index[pq.xs[j].first] = i
            pq.xs[i] = pq.xs[j]
            i = j
        else
            break
        end
    end
    pq.index[x.first] = i
    pq.xs[i] = x
end

function percolate_up!(pq::PriorityQueue, i::Integer)
    x = pq.xs[i]
    @inbounds while i > 1
        j = heapparent(i)
        if lt(pq.o, x.second, pq.xs[j].second)
            pq.index[pq.xs[j].first] = i
            pq.xs[i] = pq.xs[j]
            i = j
        else
            break
        end
    end
    pq.index[x.first] = i
    pq.xs[i] = x
end

# Equivalent to percolate_up! with an element having lower priority than any other
function force_up!(pq::PriorityQueue, i::Integer)
    x = pq.xs[i]
    @inbounds while i > 1
        j = heapparent(i)
        pq.index[pq.xs[j].first] = i
        pq.xs[i] = pq.xs[j]
        i = j
    end
    pq.index[x.first] = i
    pq.xs[i] = x
end

function getindex{K,V}(pq::PriorityQueue{K,V}, key)
    pq.xs[pq.index[key]].second
end

function get{K,V}(pq::PriorityQueue{K,V}, key, deflt)
    i = get(pq.index, key, 0)
    i == 0 ? deflt : pq.xs[i].second
end

# Change the priority of an existing element, or equeue it if it isn't present.
function setindex!{K,V}(pq::PriorityQueue{K, V}, value, key)
    if haskey(pq, key)
        i = pq.index[key]
        oldvalue = pq.xs[i].second
        pq.xs[i] = Pair{K,V}(key, value)
        if lt(pq.o, oldvalue, value)
            percolate_down!(pq, i)
        else
            percolate_up!(pq, i)
        end
    else
        enqueue!(pq, key, value)
    end
    value
end

"""
    enqueue!(pq, k, v)

Insert the a key `k` into a priority queue `pq` with priority `v`.

```jldoctest
julia> a = PriorityQueue(["a","b","c"],[2,3,1],Base.Order.Forward)
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "c" => 1
  "b" => 3
  "a" => 2

julia> enqueue!(a, "d", 4)
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 4 entries:
  "c" => 1
  "b" => 3
  "a" => 2
  "d" => 4
```
"""
function enqueue!{K,V}(pq::PriorityQueue{K,V}, key, value)
    if haskey(pq, key)
        throw(ArgumentError("PriorityQueue keys must be unique"))
    end
    push!(pq.xs, Pair{K,V}(key, value))
    pq.index[key] = length(pq)
    percolate_up!(pq, length(pq))
    pq
end

"""
    dequeue!(pq)

Remove and return the lowest priority key from a priority queue.

```jldoctest
julia> a = PriorityQueue(["a","b","c"],[2,3,1],Base.Order.Forward)
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 3 entries:
  "c" => 1
  "b" => 3
  "a" => 2

julia> dequeue!(a)
"c"

julia> a
PriorityQueue{String,Int64,Base.Order.ForwardOrdering} with 2 entries:
  "b" => 3
  "a" => 2
```
"""
function dequeue!(pq::PriorityQueue)
    x = pq.xs[1]
    y = pop!(pq.xs)
    if !isempty(pq)
        pq.xs[1] = y
        pq.index[y.first] = 1
        percolate_down!(pq, 1)
    end
    delete!(pq.index, x.first)
    x.first
end

function dequeue!(pq::PriorityQueue, key)
    idx = pq.index[key]
    force_up!(pq, idx)
    dequeue!(pq)
    key
end

# Unordered iteration through key value pairs in a PriorityQueue
start(pq::PriorityQueue) = start(pq.index)

done(pq::PriorityQueue, i) = done(pq.index, i)

function next{K,V}(pq::PriorityQueue{K,V}, i)
    (k, idx), i = next(pq.index, i)
    return (pq.xs[idx], i)
end


text:
@testset "Test Constructor" begin
    mdl1 = Linear(10, 2)
    mdl2 = Linear(2, 5)
    mdl3 = Linear(5, 3)
    mdl4 = Linear(3, 4)

    mdl2wrong = Linear(4, 5)
    mdl3wrong = Linear(5, 2)

    mdl = SeriesModel([mdl1, mdl2, mdl3, mdl4])

    @test mdl.ninputs == 10
    @test mdl.noutputs == 4
    @test mdl.nparams == 10*2+2*5+5*3+3*4
    @test mdl.nmodels == 4

    @test_throws DimensionMismatch SeriesModel([mdl1, mdl2wrong, mdl3, mdl4])
    @test_throws DimensionMismatch SeriesModel([mdl1, mdl2, mdl3wrong, mdl4])
end

@testset "Test Constructor with Multiple Series Model" begin
    mdl1 = Linear(10, 5)
    mdl2 = SeriesModel([Identity(5), Bias(5)])
    mdl3 = Linear(5, 3)
    mdl4 = SeriesModel([Identity(3), Bias(3)])

    mdl = SeriesModel([mdl1, mdl2, mdl3, mdl4])

    @test mdl.ninputs == 10
    @test mdl.noutputs == 3
    @test mdl.nparams == 10*5+5+5*3+3
    @test mdl.nmodels == 6
    @test isequal(mdl.mdls, [Linear(10, 5), Identity(5),  Bias(5), Linear(5, 3), Identity(3), Bias(3)])
end

@testset "Test Evaluation" begin
    mdl1 = Linear(2, 2)
    mdl2 = Identity(2)
    mdl3 = Linear(2, 3)
    mdl4 = Identity(3)

    mdl = SeriesModel([mdl1, mdl2, mdl3, mdl4])

    A1 = [1 2;
          3 4]
    A2 = [5 6;
          7 8;
          9 10]

    Θ = [vec(A1); vec(A2)]

    z = z_buffer(mdl)
    dx = dx_buffer(mdl)
    dΘ = dΘ_buffer(mdl)
    x = [2.3, 4]
    evaluate!(mdl, x, Θ, z; dx=dx, dΘ=dΘ)
    @test z ≈ A2*A1*x
    @test dx ≈ A2*A1
    @test dΘ ≈ [11.5 13.8 20.0 24.0 10.3 0.0 0.0 22.9 0.0 0.0;
                 16.1 18.4 28.0 32.0 0.0 10.3 0.0 0.0 22.9 0.0;
                 20.7 23.0 36.0 40.0 0.0 0.0 10.3 0.0 0.0 22.9]

    z = z_buffer(mdl)
    x = [5, 4]
    evaluate!(mdl, x, Θ, z)
    @test z ≈ A2*A1*x

    z = z_buffer(mdl)
    dx = dx_buffer(mdl)
    x = [10.2, 4]
    evaluate!(mdl, x, Θ, z; dx=dx, dΘ=dΘ)
    @test z ≈ A2*A1*x
    @test dx ≈ A2*A1

    z = z_buffer(mdl)
    dΘ = dΘ_buffer(mdl)
    x = [2.3, 4]
    evaluate!(mdl, x, Θ, z; dΘ=dΘ)
    @test z ≈ A2*A1*x
    @test dΘ ≈ [11.5 13.8 20.0 24.0 10.3 0.0 0.0 22.9 0.0 0.0;
                 16.1 18.4 28.0 32.0 0.0 10.3 0.0 0.0 22.9 0.0;
                 20.7 23.0 36.0 40.0 0.0 0.0 10.3 0.0 0.0 22.9]
end

@testset "Test Initializer" begin
    N = 1000
    srand(0)

    mdl1 = Linear(10, 5)
    mdl2 = HyperbolicTangent(5)
    mdl3 = Linear(5, 3)
    mdl4 = Bias(3)

    mdl = SeriesModel([mdl1, mdl2, mdl3, mdl4])

    init1 = CustomRandomInitializer(Distributions.Normal(0, 1))
    init2 = CustomRandomInitializer(Distributions.Normal(0, 2))
    init3 = ZeroInitializer()

    init = ComposedInitializer([init1, init2, init3])

    # First System
    mean_estimated1 = 0
    var_estimated1 = 0
    mean_estimated2 = 0
    var_estimated2 = 0
    mean_estimated3 = 0
    var_estimated3 = 0
    for i=1:N
        params = initial_guess(mdl, init)
        params1 = params[1:50]
        params2 = params[51:65]
        params3 = params[66:end]
        mean_estimated1 += sum(params1)/(N*50)
        var_estimated1 += sum(params1.^2)/(N*50)
        mean_estimated2 += sum(params2)/(N*15)
        var_estimated2 += sum(params2.^2)/(N*15)
        mean_estimated3 += sum(params3)/(N*3)
        var_estimated3 += sum(params3.^2)/(N*3)
    end

    std_estimated1 = sqrt(var_estimated1)
    std_estimated2 = sqrt(var_estimated2)
    std_estimated3 = sqrt(var_estimated3)

    @test std_estimated1 ≈ 1.0 atol=1e-1
    @test mean_estimated1 ≈ 0.0 atol=1e-1

    @test std_estimated2 ≈ 2.0 atol=1e-1
    @test mean_estimated2 ≈ 0.0 atol=1e-1

    @test std_estimated3 ≈ 0.0 atol=1e-16
    @test mean_estimated3 ≈ 0.0 atol=1e-16
end

@testset "Test Default Initializer" begin
    N = 1000
    srand(0)

    mdl1 = Linear(10, 5)
    mdl2 = HyperbolicTangent(5)
    mdl3 = Linear(5, 3)
    mdl4 = Bias(3)

    mdl = SeriesModel([mdl1, mdl2, mdl3, mdl4])

    # First System
    mean_estimated1 = 0
    var_estimated1 = 0
    mean_estimated2 = 0
    var_estimated2 = 0
    mean_estimated3 = 0
    var_estimated3 = 0
    for i=1:N
        params = initial_guess(mdl)
        params1 = params[1:50]
        params2 = params[51:65]
        params3 = params[66:end]
        mean_estimated1 += sum(params1)/(N*50)
        var_estimated1 += sum(params1.^2)/(N*50)
        mean_estimated2 += sum(params2)/(N*15)
        var_estimated2 += sum(params2.^2)/(N*15)
        mean_estimated3 += sum(params3)/(N*3)
        var_estimated3 += sum(params3.^2)/(N*3)
    end

    std_estimated1 = sqrt(var_estimated1)
    std_estimated2 = sqrt(var_estimated2)
    std_estimated3 = sqrt(var_estimated3)

    @test std_estimated1 ≈ 1/sqrt(10) atol=1e-1
    @test mean_estimated1 ≈ 0.0 atol=1e-1

    @test std_estimated2 ≈ 1/sqrt(5) atol=1e-1
    @test mean_estimated2 ≈ 0.0 atol=1e-1

    @test std_estimated3 ≈ 0.0 atol=1e-16
    @test mean_estimated3 ≈ 0.0 atol=1e-16
end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_m4_rho0_69(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_m4_rho0_69!(A_copy)
end

@inline function exp_ps_m4_rho0_69!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=5
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: B2 B3 Bb4 B4 Bb5 B5 T2k7
    # Computing B2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing B3 with operation: mult
    mul!(memslots3,memslots2,memslots1)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=0.001388888888888889
    coeff2=0.0001984126984126984
    coeff3=2.48015873015873e-5
    coeff4=2.7557319223985893e-6
    memslots4 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3
    mul!(memslots4,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots5,memslots3,memslots4)
    # Deallocating Bb4 in slot 4
    # Computing Bb5 = x*I+x*A+x*B2+x*B4
    coeff1=0.16666666666666666
    coeff2=0.041666666666666664
    coeff3=0.008333333333333333
    coeff4=1.0
    # Smart lincomb recycle B4
    memslots5 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots5
    mul!(memslots5,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots4,memslots3,memslots5)
    # Deallocating B3 in slot 3
    # Deallocating Bb5 in slot 5
    # Computing T2k7 = x*I+x*A+x*B2+x*B5
    coeff1=1.0
    coeff2=1.0
    coeff3=0.5
    coeff4=1.0
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots4
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 2
    # Deallocating B5 in slot 4
    return memslots1 # Returning T2k7
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using Base.Test

import Base.Test.@test_skip, Base.Test.@test_broken

# When testing the example solution, all tests must pass, even ones marked as skipped or broken.
# The track user will not be affected by this.
# Overwrite @test_skip, @test_broken with @test
macro test_skip(ex)
    @test eval(current_module(), ex)
end

macro test_broken(ex)
    @test eval(current_module(), ex)
end

for exercise in readdir("exercises")
    # Allow only testing specified execises
    if !isempty(ARGS) && !(exercise in ARGS)
        continue
    end

    exercise_path = joinpath("exercises", exercise)
    isdir(exercise_path) || continue

    # Create temporary directory
    temp_path = mktempdir(".")

    # Copy tests & example to the temporary directory
    cp(joinpath(exercise_path, "example.jl"), joinpath(temp_path, "$exercise.jl"))
    cp(joinpath(exercise_path, "runtests.jl"), joinpath(temp_path, "runtests.jl"))

    try
        # Run the tests
        @testset "$exercise example" begin
            # Run the tests within an anonymous module to prevent definitions from
            # one exercise leaking into another.
            eval(Module(), :(include(joinpath($temp_path, "runtests.jl"))))
        end
    finally
        # Delete the temporary directory
        rm(temp_path, recursive=true)
    end
end


text:
# TODO: switch the inferface from Vector{arms} to Set{arms}? For now, order might still be important... but it should not!
# TODO: have a look at https://github.com/KristofferC/TimerOutputs.jl for timing.
# TODO: switch to the generic solve() interface of Kombinator, THEN USE THE PACKAGE'S approximation_ratio/term!!!!!!
# TODO: refactor the instances, an LP formulation is a property of the problem, not of an algorithm. (Still need an extra abstraction layer, because there may be several good formulations for a given problem.)

module CombinatorialBandits
  using LinearAlgebra
  using SparseArrays

  using DataStructures
  using IterTools
  using Distributions

  import Munkres # Avoid clashing with Hungarian: both packages propose a munkres() function.
  import Hungarian
  using LightGraphs
  using JuMP

  import Base: push!, copy, hash, isequal
  import JuMP: value

  ## Combinatorial algorithms. TODO: put this in another package.
  include("algos/Kombinator.jl")
  using .Kombinator
  import .Kombinator: dimension, value

  ## Nonsmooth optimisation. TODO: put this in another package.
  include("nso/NonsmoothOptim.jl")
  using .NonsmoothOptim

  # General algorithm.
  abstract type Policy end
  abstract type PolicyDetails end

  """
      CombinatorialInstance{T}

  A combinatorial instance, also called an environment in reinforcement learning. An instance is the combination of a
  combinatorial structure (i.e. what the bandit is allowed to make in the environment) and of the statistical
  properties of the arms (i.e. what reward the bandit gets).

  Two properties are mandatory:

  * `n_arms`: the number of arms that are available in the instance (i.e. the dimension of the problem)
  * `optimal_average_reward`: the best reward that can be obtained at any round, in expectation
  """
  abstract type CombinatorialInstance{T} end

  """
      State{T}

  Memorise the (global) evolution of the bandit in the environment. This structure mostly contains information required
  by most bandit algorithms to decide their next step.

  * `round`: the number of rounds (time steps) the bandit has already played
  * `regret`: the total regret the bandit faced (i.e. how much more reward it would have got if it always played the
    optimum solution)
  * `reward`: the total reward the bandit managed to gather
  * `arm_counts`: the number of times each arm (and not each solution, which is a set of arms)  has been played.
    This information is used by most bandit algorithms
  * `arm_reward`: the total reward obtained by this arm
  * `arm_average_reward`: the average reward obtained by this arm. This information is used by most bandit algorithms
  * `policy_extension`: if a policy needs more information than this structure contains, it can store, in a key that
    uniquely identifies the policy, anything it might require
  """
  mutable struct State{T}
    round::Int
    regret::Float64
    reward::Float64
    arm_counts::Dict{T, Int}
    arm_reward::Dict{T, Float64}
    arm_average_reward::Dict{T, Float64}
    policy_extension::Dict{Symbol, Any}
  end

  copy(s::State{T}) where T = State{T}(s.round, s.regret, s.reward, copy(s.arm_counts), copy(s.arm_reward),
                                       copy(s.arm_average_reward), copy(s.policy_extension))

  # Define the trace of the execution throughout the rounds. The structure itself is immutable, things only get appended.
  """
      Trace{T}

  Memorise all the details about the evolution of the bandit in the environment. This object can therefore become
  very heavy for long episodes.

  The information is stored in several vectors, all of them having the same length equal to the number of rounds the
  bandit was used:

  * `states`: a list of states the bandit entered into over the episode
  * `arms`: a list of sets of arms that have been played over the episode
  * `reward`: a list of reward the policy received over the episode
  * `policy_details`: a list of policy-defined details over the episode (typically, computation times)
  * `time_choose_action`: a list of times (in milliseconds) to take a round decision over the episode
  """
  struct Trace{T}
    states::Vector{State{T}}
    arms::Vector{Vector{T}}
    reward::Vector{Vector{Float64}}
    policy_details::Vector{PolicyDetails}
    time_choose_action::Vector{Int}
  end

  """
      push!(trace::Trace{T}, state::State{T}, arms::Vector{T}, reward::Vector{Float64}, policy_details::PolicyDetails,
            time_choose_action::Int) where T

  Appends the arguments to the execution trace of the bandit algorithm. More specifically, `trace`'s data structures
  are updated to also include `state`, `arms`, `reward`, `policy_details`, and `time_choose_action` (expressed in
  milliseconds).

  All of these arguments are copied, *except* `policy_details`. (Indeed, the usual scenario is to keep updating the
  state, the chosen solution (set of arms), and the reward objects, but to build the detail object at each round from
  the ground up.)
  """
  function push!(trace::Trace{T}, state::State{T}, arms::Vector{T}, reward::Vector{Float64},
                 policy_details::PolicyDetails, time_choose_action::Int) where T
    push!(trace.states, copy(state))
    push!(trace.arms, copy(arms))
    push!(trace.reward, copy(reward))
    push!(trace.policy_details, policy_details)
    push!(trace.time_choose_action, time_choose_action)
  end

  """
      all_arm_indices(reward)

  Returns a list of arm indices for the given reward. For instance, for a vector of arms, it returns the list of indices
  in that vector: each index is associated to an arm.

      all_arm_indices(instance::CombinatorialInstance{T}) where T

  Returns a list of arm indices for the given combinatorial instance.
  """
  function all_arm_indices end

  all_arm_indices(reward::Vector{Distribution}) = eachindex(view(reward, [1:s for s in size(reward)]...))
  all_arm_indices(reward::Dict{T, Distribution}) where T = collect(keys(reward))

  function all_arm_indices(reward::Matrix{Distribution})
    reward_indices_cartesian = eachindex(view(reward, [1:s for s in size(reward)]...))
    return [Tuple(i) for i in reward_indices_cartesian]
  end

  """
      initial_state(instance::CombinatorialInstance{Int})

  Returns a new empty `State` object for the given combinatorial instance.
  """
  function initial_state(instance::CombinatorialInstance{T}) where T
    zero_counts = Dict(i => 0 for i in all_arm_indices(instance))
    zero_rewards = Dict(i => 0.0 for i in all_arm_indices(instance))
    return State{T}(0, 0.0, 0.0, zero_counts, zero_rewards, copy(zero_rewards), Dict{Symbol, Any}())
  end

  """
      initial_state(instance::CombinatorialInstance{Int})

  Returns a new empty `Trace` object for the given combinatorial instance.
  """
  function initial_trace(::CombinatorialInstance{T}) where T
    return Trace{T}(State{T}[], Vector{T}[], Vector{Float64}[], PolicyDetails[], Int[])
  end

  """
      dimension(instance::CombinatorialInstance{T}) where T

  Returns the dimension of a solution.
  """
  function dimension(instance::CombinatorialInstance{T}) where T
    return instance.n_arms
  end

  """
      is_feasible(instance::CombinatorialInstance{T}, arms::Vector{T})

  Returns whether the set of arms is a solution that can be played for the given combinatorial instance.
  """
  function is_feasible(::CombinatorialInstance{T}, ::Vector{T}) where T end

  """
      is_partially_acceptable(instance::CombinatorialInstance{T}, arms::Vector{T})

  Returns whether the set of arms is either a solution or the subset of a solution that can be played for the given
  combinatorial instance. In some cases, a subset of a solution is not a solution: the distinction between the two
  states is made by calling `is_feasible`.

  Typically, `is_partially_acceptable` returns `true` for an empty set of arms: even if this is not an acceptable
  solution, adding elements to the empty set may yield a perfectly acceptable solution.
  """
  function is_partially_acceptable(::CombinatorialInstance{T}, ::Vector{T}) where T end

  function all_arm_indices(instance::CombinatorialInstance{T}) where T
    if isa(instance.reward, MultivariateDistribution) # Correlated arms.
      # Nothing as generic as the uncorrelated case is available, due to
      # the fact that only a vector is known, for any kind of arms (unlike
      # the uncorrelated case, where there is a distinction between vectors
      # and matrices of reward distributions).
      return instance.all_arm_indices
    else # Uncorrelated arms.
      return all_arm_indices(instance.reward)
    end
  end

  """
      pull(instance::CombinatorialInstance{T}, arms::Vector{T}) where T

  For the given bandit, plays the given solution (i.e. set of arms). This function returns both the reward and the
  regret this action caused (zero if the action is optimum, greater than zero otherwise).
  """
  function pull(instance::CombinatorialInstance{T}, arms::Vector{T}) where T
    # Draw the rewards for this round. If T is a tuple, the reward distributions
    # are stored in a matrix, hence the splatting.
    arm_indices = all_arm_indices(instance)
    if isa(instance.reward, MultivariateDistribution) # Correlated arms.
      true_rewards_vector = mean(instance.reward)
      true_rewards = Dict{T, Float64}(arm => true_rewards_vector[i] for (i, arm) in enumerate(arm_indices))

      drawn_rewards_vector = rand(instance.reward)
      drawn_rewards = Dict{T, Float64}(arm => true_rewards_vector[i] for (i, arm) in enumerate(arm_indices))
    else # Uncorrelated arms.
      true_rewards = Dict{T, Float64}(i => mean(instance.reward[i...]) for i in arm_indices)
      drawn_rewards = Dict{T, Float64}(i => rand(instance.reward[i...]) for i in arm_indices)
    end

    # Select the information that will be provided back to the bandit policy.
    # Here is implemented the semi-bandit setting.
    true_reward = sum(true_rewards[arm] for arm in arms)
    reward = Float64[drawn_rewards[arm] for arm in arms]

    # Compute the incurred regret from the provided solution.
    incurred_regret = instance.optimal_average_reward - sum(true_reward)

    return reward, incurred_regret
  end

  solve_linear(instance::CombinatorialInstance{T}, rewards::Dict{T, Float64}) where T = solve_linear(instance.solver, rewards)
  solve_budgeted_linear(instance::CombinatorialInstance{T}, rewards::Dict{T, Float64}, weights::Dict{T, Int}, budget::Int) where T =
    solve_budgeted_linear(instance.solver, rewards, weights, budget)
  supports_solve_budgeted_linear(instance::CombinatorialInstance{T}) where T = supports_solve_budgeted_linear(instance.solver)
  solve_all_budgeted_linear(instance::CombinatorialInstance{T}, rewards::Dict{T, Float64}, weights::Dict{T, Int}, max_budget::Int) where T =
    solve_all_budgeted_linear(instance.solver, rewards, weights, max_budget)
  supports_solve_all_budgeted_linear(instance::CombinatorialInstance{T}) where T = supports_solve_budgeted_linear(instance.solver)
  has_lp_formulation(instance::CombinatorialInstance{T}) where T = has_lp_formulation(instance.solver)
  get_lp_formulation(instance::CombinatorialInstance{T}, rewards::Dict{T, Float64}) where T = has_lp_formulation(instance) ?
    get_lp_formulation(instance.solver, rewards) :
    error("The chosen solver uses no LP formulation.")
  approximation_ratio(instance::CombinatorialInstance{T}) where T = approximation_ratio(instance.solver)
  approximation_term(instance::CombinatorialInstance{T}) where T = approximation_term(instance.solver)
  approximation_ratio_budgeted(instance::CombinatorialInstance{T}) where T = approximation_ratio_budgeted(instance.solver)
  approximation_term_budgeted(instance::CombinatorialInstance{T}) where T = approximation_term_budgeted(instance.solver)

  # Interface for policies.
  function choose_action(::CombinatorialInstance{T}, policy::Policy, ::State{T}) where T
    error("Policy not implemented: $policy")
  end

  # Update the state before the new round.
  function update!(state::State{T}, ::CombinatorialInstance{T}, arms::Vector{T}, reward::Vector{Float64}, incurred_regret::Float64) where T
    state.round += 1

    # One reward per arm, i.e. semi-bandit feedback (not bandit feedback, where there would be only one reward for all arms).
    for i in 1:length(arms)
      state.arm_counts[arms[i]] += 1
      state.arm_reward[arms[i]] += reward[i]
      state.arm_average_reward[arms[i]] = state.arm_reward[arms[i]] / state.arm_counts[arms[i]]
      state.reward += reward[i]
    end

    state.regret += incurred_regret
  end

  # Use the bandit for the given number of steps.
  function simulate(instance::CombinatorialInstance{T}, policy::Policy, steps::Int; with_trace::Bool=false) where T
    state = initial_state(instance)
    if with_trace
      trace = initial_trace(instance)
    end

    for i in 1:steps
      t0 = time_ns()
      if with_trace
        arms, run_details = choose_action(instance, policy, state, with_trace=true)
      else
        arms = choose_action(instance, policy, state, with_trace=false)
      end
      t1 = time_ns()

      if length(arms) == 0
        error("No arms have been chosen at round $(i)!")
      end

      reward, incurred_regret = pull(instance, arms)
      update!(state, instance, arms, reward, incurred_regret)

      if with_trace
        push!(trace, state, arms, reward, run_details, round(Int, (t1 - t0) / 1_000_000_000))
      end

      if i % 100 == 0
        println(i)
      end
    end

    if ! with_trace
      return state
    else
      return state, trace
    end
  end

  ## Common estimation procedures.
  function estimate_Δmax(instance::CombinatorialInstance{T}, θ::Dict{T, Float64}) where T # TODO: expensive. Memoisation?
    x_max = solve_linear(copy(instance), θ)
    r_max = sum(θ[i] for i in x_max)
    x_min = solve_linear(copy(instance), Dict(k => -v for (k, v) in θ))
    r_min = (length(x_min) > 0) ? sum(θ[i] for i in x_min) : 0.0
    Δmax = r_max - r_min
    return Δmax
  end

  function estimate_Δmin(instance::CombinatorialInstance{T}, θ::Dict{T, Float64}) where T # TODO: expensive. Memoisation?
    x_max = solve_linear(copy(instance), θ) # TODO: redundant with estimate_Δmax. Memoisation?
    r_max = sum(θ[arm] for arm in x_max)

    Δmin = Inf
    for i in keys(θ)
      new_weights = copy(θ)
      if i in x_max
        new_weights[i] = 0.0 # Remove the element.
      else
        new_weights[i] = dimension(instance) * maximum(values(θ)) # Force the element to be taken.
      end
      new_x = solve_linear(instance, new_weights)
      new_r = sum(θ[i] for i in new_x)

      # Focus on suboptimal solutions: Δmin > 0 even if there are multiple optimal solutions.
      if r_max == new_r
        continue
      end

      # Finally, update Δmin.
      if r_max - new_r < Δmin
        Δmin = r_max - new_r
      end
    end
    return Δmin
  end

  function maximum_solution_length(instance::CombinatorialInstance{T}) where T # Common symbol: m; too common, high risk of clash.
    # TODO: expensive but very very common. Memoisation?
    weights = Dict{T, Float64}(k => 1.0 for k in all_arm_indices(instance))
    return length(solve_linear(instance, weights))
  end

  ## Helpers.
  include("budgeted.jl")

  ## Bandit policies.
  include("policies/thompson.jl")
  include("policies/llr.jl")
  include("policies/cucb.jl")
  include("policies/escb2.jl")
  include("policies/olsucb.jl")
  include("policies/ossb.jl")

  include("policies/escb2_exact.jl")
  include("policies/escb2_greedy.jl")
  include("policies/escb2_budgeted.jl")

  include("policies/ossb_helpers.jl")
  include("policies/ossb_convexcombination.jl")
  include("policies/ossb_exact.jl")
  include("policies/ossb_exactsmart.jl")
  include("policies/ossb_exactnaive.jl")
  include("policies/ossb_subgradient_budgeted.jl")

  include("policies/olsucb_greedy.jl")

  ## Potential problems to solve.
  include("instances/_lp_generic.jl")

  include("instances/perfectbipartitematching.jl")
  include("instances/perfectbipartitematching_algos.jl")
  include("instances/perfectbipartitematching_lp.jl")
  include("instances/perfectbipartitematching_munkres.jl")
  include("instances/perfectbipartitematching_hungarian.jl")
  # Not using LightGraphsMatching, due to BlossomV dependency (hard to get to work on Windows, especially…)

  include("instances/elementarypath.jl")
  include("instances/elementarypath_algos.jl")
  include("instances/elementarypath_lightgraphsdijkstra.jl")
  include("instances/elementarypath_lp.jl")

  include("instances/spanningtree.jl")
  include("instances/spanningtree_algos.jl")
  include("instances/spanningtree_lightgraphsprim.jl")
  include("instances/spanningtree_lp.jl")

  include("instances/mset.jl")
  include("instances/mset_algos.jl")
  include("instances/mset_lp.jl")

  # Export all symbols. Code copied from JuMP.
  const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include]

  for sym in names(@__MODULE__, all=true)
    sym_string = string(sym)
    if sym in _EXCLUDE_SYMBOLS || startswith(sym_string, "_")
        continue
    end
    if !(Base.isidentifier(sym) || (startswith(sym_string, "@") &&
         Base.isidentifier(sym_string[2:end])))
       continue
    end
    @eval export $sym
  end
end


text:
using SurrogateModelOptim
using PlotlyJS
using Statistics
using Distributed

# Discontinuous Forrester function
Forrester(x) = (6*x[1] - 2)^2 * sin(12*x[1] - 4)
Forrester_disc(x) = x[1]>0.15 ? Forrester(x[1])+15 : Forrester(x[1])
search_range=[(0.0,1.0)]

# Optimize the test function
result = smoptimize(Forrester_disc, search_range;
                    options=SurrogateModelOptim.Options(
                    iterations=25,
                    num_interpolants=20, #Preferably even number of added processes
                    num_start_samples=5,
                        ));

function plot_fun_1D(fun_original,fun_estimate,search_range)
    N = 101
    x = range(search_range[1][1], length=N,stop=search_range[1][2])
    y_original = fun_original.(x)
    y_estimate = [median(fun_estimate([x])) for x in x]

    trace1 = scatter(;x=x, y=y_original, mode="lines",name="Original function")
    trace2 = scatter(;x=x, y=y_estimate, mode="lines",name="Estimated function")
    p = plot([trace1,trace2])
end

# Plot the results
display(plot_fun_1D(Forrester_disc,result.sm_interpolant,search_range))


text:
toSMatrix() = SMatrix{0,0,Float64}()
toSMatrix(ss::NTuple{N,Number}...) where {N} = toSMatrix(SVector{N}.(ss)...)
toSMatrix(ss::SVector{N}...) where {N} = hcat(ss...)
toSMatrix(::Type{T}, ss...) where {T} = _toSMatrix(T, toSMatrix(ss...))
_toSMatrix(::Type{T}, s::SMatrix{N,M}) where {N,M,T} = convert(SMatrix{N,M,T}, s)
# Dynamic dispatch
toSMatrix(ss::AbstractVector...) = toSMatrix(Tuple.(ss)...)
toSMatrix(s::AbstractMatrix) = SMatrix{size(s,1), size(s,2)}(s)

toSVector(::Tuple{}) = SVector{0,Float64}()
toSVectors(vs...) = [promote(toSVector.(vs)...)...]
toSVector(v::SVector) = v
toSVector(v::NTuple{N,Number}) where {N} = SVector(v)
toSVector(x::Number) = SVector{1}(x)
toSVector(::Type{T}, v) where {T} = T.(toSVector(v))
toSVector(::Type{T}, ::Tuple{}) where {T} = SVector{0,T}()
# Dynamic dispatch
toSVector(v::AbstractVector) = SVector(Tuple(v))

# ensureSMatrix(f::Function) = f
# ensureSMatrix(m::T) where T<:Number = SMatrix{1,1,T,1}(m)
# ensureSMatrix(m::SMatrix) = m
# ensureSMatrix(m::Array) =
#     throw(ErrorException("Write all model terms using scalars or @SMatrix[matrix]"))

_rdr(r1, r2) = (0.5 * (r1 + r2), r2 - r1)

# zerotuple(::Type{T}, ::Val{L}) where {T,L} = ntuple(_ -> zero(T), Val(L))

function padright!(v::Vector, x, n::Integer)
    n0 = length(v)
    resize!(v, max(n, n0))
    for i in (n0 + 1):n
        v[i] = x
    end
    return v
end

padright(sv::StaticVector{E,T}, x::T, ::Val{E}) where {E,T} = sv
padright(sv::StaticVector{E1,T1}, x::T2, ::Val{E2}) where {E1,T1,E2,T2} =
    (T = promote_type(T1,T2); SVector{E2, T}(ntuple(i -> i > E1 ? x : convert(T, sv[i]), Val(E2))))
padright(sv::StaticVector{E,T}, ::Val{E2}) where {E,T,E2} = padright(sv, zero(T), Val(E2))
padright(sv::StaticVector{E,T}, ::Val{E}) where {E,T} = sv
padright(t::NTuple{N´,<:Any}, x, ::Val{N}) where {N´,N} = ntuple(i -> i > N´ ? x : t[i], Val(N))
padright(t::NTuple{N´,<:Any}, ::Val{N}) where {N´,N} = ntuple(i -> i > N´ ? 0 : t[i], Val(N))

filltuple(x, ::Val{L}) where {L} = ntuple(_ -> x, Val(L))

@inline padtotype(s::SMatrix{E,L}, st::Type{S}) where {E,L,E2,L2,S<:SMatrix{E2,L2}} =
    S(SMatrix{E2,E}(I) * s * SMatrix{L,L2}(I))
@inline padtotype(s::UniformScaling, st::Type{S}) where {S<:SMatrix} = S(s)
@inline padtotype(s::Number, ::Type{S}) where {S<:SMatrix} = S(s*I)
@inline padtotype(s::Number, ::Type{T}) where {T<:Number} = T(s)
@inline padtotype(s::AbstractArray, ::Type{T}) where {T<:Number} = T(first(s))
@inline padtotype(s::UniformScaling, ::Type{T}) where {T<:Number} = T(s.λ)

## Work around BUG: -SVector{0,Int}() isa SVector{0,Union{}}
negative(s::SVector{L,<:Number}) where {L} = -s
negative(s::SVector{0,<:Number}) where {L} = s

empty_sparse(::Type{M}, n, m) where {M} = sparse(Int[], Int[], M[], n, m)

display_as_tuple(v, prefix = "") = isempty(v) ? "()" :
    string("(", prefix, join(v, string(", ", prefix)), ")")

displayvectors(mat::SMatrix{E,L,<:AbstractFloat}; kw...) where {E,L} =
    ntuple(l -> round.(Tuple(mat[:,l]); kw...), Val(L))
displayvectors(mat::SMatrix{E,L,<:Integer}; kw...) where {E,L} =
    ntuple(l -> Tuple(mat[:,l]), Val(L))

# pseudoinverse of s times an integer n, so that it is an integer matrix (for accuracy)
pinvmultiple(s::SMatrix{L,0}) where {L} = (SMatrix{0,0,Int}(), 0)
function pinvmultiple(s::SMatrix{L,L´}) where {L,L´}
    L < L´ && throw(DimensionMismatch("Supercell dimensions $(L´) cannot exceed lattice dimensions $L"))
    qrfact = qr(s)
    n = det(qrfact.R)
    # Cannot check det(s) ≈ 0 because s can be non-square
    abs(n) ≈ 0 && throw(ErrorException("Supercell appears to be singular"))
    pinverse = inv(qrfact.R) * qrfact.Q'
    return round.(Int, n * inv(qrfact.R) * qrfact.Q'), round(Int, n)
end

@inline tuplejoin() = ()
@inline tuplejoin(x) = x
@inline tuplejoin(x, y) = (x..., y...)
@inline tuplejoin(x, y, z...) = (x..., tuplejoin(y, z...)...)

function isgrowing(vs::AbstractVector, i0 = 1)
    i0 > length(vs) && return true
    vprev = vs[i0]
    for i in i0 + 1:length(vs)
        v = vs[i]
        v <= vprev && return false
        vprev = v
    end
    return true
end

function ispositive(ndist)
    result = false
    for i in ndist
        i == 0 || (result = i > 0; break)
    end
    return result
end

isnonnegative(ndist) = iszero(ndist) || ispositive(ndist)

############################################################################################
######## _copy! and _add! #  Revise after #33589 is merged #################################
############################################################################################

_copy!(dest, src) = copy!(dest, src)
_copy!(dst::DenseMatrix{<:Number}, src::SparseMatrixCSC{<:Number}) = _fast_sparse_copy!(dst, src)
_copy!(dst::DenseMatrix{<:SMatrix{N,N}}, src::SparseMatrixCSC{<:SMatrix{N,N}}) where {N} = _fast_sparse_copy!(dst, src)
# _copy!(dst::AbstractMatrix{<:Number}, src::AbstractMatrix{<:SMatrix}) = _flatten_muladd!(dst, src)

_add!(dest, src, α) = _plain_muladd(dest, src, α)
_add!(dst::DenseMatrix{<:Number}, src::SparseMatrixCSC{<:Number}, α = 1) = _fast_sparse_muladd!(dst, src, α)
_add!(dst::DenseMatrix{<:SMatrix{N,N}}, src::SparseMatrixCSC{<:SMatrix{N,N}}, α = I) where {N} = _fast_sparse_muladd!(dst, src, α)
# _add!(dst::AbstractMatrix{<:Number}, src::AbstractMatrix{<:SMatrix}, α = I) = _flatten_muladd!(dst, src, α)

# Using broadcast .+= instead allocates unnecesarily
function _plain_muladd(dst, src, α)
    @boundscheck checkbounds(dst, axes(src)...)
    for i in eachindex(src)
        @inbounds dst[i] += α * src[i]
    end
    return dst
end

function _fast_sparse_copy!(dst::DenseMatrix{T}, src::SparseMatrixCSC) where {T}
    @boundscheck checkbounds(dst, axes(src)...)
    fill!(dst, zero(eltype(src)))
    for col in 1:size(src, 1)
        for p in nzrange(src, col)
            @inbounds dst[rowvals(src)[p], col] = nonzeros(src)[p]
        end
    end
    return dst
end

# Only needed for dense <- sparse (#33589), copy!(sparse, sparse) is fine
function _fast_sparse_muladd!(dst::DenseMatrix{T}, src::SparseMatrixCSC, α = I) where {T}
    @boundscheck checkbounds(dst, axes(src)...)
    for col in 1:size(src, 1)
        for p in nzrange(src, col)
            @inbounds dst[rowvals(src)[p], col] += α * nonzeros(src)[p]
        end
    end
    return dst
end

rclamp(r1::UnitRange, r2::UnitRange) = clamp(minimum(r1), extrema(r2)...):clamp(maximum(r1), extrema(r2)...)

# function _flatten_muladd!(dst::DenseMatrix{T}, src::SparseMatrixCSC{S}, α = zero(T)) where {T<:Number,N,S<:SMatrix{N,N}}
#     checkflattenaxes(dst, src)
#     iszero(α) ? fill!(dst, zero(eltype(src))) : (α != 1 && α != I && rmul!(dst, α))
#     for col in 1:size(src, 1)
#         for p in nzrange(src, col)
#             coffset = CartesianIndex(((rowvals(src)[p] - 1) * N, (col - 1) * N))
#             smatrix = nonzeros(src)[p]
#             for i in CartesianIndices((1:N, 1:N))
#                 @inbounds dst[coffset + i] += smatrix[i]
#             end
#         end
#     end
#     return dst
# end

# function _flatten_muladd!(dst::SparseMatrixCSC{T}, src::SparseMatrixCSC{S}, α = zero(T)) where {T<:Number,N,S<:SMatrix{N,N}}
#     rdst = rowvals(dst)
#     ndst = nonzeros(dst)
#     cdst = getcolptr(dst)
#     cdst[1] = 1

#     iszero(α) ? fill!(ndst, zero(eltype(S))) : (α != 1 && α != I && rmul!(ndst, α))
#     copy_structure!(dst, src)

#     p´ = col´ = 0
#     for col in 1:size(src, 2), j in 1:N
#         col´ += 1
#         for p in nzrange(src, col)
#             nz = nonzeros(src)[p]
#             row´ = (rowvals(src)[p] - 1) * N
#             for i in 1:N
#                 row´ += 1
#                 p´ += 1
#                 rdst[p´] = row´
#                 ndst[p´] += nz[i, j]
#             end
#         end
#         cdst[col´ + 1] = p´ + 1
#     end
#     return dst
# end

# function _flatten_copy!(dst::SparseMatrixCSC{T}, src::SparseMatrixCSC{S}) where {T<:Number,N,S<:SMatrix{N,N}}
#     checkflattenaxes(dst, src)
#     l = length(nonzeros(src))
#     l´ = N * N * l
#     rdst = resize!(rowvals(dst), l´)
#     ndst = resize!(nonzeros(dst), l´)
#     cdst = resize!(getcolptr(dst), size(src, 2) * N + 1)
#     cdst[1] = 1

#     fill!(ndst, zero(eltype(S)))

#     p´ = col´ = 0
#     for col in 1:size(src, 2), j in 1:N
#         col´ += 1
#         for p in nzrange(src, col)
#             row´ = (rowvals(src)[p] - 1) * N
#             nz = nonzeros(src)[p]
#             for i in 1:N
#                 row´ += 1
#                 p´ += 1
#                 rdst[p´] = row´
#                 ndst[p´] = nz[i, j]
#             end
#         end
#         cdst[col´ + 1] = p´ + 1
#     end
#     return dst
# end

# function _flatten_muladd!(dst::DenseMatrix{<:Number}, src::DenseMatrix{S}, α = zero(T)) where {T<:Number,N,S<:SMatrix{N,N}}
#     checkflattenaxes(dst, src)
#     iszero(α) ? fill!(dst, zero(eltype(S))) : (α != 1 && α != I && rmul!(dst, α))
#     c = CartesianIndices(src)
#     i0 = first(c)
#     for i in c
#         smatrix = src[i]
#         ioffset = (i - i0) * N
#         for j in 1:N, i in 1:N
#             @inbounds dst[ioffset + CartesianIndex(i, j)] += smatrix[i, j]
#         end
#     end
#     return dst
# end

# function checkflattenaxes(dst::AbstractMatrix{<:Number}, src::AbstractMatrix{S}) where {T,N,S<:SMatrix{N,N,T}}
#     Base.require_one_based_indexing(src)
#     axes(dst) == (1:(N * size(src, 1)), 1:(N * size(src, 2))) ||
#         throw(ArgumentError( "arrays must have the same axes (after flattening) for copy!"))
# end

# checkflattenaxes(dst::AbstractMatrix{<:SMatrix}, src::AbstractMatrix{<:Number}) =
#     throw(ArgumentError("unflattening not supported"))

# function checkflattenaxes(dst::AbstractMatrix, src::AbstractMatrix)
#     axes(dst) == axes(src) ||
#         throw(ArgumentError( "arrays must have the same axes for copy!"))
# end

# function copy_structure!(dst::SparseMatrixCSC{T}, src::SparseMatrixCSC) where {T}
#     @boundscheck checkflattenaxes(dst, src)
#     if !samestructure(dst, src)
#         N = blocksize(eltype(src))
#         rdst = rowvals(dst)
#         ndst = nonzeros(dst)
#         for col in 1:size(src, 2), p in nzrange(src, col)
#             row´ = (rowvals(src)[p] - 1) * N + 1
#             col´ = (col - 1) * N + 1
#             isstored(dst, row´, col´) || _sparse_insert!(dst, row´, col´, N)
#         end
#     end
#     return dst
# end

# samestructure(a, b) =
#     blocksize(eltype(a)) == blocksize(eltype(b)) && rowvals(a) == rowvals(b) && getcolptr(a) == getcolptr(b)

# function _sparse_insert!(s::SparseMatrixCSC{T}, row, col, N) where {T}
#     if N != 1         # Block insert
#         v = view(s, (row):(row + N - 1), (col):(col + N - 1))
#         v .= one(T)   # necessary to create new stored entries
#         v .= zero(T)  # here they are not removed, see #17404
#     else
#         s[row, col] = one(T)
#         s[row, col] = zero(T)
#     end
#     return s
# end

# # Taken fron Julialang/#33821
# function isstored(A::SparseMatrixCSC, i::Integer, j::Integer)
#     @boundscheck checkbounds(A, i, j)
#     rows = rowvals(A)
#     for istored in nzrange(A, j) # could do binary search if the row indices are sorted?
#         i == rows[istored] && return true
#     end
#     return false
# end

# blocksize(s::Type{<:SMatrix{N,N}}) where {N} = N
# blocksize(s::Type{<:Number}) = 1

############################################################################################

function pushapproxruns!(runs::AbstractVector{<:UnitRange}, list::AbstractVector{T},
                         offset = 0, degtol = sqrt(eps(real(T)))) where {T}
    len = length(list)
    len < 2 && return runs
    rmin = rmax = 1
    prev = list[1]
    @inbounds for i in 2:len
        next = list[i]
        if abs(next - prev) < degtol
            rmax = i
        else
            rmin < rmax && push!(runs, (offset + rmin):(offset + rmax))
            rmin = rmax = i
        end
        prev = next
    end
    rmin < rmax && push!(runs, (offset + rmin):(offset + rmax))
    return runs
end

function hasapproxruns(list::AbstractVector{T}, degtol = sqrt(eps(real(T)))) where {T}
    for i in 2:length(list)
        abs(list[i] - list[i-1]) < degtol && return true
    end
    return false
end

eltypevec(::AbstractMatrix{T}) where {T<:Number} = T
eltypevec(::AbstractMatrix{S}) where {N,T<:Number,S<:SMatrix{N,N,T}} = SVector{N,T}

# pinverse(s::SMatrix) = (qrfact = qr(s); return inv(qrfact.R) * qrfact.Q')

# padrightbottom(m::Matrix{T}, im, jm) where {T} = padrightbottom(m, zero(T), im, jm)

# function padrightbottom(m::Matrix{T}, zeroT::T, im, jm) where T
#     i0, j0 = size(m)
#     [i <= i0 && j<= j0 ? m[i,j] : zeroT for i in 1:im, j in 1:jm]
# end

tuplesort((a,b)::Tuple{<:Number,<:Number}) = a > b ? (b, a) : (a, b)
tuplesort(t::Tuple) = t
tuplesort(::Missing) = missing

# collectfirst(s::T, ss...) where {T} = _collectfirst((s,), ss...)
# _collectfirst(ts::NTuple{N,T}, s::T, ss...) where {N,T} = _collectfirst((ts..., s), ss...)
# _collectfirst(ts::Tuple, ss...) = (ts, ss)
# _collectfirst(ts::NTuple{N,System}, s::System, ss...) where {N} = _collectfirst((ts..., s), ss...)
# collectfirsttolast(ss...) = tuplejoin(reverse(collectfirst(ss...))...)

# allorderedpairs(v) = [(i, j) for i in v, j in v if i >= j]

# Like copyto! but with potentially different tensor orders (adapted from Base.copyto!)
function copyslice!(dest::AbstractArray{T1,N1}, Rdest::CartesianIndices{N1},
                    src::AbstractArray{T2,N2}, Rsrc::CartesianIndices{N2}, by = identity) where {T1,T2,N1,N2}
    isempty(Rdest) && return dest
    if length(Rdest) != length(Rsrc)
        throw(ArgumentError("source and destination must have same length (got $(length(Rsrc)) and $(length(Rdest)))"))
    end
    checkbounds(dest, first(Rdest))
    checkbounds(dest, last(Rdest))
    checkbounds(src, first(Rsrc))
    checkbounds(src, last(Rsrc))
    src′ = Base.unalias(dest, src)
    @inbounds for (Is, Id) in zip(Rsrc, Rdest)
        @inbounds dest[Id] = by(src′[Is])
    end
    return dest
end

function appendslice!(dest::AbstractArray, src::AbstractArray{T,N}, Rsrc::CartesianIndices{N}) where {T,N}
    checkbounds(src, first(Rsrc))
    checkbounds(src, last(Rsrc))
    Rdest = (length(dest) + 1):(length(dest) + length(Rsrc))
    resize!(dest, last(Rdest))
    src′ = Base.unalias(dest, src)
    @inbounds for (Is, Id) in zip(Rsrc, Rdest)
        @inbounds dest[Id] = src′[Is]
    end
    return dest
end

######################################################################
# Permutations (taken from Combinatorics.jl)
######################################################################

struct Permutations{T}
    a::T
    t::Int
end

Base.eltype(::Type{Permutations{T}}) where {T} = Vector{eltype(T)}

Base.length(p::Permutations) = (0 <= p.t <= length(p.a)) ? factorial(length(p.a), length(p.a)-p.t) : 0

"""
    permutations(a)
Generate all permutations of an indexable object `a` in lexicographic order. Because the number of permutations
can be very large, this function returns an iterator object.
Use `collect(permutations(a))` to get an array of all permutations.
"""
permutations(a) = Permutations(a, length(a))

"""
    permutations(a, t)
Generate all size `t` permutations of an indexable object `a`.
"""
function permutations(a, t::Integer)
    if t < 0
        t = length(a) + 1
    end
    Permutations(a, t)
end

function Base.iterate(p::Permutations, s = collect(1:length(p.a)))
    (!isempty(s) && max(s[1], p.t) > length(p.a) || (isempty(s) && p.t > 0)) && return
    nextpermutation(p.a, p.t ,s)
end

function nextpermutation(m, t, state)
    perm = [m[state[i]] for i in 1:t]
    n = length(state)
    if t <= 0
        return(perm, [n+1])
    end
    s = copy(state)
    if t < n
        j = t + 1
        while j <= n &&  s[t] >= s[j]; j+=1; end
    end
    if t < n && j <= n
        s[t], s[j] = s[j], s[t]
    else
        if t < n
            reverse!(s, t+1)
        end
        i = t - 1
        while i>=1 && s[i] >= s[i+1]; i -= 1; end
        if i > 0
            j = n
            while j>i && s[i] >= s[j]; j -= 1; end
            s[i], s[j] = s[j], s[i]
            reverse!(s, i+1)
        else
            s[1] = n+1
        end
    end
    return (perm, s)
end

# Taken from Combinatorics.jl
# TODO: This should really live in Base, otherwise it's type piracy
"""
    factorial(n, k)

Compute ``n!/k!``.
"""
function Base.factorial(n::T, k::T) where T<:Integer
    if k < 0 || n < 0 || k > n
        throw(DomainError((n, k), "n and k must be nonnegative with k ≤ n"))
    end
    f = one(T)
    while n > k
        f = Base.checked_mul(f, n)
        n -= 1
    end
    return f
end

Base.factorial(n::Integer, k::Integer) = factorial(promote(n, k)...)