text:
function readRetirement( mpSim::MPsim, xf::XF )

    if !XLSX.hassheet( xf, "Default retirement" )
        @warn "Configuration file is missing a 'Default retirement' sheet, cannot define a default, fallback retirement scheme for the organisation."
        return
    end  # if !XLSX.hassheet( xf, "Default retirement" )

    ws = xf[ "Default retirement" ]
    retirement = Retirement()

    # Set schedule of retirement.
    freq, offset = ws[ "B3" ], ws[ "B4" ]
    setRetirementSchedule!( retirement, freq isa Real ? freq : 0,
        offset isa Real ? offset : 0 )

    # Set retirement conditions.
    tenure, age = ws[ "B5" ], ws[ "B6" ]
    setRetirementCareerLength!( retirement,
        tenure isa Real ? tenure * 12.0 : 0.0 )
    setRetirementAge!( retirement, age isa Real ? age * 12.0 : 0.0 )
    setRetirementIsEither!( retirement, uppercase( ws[ "B7" ] ) == "EITHER" )

    setSimulationRetirement!( mpSim, retirement )

end  # readRetirement( mpSim, xf )

function readRetirement( mpSim::MPsim, configDB::SQLite.DB,
    configName::String )

    retirePars = DataFrame( DBInterface.execute( configDB,
        string( "SELECT * FROM `", configName,
        "` WHERE parType IS 'Retirement'" ) ) )

    if !isempty( retirePars )
        retirement = Retirement()
        pars = split( retirePars[ 1, :parValue ], ";" )

        if ( tryparse( Float64, pars[ 1 ] ) isa Float64 ) &&
            ( tryparse( Float64, pars[ 2 ] ) isa Float64 )
            setRetirementSchedule!( retirement, parse( Float64, pars[ 1 ] ),
                parse( Float64, pars[ 2 ] ) )
        end  # if ( tryparse( Float64, pars[ 1 ] ) isa Float64 ) && ...

        if tryparse( Float64, pars[ 3 ] ) isa Float64
            setRetirementCareerLength!( retirement,
                parse( Float64, pars[ 3 ] ) )
        end  # if tryparse( Float64, pars[ 3 ] ) isa Float64

        if tryparse( Float64, pars[ 4 ] ) isa Float64
            setRetirementAge!( retirement, parse( Float64, pars[ 4 ] ) )
        end  # if tryparse( Float64, pars[ 4 ] ) isa Float64

        if tryparse( Bool, pars[ 5 ] ) isa Bool
            setRetirementIsEither!( retirement, parse( Bool, pars[ 5 ] ) )
        end  # if tryparse( Bool, pars[ 5 ] ) isa Bool

        setSimulationRetirement!( mpSim, retirement )
    end  # if !isempty( retirePars )

end  # readRetirement( mpSim, configDB, configName )


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

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

ensuretuple(s::Tuple) = s
ensuretuple(s) = (s,)

tupletopair(s::Tuple) = Pair(s...)

tuplemost(t::NTuple{N,Any}) where {N} = ntuple(i -> t[i], Val(N-1))

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

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

tupleproduct(p1, p2) = tupleproduct(ensuretuple(p1), ensuretuple(p2))
tupleproduct(p1::NTuple{M,Any}, p2::NTuple{N,Any}) where {M,N} =
    ntuple(i -> (p1[1+fld(i-1, N)], p2[1+mod(i-1, N)]), Val(M * N))

mergetuples(ts...) = keys(merge(tonamedtuple.(ts)...))
tonamedtuple(ts::Val{T}) where {T} = NamedTuple{T}(filltuple(0,T))

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

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

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

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

# Pad element type to a "larger" type
@inline padtotype(s::SMatrix{E,L}, ::Type{S}) where {E,L,E2,L2,S<:SMatrix{E2,L2}} =
    S(SMatrix{E2,E}(I) * s * SMatrix{L,L2}(I))
@inline padtotype(x::Number, ::Type{S}) where {E,L,S<:SMatrix{E,L}} =
    S(x * (SMatrix{E,1}(I) * SMatrix{1,L}(I)))
@inline padtotype(x::Number, ::Type{T}) where {T<:Number} = T(x)
@inline padtotype(u::UniformScaling, ::Type{T}) where {T<:Number} = T(u.λ)
@inline padtotype(u::UniformScaling, ::Type{S}) where {S<:SMatrix} = S(u)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

######################################################################
# convert a matrix/number block to a matrix/inlinematrix string
######################################################################

_isreal(x) = all(o -> imag(o) ≈ 0, x)
_isimag(x) = all(o -> real(o) ≈ 0, x)

matrixstring(row, x) = string("Onsite[$row] : ", _matrixstring(x))
matrixstring(row, col, x) = string("Hopping[$row, $col] : ", _matrixstring(x))

matrixstring_inline(row, x) = string("Onsite[$row] : ", _matrixstring_inline(x))
matrixstring_inline(row, col, x) = string("Hopping[$row, $col] : ", _matrixstring_inline(x))

_matrixstring(x::Number) = numberstring(x)
_matrixstring_inline(x::Number) = numberstring(x)
function _matrixstring(s::SMatrix)
    ss = repr("text/plain", s)
    pos = findfirst(isequal('\n'), ss)
    return pos === nothing ? ss : ss[pos:end]
end

function _matrixstring_inline(s::SMatrix{N}) where {N}
    stxt = numberstring.(transpose(s))
    stxt´ = vcat(stxt, SMatrix{1,N}(ntuple(_->";", Val(N))))
    return string("[", stxt´[1:end-1]..., "]")
end

numberstring(x) = _isreal(x) ? string(" ", real(x)) : _isimag(x) ? string(" ", imag(x), "im") : string(" ", x)

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

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

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

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

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

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

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

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

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

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

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

# ######################################################################
# # SparseMatrixIJV
# ######################################################################

# struct SparseMatrixIJV{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
#     I::Vector{Ti}
#     J::Vector{Ti}
#     V::Vector{Tv}
#     m::Ti
#     n::Ti
#     klasttouch::Vector{Ti}
#     csrrowptr::Vector{Ti}
#     csrcolval::Vector{Ti}
#     csrnzval::Vector{Tv}
#     csccolptr::Vector{Ti}
#     cscrowval::Vector{Ti}
#     cscnzval::Vector{Tv}
# end

# SparseMatrixIJV{Tv}(m::Ti, n::Ti) where {Tv,Ti} = SparseMatrixIJV{Tv,Ti}(m,n)

# function SparseMatrixIJV{Tv,Ti}(m::Integer, n::Integer; hintnnz = 0) where {Tv,Ti}
#     I = Ti[]
#     J = Ti[]
#     V = Tv[]
#     klasttouch = Vector{Ti}(undef, n)
#     csrrowptr = Vector{Ti}(undef, m + 1)
#     csrcolval = Vector{Ti}()
#     csrnzval = Vector{Tv}()
#     csccolptr = Vector{Ti}(undef, n + 1)
#     cscrowval = Vector{Ti}()
#     cscnzval = Vector{Tv}()

#     if hintnnz > 0
#         sizehint!(I, hintnnz)
#         sizehint!(J, hintnnz)
#         sizehint!(V, hintnnz)
#         sizehint!(csrcolval, hintnnz)
#         sizehint!(csrnzval, hintnnz)
#         sizehint!(cscrowval, hintnnz)
#         sizehint!(cscnzval, hintnnz)
#     end

#     return SparseMatrixIJV{Tv,Ti}(I, J, V, m, n, klasttouch, csrrowptr, csrcolval, csrnzval,
#                                                              csccolptr, cscrowval, cscnzval)
# end

# Base.summary(::SparseMatrixIJV{Tv,Ti}) where {Tv,Ti} =
#     "SparseMatrixIJV{$Tv,$Ti} : Sparse matrix builder using the IJV format"

# function Base.show(io::IO, ::MIME"text/plain", s::SparseMatrixIJV)
#     i = get(io, :indent, "")
#     print(io, i, summary(s), "\n", "$i  Nonzero elements : $(length(s.I))")
# end

# function Base.push!(s::SparseMatrixIJV, (i, j, v))
#     push!(s.I, i)
#     push!(s.J, j)
#     push!(s.V, v)
#     return s
# end

# function SparseArrays.sparse(s::SparseMatrixIJV)
#     numnz = length(s.I)
#     resize!(s.csrcolval, numnz)
#     resize!(s.csrnzval,  numnz)
#     resize!(s.cscrowval, numnz)
#     resize!(s.cscnzval,  numnz)
#     return SparseArrays.sparse!(s.I, s.J, s.V, s.m, s.n, +, s.klasttouch,
#         s.csrrowptr, s.csrcolval, s.csrnzval, s.csccolptr, s.cscrowval, s.cscnzval)
# end

# Base.size(s::SparseMatrixIJV) = (s.m, s.n)

############################################################################################
######## fast sparse copy #  Revise after #33589 is merged #################################
############################################################################################

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

# Only needed for dense <- sparse (#33589), copy!(sparse, sparse) is fine in v1.5+
function _fast_sparse_copy!(dst::AbstractMatrix{T}, src::SparseMatrixCSC) where {T}
    @boundscheck checkbounds(dst, axes(src)...)
    fill!(dst, zero(eltype(src)))
    for col in 1:size(src, 1)
        for p in nzrange(src, col)
            @inbounds dst[rowvals(src)[p], col] = nonzeros(src)[p]
        end
    end
    return dst
end

function _fast_sparse_muladd!(dst::AbstractMatrix{T}, src::SparseMatrixCSC, α = I) where {T}
    @boundscheck checkbounds(dst, axes(src)...)
    for col in 1:size(src, 1)
        for p in nzrange(src, col)
            @inbounds dst[rowvals(src)[p], col] += α * nonzeros(src)[p]
        end
    end
    return dst
end

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

ext: jl


text:
function separate_columns(data::Vector{Vector{T}}) where T <: Number
    if length(data) > 0
        columns_number = length(data[1])
        columns = Vector{Vector{Float64}}(undef, columns_number)
        for i = 1:columns_number
            columns[i] = getindex.(data, i)
        end
        return columns
    else
        error("Null data vector.")
    end
