text:
module ManifoldConstructors

using LinearAlgebra
using Random
using SparseArrays
using StaticArrays

using ..Manifolds
using ..Meshing
using ..SparseOps
using ..ZeroOrOne

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

export empty_manifold
"""
The empty manifold
"""
function empty_manifold(::Val{D}, ::Type{S}; options...) where {D,S}
    return Manifold("empty manifold", zero(SparseOp{0,D,One}, 0, 0),
                    IDVector{0,SVector{D,S}}(), IDVector{0,S}(); options...)
end

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

export simplex_manifold
"""
Manifold with one standard simplex
"""
function simplex_manifold(::Val{D}, ::Type{S}; options...) where {D,S}
    nvertices = D + 1
    coords = IDVector{0}(regular_simplex(Val(D), S))
    weights = IDVector{0}(zeros(S, length(coords)))
    simplices0 = MakeSparse{One}(nvertices, 1)
    for i in 1:nvertices
        simplices0[i, 1] = One()
    end
    simplices = SparseOp{0,D,One}(sparse(simplices0))
    return Manifold("simplex manifold", simplices, coords, weights; options...)
end

"""
Generate the coordinate positions for a regular D-simplex with edge
length 1.

The algorithm proceeds recursively. A 0-simplex is a point. A
D-simplex is a (D-1)-simplex that is shifted down along the new axis,
plus a new point on the new axis.
"""
function regular_simplex(::Val{D}, ::Type{S}) where {D,S}
    @assert D ≥ 0
    N = D + 1
    s = Array{SVector{D,S}}(undef, N)
    if D > 0
        s0 = regular_simplex(Val(D - 1), S)
        # Choose height so that edge length is 1
        if D == 1
            z = S(1)
        else
            z0 = sqrt(1 - sum(s0[1] .^ 2))
            if S <: Rational
                z = rationalize(typeof(zero(S).den), z0; tol=sqrt(eps(z0)))
            else
                z = z0::S
            end
        end
        z0 = -z / (D + 1)
        for n in 1:(N - 1)
            s[n] = SVector{D,S}(s0[n]..., z0)
        end
        s[N] = SVector{D,S}(zero(SVector{D - 1,S})..., z + z0)
    end
    return s
end

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

export orthogonal_simplex_manifold
"""
Manifold with one orthogonal simplex
"""
function orthogonal_simplex_manifold(::Val{D}, ::Type{S};
                                     options...) where {D,S}
    N = D + 1
    coords, weights = orthogonal_simplex(Val(D), S)
    coords = IDVector{0}(coords)
    weights = IDVector{0}(weights)
    simplices0 = MakeSparse{One}(N, 1)
    for i in 1:N
        simplices0[i, 1] = One()
    end
    simplices = SparseOp{0,D,One}(sparse(simplices0))
    return Manifold("orthogonal simplex manifold", simplices, coords, weights;
                    options...)
end

"""
Generate the coordinate positions for an orthogonal D-simplex with
edge lengths 1.
"""
function orthogonal_simplex(::Val{D}, ::Type{S}) where {D,S}
    @assert D ≥ 0
    s = zeros(SVector{D,S}, D + 1)
    w = zeros(S, D + 1)
    w[1] = Dict(0 => 0, 1 => 0, 2 => -S(1) / 2,
                # 2 => -S(1) / 8,
                3 => -S(1) / 2, 4 => -S(2) / 3, 5 => -S(2) / 3)[D]
    w[1] = -1 + S(1) / D
    for d in 1:D
        s[d + 1] = setindex(zero(SVector{D,S}), 1, d)
    end
    return s, w
end

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

export hypercube_manifold
"""
Standard tesselation of a hypercube
"""
function hypercube_manifold(::Val{D}, ::Type{S}; options...) where {D,S}
    @assert D ≥ 0
    N = D + 1

    # Find simplices
    simplices = SVector{N,Int}[]
    corner = zero(SVector{D,Bool})
    vertices = [corner2vertex(corner)]
    next_corner!(simplices, vertices, corner)
    nsimplices = length(simplices)
    @assert nsimplices == factorial(D)

    # Set up coordinates
    coords = IDVector{0,SVector{D,S}}()
    imin = CartesianIndex(ntuple(d -> 0, D))
    imax = CartesianIndex(ntuple(d -> 1, D))
    for i in imin:imax
        push!(coords, SVector{D,S}(i.I))
    end
    nvertices = length(coords)
    @assert nvertices == 2^D
    weights = IDVector{0}(zeros(S, nvertices))

    simplices′ = MakeSparse{One}(nvertices, nsimplices)
    for (j, sj) in enumerate(simplices)
        for i in sj
            simplices′[Int(i), Int(j)] = One()
        end
    end
    simplices″ = SparseOp{0,D,One}(sparse(simplices′))

    return Manifold("hypercube manifold", simplices″, coords, weights;
                    options...)
end

"""
- Accumulate the simplices in `simplices`.
- `vertices` is the current set of vertices as we sweep from the
  origin to diagonally opposide vertex.
- `corner` is the current corner.
"""
function next_corner!(simplices::Vector{SVector{N,Int}}, vertices::Vector{Int},
                      corner::SVector{D,Bool})::Nothing where {N,D}
    @assert N == D + 1
    if D == 0
        @assert length(vertices) == 1
    end
    @assert count(corner) == length(vertices) - 1
    if length(vertices) == D + 1
        # We have all vertices; build the simplex
        push!(simplices, SVector{N}(vertices))
        return
    end
    # Loop over all neighbouring corners
    for d in 1:D
        if !corner[d]
            new_corner = setindex(corner, true, d)
            new_vertex = corner2vertex(new_corner)
            new_vertices = [vertices; new_vertex]
            next_corner!(simplices, new_vertices, new_corner)
        end
    end
end

function corner2vertex(c::SVector{D,Bool})::Int where {D}
    D == 0 && return 1
    return 1 + sum(c[d] << (d - 1) for d in 1:D)
end

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

export large_hypercube_manifold
"""
Triangulate a large hypercube by splitting it into smaller hypercubes
"""
function large_hypercube_manifold(::Val{D}, ::Type{S};
                                  nelts::Union{Nothing,NTuple{D,Int}}=nothing,
                                  torus::Bool=false, options...) where {D,S}
    @assert D ≥ 0
    N = D + 1

    if nelts ≡ nothing
        nelts = ntuple(d -> 2, D)
    end

    @assert all(nelts .≥ 0)
    # We could relax this
    @assert all(nelts .% 2 .== 0)

    symmetric_hypercube = create_hypercube(Val(D))
    for dir in 1:D
        append!(symmetric_hypercube, mirror_hypercube(symmetric_hypercube, dir))
    end
    @assert length(symmetric_hypercube) == factorial(D) * 2^D

    stride = SVector{D,Int}(d == 1 ? 1 : prod(nelts[1:(d - 1)] .+ 1)
                            for d in 1:D)
    function vertex2ind(i::SVector{D,Int})
        if !torus
            return 1 + sum(i .* stride)
        else
            return 1 + sum(i .% nelts .* stride)
        end
    end

    nsimplices = factorial(D) * (D == 0 ? 1 : prod(nelts))
    simplices = Array{SVector{N,Int}}(undef, nsimplices)
    i = 0
    for elt in CartesianIndices(nelts .÷ 2)
        offset = 2 * (SVector{D,Int}(elt.I) .- 1)
        for s in symmetric_hypercube
            simplices[i += 1] = map(v -> vertex2ind(v .+ offset), s)
        end
    end
    @assert i == length(simplices)

    nvertices = D == 0 ? 1 : prod(nelts .+ 1)
    coords0 = Array{SVector{D,S}}(undef, nvertices)
    i = 0
    for elt in CartesianIndices(nelts .+ 1)
        coords0[i += 1] = (SVector{D,S}(elt.I) .- 1) ./ nelts
    end
    @assert i == length(coords0)
    coords = IDVector{0}(coords0)

    # fact = SVector{D,S}(a == D ? 1 / S(2) : 1 for a in 1:D)
    # map!(x -> fact .* x, coords, coords)

    if false
        weights = IDVector{0}(zeros(S, nvertices))
    else
        # These weights put the weighted circumcentre onto the barycentre
        levelweights = D == 0 ? S[0] :
                       D == 1 ? S[0, 0] :
                       D == 2 ? S[0, -S(1) / 6, 0] :
                       D == 3 ? S[0, -S(1) / 6, -S(1) / 6, 0] :
                       D == 4 ? S[0, -S(3) / 20, -S(1) / 5, -S(3) / 20, 0] :
                       D == 5 ?
                       S[0, -S(2) / 15, -S(1) / 5, -S(1) / 5, -S(2) / 15, 0] :
                       nothing
        if D > 0
            @assert all(==(nelts[1]), nelts)
            levelweights ./= nelts[1]^2
        end

        weights = IDVector{0}(Vector{S}(undef, nvertices))
        for vi in CartesianIndex(ntuple(d -> 0, D)):CartesianIndex(nelts)
            v = SVector{D,Int}(vi.I)
            l = D == 0 ? 0 : sum(v[d] % 2 for d in 1:D)
            weights[ID{0}(vertex2ind(v))] = levelweights[l + 1]
        end
    end

    simplices′ = MakeSparse{One}(nvertices, nsimplices)
    for (j, sj) in enumerate(simplices)
        for i in sj
            simplices′[Int(i), Int(j)] = One()
        end
    end
    simplices″ = SparseOp{0,D,One}(sparse(simplices′))

    return Manifold("large hypercube manifold", simplices″, coords, weights;
                    options...)
end

function create_hypercube(::Val{D}) where {D}
    @assert D ≥ 0
    N = D + 1

    # Find simplices
    simplices = SVector{N,SVector{D,Int}}[]
    corner = zero(SVector{D,Bool})
    vertices = SVector{D,Int}[corner]
    next_corner!(simplices, vertices, corner)
    nsimplices = length(simplices)
    @assert nsimplices == factorial(D)

    # Ensure we have all 2^D vertices
    cvertices = SVector{D,Int}[]
    for si in simplices
        for j in si
            push!(cvertices, j)
        end
    end
    unique!(sort!(cvertices))
    @assert length(cvertices) == 2^D

    return simplices
end

"""
- Accumulate the simplices in `simplices`.
- `vertices` is the current set of vertices as we sweep from the
  origin to diagonally opposide vertex.
- `corner` is the current corner.
"""
function next_corner!(simplices::Vector{SVector{N,SVector{D,Int}}},
                      vertices::Vector{SVector{D,Int}},
                      corner::SVector{D,Bool}) where {N,D}
    @assert N == D + 1
    if D == 0
        @assert length(vertices) == 1
    end
    @assert sum(corner) == length(vertices) - 1
    if length(vertices) == D + 1
        # We have all vertices; build the simplex
        push!(simplices, SVector{N}(vertices))
        return
    end
    # Loop over all neighbouring corners
    for d in 1:D
        if !corner[d]
            new_corner = setindex(corner, true, d)
            new_vertex = SVector{D,Int}(new_corner)
            new_vertices = copy(vertices)
            push!(new_vertices, new_vertex)
            next_corner!(simplices, new_vertices, new_corner)
        end
    end
    return
end

function mirror_hypercube(simplices::Vector{SVector{N,SVector{D,Int}}},
                          dir::Int) where {D,N}
    @assert D ≥ 0
    @assert N == D + 1
    @assert 1 ≤ dir ≤ D

    flip = SVector{D,Int}(d == dir ? -1 : 1 for d in 1:D)
    offset = SVector{D,Int}(d == dir ? 2 : 0 for d in 1:D)
    simplices′ = map(s -> map(v -> flip .* v .+ offset, s), simplices)

    return simplices′
end

function shift_hypercube(simplices::Vector{SVector{N,SVector{D,Int}}},
                         offset::SVector{D,Int}) where {D,N}
    @assert D ≥ 0
    @assert N == D + 1

    simplices′ = map(s -> map(v -> v .+ offset, s), simplices)

    return simplices′
end

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

export delaunay_hypercube_manifold
"""
Delaunay triangulation of a hypercube
"""
function delaunay_hypercube_manifold(::Val{D}, ::Type{S};
                                     options...) where {D,S}
    @assert D ≥ 0
    N = D + 1

    # Set up coordinates
    coords = IDVector{0,SVector{D,S}}()
    imin = CartesianIndex(ntuple(d -> 0, D))
    imax = CartesianIndex(ntuple(d -> 1, D))
    for i in imin:imax
        # Re-map coordinates to avoid round-off errors
        push!(coords, SVector{D,S}(i.I) .+ S(5) / 2)
    end
    nvertices = length(coords)
    @assert nvertices == 2^D
    weights = IDVector{0}(zeros(S, length(coords)))

    simplices = delaunay_mesh(coords)

    return Manifold("delaunay hypercube manifold", simplices, coords, weights;
                    options...)
