text:
import Base: Forward, ForwardOrdering, Reverse, ReverseOrdering, Lexicographic, LexicographicOrdering, sortperm, Ordering,
            setindex!, getindex, similar, Algorithm
import SortingAlgorithms: RadixSort, RadixSortAlg

struct StringRadixSortAlg <: Algorithm end
const StringRadixSort = StringRadixSortAlg() # this is needed to make sortperm work

import Base: Forward, ForwardOrdering, Reverse, ReverseOrdering, sortperm, Ordering,
            setindex!, getindex, similar

# Radix sort for strings
function sort!(svec::AbstractVector{String}, ::StringRadixSortAlg, o::Perm)
    sort!(svec, 1, length(svec), StringRadixSort, o)
end

function sort!(svec::AbstractVector, lo::Int, hi::Int, ::StringRadixSortAlg, o::O) where O <: Union{ForwardOrdering, ReverseOrdering, Perm}
    if isa(o, Perm)
        if eltype(o.data) != String
            throw(ArgumentError("Cannot use StringRadixSort on type $(eltype(o.data))"))
        end
        return sortperm_radixsort(o.data, order = o.order)
    else
        sort!(svec, lo, hi, StringRadixSort, o)
    end
end

function sort!(svec::AbstractVector{String}, lo::Int, hi::Int, ::StringRadixSortAlg, o::O) where O <: Union{ForwardOrdering, ReverseOrdering}
    if lo >= hi;  return svec;  end
    # the length subarray to sort
    l = hi - lo + 1

    # find the maximum string length
    lens = maximum(sizeof, svec)
    skipbytes = lens
    if lens > 0
        # while lens > 8
        #     skipbytes = max(0, skipbytes - 16)
        #     bits128 = zeros(UInt128, l)
        #     if o == Reverse
        #         bits128[lo:hi] .= .~load_bits.(UInt128, @view(svec[lo:hi]), skipbytes)
        #         sorttwo!(bits128, svec, lo, hi)
        #     else
        #         bits128[lo:hi] .= load_bits.(UInt128, @view(svec[lo:hi]), skipbytes)
        #         sorttwo!(bits128, svec, lo, hi)
        #     end
        #     lens -= 16
        # end
        while lens > 4
            skipbytes = max(0, skipbytes - 8)
            bits64 = zeros(UInt64, l)
            if o == Reverse
                bits64[lo:hi] .= .~load_bits.(UInt64, @view(svec[lo:hi]), skipbytes)
                sorttwo!(bits64, svec, lo, hi)
            else
                bits64[lo:hi] .= load_bits.(UInt64, @view(svec[lo:hi]), skipbytes)
                sorttwo!(bits64, svec, lo, hi)
            end
            lens -= 8
        end
        if lens > 0
            skipbytes = max(0, skipbytes - 4)
            bits32 = zeros(UInt32, l)
            if o == Reverse
                bits32[lo:hi] .= .~load_bits.(UInt32, @view(svec[lo:hi]), skipbytes)
                sorttwo!(bits32, svec, lo, hi)
            else
                bits32[lo:hi] .= load_bits.(UInt32, @view(svec[lo:hi]), skipbytes)
                sorttwo!(bits32, svec, lo, hi)
            end
            lens -= 4
        end
    end
    svec
end

"""
    sorttwo!(vs, index)

Sort both the `vs` and reorder `index` at the same. This allows for faster sortperm
for radix sort.
"""
function sorttwo!(vs::AbstractVector{T}, index, lo::Int = 1, hi::Int=length(vs)) where T
    # Input checking
    if lo >= hi;  return (vs, index);  end

    # Make sure we're sorting a bits type
    # T = Base.Order.ordtype(o, vs)
    if !isbits(T)
        error("Radix sort only sorts bits types (got $T)")
    end
    o = Forward

    # Init
    iters = ceil(Integer, sizeof(T)*8/RADIX_SIZE)
    # number of buckets in the counting step
    nbuckets = 2^RADIX_SIZE
    bin = zeros(UInt32, nbuckets, iters)
    if lo > 1;  bin[1,:] = lo-1;  end

    # Histogram for each element, radix
    for i = lo:hi
        v = uint_mapping(o, vs[i])
        for j = 1:iters
            idx = Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK) + 1
            @inbounds bin[idx,j] += 1
        end
    end

    # Sort!
    swaps = 0
    len = hi-lo+1

    index1 = similar(index)
    ts=similar(vs)
    for j = 1:iters
        # Unroll first data iteration, check for degenerate case
        v = uint_mapping(o, vs[hi])
        idx = Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK) + 1

        # are all values the same at this radix?
        if bin[idx,j] == len;  continue;  end

        # cbin = cumsum(bin[:,j])
        # tries to achieve the above one-liner with more efficiency
        cbin = zeros(UInt32, nbuckets)
        cbin[1] = bin[1,j]
        for i in 2:nbuckets
            cbin[i] = cbin[i-1] + bin[i,j]
        end

        ci = cbin[idx]
        ts[ci] = vs[hi]
        index1[ci] = index[hi]
        cbin[idx] -= 1

        # Finish the loop...
        @inbounds for i in hi-1:-1:lo
            v = uint_mapping(o, vs[i])
            idx = Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK) + 1
            ci = cbin[idx]
            ts[ci] = vs[i]
            index1[ci] = index[i]
            cbin[idx] -= 1
        end
        vs,ts = ts,vs
        index, index1 = index1, index
        swaps += 1
    end

    if isodd(swaps)
        vs,ts = ts,vs
        index, index1 = index1, index
        for i = lo:hi
            @inbounds vs[i] = ts[i]
            @inbounds index[i] = index1[i]
        end
    end
    (vs, index)
end

"""
    radixsort(svec, rev = false)
"""
radixsort(svec, rev = false) = radixsort!(copy(svec), rev)
radixsort!(svec, rev = false) =  sort!(svec, 1, length(svec), StringRadixSort, rev? Base.Reverse:Base.Forward)

"""
    sortperm_radixsort(svec, rev = nothing, order = Forward)

To return a `String` vector using LSD radixsort
"""
function sortperm_radixsort(svec::AbstractVector{String}; rev::Union{Bool,Void}=nothing, order::Ordering=Forward)
    sortperm_radixsort!(copy(svec), rev = rev, order =order)
end

function sortperm_radixsort!(svec::AbstractVector{String}; rev::Union{Bool,Void}=nothing, order::Ordering=Forward)
    siv = StringIndexVector(svec, collect(1:length(svec)))

    # find the maximum string length
    lens = reduce((x,y) -> max(x,sizeof(y)),0, svec)
    skipbytes = lens
    while lens > 0
       if lens > 8
            skipbytes = max(0, skipbytes - 16)
            if order == Reverse
                sorttwo!(.~load_bits.(UInt128, siv.svec, skipbytes), siv)
            else
                sorttwo!(load_bits.(UInt128, siv.svec, skipbytes), siv)
            end
            lens -= 16
        elseif lens > 4
            skipbytes = max(0, skipbytes - 8)
            if order == Reverse
                sorttwo!(.~load_bits.(UInt64, siv.svec, skipbytes), siv)
            else
                sorttwo!(load_bits.(UInt64, siv.svec, skipbytes), siv)
            end
            lens -= 8
        else
            skipbytes = max(0, skipbytes - 4)
            if order == Reverse
                sorttwo!(.~load_bits.(UInt32, siv.svec, skipbytes), siv)
            else
                sorttwo!(load_bits.(UInt32, siv.svec, skipbytes), siv)
            end
            lens -= 4
        end
    end
    siv.index
end

"Simple data structure for carrying a string vector and its index; this allows
`sorttwo!` to sort the radix of the string vector and reorder the string and its
index at the same time opening the way for faster sort_perm"
struct StringIndexVector
    svec::Vector{String}
    index::Vector{Int}
end

function setindex!(siv::StringIndexVector, X::StringIndexVector, inds)
    siv.svec[inds] = X.svec
    siv.index[inds] = X.index
end

function setindex!(siv::StringIndexVector, X, inds)
    siv.svec[inds] = X[1]
    siv.index[inds] = X[2]
end

getindex(siv::StringIndexVector, inds::Integer) = siv.svec[inds], siv.index[inds]
getindex(siv::StringIndexVector, inds...) = StringIndexVector(siv.svec[inds...], siv.index[inds...])
similar(siv::StringIndexVector) = StringIndexVector(similar(siv.svec), similar(siv.index))


text:
@enum TypePoisson StdPoisson = 1 StdPoisson2d = 2  StdABp = 3

function _get_fctv_k(adv::Advection{T,N,timeopt}) where {T,N,timeopt}
    fct_k(v) = im / sum(v .^ 2)
    Nsp = div(N, 2)
    v_k = vec_k_fft.(adv.t_mesh[1:Nsp])
    sz = length.(adv.t_mesh[1:Nsp])
    fctv_k_gen = fct_k.(collect(Iterators.product(v_k...)))
    fctv_k_gen[1] = 0
    return ntuple(x -> reshape(v_k[x], tupleshape(x, Nsp, sz[x])) .* fctv_k_gen, Nsp)
end

# using SHA
# cod(a::Array) = bytes2hex(sha256(reinterpret(UInt8, collect(Iterators.flatten(a)))))

"""
    PoissonConst{T, Nsp, Nv}
    PoissonConst(adv::Advection{T, Nsp, Nv, Nsum}; isfftbig=true)

Constant data for the computation of poisson coefficients

# Arguments
- `adv::Advection{T, Nsp, Nv, Nsum, timeopt}` : Advection constant data
- `isfftbig=true`: if true compute the fttbig structure

# Implementation
- `adv` : Advection constant data
- `v_k' : vector of vector of fourier coefficents for integration for each space dimension
- `fctv_k` : Array of space dimensions of the inverse of the norm of fourier coefficients
- `pfftbig` : Fourier data for space dimensions
"""
struct PoissonConst{T,N,Nsp,Nv,type,typeadd}
    adv::Advection
    fctv_k::Any
    v_square::Array{T,Nv}
    pfftbig::Any
    abcoef::ABcoef
    function PoissonConst(
        adv::Advection{T,N,timeopt};
        isfftbig = true,
        type::TypePoisson = StdPoisson,
        typeadd = 0,
    ) where {T,N,timeopt}
        N % 2 == 0 || thrown(ArgumentError("N=$N must be a multiple of 2"))

        Nsp = Nv = div(N, 2)
        fctv_k = _get_fctv_k(adv)
        v_square = dotprod(points.(adv.t_mesh[(Nsp+1):N])) .^ 2 # precompute for ke
        pfftbig = if isfftbig && T != Float64
            PrepareFftBig(sizeall(adv)[1:Nsp], T; numdims = Nsp, dims = ntuple(x -> x, Nsp))
        else
            missing
        end

        return new{T,N,Nsp,Nv,type,typeadd}(adv, fctv_k, v_square, pfftbig, ABcoef(typeadd))
    end
end
getNspNv(pc::PoissonConst{T,N,Nsp,Nv}) where {T,N,Nsp,Nv} = (Nsp, Nv)

"""
    PoissonVar{T, N, Nsp, Nv} <: AbstractExtDataAdv{T}
    PoissonVar(pc::PoissonConst{T, N, Nsp, Nv})

mutable structure of variable data for the poisson computation

# Arguments
- `pc::PoissonConst{T, N, Nsp, Nv}` : poisson constant data

# Implementation
- `pc::PoissonConst{T, Nsp, Nv}` : poisson constant data
- `rho::Array{T, Nsp}` : result of the compute_charge that is the sum along velocity dimensions
- `t_elfield::NTuple{Nsp,Array{Complex{T}, Nsp}}` : electric fields initialized at each beginning of velocity advection subseries
"""
mutable struct PoissonVar{T,N,Nsp,Nv,type,typeadd} <: AbstractExtDataAdv
    pc::PoissonConst{T,N,Nsp,Nv,type,typeadd}
    rho::Array{T,Nsp}
    t_elfield::Union{NTuple{Nsp,Array{T,Nsp}},Missing}
    bufcur_sp::Any
    bufcur_v::Any
    tupleind::Any

    function PoissonVar(
        pc::PoissonConst{T,N,Nsp,Nv,type,typeadd},
    ) where {T,N,Nsp,Nv,type,typeadd}
        sz = length.(pc.adv.t_mesh[1:Nsp])
        rho = Array{T,Nsp}(undef, sz)

        return new{T,N,Nsp,Nv,type,typeadd}(
            pc,
            rho,
            missing,
            missing,
            missing,
            missing
        )
    end
end
function getpoissonvar(adv::Advection; type::TypePoisson = StdPoisson, typeadd = 0)
    pc = PoissonConst(adv, type = type, typeadd = typeadd)
    return PoissonVar(pc)
end