end

function separate_columns(data::Vector{Vector{T}}) where T <: Any
    if length(data) > 0
        columns_number = length(data[1])
        type = typeof(data[1])
        columns = Vector{type}(undef, columns_number)
        for i = 1:columns_number
            columns[i] = getindex.(data, i)
        end
        return columns
    else
        error("Null data vector.")
    end
end

function correct_upper_bounds!(vector, bounds)
    if length(vector) == length(bounds)
        for i = 1:length(vector)
            if vector[i] > bounds[i]
                vector[i] = bounds[i]
            end
        end
    else
        error("Bounds were not set to all elements.")
    end
end

function correct_upper_bounds(vector, bounds)
    if length(vector) == length(bounds)
        new_vector = typeof(vector)(undef, length(vector))
        for i = 1:length(vector)
            if vector[i] > bounds[i]
                new_vector[i] = bounds[i]
            else
                new_vector[i] = vector[i]
            end
        end
        return new_vector
    else
        error("Bounds were not set to all elements.")
    end
end

function get_bin_from_indices(vector, bounds)
    if length(vector) == length(bounds)
        index = 0
        len = length(vector)
        for i = 1:(length(vector)-1)
            index += (vector[i]-1)*prod(bounds[(i+1):end])
        end
        if len >= 1 index += vector[end] end
        return index
    else
        error("Vector and bound vector have different lengths.")
    end
end

get_weights_length(h::Histogram) = length(h.weights)

function flat_vector(vector; nbins = 5, lower_lim=-1, upper_lim=1)
    step = (upper_lim-lower_lim)/nbins
    columns = separate_columns(vector)
    hists = fit.(Histogram, columns, [lower_lim:step:upper_lim])
    bounds = (get_weights_length.(hists)...,)
    indices = [getindex.(get_bin.(hists, v),1) for v in vector]
    indices = [correct_upper_bounds(ind, bounds) for ind in indices]
    get_bin_from_indices.(indices, [bounds])
end

#=
function get_layer_bins_vector(vector; nbins = 5, lower_lim=-1, upper_lim=1)
    columns = separate_columns(vector)
    hists = create_bin_histogram.(columns, [nbins], [lower_lim], [upper_lim])
    bounds = (get_weights_length.(hists)...,)
    indices = [getindex.(get_bin.(hists, v),1) for v in vector]
    indices = [correct_upper_bounds(ind, bounds) for ind in indices]
    combination_values = get_bin_from_indices.(indices, [bounds])
    max_value = foldr(*, bounds)
    create_bin_vector(combination_values, nbins, 1, max_value)
end
=#

function get_layer_bins_vector(vector; neuron_nbins = 5, nbins = 5, lower_lim=-1, upper_lim=1)
    columns = separate_columns(vector)
    hists = create_bin_histogram.(columns, [neuron_nbins], [lower_lim], [upper_lim])
    bounds = (get_weights_length.(hists)...,)
    indices = [getindex.(get_bin.(hists, v),1) for v in vector]
    indices = [correct_upper_bounds(ind, bounds) for ind in indices]
    combination_values = get_bin_from_indices.(indices, [bounds])
    max_value = foldr(*, bounds)
    create_bin_vector(combination_values, nbins, 1, max_value)
end

#=
function fit_into_histogram(vector, nbins::Int, lower_lim, upper_lim; ϵ = 0.01)
    hist_step = (upper_lim-lower_lim)/nbins
    ϵ = 0.01
    edges = collect(lower_lim:hist_step:upper_lim)
    edges[1] -= ϵ; edges[end] += ϵ
    hist = fit(Histogram, vector, edges)
    [getindex.(get_bin(hist, v),1) for v in vector]
end
=#

function create_bin_histogram(vector::Vector{T}, nbins::Int, lower_lim, upper_lim; ϵ = 0.01) where T <: Number
    hist_step = abs(upper_lim-lower_lim)/nbins
    edges = collect(lower_lim:hist_step:upper_lim)
    edges[1] -= ϵ; edges[end] += ϵ
    hist = fit(Histogram, vector, edges)
end

function create_bin_vector(vector::Vector{T}, nbins::Int, lower_lim, upper_lim; ϵ = 0.01) where T <: Number
    hist = create_bin_histogram(vector, nbins, lower_lim, upper_lim, ϵ = ϵ)
    [getindex.(get_bin(hist, v),1) for v in vector]
end

function get_input_bin_vector(vector, nfeatures; nbins = 5, return_bins::Int=0)
    hists = Vector{Vector{Int}}(undef, nfeatures)
    columns = separate_columns(vector)
    for feature in 1:nfeatures
        feature_values = columns[feature]
        min = minimum(feature_values)
        max = maximum(feature_values)
        hists[feature] = create_bin_vector(feature_values, nbins, min, max)
    end
    bounds = repeat([nbins], nfeatures)
    indices = separate_columns(hists)
    combination_values = get_bin_from_indices.(indices, [bounds])
    if return_bins > 0
        max_value = foldr(*, bounds)
        return create_bin_vector(combination_values, return_bins, 1, max_value)
    else
        return floor.(Int, combination_values)
    end
end

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
import Base: convert
import Base: ==

# abstract type InferenceType end
# abstract type PackedInferenceType end
#
# abstract type FunctorInferenceType <: Function end
#
# abstract type InferenceVariable end
# abstract type ConvolutionObject <: Function end

# abstract type FunctorSingleton <: FunctorInferenceType end
# abstract type FunctorPairwise <: FunctorInferenceType end
# abstract type FunctorPairwiseMinimize <: FunctorInferenceType end

# TODO been replaced by Functor types, but may be reused for non-numerical cases
abstract type Pairwise <: InferenceType end
abstract type Singleton <: InferenceType end

# TODO deprecate with standard null hypothesis only
abstract type FunctorSingletonNH <: FunctorSingleton end
abstract type FunctorPairwiseNH <: FunctorPairwise end
# abstract type FunctorPairwiseNHMinimize <: FunctorPairwiseMinimize end # TODO

const FGG = Graphs.GenericIncidenceList{TreeClique,Graphs.Edge{TreeClique},Array{TreeClique,1},Array{Array{Graphs.Edge{TreeClique},1},1}}
const FGGdict = Graphs.GenericIncidenceList{TreeClique,Graphs.Edge{TreeClique},Dict{Int,TreeClique},Dict{Int,Array{Graphs.Edge{TreeClique},1}}}

"""
$(TYPEDEF)

Solver parameters for the DistributedFactoGraph.

Dev Notes
- TODO remove NothingUnion
"""
mutable struct SolverParams <: DFG.AbstractParams
  dimID::Int
  # TODO remove NothingUnion
  registeredModuleFunctions::NothingUnion{Dict{Symbol, Function}} # remove from
  reference::NothingUnion{Dict{Symbol, Tuple{Symbol, Vector{Float64}}}}
  stateless::Bool
  qfl::Int # Quasi fixed length
  isfixedlag::Bool # true when adhering to qfl window size for solves
  limitfixeddown::Bool # if true, then fixed lag will also not update marginalized during down (default false)
  # new functions
  incremental::Bool
  upsolve::Bool
  downsolve::Bool
  drawtree::Bool
  showtree::Bool
  dbg::Bool
  async::Bool
  limititers::Int
  N::Int
  multiproc::Bool
  logpath::String
  devParams::Dict{Symbol,String}
  SolverParams(;dimID::Int=0,
                registeredModuleFunctions=nothing,
                reference=nothing,
                stateless::Bool=false,
                qfl::Int=99999999999,
                isfixedlag::Bool=false,
                limitfixeddown::Bool=false,
                incremental::Bool=true,
                upsolve::Bool=true,
                downsolve::Bool=true,
                drawtree::Bool=false,
                showtree::Bool=false,
                dbg::Bool=false,
                async::Bool=false,
                limititers::Int=500,
                N::Int=100,
                multiproc::Bool=true,
                logpath::String="/tmp/caesar/$(now())",
                devParams::Dict{Symbol,String}=Dict{Symbol,String}()) = new(dimID,
                                                                            registeredModuleFunctions,
                                                                            reference,
                                                                            stateless,
                                                                            qfl,
                                                                            isfixedlag,
                                                                            limitfixeddown,
                                                                            incremental,
                                                                            upsolve,
                                                                            downsolve,
                                                                            drawtree,
                                                                            showtree,
                                                                            dbg,
                                                                            async,
                                                                            limititers,
                                                                            N,
                                                                            multiproc,
                                                                            logpath,
                                                                            devParams)
  #
end

"""
$(TYPEDEF)

NOTE: Deprecated by DistributedFactorGraphs.
"""
mutable struct FactorGraph
  g::FGGdict
  bn
  IDs::Dict{Symbol,Int}
  fIDs::Dict{Symbol,Int}
  id::Int
  nodeIDs::Array{Int,1} # TODO -- ordering seems improved to use adj permutation -- pending merge JuliaArchive/Graphs.jl/#225
  factorIDs::Array{Int,1}
  bnverts::Dict{Int,TreeClique} # TODO -- not sure if this is still used, remove
  bnid::Int # TODO -- not sure if this is still used
  dimID::Int
  cg
  cgIDs::Dict{Int,Int} # cgIDs[exvid] = neoid
  sessionname::String
  robotname::String
  username::String
  registeredModuleFunctions::NothingUnion{Dict{Symbol, Function}}
  reference::NothingUnion{Dict{Symbol, Tuple{Symbol, Vector{Float64}}}}
  stateless::Bool
  fifo::Vector{Symbol}
  qfl::Int # Quasi fixed length
  isfixedlag::Bool # true when adhering to qfl window size for solves
  FactorGraph(;reference::NothingUnion{Dict{Symbol, Tuple{Symbol, Vector{Float64}}}}=nothing, is_directed::Bool=true ) = new(Graphs.incdict(TreeClique,is_directed=false),
                      Graphs.incdict(TreeClique,is_directed=is_directed),
                      #  Dict{Int,TreeClique}(),
                      #  Dict{Int,TreeClique}(),
                      Dict{Symbol,Int}(),
                      Dict{Symbol,Int}(),
                      0,
                      [],
                      [],
                      Dict{Int,TreeClique}(),
                      0,
                      0,
                      nothing,
                      Dict{Int,Int}(),
                      "",
                      "",
                      "",
                      Dict{Symbol, Function}(:IncrementalInference=>IncrementalInference.getSample), # TODO likely to be removed
                      reference,
                      false,
                      Symbol[],
                      0,
                      false  )