end

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

export large_delaunay_hypercube_manifold
"""
Delaunay triangulation of a large hypercube
"""
function large_delaunay_hypercube_manifold(::Val{D}, ::Type{S};
                                           nelts::Union{Nothing,Int}=nothing,
                                           options...) where {D,S}
    @assert D ≥ 0
    N = D + 1

    if nelts ≡ nothing
        ns = Dict{Int,Int}(0 => 1, 1 => 1024, 2 => 32, 3 => 16, 4 => 4, 5 => 2)
        nelts = ns[D]
    end
    # All n should be powers of 2 to avoid round-off below
    @assert ispow2(nelts)

    rng = MersenneTwister(1)

    # Set up coordinates
    coords = IDVector{0,SVector{D,S}}()
    # coords′ = SVector{D,S}[]
    imin = CartesianIndex(ntuple(d -> 0, D))
    imax = CartesianIndex(ntuple(d -> nelts, D))
    for i in imin:imax

        # # We also set up slightly different coordinates to help the
        # # Delaunay triangulation
        # # 1. Perturb points in the interior randomly to avoid
        # #    degeneracies
        # # 2. Make points on the boundary bulge outwards to avoid
        # #    degeneracies
        #
        # x = SVector{D,S}(i[d] for d in 1:D) / n
        #
        # # # Avoid floating-point round-off errors
        # # x = x .+ S(5) / 2
        #
        # # Add noise
        # isbnd(d) = i[d] ∈ (0, n)
        # x += SVector{D,S}(isbnd(d) ? 0 : (2 * rand(rng, S) - 1) / (256 * n)
        #                   for d in 1:D)
        #
        # # These are the coordinates we want
        # push!(coords, x)
        #
        # # x′ = x
        # #
        # # # Barrel transformation
        # # r = norm(2 * x′ .- 1)
        # # x′ = ((2 * x′ .- 1) * (1 - r / (256^3 * n)) .+ 1) / 2
        # #
        # # push!(coords′, x′)

        # Stagger every second point, except boundary points
        x = SVector{D,S}(i[d] for d in 1:D) / nelts
        for d in 1:D
            for d′ in 1:(d - 1)
                i[d′] ∈ (0, nelts) && continue
                di = (S(1) / 2) / sqrt(S(d - 1))
                x1 = x[d′]
                x1 = x1 * nelts / (nelts + di)
                if isodd(i[d])
                    x1 = x1 + di / (nelts + di)
                end
                x = Base.setindex(x, x1, d′)
            end
        end

        push!(coords, x)
    end
    # @assert length(coords′) == length(coords)
    nvertices = length(coords)
    @assert nvertices == (nelts + 1)^D
    weights = IDVector{0}(zeros(S, length(coords)))

    simplices = delaunay_mesh(coords)

    return Manifold("large delaunay hypercube manifold", simplices, coords,
                    weights; options...)
end

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

export random_hypercube_manifold
"""
Hypercube manifold with random vertices
"""
function random_hypercube_manifold(xmin::SVector{D,S}, xmax::SVector{D,S},
                                   npoints_per_dim::Int; options...) where {D,S}
    rng = MersenneTwister(1)

    function coord(f, d)
        f == -1 && return xmin[d]
        f == +1 && return xmax[d]
        local n = 2^20
        local i = rand(rng, 1:(n - 1))
        local x = (n - i) / S(n) * xmin[d] + i / S(n) * xmax[d]
        return x
    end

    coords = IDVector{0}(Vector{SVector{D,S}}())
    for face in
        CartesianIndex(ntuple(d -> -1, D)):CartesianIndex(ntuple(d -> 1, D))
        dim = count(==(0), face[d] for d in 1:D)
        for i in 1:(npoints_per_dim^dim)
            x = SVector{D,S}(coord(face[d], d) for d in 1:D)
            push!(coords, x)
        end
    end
    weights = IDVector{0}(zeros(S, length(coords)))

    simplices = delaunay_mesh(coords)

    return Manifold("random hypercube manifold", simplices, coords, weights;
                    options...)
end

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

export refined_manifold
"""
Refined manifold
"""
function refined_manifold(mfd0::Manifold{D,C,S}; options...) where {D,C,S}
    D == 0 && return mfd0
    coords = refine_coords(get_lookup(mfd0, 0, 1), get_coords(mfd0))
    weights = IDVector{0}(zeros(S, length(coords)))
    simplices = delaunay_mesh(coords)
    mfd = Manifold("refined $(mfd0.name)", simplices, coords, weights;
                   options...)
    return mfd
end

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

export refined_simplex_manifold
"""
Refined simplex manifold
"""
function refined_simplex_manifold(::Val{D}, ::Type{S}; options...) where {D,S}
    mfd0 = simplex_manifold(Val(D), S; options...)
    mfd = refined_manifold(mfd0; options...)
    return mfd
end

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

export boundary_manifold
"""
Boundar manifold
"""
function boundary_manifold(mfd0::Manifold{D,C,S}; options...) where {D,C,S}
    @assert D > 0
    # Find boundary faces
    B = get_boundaries(mfd0, D)
    volume_weights = IDVector{D}(ones(Int8, nsimplices(mfd0, D)))
    keep_face = B * volume_weights
    keep_face::IDVector{D - 1,Int8}
    @assert all(x -> abs(x) ≤ 1, keep_face)
    faces_old2new = IDVector{D - 1}(zeros(ID{D - 1}, nsimplices(mfd0, D - 1)))
    faces_new2old = IDVector{D - 1,ID{D - 1}}()
    for j in axes(get_simplices(mfd0, D - 1), 2)
        if keep_face[j] ≠ 0
            push!(faces_new2old, j)
            faces_old2new[j] = ID{D - 1}(length(faces_new2old))
        end
    end
    nfaces_new = length(faces_new2old)
    # Find boundary vertices
    F = get_simplices(mfd0, D - 1)
    keep_vertex = IDVector{0}(zeros(Bool, nsimplices(mfd0, 0)))
    for j in faces_new2old
        for i in sparse_column_rows(F, j)
            keep_vertex[i] = true
        end
    end
    vertices_old2new = IDVector{0}(zeros(ID{0}, nsimplices(mfd0, 0)))
    vertices_new2old = IDVector{0,ID{0}}()
    for i in axes(get_simplices(mfd0, 0), 2)
        if keep_vertex[i]
            push!(vertices_new2old, i)
            vertices_old2new[i] = ID{0}(length(vertices_new2old))
        end
    end
    nvertices_new = length(vertices_new2old)
    # Define new manifold
    name = "boundary $(mfd0.name)"
    simplices′ = MakeSparse{One}(nvertices_new, nfaces_new)
    for jnew in axes(faces_new2old, 1)
        jold = faces_new2old[jnew]
        for iold in sparse_column_rows(F, jold)
            inew = vertices_old2new[iold]
            simplices′[Int(inew), Int(jnew)] = One()
        end
    end
    simplices″ = SparseOp{0,D - 1,One}(sparse(simplices′))
    coords = get_coords(mfd0)[vertices_new2old]
    weights = get_weights(mfd0)[vertices_new2old]
    mfd = Manifold(name, simplices″, coords, weights; options...)
    return mfd
end

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

export boundary_simplex_manifold
"""
Boundary simplex manifold
"""
function boundary_simplex_manifold(::Val{D}, ::Type{S}; options...) where {D,S}
    mfd0 = simplex_manifold(Val(D + 1), S; options...)
    mfd = boundary_manifold(mfd0; options...)
    return mfd
end

end


text:
# max x₁+x₂ s.t.
#  x₁ + 2x₂ ≤ 4
# 4x₁ + 2x₂ ≤ 12
# -x₁ +  x₂ ≤ 1
# and xᵢ ≥ 0 ∀ i
a = [
     1  2  1  0  0  4
     4  2  0  1  0 12
    -1  1  0  0  1  1
     1  1  0  0  0  0
]
# should give 10/3 with x = (8/3, 2/3)

simplex_max(a::Array{<:Integer}) =
    simplex_max(Rational.(a))

simplex_max(a) = _simplex(a, findmax, x -> x .> 0)
simplex_min(a) = _simplex(a, findmin, x -> x .< 0)

"""
solve a linear problem with the basic simplex algorithm
TODO doesn't handle unbounded solution or empty solution spaces
TODO currently the initial matrix must represent a feasible solution
    => use simplex_dualpass
# Arguments
- `a`: the problem as a matrix
"""
function _simplex(a, objective, filter)
    nrow,ncol = size(a) .- 1
    base = collect(ncol-nrow+1:ncol)
    # assume qu'on a l'identité à droite du tableau, donc la base c'est juste
    # les n dernières variables avec n = nb de contraintes (oui on assume aussi
    # que c'est que des contraintes ≤)
    while any(filter(getcosts(a)))
        _,new,old = simplex_pass!(a, objective, filter)
        base[old] = new # maj la base courante
        println("colonne #$new rentre, pivot sur ligne #$old, nouvelle base=$base")
    end
    println("Coût final = $(-a[end])")
    values = a[1:end-1,end]
    solution = zeros(eltype(a), ncol)
    solution[base] = values
    print("Solution = ") ; flush(STDOUT)
    display(solution)
    -a[end], base, values, a
end

getcosts(a) = a[end,1:end-1]
debug(obj::AbstractString) = print(obj)
debug(obj) = (flush(STDOUT); display(obj); println())
debug(objs...) = for obj in objs debug(obj) end
indf(indfunc, arr, f) = indfunc(f(e)?e:typemax(e) for e in arr)
indminpos(itr) = indmin(e>0?e:NaN for e in itr)

simplex_pass(a, objective, filter) = simplex_pass!(copy(a), objective, filter)
# marche avec findmax TODO tester avec findmin
function simplex_pass!(a::Array{<:Real}, heuristic, filter)
    costs = getcosts(a)
    extr, idx = findmax(costs) ; @assert filter(extr) # variable à entrer dans la base
    p = a[1:end-1,end] ./ a[1:end-1,idx]        # pour trouver le pivot
    ip = indf(indmin, p, filter)                # pivot, ⟹ variable à sortir
    a[ip,:] ./= a[ip,idx]                       # normalise la ligne du pivot
    temp = a[:,idx] * a[ip,:]'                  # on fait entrer la nouvelle variable
    temp[ip,:] = 0                              #
    a .-= temp                                  #
    a, idx, ip
end