getNspNv(_::PoissonVar{T,N,Nsp,Nv}) where {T,N,Nsp,Nv} = (Nsp, Nv)
"""
    compute_charge!( self::PoissonVar, advd::AdvectionData)

 Compute charge density

 ρ(x,t) = ∫ f(x,v,t) dv

 # Argument
 - `self::AdvectionData` : mutable structure of variables data.
"""
function compute_charge!(self::PoissonVar{T,N, Nsp,Nv}, advd::AdvectionData) where {T,N, Nsp,Nv}
    compute_charge!(self.rho, self.pc.adv.t_mesh[end-Nv+1:end], advd.data)
    missing
end

"""
    compute_elfield!( self:PoissonVar)

computation of electric field
    ∇.e = - ρ

# Argument
 - `self::PoissonVar` : mutable structure of variables data.

"""
function compute_elfield!(self::PoissonVar{T,N, Nsp,Nv}) where {T,N, Nsp,Nv}
    pfft = self.pc.pfftbig
    buf = fftgenall(pfft, self.rho)
    self.t_elfield = ntuple(x -> real(ifftgenall(pfft, self.pc.fctv_k[x] .* buf)), Nsp)
    missing
end

"""
    initcoef!(pv::PoissonVar{T, N,Nsp, Nv}, self::AdvectionData{T, N})

Implementation of the interface function that is called at the begining of each advection
    This is implementation for Vlasov-Poisson equation

"""
function isvelocity(pv::PoissonVar{T,N,Nsp,Nv}, advd::AdvectionData{T,N}) where {T,N,Nsp,Nv}
    st = getst(advd)
    return st.perm[1] > Nsp
end
function initcoef!(
    pv::PoissonVar{T,N,Nsp,Nv,StdPoisson},
    advd::AdvectionData{T,N},
) where {T,N,Nsp,Nv}
    st = getst(advd)
    adv = advd.adv
    dt = getcur_t(advd)
    @show dt, isvelocity(pv, advd)
    if isvelocity(pv, advd)
        if (Nsp + 1) in st.perm[1:st.ndims]
            compute_charge!(pv, advd)
            compute_elfield!(pv)
        end
#        @show dt
        bc_v = ntuple(
            x ->
                (dt / step(adv.t_mesh[st.perm[x]])) *
                pv.t_elfield[st.perm[x]-Nsp],
            st.ndims,
        )

        if ismissing(pv.bufcur_v)
            pv.bufcur_v = bc_v
        else
            for x = 1:length(bc_v)
                pv.bufcur_v[x] .= bc_v[x]
            end
        end
    else
        pv.tupleind = ntuple(x -> st.perm[st.invp[x]+Nsp] - st.ndims, st.ndims)
        bc_sp = ntuple(
            x ->
                (-dt / step(adv.t_mesh[st.invp[x]])) *
                adv.t_mesh[st.invp[x]+Nsp].points,
            st.ndims,
        )
        if ismissing(pv.bufcur_sp)
            pv.bufcur_sp = bc_sp
        else
            for x = 1:length(bc_sp)
                pv.bufcur_sp[x] .= bc_sp[x]
            end
        end
    end
end
@inline function getalpha(
    pv::PoissonVar{T,N,Nsp,Nv,StdPoisson},
    advd::AdvectionData{T},
    ind,
) where {T,N,Nsp,Nv}
    st = getst(advd)
    return if isvelocity(pv, advd)
        ntuple(
            x -> pv.bufcur_v[x][CartesianIndex(ind.I[end-Nsp+1:end])],
            size(pv.bufcur_v, 1),
        )
    else
        ntuple(x -> pv.bufcur_sp[x][ind.I[pv.tupleind[x]]], size(pv.bufcur_sp, 1))
    end
end

function initcoef!(
    pv::PoissonVar{T,N,Nsp,Nv,StdPoisson2d},
    advd::AdvectionData{T,N, timeopt},
) where {T,N,Nsp,Nv, timeopt}
    #    st = getst(self)
    adv = advd.adv
    compute_charge!(pv, advd)
    compute_elfield!(pv)

    bufc_v = (getcur_t(advd) / step(adv.t_mesh[2])) * pv.t_elfield[1]
    bufc_sp = (-getcur_t(advd) / step(adv.t_mesh[1])) * adv.t_mesh[2].points

    if ismissing(advd.bufcur)
        advd.bufcur = zeros(OpTuple{N,T}, sizeall(adv))
    end
    for ind in CartesianIndices(sizeall(adv))
        advd.bufcur[ind] = OpTuple((bufc_sp[ind.I[2]], bufc_v[ind.I[1]]))
    end
end

ext: jl


text: tonamedtuple(d::Dict) = NamedTuple{Tuple(keys(d))}(values(d))


text:
"""
$(TYPEDEF)

A `Body` is a component of a [`Mechanism`](@ref).
# Important attributes
* `id`:    The unique ID of a body. Assigned when added to a `Mechanism`.
* `name`:  The name of a body. The name is taken from a URDF or can be assigned by the user.
* `m`:     The mass of a body.
* `J`:     The inertia of a body.
* `state`: The state of a body. Contains all position and velocity information (see [`State`](@ref)).
* `shape`: The visualization shape of a body (see [`Shape`](@ref)).

# Constructors
    Body(m, J; name, shape)
    Mesh(path, m, J; scale, kwargs...)
    Box(x, y, z, m; kwargs...)
    Cylinder(r, h, m; kwargs...)
    Sphere(r, m; kwargs...)
    Pyramid(w, h, m; kwargs...)
"""
mutable struct Body{T} <: AbstractBody{T}
    id::Int64
    name::String
    active::Bool

    m::T
    J::SMatrix{3,3,T,9}

    state::State{T}

    shape::Shape{T}

    function Body(m::Real, J::AbstractArray; name::String="", shape::Shape=EmptyShape())
        T = promote_type(eltype.((m, J))...)
        new{T}(getGlobalID(), name, true, m, J, State{T}(), shape)
    end

    function Body{T}(contents...) where T
        new{T}(getGlobalID(), contents...)
    end
end

"""
$(TYPEDEF)

The `Origin` is the root of a [`Mechanism`](@ref).
# Important attributes
* `id`:    The unique ID of the origin. Assigned when added to a `Mechanism`.
* `name`:  The name of the origin. The name is taken from a URDF or can be assigned by the user.
* `shape`: The visualization shape of the origin (see [`Shape`](@ref)).

# Constructors
    Origin(; name, shape)
    Origin{Type}(; name, shape)
    Origin(body)
"""
mutable struct Origin{T} <: AbstractBody{T}
    id::Int64
    name::String

    shape::Shape{T}

    function Origin{T}(; name::String="", shape::Shape=EmptyShape()) where T
        new{T}(getGlobalID(), name, shape)
    end

    function Origin(; name::String="", shape::Shape=EmptyShape())
        Origin{Float64}(; name = name, shape = shape)
    end

    function Origin(body::Body{T}) where T
        new{T}(body.id, body.name, body.shape)
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, body::Body{T}) where {T}
    summary(io, body)
    println(io,"")
    println(io, " id:     "*string(body.id))
    println(io, " name:   "*string(body.name))
    println(io, " active: "*string(body.active))
    println(io, " m:      "*string(body.m))
    println(io, " J:      "*string(body.J))
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, origin::Origin{T}) where {T}
    summary(io, origin)
    println(io, " id:   "*string(origin.id))
    println(io, " name: "*string(origin.name))
end

function Base.deepcopy(b::Body{T}) where T
    contents = []
    for i = 2:getfieldnumber(b)
        push!(contents, deepcopy(getfield(b, i)))
    end

    return Body{T}(contents...)
end

Base.length(::Body) = 6

@inline getM(body::Body{T}) where T = [[I*body.m;szeros(T,3,3)] [szeros(T,3,3);body.J]]

@inline function ∂dyn∂vel(mechanism, body::Body{T}) where T
    state = body.state
    Δt = mechanism.Δt
    J = body.J
    ω2 = state.ωsol[2]
    sq = sqrt(4 / Δt^2 - ω2' * ω2)

    dynT = I * body.m / Δt
    dynR = skewplusdiag(ω2, sq) * J - J * ω2 * (ω2' / sq) - skew(J * ω2)

    Z = szeros(T, 3, 3)

    state.D = [[dynT; Z] [Z; dynR]]

    for connectionid in damperconnections(mechanism.graph, body.id)
        damperToD!(mechanism, body, geteqconstraint(mechanism, connectionid))
    end

    return state.D
end

# Derivatives for linearizations

function ∂F∂z(body::Body{T}, Δt) where T
    state = body.state
    Z3 = szeros(T,3,3)
    Z76 = szeros(T,7,6)
    Z6 = szeros(T,6,6)
    E3 = SMatrix{3,3,T,9}(I)

    AposT = [-I -E3*Δt] # TODO is there a UniformScaling way for this instead of E3?

    AvelT = [Z3 -I*body.m/Δt]

    AposR = [-Rmat(ωbar(state.ωc, Δt)*Δt/2)*LVᵀmat(state.qc) -Lmat(state.qc)*derivωbar(state.ωc, Δt)*Δt/2]

    J = body.J
    ω1 = state.ωc
    sq1 = sqrt(4 / Δt^2 - ω1' * ω1)
    ω1func = skewplusdiag(ω1, -sq1) * J + J * ω1 * (ω1' / sq1) - skew(J * ω1)
    AvelR = [Z3 ω1func]

    return [[AposT;AvelT] Z6;Z76 [AposR;AvelR]]
end

function ∂F∂u(body::Body{T}, Δt) where T
    Z3 = szeros(T,3,3)
    Z43 = szeros(T,4,3)

    BposT = [Z3 Z3] # TODO is there a UniformScaling way for this instead of E3?

    BvelT = [-I Z3]

    BposR = [Z43 Z43]

    BvelR = [Z3 -2.0*I]

    return [BposT;BvelT;BposR;BvelR]
end

function ∂F∂fz(body::Body{T}, Δt) where T
    state = body.state
    Z3 = szeros(T,3,3)
    Z34 = szeros(T,3,4)
    Z43 = szeros(T,4,3)
    Z67 = szeros(T,6,7)
    Z76 = szeros(T,7,6)

    AposT = [I Z3]

    AvelT = [Z3 I*body.m/Δt]

    AposR = [I Z43]

    J = body.J
    ω2 = state.ωsol[2]
    sq2 = sqrt(4 / Δt^2 - ω2' * ω2)
    ω2func = skewplusdiag(ω2, sq2) * J - J * ω2 * (ω2' / sq2) - skew(J * ω2)
    AvelR = [Z34 ω2func]

    return [[AposT;AvelT] Z67;Z76 [AposR;AvelR]], [I Z3 Z34 Z3;Z3 I Z34 Z3;Z3 Z3 VLᵀmat(state.qsol[2]) Z3;Z3 Z3 Z34 I]
end


text:
isdefined(Base, :__precompile__) && __precompile__()

module DiffBase

using Compat

include("results.jl")
include("testfuncs.jl")

######################
# Test Function Sets #
######################
#=
These functions are organized in sets based on input/output type. They are unary and not
in-place unless otherwised specified. These functions have been written with the following
assumptions:

- Array input is of length >= 5
- Matrix input is square
- Matrix inputs for n-ary functions are of equal shape
=#

const NUMBER_TO_NUMBER_FUNCS = (num2num_1, num2num_2, num2num_3,
                                num2num_4, num2num_5, num2num_6)

const NUMBER_TO_ARRAY_FUNCS = (num2arr_1,)

const VECTOR_TO_NUMBER_FUNCS = (vec2num_1, vec2num_2,  vec2num_3, vec2num_4, vec2num_5,
                                rosenbrock_1, rosenbrock_2, rosenbrock_3, rosenbrock_4,
                                ackley, self_weighted_logit, first)

const MATRIX_TO_NUMBER_FUNCS = (det, mat2num_1, mat2num_2, mat2num_3, mat2num_4, softmax)

const INPLACE_ARRAY_TO_ARRAY_FUNCS = (chebyquad!, brown_almost_linear!, trigonometric!,
                                      mutation_test_1!, mutation_test_2!)

const ARRAY_TO_ARRAY_FUNCS = (-, chebyquad, brown_almost_linear, trigonometric, arr2arr_1,
                              arr2arr_2, mutation_test_1, mutation_test_2, identity)

const MATRIX_TO_MATRIX_FUNCS = (inv,)

if VERSION >= v"0.6.0-dev.1614"
    const BINARY_BROADCAST_OPS = ((a, b) -> broadcast(+, a, b),
                                  (a, b) -> broadcast(-, a, b),
                                  (a, b) -> broadcast(*, a, b),
                                  (a, b) -> broadcast(/, a, b),
                                  (a, b) -> broadcast(\, a, b),
                                  (a, b) -> broadcast(^, a, b))
else
    const BINARY_BROADCAST_OPS = (.+, .-, .*, ./, .\, .^)
end

const BINARY_MATRIX_TO_MATRIX_FUNCS = (+, -, *, /, \,
                                       BINARY_BROADCAST_OPS...,
                                       A_mul_Bt, At_mul_B, At_mul_Bt,
                                       A_mul_Bc, Ac_mul_B, Ac_mul_Bc)

const TERNARY_MATRIX_TO_NUMBER_FUNCS = (neural_step,)

end # module


text:
module Thermo

    # Write your package code here.

end


text:
module TestElixirs

using LinearAlgebra
using Test
using Trixi

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
const EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples")

@testset "Special elixirs" begin
  @testset "Convergence test" begin
    mean_values = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "2d", "elixir_advection_extended.jl"), 3)
    @test isapprox(mean_values[:l2], [4.0], rtol=0.01)

    mean_values = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "paper-self-gravitating-gas-dynamics", "elixir_eulergravity_eoc_test.jl"), 2, tspan=(0.0, 0.1))
    @test isapprox(mean_values[:l2], 4 * ones(4), atol=0.4)
  end

  @testset "Test linear structure (2D)" begin
    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "2d", "elixir_advection_extended.jl"),
                  tspan=(0.0, 0.0), initial_refinement_level=2)
    A, b = linear_structure(semi)
    λ = eigvals(Matrix(A))
    @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))

    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "2d", "elixir_hypdiff_lax_friedrichs.jl"),
                  tspan=(0.0, 0.0), initial_refinement_level=2)
    A, b = linear_structure(semi)
    λ = eigvals(Matrix(A))
    @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))

    # check whether the user can modify `b` without changing `A`
    x = vec(ode.u0)
    Ax = A * x
    @. b = 2 * b + x
    @test A * x ≈ Ax
  end

  @testset "Test Jacobian of DG (2D)" begin
    @testset "Linear advection" begin
      trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "2d", "elixir_advection_extended.jl"),
                    tspan=(0.0, 0.0), initial_refinement_level=2)
      A, _ = linear_structure(semi)

      J = jacobian_ad_forward(semi)
      @test Matrix(A) ≈ J
      λ = eigvals(J)
      @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))

      J = jacobian_fd(semi)
      @test Matrix(A) ≈ J
      λ = eigvals(J)
      @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))
    end

    @testset "Compressible Euler equations" begin
      trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "2d", "elixir_euler_density_wave.jl"),
                    tspan=(0.0, 0.0), initial_refinement_level=2)

      J = jacobian_ad_forward(semi)
      λ = eigvals(J)
      @test maximum(real, λ) < 7.0e-7

      J = jacobian_fd(semi)
      λ = eigvals(J)
      @test maximum(real, λ) < 7.0e-3

      trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "2d", "elixir_euler_shockcapturing.jl"),
                    tspan=(0.0, 0.0), initial_refinement_level=1)
      # This does not work yet because of the indicators...
      @test_skip jacobian_ad_forward(semi)
    end

    @testset "MHD" begin
      trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "2d", "elixir_mhd_alfven_wave.jl"),
                    tspan=(0.0, 0.0), initial_refinement_level=1)
      @test_nowarn jacobian_ad_forward(semi)

      trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "2d", "elixir_mhd_alfven_wave_mortar.jl"),
                    tspan=(0.0, 0.0), initial_refinement_level=1)
      @test_nowarn jacobian_ad_forward(semi)
    end
  end

  @testset "Test linear structure (3D)" begin
    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "3d", "elixir_advection_extended.jl"),
                  tspan=(0.0, 0.0), initial_refinement_level=1)
    A, b = linear_structure(semi)
    λ = eigvals(Matrix(A))
    @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))
  end

  @testset "Test Jacobian of DG (3D)" begin
    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "3d", "elixir_advection_extended.jl"),
                  tspan=(0.0, 0.0), initial_refinement_level=1)
    A, _ = linear_structure(semi)

    J = jacobian_ad_forward(semi)
    @test Matrix(A) ≈ J
    λ = eigvals(J)
    @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))

    J = jacobian_fd(semi)
    @test Matrix(A) ≈ J
    λ = eigvals(J)
    @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end #module