end

"""
    $SIGNATURES

Initialize an empty in-memory DistributedFactorGraph `::DistributedFactorGraph` object.
"""
function initfg(dfg::T=InMemDFGType(params=SolverParams());
                                               sessionname="NA",
                                               robotname="",
                                               username="",
                                               cloudgraph=nothing)::T where T <: AbstractDFG
  #
  return dfg
end

#init an empty fg with a provided type and SolverParams
function initfg(::Type{T}; params=SolverParams(),
                           sessionname="NA",
                           robotname="",
                           username="",
                           cloudgraph=nothing)::AbstractDFG where T <: AbstractDFG
  return T(params=params)
end

function initfg(::Type{T}, params::SolverParams;
                           sessionname="NA",
                           robotname="",
                           username="",
                           cloudgraph=nothing)::AbstractDFG where T <: AbstractDFG
  return T{SolverParams}(params=params)
end

"""
$(TYPEDEF)

TODO remove Union types -- issue #383
"""
mutable struct FactorMetadata
  factoruserdata
  variableuserdata::Union{Vector, Tuple}
  variablesmalldata::Union{Vector, Tuple}
  solvefor::Union{Symbol, Nothing}
  variablelist::Union{Nothing, Vector{Symbol}}
  dbg::Bool
  FactorMetadata() = new() # [], []
  FactorMetadata(x1, x2::Union{Vector,Tuple},x3) = new(x1, x2, x3, nothing, nothing, false)
  FactorMetadata(x1, x2::Union{Vector,Tuple},x3,x4::Symbol) = new(x1, x2, x3, x4, nothing, false)
  FactorMetadata(x1, x2::Union{Vector,Tuple},x3,x4::Symbol,x5::Vector{Symbol};dbg::Bool=false) = new(x1, x2, x3, x4, x5, dbg)
end

"""
$(TYPEDEF)
"""
struct SingleThreaded
end
"""
$(TYPEDEF)
"""
struct MultiThreaded
end

"""
$(TYPEDEF)
"""
mutable struct ConvPerThread
  thrid_::Int
  particleidx::Int # the actual particle being solved at this moment
  factormetadata::FactorMetadata # additional data passed to user function -- optionally used by user function
  activehypo::Union{UnitRange{Int},Vector{Int}} # subsection indices to select which params should be used for this hypothesis evaluation
  p::Vector{Int} # a permutation vector for low-dimension solves (FunctorPairwise only)
  perturb::Vector{Float64} # slight numerical perturbation for degenerate solver cases such as division by zero
  X::Array{Float64,2}
  Y::Vector{Float64}
  res::Vector{Float64}
  ConvPerThread() = new()
end

function ConvPerThread(X::Array{Float64,2},
                       zDim::Int;
                       factormetadata::FactorMetadata=FactorMetadata(),
                       particleidx::Int=1,
                       activehypo= 1:length(params),
                       p=collect(1:size(X,1)),
                       perturb=zeros(zDim),
                       Y=zeros(size(X,1)),
                       res=zeros(zDim)  )
  #
  cpt = ConvPerThread()
  cpt.thrid_ = 0
  cpt.X = X
  cpt.factormetadata = factormetadata
  cpt.particleidx = particleidx
  cpt.activehypo = activehypo
  cpt.p = p
  cpt.perturb = perturb
  cpt.Y = Y
  cpt.res = res
  return cpt
end

"""
$(TYPEDEF)
"""
mutable struct CommonConvWrapper{T} <: ConvolutionObject where {T<:FunctorInferenceType}
  ### Values consistent across all threads during approx convolution
  usrfnc!::T # user factor / function
  # general setup
  xDim::Int
  zDim::Int
  # special case settings
  specialzDim::Bool # is there a special zDim requirement -- defined by user
  partial::Bool # is this a partial constraint -- defined by user
  # multi hypothesis settings
  hypotheses::Union{Nothing, Distributions.Categorical} # categorical to select which hypothesis is being considered during convolugtion operation
  certainhypo::Union{Nothing, Vector{Int}}
  # values specific to one complete convolution operation
  params::Vector{Array{Float64,2}} # parameters passed to each hypothesis evaluation event on user function
  varidx::Int # which index is being solved for in params?
  measurement::Tuple # user defined measurement values for each approxConv operation
  threadmodel::Union{Type{SingleThreaded}, Type{MultiThreaded}}

  ### particular convolution computation values per particle idx (varies by thread)
  cpt::Vector{ConvPerThread}
  # varidx::Int # which index is being solved for in params?
  # factormetadata::FactorMetadata # additional data passed to user function -- optionally used by user function
  # activehypo::Union{UnitRange{Int},Vector{Int}} # subsection indices to select which params should be used for this hypothesis evaluation
  # particleidx::Int # the actual particle being solved at this moment
  # p::Vector{Int} # a permutation vector for low-dimension solves (FunctorPairwise only)
  # perturb::Vector{Float64} # slight numerical perturbation for degenerate solver cases such as division by zero
  # X::Array{Float64,2}
  # Y::Vector{Float64}
  # res::Vector{Float64}

  CommonConvWrapper{T}() where {T<:FunctorInferenceType} = new{T}()
end

function CommonConvWrapper(fnc::T,
                           X::Array{Float64,2},
                           zDim::Int,
                           params::Vector{Array{Float64,2}};
                           factormetadata::FactorMetadata=FactorMetadata(),
                           specialzDim::Bool=false,
                           partial::Bool=false,
                           hypotheses=nothing,
                           certainhypo=nothing,
                           activehypo= 1:length(params),
                           varidx::Int=1,
                           measurement::Tuple=(zeros(0,1),),
                           particleidx::Int=1,
                           p=collect(1:size(X,1)),
                           perturb=zeros(zDim),
                           Y=zeros(size(X,1)),
                           xDim=size(X,1),
                           res=zeros(zDim),
                           threadmodel=MultiThreaded  ) where {T<:FunctorInferenceType}
  #
  ccw = CommonConvWrapper{T}()

  ccw.usrfnc! = fnc
  ccw.xDim = xDim
  ccw.zDim = zDim
  ccw.specialzDim = specialzDim
  ccw.partial = partial
  ccw.hypotheses = hypotheses
  ccw.certainhypo=certainhypo
  ccw.params = params
  ccw.varidx = varidx
  ccw.threadmodel = threadmodel
  ccw.measurement = measurement

  # thread specific elements
  ccw.cpt = Vector{ConvPerThread}(undef, Threads.nthreads())
  for i in 1:Threads.nthreads()
    ccw.cpt[i] = ConvPerThread(X, zDim,
                    factormetadata=factormetadata,
                    particleidx=particleidx,
                    activehypo=activehypo,
                    p=p,
                    perturb=perturb,
                    Y=Y,
                    res=res )
  end

  return ccw
end

## moved to DFG

# # where {T <: Union{InferenceType, FunctorInferenceType}}
# const FunctionNodeData{T} = GenericFunctionNodeData{T, Symbol}
# FunctionNodeData(x1, x2, x3, x4, x5::Symbol, x6::T, x7::String="", x8::Vector{Int}=Int[]) where {T <: Union{FunctorInferenceType, ConvolutionObject}}= GenericFunctionNodeData{T, Symbol}(x1, x2, x3, x4, x5, x6, x7, x8)

# # where {T <: PackedInferenceType}
# const PackedFunctionNodeData{T} = GenericFunctionNodeData{T, <: AbstractString}
# PackedFunctionNodeData(x1, x2, x3, x4, x5::S, x6::T, x7::String="", x8::Vector{Int}=Int[]) where {T <: PackedInferenceType, S <: AbstractString} = GenericFunctionNodeData(x1, x2, x3, x4, x5, x6, x7, x8)

# excessive function, needs refactoring
function updateFullVertData!(fgl::AbstractDFG,
                             nv::DFGNode;
                             updateMAPest::Bool=false )
  #
  @warn "Deprecated updateFullVertData!, need alternative"

  sym = Symbol(nv.label)
  isvar = isVariable(fgl, sym)

  lvert = isvar ? DFG.getVariable(fgl, sym) : DFG.getFactor(fgl, sym)
  lvd = solverData(lvert)
  nvd = solverData(nv)

  if isvar
    lvd.val[:,:] = nvd.val[:,:]
    lvd.bw[:] = nvd.bw[:]
    lvd.initialized = nvd.initialized
    lvd.inferdim = nvd.inferdim
  else
    # assuming nothing to be done
  end

  nothing
end


text:
"""
    test_variable_add_variable(model::MOI.ModelLike, config::Config)

Test adding a single variable.
"""
function test_variable_add_variable(model::MOI.ModelLike, ::Config)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    v = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1
    _test_attribute_value_type(model, MOI.NumberOfVariables())
    return
end

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

Test adding multiple variables.
"""
function test_variable_add_variables(model::MOI.ModelLike, ::Config)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    v = MOI.add_variables(model, 2)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2
    return
end

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

Tess adding, and then deleting, a single variable.
"""
function test_variable_delete(model::MOI.ModelLike, ::Config)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    v = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1
    MOI.delete(model, v)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    return
end

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

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

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

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

"""
    test_variable_delete_Nonnegatives_row(model::MOI.ModelLike, ::Config)

Test adding, and then deleting one by one, nonnegative variables.
"""
function test_variable_delete_Nonnegatives_row(model::MOI.ModelLike, ::Config)
    @requires MOI.supports_add_constrained_variables(model, MOI.Nonnegatives)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    v, cv = MOI.add_constrained_variables(model, MOI.Nonnegatives(2))
    @test MOI.get(model, MOI.NumberOfVariables()) == 2
    MOI.delete(model, v[1])
    @test !MOI.is_valid(model, v[1])
    @test_throws MOI.InvalidIndex(v[1]) MOI.delete(model, v[1])
    @test MOI.is_valid(model, cv)
    @test MOI.is_valid(model, v[2])
    MOI.delete(model, v[2])
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    @test !MOI.is_valid(model, v[1])
    @test_throws MOI.InvalidIndex(v[1]) MOI.delete(model, v[1])
    @test !MOI.is_valid(model, v[2])
    @test_throws MOI.InvalidIndex(v[2]) MOI.delete(model, v[2])
    @test !MOI.is_valid(model, cv)
    return