function simplex_dualpass!(a, f)
    a,idx,ip = simplex_pass!(a', f)
    a',idx,ip
end

@time x,zup,zop = simplex_max(a)
println()
@time simplex_min(a)

ext: jl

max_stars_repo_name: artemisart/julia-simplex

max_issues_repo_name: artemisart/julia-simplex

max_forks_repo_name: artemisart/julia-simplex


text:
# SPDX-License-Identifier: MIT

#======================================================================
This file wrapped a few basic linear algebra operations in order to
have more fine-grained control between CPU/GPU code and different
semifield algebra operations.
======================================================================#

## Element-wise division.

eldiv!(out, x::AbstractArray, y::AbstractArray) = broadcast!(/, out, x, y)

## Element-wise multiplication.

elmul!(out, x::AbstractArray, y::AbstractArray) =  broadcast!(*, out, x, y)
elmul!(out, x::CuSparseVector, y::AbstractVector) = elmul_svdv!(out, x, y)
elmul!(out, x::CuSparseVector, Y::AbstractMatrix) = elmul_svdm!(out, x, Y)
function elmul!(out, x::AbstractSparseArray{T1}, y::AbstractArray{T2}) where {T1,T2}
    T = promote_type(T1, T2)
    fill!(out, zero(T))
    broadcast!(*, out, x, y)
end
elmul!(out, x::AbstractSparseMatrix{T1}, y::AbstractArray{T2}) where {T1,T2} = _elmul_sm_dm!(out, x, y)
elmul!(out, x::AbstractSparseMatrix{T1}, y::AbstractSparseMatrix{T2}) where {T1,T2} = _elmul_sm_dm!(out, x, y)
elmul!(out, x::AbstractArray{T1}, y::AbstractSparseMatrix{T2}) where {T1,T2} = _elmul_sm_dm!(out, y, x)

function _elmul_sm_dm!(out, X::AbstractSparseMatrix{T1}, y::AbstractArray{T2}) where {T1,T2}
    T = promote_type(T1, T2)
    fill!(out, zero(T))

    if sparsity_level(X) > 0.1
        return broadcast!(*, out, X, y)
    end

    m, n = size(X)
    length(y) in [m, n, m*n] || throw(DimensionMismatch("arrays could not be broadcast to a common size; got a matrix ($m, $n) and $(length(y))"))
    rows = rowvals(X)
    vals = nonzeros(X)
    for j = 1:n
       for i in nzrange(X, j)
          row = rows[i]
          val = vals[i]
          out[row, j] = val * _getindex(y, row, j)
       end
    end
    return out
end

sparsity_level(x::AbstractSparseArray{T}) where T = begin
    nnz(x) / reduce(*, size(x))
end
_getindex(x::AbstractMatrix{T}, i, j) where T = begin
    if size(x, 1) == 1
        x[1, j]
    elseif size(x, 2) == 1
        x[i, 1]
    else
        x[i, j]
    end
end
_getindex(x::AbstractVector{T}, i, j) where T = x[i]

## Matrix multiplication.

function matmul!(out, X::AbstractMatrix{T1}, y::AbstractArray{T2}) where {T1,T2}
    T = promote_type(T1, T2)
    mul!(out, X, y, one(T), zero(T))
end

function matmul!(out, X::AbstractSparseMatrix{T1}, y::AbstractArray{T2}) where {T1,T2}
    T = promote_type(T1, T2)
    fill!(out, zero(T))
    mul!(out, X, y, one(T), zero(T))
end

matmul!(out, X::CuMatrix, Y::CuMatrix) = mul_dmdm!(out, X, Y)
matmul!(out, X::CuSparseMatrixCSR, y::AbstractVector) = mul_smdv!(out, X, y)
matmul!(out, X::CuSparseMatrixCSR, Y::AbstractMatrix) = mul_smdm!(out, X, Y)

max_stars_repo_name: MartinKocour/HiddenMarkovModel

max_issues_repo_name: MartinKocour/HiddenMarkovModel

max_forks_repo_name: MartinKocour/HiddenMarkovModel


text:
using Documenter, EAGOBranchBound

makedocs(modules=[EAGOBranchBound],
         doctest=false, clean=true,
         format =:markdown,
         sitename="EAGO.jl",
         authors="Matthew Wilhelm",
         pages = Any[
         "Home" => "index.md",
         "Tutorials" => Any[
           "tutorials/page1.md",
           "tutorials/page2.md",
           "tutorials/page3.md"
         ],
         "EAGOBranchBound" => Any[
           "userdocs/page1.md",
         ]
         ])


text:
"""
    Element(session::Session, location_strategy::AbstractString, value::AbstractString)::Element
    Element(element::Element, location_strategy::AbstractString, value::AbstractString)::Element

The Find Element command is used to find an element in the current browsing context that can be used as the web element context for future element-centric commands.

May use a location strategy from "css selector", "link text", "partial link text", "tag name" or "xpath".

```jldoctest
julia> capabilities = Capabilities("chrome")
Remote WebDriver Capabilities
browserName: chrome
julia> wd = RemoteWebDriver(capabilities, host = ENV["WEBDRIVER_HOST"], port = parse(Int, ENV["WEBDRIVER_PORT"]))
Remote WebDriver
julia> session = Session(wd)
Session
julia> isa(session, Session)
true
julia> delete!(session);
```
"""
struct Element
    session::Session
    id::String
    Element(session::Session, id::AbstractString) = new(session, id)
    # Return Element
    Element(element::Element, id::AbstractString) = new(element.session, id)
    # Command: Find Element
    function Element(
        session::Session,
        location_strategy::AbstractString,
        value::AbstractString,
    )::Element
        location_strategy ∈
        ["css selector", "link text", "partial link text", "tag name", "xpath"] ||
        throw(ArgumentError("""location_strategy must be "css selector", "link text", "partial link text", "tag name" or "xpath"."""))
        @unpack addr, id = session
        # location_strategy = "xpath"
        # value = "//select[@id='selecttype']"
        response = HTTP.post(
            "$addr/session/$id/element",
            [("Content-Type" => "application/json; charset=utf-8")],
            JSON3.write(Dict("using" => location_strategy, "value" => value)),
        )
        @assert response.status == 200
        new(session, first(values(JSON3.read(response.body).value)))
    end
    # Command: Find Element from Element
    function Element(
        element::Element,
        location_strategy::AbstractString,
        value::AbstractString,
    )::Element
        location_strategy ∈
        ["css selector", "link text", "partial link text", "tag name", "xpath"] ||
        throw(ArgumentError("""location_strategy must be "css selector", "link text", "partial link text", "tag name" or "xpath"."""))
        @unpack addr, id = element.session
        element_id = element.id
        response = HTTP.post(
            "$addr/session/$id/element/$element_id/element",
            [("Content-Type" => "application/json; charset=utf-8")],
            JSON3.write(Dict("using" => location_strategy, "value" => value)),
        )
        @assert response.status == 200
        new(element.session, first(values(JSON3.read(response.body).value)))
    end
end
broadcastable(obj::Element) = Ref(obj)
summary(io::IO, obj::Element) = println(io, "Element")
function show(io::IO, obj::Element)
    print(io, summary(obj))
end


text:
using PyPlot

function plot(vector::Array{Float64, 1}, dofmap::DofMap, filename::String)
    mesh = dofmap.mesh

    if dofmap.element.cell() == ufl.triangle

        # Prepare handles to vertex coords
        xs = mesh.vertices[:, 1]
        ys = mesh.vertices[:, 2]

        gca(aspect="equal")
        # -1 in topology for Python indexing
        triplot(xs, ys, mesh.topology[(2,0)] .- 1, color="black", linewidth=0.1)
        ax = tripcolor(xs, ys, mesh.topology[(2,0)] .- 1, vector, shading="gouraud")
        colorbar(ax)
        savefig(filename)
    else
        throw("Unsupported dofmap")
    end
end


text:
using JuMP, GLPK
using DataFrames, CSV
using LinearAlgebra
using JSON
using LightGraphs

#cd("/home/xiaotong/Dropbox (MIT)/13_Amazon_LM_competition/submission_new")
cd(@__DIR__)

tsp_opt_seq = JSON.parsefile("../data/model_apply_outputs/model_apply_output/mean_dist/opt_complete_seq_tour_post_process.json")

# Output and store sequences
model_apply_output = Dict()
for i in keys(tsp_opt_seq)
    model_apply_output[i] = Dict()
    model_apply_output[i]["proposed"] = Dict()
    for j in 1:length(tsp_opt_seq[i])
        model_apply_output[i]["proposed"][tsp_opt_seq[i][j]] = j-1
    end
end

open("../data/model_apply_outputs/proposed_sequences.json","w") do f
    JSON.print(f, model_apply_output)
end

println("Successfully store proposed sequences")


text:
#=

  Four Islands puzzle (Dell Logic Puzzles) in ConstraintSolver.jl

  http://brownbuffalo.sourceforge.net/FourIslandsClues.html
  """
  Title: Four Islands
  Author: Humphrey Dudley
  Publication: Dell Logic Puzzles
  Issue: April, 1998
  Page: 9
  Stars: 1

  A tiny nation in the South Pacific contains four islands connected by bridges
  as shown (see below). Each of the four islands (Pwana, Quero, Rayou, and Skern)
  boasts a different primary export (alabaster, bananas, coconuts, and durian
  fruit) and a different tourist attraction (hotel, ice skating rink, jai alai
  stadium, and koala preserve). Can you find the name, export, and tourist
  attraction of each island on the map?

    N
  W   E     *compass directions
    S

  A, B, C, D are the islands

  (A) -- (B)
   |      |
   |      |
  (C) -- (D)

  1. The island noted for its koala preserve is due south of Pwana.
  2. The island with the largest alabaster quarry is due west of Quero.
  3. The island with the resort hotel is due east of the one that exports
     durian fruit.
  4. Skern and the island with the jai alai stadium are connected by a
     north-south bridge.
  5. Rayou and the island that exports bananas are connected by an east-west
     bridge.
  6. The islands noted for the South Pacific's largest ice skating rink and
     for the jai alai stadium are not connected by a bridge.

  Determine: Island island -- Island name -- Export -- Tourist Attraction
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function four_islands(print_solutions=true,all_solutions=true,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    n = 4
    A,B,C,D = 1:4
    @variable(model, 1 <= Island[1:n] <= n, Int)
    Pwana, Quero, Rayou, Skern = Island

    @variable(model, 1 <= Export[1:n] <= n, Int)
    Alabaster, Bananas, _Coconuts, DurianFruit = Export

    @variable(model, 1 <= Attraction[1:n] <= n, Int)
    ResortHotel, IceSkatingRink, JaiAlaiStadium, KoalaPreserve = Attraction

    @constraint(model, Island in CS.AllDifferent())
    @constraint(model, Export in CS.AllDifferent())
    @constraint(model, Attraction in CS.AllDifferent())

    # 1. The island noted for its koala preserve is due south of Pwana.
    @constraint(model, Pwana == A && KoalaPreserve == C ||
                       Pwana == B && KoalaPreserve == D)

    # 2. The island with the largest alabaster quarry is due west of Quero.
    @constraint(model, Alabaster == A && Quero == B ||
                       Alabaster == C && Quero == D)

    # 3. The island with the resort hotel is due east of the one
    #    that exports durian fruit.
    @constraint(model, DurianFruit == A && ResortHotel == B  ||
                       DurianFruit == C && ResortHotel == D)

    # 4. Skern the island with the jai alai stadium are connected by a
    #    north-south bridge.
    @constraint(model, Skern == A && JaiAlaiStadium == C ||
                       Skern == C && JaiAlaiStadium == A ||
                       Skern == B && JaiAlaiStadium == D ||
                       Skern == D && JaiAlaiStadium == B)

    # 5. Rayou/\the island that exports bananas are connected by an
    #    east-west bridge.
    @constraint(model, (Rayou == A && Bananas == B) ||
                       (Rayou == B && Bananas == A) ||
                       (Rayou == C && Bananas == D) ||
                       (Rayou == D && Bananas == C))

    # 6. The islands noted for the South Pacific's largest ice skating rink
    #    for the jai alai stadium are not connected by a bridge.
    @constraint(model, IceSkatingRink == A && JaiAlaiStadium == D ||
                       IceSkatingRink == D && JaiAlaiStadium == A ||
                       IceSkatingRink == B && JaiAlaiStadium == C ||
                       IceSkatingRink == C && JaiAlaiStadium == B)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                island_val = convert.(Integer,JuMP.value.(Island; result=sol))
                export_val = convert.(Integer,JuMP.value.(Export; result=sol))
                attraction_val = convert.(Integer,JuMP.value.(Attraction; result=sol))
                println("island    :$island_val")
                println("export    :$export_val")
                println("attraction:$attraction_val")

            end
        end
    else
        println("status:$status")
    end

    return status, num_sols
end

@time four_islands()

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct Ipv6CircuitConnectionConfig <: SwaggerModel
    addressPrefix::Any # spec type: Union{ Nothing, String } # spec name: addressPrefix
    circuitConnectionStatus::Any # spec type: Union{ Nothing, CircuitConnectionStatus } # spec name: circuitConnectionStatus

    function Ipv6CircuitConnectionConfig(;addressPrefix=nothing, circuitConnectionStatus=nothing)
        o = new()
        validate_property(Ipv6CircuitConnectionConfig, Symbol("addressPrefix"), addressPrefix)
        setfield!(o, Symbol("addressPrefix"), addressPrefix)
        validate_property(Ipv6CircuitConnectionConfig, Symbol("circuitConnectionStatus"), circuitConnectionStatus)
        setfield!(o, Symbol("circuitConnectionStatus"), circuitConnectionStatus)
        o
    end
end # type Ipv6CircuitConnectionConfig

const _property_map_Ipv6CircuitConnectionConfig = Dict{Symbol,Symbol}(Symbol("addressPrefix")=>Symbol("addressPrefix"), Symbol("circuitConnectionStatus")=>Symbol("circuitConnectionStatus"))
const _property_types_Ipv6CircuitConnectionConfig = Dict{Symbol,String}(Symbol("addressPrefix")=>"String", Symbol("circuitConnectionStatus")=>"CircuitConnectionStatus")
Base.propertynames(::Type{ Ipv6CircuitConnectionConfig }) = collect(keys(_property_map_Ipv6CircuitConnectionConfig))
Swagger.property_type(::Type{ Ipv6CircuitConnectionConfig }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_Ipv6CircuitConnectionConfig[name]))}
Swagger.field_name(::Type{ Ipv6CircuitConnectionConfig }, property_name::Symbol) =  _property_map_Ipv6CircuitConnectionConfig[property_name]

function check_required(o::Ipv6CircuitConnectionConfig)
    true
end

function validate_property(::Type{ Ipv6CircuitConnectionConfig }, name::Symbol, val)
end


text:
using Pkg
Pkg.activate("rl-test")
includet("identity-skip.jl")
includet("code.jl")

using ReinforcementLearningEnvironments, ReinforcementLearningBase, Flux
using StatsBase, Statistics
const γ = 1 # discount rate
const α = 0.0001 # learning rate

function update_params!(ps, state, action, reward)
    grad = gradient(()->reward*log(eps(Float32) + policy(state)[action]), ps)
    for p in ps
        Flux.update!(p, -α*grad[p])
    end
end

env = CartPoleEnv(max_steps=1000)

action_space = get_action_space(env)

policy = Chain(
    Dense(4, 16, relu),
    Dense(16, 16, relu),
    Dense(16, length(action_space)),
    softmax
)

ps = Flux.params(policy)

obs = observe(env)
policy(obs.state)
update_params!(ps, obs.state, 1, 1)
policy(obs.state)
update_params!(ps, obs.state, 1, -1)
policy(obs.state)
update_params!(ps, obs.state, 2, 1)
policy(obs.state)
update_params!(ps, obs.state, 2, -1)

function ok(n)
    all_states = []
    all_rewards = []
    all_actions = []
    epsilon = 1.0
    epsilon_min = 0.01
    epsilon_decay = 0.995
    for i in 1:n
        policy(obs.state)
        states, rewards, actions = rollout(policy; epsilon=epsilon, max_steps=1000)

        length(rewards)

        # train the network
        for i in length(rewards)-1:-1:1
            rewards[i] += rewards[i+1]*γ
        end

        all_actions = vcat(all_actions, actions)
        all_rewards = vcat(all_rewards, rewards)
        all_states = vcat(all_states, states)

        epsilon = max(epsilon_min, epsilon*epsilon_decay)
    end

    all_rewards .= (all_rewards .- mean(all_rewards))./std(all_rewards)

    # randomise the order
    s = sample(1:length(all_rewards), length(all_rewards), replace=false)

    all_states .= all_states[s]
    all_rewards .= all_rewards[s]
    all_actions .= all_actions[s]

    for (s, a, r) in zip(all_states, all_actions, all_rewards)
        update_params!(ps, s, a, r)
    end

    # play it  to show how well it does
    # [length(rollout(policy; epsilon=0.0, max_steps=1000)[3]) for j in 1:1000] |> mean |> println
    # [length(rollout(policy; epsilon=0.0, max_steps=1000)[3]) for j in 1:1000] |> extrema |> println

end

function playone()
    env = CartPoleEnv(max_steps=1000)
    action_space = get_action_space(env) |> collect
    while true
        local action
        _, action = findmax(policy(env.state))
        env(action)
        render(env)
        if get_terminal(observe(env))
            break
        end
        sleep(0.001)
    end
end

@time ok(50)

while true
    @time ok(50)
    if rand() < 0.1
        playone()
    end
end

max_stars_repo_name: xiaodaigh/my-rl-gym

max_issues_repo_name: xiaodaigh/my-rl-gym

max_forks_repo_name: xiaodaigh/my-rl-gym


text:
#=

  Golomb ruler in Julia ConstraintSolver.jl

  A Golomb ruler is a set of integers (marks) a(1) < ...  < a(n) such
  that all the differences a(i)-a(j) (i > j) are distinct.  Clearly we
  may assume a(1)=0.  Then a(n) is the length of the Golomb ruler.
  For a given number of marks, n, we are interested in finding the
  shortest Golomb rulers.  Such rulers are called optimal.

  See http://www.research.ibm.com/people/s/shearer/grule.html

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function golomb_ruler(n=6,print_solutions=true,all_solutions=true,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            # "traverse_strategy"=>:BFS,
                                                            "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 6, # default, 10
                                                            # "activity.max_confidence_deviation" => 10, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=> timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    # m = 2^(n-1) -1
    m = n^2
    @variable(model, 0 <= x[1:n] <= m, Int)

    @constraint(model, x[1] == 0)
    @constraint(model, x in CS.AllDifferent())
    increasing(model,x)

    # This don't work:
    # """Each variable must be an integer and bounded."""
    # diffs = [x[i]-x[j] for i in 1:n, j in 1:n if i < j]
    # println("diffs: $diffs")
    # This works:
    diffs = []
    for i in 1:n, j in 1:n
        if i != j
            d = @variable(model,integer=true, lower_bound = -m, upper_bound = m )
            @constraint(model, d == x[i]-x[j])
            push!(diffs,d)
        end
    end

    @constraint(model, diffs in CS.AllDifferent())

    # Symmetry breaking
    if n > 2
        @constraint(model, x[2] - x[1] < x[n] - x[n-1] )
    end
    @constraint(model, diffs[1] < diffs[end])

    @objective(model,Min,x[n])

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        # println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                diffs_val = convert.(Integer,JuMP.value.(diffs; result=sol))
                println("x:$x_val")
                println("diffs:$diffs_val")
                println()

            end
        end
    else
        println("status:$status")
    end

    return status
end # end golomb_ruler

# @time golomb_ruler(8,true,false) # 1.5

# n  time
# -------
#  8   1.3s
#  9  13.3s
# 10  > 60s
for n in 2:10
    println("\nn:$n")
    @time status = golomb_ruler(n,true,false,60)
    if status == MOI.TIME_LIMIT
        break
    end
end


text:
function nlp_mi_004_011(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - non-linear objective and linear, quadratic constraints as NL

    model = Model(optimizer)

    @variable(model, -1 <= x <= 1)
    @variable(model, y)
    @variable(model, z, Int)

    @NLobjective(model, Min, tan(x) + y + x * z + 0.5 * abs(y))
    @NLconstraint(model, x^2 + y^2 + z^2 <= 10)
    @NLconstraint(model, -1.2 * x - y <= z / 1.35)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, -4.5768813091909015, tol = objective_tol)
    return check_solution(
        [x, y, z],
        [-0.9969558935432885, -0.07796881265504599, 3],
        tol = primal_tol,
    )
end


text:
using DictDotDot
using Test

@testset "Set" begin
    @testset "Set single property" begin
        d = Dict()
        DictDotDot.dotset!(d, "hello", "world")
        DictDotDot.dotset!(d, "hi", Dict("im" => "DictDotDot"))
        @test d["hello"] == "world"
        @test d["hi"]["im"] == "DictDotDot"
    end
end

max_stars_repo_name: nampdn/DictDotDot

max_issues_repo_name: nampdn/DictDotDot

max_forks_repo_name: nampdn/DictDotDot


text:
using TightlyBound
using Plots

#setup chosen Plots backend, i recommend pyplot
#pyplot()
#gr()

#make the crystal object
#here we choose fcc Al

types=["Al"];

#positions, crystal units
pos=zeros((1,3));

#lattice vectors, in Angstrom by default
A=[ [4.0 4.0 0]; [4.0 0 4.0 ]; [ 0 4.0 4.0]] * 0.529177;

#makes the crystal
c=makecrys(A, pos, types)

println("Initial Crystal:")
println()
println(c)
println()

energy, forces, stress, tbc = scf_energy_force_stress(c)

println("Initial Forces--")
println(TightlyBound.Utility.arr2str(forces))
println()
println("Initial Stress--")
println(TightlyBound.Utility.arr2str(stress))
println()

cfinal, tbc, energy, forces, stress = relax_structure(c)

println("Final crystal:")
println(cfinal)
println()
println("Final Forces--")
println(TightlyBound.Utility.arr2str(forces))
println()
println("Finial Stress--")
println(TightlyBound.Utility.arr2str(stress))
println()

#plot band structure
#kpath=[0.0 0.0 0.0;0.5 0 0; 0.5 0.25 0.75; 0.0 0.0 0.0]
#knames=["Γ", "X", "W", "Γ"]
#plot_bandstr(tbc, kpath=kpath, names=knames)


text:
mutable struct Particle
    values :: Array{Float64,1}
    currentValue :: Float64
    velocity :: Array{Float64,1}
    bestValues :: Array{Float64,1}
    bestValue :: Float64
end

function updateValue!(particle)
    x = particle.values[1];
    y = particle.values[2];
    z = particle.values[3];
    w = particle.values[4];
    u = particle.values[5];
    v = particle.values[6];
    particle.currentValue = 4*x^2 - 6*x - 3*y^3 + 0.5*y + 3*z + 8*w - 6.1*u + 2*v -10;
end

function generateParticle(numberOfOperands,minRang,maxRang)
    particle = Particle(rand(minRang:maxRang,numberOfOperands),0,
                        fill(0.0, numberOfOperands),fill(0,numberOfOperands),0);
    updateValue!(particle);
    particle.bestValues = copy(particle.values);
    particle.bestValue = particle.currentValue;
    return particle;
end

function printParticle(particle)
    x = particle.values[1];
    y = particle.values[2];
    z = particle.values[3];
    w = particle.values[4];
    u = particle.values[5];
    v = particle.values[6];
    print("4*$x^2-6*$x-3*$y^3+0.5*$y+3*$z+8*$w-6.1*$u+2*$v-10\n");
    print(particle.currentValue,"\n");
end

# cv*iv + cp*rp*(pi-xi) + cg*rg*(g-xi)
# Cv = 1 Cp = 2 Cg = 2
function updateVelocity!(particle, globalBestParticle, maxVelocity)
    for i in 1:length(particle.velocity)
        particle.velocity[i] += 2.0*particle.velocity[i] + 2.0*rand(0:1)*(particle.bestValues[i]
                             - particle.values[i]) + 2.0*rand(0:1)*(globalBestParticle.values[i]
                             - particle.values[i]);
        if particle.velocity[i] > maxVelocity
            particle.velocity[i] = maxVelocity;
        elseif particle.velocity[i] <  -maxVelocity
            particle.velocity[i] = -maxVelocity;
        end
    end
end

function updatePosition!(particle, targetValue)
    for i in 1:length(particle.values)
        particle.values[i] += particle.velocity[i];
        if particle.values[i] > 1
            particle.values[i] = 1
        elseif particle.values[i] < -1
            particle.values[i] = -1;
        end
    end
    updateValue!(particle);
    if abs(particle.currentValue - targetValue) < abs(particle.bestValue - targetValue)
        particle.bestValues = copy(particle.bestValues);
        particle.bestValue = particle.currentValue;
    end
end

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
"""
    value_at_risk(returns, confidence, method; multiplier=1.0)

Computes the Value-at-Risk (VaR) for a given significance level `α` based on the chosen estimation method. The VaR value represents the maximum expected loss at a certain significance level `α`. For a more tail-risk focused measure, see `expected_shortfall`.

# Arguments
- `returns`:     Vector of asset returns.
- `α`:           Significance level, e.g. use `0.05` for 95% confidence, or `0.01` for 99% confidence.
- `method`:      Distribution estimation method: `:historical`, `:gaussian` or `:cornish_fisher`.
- `multiplier`:  Optional scalar multiplier, i.e. use `12` to annualize monthly returns, and use `252` to annualize daily returns.

# Methods
- `:historical`:        Historical based on empirical distribution of returns.
- `:gaussian`:          Gaussian distribution based on parametric fit (mean, variance).
- `:cornish_fisher`:    Cornish-Fisher based on Gaussian parametric distribution fit adjusted for third and fourth moments (skewness, kurtosis). Cornish-Fisher expansion aims to approximate the quantile of a true distribution by using higher moments (skewness and kurtosis) of that distribution to adjust for its non-normality. See https://thema.u-cergy.fr/IMG/pdf/2017-21.pdf for details.

# Sources
- Favre, Laurent and Galeano, Jose-Antonio (2002). Mean-Modified Value-at-Risk Optimization with Hedge Funds. Journal of Alternative Investment.
- Amédée-Manesme, Charles-Olivier and Barthélémy, Fabrice and Maillard, Didier (2017). Computation of the Corrected Cornish–Fisher Expansion using the Response Surface Methodology: Application to VaR and CVaR. THEMA Working Paper n°2017-21, Université de Cergy-Pontoise, France.
"""
function value_at_risk(returns, α, method::Symbol; multiplier=1.0)
    if method == :historical
        # empirical quantile for VaR estimation
        return quantile(returns, α) * sqrt(multiplier)
    elseif method == :gaussian
        # parametric Gaussian distribution fit
        μ = mean(returns)
        σ = std(returns; corrected=false)
        return quantile(Normal(μ, σ), α)
    elseif method == :cornish_fisher
        # third/fourth moment adjusted Gaussian distribution fit
        # http://www.diva-portal.org/smash/get/diva2:442078/FULLTEXT01.pdf
        # https://papers.ssrn.com/sol3/papers.cfm?abstract_id=1024151
        q = quantile(Normal(), α)
        S = skewness(returns)
        K = kurtosis(returns; method=:excess)
        z = q + 1/6*(q^2-1)S + 1/24*(q^3-3q)*K - 1/36*(2q^3-5q)*S^2
        μ = mean(returns)
        σ = std(returns; corrected=false)
        return (μ + z*σ) * sqrt(multiplier)
    end

    throw(ArgumentError("Passed method parameter '$(method)' is invalid, must be one of :historical, :gaussian, :cornish_fisher."))
end


text:
# visualize various lines

function drawLine!(vispath, from::Vector{Float64}, to::Vector{Float64}; scale=0.01,color=RGBA(0,1.0,0,0.5))
  vector = to-from
  len = norm(vector)
  buildline = Float64[len, 0, 0]

  v = norm(buildline-vector) > 1e-10 ? LinearAlgebra.cross(buildline, vector)  : [0,0,1.0]
  axis = v/norm(v)
  angle = acos(dot(vector, buildline)/(len^2) )
  rot = LinearMap( CoordinateTransformations.AngleAxis(angle, axis...) )

  mol = HyperRectangle(Vec(0.0,-scale,-scale), Vec(len,scale,scale))
  molbox = GeometryData(mol, color)

  setgeometry!(vispath, molbox)
  settransform!(vispath, Translation(from...) ∘ rot )
  nothing
end


text:
using CSV, DataFrames, RCall

data1 = CSV.read("../data/machine1.csv", header=false)[:,1]
data2 = CSV.read("../data/machine2.csv", header=false)[:,1]
data3 = CSV.read("../data/machine3.csv", header=false)[:,1]

function R_ANOVA(allData)
    data = vcat([ [x fill(i, length(x))] for (i, x) in
    			enumerate(allData) ]...)
    df = DataFrame(data, [:Diameter, :MachNo])
    @rput df

    R"""
    df$MachNo <- as.factor(df$MachNo)
    anova <- summary(aov( Diameter ~ MachNo, data=df))
    fVal <- anova[[1]]["F value"][[1]][1]
    pVal <- anova[[1]]["Pr(>F)"][[1]][1]
    """
    println("R ANOVA f-value: ", @rget fVal)
    println("R ANOVA p-value: ", @rget pVal)
end

R_ANOVA([data1, data2, data3])

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using WAV
using PyPlot

y, fs = wavread("sample.wav")
plot(y, ".")
fy = fft(y)
plot(fy)


text:
using CLI

add(x::Int, y::Int) = x + y
mul(x::Int, y::Int) = x * y

leaf1 = LeafCommand(:add, add, Signature([Int, Int], []))
leaf2 = LeafCommand(:mul, mul, Signature([Int, Int], []))
node = NodeCommand("math", [leaf1, leaf2])

maincmd = MainCommand("cal", [], [node])

stream = ARGStream(ARGS)
resolve(stream, maincmd; debug=false)


text:
export
    Camera,
    StaticCamera,
    FitToContentCamera,
    CarFollowCamera,
    SceneFollowCamera

abstract type Camera end
camera_set!(::RenderModel, cam::Camera, ::EntityFrame{S,D,I}, ::R, canvas_width::Int, canvas_height::Int) where {S,D,I,R} = error("camera_set! not implemented for Camera $cam")

mutable struct StaticCamera <: Camera
    pos::VecE2
    zoom::Float64 # [pix/meter]
    StaticCamera(pos::VecE2, zoom::Float64=3.0) = new(pos, zoom)
    StaticCamera(x::Float64, y::Float64, zoom::Float64=3.0) = new(VecE2(x,y), zoom)
end
function camera_set!(rendermodel::RenderModel, cam::StaticCamera, canvas_width::Int, canvas_height::Int)

    camera_set_pos!(rendermodel, cam.pos)
    camera_setzoom!(rendermodel, cam.zoom)

    rendermodel
end
camera_set!(rendermodel::RenderModel, cam::StaticCamera, scene::EntityFrame{S,D,I}, roadway::R, canvas_width::Int, canvas_height::Int) where {S,D,I,R} = camera_set!(rendermodel, cam, canvas_width, canvas_height)

# method for new interface
camera_set!(rm::RenderModel, cam::StaticCamera, scene, canvas_width::Int, canvas_height::Int) = camera_set!(rm, cam, canvas_width, canvas_height)

mutable struct FitToContentCamera <: Camera
    percent_border::Float64
    FitToContentCamera(percent_border::Float64=0.1) = new(percent_border)
end
function camera_set!(rendermodel::RenderModel, cam::FitToContentCamera, canvas_width::Int, canvas_height::Int)
    camera_fit_to_content!(rendermodel, canvas_width, canvas_height, percent_border=cam.percent_border)
    rendermodel
end
camera_set!(rendermodel::RenderModel, cam::FitToContentCamera, scene::EntityFrame{S,D,I}, roadway::R, canvas_width::Int, canvas_height::Int) where {S,D,I,R} = camera_set!(rendermodel, cam, canvas_width, canvas_height)

# method for new interface
camera_set!(rendermodel::RenderModel, cam::FitToContentCamera, scene, canvas_width::Int, canvas_height::Int) = camera_set!(rendermodel, cam, canvas_width, canvas_height)

mutable struct CarFollowCamera{I} <: Camera
    targetid::I
    zoom::Float64 # [pix/meter]
end
CarFollowCamera(targetid::I) where {I} = CarFollowCamera{I}(targetid, 3.0)

function camera_set!(rendermodel::RenderModel, cam::CarFollowCamera{I}, scene::EntityFrame{S,D,I}, roadway::R, canvas_width::Int, canvas_height::Int) where {S<:State1D,D,I,R}

    veh_index = findfirst(cam.targetid, scene)
    if veh_index != 0
        camera_set_pos!(rendermodel, VecE2(scene[veh_index].state.s, 0.0))
        camera_setzoom!(rendermodel, cam.zoom)
    else
        add_instruction!( rendermodel, render_text, (@sprintf("CarFollowCamera did not find id %d", cam.targetid), 10, 15, 15, colorant"white"), incameraframe=false)
        camera_fit_to_content!(rendermodel, canvas_width, canvas_height)
    end

    rendermodel
end
function camera_set!(rendermodel::RenderModel, cam::CarFollowCamera{I}, scene::EntityFrame{S,D,I}, roadway::R, canvas_width::Int, canvas_height::Int) where {S<:VehicleState,D,I,R}

    veh_index = findfirst(cam.targetid, scene)
    if veh_index != 0
        camera_set_pos!(rendermodel, scene[veh_index].state.posG)
        camera_setzoom!(rendermodel, cam.zoom)
    else
        add_instruction!( rendermodel, render_text, (@sprintf("CarFollowCamera did not find id %d", cam.targetid), 10, 15, 15, colorant"white"), incameraframe=false)
        camera_fit_to_content!(rendermodel, canvas_width, canvas_height)
    end

    rendermodel
end

# method for new interface
function camera_set!(rendermodel::RenderModel, cam::CarFollowCamera, scene, canvas_width::Int, canvas_height::Int)

    inds = findall(x -> x isa ArrowCar && id(x) == cam.targetid, scene)
    if isempty(inds)
        ids = [c.id for c in scene if c isa ArrowCar]
        add_instruction!( rendermodel, render_text, ("CarFollowCamera did not find an ArrowCar with id $(cam.targetid) (found ids: $ids)", 10, 15, 15, colorant"white"), incameraframe=false)
        camera_fit_to_content!(rendermodel, canvas_width, canvas_height)
    else
        veh_index = first(inds)
        camera_set_pos!(rendermodel, pos(scene[veh_index])...)
        camera_setzoom!(rendermodel, cam.zoom)
    end

    rendermodel
end

mutable struct SceneFollowCamera <: Camera
    zoom::Float64 # [pix/meter]
    SceneFollowCamera(zoom::Float64=3.0) = new(zoom)
end
function camera_set!(rendermodel::RenderModel, cam::SceneFollowCamera, scene::EntityFrame{S,D,I}, roadway::R, canvas_width::Int, canvas_height::Int) where {S<:State1D,D,I,R}

    if length(scene) > 0

        # get camera center
        C = 0.0
        for veh in scene
            C += veh.state.s
        end
        C = C / length(scene)

        camera_set_pos!(rendermodel, VecE2(C, 0.0))
        camera_setzoom!(rendermodel, cam.zoom)
    else
        add_instruction!( rendermodel, render_text, ("SceneFollowCamera did not find any vehicles", 10, 15, 15, colorant"white"), incameraframe=false)
        camera_fit_to_content!(rendermodel, canvas_width, canvas_height, 0.1)
    end

    rendermodel
end
function camera_set!(rendermodel::RenderModel, cam::SceneFollowCamera, scene::EntityFrame{S,D,I}, roadway::R, canvas_width::Int, canvas_height::Int) where {S<:VehicleState,D,I,R}

    if length(scene) > 0

        # get camera center
        C = VecE2(0.0,0.0)
        for veh in scene
            C += convert(VecE2, veh.state.posG)
        end
        C = C / length(scene)

        camera_set_pos!(rendermodel, C)
        camera_setzoom!(rendermodel, cam.zoom)
    else
        add_instruction!( rendermodel, render_text, ("SceneFollowCamera did not find any vehicles", 10, 15, 15, colorant"white"), incameraframe=false)
        camera_fit_to_content!(rendermodel, canvas_width, canvas_height, 0.1)
    end

    rendermodel
end


text:
# Autogenerated wrapper script for ZeroMQ_jll for x86_64-w64-mingw32-cxx03
export libzmq

using libsodium_jll
JLLWrappers.@generate_wrapper_header("ZeroMQ")
JLLWrappers.@declare_library_product(libzmq, "libzmq.dll")
function __init__()
    JLLWrappers.@generate_init_header(libsodium_jll)
    JLLWrappers.@init_library_product(
        libzmq,
        "bin\\libzmq.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

if !isdefined(@__MODULE__, :DefsTeam)
    const DefsTeam = String
else
    @warn("Skipping redefinition of DefsTeam to String")
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ONNX"
version = v"1.10.2"

sources = [
    ArchiveSource(
        "https://github.com/onnx/onnx/archive/refs/tags/v$(version).tar.gz",
        "520b3aa34272cc215e2eb41385f58adf01750d88858d4722563edca8410c5dc9",
    ),
]

script = raw"""
cd onnx*
mkdir build
cd build
cmake \
    -DBUILD_ONNX_PYTHON=OFF \
    -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DONNX_USE_LITE_PROTO=ON \
    -DONNX_USE_PROTOBUF_SHARED_LIBS=ON \
    -DProtobuf_PROTOC_EXECUTABLE=${host_bindir}/protoc \
    ..
cmake --build . --config Release --target install -- -j${nproc}
"""

#=
Sadly, no windows support yet. It looks to me like the upstream treats
`defined(_WIN32)` as synonymous with `build with MSVC`.
=#
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true, exclude=Sys.iswindows))

products = [
    FileProduct("lib/libonnx.a", :libonnx),
    LibraryProduct("libonnxifi", :libonnxifi),
    LibraryProduct("libonnxifi_dummy", :libonnxifi_dummy),
]

dependencies = [HostBuildDependency("protoc_jll"), Dependency("protoc_jll")]

build_tarballs(
    ARGS,
    name,
    version,
    sources,
    script,
    platforms,
    products,
    dependencies;
    preferred_gcc_version = v"9",
    julia_compat = "1.6",
)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Plots
using Random
Random.seed!(1)

include_implicit_dynamics()
include_ddp()
include(joinpath(pwd(), "examples/implicit_dynamics/models/cartpole/model.jl"))

# build implicit dynamics
h = 0.05
T = 51

# friction model
data = dynamics_data(model, h, T,
        r_func, rz_func, rθ_func, rz_array, rθ_array;
		idx_soc = idx_soc,
		z_subset_init = z_subset_init,
        diff_idx = 3,
        # θ_params = [1.0; 1.0], # fails to reach constraint tolerance
        # θ_params = [0.5; 0.5], # fails to reach constraint tolerance
        θ_params = [0.35; 0.35],
        # θ_params = [0.25; 0.25],
        # θ_params = [0.1; 0.1],
        # θ_params = [0.01; 0.01],
        opts =  InteriorPointOptions{Float64}(
						r_tol = 1.0e-8,
						κ_tol = 1.0e-4,
						κ_init = 1.0,
						diff_sol = false))

model_implicit = ImplicitDynamics{Midpoint, FixedTime}(2 * model.dim.q, model.dim.u, 0, data)

# no friction model
# data = dynamics_data(model, h,
#         r_no_friction_func, rz_no_friction_func, rθ_no_friction_func,
# 		rz_no_friction_array, rθ_no_friction_array,
#         θ_params = [0.0; 0.0],
#         dyn_opts =  InteriorPointOptions{Float64}(
# 						r_tol = 1.0e-8,
# 						κ_tol = 0.1,
# 						κ_init = 0.1,
# 						diff_sol = true),
# 		jac_opts =  InteriorPointOptions{Float64}(
# 						r_tol = 1.0e-8,
# 						κ_tol = 0.1,
# 						κ_init = 0.1,
# 						diff_sol = true))
#
# model_implicit = ImplicitDynamics{Midpoint, FixedTime}(2 * model.dim.q, model.dim.u, 0, data)

n = model_implicit.n
nq = model_implicit.dynamics.m.dim.q
m = model_implicit.m

q0 = [0.0; 0.0]
q1 = [0.0; 0.0]
qT = [0.0; π]
q_ref = [0.0; π]

x1 = [q1; q1]
xT = [qT; qT]

# Objective
V = 1.0 * Diagonal(ones(nq))
Q_velocity = [V -V; -V V] ./ h^2.0
Q_track = 1.0 * Diagonal(ones(2 * nq))

Q = [t < T ? Q_velocity + Q_track : Q_velocity + 1.0 * Q_track for t = 1:T]
q = [-2.0 * (t == T ? 1.0 : 1.0) * Q_track * xT for t = 1:T]
R = [Diagonal(1.0 * ones(m)) for t = 1:T-1]
r = [zeros(m) for t = 1:T-1]

obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return x' * Q * x + q' * x + u' * R * u + r' * u
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return x' * Q * x + q' * x
    else
        return 0.0
    end
end

# Constraints
ul = [-10.0]
uu = [10.0]
p = [t < T ? 2 * m : n for t = 1:T]
info_t = Dict(:ul => ul, :uu => uu, :inequality => (1:2 * m))
info_T = Dict(:xT => xT)
con_set = [StageConstraint(p[t], t < T ? info_t : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	p = cons.con[t].p

	if t < T
		ul = cons.con[t].info[:ul]
		uu = cons.con[t].info[:uu]
		# c .= [ul - u; u - uu]
	else
		c .= x - cons.con[T].info[:xT]
	end
end

# ū = [(t == 1 ? -1.0 : 0.0) * ones(m) for t = 1:T-1]
ū = [(t == 1 ? 1.0 : 0.0) * ones(m) for t = 1:T-1] # initialize no friciton model this direction (not sure why it goes the oppostive direction...)
w = [zeros(model_implicit.d) for t = 1:T-1]

# Rollout
x̄ = rollout(model_implicit, x1, ū, w, h, T)

q̄ = state_to_configuration(x̄)
# visualize!(vis, s.model, q̄, Δt = h)

prob = problem_data(model_implicit, obj, con_set, copy(x̄), copy(ū), w, h, T,
	analytical_dynamics_derivatives = true)

# Solve
@time stats = constrained_ddp_solve!(prob,
	max_iter = 1000,
    max_al_iter = 10,
	ρ_init = 1.0,
    ρ_scale = 10.0,
    grad_tol = 1.0e-3,
	con_tol = 0.001)

@show ilqr_iterations(stats)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)
b = [data.z_cache[t][collect([model.dim.q + 2, model.dim.q + nc + 2])] for t = 1:T-1]

# x_01, u_01, b_01 = x, u, b
# @save joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/cartpole_01.jld2") x_01 u_01 b_01
@load joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/cartpole_01.jld2") x_01 u_01 b_01

# x_1, u_1, b_1 = x, u, b
# @save joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/cartpole_1.jld2") x_1 u_1 b_1
@load joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/cartpole_1.jld2") x_1 u_1 b_1

# x_25, u_25, b_25 = x, u, b
# @save joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/cartpole_25.jld2") x_25 u_25 b_25
@load joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/cartpole_25.jld2") x_25 u_25 b_25

# x_35, u_35, b_35 = x, u, b
# @save joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/cartpole_35.jld2") x_35 u_35 b_35
@load joinpath(pwd(), "examples/implicit_dynamics/examples/trajectories/cartpole_35.jld2") x_35 u_35 b_35

# using Plots
# t = range(0, stop = h * (T - 1), length = T)
# plot(t, hcat(b_01..., b_01[end])[1, :], linetype = :steppost, color = :magenta, width = 2.0, labels = "")
# plot!(t, hcat(b_1..., b_1[end])[1, :], linetype = :steppost, color = :orange, width = 2.0, labels = "")
# plot!(t, hcat(b_25..., b_25[end])[1, :], linetype = :steppost, color = :green, width = 2.0, labels = "")
# plot!(t, hcat(b_35..., b_35[end])[1, :], linetype = :steppost, color = :cyan, width = 2.0, labels = "")
#
# plot(t, hcat(b_01..., b_01[end])[2, :], linetype = :steppost, color = :magenta, width = 2.0, labels = "")
# plot!(t, hcat(b_1..., b_1[end])[2, :], linetype = :steppost, color = :orange, width = 2.0, labels = "")
# plot!(t, hcat(b_25..., b_25[end])[2, :], linetype = :steppost, color = :green, width = 2.0, labels = "")
# plot!(t, hcat(b_35..., b_35[end])[2, :], linetype = :steppost, color = :cyan, width = 2.0, labels = "")
#
# using PGFPlots
# const PGF = PGFPlots
#
# plt_b_01_s = PGF.Plots.Linear(t, hcat(b_01..., b_01[end])[1, :],
# 	mark="none",style="const plot, color=magenta, line width = 2pt",legendentry="0.01")
#
# plt_b_1_s = PGF.Plots.Linear(t, hcat(b_1..., b_1[end])[1, :],
# 	mark="none",style="const plot, color=orange, line width = 2pt",legendentry="0.1")
#
# plt_b_25_s = PGF.Plots.Linear(t, hcat(b_25..., b_25[end])[1, :],
# 	mark="none",style="const plot, color=green, line width = 2pt",legendentry="0.25")
#
# plt_b_35_s = PGF.Plots.Linear(t, hcat(b_35..., b_35[end])[1, :],
# 	mark="none",style="const plot, color=cyan, line width = 2pt",legendentry="0.35")
#
# plt_b_01_a = PGF.Plots.Linear(t, hcat(b_01..., b_01[end])[2, :],
# 	mark="none",style="const plot, color=magenta, line width = 2pt",legendentry="0.01")
#
# plt_b_1_a = PGF.Plots.Linear(t, hcat(b_1..., b_1[end])[2, :],
# 	mark="none",style="const plot, color=orange, line width = 2pt",legendentry="0.1")
#
# plt_b_25_a = PGF.Plots.Linear(t, hcat(b_25..., b_25[end])[2, :],
# 	mark="none",style="const plot, color=green, line width = 2pt",legendentry="0.25")
#
# plt_b_35_a = PGF.Plots.Linear(t, hcat(b_35..., b_35[end])[2, :],
# 	mark="none",style="const plot, color=cyan, line width = 2pt",legendentry="0.35")
#
# a_s = Axis([plt_b_01_s, plt_b_1_s, plt_b_25_s, plt_b_35_s],
#     axisEqualImage=false,
#     hideAxis=false,
# 	ylabel="friction",
# 	xlabel="time (s)",
# 	# xlim=(0.0, 5.0),
# 	legendStyle="{at={(0.01,0.99)},anchor=north west}")
#
# a_a = Axis([plt_b_01_a, plt_b_1_a, plt_b_25_a, plt_b_35_a],
#     axisEqualImage=false,
#     hideAxis=false,
# 	ylabel="configuration",
# 	xlabel="time (s)",
# 	# xlim=(0.0, 5.0),
# 	legendStyle="{at={(0.01,0.99)},anchor=north west}")
#
# PGF.save("/home/taylor/Research/implicit_dynamics_manuscript/figures/cartpole_friction_slider.tikz", a_s, include_preamble=false)
# PGF.save("/home/taylor/Research/implicit_dynamics_manuscript/figures/cartpole_friction_angle.tikz", a_a, include_preamble=false)

plot(hcat(b..., b[end])', linetype = :steppost)
data.θ_params[1] * (model.mp + model.mc) * model.g * h[1]
data.θ_params[1] * (model.mp * model.g * model.l) * h[1]

q̄ = state_to_configuration(x̄)
v̄ = [(q̄[t+1] - q̄[t]) ./ h for t = 1:length(q̄)-1]

# t = range(0, stop = h * (length(q̄) - 1), length = length(q̄))
# plt = plot();
# plt = plot!(t, hcat(q̄...)', width = 2.0,
# 	color = [:magenta :orange],
# 	labels = ["q1" "q2"],
# 	legend = :topleft,
# 	xlabel = "time (s)",
# 	ylabel = "configuration",
# 	# title = "cartpole (w / o friction)")
# 	title = "cartpole (w/ friction)")
#
# plt = plot();
# plt = plot!(t, hcat(v̄..., v̄[end])', width = 2.0,
# 	color = [:magenta :orange],
# 	labels = ["q1" "q2"],
# 	legend = :topleft,
# 	xlabel = "time (s)",
# 	ylabel = "velocity",
# 	linetype = :steppost,
# 	# title = "cartpole (w / o friction)")
# 	title = "cartpole (w/ friction)")
# 	# title = "acrobot (w/ joint limits)")

# show(plt)
# savefig(plt, "/home/taylor/Research/implicit_dynamics_manuscript/figures/cartpole_friction.png")
# savefig(plt, "/home/taylor/Research/implicit_dynamics_manuscript/figures/cartpole_no_friction.png")

# plot(hcat(ū..., ū[end])', linetype = :steppost)
# #
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model, q̄, Δt = h)
open(vis)
default_background!(vis)
settransform!(vis["/Cameras/default"],
        compose(Translation(0.0, -95.0, -1.0), LinearMap(RotY(0.0 * π) * RotZ(-π / 2.0))))
setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 50)
setvisible!(vis["/Grid"], false)
#
# # q̄ = state_to_configuration(x̄)
# # q_anim = [[q̄[1] for t = 1:20]..., q̄..., [q̄[end] for t = 1:20]...]
# # visualize!(vis, model, q_anim, Δt = h)
#
# # ghost
# limit_color = [0.0, 0.0, 0.0]
# # limit_color = [0.0, 1.0, 0.0]
#
# t = 1
# id = t
# tl = 0.05
# _create_cartpole!(vis, model;
#         tl = tl,
#         color = RGBA(limit_color..., tl),
#         i = id)
# _set_cartpole!(vis, model, x[t], i = id)
#
# t = 5
# id = t
# tl = 0.15
# _create_cartpole!(vis, model;
#         tl = tl,
#         color = RGBA(limit_color..., tl),
#         i = id)
# _set_cartpole!(vis, model, x[t], i = id)
#
# t = 10
# id = t
# tl = 0.25
# _create_cartpole!(vis, model;
#         tl = tl,
#         color = RGBA(limit_color..., tl),
#         i = id)
# _set_cartpole!(vis, model, x[t], i = id)
#
# t = 15
# id = t
# tl = 0.35
# _create_cartpole!(vis, model;
#         tl = tl,
#         color = RGBA(limit_color..., tl),
#         i = id)
# _set_cartpole!(vis, model, x[t], i = id)
#
# t = 20
# id = t
# tl = 0.45
# _create_cartpole!(vis, model;
#         tl = tl,
#         color = RGBA(limit_color..., tl),
#         i = id)
# _set_cartpole!(vis, model, x[t], i = id)
#
# t = 25
# id = t
# tl = 0.55
# _create_cartpole!(vis, model;
#         tl = tl,
#         color = RGBA(limit_color..., tl),
#         i = id)
# _set_cartpole!(vis, model, x[t], i = id)
#
# t = 30
# id = t
# tl = 0.65
# _create_cartpole!(vis, model;
#         tl = tl,
#         color = RGBA(limit_color..., tl),
#         i = id)
# _set_cartpole!(vis, model, x[t], i = id)
#
# t = 35
# id = t
# tl = 0.75
# _create_cartpole!(vis, model;
#         tl = tl,
#         color = RGBA(limit_color..., tl),
#         i = id)
# _set_cartpole!(vis, model, x[t], i = id)
#
# t = 40
# id = t
# tl = 0.85
# _create_cartpole!(vis, model;
#         tl = tl,
#         color = RGBA(limit_color..., tl),
#         i = id)
# _set_cartpole!(vis, model, x[t], i = id)
#
# t = 45
# id = t
# tl = 0.95
# _create_cartpole!(vis, model;
#         tl = tl,
#         color = RGBA(limit_color..., tl),
#         i = id)
# _set_cartpole!(vis, model, x[t], i = id)
#
#
# t = 51
# id = t
# tl = 1.0
# _create_cartpole!(vis, model;
#         tl = tl,
#         color = RGBA(limit_color..., tl),
#         i = id)
# _set_cartpole!(vis, model, x[t], i = id)
#
# line_mat = LineBasicMaterial(color=color=RGBA(1.0, 153.0 / 255.0, 51.0 / 255.0, 1.0), linewidth=10.0)
# # line_mat = LineBasicMaterial(color=color=RGBA(51.0 / 255.0, 1.0, 1.0, 1.0), linewidth=10.0)
#
# points = Vector{Point{3,Float64}}()
# for (i, xt) in enumerate(x̄)
#     k = kinematics(model, xt)
# 	push!(points, Point(k[1], 0.0, k[2]))
#
#     setobject!(vis["ee_vertex_$i"], Sphere(Point3f0(0),
#         convert(Float32, 0.001)),
#         MeshPhongMaterial(color = RGBA(1.0, 153.0 / 255.0, 51.0 / 255.0, 1.0)))
#         settransform!(vis["ee_vertex_$i"], Translation(points[i]))
# end
# setobject!(vis[:ee_traj], MeshCat.Line(points, line_mat))
#
#
# # using PGFPlots
# # const PGF = PGFPlots
# #
# # plt_q1_smooth = PGF.Plots.Linear(t, hcat(q̄_smooth...)[1,:],
# # 	mark="none",style="color=cyan, line width = 2pt, dashed",legendentry="q1")
# #
# # plt_q2_smooth = PGF.Plots.Linear(t, hcat(q̄_smooth...)[2,:],
# # 	mark="none",style="color=orange, line width = 2pt, dashed",legendentry="q2")
# #
# # plt_qd1_smooth = PGF.Plots.Linear(t, hcat(v̄_smooth..., v̄_smooth[end])[1,:],
# # 	mark="none",style="const plot, color=magenta, line width = 2pt, dashed",legendentry="q1")
# #
# # plt_qd2_smooth = PGF.Plots.Linear(t, hcat(v̄_smooth..., v̄_smooth[end])[2,:],
# # 	mark="none",style="const plot, color=green, line width = 2pt, dashed",legendentry="q2")
# #
# # plt_q1_friction = PGF.Plots.Linear(t, hcat(q̄_friction...)[1,:],
# # 	mark="none",style="color=cyan, line width = 2pt",legendentry="q1 (friction)")
# #
# # plt_q2_friction = PGF.Plots.Linear(t, hcat(q̄_friction...)[2,:],
# # 	mark="none",style="color=orange, line width = 2pt",legendentry="q2 (friction)")
# #
# # plt_qd1_friction = PGF.Plots.Linear(t, hcat(v̄_friction..., v̄_friction[end])[1,:],
# # 	mark="none",style="const plot, color=magenta, line width = 2pt",legendentry="q1 (friction)")
# #
# # plt_qd2_friction = PGF.Plots.Linear(t, hcat(v̄_friction..., v̄_friction[end])[2,:],
# # 	mark="none",style="const plot, color=green, line width = 2pt",legendentry="q2 (friction)")
# #
# # aq = Axis([plt_q1_friction; plt_q2_friction; plt_q1_smooth; plt_q2_smooth],#; plt_qd1_smooth; plt_qd1_friction; plt_qd2_smooth; plt_qd2_friction],
# #     axisEqualImage=false,
# #     hideAxis=false,
# # 	ylabel="configuration",
# # 	xlabel="time (s)",
# # 	xlims=(0.0, 2.5),
# # 	legendStyle="{at={(0.01,0.99)},anchor=north west}")
# #
# # av = Axis([plt_qd1_friction; plt_qd2_friction; plt_qd1_smooth; plt_qd2_smooth],
# #     axisEqualImage=false,
# #     hideAxis=false,
# # 	ylabel="velocity",
# # 	xlabel="time (s)",
# # 	legendStyle="{at={(0.01,0.99)},anchor=north west}")
# #
# # # Save to tikz format
# # PGF.save("/home/taylor/Research/implicit_dynamics_manuscript/figures/cartpole_friction_configuration.tikz", aq, include_preamble=false)
# # PGF.save("/home/taylor/Research/implicit_dynamics_manuscript/figures/cartpole_friction_velocity.tikz", av, include_preamble=false)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
## hclust.jl (c) 2014 David A. van Leeuwen
## Hierarchical clustering, similar to R's hclust()

## Algorithms are based upon C. F. Olson, Parallel Computing 21 (1995) 1313--1325.

## This is also in types.jl, but that is not read...
## Mostly following R's hclust class
type Hclust{T<:Real}
    merge::Matrix{Int}
    height::Vector{T}
    order::Vector{Int}
    labels::Vector
    method::Symbol
end

function assertdistancematrix(d::AbstractMatrix)
    nr, nc = size(d)
    nr == nc || throw(DimensionMismatch("Distance matrix should be square."))
    issym(d) || error("Distance matrix should be symmetric.")
end

## This seems to work like R's implementation, but it is extremely inefficient
## This probably scales O(n^3) or worse. We can use it to check correctness
function hclust_n3{T<:Real}(d::AbstractMatrix{T}, method::Function)
    assertdistancematrix(d)
    mr = Int[]                  # min row
    mc = Int[]                  # min col
    h = T[]                     # height
    nc = size(d,1)              # number of clusters
    cl = -[1:nc]                # segment to cluster attribution, initially negative
    next = 1                    # next cluster label
    while next < nc
        mindist = Inf
        mi = mj = 0
        cli = unique(cl)
        mask = BitVector(nc)
        for j in 1:length(cli)           # loop over for lower triangular indices, i>j
            cols = cl .== cli[j]
            for i in (j+1):length(cli)
                rows = cl.==cli[i]
                distance = method(d[rows,cols]) # very expensive
                if distance < mindist
                    mindist = distance
                    mi = cli[i]
                    mj = cli[j]
                    mask = cols | rows
                end
            end
        end
        ## simulate R's order
        if mi < 0 && mj < 0 && mi > mj ||
            mi > 0 && mj > 0 && mi < mj ||
            mi < 0 && mj > 0
            push!(mr, mi)
            push!(mc, mj)
        else
            push!(mr, mj)
            push!(mc, mi)
        end
        push!(h, mindist)
        cl[mask] = next
        next += 1
    end
    hcat(mr, mc), h
end

## Efficient single link algorithm, according to Olson, O(n^2), fig 2.
## Verified against R's implementation, correct, and about 2.5 x faster
## For each i < j compute D(i,j) (this is already given)
## For each 0 < i <= n compute Nearest Neighbor N(i)
## Repeat n-1 times
##   find i,j that minimize D(i,j)
##   merge clusters i and j
##   update D(i,j) and N(i) accordingly
function hclust_minimum{T<:Real}(ds::Symmetric{T})
    ## For each i < j compute d[i,j] (this is already given)
    d = full(ds)                #  we need a local copy
    nc = size(d,1)
    mr = Array(Int, nc-1)       # min row
    mc = Array(Int, nc-1)       # min col
    h = Array(T, nc-1)          # height
    merges = -[1:nc]
    next = 1
    ## For each 0 < i <= n compute Nearest Neighbor N[i]
    N = zeros(Int, nc)
    for k = 1:nc
        mindist = Inf
        mk = 0
        for i = 1:(k-1)
            if d[i,k] < mindist
                mindist = d[i,k]
                mk = i
            end
        end
        for j = (k+1):nc
            if d[k,j] < mindist
                mindist = d[k,j]
                mk = j
            end
        end
        N[k] = mk
    end
    ## the main loop
    while nc > 1                # O(n)
        mindist = d[1,N[1]]
        i = 1
        for k in 2:nc           # O(n)
            if k < N[k]
                distance = d[k,N[k]]
            else
                distance = d[N[k],k]
            end
            if distance < mindist
                mindist = distance
                i = k
            end
        end
        j = N[i]
        if i > j
            i, j = j, i     # make sure i < j
        end
        ## update result, compatible to R's order.  It must be possible to do this simpler than this...
        if merges[i] < 0 && merges[j] < 0 && merges[i] > merges[j] ||
            merges[i] > 0 && merges[j] > 0 && merges[i] < merges[j] ||
            merges[i] < 0 && merges[j] > 0
            mr[next] = merges[i]
            mc[next] = merges[j]
        else
            mr[next] = merges[j]
            mc[next] = merges[i]
        end
        h[next] = mindist
        merges[i] = next
        merges[j] = merges[nc]
        ## update d, split in ranges k<i, i<k<j, j<k<=nc
        for k = 1:(i-1)         # k < i
            if d[k,i] > d[k,j]
                d[k,i] = d[k,j]
            end
        end
        for k = (i+1):(j-1)     # i < k < j
            if d[i,k] > d[k,j]
                d[i,k] = d[k,j]
            end
        end
        for k = (j+1):nc        # j < k <= nc
            if d[i,k] > d[j,k]
                d[i,k] = d[j,k]
            end
        end
        ## move the last row/col into j
        for k=1:(j-1)           # k==nc,
            d[k,j] = d[k,nc]
        end
        for k=(j+1):(nc-1)
            d[j,k] = d[k,nc]
        end
        ## update N[k], k !in (i,j)
        for k=1:nc
            if N[k] == j       # update nearest neigbors != i
                N[k] = i
            elseif N[k] == nc
                N[k] = j
            end
        end
        N[j] = N[nc]            # update N[j]
        ## update nc, next
        nc -= 1
        next += 1
        ## finally we need to update N[i], because it was nearest to j
        mindist = Inf
        mk = 0
        for k=1:(i-1)
            if d[k,i] < mindist
                mindist = d[k,i]
                mk = k
            end
        end
        for k = (i+1):nc
            if d[i,k] < mindist
                mindist = d[i,k]
                mk = k
            end
        end
        N[i] = mk
#        for n in N[1:nc] print(n, " ") end; println()
    end
    return hcat(mr, mc), h
end

## functions to compute maximum, minimum, mean for just a slice of an array

function slicemaximum{T<:Real}(d::AbstractMatrix{T}, cl1::Vector{Int}, cl2::Vector{Int})
    maxdist = -Inf
    for i in cl1 for j in cl2
        if d[i,j] > maxdist
            maxdist = d[i,j]
        end
    end end
    maxdist
end

function sliceminimum{T<:Real}(d::AbstractMatrix{T}, cl1::Vector{Int}, cl2::Vector{Int})
    mindist = Inf
    for i in cl1 for j in cl2
        if d[i,j] < mindist
            mindist = d[i,j]
        end
    end end
    mindist
end

function slicemean{T<:Real}(d::AbstractMatrix{T}, cl1::Vector{Int}, cl2::Vector{Int})
    s = zero(T)
    for i in cl1 for j in cl2
        s += d[i,j]
    end end
    s / (length(cl1)*length(cl2))
end

## This reorders the pairs to be compatible with R's hclust()
function rorder!(mr, mc, h)
    o = sortperm(h)
    io = invperm(o)
    for i in 1:length(mr)
        if mr[i] > 0
            mr[i] = io[mr[i]]
        end
        if mc[i] > 0
            mc[i] = io[mc[i]]
        end
        ## R's order of pairs
        if ! (mr[i] < 0 && mc[i] < 0 && mr[i] > mc[i] ||
              mr[i] > 0 && mc[i] > 0 && mr[i] < mc[i] ||
              mr[i] < 0 && mc[i] > 0)
            mr[i], mc[i] = mc[i], mr[i]
        end
    end
    return o
end

## Another nearest neighbor algorithm, for reducible metrics
## From C. F. Olson, Parallel Computing 21 (1995) 1313--1325, fig 5
## Verfied against R implementation for mean and maximum, correct but ~ 5x slower
## Pick c1: 0 <= c1 <= n random
## i <- 1
## repeat n-1 times
##   repeat
##     i++
##     c[i] = nearest neigbour c[i-1]
##   until c[i] = c[i-2] ## nearest of nearest is cluster itself
##   merge c[i] and nearest neigbor c[i]
##   if i>3 i -= 3 else i <- 1
function hclust2{T<:Real}(d::Symmetric{T}, method::Function)
    nc = size(d,1)                      # number of clusters
    mr = Array(Int, nc-1)               # min row
    mc = Array(Int, nc-1)               # min col
    h = Array(T, nc-1)                  # height
    cl = [[x] for x in 1:nc]            # clusters
    merges = -[1:nc]
    next = 1
    i = 1
    N = Array(Int, nc+1)
    N[1] = 1                            # arbitrary choice
    while nc > 1
        found=false
        mindist = Inf
        while !found
            i += 1
            mi = 0
            mindist = Inf
            Nim1 = N[i-1]
            ## c[i] = nearest neigbour c[i-1]
            for j = 1:nc if Nim1 != j
                distance = method(d, cl[Nim1], cl[j])
                if distance < mindist
                    mindist = distance
                    mi = j
                end
            end end
            N[i] = mi           # N[i+1] is nearest neigbor to N[i]
            found = i > 2 && N[i] == N[i-2]
        end
        ## merge c[i] and nearest neigbor c[i], i.e., c[i-1]
        if N[i-1] < N[i]
            lo, high = N[i-1], N[i]
        else
            lo, high = N[i], N[i-1]
        end
        ## first, store the result
        mr[next] = merges[lo]
        mc[next] = merges[high]
        h[next] = mindist
        merges[lo] = next
        merges[high] = merges[nc]
        next += 1
        ## then perform the actual merge
        cl[lo] = vcat(cl[lo], cl[high])
        cl[high] = cl[nc]
        if i>3
            i -= 3
        else
            i = 1
        end
        ## replace any nearest neighbor referring to nc
        for k=1:i
            if N[k] == nc
                N[k] = high
            end
        end
        nc -= 1
    end
    ## fix order for presenting result
    o = rorder!(mr, mc, h)
    hcat(mr[o], mc[o]), h[o]
end

## this calls the routine that gives the correct answer, fastest
## method names are inspired by R's hclust
function hclust{T<:Real}(d::Symmetric{T}, method::Symbol)
    nc = size(d,1)
    if method == :single
        h = hclust_minimum(d)
    elseif method == :complete
        h = hclust2(d, slicemaximum)
    elseif method == :average
        h = hclust2(d, slicemean)
    else
        error("Unsupported method ", method)
    end
    ## order and label are placeholders for the moment
    Hclust(h..., [1:nc], [1:nc], method)
end

## uplo may be Char for v0.3, Symbol for v0.4
hclust{T<:Real}(d::AbstractMatrix{T}, method::Symbol, uplo) = hclust(Symmetric(d, uplo), method)

function hclust{T<:Real}(d::AbstractMatrix{T}, method::Symbol)
    assertdistancematrix(d)
    hclust(Symmetric(d), method)
end

## cut a tree at height `h' or to `k' clusters
function cutree(hclust::Hclust; k::Int=1,
                h::Real=maximum(hclust.height))
    clusters = Vector{Int}[]
    nnodes = length(hclust.labels)
    nodes = [[i::Int] for i=1:nnodes]
    N = nnodes - k
    i = 1
    while i<=N && hclust.height[i] <= h
        both = vec(hclust.merge[i,:])
        new = Int[]
            for x in both
                if x<0
                    push!(new, -x)
                    nodes[-x] = []
                else
                    append!(new, clusters[x])
                    clusters[x] = []
                end
            end
        push!(clusters, new)
        i += 1
    end
    all = vcat(clusters, nodes)
    all = all[map(length, all) .> 0]
    ## convert to a single array of cluster indices
    res = Array(Int, nnodes)
    for (i,cl) in enumerate(all)
        res[cl] = i
    end
    res
end

## some diagnostic functions, not exported
function printupper(d::Matrix)
    n = size(d,1)
    for i = 1:(n-1)
        print(" " ^ ((i-1) * 6))
        for j = (i+1):n
            print(@sprintf("%5.2f ", d[i,j]))
        end
        println()
    end
end


text:
# size_tools.jl

export collect_dim, select_sizes, select_sizes_squeeze, center_position, center_value
export expand_add, expand_size, optional_pos_z, reorient, single_dim_size

"""
    contains all functions that operate on Tuples, sizes and alike
"""

# Functions from ROIViews:

"""
    expand_add(t1,t2)  # adds t1 to t2 as a tuple and returns t2[n] for n > length(t1)

Adds the elements of the tuple `t1` to `t2`.
If `t1` is shorter than `t2`, we take the tail of `t2` as elements.
Output length is always the same length as `t2`

Example:
```jldoctest
julia> expand_add((1,2,3),(4,5,6,7,8,9))
(5, 7, 9, 7, 8, 9)

julia> expand_add((1,1), (0,))
(1,)

julia> expand_add((1,2), (0,0))
(1, 2)

julia> expand_add((1,2), (0,0,0))
(1, 2, 0)
```
"""
function expand_add(t1::NTuple{M, T1},t2::NTuple{N, T2}) where {M, N, T1, T2}
    return ntuple(i -> i ≤ M ? t1[i] + t2[i] : t2[i], Val(N))
end

"""
    expand_size(sz,sz2)

Expands a size tuple `sz` with the sizes as given in the tuple `sz2` for positions which do not exist in `sz`.
Typically one wants to

Example:
```jldoctest
julia> sz = expand_size((1,2,3),(4,5,6,7,8,9))
Base.Generator{UnitRange{Int64}, var"#7#8"{Tuple{Int64, Int64, Int64}, NTuple{6, Int64}, Int64}}(var"#7#8"{Tuple{Int64, Int64, Int64}, NTuple{6, Int64}, Int64}((1, 2, 3), (4, 5, 6, 7, 8, 9), 3), 1:6)

julia> Tuple(sz)
(1, 2, 3, 7, 8, 9)
```
"""
function expand_size(sz::NTuple{M, T1},sz2::NTuple{N, T2}) where {M, N, T1, T2}
    return ntuple(i -> i ≤ M ? sz[i] : sz2[i], Val(N))
end

# Functions from IndexFunArrays:

"""
    optional_pos_z(x, offset)

Returns a `z` position `posZ = x[3]-offset[3]` if the tuple is long enough and `posZ=1` otherwise.
This is useful for 3D routines, which should also work for 2d data.

Example:
```jldoctest
julia> optional_pos_z((5,5),(2,3))
1

julia> optional_pos_z((5,5,5),(3,2,1))
4
```
"""
optional_pos_z(x::NTuple{1,T}, offset::NTuple{1,T}) where {T,N} = 1
optional_pos_z(x::NTuple{2,T}, offset::NTuple{2,T}) where {T,N} = 1
optional_pos_z(x::NTuple{N,T}, offset::NTuple{N,T}) where {T,N} = x[3]-offset[3]

"""
    curry(f, x) = (xs...) -> f(x, xs...)   # just a shorthand to remove x

allows to remove the fix the first argument x in a function f with any number of arguments

Example:
```jldoctest
julia> g = curry(+,10.0)
#9 (generic function with 1 method)

julia> g(3)
13.0
```
"""
curry(f, x) = (xs...) -> f(x, xs...)   # just a shorthand to remove x

# Functions for IndexFunArrays.utils

"""
    single_dim_size(dim::Int, dim_size::Int, tdim=dim)

Returns a tuple (length `tdim`, which by default is `dim`) of singleton sizes
except at the final position `dim`, which contains `dim_size`.

!!! warning "Not type-stable!"
    Is not type stable!

Arguments:
+ dim: non-zero position
+ dim_size: the value this non-zero position is given in the returned NTuple
+ tdim: total length of the returned NTuple

Example
```jldoctest
julia> single_dim_size(4, 3)
(1, 1, 1, 3)

julia> single_dim_size(4, 5)
(1, 1, 1, 5)

julia> single_dim_size(2, 5)
(1, 5)

julia> single_dim_size(3,5, 4)
(1, 1, 5, 1)
```
"""
function single_dim_size(dim::Int, dim_size::Int, tdim=dim)
    Base.setindex(ntuple(i -> 1, Val(tdim)), dim_size, dim)::NTuple{tdim, Int}
end

"""
    single_dim_size(::Val{dim}, dim_size::Int, tdim=::Val{dim})

Same `single_dim_size` but type stable with `Val`.
"""
function single_dim_size(::Val{dim}, dim_size::Int, tdim=Val(dim)) where dim
    Base.setindex(ntuple(i -> 1, tdim), dim_size, dim)
end

"""
    reorient(vec, dim)

Reorients a 1D vector `vec` along dimension `dim`.

!!! warning "Not type-stable!"
    Is not type stable!

```julia
julia> reorient([1,2,3,4], 3)
1×1×4 Array{Int64, 3}:
[:, :, 1] =
 1

[:, :, 2] =
 2

[:, :, 3] =
 3

[:, :, 4] =
 4
```
"""
function reorient(vec, dim::Int)
    reshape(vec, single_dim_size(dim, length(vec)))
end

"""
    reorient(vec, ::Val{dim})

Type stable version of `reorient`!

```jldoctest
reorient([1,2,3,4], Val(3))
1×1×4 Array{Int64, 3}:
[:, :, 1] =
 1

[:, :, 2] =
 2

[:, :, 3] =
 3

[:, :, 4] =
 4
```
"""
function reorient(vec, ::Val{dim}) where dim
    reshape(vec, single_dim_size(Val(dim), length(vec)))
end

"""
    collect_dim(col, dim::Int)

collects a collection `col` and reorients it into direction `dim`.

Example:
 ```jldoctest
julia> collect_dim(1:5,2)
1×5 Matrix{Int64}:
 1  2  3  4  5
```
"""
function collect_dim(col, dim::Int)
    reorient(collect(col), dim)
end

"""
    collect_dim(col, dim::Val)

Type stable version!

```jldoctest
julia> collect_dim(1:5, Val(2))
1×5 Matrix{Int64}:
 1  2  3  4  5
```

"""
function collect_dim(col, dim::Val)
    reorient(collect(col), dim)
end

## Functions from FourierTools utils:

"""
    select_sizes(x::AbstractArray, dim)

Additional size method to access the size at several dimensions
in one call.
Keep singleton dimensions.

Examples
```jldoctest
julia> x = ones((2,4,6,8, 10));

julia> select_sizes(x, (2,3))
(1, 4, 6, 1, 1)

julia> select_sizes(x, 5)
(1, 1, 1, 1, 10)

julia> select_sizes(x, (5,))
(1, 1, 1, 1, 10)
```
"""
function select_sizes(x::AbstractArray{T, M}, dim::NTuple{N,Int}) where {T,N,M}
    sz = ntuple(i -> 1, Val(M))
    for n in dim
        sz = Base.setindex(sz, size(x, n), n)
    end
    return Tuple(sz)
end

function select_sizes(x::AbstractArray, dim::Integer)
    select_sizes(x, (dim,))
end

"""
    select_sizes_squeeze(x::AbstractArray, dim)

Additional size method to access the size at several dimensions
in one call.
Remove singleton dimensions.

See also [`select_sizes`](@ref select_sizes) which does not remove singleton dimensions.

Examples
```jldoctest
julia> select_sizes_squeeze(randn((5,6,7)), (2,3))
(6, 7)

julia> select_sizes_squeeze(randn((5,6,7)), 2)
(6,)
```
"""
function select_sizes_squeeze(x::AbstractArray{T, M}, dim::NTuple{N,Int}) where {T,N,M}
    return map(n -> size(x,n), dim)
end

function select_sizes_squeeze(x::AbstractArray, dim::Integer)
    select_sizes_squeeze(x, (dim,))
end

"""
    center_position(field)

position of the center of the `field` according to the Fourier-space nomenclature (pixel right of center for even-sized arrays).
returns a tuple of ints which can be used for indexing.
"""
function center_position(field)
    (size(field) .÷ 2).+1
end

"""
    center_value(field)

Value of the `field` at the position of the center of the `field` according to the Fourier-space nomenclature (pixel right of center for even-sized arrays).
"""
function center_value(field)
    field[center_position(field)...]
end

"""
    ft_center_diff(s [, dims])

Calculates how much each dimension must be shifted that the
center frequency is at the Fourier center.

"""
function ft_center_diff(s::NTuple{N, T}, dims=ntuple(identity, Val(N))) where {N, T}
    ntuple(i -> i ∈ dims ?  s[i] ÷ 2 : 0 , N)
end