text:
function SelectNode(Tree;VariableDef = 0, MU = 0, RO1 = 0, RO2 = 0, Type = 0, Type2Flag = 0)

    const  EPS = 0.000001

    # Type = (Type_of_selection, Key_Having_Feasible_Sol=1)
    if Type == 0 || isempty(Type)
        return length(Tree)
    elseif Type == 1
        return MinimalT(Tree)[1]
    elseif Type == 2
        return SelectNode(Tree,Type = Type2Flag)
    elseif Type == 3
        ExpectedImprove = zeros(length(Tree))
        for i = 1:length(Tree)
            delta1 = (ceil.(Tree[i].X[VariableDef[1]+1:end])-Tree[i].X[VariableDef[1]+1:end]).*RO1./(EPS+n1)
            delta2 = (Tree[i].X[VariableDef[1]+1:end]-floor.(Tree[i].X[VariableDef[1]+1:end])).*RO2./(EPS+n2)
            NonIntegers = find(.!isinteger.(Tree[i].X[VariableDef[1]+1:end]))
            S = MU*min.(delta1,delta2) + (1-MU)*max.(delta1,delta2)
            ind = findmax(S[NonIntegers])[2]
            ExpectedImprove[i] = max.(delta1[ind]+Tree[i].LB, delta2[ind]+Tree[i].LB)
        end
        ind = findmax(ExpectedImprove)[2]
        return ind
    elseif Type == 4
        ExpectedImprove = zeros(length(Tree));
        for i = 1:length(Tree)
            delta1 = (ceil.(Tree[i].X[VariableDef[1]+1:end])-Tree[i].X[VariableDef[1]+1:end]).*RO1./(EPS+n1)
            delta2 = (Tree[i].X[VariableDef[1]+1:end]-floor.(Tree[i].X[VariableDef[1]+1:end])).*RO2./(EPS+n2)
            delta = min.(delta1, delta2);
            NonIntegers = find(.!isinteger.(Tree[i].X[VariableDef[1]+1:end]))
            ExpectedImprove[i] = sum(delta[NonIntegers])+Tree[i].LB
        end
        ind = findmax(ExpectedImprove)[2]
        return ind
    end
end


text:
# Exports
export construct_graph
export solve_graph
export satellite_plan_graph

"""
Construct decision graph for solving network
"""
function construct_graph(problem::SatPlanningProblem; sat_id::Integer=1)
    # Sort opportunities to ensure time ascending order
    sort!(problem.opportunities, by = x -> x.t_start)

    # Get max slew time
    sc = problem.spacecraft[sat_id]

    graph = Dict{Opportunity, Array{<:Opportunity, 1}}()

    for start_opp in problem.opportunities
        # Create list of valid transitions for start_collect
        graph[start_opp] = Opportunity[]

        for end_opp in problem.opportunities

            if (start_opp == end_opp) || (start_opp.t_start >= end_opp.t_end)
               # Conditions to skip insertion
            else
                # Valuate transition to see if valid
                valid = true

                # Check to see if each transition constraint is valid
                for constraint in problem.constraints
                    if valid == false
                        break
                    end

                    valid = valid && constraint(start_opp, end_opp)
                end

                # If valid transition add to edges
                if valid == true
                    push!(graph[start_opp], end_opp)
                end

            end
        end
    end

    return graph
end

"""
Solve a graph to find optimal transition path
"""
function solve_graph(problem::SatPlanningProblem, graph::Dict{Opportunity, Array{<:Opportunity, 1}};
            contact_reward::Real=0.0)

    # Create graph based on incoming edges for each node
    igraph = Dict{Opportunity, Array{<:Opportunity, 1}}()

    for opp in problem.opportunities
        igraph[opp] = Opportunity[]
    end

    # For each node in the original graph add the source node as a possible
    # predecessor to each outgoing edge it has
    for node in sort(collect(keys(graph)), by = x -> x.t_start)
        for edge in graph[node]
            push!(igraph[edge], node)
        end
    end

    # Initialize structure to store optimal decision at each node
    optimal_path = Dict{Opportunity, Tuple{Union{Opportunity, Nothing}, Float64, Array{Request, 1}}}()

    for opp in problem.opportunities
        # Initialize node weight
        if typeof(opp) == Collect
            if problem.solve_allow_repeats == false
                optimal_path[opp] = (nothing, opp.reward, Request[opp.location])
            else
                optimal_path[opp] = (nothing, opp.reward, Request[])
            end
        else
            optimal_path[opp] = (nothing, 0.0, Request[])
        end
    end

    # For each node iterate over all incoming nodes
    for nodej in sort(collect(keys(igraph)), by = x -> x.t_start)
        # Iterate over all possible incoming nodes to node j
        for nodei in sort(igraph[nodej], by = x -> x.t_start)
            # Compute reward for current action
            rij = optimal_path[nodei][2]
            if typeof(nodej) == Collect
                # Calculate reward for taking j after the node i path
                if problem.solve_allow_repeats == true || !(nodej.location in optimal_path[nodei][3])
                    rij += nodej.reward
                end
            elseif typeof(nodej) == Contact
                rij += contact_reward
            end

            # If incoming node i is higher than anything seen previously set it
            # as new predecessory
            if rij > optimal_path[nodej][2]
                request_history = copy(optimal_path[nodei][3])
                if problem.solve_allow_repeats == false && typeof(nodej.location) == Request
                    push!(request_history, nodej.location)
                end
                optimal_path[nodej] = (nodei, rij, request_history)
            end
        end
    end

    return optimal_path
end

"""
Schedule satellite activity using graph planning.

Returns:
- `plan::Array{<:Opportunity, 1}`
- `reward::Float64`
- `requests::{<:Request, 1}`
"""
function satellite_plan_graph(problem::SatPlanningProblem; sat_id::Integer=1,
            contact_reward::Real=0)
    # Construct graph to solve
    graph = construct_graph(problem, sat_id=sat_id)

    # Solve graph
    optimal_path = solve_graph(problem, graph,
                        contact_reward=contact_reward)

    # Extract prameters of optimal path
    path     = Opportunity[]
    reward   = 0.0

    # Find node with highest reward:
    current_max = -Inf
    end_node = nothing
    for node in keys(optimal_path)
        if optimal_path[node][2] > current_max
            current_max = optimal_path[node][2]
            end_node    = node
        end
    end
    reward = current_max

    # Backtrack from end node to get optimal path
    push!(path, end_node)
    current_node = end_node
    while optimal_path[current_node][1] != nothing
        previous_node = optimal_path[current_node][1]
        push!(path, previous_node)
        current_node = previous_node
    end

    # Reverse path so it's in chronological order
    path = reverse!(path)

    return path, reward
end


text:
using Test
using Statistics
using LESbrary.Utils
using LESbrary.Utils: ∂z, ∂t

@testset "InterpolatedProfile" begin

    # interpolation

    zs = 0:5:50
    Δz = zs.step
    N = length(zs)
    data = randn(N)

    prof = InterpolatedProfile(data, zs, Δz, zs[end])

    @test prof(0) == prof.data[1]
    @test prof(5) == prof.data[2]

    @test prof(12.5) ≈ (prof.data[3] + prof.data[4]) / 2

    # derivatives

    zs = 0:5:50
    Δz = zs.step
    N = length(zs)

    f(z) = z^2
    data = f.(zs)

    prof = InterpolatedProfile(data, zs, Δz, zs[end])

    @test ∂z(prof, 2) ≈ (f(5) - f(0)) / Δz
    @test ∂z(prof, 7.4) ≈ (f(10) - f(5)) / Δz
    @test ∂z(prof, 11.8) ≈ (f(15) - f(10)) / Δz

    @test ∂z(prof, 2) == ∂z(prof, 4)
    @test ∂z(prof, 7.4) == ∂z(prof, 9.99)
    @test ∂z(prof, 11.8) == ∂z(prof, 10.01)
end