end

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

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

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

Test getting variables by name.
"""
function test_variable_get_VariableIndex(model::MOI.ModelLike, ::Config)
    @requires MOI.supports(model, MOI.VariableName(), MOI.VariableIndex)
    variable = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), variable, "x")
    x = MOI.get(model, MOI.VariableIndex, "x")
    @test MOI.is_valid(model, x)
    @test MOI.get(model, MOI.VariableIndex, "y") === nothing
    return
end

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

Test getting and setting variable names.
"""
function test_variable_VariableName(model::MOI.ModelLike, ::Config)
    @requires MOI.supports(model, MOI.VariableName(), MOI.VariableIndex)
    v = MOI.add_variable(model)
    @test MOI.get(model, MOI.VariableName(), v) == ""
    MOI.set(model, MOI.VariableName(), v, "x")
    @test MOI.get(model, MOI.VariableName(), v) == "x"
    MOI.set(model, MOI.VariableName(), v, "y")
    @test MOI.get(model, MOI.VariableName(), v) == "y"
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    @test MOI.get(model, MOI.VariableName(), x) == "x"
    return
end

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

Test setting the upper bound of a variable, confirm that it solves correctly.
"""
function test_variable_solve_with_upperbound(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    c1 = MOI.add_constraint(model, x, MOI.LessThan(T(1)))
    @test x.value == c1.value
    c2 = MOI.add_constraint(model, x, MOI.GreaterThan(T(0)))
    @test x.value == c2.value
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(2), x)], T(0))
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    _test_model_solution(
        model,
        config;
        objective_value = T(2),
        variable_primal = [(x, T(1))],
        constraint_primal = [(c1, T(1)), (c2, T(1))],
        constraint_dual = [(c1, T(-2)), (c2, T(0))],
    )
    return
end

function setup_test(
    ::typeof(test_variable_solve_with_upperbound),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1]),
            MOI.FEASIBLE_POINT,
            (MOI.VariableIndex, MOI.LessThan{T}) => T[-2],
            (MOI.VariableIndex, MOI.GreaterThan{T}) => T[0],
        ),
    )
    model.eval_variable_constraint_dual = false
    return () -> model.eval_variable_constraint_dual = true
end

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

Test setting the lower bound of a variable, confirm that it solves correctly.
"""
function test_variable_solve_with_lowerbound(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    c1 = MOI.add_constraint(model, x, MOI.GreaterThan(T(1)))
    c2 = MOI.add_constraint(model, x, MOI.LessThan(T(2)))
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(2), x)], T(0))
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    _test_model_solution(
        model,
        config;
        objective_value = T(2),
        variable_primal = [(x, T(1))],
        constraint_primal = [(c1, T(1)), (c2, T(1))],
        constraint_dual = [(c1, T(2)), (c2, T(0))],
    )
    return
end

function setup_test(
    ::typeof(test_variable_solve_with_lowerbound),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, T[1]),
            MOI.FEASIBLE_POINT,
            (MOI.VariableIndex, MOI.GreaterThan{T}) => T[2],
            (MOI.VariableIndex, MOI.LessThan{T}) => T[0],
        ),
    )
    model.eval_variable_constraint_dual = false
    return () -> model.eval_variable_constraint_dual = true
end

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

Test an integer variable with fractional lower bound.
"""
function test_variable_solve_Integer_with_lower_bound(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.GreaterThan(T(3 // 2)))
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(2), x)], T(0))
    MOI.add_constraint(model, x, MOI.Integer())
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    _test_model_solution(
        model,
        config;
        objective_value = T(4),
        variable_primal = [(x, T(2))],
    )
    return
end

function setup_test(
    ::typeof(test_variable_solve_Integer_with_lower_bound),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, T[2])),
    )
    return
end

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

Test an integer variable with fractional upper bound.
"""
function test_variable_solve_Integer_with_upper_bound(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.LessThan(T(3 // 2)))
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(-2), x)], T(0))
    MOI.add_constraint(model, x, MOI.Integer())
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    _test_model_solution(
        model,
        config;
        objective_value = T(-2),
        variable_primal = [(x, T(1))],
    )
    return
end

function setup_test(
    ::typeof(test_variable_solve_Integer_with_upper_bound),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, T[1])),
    )
    return
end

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

Test a binary variable `<= 2`.
"""
function test_variable_solve_ZeroOne_with_upper_bound(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.LessThan(T(2)))
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(-2), x)], T(0))
    MOI.add_constraint(model, x, MOI.ZeroOne())
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    _test_model_solution(
        model,
        config;
        objective_value = T(-2),
        variable_primal = [(x, T(1))],
    )
    return
end

function setup_test(
    ::typeof(test_variable_solve_ZeroOne_with_upper_bound),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, T[1])),
    )
    return
end

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

Test a binary variable `<= 0`.
"""
function test_variable_solve_ZeroOne_with_0_upper_bound(
    model::MOI.ModelLike,
    config::Config{T},
) where {T}
    x = MOI.add_variable(model)
    MOI.add_constraint(model, x, MOI.LessThan(T(0)))
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(T(1), x)], T(0))
    MOI.add_constraint(model, x, MOI.ZeroOne())
    MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    _test_model_solution(
        model,
        config;
        objective_value = T(0),
        variable_primal = [(x, T(0))],
    )
    return
end

function setup_test(
    ::typeof(test_variable_solve_ZeroOne_with_0_upper_bound),
    model::MOIU.MockOptimizer,
    ::Config{T},
) where {T}
    MOIU.set_mock_optimize!(
        model,
        (mock::MOIU.MockOptimizer) ->
            MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, T[0])),
    )
    return
end

"""
    test_variable_VariablePrimalStart(model::MOI.ModelLike, ::Config)

Test VariablePrimalStart
"""
function test_model_VariablePrimalStart(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex)
    x = MOI.add_variable(model)
    @test MOI.get(model, MOI.VariablePrimalStart(), x) === nothing
    MOI.set(model, MOI.VariablePrimalStart(), x, T(1))
    @test MOI.get(model, MOI.VariablePrimalStart(), x) == T(1)
    MOI.set(model, MOI.VariablePrimalStart(), x, nothing)
    @test MOI.get(model, MOI.VariablePrimalStart(), x) === nothing
    return
end

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

Test vector method of `add_constrained_variables`.
"""
function test_add_constrained_variables_vector(
    model::MOI.ModelLike,
    ::Config{T},
) where {T}
    @requires MOI.supports_add_constrained_variable(model, MOI.GreaterThan{T})
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    sets = [MOI.GreaterThan(zero(T)), MOI.GreaterThan(one(T))]
    v, cv = MOI.add_constrained_variables(model, sets)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2
    @test length(v) == 2
    @test length(cv) == 2
    @test MOI.get(model, MOI.ConstraintSet(), cv[1]) == sets[1]
    @test MOI.get(model, MOI.ConstraintSet(), cv[2]) == sets[2]
    return
end


text:
module Write

import ..Git, ..Cache, ..Read

function prefetch(pkg::String, sha1::String)
    isempty(Cache.prefetch(pkg, Read.url(pkg), sha1)) && return
    error("$pkg: couldn't find commit $(sha1[1:10])")
end

function fetch(pkg::String, sha1::String)
    refspec = "+refs/*:refs/remotes/cache/*"
    Git.run(`fetch -q $(Cache.path(pkg)) $refspec`, dir=pkg)
    Git.iscommit(sha1, dir=pkg) && return
    f = Git.iscommit(sha1, dir=Cache.path(pkg)) ? "fetch" : "prefetch"
    error("$pkg: $f failed to get commit $(sha1[1:10]), please file a bug")
end

function checkout(pkg::String, sha1::String)
    Git.set_remote_url(Read.url(pkg), dir=pkg)
    Git.run(`checkout -q $sha1`, dir=pkg)
end

function install(pkg::String, sha1::String)
    prefetch(pkg, sha1)
    if isdir(".trash/$pkg")
        run(`mv .trash/$pkg ./`)
    else
        Git.run(`clone -q $(Cache.path(pkg)) $pkg`)
    end
    fetch(pkg, sha1)
    checkout(pkg, sha1)
end

function update(pkg::String, sha1::String)
    prefetch(pkg, sha1)
    fetch(pkg, sha1)
    checkout(pkg, sha1)
end

function remove(pkg::String)
    isdir(".trash") || mkdir(".trash")
    ispath(".trash/$pkg") && run(`rm -rf .trash/$pkg`)
    run(`mv $pkg .trash/`)
end

end # module

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
"""
	transferoperator_exact(t::AbstractTriangulation)

Compute the transfer operator from a triangulation.
"""
function transferoperator_exact(t::AbstractTriangulation)
    n_simplices = size(t.simplex_inds, 1)
    dim = size(t.points, 2)

    #=
    Mismatch for the markovity of the Markov matrix is at most δ if ϵ tolerance
	as below
    =#
    δ::Float64 = 1/10^5
    ϵ::Float64 = δ/n_simplices

    #=
    Tolerance for similary of convex expansion coefficients of simplex vertices in
    simplexintersection function.
    =#
    convex_params_tol::Float64 = 1/10^12

    TO = zeros(Float64, n_simplices, n_simplices) # intersecting volumes

    for i in 1:n_simplices
        imvol = t.volumes_im[i]
        for j in 1:n_simplices
            vol = t.volumes[j]
            if vol * imvol > 0 && (vol/imvol) > ϵ
                # Intersecting volume between these two simplices
                TO[i, j] = simplexintersection(
                                transpose(t.points[t.simplex_inds[j, :], :]),
                                transpose(t.impoints[t.simplex_inds[i, :], :])
                            ) / imvol
            end
        end
    end
    return ExactSimplexTransferOperator(TO)
end

"""
    transferoperator_exact_p(t::AbstractTriangulation)

Estimate the transfer operator from a triangulation of a state space.
Will run in parallel if `nprocs() > 1`.
"""
function transferoperator_exact_p(t::AbstractTriangulation)

    n_simplices = size(t.simplex_inds, 1)
    dim = size(t.points, 2)

    #=
    Mismatch for the markovity of the Markov matrix is at most δ if ϵ
    tolerance as below.
    =#
    δ::Float64 = 1/10^5
    ϵ::Float64 = δ/n_simplices

    #=
    Tolerance for similary of convex expansion coefficients of simplex vertices
    in simplexintersection function.
    =#
    convex_params_tol::Float64 = 1/10^12

    TO = SharedArray{Float64}(n_simplices, n_simplices)

    @sync @distributed for i in 1:n_simplices
        imvol = t.volumes_im[i]
        for j in 1:n_simplices
            vol = t.volumes[j]
            if vol * imvol > 0 && (vol/imvol) > ϵ
                # Intersecting volume between these two simplices
                TO[i, j] = simplexintersection(
                                transpose(t.points[t.simplex_inds[j, :], :]),
                                transpose(t.impoints[t.simplex_inds[i, :], :])
                            ) / imvol
            end
        end
    end

    return ExactSimplexTransferOperator(Array(TO))
end

transferoperator_triang_exact = transferoperator_exact
export transferoperator_triang_exact


text:
module TEST_NODE
using NRES
using Test

throw("error: node.jl (og test_node.jl): This file is not maintained: Planned for destruction.. ")
# TODO Sjå førre linje
#
#
#
#
#
#
#
#
#
#
#
#
struct DUMMY_container <: NRES.SAT_container end

@testset "SAT init" begin
    demoSet = DUMMY_container()
    case = Node.SAT(member_of_set=demoSet)
    @test isa(case, Node.SAT)
    " Node.SAT initializes an instance of SAT "

    #@test isa(case._member_of_group, NRES.Node.SAT_container)

    @test !Node.is_active(case)
    @test  Node.not_active(case)
    " A node is initially inactive // default: not_active "

    Node.activate!(case)
    @test Node.is_active(case)
    " Node.activate! "

    Node.deactivate!(case)
    @test Node.not_active(case)
    " Node.deactivate! "
end

end #module TEST_NODE

max_stars_repo_name: leikanger/NRES

max_issues_repo_name: leikanger/NRES

max_forks_repo_name: leikanger/NRES


text:
using TexExtensions
using Compat

using Base.Test

@test istextmime(MIME("text/mathtex+latex"))


text:
double printInt() {
   return 1;
}
int main() {
   printInt(1.0);
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
"""siceのシミュレーション"""

using Plots
using StaticArrays
using ArraysOfArrays
using Parameters
using LinearAlgebra
using Random
using LaTeXStrings

include("./sice_kinematics.jl")
include("./sice_robot_dynamics.jl")
include("./rmp.jl")

using .SiceKinematics
using .SiceDynamics

const q_min = [(-3/4)π+(1/2)π, (-3/4)π, (-3/4)π, (-3/4)π]
const q_max = [(3/4)π+(1/2)π, (3/4)π, (3/4)π, (3/4)π]
const q_neutral = [(1/2)π, 0.0, 0.0, 0.0]

function split_vec_of_arrays(u)
    vec.(u) |>
    x ->
    VectorOfSimilarVectors(x).data |>
    transpose |>
    VectorOfSimilarVectors
end

function calc_min_dis_to_obs(xs::Vector{Vector{T}}, xos::Vector{Vector{T}}) where T

    d = Matrix{T}(undef, length(xs), length(xos))
    for (i, x) in enumerate(xs)
        for (j, xo) in enumerate(xos)
            d[i, j] = norm(x - xo)
        end
    end

    return minimum(d)
end

"""ジョイント制限を守れてるかチェック"""
function check_JointLimitation(q::Vector{T}) where T
    return (q .< q_min) .| (q_max .< q)
end

"""円物体が与える外力"""
function Fc_from_cicle(
    x::Vector{T}, dx::Vector{T},
    center::Vector{T}, r::T,
    K::T, D::T
    ) where T
    s = norm(x .- center)

    if s >= r
        return zero(x)
    else
        n = (x - center) / s  # 力の方向
        return @. K * (r-s) .* n .+ D * dx
    end
end

function multi_avoidance_rmp(x::Vector{T}, dx::Vector{T}, o::Vector{T}) where T
    R = 0.3
    alpha = 0.01#1e-5
    eta = 0.0
    epsilon = 0.2

    z_vec = (x - o)
    dz_vec = dx

    # ## マルチロボットの安全距離付きタスク写像 ###
    # z = norm(z_vec)/R - 1
    # J = 1/norm(z_vec) * (z_vec)' / R
    # dJ = (dz_vec' * (-1/norm(z_vec)^3 * (z_vec) * (z_vec)' + 1/norm(z_vec) * Matrix{T}(I, 2, 2))) / R

    ### 普通の障害物との距離関数のタスク写像 ###
    z = norm(z_vec)
    dz = (1/z .* dot((z_vec), (dz_vec)))[1]
    J = 1/z * (z_vec)'
    dJ = -z^(-2) .* (dz_vec' .- ((z_vec)' .* dz))

    dz = J * dx

    if z < 0
        w = 1.0e10
        grad_w = 0.0
    else
        w = 1 / z^4
        grad_w = -4 / z^5
    end

    u = epsilon + min(0.0, dz) * dz
    g = w * u

    grad_u = 2 * min(0.0, dz)
    grad_phi = alpha * w * grad_w
    xi = 1/2 * dz^2 * u * grad_w

    _M = g + 1/2 * dz * w * grad_u
    M = min(max(_M, -1e5), 1e5)

    Bx_dot = eta * g * dz

    _f = - grad_phi - xi -Bx_dot
    f = min(max(_f, -1e10), 1e10)

    pulled_f = J' * (f .- M * dJ * dx)[1]

    pulled_M = J' * M * J

    return pulled_f, pulled_M
end

"""実験データ"""
struct Data{T}
    t::StepRangeLen{T}
    q::Vector{Vector{T}}
    dq::Vector{Vector{T}}
    ddq::Vector{Vector{T}}
    desired_ddq::Vector{Vector{T}}
    u::Vector{Vector{T}}
    x1::Vector{Vector{T}}
    x2::Vector{Vector{T}}
    x3::Vector{Vector{T}}
    x4::Vector{Vector{T}}
    dx1::Vector{Vector{T}}
    dx2::Vector{Vector{T}}
    dx3::Vector{Vector{T}}
    dx4::Vector{Vector{T}}
    J1::Vector{Matrix{T}}
    J2::Vector{Matrix{T}}
    J3::Vector{Matrix{T}}
    J4::Vector{Matrix{T}}
    dJ1::Vector{Matrix{T}}
    dJ2::Vector{Matrix{T}}
    dJ3::Vector{Matrix{T}}
    dJ4::Vector{Matrix{T}}
    f0::Vector{Vector{T}}
    f1::Vector{Vector{T}}
    f2::Vector{Vector{T}}
    f3::Vector{Vector{T}}
    f4::Vector{Vector{T}}
    M0::Vector{Matrix{T}}
    M1::Vector{Matrix{T}}
    M2::Vector{Matrix{T}}
    M3::Vector{Matrix{T}}
    M4::Vector{Matrix{T}}
    error::Vector{T}
    min_dit_to_obs::Vector{T}
    jl::Vector{BitVector}
    F_distur::Vector{Vector{T}}  # トルクに入る外乱
    Fc::Vector{Vector{T}}  # 対象に与えた力（i.e. 対象が受ける力 = -Fc）
end

"""全部実行"""
function run_simulation(;
    TIME_SPAN::T=20.0, Δt::T=0.0001, isImpedance::Bool=false
    ) where T

    # 初期値
    q₀ = q_neutral
    dq₀ = zero(q₀)
    ddq₀ = zero(q₀)
    desired_ddq₀ = zero(q₀)
    u₀ = zero(q₀)
    F_distur₀ = zero(q₀)

    circle = (
        r = 0.5, x = 2.0, y = 0.5, K = 1.0, D = 1.0,
    )  # 円の物体の情報

    # 発生させたい力
    Fc_desired = 0.05  # [N]
    Δx = Fc_desired / circle.K

    # 目標値とコンプライアンス中心
    xd_true = [2.0, 1.0]

    if isImpedance
        xd = xd_true
        xe = xd .- [0.0, Δx]
    else
        xd = xd_true .- [0.0, Δx]
        xe = nothing
    end

    # 障害物
    # xo = [
    #     [1.5, 0.5],
    #     [2.5, 0.5]
    # ] #.* 1000000

    xo = [[2.0, 2.5]] .* 1000

    # 対象物
    # box_l = 2.0
    # box_h = 1.0
    # box_center = [2.0, 0.5]

    # rmpのパラメータ

    obs_avoidance = RMPfromGDSCollisionAvoidance(
        rw = 1.0,
        sigma = 1.0,
        alpha = 1.0,
    )

    jl_avoidance = RMPfromGDSJointLimitAvoidance(
        gamma_p = 0.01,
        gamma_d = 0.05,
        lambda = 1.0,
        sigma = 0.1,
    )

    if isImpedance
        # インピーダンス特性の決定
        zeta_d = 0.9  # 所望の減衰係数
        omega_d = 2.0
        dd = 4.0

        de = circle.D
        ke = circle.K

        md = (de + dd)/(2*omega_d*zeta_d)
        kd = (-ke*zeta_d + omega_d*(de + dd)/2)/zeta_d

        println("md = ", md)
        println("kd = ", kd)
        println("dd = ", dd)

        impedance = RMPfromGDSImpedance(
            M_d = Matrix{T}(I, 2, 2) * md,
            D_d = Matrix{T}(I, 2, 2) * dd,
            P_d = Matrix{T}(I, 2, 2) * kd,
            D_e = Matrix{T}(I, 2, 2) * de,
            P_e = Matrix{T}(I, 2, 2) * ke,
            a=1.0,
            eta_d=1.0,
            eta_e=1.0,
            f_alpha=0.15,
            sigma_alpha=1.0,
            sigma_gamma=1.0,
            wu=10.0,
            wl=0.1,
            alpha=0.15,
            epsilon=0.5,
        )

    else
        attractor = RMPfromGDSAttractor(
            max_speed = 3.3,
            gain = 10.0,
            f_alpha = 0.15,
            sigma_alpha = 1.0,
            sigma_gamma = 1.0,
            wu = 10.0,
            wl = 0.1,
            alpha = 0.15,
            epsilon = 0.5,
        )

    end

    t = range(0.0, TIME_SPAN, step=Δt)  # 時間列

    data = Data(
        t,
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Matrix{T}}(undef, length(t)),
        Vector{Matrix{T}}(undef, length(t)),
        Vector{Matrix{T}}(undef, length(t)),
        Vector{Matrix{T}}(undef, length(t)),
        Vector{Matrix{T}}(undef, length(t)),
        Vector{Matrix{T}}(undef, length(t)),
        Vector{Matrix{T}}(undef, length(t)),
        Vector{Matrix{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Matrix{T}}(undef, length(t)),
        Vector{Matrix{T}}(undef, length(t)),
        Vector{Matrix{T}}(undef, length(t)),
        Vector{Matrix{T}}(undef, length(t)),
        Vector{Matrix{T}}(undef, length(t)),
        Vector{T}(undef, length(t)),
        Vector{T}(undef, length(t)),
        Vector{BitVector}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
        Vector{Vector{T}}(undef, length(t)),
    )

    # 初期値代入
    data.q[1] = q₀
    data.dq[1] = dq₀
    data.ddq[1] = ddq₀
    data.desired_ddq[1] = desired_ddq₀
    data.u[1] = u₀

    data.x1[1], data.x2[1], data.x3[1], data.x4[1] = calc_x(q₀)
    data.J1[1], data.J2[1], data.J3[1], data.J4[1] = calc_J(q₀)
    data.dJ1[1], data.dJ2[1], data.dJ3[1], data.dJ4[1] = calc_dJ(q₀, dq₀)

    data.dx1[1] = data.J1[1] * data.dq[1]
    data.dx2[1] = data.J2[1] * data.dq[1]
    data.dx3[1] = data.J3[1] * data.dq[1]
    data.dx4[1] = data.J4[1] * data.dq[1]

    data.f0[1] = zero(q₀)
    data.f1[1] = zeros(T, 2)
    data.f2[1] = zeros(T, 2)
    data.f3[1] = zeros(T, 2)
    data.f4[1] = zeros(T, 2)

    data.M0[1] = zeros(T, 4, 4)
    data.M1[1] = zeros(T, 2, 2)
    data.M2[1] = zeros(T, 2, 2)
    data.M3[1] = zeros(T, 2, 2)
    data.M4[1] = zeros(T, 2, 2)

    #data.error[1] = norm(data.x4[1] - xd)
    data.error[1] = norm(data.x4[1] - [circle.x, circle.y]) - circle.r  # 円物体のとき

    data.min_dit_to_obs[1] = calc_min_dis_to_obs(
        [data.x1[1], data.x2[1], data.x3[1], data.x4[1]],
        xo
    )
    data.jl[1] = check_JointLimitation(q₀)
    data.F_distur[1] = zero(u₀)
    data.Fc[1] = Fc_from_cicle(
        data.x4[1], data.dx4[1],
        [circle.x, circle.y], circle.r, circle.K, circle.D
    )

    # ぐるぐる回す
    for i in 1:length(data.t)-1
        #println("t = ", data.t[i])

        # 状態を更新
        data.q[i+1] = data.q[i] .+ data.dq[i] .* Δt
        data.dq[i+1] = data.dq[i] .+ data.ddq[i] .*Δt

        data.x1[i+1], data.x2[i+1], data.x3[i+1], data.x4[i+1] = calc_x(data.q[i+1])
        data.J1[i+1], data.J2[i+1], data.J3[i+1], data.J4[i+1] = calc_J(data.q[i+1])
        data.dJ1[i+1], data.dJ2[i+1], data.dJ3[i+1], data.dJ4[i+1] = calc_dJ(data.q[i+1], data.dq[i+1])

        data.dx1[i+1] = data.J1[i+1] * data.dq[i+1]
        data.dx2[i+1] = data.J2[i+1] * data.dq[i+1]
        data.dx3[i+1] = data.J3[i+1] * data.dq[i+1]
        data.dx4[i+1] = data.J4[i+1] * data.dq[i+1]

        # rmpを計算

        # rmpをセット
        # root
        data.f0[i+1], data.M0[i+1] = get_natural(
            jl_avoidance, data.q[i+1], data.dq[i+1], q_max, q_min, q_neutral
        )
        data.f1[i+1] = zeros(T, 2)
        data.f2[i+1] = zeros(T, 2)
        data.f3[i+1] = zeros(T, 2)
        data.f4[i+1] = zeros(T, 2)
        data.M1[i+1] = zeros(T, 2, 2)
        data.M2[i+1] = zeros(T, 2, 2)
        data.M3[i+1] = zeros(T, 2, 2)
        data.M4[i+1] = zeros(T, 2, 2)

        for o in xo

            # multi robotのやつ
            f, M = multi_avoidance_rmp(data.x1[i+1], data.dx1[i+1], o)
            @. data.f1[i+1] += f
            @. data.M1[i+1] += M

            f, M = multi_avoidance_rmp(data.x2[i+1], data.dx2[i+1], o)
            @. data.f2[i+1] += f
            @. data.M2[i+1] += M

            f, M = multi_avoidance_rmp(data.x3[i+1], data.dx3[i+1], o)
            @. data.f3[i+1] += f
            @. data.M3[i+1] += M

            f, M = multi_avoidance_rmp(data.x4[i+1], data.dx4[i+1], o)
            @. data.f4[i+1] += f
            @. data.M4[i+1] += M

            # f, M = get_natural(obs_avoidance, data.x1[i+1], data.dx1[i+1], o)
            # @. data.f1[i+1] += f
            # @. data.M1[i+1] += M

            # f, M = get_natural(obs_avoidance, data.x2[i+1], data.dx2[i+1], o)
            # @. data.f2[i+1] += f
            # @. data.M2[i+1] += M

            # f, M = get_natural(obs_avoidance, data.x3[i+1], data.dx3[i+1], o)
            # @. data.f3[i+1] += f
            # @. data.M3[i+1] += M

            # f, M = get_natural(obs_avoidance, data.x4[i+1], data.dx4[i+1], o)
            # @. data.f4[i+1] += f
            # @. data.M4[i+1] += M
        end

        if isImpedance
            f, M = get_natural(impedance, data.x4[i+1], xd, xe, data.dx4[i+1], [circle.x, circle.y], circle.r)  # インピーダンス

            # 対象物表面との距離を測る
            n_surface = hat(data.x4[i+1] .- [circle.x, circle.y])  # 表面との法線ベクトル
            d_surface = norm(data.x4[i+1] .- [circle.x, circle.y]) - circle.r  # 表面との距離

            if d_surface < 0
                f_damp = zero(f)
                M_damp = zero(f)
            else
                #println("d = ", d_surface)
                a_damp = 0.8 / (d_surface/1.0 + 1e-3) * norm(data.dx4[i+1])
                f_damp = a_damp * n_surface
                #println("f = ", norm(f_damp))
                M_damp = Matrix{T}(I, 2, 2)

                @. f += f_damp
                @. M += M_damp
            end

        else
            f, M = get_natural(attractor, data.x4[i+1], data.dx4[i+1], xd)  # アトラクタ
        end
        @. data.f4[i+1] += f
        @. data.M4[i+1] += M
        #println("attractor = ", f)

        # #println("attractor = ", f)

        # pullback演算
        root_f = data.f0[i+1]
        root_M = data.M0[i+1]

        _rf, _rM = pullbacked_rmp(data.f1[i+1], data.M1[i+1], data.J1[i+1], data.dJ1[i+1], data.dq[i+1])
        @. root_f += _rf
        @. root_M += _rM

        _rf, _rM = pullbacked_rmp(data.f2[i+1], data.M2[i+1], data.J2[i+1], data.dJ2[i+1], data.dq[i+1])
        @. root_f += _rf
        @. root_M += _rM

        _rf, _rM = pullbacked_rmp(data.f3[i+1], data.M3[i+1], data.J3[i+1], data.dJ3[i+1], data.dq[i+1])
        @. root_f += _rf
        @. root_M += _rM

        _rf, _rM = pullbacked_rmp(data.f4[i+1], data.M4[i+1], data.J4[i+1], data.dJ4[i+1], data.dq[i+1])
        @. root_f += _rf
        @. root_M += _rM

        # resolve演算
        data.desired_ddq[i+1] = pinv(root_M) * root_f

        data.u[i+1] = calc_torque(data.q[i+1], data.dq[i+1], data.desired_ddq[i+1])

        data.F_distur[i+1] = zero(u₀)  # 外乱無し
        #data.F_distur[i+1] = rand(T, 4)*0.1

        data.Fc[i+1] = Fc_from_cicle(
            data.x4[i+1], data.dx4[i+1],
            [circle.x, circle.y], circle.r, circle.K, circle.D
        )

        data.ddq[i+1] = calc_real_ddq(
            u = data.u[i+1],
            q = data.q[i+1],
            dq = data.dq[i+1],
            F = data.F_distur[i+1],
            Fc = data.Fc[i+1],
            Jend = data.J4[i+1],
        )

        #data.error[i+1] = norm(data.x4[i+1] - xd)
        data.error[i+1] = norm(data.x4[i+1] - [circle.x, circle.y]) - circle.r  # 円物体のとき
        #println("error = ", data.error[i+1])
        data.min_dit_to_obs[i+1] = calc_min_dis_to_obs(
            [data.x1[i+1], data.x2[i+1], data.x3[i+1], data.x4[i+1]],
            xo,
        )
        data.jl[i+1] = check_JointLimitation(data.q[i+1])

    end

    # グラフ化

    x, y, z, w = split_vec_of_arrays(data.q)
    fig_q = plot(xlim=(0, TIME_SPAN), legend=:outerright)
    plot!(fig_q, data.t, x, label="_q1")
    plot!(fig_q, data.t, y, label="_q2")
    plot!(fig_q, data.t, z, label="_q3")
    plot!(fig_q, data.t, w, label="_q4")

    x, y, z, w = split_vec_of_arrays(data.dq)
    fig_dq = plot(xlim=(0, TIME_SPAN), legend=:outerright)
    plot!(fig_dq, data.t, x, label="_w1")
    plot!(fig_dq, data.t, y, label="_w2")
    plot!(fig_dq, data.t, z, label="_w3")
    plot!(fig_dq, data.t, w, label="_w4")

    x, y, z, w = split_vec_of_arrays(data.ddq)
    fig_ddq = plot(xlim=(0, TIME_SPAN), legend=:outerright)
    plot!(fig_ddq, data.t, x, label="_a1")
    plot!(fig_ddq, data.t, y, label="_a2")
    plot!(fig_ddq, data.t, z, label="_a3")
    plot!(fig_ddq, data.t, w, label="_a4")

    x, y, z, w = split_vec_of_arrays(data.desired_ddq)
    fig_desired_ddq = plot(xlim=(0, TIME_SPAN), legend=:outerright)
    plot!(fig_desired_ddq, data.t, x, label="a*1")
    plot!(fig_desired_ddq, data.t, y, label="a*2")
    plot!(fig_desired_ddq, data.t, z, label="a*3")
    plot!(fig_desired_ddq, data.t, w, label="a*4")

    x, y, z, w = split_vec_of_arrays(data.u)
    fig_u = plot(xlim=(0, TIME_SPAN), legend=:outerright)
    plot!(fig_u, data.t, x, label="_u1")
    plot!(fig_u, data.t, y, label="_u2")
    plot!(fig_u, data.t, z, label="_u3")
    plot!(fig_u, data.t, w, label="_u4")

    f0 = norm.(data.f0)
    f1 = norm.(data.f1)
    f2 = norm.(data.f2)
    f3 = norm.(data.f3)
    f4 = norm.(data.f4)
    _fs = [0.0]
    f_max = maximum(append!(_fs, f0, f1, f2, f3, f4))
    fig_f = plot(
        xlim=(0, TIME_SPAN), ylim=(0, f_max),
        legend=:outerright,
    )
    plot!(fig_f, data.t, f0, label="roo")
    plot!(fig_f, data.t, f1, label="_f1")
    plot!(fig_f, data.t, f2, label="_f2")
    plot!(fig_f, data.t, f3, label="_f3")
    plot!(fig_f, data.t, f4, label="_ee")

    fig_error = plot(
        data.t, data.error,
        label="err",
        xlim=(0, TIME_SPAN),
        #ylim=(0, maximum(data.error)),
        ylim=(minimum(data.error), 0.005),
        legend=:outerright,
    )

    # # 物体に与えた力をノルムでplot
    # Fcs = norm.(data.Fc)
    # fig_Fc = plot(
    #     data.t, Fcs,
    #     label="_Fc", ylabel="Fc",
    #     xlim=(0, TIME_SPAN), ylim=(0, maximum(Fcs)),
    #     legend=:outerright,
    # )

    x, y = split_vec_of_arrays(data.Fc)
    fig_Fc = plot(
        xlim=(0, TIME_SPAN),
        legend=:outerright,
    )
    plot!(fig_Fc, data.t, norm.(data.Fc), label="nFc")
    plot!(fig_Fc, data.t, x, label="Fcx")
    plot!(fig_Fc, data.t, y, label="Fcy")
    plot!(fig_Fc, ylim=(-0.06, 0.08))

    s = Vector{T}(undef, length(t))
    ds = Vector{T}(undef, length(t))
    for i in 1:length(t)
        s[i] = norm(data.x4[i] .- [circle.x, circle.y]) - circle.r
        ds[i] = 1/s[i] * dot((data.x4[i] .- [circle.x, circle.y]), data.dx4[i])
    end

    fig_s = plot(
        data.t, s, label="__s",
        xlim=(0, TIME_SPAN),
        legend=:outerright,
    )
    fig_ds = plot(
        data.t, ds, label="_ds",
        xlim=(0, TIME_SPAN),
        legend=:outerright,
    )

    fig_ee_dx = plot(
        data.t, norm.(data.dx4),
        label="dx4",
        xlim=(0, TIME_SPAN),
        legend=:outerright,
    )
    x, y = split_vec_of_arrays(data.dx4)
    plot!(fig_ee_dx, data.t, x, label="_dx")
    plot!(fig_ee_dx, data.t, y, label="_dy")

    fig_dis_to_obs = plot(
        data.t, data.min_dit_to_obs,
        label="obs",
        xlim=(0, TIME_SPAN), ylim=(0, maximum(data.min_dit_to_obs)),
        legend=:outerright,
    )

    fig = plot(
        fig_q, fig_dq, fig_ddq,
        fig_desired_ddq, fig_u, fig_f,
        fig_error, fig_Fc, fig_s, fig_ds, fig_ee_dx,
        fig_dis_to_obs,
        layout=(6, 2),
        size=(1200, 2000)
    )

    if isImpedance
        name = "fig/sice_simple_proposed.png"
    else
        name = "fig/sice_simple_conventional.png"
    end

    savefig(fig, name)

    # アニメ作成

    # 少し準備
    """箱の形"""
    function rectangle(w, h, x, y)
        Shape(x .+ [-w/2,w/2,w/2,-w/2], y .+ [-h/2,-h/2,h/2,h/2])
    end

    """楕円"""
    function ellipse(x, y, w, h)
        ang = range(0, 2π, length = 60)
        Shape(w*sin.(ang).+x, h*cos.(ang).+y)
    end

    """1フレームを描写"""
    function draw_frame(i)
        arm = [
            [0.0, 0.0],
            data.x1[i],
            data.x2[i],
            data.x3[i],
            data.x4[i]
        ]

        x, y = split_vec_of_arrays(arm)
        fig = plot(
            x, y,
            marker=:circle,
            aspect_ratio = 1,
            xlabel = "X[m]", ylabel = "Y[m]",
        )

        scatter!(
            fig,
            [xd_true[1]], [xd_true[2]],
            label="xd_true",
            markershape=:star6,
        )

        scatter!(
            fig,
            [xd[1]], [xd[2]],
            label="xd_true",
            markershape=:star6,
        )

        if !isnothing(xe)
            scatter!(
                fig,
                [xe[1]], [xe[2]],
                label="xe",
                markershape=:star6,
            )
        end

        x, y = split_vec_of_arrays(xo)
        scatter!(
            fig,
            x, y
        )

        # 対象物
        #plot!(rectangle(box_l, box_h, box_center[1], box_center[2]), opacity=.5)
        plot!(ellipse(circle.x, circle.y, circle.r, circle.r), opacity=.5)

        x_max = 4.1
        x_min = -4.1
        y_max = 4.1
        y_min = -4.1
        max_range = max(x_max-x_min, y_max-y_min)*0.5
        x_mid = (x_max + x_min) / 2
        y_mid = (y_max + y_min) / 2

        plot!(
            fig,
            xlims=(x_mid-max_range, x_mid+max_range),
            ylims=(y_mid-max_range, y_mid+max_range),
            legend = false,
            size=(600, 600),
        )

        # rmpも描画

        scale = 10.0

        f = [data.x1[i], data.x1[i] .+ data.f1[i] .* scale]
        x, y = split_vec_of_arrays(f)
        plot!(
            fig,
            x, y
        )

        f = [data.x2[i], data.x2[i] .+ data.f2[i] .* scale]
        x, y = split_vec_of_arrays(f)
        plot!(
            fig,
            x, y
        )

        f = [data.x3[i], data.x3[i] .+ data.f3[i] .* scale]
        x, y = split_vec_of_arrays(f)
        plot!(
            fig,
            x, y
        )

        f = [data.x4[i], data.x4[i] .+ data.f4[i] .* scale]
        x, y = split_vec_of_arrays(f)
        plot!(
            fig,
            x, y
        )

        # function make_metric_cicle(i, M, q_global)
        #     e = eigen(M)
        #     axis_len = 1 / sqrt.(e.Vector) * 0.1
        #     max_len = maximum(axis_len)
        #     scale = minimum(2/max_len, 1)
        #     M_xis_len = axis_len .* scale
        #     angle = atan(e.Vector[2], e.Vector[1]) + q_global |>
        #     rad2deg

        # end

        fig_2 = plot(fig, xlim=(1.8, 2.2), ylim=(0.8, 1.2))

        fig_i = plot(
            fig, fig_2,
            layout=(2,1),
            size=(600,1200)
        )

        id = findall(x->x==true, data.jl[i])
        plot!(fig_i, title = string(round(data.t[i], digits=2)) * "[s]" * ", " * string(id))
        return fig_i
    end

    # println("アニメ作成中...")
    # # 枚数決める
    # #println(data.t)
    # epoch_max = 100
    # epoch = length(data.t)
    # if epoch < epoch_max
    #     step = 1
    # else
    #     step = div(epoch, epoch_max)
    # end

    # #println(step)

    # anim = Animation()
    # @gif for i in 1:step:length(data.q)
    #     _fig = draw_frame(i)
    #     frame(anim, _fig)
    # end

    # if isImpedance
    #     name = "fig/sice_animation_proposed.gif"
    # else
    #     name = "fig/sice_animation_conventional.gif"
    # end

    # gif(anim, name, fps = 60)

    # println("アニメ作成完了")

    ### sice用のグラフ作成 ###
    fig_error = plot(
        data.t, data.error,
        label=L"||x - x_g||",
        xlabel=L"Time [s]",
        ylabel=L"Position Error [m]",
        xlim=(0, TIME_SPAN),
        ylim=(0.00, maximum(data.error)),
    )
    plot!(size=(400, 200))

    savefig(fig_error, "error.png")

    x, y = split_vec_of_arrays(data.Fc)
    fig_Fc = plot(
        xlim=(0, TIME_SPAN),
    )
    plot!(fig_Fc, data.t, norm.(data.Fc), label=L"F_c")
    plot!(ylabel = L"Applied Force [N]")
    plot!(xlabel = L"Time [s]")
    plot!(fig_Fc, ylim=(0.0, 0.10))
    plot!(size=(400, 200))

    savefig(fig_Fc, "Fc.png")

    return data
    # fig_q, fig_dq, fig_ddq,
    # fig_desired_ddq, fig_u, fig_f,
    # fig_error, fig_Fc, fig_dis_to_obs
end

println("実行中...")
# @time data,
# fig_q, fig_dq, fig_ddq,
# fig_desired_ddq, fig_u, fig_f,
# fig_error, fig_Fc, fig_dis_to_obs = run_simulation()

@time data_pro = run_simulation(isImpedance=true)
@time data_con = run_simulation(isImpedance=false)

println("実行終了!!")


text:
"""
    DofHandler(grid::Grid)