@testset "InterpolatedProfileTimeSeries" begin

    # interpolation

    zs = 0:5:50
    ts = 0:60:600
    Δz = zs.step
    Δt = ts.step
    Nz = length(zs)
    Nt = length(ts)
    data = randn(Nz, Nt)

    prof = InterpolatedProfileTimeSeries(data, zs, ts, Δz, Δt, zs[end], ts[end])

    @test prof(0, 0) == prof.data[1, 1]
    @test prof(5, 0) == prof.data[2, 1]
    @test prof(0, 60) == prof.data[1, 2]
    @test prof(5, 60) == prof.data[2, 2]

    @test prof(27.5, 330) ≈ mean(prof.data[6:7, 6:7])

    # derivatives

    zs = 0:5:50
    ts = 0:60:600
    Δz = zs.step
    Δt = ts.step
    Nz = length(zs)
    Nt = length(ts)

    f(z, t) = 3z^2 + 0.1t^2
    data = f.(zs, ts')

    prof = InterpolatedProfileTimeSeries(data, zs, ts, Δz, Δt, zs[end], ts[end])

    @test ∂z(prof, 3, 0) ≈ (f(5, 0) - f(0, 0)) / Δz
    @test ∂z(prof, 7.5, 10) ≈ (f(10, 10) - f(5, 10)) / Δz
    @test ∂z(prof, 10.01, 420.69) ≈ (f(15, 420.69) - f(10, 420.69)) / Δz

    @test ∂t(prof, 0, 45) ≈ (f(0, 60) - f(0, 0)) / Δt
    @test ∂t(prof, 6, 100) ≈ (f(6, 120) - f(6, 60)) / Δt
    @test ∂t(prof, 22.2, 310) ≈ (f(22.2, 360) - f(22.2, 300)) / Δt
end


text:
function geom(x::Blade, y::Blade)
    i, ρ = TABLE[x.index, y.index]
    Blade(ρ * x.coef * y.coef, i)
end


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    @metasolver solver solvertype body

A helper macro to create a solver named `solver` of type `solvertype`
with parameters specified in `body`.

## Examples

Create a solver with parameters `mean` and `variogram` for each variable
of the problem, and a global parameter that specifies whether or not
to use the GPU:

```julia
julia> @metasolver MySolver AbstractSimulationSolver begin
  @param mean = 0.
  @param variogram = GaussianVariogram()
  @global gpu = false
end
```

### Notes

This macro is not intended to be used directly, see other macros defined
below for estimation and simulation solvers.
"""
macro metasolver(solver, solvertype, body)
  # discard any content that doesn't start with @param or @global
  content = filter(arg -> arg isa Expr, body.args)

  # lines starting with @param refer to variable parameters
  vparams = filter(p -> p.args[1] == Symbol("@param"), content)
  vparams = map(p -> p.args[3], vparams)

  # lines starting with @global refer to global solver parameters
  gparams = filter(p -> p.args[1] == Symbol("@global"), content)
  gparams = map(p -> p.args[3], gparams)

  # add default value of `nothing` if necessary
  gparams = map(p -> p isa Symbol ? :($p = nothing) : p, gparams)

  # replace Expr(:=, a, 2) by Expr(:kw, a, 2) for valid kw args
  gparams = map(p -> Expr(:kw, p.args...), gparams)

  # keyword names
  gkeys = map(p -> p.args[1], gparams)

  # solver parameter type
  solverparam = Symbol(solver,"Param")

  esc(quote
    $Parameters.@with_kw_noshow struct $solverparam
      __dummy__ = nothing
      $(vparams...)
    end

    @doc (@doc $solverparam) (
    struct $solver <: $solvertype
      params::Dict{Symbol,$solverparam}

      $(gkeys...)

      function $solver(params::Dict{Symbol,$solverparam}, $(gkeys...))
        new(params, $(gkeys...))
      end
    end)

    function $solver(params...; $(gparams...))
      # build dictionary for inner constructor
      dict = Dict{Symbol,$solverparam}()

      # convert named tuples to solver parameters
      for (varname, varparams) in params
        kwargs = [k => v for (k,v) in zip(keys(varparams), varparams)]
        push!(dict, varname => $solverparam(; kwargs...))
      end

      $solver(dict, $(gkeys...))
    end

    # ------------
    # IO methods
    # ------------
    function Base.show(io::IO, solver::$solver)
      print(io, $solver)
    end

    function Base.show(io::IO, ::MIME"text/plain", solver::$solver)
      println(io, solver)
      for (var, varparams) in solver.params
        println(io, "  └─$var")
        pnames = setdiff(fieldnames(typeof(varparams)), [:__dummy__])
        for pname in pnames
          pval = getfield(varparams, pname)
          if pval ≠ nothing
            print(io, "    └─$pname = ")
            show(IOContext(io, :compact => true), pval)
            println(io, "")
          end
        end
      end
    end
  end)
end

"""
    @estimsolver solver body

A helper macro to create a estimation solver named `solver` with parameters
specified in `body`. For examples, please check the documentation for
`@metasolver`.
"""
macro estimsolver(solver, body)
  esc(quote
    GeoStatsDevTools.@metasolver $solver GeoStatsBase.AbstractEstimationSolver $body
  end)
end

"""
    @estimsolver solver body

A helper macro to create a simulation solver named `solver` with parameters
specified in `body`. For examples, please check the documentation for
`@metasolver`.
"""
macro simsolver(solver, body)
  esc(quote
    GeoStatsDevTools.@metasolver $solver GeoStatsBase.AbstractSimulationSolver $body
  end)
end


text:
using Optim
include("NN_util.jl")

mutable struct NNModel
    Nvar     # number of ODE equations
    f_list   # f(t, y) on the right hand side
    y0_list  # initial conditions
    t        # training points
    t0       # initial t, i.e. t[1]
    n_hidden # number of hidden units

    params_list    # a list of NN parameters
    sizes   # size of NN parameters, for flatten operation
    lengths # length of NN parameters, for flatten operation
    total_l  # total length of NN parameters, for flatten operation
    all_flattened_params # flattenned parameters, for Optim.jl
end

function show(nn::NNModel)
    println("Neural ODE Solver")
    println("Number of equations:       ", nn.Nvar)
    println("Initial condition y0:      ", nn.y0_list)
    println("Numnber of hidden units:   ", nn.n_hidden)
    println("Number of training points: ", length(nn.t))
end

function pre_init_nn(f!, t, y0_list; n_hidden=10)
    Nvar = length(y0_list)
    t0 = t[1]

    function f_list(t, y)
        dy = Vector{}(length(y))
        f!(t, y, dy)
        return dy
    end

    return NNModel(Nvar, f_list, y0_list, t, t0, n_hidden,
                   nothing, nothing, nothing, nothing, nothing)
end

function reset_weights!(nn::NNModel)
    # update weights
    nn.params_list = [init_weights(n_hidden = nn.n_hidden)
                      for _ in 1:nn.Nvar]
    nn.sizes, nn.lengths, nn.total_l = sizes_and_length(nn.params_list[1])

    # update flattened weights
    nn.all_flattened_params = vcat(map(flat_opt, nn.params_list)...)
    return nothing
end

function init_nn(f!, t, y0_list; n_hidden=10)
    nn = pre_init_nn(f!, t, y0_list, n_hidden=n_hidden)
    reset_weights!(nn)
    return nn
end

function loss_func(params_list, nn::NNModel)
    # need to expose params_list to optim.jl

    # shortcut
    Nvar = nn.Nvar
    f_list = nn.f_list
    t = nn.t
    t0 = nn.t0

    y_pred_list = []
    dydt_pred_list = []

    for i = 1:Nvar
        params = params_list[i]
        y0 = y0_list[i]
        y_pred =  predict(params, t, y0, t0)
        dydt_pred = predict_dt(params, t, y0, t0)

        push!(y_pred_list, y_pred)
        push!(dydt_pred_list, dydt_pred)
    end

    f_pred_list = f_list(t, y_pred_list)

    loss_total = 0.0
    for i = 1:Nvar
        f_pred = f_pred_list[i]
        dydt_pred = dydt_pred_list[i]
        loss = mean(abs2, dydt_pred - f_pred)
        loss_total += loss
    end
    return loss_total
end

function get_unflat(all_flattened_params, nn::NNModel)
    # don't modify nn itself
    return unflatten_all(all_flattened_params,
        nn.Nvar, nn.total_l, nn.sizes, nn.lengths)
end

function train!(nn::NNModel, method; kwargs...)
    # for optim.jl
    function loss_wrap(all_flattened_params)
        params_list = get_unflat(all_flattened_params, nn)
        return loss_func(params_list, nn)
    end

    # configuration optimazation options
    od = OnceDifferentiable(loss_wrap, nn.all_flattened_params; autodiff =:forward);
    println(kwargs)
    option = Optim.Options(;store_trace=true, extended_trace=true, kwargs...)

    # training
    res = optimize(od, nn.all_flattened_params, method, option)

    # update weights
    nn.all_flattened_params = res.minimizer # flattened weights
    nn.params_list = get_unflat(res.minimizer, nn) # original weights

    return res
end

function predict(nn::NNModel; t=nothing)
    if t == nothing
        t = nn.t # predict on training points by default
    end

    # shortcut
    Nvar = nn.Nvar
    f_list = nn.f_list
    params_list = nn.params_list
    t0 = nn.t0

    y_pred_list = []
    dydt_pred_list = []

    for i = 1:Nvar
        params = params_list[i]
        y0 = y0_list[i]
        y_pred =  predict(params, t, y0, t0)
        dydt_pred = predict_dt(params, t, y0, t0)

        push!(y_pred_list, y_pred)
        push!(dydt_pred_list, dydt_pred)
    end

    return y_pred_list, dydt_pred_list
end


text:
insts = OrderedDict()

insts["test"] = [
    ((2, 3, false, VecSpecExt[]),),
    ((3, 3, true, VecSpecExt[]),),
    ((5, 5, false, [VecNegRtdet()]),),
    ((6, 4, true, [VecNegEntropy(), VecNegSqrtConj()]),),
]

function vectorregression_insts(use_extended_form::Bool, max_n::Int)
    options = (; use_extended_form = use_extended_form)
    hs = [VecNegEntropy(), VecNegSqrtConj()]
    return [((n, n, true, hs), options) for n in vcat(10, 10:10:max_n)]
end

insts["nat"] = vectorregression_insts(true, 100)
insts["noext"] = vectorregression_insts(false, 60)

return (VectorRegression, insts)


text:
"""
    Adam(;lr=0.001, gclip=0, beta1=0.9, beta2=0.999, eps=1e-8)
    update!(w,g,p::Adam)

Container for parameters of the Adam optimization algorithm used by
[`update!`](@ref).

Adam is one of the methods that compute the adaptive learning rate. It
stores accumulated gradients (first moment) and the sum of the squared
of gradients (second).  It scales the first and second moment as a
function of time. Here is the update formulas:

    m = beta1 * m + (1 - beta1) * g
    v = beta2 * v + (1 - beta2) * g .* g
    mhat = m ./ (1 - beta1 ^ t)
    vhat = v ./ (1 - beta2 ^ t)
    w = w - (lr / (sqrt(vhat) + eps)) * mhat

where `w` is the weight, `g` is the gradient of the objective function
w.r.t `w`, `lr` is the learning rate, `m` is an array with the same
size and type of `w` and holds the accumulated gradients. `v` is an
array with the same size and type of `w` and holds the sum of the
squares of the gradients. `eps` is a small constant to prevent a zero
denominator. `beta1` and `beta2` are the parameters to calculate bias
corrected first and second moments. `t` is the update count.

If `vecnorm(g) > gclip > 0`, `g` is scaled so that its norm is equal
to `gclip`.  If `gclip==0` no scaling takes place.

Reference: [Kingma, D. P., & Ba,
J. L. (2015)](https://arxiv.org/abs/1412.6980). Adam: a Method for
Stochastic Optimization. International Conference on Learning
Representations, 1–13.

"""
type Adam
    lr::AbstractFloat
    gclip::AbstractFloat
    beta1::AbstractFloat
    beta2::AbstractFloat
    eps::AbstractFloat
    t::Int
    fstm
    scndm
end

Adam(; lr=0.001, gclip=0, beta1=0.9, beta2=0.999, eps=1e-8)=Adam(lr, gclip, beta1, beta2, eps, 0, nothing, nothing)

function update!(w, g, p::Adam)
    gclip!(g, p.gclip)
    if p.fstm===nothing; p.fstm=zeros(w); p.scndm=zeros(w); end
    p.t += 1
    scale!(p.beta1, p.fstm)
    axpy!(1-p.beta1, g, p.fstm)
    scale!(p.beta2, p.scndm)
    axpy!(1-p.beta2, g .* g, p.scndm)
    fstm_corrected = p.fstm / (1 - p.beta1 ^ p.t)
    scndm_corrected = p.scndm / (1 - p.beta2 ^ p.t)
    axpy!(-p.lr, (fstm_corrected ./ (sqrt.(scndm_corrected) + p.eps)), w)
end


text:
    """
Sample a value For each parameter uniformly at random from the candidate vectors. Log-uniform sampling available by providing a log-spaced candidate vector.
"""
struct RandomSampler <: Sampler end

function (s::RandomSampler)(ho, iter)
    [list[rand(HO_RNG[threadid()], 1:length(list))] for list in ho.candidates]
end

# Latin Hypercube Sampler ======================================================

"""
Sample from a latin hypercube
"""
Base.@kwdef mutable struct LHSampler <: Sampler
    samples = zeros(0,0)
    iters = -1
end
function init!(s::LHSampler, ho)
    s.samples != zeros(0,0) && return # Already initialized
    ndims = length(ho.candidates)
    all(length(c) == ho.iterations for c in ho.candidates) ||
        throw(ArgumentError("Latin hypercube sampling requires all candidate vectors to have the same length as the number of iterations, got lengths $(repr.(collect(zip(ho.params, length.(ho.candidates))))) with $(ho.iterations) iterations"))
        s.iters == -1 && (s.iters = (1000*100*2)÷ho.iterations÷ndims)
    X, fit = LHCoptim(ho.iterations,ndims,s.iters)
    s.samples = copy(X')
end

function (s::LHSampler)(ho, iter)
    [list[s.samples[dim,iter]] for (dim,list) in enumerate(ho.candidates)]
end

"""
    CLHSampler(dims=[Continuous(), Categorical(2), ...])
Sample from a categorical/continuous latin hypercube. All continuous variables must have the same length of the candidate vectors.
"""
Base.@kwdef mutable struct CLHSampler <: Sampler
    samples = zeros(0,0)
    dims = []
end
function init!(s::CLHSampler, ho)
    s.samples != zeros(0,0) && return # Already initialized
    all(zip(s.dims, ho.candidates)) do (d,c)
        d isa Categorical || length(c) == ho.iterations
    end || throw(ArgumentError("Latin hypercube sampling requires all candidate vectors for Continuous variables to have the same length as the number of iterations, got lengths $(repr.(collect(zip(ho.params, length.(ho.candidates)))))"))
    ndims = length(ho.candidates)
    initialSample = randomLHC(ho.iterations,s.dims)
    X,_ = LHCoptim!(initialSample, 500, dims=s.dims)
    s.samples = copy(X')
end

function (s::CLHSampler)(ho, iter)
    [list[s.samples[dim,iter]] for (dim,list) in enumerate(ho.candidates)]
end

# GaussianProcesses sampler ====================================================

"""
Sample using Bayesian optimization. `GPSampler(Min)/GPSampler(Max)` fits a Gaussian process to the data and tries to use this model to figure out where the best point to sample next is (using expected improvement). Underneath, the package [BayesianOptimization.jl](https://github.com/jbrea/BayesianOptimization.jl/) is used. We try to provide reasonable defaults for the underlying model and optimizer and we do not provide any customization options for this sampler. If you want advanced control, use BayesianOptimization.jl directly.
"""
mutable struct GPSampler <: Sampler
    sense
    model
    # opt
    modeloptimizer
    logdims
    candidates
end

GPSampler() = error("You must specify GPSampler(Min)/GPSampler(Max) for minimization or maximization.")
GPSampler(sense) = GPSampler(sense,nothing,nothing,nothing,nothing)

function islogspace(x)
    all(x->x > 0, x) || return false
    length(x) > 4    || return false
    std(diff(log.(x))) < sqrt(eps(eltype(x)))
end

function to_logspace(x,logdims)
    map(x,logdims) do x,l
        l ? log10.(x) : x
    end
end

function from_logspace(x,logdims)
    map(x,logdims) do x,l
        l ? exp10.(x) : x
    end
end

function init!(s::GPSampler, ho)
    # s.model === nothing || return # Already initialized
    ndims                = length(ho.candidates)
    logdims              = islogspace.(ho.candidates)
    candidates           = to_logspace(ho.candidates, logdims)
    lower_bounds         = [minimum.(candidates)...]
    upper_bounds         = [maximum.(candidates)...]
    widths               = upper_bounds - lower_bounds
    kernel_widths        = widths/10
    log_function_noise   = 0.
    hypertuning_interval = max(9, max(ho.iterations, length(ho.history)) ÷ 9)
    model = ElasticGPE(ndims, mean = MeanConst(0.),
                       kernel = SEArd(log.(kernel_widths), log_function_noise), logNoise = 0., capacity=max(ho.iterations, length(ho.history))+1)
    # set_priors!(model.mean, [GaussianProcesses.Normal(0, 100)])

    modeloptimizer = MAPGPOptimizer(every = hypertuning_interval, noisebounds = [-4, 3], # log bounds on the function noise?
                                    maxeval = 40, kernel_bounds=fill([0,0,0], length(kernel_widths))) # max iters for optimization of the GP hyperparams

    s.model = model
    s.modeloptimizer = modeloptimizer
    s.logdims = logdims
    s.candidates = candidates
end

function train_model!(s, ho)
    xv = [reshape(to_logspace(float.(ho.history[j]), s.logdims), :, 1) for j in 1:length(ho.history)]
    input = reduce(hcat, xv)
    targets = Int(s.sense)*ho.results
    BayesianOptimization.update!(s.model, input, targets)
    BayesianOptimization.optimizemodel!(s.modeloptimizer, s.model) # This determines whether to run or not internally?
end

function Base.iterate(ho::Hyperoptimizer{<:GPSampler}, iter=1)
    iter > ho.iterations && return nothing
    s = ho.sampler
    init!(s, ho)
    th = max(4, length(ho.candidates)+1)
    if length(ho.history) < th
        samples = [rand(HO_RNG[threadid()], list) for (dim,list) in enumerate(ho.candidates)]
        nt = (; Pair.((:i, ho.params...), (iter, samples...))...)
        # res = ho.objective(iter, samples...)
        # push!(ho.results, res)
        push!(ho.history, samples)
        return nt, iter+1
    else
        try
            train_model!(s, ho)
        catch ex
            @warn("BayesianOptimization failed with optimizemodel! at iter $iter: error: ", ex, maxlog=1)
        end
    end

    # We now optimize the acq function using the random sampler. This could potentially be improved upon
    @assert length(s.model.y) >= th "Hyperoptimizer does not contain enough stored samples. If you iterate the optimizer manually, do not forget to push the results to the internal storage, see https://github.com/baggepinnen/Hyperopt.jl#full-example for an example."
    acqfunc = BayesianOptimization.acquisitionfunction(ExpectedImprovement(maximum(s.model.y)), s.model)
    # plot(s) |> display
    # plot!(x->acqfunc([x]), 1, 5, sp=2)
    # sleep(0.2)
    iters = min(80, prod(length, s.candidates))
    ho2 = Hyperoptimizer(iterations=iters, params=ho.params, candidates=s.candidates, objective=acqfunc)
    for params in ho2
        params2 = collect(params)[2:end]
        res = -Inf
        try
            res = acqfunc(params2)
            if isnan(res)
                init!(s, ho); train_model!(s, ho)
                acqfunc = BayesianOptimization.acquisitionfunction(ExpectedImprovement(maximum(s.model.y)), s.model)
            end
        catch ex
            @warn("BayesianOptimization acqfunc failed at iter $iter: error: ", ex, maxlog=1)
        end
        push!(ho2.results, res)
    end

    # @show ho2.maximizer
    # xpl = reduce(vcat,ho2.history)
    # @show length(xpl)
    # scatter!(xpl, zeros(3000), sp=2) |> display

    samples = from_logspace(ho2.maximizer, s.logdims)
    push!(ho.history, samples)
    # push!(ho.results, res)
    nt = (; Pair.((:i, ho.params...), (iter, samples...))...)
    return nt, iter+1
end

# Hyperband ====================================================================

"""

"""
Base.@kwdef mutable struct Hyperband <: Sampler
    R
    η::Int = 3
    minimum = (Inf,)
    inner = RandomSampler()
end
Hyperband(R) = Hyperband(R=R)

function hyperband_costfun(ex, params, candidates, sampler, ho_, objective)
    quote
        iters = $(esc(candidates[1]))
        if $(esc(sampler)).inner isa LHSampler
            smax = floor(Int, log($(esc(sampler)).η,$(esc(sampler)).R))
            B = (smax + 1)*$(esc(sampler)).R
            iters = floor(Int,$(esc(sampler)).R*$(esc(sampler)).η^smax)
            ss = string($(esc(sampler)).inner)
            @info "Starting Hyperband with inner sampler $(ss). Setting the number of iterations to R*η^log(η,R)=$(iters), make sure all candidate vectors have this length as well!"
        end
        costfun = $(objective)
        (::$typeof(costfun))($(esc(params[1])), $(esc(:state))) = $(esc(ex.args[2]))
        costfun, iters
    end
end

function hyperband(ho::Hyperoptimizer{Hyperband})
    hb = ho.sampler
    R,η = hb.R, hb.η
    hb.minimum = (Inf,)
    smax = floor(Int, log(η,R))
    B = (smax + 1)*R
    # ho.iterations >= R*η^smax || error("The number of iterations must be larger than R*η^log(η,R) = $(R*η^smax)")
    Juno.progress() do id
        for s in smax:-1:0
            n = ceil(Int, (B/R)*((η^s)/(s+1)))
            r = R / (η^s)
            minᵢ = successive_halving(ho, n, r, s)
            if minᵢ[1] < hb.minimum[1]
                hb.minimum = minᵢ
            end
            Base.CoreLogging.@logmsg -1 "Hyperband" progress=(smax-s)+1/(smax+1)  _id=id
        end
    end
    return hb.minimum
end

function successive_halving(ho, n, r=1, s=round(Int, log(hb.η, n)))
    hb = ho.sampler
    costfun = ho.objective
    η = hb.η
    minimum = Inf
    T = [ hb.inner(ho, i) for i=1:n ]
    # append!(ho.history, T)
    Juno.progress() do id
        for i in 0:s
            nᵢ = floor(Int,n/(η^i))
            rᵢ = floor(Int,r*(η^i))
            if i == 0
                LTend = [ costfun(rᵢ, t...) for t in T ]
            else
                LTend = [ costfun(rᵢ, t) for t in T ]
            end
            L, T = first.(LTend), last.(LTend)
            # if i == 0
            #     append!(ho.results, L)
            # else
            #     for t in eachindex(T) # Update results for those that were continued
            #         ho.results[findlast(x->x==T[t], ho.history)] = L[t]
            #     end
            # end
            append!(ho.history, T)
            append!(ho.results, L)

            perm = sortperm(L)
            besti = perm[1]
            if L[besti] < minimum[1]
                minimum = (L[besti], rᵢ, T[besti])
            end
            T = T[perm[1:floor(Int,nᵢ/η)]]
            # T, minimum
            # T, minimum = top_k(hb,T,L,nᵢ,minimum)
            Base.CoreLogging.@logmsg -1 "successive_halving" progress=i/s  _id=id
        end
    end
    return minimum
end


text:
function factors(n)
  f = [one(n)]
  for (p,e) in factor(n)
    f = reduce(vcat, f, [f*p^j for j in 1:e])
  end
  return length(f) == 1 ? [one(n), n] : sort!(f)
end


text:
using IterativeSolvers
using FactCheck

for rsvd in [IterativeSolvers.rsvdfact]
    facts(rsvd) do

        context("Small wide rectangular") do
            A = [1. 2 3 4; 5 6 7 8]
            S1 = svdfact(A)
            S2 = rsvd(A, 2, 0)

            @fact vecnorm(abs(S1[:U]) - abs(S2[:U])) --> less_than( √(eps()))
            @fact vecnorm(abs(S1[:Vt]) - abs(S2[:Vt])) --> less_than(√(eps()))
            @fact norm(S1[:S] - S2[:S]) --> less_than(√(eps()))
        end

        context("Small tall rectangular") do
            A = [1. 2; 3 4; 5 6; 7 8]
            S1 = svdfact(A)
            S2 = rsvd(A, 2, 0)

            @fact vecnorm(abs(S1[:U]) - abs(S2[:U])) --> less_than(√(eps()))
            @fact vecnorm(abs(S1[:Vt]) - abs(S2[:Vt])) --> less_than(√(eps()))
            @fact norm(S1[:S] - S2[:S]) --> less_than(√(eps()))
        end

        context("Small square") do
            A = [1. 2 3; 4 5 6; 7 8 9]
            S1 = svdfact(A)
            S2 = rsvd(A, 3, 0)

            @fact vecnorm(abs(S1[:U]) - abs(S2[:U]))  --> less_than(√(eps()))
            @fact vecnorm(abs(S1[:Vt]) - abs(S2[:Vt]))  --> less_than(√(eps()))
            @fact norm(S1[:S] - S2[:S])  --> less_than(√(eps()))
        end

        context("Low rank") do #Issue #42
            n = 10
            r = 2
            A = randn(n, r)*randn(r, n)
            S = svdvals(A)
            for nvals = 1:r
                S1= IterativeSolvers.rsvdvals(A, nvals, r-nvals)
                for i = 1:nvals
                    @fact abs(S[i] - S1[i]) --> less_than(n^2*r*eps())
                end
            end
        end

        context("rrange_adaptive") do
            A = [1. 2 3; 4 5 6; 7 8 9]
            @fact size(IterativeSolvers.rrange_adaptive(A, 3, 1e-3)) --> (3,2)
        end
    end
end


text:
export DataNeuron
type DataNeuron <: Neuron
    value :: Float32
    ∇     :: Float32
end

function DataNeuron(value::AbstractFloat)
    DataNeuron(value, 0.0)
end

export WeightedNeuron
@neuron type WeightedNeuron
    weights  :: DenseArray{Float32}
    ∇weights :: DenseArray{Float32}

    bias     :: DenseArray{Float32}
    ∇bias    :: DenseArray{Float32}
end

@neuron forward(neuron::WeightedNeuron) do
    for i in 1:length(neuron.inputs[1])
        neuron.value += neuron.weights[i] * neuron.inputs[1][i]
    end
    neuron.value += neuron.bias[1]
end

@neuron backward(neuron::WeightedNeuron) do
    # Only backpropogate if input layer requests it
    # Length 0 indicates no buffer was created to receive backpropogated ∇s
    # if length(neuron.∇inputs[1]) > 0
        for i in 1:length(neuron.inputs[1])
            neuron.∇inputs[1][i] += neuron.weights[i] * neuron.∇
        end
    # end
    for i in 1:length(neuron.inputs[1])
        neuron.∇weights[i] += neuron.inputs[1][i] * neuron.∇
    end
    neuron.∇bias[1] += neuron.∇
end


text:
using Test
@testset "JUDI Unit tests" begin
    include("test_abstract_vectors.jl")
    include("test_geometry.jl")
    include("test_judiVector.jl")
    include("test_linear_operators.jl")
end


text:
function solve(Op::BEMOperator, b)
    A = matrix(Op)
    u = A \ b
    Expansion(src(A), u)
end


text:
for N in [Float64, Rational{Int}, Float32]
    # random line segment
    rand(LineSegment)

    # construction
    p, q = N[1, 1], N[2, 2]
    l = LineSegment(p, q)

    # dimension
    @test dim(l) == 2

    # support vector
    @test σ(N[1, 1], l) == q
    @test σ(N[0, 1], l) == q
    @test σ(N[-1, 1], l) == q
    @test σ(N[-1, 0], l) == p
    @test σ(N[-1, -1], l) == p
    @test σ(N[0, -1], l) == p
    @test σ(N[1, -1], l) == q
    @test σ(N[1, 0], l) == q

    # boundedness
    @test isbounded(l)

    # membership
    @test N[0, 0] ∉ l
    @test N[1, 1] ∈ l
    @test N[1.5, 1.5] ∈ l
    @test N[3, 4] ∉ l
    @test N[7, 4] ∉ l
    @test N[1.5, 1.6] ∉ l

    # approximate membership test
    if N == Float64
        @test N[1.5, 1.5] ∈ LineSegment(N[1.5, 1.50000000000001], N[1.5, 2.0])
        r = LazySets._rtol(N)
        LazySets.set_rtol(N, 1e-20)
        @test !(N[1.5, 1.5] ∈ LineSegment(N[1.5, 1.50000000000001], N[1.5, 2.0]))
        LazySets.set_rtol(N, r)
    end

    # center/generators
    @test center(l) == N[1.5, 1.5]
    @test collect(generators(l)) ∈ [[N[1/2, 1/2]], [N[-1/2, -1/2]]]
    @test genmat(l) ∈ [hcat(N[1/2 1/2]), hcat(N[-1/2, -1/2])]
    @test ngens(l) == 1
    l_degenerate = LineSegment(p, p)
    @test collect(generators(l_degenerate)) == Vector{N}()
    @test genmat(l_degenerate) == Matrix{N}(undef, 2, 0)
    @test ngens(l_degenerate) == 0

    # isempty
    @test !isempty(l)

    # isuniversal
    answer, w = isuniversal(l, true)
    @test !isuniversal(l) && !answer && w ∉ l

    # an_element function
    @test an_element(l) ∈ l

    # vertices_list function
    vl = vertices_list(l)
    @test ispermutation(vl, [l.p, l.q])

    # translation
    @test translate(l, N[1, 2]) == LineSegment(N[2, 3], N[3, 4])

    # intersection emptiness
    l1 = LineSegment(N[1, 1], N[2, 2])
    l2 = LineSegment(N[2, 1], N[1, 2])
    l3 = LineSegment(N[0, 1], N[0, 2])
    l4 = LineSegment(N[1, 1], N[1, 1])
    l5 = LineSegment(N[0, 0], N[0, 0])
    l6 = LineSegment(N[13//10, 13//10], N[23//10, 23//10])
    l7 = LineSegment(N[3, 3], N[4, 4])
    l8 = LineSegment(N[1, 2], N[2, 3])
    l1_copy = LineSegment(copy(l1.p), copy(l1.q))
    intersection_empty, point = is_intersection_empty(l1, l1, true)
    @test !is_intersection_empty(l1, l1) && !intersection_empty && point ∈ l1
    intersection_empty, point = is_intersection_empty(l1, l2, true)
    @test !is_intersection_empty(l1, l2) && !intersection_empty && point ∈ l1 && point ∈ l2
    @test is_intersection_empty(l1, l3) && is_intersection_empty(l1, l3, true)[1]
    intersection_empty, point = is_intersection_empty(l1, l4, true)
    @test !is_intersection_empty(l1, l4) && !intersection_empty && point ∈ l1 && point ∈ l4
    @test is_intersection_empty(l1, l5) && is_intersection_empty(l1, l5, true)[1]
    intersection_empty, point = is_intersection_empty(l1, l6, true)
    @test !is_intersection_empty(l1, l6) && !intersection_empty && point ∈ l1 && point ∈ l6
    intersection_empty, point = is_intersection_empty(l2, l2, true)
    @test !is_intersection_empty(l2, l2) && !intersection_empty && point ∈ l2
    @test is_intersection_empty(l2, l3) && is_intersection_empty(l2, l3, true)[1]
    @test is_intersection_empty(l2, l4) && is_intersection_empty(l2, l4, true)[1]
    @test is_intersection_empty(l2, l5) && is_intersection_empty(l2, l5, true)[1]
    intersection_empty, point = is_intersection_empty(l2, l6, true)
    @test !is_intersection_empty(l2, l6) && !intersection_empty && point ∈ l2 && point ∈ l6
    intersection_empty, point = is_intersection_empty(l3, l3, true)
    @test !is_intersection_empty(l3, l3) && !intersection_empty && point ∈ l3
    @test is_intersection_empty(l3, l4) && is_intersection_empty(l3, l4, true)[1]
    @test is_intersection_empty(l3, l5) && is_intersection_empty(l3, l5, true)[1]
    @test is_intersection_empty(l3, l6) && is_intersection_empty(l3, l6, true)[1]
    intersection_empty, point = is_intersection_empty(l4, l4, true)
    @test !is_intersection_empty(l4, l4) && !intersection_empty && point ∈ l4
    @test is_intersection_empty(l4, l5) && is_intersection_empty(l4, l5, true)[1]
    @test is_intersection_empty(l4, l6) && is_intersection_empty(l4, l6, true)[1]
    intersection_empty, point = is_intersection_empty(l5, l5, true)
    @test !is_intersection_empty(l5, l5) && !intersection_empty && point ∈ l5
    @test is_intersection_empty(l5, l6) && is_intersection_empty(l5, l6, true)[1]
    intersection_empty, point = is_intersection_empty(l6, l6, true)
    @test !is_intersection_empty(l6, l6) && !intersection_empty && point ∈ l6
    @test is_intersection_empty(l1, l7) && is_intersection_empty(l1, l7, true)[1]
    @test is_intersection_empty(l1, l8) && is_intersection_empty(l1, l8, true)[1]
    intersection_empty, point = is_intersection_empty(l1, l1_copy, true)
    @test !is_intersection_empty(l1, l1_copy) && !intersection_empty && point ∈ l1

    # intersection
    s1 = LineSegment(N[-5.0,-5.0], N[5.0,5.0])
    # colinear shifted down
    s2 = LineSegment(N[-6.0,-6.0], N[4.0,4.0])
    @test intersection(s1, s2) isa LineSegment
    @test isapprox(intersection(s1, s2).p, N[-5,-5])
    @test isapprox(intersection(s1, s2).q, N[4,4])
    # parallel, not intersect
    s3 = LineSegment(N[-5.0,-4.0], N[4.0,5.0])
    @test intersection(s1, s3) isa EmptySet
    # intersect outside of segment
    s4 = LineSegment(N[0.0,10.0], N[6.0,5.0])
    @test intersection(s1, s4) isa EmptySet
    # intersect in segment
    s5 = LineSegment(N[5.0,-5.0], N[-5.0,5.0])
    @test intersection(s1, s5) isa Singleton
    @test isapprox(intersection(s1,s5).element, N[0,0])
    # parallel, no points in common
    s6 = LineSegment(N[10.0, 10.0], N[11.0,11.0])
    @test intersection(s1, s6) isa EmptySet
    # parallel one point in common
    s7 = LineSegment(N[5.0,5.0], N[6.0,6.0])
    @test intersection(s1, s7) isa Singleton
    @test isapprox(intersection(s1,s7).element, N[5,5])
    s8 = LineSegment(N[0.0,0.0], N[1.0,0.0])
    s9 = LineSegment(N[0.0,0.0], N[2.0,0.0])
    @test intersection(s9, s8) isa LineSegment
    @test isapprox(intersection(s9, s8).p, N[0,0])
    @test isapprox(intersection(s9, s8).q, N[1,0])

    # subset
    l = LineSegment(N[1, 1], N[2, 2])
    b1 = Ball1(N[1.5, 1.5], N(1.1))
    b2 = Ball1(N[1.5, 1.5], N(0.4))
    subset, point = ⊆(l, b1, true)
    @test l ⊆ b1 && subset && point == N[]
    subset, point = ⊆(l, b2, true)
    @test l ⊈ b2 && !subset && point ∈ l && point ∉ b2
    h1 = Hyperrectangle(N[1.5, 1.5], N[0.6, 0.8])
    h2 = Hyperrectangle(N[1.5, 1.5], N[0.4, 0.8])
    subset, point = ⊆(l, h1, true)
    @test l ⊆ h1 && subset && point == N[]
    subset, point = ⊆(l, h2, true)
    @test l ⊈ h2 && !subset && point ∈ l && point ∉ h2

    # halfspace_left & halfspace_right
    @test N[1, 2] ∈ halfspace_left(l)
    @test N[2, 1] ∈ halfspace_right(l)

    # constraints list
    l = LineSegment(N[0, 0], N[1, 1])
    clist = constraints_list(l)
    @test ispermutation(clist, [HalfSpace(N[1, -1], N(0)),  # x <= y
                                HalfSpace(N[-1, 1], N(0)),  # x >= y
                                HalfSpace(N[-1, -1], N(0)), # y >= -x
                                HalfSpace(N[1, 1], N(2))])  # y <= 2-x

    # conversion
    z = convert(Zonotope, l)
    G = hcat(N[1//2, 1//2])
    @test center(z) == N[1//2, 1//2] && genmat(z) ∈ [G, -G]
end


text:
#
# from and to_list for Sequence
#

a = [1, 2, 3, 4, 5]
t = from(a)
@test consume(t.chain) == Subject(1)
@test consume(t.chain) == Subject(2)
@test consume(t.chain) == Subject(3)
@test consume(t.chain) == Subject(4)
@test consume(t.chain) == Subject(5)
@test consume(t.chain) == :done

#
# map, filter and reduce
#

t = from(a)(map_, x -> x*2)
@test consume(t.chain) == Subject(2)
@test consume(t.chain) == Subject(4)
@test consume(t.chain) == Subject(6)
@test consume(t.chain) == Subject(8)
@test consume(t.chain) == Subject(10)

# (reduce, (x, y) -> x + y)

#
# Assembly line
#

a = [1, 2, 3, 4, 5]
t = from(a)(map_, x -> x+5)(map_, x -> x*2)
@test to_list(t) == [12, 14, 16, 18, 20]


text:
module TestDynaPlanning_Instances

include("../src/DynaPlanning.jl")
include("../src/DynaUtils.jl")

using LightGraphs, LightGraphsFlows, GraphPlot, Compose
using SparseArrays
using DataStructures
using Random
using Test
using Logging

import .DynaPlanning
import .DynaUtils

# Logging
logger = LogLevel(Logging.Debug)

mutable struct SimulationParameters
    episodes::Int64
end

parameters = SimulationParameters(100)

################################################################################
# Pigou's example
################################################################################

sources = [1, 2, 1, 4]
sinks = [2, 3, 4, 3]

Pigou = SimpleDiGraph(length(sources))

for i = 1:length(sources)
    add_edge!(Pigou, sources[i], sinks[i])
end

a_e = [0.0, 0.0, 1.0, 0.0]
b_e = [1.0, 0.0, 0.0, 0.0]

# Vector of linear cost functions
c_e = Array{Function}(undef, 4)
for i = 1:length(sources)
    c_e[i] = x -> a_e[i] * x + b_e[i]
end

################################################################################
# # Create SimpleWeightedDiGraph from a Barabási–Albert model random graph

G = barabasi_albert(12, 2, 2, is_directed=true, seed=23)
A = adjacency_matrix(G)

# Cost function
rng = MersenneTwister(1234);
a_e_G = rand!(rng, zeros(12))

rng = MersenneTwister(5678);
b_e_G = rand!(rng, zeros(12))

locs_x, locs_y = spring_layout(G)

@info("Plotting Barabási–Albert network (G)")
DynaUtils.plot_graph(G, locs_x, locs_y)

# DynaAgents
N_DynaAgents = 4

@info "Setting Dyna_agentPlus and StackelbergInstance"  N_DynaAgents

srcDst = Array{Tuple,1}()
push!(srcDst, (11, 1))
push!(srcDst, (12, 1))
push!(srcDst, (12, 1))
push!(srcDst, (10, 1))

agents = Array{DynaPlanning.DynaAgentPlus,1}()
for n in range(1, stop=N_DynaAgents)
    testAgent = DynaPlanning.DynaAgentPlus()
    testSrcDst = pop!(srcDst)

    # A* shortest path
    testShortestPath = a_star(SimpleDiGraph(A), testSrcDst[1], testSrcDst[2])

    testAgent.path = Queue{Edge}()
    for i in 1:length(testShortestPath)
        enqueue!(testAgent.path, testShortestPath[i])
    end

    push!(agents, testAgent)
end

# Stackelberg instance
testStackelbergInstance = DynaPlanning.StackelbergInstance()
testStackelbergInstance.network = G
testStackelbergInstance.agents = agents

end # module

max_stars_repo_name: mgoar/dynaplanning

max_issues_repo_name: mgoar/dynaplanning

max_forks_repo_name: mgoar/dynaplanning


text:
"""
    weights, wids, nvrtx = signed_volume(simplex, wids, nvrtx)

This function implements the Signed Volume distance sub-algorithm from
Montanari, Petrinijc: Improving the GJK Algorithm for faster and more
reliable distance queries between convex objects.  The advantage is to avoid
Johnson's sub-algorithm numerical instability in the Ax=b system for calculating
the next v_k vector.

M is number of points in simplex, and SVector{N,T} is size 3 and type Float64.
In the main algorithm, M is always fully populated (size = 4), and the weights
are set to zero for non-active simplex points. wids tracks the top nvrtx active vertices.
"""
function signed_volume(simplex::SVector{M,SVector{N,T}}, wids::SVector{M,P}, nvrtx::P) where {M,N,T,P}
    # Calculate weights and mininum active support set
    if nvrtx == 4
        weights, wids, nvrtx = S3D(simplex, wids)
        return weights, wids, nvrtx
    elseif nvrtx == 3
        weights, wids, nvrtx = S2D(simplex, wids)
        return weights, wids, nvrtx
    elseif nvrtx == 2
        weights, wids, nvrtx = S1D(simplex, wids)
        return weights, wids, nvrtx
    elseif nvrtx == 1
        weights = SVector{M,T}(1, 0, 0, 0)
        weights = order_weights(weights, wids)
        return weights, wids, nvrtx
    end
end

"""
S3D function for determining minimum active support set for tetrahedron and associated weights.
"""
function S3D(simplex::SVector{M,SVector{N,T}}, wids::SVector{M,P}) where {M,N,T,P}
    # Unpack simplex in wids order
    _s = view(simplex,1:4)

    # Calculate signed determinants to compare tetrahedron orientations
    M_41 = SMatrix{3,3,T}(_s[2][1],_s[2][2],_s[2][3],_s[3][1],_s[3][2],_s[3][3],_s[4][1],_s[4][2],_s[4][3])
    M_42 = SMatrix{3,3,T}(_s[1][1],_s[1][2],_s[1][3],_s[3][1],_s[3][2],_s[3][3],_s[4][1],_s[4][2],_s[4][3])
    M_43 = SMatrix{3,3,T}(_s[1][1],_s[1][2],_s[1][3],_s[2][1],_s[2][2],_s[2][3],_s[4][1],_s[4][2],_s[4][3])
    M_44 = SMatrix{3,3,T}(_s[1][1],_s[1][2],_s[1][3],_s[2][1],_s[2][2],_s[2][3],_s[3][1],_s[3][2],_s[3][3])
    C = SVector{4,T}(-det(M_41), det(M_42), -det(M_43), det(M_44))
    detM = sum(C)
    facets_test = compare_signs.(detM, C)

    # Calculate minimum support set and associated weights
    if facets_test[1] == facets_test[2] == facets_test[3] == facets_test[4] == true
        weights = C / detM
        new_wids = SVector{4,Int}(1, 2, 3, 4)
        return weights, new_wids, 4
    else
        d_min = Inf
        best_weights, best_wids, best_nvrtx = undef, undef, undef
        tmp_weights, tmp_wids, tmp_nvrtx = undef, undef, undef
        @inbounds for i = 1:4
            if facets_test[i] == false  # wrt s ordering
                if i == 1
                    new_wids = SVector{4,Int}(2, 3, 4, 1)
                elseif i == 2
                    new_wids = SVector{4,Int}(1, 3, 4, 2)
                elseif i == 3
                    new_wids = SVector{4,Int}(1, 2, 4, 3)
                elseif i == 4
                    new_wids = SVector{4,Int}(1, 2, 3, 4)
                end
                tmp_weights, tmp_wids, tmp_nvrtx = S2D(simplex, new_wids)  # weights are in order
                v = linear_combination(tmp_weights, simplex)
                v_len = dot(v, v)
                if v_len < d_min
                    best_weights, d_min, best_wids, best_nvrtx = tmp_weights, v_len, tmp_wids, tmp_nvrtx
                end
            end
        end
        return best_weights, best_wids, best_nvrtx
    end
end

"""
S2D function for determining minimum active support set for plane (triangle) and associated weights.
Modification to the original algorithm by rotating the plane to a flat orientation.
"""
function S2D(simplex::SVector{M,SVector{N,T}}, wids::SVector{M,P}) where {M,N,T,P}
    # Unpack indices
    _wids = view(wids,1:4)

    # Unpack simplex in wids order
    _s = view(simplex[wids],1:3)

    # Setup normals and origin projection
    AB = _s[2] - _s[1]
    AC = _s[3] - _s[1]
    n = cross(AB, AC)
    p = (dot(_s[1], n)*n) / dot(n, n)
    _p = view(p,1:3)

    # Calculated signed determinants
    # ABC
    xyABC = SMatrix{3,3,T}(_s[1][1], _s[1][2], 1, _s[2][1], _s[2][2], 1, _s[3][1], _s[3][2], 1)
    yzABC = SMatrix{3,3,T}(_s[1][2], _s[1][3], 1, _s[2][2], _s[2][3], 1, _s[3][2], _s[3][3], 1)
    xzABC = SMatrix{3,3,T}(_s[1][1], _s[1][3], 1, _s[2][1], _s[2][3], 1, _s[3][1], _s[3][3], 1)
    signed_area = SVector{3,T}(det(yzABC), det(xzABC), det(xyABC))
    _, I3 = findmax(abs.(signed_area))  # index to exclude
    if I3 == 1
        I1, I2 = 2, 3
    elseif I3 == 2
        I1, I2 = 1, 3
    elseif I3 == 3
        I1, I2 = 1, 2
    end
    # ABC
    nu_max = signed_area[I3]
    # ABP vs. ABC
    ABP = SMatrix{3,3,T}(_s[1][I1], _s[1][I2], 1, _s[2][I1], _s[2][I2], 1, _p[I1], _p[I2], 1)
    # BCP vs. BCA
    BCP = SMatrix{3,3,T}(_s[2][I1], _s[2][I2], 1, _s[3][I1], _s[3][I2], 1, _p[I1], _p[I2], 1)
    # ACP vs. ACB
    ACP = SMatrix{3,3,T}(_s[1][I1], _s[1][I2], 1, _s[3][I1], _s[3][I2], 1, _p[I1], _p[I2], 1)
    C = SVector{4,T}(det(BCP), -det(ACP), det(ABP), 0)  # wrt s ordering
    facets_test = compare_signs.(nu_max, C)

    # Calculate minimum support set and associated weights
    if facets_test[1] == facets_test[2] == facets_test[3] == true
        weights = C / nu_max
        weights = order_weights(weights, wids)
        return weights, wids, 3
    else
        d_min = Inf
        best_weights, best_wids, best_nvrtx = undef, undef, undef
        tmp_weights, tmp_wids, tmp_nvrtx = undef, undef, undef
        @inbounds for i = 1:3
            if facets_test[i] == false
                if i == 1
                    new_wids = SVector{4,P}(_wids[2], _wids[3], _wids[4], _wids[1])
                elseif i == 2
                    new_wids = SVector{4,P}(_wids[1], _wids[3], _wids[4], _wids[2])
                elseif i == 3
                    new_wids = SVector{4,P}(_wids[1], _wids[2], _wids[4], _wids[3])
                end
                tmp_weights, tmp_wids, tmp_nvrtx = S1D(simplex, new_wids)  # weights are in order
                v = linear_combination(tmp_weights, simplex)
                v_len = dot(v, v)
                if v_len < d_min
                    best_weights, d_min, best_wids, best_nvrtx = tmp_weights, v_len, tmp_wids, tmp_nvrtx
                end
            end
        end
        return best_weights, best_wids, best_nvrtx
    end
end

"""
S1D function for determining minimum active support set for line and associated weights.
Modification to the original algorithm with the calculation of the weights.
"""
function S1D(simplex::SVector{M,SVector{N,T}}, wids::SVector{M,P}) where {M,N,T,P}

    # Unpack indices
    _wids = view(wids,1:2)

    # Unpack simplex in wids order
    _s = view(simplex[_wids],:)

    # Calculate projection
    AB = _s[2] - _s[1]
    t = - dot(_s[1], AB) / dot(AB, AB)  # s(t) = (1-t)*s1 + t*s2

    if t < 0
        weights = SVector{4,T}(1, 0, 0, 0)
        weights = order_weights(weights, wids)
        return weights, wids, 1
    elseif t < 1
        weights = SVector{4,T}(1-t, t, 0, 0)
        weights = order_weights(weights, wids)
        return weights, wids, 2
    else
        weights = SVector{4,T}(0, 1, 0, 0)
        weights = order_weights(weights, wids)
        return weights, wids, 1
    end
end

"""
Helper function for sign comparison used in orientation determination.
"""
function compare_signs(a::T, b::T) where {T}
    if (a>0) && (b>0) || (a<0) && (b<0)
        return true
    else
        return false
    end
end

"""
Helper function to order weights based on wids.
"""
function order_weights(weights::SVector{M,T}, wids::SVector{M,P}) where {M,T,P}
    ordered_weights = MVector{M,T}(undef)
    ordered_weights[wids] = weights
    return SVector(ordered_weights)
end

num_johnson_subsets(simplex_length::Integer) = 2^simplex_length - 1

"""
Compute all subsets of the points in the simplex in a reliable order. The
order is arbitrary, but was chosen to match the implemenation in
S. Cameron, “Enhancing GJK: computing mininum and penetration distances between
convex polyhedra,”. Specifically, subset i contains point j iff the binary
representation of i has a one at the jth bit.
"""
function johnson_subsets(simplex_length::Integer)
    num_subsets = num_johnson_subsets(simplex_length)
    subsets = falses(simplex_length, num_subsets)
    for i in 1:num_subsets
        for j in 1:simplex_length
            subsets[j, i] = (i >> (j - 1)) & 1
        end
    end
    subsets
end

"""
    weights = projection_weights(simplex)

This function implements Johnson's distance subalgorithm, as described in
E. G. Gilbert, D. W. Johnson, and S. S. Keerthi, “A fast procedure for
computing the distance between complex objects in three-dimensional space,”
1988. Given a simplex (a length N+1 vector of points of dimension N), it
returns weights such that LinearAlgebra.dot(weights, simplex) yields the point in the convex
hull of the simplex which is closest to the origin.

This is the critical loop of the GJK algorithm, so it has been heavily optimized
to precompute, inline, and unroll as nuch of the algorithm as possible. For a
more readable (and nuch slower) implementation, see
projection_weights_reference()
"""
@generated function projection_weights(simplex::SVector{M, SVector{N, T}}; atol=sqrt(eps(T))) where {M,N,T}
    simplex_length = M
    num_subsets = num_johnson_subsets(M)
    subsets = johnson_subsets(M)
    complements = .!subsets

    expr = quote
        deltas = SVector(tuple($([zero(SVector{simplex_length, T}) for i = 1 : num_subsets]...)))
    end

    # Set the weight of every singleton subset to 1.
    for i in 1:simplex_length
        elements = [:(zero($T)) for j in 1:simplex_length]
        elements[i] = :(one($T))
        arg = :(SVector{$simplex_length, $T}(tuple($(elements...))))

        push!(expr.args, quote
            deltas = setindex(deltas, $(arg), $(2^(i - 1)))
        end)
    end

    for s in 1:(num_subsets - 1)
        k = first((1:M)[subsets[:,s]])
        push!(expr.args, quote
            viable = true
        end)
        for j in (1:M)[complements[:,s]]
            s2 = s + (1 << (j - 1))
            push!(expr.args, quote
                d = deltas[$s2][$j]
            end)
            for i in (1:M)[subsets[:,s]]
                push!(expr.args, quote
                    d += deltas[$s][$i] *
                    (LinearAlgebra.dot(simplex[$i], simplex[$k]) - LinearAlgebra.dot(simplex[$i], simplex[$j]))
                end)
            end
            push!(expr.args, quote
                deltas = setindex(deltas, setindex(deltas[$s2], d, $j), $s2)
            end)
            push!(expr.args, quote
                if d > atol # See #17
                    viable = false
                end
            end)
        end
        push!(expr.args, quote
            viable = viable && all($(Expr(:tuple,
            [:(deltas[$s][$i] >= 0) for i in (1:M)[subsets[:,s]]]...)))
            if viable
                return deltas[$s] ./ sum(deltas[$s])
            end
        end)
    end
    push!(expr.args, quote
        return deltas[$num_subsets] ./ sum(deltas[$num_subsets])
    end)
    return expr
end

function projection_weights_reference(simplex::SVector{M, SVector{N, T}}) where {M,N,T}
    num_subsets = 2^M - 1
    complements = falses(M, num_subsets)
    subsets = falses(M, num_subsets)
    for i in 1:num_subsets
        for j in 1:M
            subsets[j, i] = (i >> (j - 1)) & 1
            complements[j, i] = !subsets[j, i]
        end
    end

    deltas = zeros(MMatrix{M, num_subsets, T})
    viable = ones(MVector{num_subsets, Bool})

    for i in 1:M
        deltas[i, 2^(i - 1)] = 1
    end

    for s in 1:(num_subsets - 1)
        k = first((1:M)[subsets[:,s]])
        for j in (1:M)[complements[:,s]]
            s2 = s + (1 << (j - 1))
            delta_j_s2 = 0
            for i in (1:M)[subsets[:,s]]
                delta_j_s2 += deltas[i, s] * (dot(simplex[i], simplex[k]) - dot(simplex[i], simplex[j]))
            end
            if delta_j_s2 > 0
                viable[s] = false
            elseif delta_j_s2 < 0
                viable[s2] = false
            end
            deltas[j, s2] = delta_j_s2
        end
        if viable[s]
            return deltas[:,s] ./ sum(deltas[:,s])
        end
    end
    return deltas[:,end] ./ sum(deltas[:,end])
end


text:
using BEAST, CompScienceMeshes, MieSeries

Θ = range(0, stop=π, length=100)
Φ = 0.0
P = [ [cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)] for θ in Θ for ϕ in Φ]

## Solve the scatering problem using BEM
Γ = readmesh(joinpath(dirname(@__FILE__),"sphere2.in"))
X = raviartthomas(Γ)

ϵ = μ = 1.0
ω = 1.0; κ = ω * √(ϵ*μ)
t = Maxwell3D.singlelayer(wavenumber=κ)
E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
e = (n × E) × n

@hilbertspace j; @hilbertspace k
efie = @discretise t[k,j]==e[k]  j∈X k∈X
u = gmres(efie)

A = Maxwell3D.farfield(wavenumber=κ)
farbem = 1/(4pi)* norm.(potential(A, P, u, X))

## Solve the scattering problem by computing the Mie series
radius = 1.0; ϵ = μ = 1.0
sphere = PECSphere(radius, ϵ, μ)
exc = excite(sphere, ω, :x)
farmie = [norm(collect(farfield(sphere, exc, p)[1:3])) for p in P]

# Compare
@show norm(farmie - farbem) / norm(farmie)


text:
"""
    mutations_per_cell(population)

calculate the number of mutations per cell in each module
"""
mutations_per_cell(population::MultiSimulation) = map(mutations_per_cell, population)

"""
    mutations_per_cell(simulation::Simulation)

calculate the number of mutations per cell
"""
mutations_per_cell(simulation::Simulation) = mutations_per_cell(simulation.output)

"""
    mutations_per_cell(moduletracker::ModuleTracker)
"""
mutations_per_cell(moduletracker::ModuleTracker) = mutations_per_cell(moduletracker.cells)

mutations_per_cell(cells::Array{Cell, 1}) = map(cell -> length(cell.mutations), cells)

mutation_ids_by_cell(moduletracker::ModuleTracker, idx=nothing) = mutation_ids_by_cell(moduletracker.cells, idx)

function mutation_ids_by_cell(cells::Array{Cell, 1}, idx=nothing)
    if isnothing(idx)
        return map(cell -> cell.mutations, cells)
    else
        return map(cell --> cell.mutations, cells[idx])
    end
end
"""
    average_mutations_per_module(population)

calculate the mean number of mutations in each module
"""
function average_mutations_per_module(population)
    return map(average_mutations, population)
end

"""
    average_mutations(population, var=false)

calculate the mean number of mutations in whole population (and variance if variance is true)
"""
function average_mutations(population, variance=false)
    mutations = (length(cell.mutations)
        for moduletracker in population for cell in moduletracker.cells)
    if variance
        return mean(mutations), var(mutations)
    else
        return mean(mutations)
    end
end

"""
    average_mutations(simulation::Simulation)
"""
average_mutations(simulation::Simulation) = average_mutations(simulation.output)
"""
    average_mutations(moduletracker::ModuleTracker)
"""
function average_mutations(moduletracker::ModuleTracker)
    return mean(length(cell.mutations) for cell in moduletracker.cells)
end

"""
    clonal_mutations(simulation::Simulation)

returns the number of clonal mutations, i.e. mutations present in every cell
"""
function clonal_mutations(simulation::Simulation)
    return clonal_mutations(simulation.output)
end

"""
    clonal_mutations(multisim::MultiSimulation)

returns the number of clonal mutations in each module/simulated population
"""
function clonal_mutations(population::MultiSimulation)
    return map(clonal_mutations, population)
end

"""
    clonal_mutations(output::SimulationResult)
"""
function clonal_mutations(moduletracker::ModuleTracker)
    return length(clonal_mutation_ids(moduletracker))
end

"""
    clonal_mutation_ids(population, idx=nothing)
"""
function clonal_mutation_ids(population, idx=nothing)
    if isnothing(idx)
        return map(clonal_mutation_ids, population)
    else
        return map(clonal_mutation_ids, population[idx])
    end
end

"""
    clonal_mutation_ids(moduletracker)
"""
function clonal_mutation_ids(moduletracker::ModuleTracker)
    return intersect([cell.mutations for cell in moduletracker.cells]...)
end

"""
    pairwise_fixed_differences(simulation::Simulation[, idx])

Calculate the number of pairwise fixed differences between every pair of cells and return
as a dictionary (number differneces => frequency). If `idx` is given, only include listed
cells.
"""

function pairwise_fixed_differences(simulation::Simulation, idx=nothing)
    return pairwise_fixed_differences(simulation.output, idx)
end

"""
    pairwise_fixed_differences(moduletracker::ModuleTracker[, idx])

See pairwise_fixed_differences(simulation::Simulation)
"""

function pairwise_fixed_differences(moduletracker::ModuleTracker, idx=nothing)
    muts = mutation_ids_by_cell(moduletracker, idx)
    return pairwise_fixed_differences(muts)
end
"""
    pairwise_fixed_differences(population[, idx])

Calculate the number of pairwise fixed differences between every pair of modules and return
as a dictionary (number differneces => frequency). If `idx` is given, only include listed
modules.
"""

function pairwise_fixed_differences(population, idx=nothing)
    clonalmuts = clonal_mutation_ids(population, idx)
    return pairwise_fixed_differences(clonalmuts)
end

function pairwise_fixed_differences(muts::Vector{Vector{Int64}})
    n = length(muts)
    pfd_vec = Int64[]
    for i in 1:n
        for j in i+1:n
            push!(pfd_vec, length(symdiff(muts[i], muts[j])))
        end
    end
    return countmap(pfd_vec)
end

"""
    pairwise_fixed_differences_matrix(population[, idx], diagonals=false)

Calculate the number of pairwise fixed differences between modules. Return an n x n matrix,
such that the value at (i,j) is the pairwise fixed differences between modules i and j, and
there are n total modules.

If idx is given only compare specified modules, otherwise compare all modules in the
population. If diagonals is true include comparison with self (i.e. number of fixed clonal
mutations in each module).
"""
function pairwise_fixed_differences_matrix(population, idx=nothing; diagonals=false)
    clonalmuts = clonal_mutation_ids(population, idx)
    return pairwise_fixed_differences_matrix(clonalmuts, diagonals=diagonals)
end

function pairwise_fixed_differences_matrix(simulation::Simulation, idx=nothing; diagonals=false)
    return pairwise_fixed_differences_matrix(simulation.output, idx, diagonals=diafgonals)
end

function pairwise_fixed_differences_matrix(moduletracker::ModuleTracker, idx=nothing; diagonals=false)
    muts = mutation_ids_by_cell(moduletracker, idx)
    return pairwise_fixed_differences_matrix(muts, diagonals=diagonals)
end

function pairwise_fixed_differences_matrix(muts::Vector{Vector{Int64}}; diagonals=false)
    n = length(muts)
    pfd = zeros(Int64, n, n)
    for i in 1:n
        if diagonals pfd[i,i] = length(muts[i]) end
        for j in i+1:n
            pfd[j,i] = length(symdiff(muts[i], muts[j]))
        end
    end
    return pfd
end

"""
    pairwise_fixed_differences_statistics(population[, idx], clonal=false)

Calculate the mean and variance of the number of pairwise fixed differences between modules.
If idx is given only compare specified modules, otherwise compare all modules in the
population. If clonal is true also calculate mean and variance of number of clonal mutations
in each module.
"""
function pairwise_fixed_differences_statistics(population, idx=nothing; clonal=true)
    clonalmuts = clonal_mutation_ids(population, idx)
    return pairwise_fixed_differences_statistics(clonalmuts; clonal=clonal)
end

function pairwise_fixed_differences_statistics(clonalmuts::Vector{Vector{Int64}}; clonal=true)
    n = length(clonalmuts)
    pfd = Int64[]
    for i in 1:n
        for j in i+1:n
            push!(pfd, length(symdiff(clonalmuts[i], clonalmuts[j])))
        end
    end
    if clonal
        nclonalmuts = (length(cmut) for cmut in clonalmuts)
        return mean(pfd), var(pfd), mean(nclonalmuts), var(nclonalmuts)
    else
        return mean(pfd), var(pfd)
    end
end

"""
    shared_fixed_mutations(population[, idx])
"""
function shared_fixed_mutations(population, idx=nothing)
    clonalmuts = clonal_mutation_ids(population, idx)
    return shared_fixed_mutations(clonalmuts)
end

function shared_fixed_mutations(clonalmuts::Vector{Vector{Int64}})
    clonalmuts_vec = reduce(union, clonalmuts)
    nclonalmuts = map(
        x -> number_modules_with_mutation(clonalmuts, x), clonalmuts_vec
    )
    return countmap(nclonalmuts)
end

function number_modules_with_mutation(clonalmuts_by_module, mutationid)
    n = 0
    for muts in clonalmuts_by_module
        if mutationid in muts
            n += 1
        end
    end
    return n
end

"""
    newmoduletimes(population)

Return a vector of times at which new modules arose in the population
"""
newmoduletimes(population) = sort([moduletracker.tvec[1] for moduletracker in population])

"""
    cellpopulationsize(population, tstep)

Return number of cells in the population at times in 1:tstep:tend
"""
function cellpopulationsize(population, tstep)
    tend = maximum(moduletracker.tvec[end] for moduletracker in population)
    popvec = Int64[]
    for time in 0:tstep:tend
        pop = 0
        for moduletracker in population
            N0 = 0
            for (N, t) in zip(moduletracker.Nvec, moduletracker.tvec)
                if t > time
                    pop += N0
                    break
                elseif t == moduletracker.tvec[end]
                    pop += N
                    break
                else
                    N0 = N
                end
            end
        end
        push!(popvec, pop)
    end
    return collect(0:tstep:tend), popvec
end


text:
using SlackSDK
using Test
using Random

@testset "construct webclient" begin
    @test_throws ErrorException WebClient()
    t = "xxob-12345-abcde"
    client = WebClient(t)
    @test client !== nothing
    b = IOBuffer()
    show(b, client)
    @test String(take!(b)) == "Slack Web API client, with base url 'https://www.slack.com/api/' and token 'xxob-123....de'"
    ENV["SLACK_API_TOKEN"] = t
    client = WebClient(t)
    @test client !== nothing
    b = IOBuffer()
    show(b, client)
    @test String(take!(b)) == "Slack Web API client, with base url 'https://www.slack.com/api/' and token 'xxob-123....de'"
end

@testset "method availability" begin
    client = WebClient("xxob-12345-abcde")  # dummy token
    @test getproperty(client, :chat_delete) !== nothing
    @test getproperty(client, :chat_postMessage) !== nothing
    @test getproperty(client, :admin_apps_approved_list) !== nothing
    @test_throws UndefVarError getproperty(client, :chat_doesnotexist)
    @test_throws UndefVarError getproperty(client, :doesnotexist_delete)
end

## TODO : Actually exercise the api. Needs workspace and token