Construct a `DofHandler` based on the grid `grid`.
"""
struct DofHandler{dim,T}
    variables::Vector{TrialFunction}
    cell_dofs::Vector{Int}
    cell_dofs_offset::Vector{Int}
    mesh::AbstractPolytopalMesh{dim,T}
end

function DofHandler(mesh::AbstractPolytopalMesh{dim,T}, trialfuncs::Vector{TrialFunction{dim,T}}) where {dim,T}
    dofhandler = DofHandler{dim,T}(trialfuncs, Int[], Int[], mesh)
    _distribute_dofs(dofhandler)
end

function Base.show(io::IO, dh::DofHandler)
    println(io, "DofHandler")
    println(io, "  Number of Fields: ",nvariables(dh))
    println(io, "  Dofs per cell: ", ndofs_per_cell(dh))
    print(io, "  Total dofs: ", ndofs(dh))
end

ndofs(dh::DofHandler) = maximum(dh.cell_dofs)
ndofs_per_cell(dh::DofHandler, cell::Int=1) = dh.cell_dofs_offset[cell+1] - dh.cell_dofs_offset[cell]
@inline nvariables(dh::DofHandler) = length(dh.variables)

function find_field(dh::DofHandler, field::TrialFunction)
    j = findfirst(i->i == field, dh.variables)
    j == 0 && error("did not find field $field")
    return j
end

# Calculate the offset to the first local dof of a field
function field_offset(dh::DofHandler, field::TrialFunction, ci::Int)
    offset = 0
    for i in 1:find_field(dh, field)-1
        offset += getnlocaldofs(getfunctionspace(dh.variables[i]), getcells(dh.cells)[ci])::Int
    end
    return offset
end

"""
    dof_range(dh::DofHandler, u::Symbol)

Return the local dof range for `u`.
"""
function dof_range(dh::DofHandler, field::TrialFunction, ci::Int)
    f = find_field(dh, field)
    offset = field_offset(dh, field, ci)
    n_field_dofs = getnlocaldofs(getfunctionspace(dh.variables[f]),getcells(dh.mesh)[ci])
    return (offset+1):(offset+n_field_dofs)
end

# close the DofHandler and distribute all the dofs
#TODO: This only work for D < 3
function _distribute_dofs(dh::DofHandler{dim,T}) where {dim,T}
    # `vertexdict` keeps track of the visited vertices. We store the global vertex
    # number and the first dof we added to that vertex.
    vertexdicts = [Dict{entityeltype(dh.mesh,0),Int}() for _ in 1:nvariables(dh)]

    # `edgedict` keeps track of the visited edges, this will only be used for a 3D problem
    # We also need to store the direction
    # of the first edge we encounter and add dofs too. When we encounter the same edge
    # the next time we check if the direction is the same, otherwise we reuse the dofs
    # in the reverse order
    edgedicts = [Dict{entityeltype(dh.mesh,1),Int}() for _ in 1:nvariables(dh)]

    # `facedict` keeps track of the visited faces. We only need to store the first dof we
    # added to the face; if we encounter the same face again we *always* reverse the order
    # In 2D a face (i.e. a line) is uniquely determined by 2 vertices, and in 3D a
    # face (i.e. a surface) is uniquely determined by 3 vertices.
    facedicts = dim < 3 ? [] : [Dict{entityeltype(dh.mesh,2),Int}() for _ in 1:nvariables(dh)]

    topologyDicts = dim == 2 ? Dict(0=>vertexdicts, 1=>edgedicts) : Dict(0=>vertexdicts, 1=>edgedicts, 2 => facedicts)

    nextdof = 1 # next free dof to distribute
    push!(dh.cell_dofs_offset, 1) # dofs for the first cell start at 1

    # loop over all the cells, and distribute dofs for all the fields
    for (ci, cell) in enumerate(getcells(dh.mesh))
        # Get topologies
        n_max_topology_elements = maximum(keys(gettopology(dh.mesh,cell)))
        geometric_cell_topology = gettopology(dh.mesh,cell)
        for fi in 1:nvariables(dh)
            element = getelement(getfunctionspace(dh.variables[fi]))
            cell_topology = gettopology(dh.mesh, cell, element)
            for n_element in 0:n_max_topology_elements-1
                n_el_dofs_cell = cell_topology[n_element]
                n_el_cell = geometric_cell_topology[n_element]
                @assert mod(n_el_dofs_cell, n_el_cell) == 0
                nelementdofs = Int(n_el_dofs_cell/n_el_cell)
                if nelementdofs > 0
                    for element in getcellsubentities(dh.mesh,ci,n_element)
                        token = ht_keyindex2!(topologyDicts[n_element][fi], element)
                        if token > 0 # reuse dofs
                            reuse_dof = topologyDicts[n_element][fi].vals[token]
                            for d in 1:getncomponents(getfunctionspace(dh.variables[fi]))
                                push!(dh.cell_dofs, reuse_dof + (d-1))
                            end
                        else # token <= 0, distribute new dofs
                            for elementdof in 1:nelementdofs
                                Base._setindex!(topologyDicts[n_element][fi], nextdof, element, -token)
                                for d in 1:getncomponents(getfunctionspace(dh.variables[fi]))
                                    push!(dh.cell_dofs, nextdof)
                                    nextdof += 1
                                end
                            end
                        end
                    end # element loop
                end
            end #elements loop
            ncelldofs = cell_topology[dim]
            if ncelldofs > 0 # always distribute new dofs for cell
                for celldof in 1:ncelldofs
                    for d in 1:getncomponents(getfunctionspace(dh.variables[fi]))
                        push!(dh.cell_dofs, nextdof)
                        nextdof += 1
                    end
                end # cell loop
            end
        end # field loop
        # push! the first index of the next cell to the offset vector
        push!(dh.cell_dofs_offset, length(dh.cell_dofs)+1)
    end # cell loop
    return dh
end

function celldofs!(global_dofs::Vector{Int}, dh::DofHandler, i::Int)
    @assert length(global_dofs) == ndofs_per_cell(dh, i)
    unsafe_copyto!(global_dofs, 1, dh.cell_dofs, dh.cell_dofs_offset[i], length(global_dofs))
    return global_dofs
end

function vertexdofs(dh::DofHandler, u::TrialFunction)
    element = getelement(getfunctionspace(u))
    dofs = zeros(Int,getnvertices(dh.mesh))
    for (ci, cell) in enumerate(getcells(dh.mesh))
        nv = gettopology(dh.mesh,cell, element)[0]*getncomponents(getfunctionspace(u))
        for i in dh.cell_dofs_offset[ci]:(dh.cell_dofs_offset[ci]+nv-1)
            dof = dh.cell_dofs[i]
            if !(dof ∈ dofs)
                dofs[getverticesindices(dh.mesh,cell)[i-dh.cell_dofs_offset[ci]+1]] = dof
            end
        end
    end
    return dofs
end

# Creates a sparsity pattern from the dofs in a DofHandler.
# Returns a sparse matrix with the correct storage pattern.
"""
    create_sparsity_pattern(dh::DofHandler)

Create the sparsity pattern corresponding to the degree of freedom
numbering in the [`DofHandler`](@ref). Return a `SparseMatrixCSC`
with stored values in the correct places.

See the [Sparsity Pattern](@ref) section of the manual.
"""
@inline create_sparsity_pattern(dh::DofHandler) = _create_sparsity_pattern(dh, false)

"""
    create_symmetric_sparsity_pattern(dh::DofHandler)

Create the symmetric sparsity pattern corresponding to the degree of freedom
numbering in the [`DofHandler`](@ref) by only considering the upper
triangle of the matrix. Return a `Symmetric{SparseMatrixCSC}`.

See the [Sparsity Pattern](@ref) section of the manual.
"""
@inline create_symmetric_sparsity_pattern(dh::DofHandler) = Symmetric(_create_sparsity_pattern(dh, true), :U)

function _create_sparsity_pattern(dh::DofHandler, sym::Bool)
    ncells = getncells(dh.mesh)
    n = ndofs_per_cell(dh)
    N = sym ? div(n*(n+1), 2) * ncells : n^2 * ncells
    N += ndofs(dh) # always add the diagonal elements
    I = Int[]; resize!(I, N)
    J = Int[]; resize!(J, N)
    global_dofs = zeros(Int, n)
    cnt = 0
    for element_id in 1:ncells
        celldofs!(global_dofs, dh, element_id)
        @inbounds for j in 1:n, i in 1:n
            dofi = global_dofs[i]
            dofj = global_dofs[j]
            sym && (dofi > dofj && continue)
            cnt += 1
            if cnt > length(J)
                resize!(I, trunc(Int, length(I) * 1.5))
                resize!(J, trunc(Int, length(J) * 1.5))
            end
            I[cnt] = dofi
            J[cnt] = dofj

        end
    end
    @inbounds for d in 1:ndofs(dh)
        cnt += 1
        if cnt > length(J)
            resize!(I, trunc(Int, length(I) + ndofs(dh)))
            resize!(J, trunc(Int, length(J) + ndofs(dh)))
        end
        I[cnt] = d
        J[cnt] = d
    end
    resize!(I, cnt)
    resize!(J, cnt)
    V = zeros(length(I))
    K = sparse(I, J, V)
    return K
end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools